# DFT is too slow for regtest, but code is here for show
@SET QS_METHOD SE
# uncomment for DFT - too slow for regtest, though
# @SET QS_METHOD DFT
# USER set the DFT functional, e.g. BLYP or PBE
@SET DFT_FUNCTIONAL BLYP

&FORCE_EVAL
  METHOD QMMM

  &DFT
    WFN_RESTART_FILE_NAME wfn.restart
    CHARGE 0

@IF ${QS_METHOD} == SE
    &QS
      METHOD AM1
      &SE
      &END SE
    &END QS
    &SCF
      SCF_GUESS ATOMIC
      EPS_SCF 1.0E-7
      &MIXING
        ALPHA 0.1
      &END MIXING
      &PRINT
        &RESTART off
        &END RESTART
      &END PRINT
    &END SCF
@ENDIF

@IF ${QS_METHOD} == DFT
    BASIS_SET_FILE_NAME BASIS_SET
    POTENTIAL_FILE_NAME POTENTIAL
    &MGRID
      COMMENSURATE
      CUTOFF 280.0
      NGRIDS 4
      REL_CUTOFF 30.0
    &END MGRID
    &POISSON
      PERIODIC XYZ
      POISSON_SOLVER PERIODIC
      &EWALD
        EWALD_TYPE EWALD
        EWALD_ACCURACY 1.E-6
        ALPHA 0.35
  GMAX 53 61 67
      &END EWALD
    &END POISSON
    &QS
      EXTRAPOLATION PS
      EXTRAPOLATION_ORDER 1
      EPS_DEFAULT 1.0E-12
    &END QS
    &SCF
      SCF_GUESS restart
      EPS_SCF 1.0E-6
      MAX_SCF 20
      &OUTER_SCF
         EPS_SCF 1.0E-6
         MAX_SCF 100
      &END
      &OT
         MINIMIZER CG
         PRECONDITIONER FULL_ALL
      &END OT
    &END SCF
    &XC
      &XC_FUNCTIONAL ${DFT_FUNCTIONAL}
      &END XC_FUNCTIONAL
      &XC_GRID
        XC_DERIV NN10_SMOOTH
        XC_SMOOTH_RHO NN10
      &END XC_GRID
    &END XC

@ENDIF
  &END DFT

  &MM
    &FORCEFIELD
      PARMTYPE AMBER
      PARM_FILE_NAME ../sample_top/Lysozyme_small.top
      VDW_SCALE14 0.5
      EI_SCALE14  0.8333333
      SHIFT_CUTOFF .FALSE.
      &SPLINE
   EMAX_SPLINE 1.0
      &END SPLINE
    &END FORCEFIELD
    &POISSON
@IF ${QS_METHOD} == DFT
      PERIODIC XYZ
      POISSON_SOLVER PERIODIC
