Dear CP2K group:<div>    I have a question about applying restraint in QM/MM calculation.I want to add harmonic restraints to some atoms during MD simulation and the reference points are the initial atomic positions. I modified the test input file. The system contains 3 water molecules, the first being QM water, the second and third being MM water. I tried to add harmonic restraint on the MM water molecules and run NVT simulation using Nose Hoover thermostat. Below is my input:</div><div><br></div><div><br></div><div><div>&FORCE_EVAL</div><div>  METHOD QMMM</div><div>  &DFT</div><div>    BASIS_SET_FILE_NAME GTH_BASIS_SETS</div><div>    POTENTIAL_FILE_NAME POTENTIAL</div><div>    &MGRID</div><div>      COMMENSURATE</div><div>      CUTOFF 300</div><div>    &END MGRID</div><div>    &QS</div><div>    &END QS</div><div>    &SCF</div><div>      SCF_GUESS atomic</div><div>    &END SCF</div><div>    &XC</div><div>      &XC_FUNCTIONAL pade</div><div>      &END XC_FUNCTIONAL</div><div>    &END XC</div><div>  &END DFT</div><div>  &MM</div><div>    &FORCEFIELD</div><div>      &BEND</div><div>        ATOMS H O H</div><div>        K 0.002</div><div>        THETA0 1.823</div><div>      &END BEND</div><div>      &BOND</div><div>        ATOMS O H</div><div>        K 0.02</div><div>        R0 1.89</div><div>      &END BOND</div><div>      &CHARGE</div><div>        ATOM O</div><div>        CHARGE -0.8476</div><div>      &END CHARGE</div><div>      &CHARGE</div><div>        ATOM H</div><div>        CHARGE 0.4238</div><div>      &END CHARGE</div><div>      &NONBONDED</div><div>        &LENNARD-JONES</div><div>          atoms O O</div><div>          EPSILON 78.198</div><div>          SIGMA 3.166</div><div>          RCUT 11.4</div><div>        &END LENNARD-JONES</div><div>        &LENNARD-JONES</div><div>          atoms O H</div><div>          EPSILON 0.0</div><div>          SIGMA 3.6705</div><div>          RCUT 11.4</div><div>        &END LENNARD-JONES</div><div>        &LENNARD-JONES</div><div>          atoms H H</div><div>          EPSILON 0.0</div><div>          SIGMA 3.30523</div><div>          RCUT 11.4</div><div>        &END LENNARD-JONES</div><div>      &END NONBONDED</div><div>    &END FORCEFIELD</div><div>    &POISSON</div><div>      &EWALD</div><div>        EWALD_TYPE ewald</div><div>        ALPHA .44</div><div>        GMAX 21</div><div>      &END EWALD</div><div>    &END POISSON</div><div>  &END MM</div><div>  &QMMM</div><div>    &CELL</div><div>      ABC 8.0 8.0 8.0</div><div>    &END CELL</div><div>    ECOUPL GAUSS</div><div>    &MM_KIND H</div><div>      RADIUS 0.44</div><div>    &END MM_KIND</div><div>    &MM_KIND O</div><div>      RADIUS 0.78</div><div>    &END MM_KIND</div><div>    &QM_KIND H</div><div>      MM_INDEX 2 3</div><div>    &END QM_KIND</div><div>    &QM_KIND O</div><div>      MM_INDEX 1</div><div>    &END QM_KIND</div><div>  # </div><div>  # QM_KINDS</div><div>  # </div><div>  # </div><div>  # MM_KINDS</div><div>  # </div><div>  # </div><div>  &END QMMM</div><div>  &SUBSYS</div><div>    &CELL</div><div>      ABC 24.955 24.955 24.955</div><div>    &END CELL</div><div>    &COORD</div><div>   O     0.000000     0.000000     0.000000    H2O</div><div>   H     0.000000     0.000000     1.000000    H2O</div><div>   H     0.942809     0.000000    -0.333333    H2O</div><div>   O    -1.617979    -0.948062    -2.341650    H2O</div><div>   H    -2.529195    -1.296822    -2.122437    H2O</div><div>   H    -1.534288    -0.833088    -3.331486    H2O</div><div>   O    -1.447990     2.117783     1.555094    H2O</div><div>   H    -1.501128     2.645178     2.403050    H2O</div><div>   H    -2.090603     1.352766     1.597519    H2O</div><div>    &END COORD</div><div>    &KIND H</div><div>      BASIS_SET SZV-GTH</div><div>      POTENTIAL GTH-PADE-q1</div><div>    &END KIND</div><div>    &KIND O</div><div>      BASIS_SET SZV-GTH</div><div>      POTENTIAL GTH-PADE-q6</div><div>    &END KIND</div><div>    &TOPOLOGY</div><div>    &END TOPOLOGY</div><div>  &END SUBSYS</div><div>  &PRINT</div><div>   &FORCES</div><div>     &EACH</div><div>       MD 1</div><div>     &END</div><div>     FILENAME ./force           #print force out</div><div>   &END</div><div>  &END</div><div>&END FORCE_EVAL</div><div>&GLOBAL</div><div><br></div><div>  PROJECT water3</div><div>  RUN_TYPE MD</div><div>&END GLOBAL</div><div>&MOTION</div><div>  &CONSTRAINT</div><div>    &FIXED_ATOMS</div><div>     LIST 4 5 6 7 8 9                 #add restraint to MM water.</div><div>     &RESTRAINT</div><div>      K 0.3</div><div>     &END</div><div>    &END FIXED_ATOMS</div><div>  &END CONSTRAINT     </div><div>     </div><div>  &MD</div><div>    ENSEMBLE NVT</div><div>    STEPS 10</div><div>    TIMESTEP 1.0</div><div>    TEMPERATURE 300</div><div>    &THERMOSTAT</div><div>      TYPE NOSE</div><div>      &NOSE                           #Nose Hoover thermostat</div><div>        TIMECON 100</div><div>      &END NOSE</div><div>    &END THERMOSTAT</div><div>  &END MD</div><div>&END MOTION</div></div><div><br></div><div><br></div><div>The problem is the energies and temperature are wrong. Below is the energy output:</div><div><br></div><div><div>#     Step Nr.          Time[fs]        Kin.[a.u.]          Temp[K]            Pot.[a.u.]        Cons Qty[a.u.]        UsedTime[s]</div><div>         0            0.000000         0.011400535       300.000000000       252.270686409       252.282561967         0.000000000</div><div>         1            1.000000        64.715425231   1702957.551343301       162.918745117       238.987526611         2.638700000</div><div>         2            2.000000        83.864213570   2206849.375982836        70.736831365       224.666805827         1.220200000</div><div>         3            3.000000        34.934722677    919291.646080756        78.458120574       226.224275908         1.224400000</div><div>         4            4.000000        14.032290213    369253.458444565        85.560270123       228.443716237         1.134200000</div><div>         5            5.000000        28.967769113    762273.924353682        51.849454675       224.604916497         1.134400000</div><div>         6            6.000000        31.469118622    828095.821048970        22.515351209       220.561280390         1.141300000</div><div>         7            7.000000        18.256265809    480405.492313880        17.640008110       219.578630756         1.153600000</div><div>         8            8.000000        11.471213005    301859.853983926        15.162278099       218.953890894         1.144900000</div><div>         9            9.000000        13.316570694    350419.618504665         4.711198408       217.260769199         1.133700000</div><div>        10           10.000000        13.377280141    352017.160530907        -4.491490564       215.534388083         1.209100000</div></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div>Below is the data file containing force evaluation in the initial step (0th step) in MD simulation:</div><div><br></div><div><div>ATOMIC FORCES in [a.u.]</div><div><br></div><div> # Atom   Kind   Element          X              Y              Z</div><div>      1      1      O           0.00260969    -0.00099565     0.00066949</div><div>      2      2      H           0.00000000     0.00000000     0.00000000</div><div>      3      2      H           0.00000000     0.00000000     0.00000000</div><div>      4      3      O          -0.00182578    -0.00164102    -0.00180926</div><div>      5      4      H           0.00027072     0.00026737    -0.00038465</div><div>      6      4      H           0.00002536     0.00025623     0.00009876</div><div>      7      3      O          -0.00157474     0.00277115     0.00209313</div><div>      8      4      H           0.00003024    -0.00036715    -0.00022287</div><div>      9      4      H           0.00046452    -0.00029094    -0.00044461</div><div> SUM OF ATOMIC FORCES           0.00000000     0.00000000     0.00000000     0.00000000</div><div><br></div><div> ATOMIC FORCES in [a.u.]</div><div><br></div><div> # Atom   Kind   Element          X              Y              Z</div><div>      1      1      H           0.00649874     0.00414139     0.07816522</div><div>      2      1      H           0.07417144     0.00048362    -0.02114914</div><div>      3      2      O          -0.08150537    -0.00333247    -0.05532788</div><div> SUM OF ATOMIC FORCES          -0.00083518     0.00129255     0.00168820     0.00228434</div><div><br></div><div> ATOMIC FORCES in [a.u.]</div><div><br></div><div> # Atom   Kind   Element          X              Y              Z</div><div>      1      1      O          -0.07889568    -0.00432811    -0.05465838</div><div>      2      2      H           0.00649874     0.00414139     0.07816522</div><div>      3      2      H           0.07417144     0.00048362    -0.02114914</div><div>      4      3      O          -4.00472633    -4.53949313    -4.16502294</div><div>      5      4      H          -3.99913852    -4.53403585    -4.15640920</div><div>      6      4      H          -4.00072579    -4.53464884    -4.15564218</div><div>      7      3      O          -4.00042985    -4.53333413    -4.15089790</div><div>      8      4      H          -4.00149003    -4.53531132    -4.15839300</div><div>      9      4      H          -4.00121088    -4.53544321    -4.16032995</div><div> SUM OF ATOMIC FORCES         -24.00594690   -27.21196959   -24.94433748    44.03404079</div></div><div><br></div><div>Looking at the last force evaluation (QMMM total force evaluation), obviously, the energy error is caused by the harmonic restraint, because it puts too much forces on the MM atoms and causes huge acceleration. But remember it is the initial step and the positions should be used as equilibrium position in harmonic restraint, thus the forces should approach zero.</div><div><br></div><div> </div><div>As a reference, the output files without restraint are listed below:</div><div>energy output:</div><div><br></div><div><div>#     Step Nr.          Time[fs]        Kin.[a.u.]          Temp[K]            Pot.[a.u.]        Cons Qty[a.u.]        UsedTime[s]</div><div>         0            0.000000         0.011400535       300.000000000       -17.029645853       -17.017770295         0.000000000</div><div>         1            1.000000         0.016402713       431.630082719       -17.035306615       -17.018372321         2.641500000</div><div>         2            2.000000         0.022246916       585.417660357       -17.041632204       -17.018773197         1.147400000</div><div>         3            3.000000         0.021978322       578.349735268       -17.041243389       -17.018556970         1.146100000</div><div>         4            4.000000         0.016973328       446.645559080       -17.035929015       -17.018159618         1.145500000</div><div>         5            5.000000         0.012217943       321.509721676       -17.031014978       -17.017933498         1.151800000</div><div>         6            6.000000         0.011191631       294.502787184       -17.030016007       -17.017906376         1.149900000</div><div>         7            7.000000         0.014440145       379.985982395       -17.033494688       -17.018076745         1.148600000</div><div>         8            8.000000         0.019534172       514.033031159       -17.039012417       -17.018420075         1.161700000</div><div>         9            9.000000         0.021857814       575.178635154       -17.041751241       -17.018734447         1.148400000</div><div>        10           10.000000         0.018116905       476.738272741       -17.037946302       -17.018569700         1.154100000</div><div><br></div><div><br></div><div>This seems more reasonable</div><div><br></div><div><br></div><div><br></div><div>force output:</div><div> ATOMIC FORCES in [a.u.]</div><div><br></div><div> # Atom   Kind   Element          X              Y              Z</div><div>      1      1      O           0.00260969    -0.00099565     0.00066949</div><div>      2      2      H           0.00000000     0.00000000     0.00000000</div><div>      3      2      H           0.00000000     0.00000000     0.00000000</div><div>      4      3      O          -0.00182578    -0.00164102    -0.00180926</div><div>      5      4      H           0.00027072     0.00026737    -0.00038465</div><div>      6      4      H           0.00002536     0.00025623     0.00009876</div><div>      7      3      O          -0.00157474     0.00277115     0.00209313</div><div>      8      4      H           0.00003024    -0.00036715    -0.00022287</div><div>      9      4      H           0.00046452    -0.00029094    -0.00044461</div><div> SUM OF ATOMIC FORCES           0.00000000     0.00000000     0.00000000     0.00000000</div><div><br></div><div> ATOMIC FORCES in [a.u.]</div><div><br></div><div> # Atom   Kind   Element          X              Y              Z</div><div>      1      1      H           0.00649874     0.00414139     0.07816522</div><div>      2      1      H           0.07417144     0.00048362    -0.02114914</div><div>      3      2      O          -0.08150537    -0.00333247    -0.05532788</div><div> SUM OF ATOMIC FORCES          -0.00083518     0.00129255     0.00168820     0.00228434</div><div><br></div><div> ATOMIC FORCES in [a.u.]</div><div><br></div><div> # Atom   Kind   Element          X              Y              Z</div><div>      1      1      O          -0.07889568    -0.00432811    -0.05465838</div><div>      2      2      H           0.00649874     0.00414139     0.07816522</div><div>      3      2      H           0.07417144     0.00048362    -0.02114914</div><div>      4      3      O          -0.00387886    -0.00415042    -0.00762564</div><div>      5      4      H           0.00170896     0.00130686     0.00098811</div><div>      6      4      H           0.00012168     0.00069388     0.00175512</div><div>      7      3      O           0.00041762     0.00200859     0.00649940</div><div>      8      4      H          -0.00064255     0.00003140    -0.00099570</div><div>      9      4      H          -0.00036341    -0.00010049    -0.00293265</div><div> SUM OF ATOMIC FORCES          -0.00086204     0.00008672     0.00004634     0.00086763</div><div>                                                                                                                                  </div></div><div><br></div><div>The key point is that the force in the later test approaches zero. This is what should be expected, because no external forces are added.</div><div><br></div><div>In addition, pure MM simulation (FIST) for the same system with restraint added works fine. Also in QMMM simulation with the MM water fixed (frozen), everything works fine.</div><div><br></div><div>Could anyone in the group help me figuring out what happened? </div><div>Thanks in advance!</div><div><br></div><div><br></div>