@ENDIF
      &EWALD
        EWALD_TYPE SPME
        EWALD_ACCURACY 1.E-6
        ALPHA 0.35
  GMAX 53 61 67
      &END EWALD
    &END POISSON
  &END MM
  &QMMM

      &QM_KIND C
   MM_INDEX 544
      &END QM_KIND

      &QM_KIND H
   MM_INDEX 547
      &END QM_KIND

      &QM_KIND O
   MM_INDEX 545
   MM_INDEX 546
      &END QM_KIND

      &LINK
   QM_INDEX 546
   MM_INDEX 543
      &END LINK

     &FORCE_MIXING

   MAX_N_QM 500

   R_CORE 1.3 1.5
   R_QM 3.0 3.5
   R_BUF 3.0 3.5
   QM_KIND_ELEMENT_MAPPING H H
   QM_KIND_ELEMENT_MAPPING O O
   QM_KIND_ELEMENT_MAPPING N N
   QM_KIND_ELEMENT_MAPPING C C

   EXTENDED_DELTA_CHARGE -1

   ADAPTIVE_EXCLUDE_MOLECULES MOL0

   &QM_NON_ADAPTIVE

      &QM_KIND C
         MM_INDEX 536
         MM_INDEX 538
         MM_INDEX 541
      &END QM_KIND

      &QM_KIND H
         MM_INDEX 537
         MM_INDEX 539
         MM_INDEX 540
         MM_INDEX 542
         MM_INDEX 543
      &END QM_KIND

   &END QM_NON_ADAPTIVE

   &BUFFER_NON_ADAPTIVE

      &QM_KIND C
         MM_INDEX 474
         MM_INDEX 480
         MM_INDEX 484
         MM_INDEX 516
         MM_INDEX 532
         MM_INDEX 548
         MM_INDEX 552
         MM_INDEX 778
         MM_INDEX 781
         MM_INDEX 835
         MM_INDEX 850
         MM_INDEX 854
         MM_INDEX 856
         MM_INDEX 859
         MM_INDEX 862
         MM_INDEX 867
         MM_INDEX 871
         MM_INDEX 1592
         MM_INDEX 1594
         MM_INDEX 1597
         MM_INDEX 1598
         MM_INDEX 1602
         MM_INDEX 1603
         MM_INDEX 1605
         MM_INDEX 1607
         MM_INDEX 1609
         MM_INDEX 1611
         MM_INDEX 1612
         MM_INDEX 1616
         MM_INDEX 1628
         MM_INDEX 1632
         MM_INDEX 1634
      &END QM_KIND

      &QM_KIND H
         MM_INDEX 475
         MM_INDEX 483
         MM_INDEX 485
         MM_INDEX 517
         MM_INDEX 535
         MM_INDEX 551
         MM_INDEX 553
         MM_INDEX 779
         MM_INDEX 780
         MM_INDEX 836
         MM_INDEX 853
         MM_INDEX 855
         MM_INDEX 857
         MM_INDEX 858
         MM_INDEX 860
         MM_INDEX 861
         MM_INDEX 865
         MM_INDEX 866
         MM_INDEX 870
         MM_INDEX 872
         MM_INDEX 1593
         MM_INDEX 1595
         MM_INDEX 1596
         MM_INDEX 1599
         MM_INDEX 1601
         MM_INDEX 1604
         MM_INDEX 1606
         MM_INDEX 1608
         MM_INDEX 1610
         MM_INDEX 1615
         MM_INDEX 1617
         MM_INDEX 1631
         MM_INDEX 1633
         MM_INDEX 1635
         MM_INDEX 1636
         MM_INDEX 1637
      &END QM_KIND

      &QM_KIND N
         MM_INDEX 482
         MM_INDEX 534
         MM_INDEX 550
         MM_INDEX 852
         MM_INDEX 864
         MM_INDEX 869
         MM_INDEX 1600
         MM_INDEX 1614
         MM_INDEX 1630
      &END QM_KIND

      &QM_KIND O
         MM_INDEX 481
         MM_INDEX 533
         MM_INDEX 549
         MM_INDEX 782
         MM_INDEX 783
         MM_INDEX 851
         MM_INDEX 863
         MM_INDEX 868
         MM_INDEX 1613
         MM_INDEX 1629
      &END QM_KIND

      &LINK
         QM_INDEX 474
         MM_INDEX 472
      &END LINK
      &LINK
         QM_INDEX 474
         MM_INDEX 476
      &END LINK
      &LINK
         QM_INDEX 484
         MM_INDEX 490
      &END LINK
      &LINK
         QM_INDEX 484
         MM_INDEX 486
      &END LINK
      &LINK
         QM_INDEX 516
         MM_INDEX 514
      &END LINK
      &LINK
         QM_INDEX 516
         MM_INDEX 518
      &END LINK
      &LINK
         QM_INDEX 552
         MM_INDEX 559
      &END LINK
      &LINK
         QM_INDEX 552
         MM_INDEX 554
      &END LINK
      &LINK
         QM_INDEX 778
         MM_INDEX 776
      &END LINK
      &LINK
         QM_INDEX 835
         MM_INDEX 833
      &END LINK
      &LINK
         QM_INDEX 835
         MM_INDEX 837
      &END LINK
      &LINK
         QM_INDEX 871
         MM_INDEX 886
      &END LINK
      &LINK
         QM_INDEX 871
         MM_INDEX 873
      &END LINK
      &LINK
         QM_INDEX 1592
         MM_INDEX 1590
      &END LINK
      &LINK
         QM_INDEX 1616
         MM_INDEX 1618
      &END LINK
      &LINK
         QM_INDEX 1632
         MM_INDEX 1638
      &END LINK

   &END BUFFER_NON_ADAPTIVE

     &END FORCE_MIXING

     &CELL
  # 19x16x18 + 4 vacuum on each side
        ABC 27.0 24.0 26.0
@IF ${QS_METHOD} == DFT
  PERIODIC XYZ
@ENDIF
     &END CELL
@IF ${QS_METHOD} == SE
# SE only supports COULOMB
    ECOUPL COULOMB
@ENDIF
@IF ${QS_METHOD} == DFT
#GAUSS is from Laino 2005 using Blochl 1995, Gaussian expansion of the electrostatic potential
    ECOUPL GAUSS
@ENDIF
    USE_GEEP_LIB 10
    &MM_KIND OX
      RADIUS 1.2
    &END MM_KIND
    &MM_KIND HX
      RADIUS 0.44
    &END MM_KIND
  &END QMMM
  &SUBSYS
    &CELL
@IF ${QS_METHOD} == DFT
      PERIODIC XYZ
@ENDIF
      ABC [angstrom] 42.4780006  48.2669983  60.3699989
    &END CELL
    &TOPOLOGY
      COORDINATE CRD
      COORD_FILE_NAME ../sample_rst/Lysozyme_small.rst
      CONNECTIVITY AMBER
      CONN_FILE_NAME ../sample_top/Lysozyme_small.top
      &DUMP_PDB off
      &END DUMP_PDB
    &END TOPOLOGY
#USER: BASIS_SET - double zeta in valence + polarization (DZVP)
#USER: try to use basis functions consistent with POTENTIAL (e.g. GTH) and XC_FUNCTIONAL (e.g. BLYP) for all
    &KIND H
      BASIS_SET DZVP-GTH-${DFT_FUNCTIONAL}
      POTENTIAL GTH-${DFT_FUNCTIONAL}-q1
    &END KIND
    &KIND C
      BASIS_SET DZVP-GTH-${DFT_FUNCTIONAL}
      POTENTIAL GTH-${DFT_FUNCTIONAL}-q4
    &END KIND
    &KIND N
      BASIS_SET DZVP-GTH-${DFT_FUNCTIONAL}
      POTENTIAL GTH-${DFT_FUNCTIONAL}-q5
    &END KIND
    &KIND O
      BASIS_SET DZVP-GTH-${DFT_FUNCTIONAL}
      POTENTIAL GTH-${DFT_FUNCTIONAL}-q6
    &END KIND
    &KIND Cl
      BASIS_SET DZVP-GTH-${DFT_FUNCTIONAL}
      POTENTIAL GTH-${DFT_FUNCTIONAL}-q7
    &END KIND
    &KIND S
      BASIS_SET DZVP-GTH-${DFT_FUNCTIONAL}
      POTENTIAL GTH-${DFT_FUNCTIONAL}-q6
    &END KIND

    &KIND OX
      ELEMENT O
    &END KIND
    &KIND HX
      ELEMENT H
    &END KIND

  &END SUBSYS
&END FORCE_EVAL


&GLOBAL
  PROJECT Lysozyme_small_NVT
# RUN_TYPE GEO_OPT should also work (no guarantees about convergence without a well defined total energy, though)
  RUN_TYPE MD
  PRINT_LEVEL low
&END GLOBAL

&MOTION
   &PRINT
      &FORCES off
         COMMON_ITERATION_LEVELS 0
         FORMAT XYZ
         &EACH
            MD 1
         &END EACH
      &END FORCES
      &FORCE_MIXING_LABELS off
         COMMON_ITERATION_LEVELS 0
         FORMAT XYZ
         &EACH
            MD 1
         &END EACH
      &END FORCE_MIXING_LABELS
      &RESTART off
      &END RESTART
      &TRAJECTORY off
         COMMON_ITERATION_LEVELS 0
         FORMAT XYZ
         &EACH
            MD 1
         &END EACH
      &END TRAJECTORY
   &END PRINT

  &MD
    ENSEMBLE NVT
    STEPS 2
    TIMESTEP 0.5
    TEMPERATURE 300.0
    &THERMOSTAT
       REGION MASSIVE
       TYPE AD_LANGEVIN
       &AD_LANGEVIN
         TIMECON_NH       74.0
         TIMECON_LANGEVIN 370.0
       &END AD_LANGEVIN
    &END
  &END MD

# not active in this regtest, just to show how it to make it work for force mixing
   &GEO_OPT
      OPTIMIZER CG
      &CG
   &LINE_SEARCH
      TYPE 2PNT
      &2PNT
         LINMIN_GRAD_ONLY
      &END 2PNT
   &END LINE_SEARCH
      &END CG
   &END GEO_OPT

&END MOTION
