Restraint in QMMM
Ruibin
ruibin... at gmail.com
Thu Jul 5 18:56:12 UTC 2012
Dear CP2K group:
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:
&FORCE_EVAL
METHOD QMMM
&DFT
BASIS_SET_FILE_NAME GTH_BASIS_SETS
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
COMMENSURATE
CUTOFF 300
&END MGRID
&QS
&END QS
&SCF
SCF_GUESS atomic
&END SCF
&XC
&XC_FUNCTIONAL pade
&END XC_FUNCTIONAL
&END XC
&END DFT
&MM
&FORCEFIELD
&BEND
ATOMS H O H
K 0.002
THETA0 1.823
&END BEND
&BOND
ATOMS O H
K 0.02
R0 1.89
&END BOND
&CHARGE
ATOM O
CHARGE -0.8476
&END CHARGE
&CHARGE
ATOM H
CHARGE 0.4238
&END CHARGE
&NONBONDED
&LENNARD-JONES
atoms O O
EPSILON 78.198
SIGMA 3.166
RCUT 11.4
&END LENNARD-JONES
&LENNARD-JONES
atoms O H
EPSILON 0.0
SIGMA 3.6705
RCUT 11.4
&END LENNARD-JONES
&LENNARD-JONES
atoms H H
EPSILON 0.0
SIGMA 3.30523
RCUT 11.4
&END LENNARD-JONES
&END NONBONDED
&END FORCEFIELD
&POISSON
&EWALD
EWALD_TYPE ewald
ALPHA .44
GMAX 21
&END EWALD
&END POISSON
&END MM
&QMMM
&CELL
ABC 8.0 8.0 8.0
&END CELL
ECOUPL GAUSS
&MM_KIND H
RADIUS 0.44
&END MM_KIND
&MM_KIND O
RADIUS 0.78
&END MM_KIND
&QM_KIND H
MM_INDEX 2 3
&END QM_KIND
&QM_KIND O
MM_INDEX 1
&END QM_KIND
#
# QM_KINDS
#
#
# MM_KINDS
#
#
&END QMMM
&SUBSYS
&CELL
ABC 24.955 24.955 24.955
&END CELL
&COORD
O 0.000000 0.000000 0.000000 H2O
H 0.000000 0.000000 1.000000 H2O
H 0.942809 0.000000 -0.333333 H2O
O -1.617979 -0.948062 -2.341650 H2O
H -2.529195 -1.296822 -2.122437 H2O
H -1.534288 -0.833088 -3.331486 H2O
O -1.447990 2.117783 1.555094 H2O
H -1.501128 2.645178 2.403050 H2O
H -2.090603 1.352766 1.597519 H2O
&END COORD
&KIND H
BASIS_SET SZV-GTH
POTENTIAL GTH-PADE-q1
&END KIND
&KIND O
BASIS_SET SZV-GTH
POTENTIAL GTH-PADE-q6
&END KIND
&TOPOLOGY
&END TOPOLOGY
&END SUBSYS
&PRINT
&FORCES
&EACH
MD 1
&END
FILENAME ./force #print force out
&END
&END
&END FORCE_EVAL
&GLOBAL
PROJECT water3
RUN_TYPE MD
&END GLOBAL
&MOTION
&CONSTRAINT
&FIXED_ATOMS
LIST 4 5 6 7 8 9 #add restraint to MM water.
&RESTRAINT
K 0.3
&END
&END FIXED_ATOMS
&END CONSTRAINT
&MD
ENSEMBLE NVT
STEPS 10
TIMESTEP 1.0
TEMPERATURE 300
&THERMOSTAT
TYPE NOSE
&NOSE #Nose Hoover thermostat
TIMECON 100
&END NOSE
&END THERMOSTAT
&END MD
&END MOTION
The problem is the energies and temperature are wrong. Below is the energy
output:
# Step Nr. Time[fs] Kin.[a.u.] Temp[K]
Pot.[a.u.] Cons Qty[a.u.] UsedTime[s]
0 0.000000 0.011400535 300.000000000
252.270686409 252.282561967 0.000000000
1 1.000000 64.715425231 1702957.551343301
162.918745117 238.987526611 2.638700000
2 2.000000 83.864213570 2206849.375982836
70.736831365 224.666805827 1.220200000
3 3.000000 34.934722677 919291.646080756
78.458120574 226.224275908 1.224400000
4 4.000000 14.032290213 369253.458444565
85.560270123 228.443716237 1.134200000
5 5.000000 28.967769113 762273.924353682
51.849454675 224.604916497 1.134400000
6 6.000000 31.469118622 828095.821048970
22.515351209 220.561280390 1.141300000
7 7.000000 18.256265809 480405.492313880
17.640008110 219.578630756 1.153600000
8 8.000000 11.471213005 301859.853983926
15.162278099 218.953890894 1.144900000
9 9.000000 13.316570694 350419.618504665
4.711198408 217.260769199 1.133700000
10 10.000000 13.377280141 352017.160530907
-4.491490564 215.534388083 1.209100000
Below is the data file containing force evaluation in the initial step (0th
step) in MD simulation:
ATOMIC FORCES in [a.u.]
# Atom Kind Element X Y Z
1 1 O 0.00260969 -0.00099565 0.00066949
2 2 H 0.00000000 0.00000000 0.00000000
3 2 H 0.00000000 0.00000000 0.00000000
4 3 O -0.00182578 -0.00164102 -0.00180926
5 4 H 0.00027072 0.00026737 -0.00038465
6 4 H 0.00002536 0.00025623 0.00009876
7 3 O -0.00157474 0.00277115 0.00209313
8 4 H 0.00003024 -0.00036715 -0.00022287
9 4 H 0.00046452 -0.00029094 -0.00044461
SUM OF ATOMIC FORCES 0.00000000 0.00000000 0.00000000
0.00000000
ATOMIC FORCES in [a.u.]
# Atom Kind Element X Y Z
1 1 H 0.00649874 0.00414139 0.07816522
2 1 H 0.07417144 0.00048362 -0.02114914
3 2 O -0.08150537 -0.00333247 -0.05532788
SUM OF ATOMIC FORCES -0.00083518 0.00129255 0.00168820
0.00228434
ATOMIC FORCES in [a.u.]
# Atom Kind Element X Y Z
1 1 O -0.07889568 -0.00432811 -0.05465838
2 2 H 0.00649874 0.00414139 0.07816522
3 2 H 0.07417144 0.00048362 -0.02114914
4 3 O -4.00472633 -4.53949313 -4.16502294
5 4 H -3.99913852 -4.53403585 -4.15640920
6 4 H -4.00072579 -4.53464884 -4.15564218
7 3 O -4.00042985 -4.53333413 -4.15089790
8 4 H -4.00149003 -4.53531132 -4.15839300
9 4 H -4.00121088 -4.53544321 -4.16032995
SUM OF ATOMIC FORCES -24.00594690 -27.21196959 -24.94433748
44.03404079
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.
As a reference, the output files without restraint are listed below:
energy output:
# Step Nr. Time[fs] Kin.[a.u.] Temp[K]
Pot.[a.u.] Cons Qty[a.u.] UsedTime[s]
0 0.000000 0.011400535 300.000000000
-17.029645853 -17.017770295 0.000000000
1 1.000000 0.016402713 431.630082719
-17.035306615 -17.018372321 2.641500000
2 2.000000 0.022246916 585.417660357
-17.041632204 -17.018773197 1.147400000
3 3.000000 0.021978322 578.349735268
-17.041243389 -17.018556970 1.146100000
4 4.000000 0.016973328 446.645559080
-17.035929015 -17.018159618 1.145500000
5 5.000000 0.012217943 321.509721676
-17.031014978 -17.017933498 1.151800000
6 6.000000 0.011191631 294.502787184
-17.030016007 -17.017906376 1.149900000
7 7.000000 0.014440145 379.985982395
-17.033494688 -17.018076745 1.148600000
8 8.000000 0.019534172 514.033031159
-17.039012417 -17.018420075 1.161700000
9 9.000000 0.021857814 575.178635154
-17.041751241 -17.018734447 1.148400000
10 10.000000 0.018116905 476.738272741
-17.037946302 -17.018569700 1.154100000
This seems more reasonable
force output:
ATOMIC FORCES in [a.u.]
# Atom Kind Element X Y Z
1 1 O 0.00260969 -0.00099565 0.00066949
2 2 H 0.00000000 0.00000000 0.00000000
3 2 H 0.00000000 0.00000000 0.00000000
4 3 O -0.00182578 -0.00164102 -0.00180926
5 4 H 0.00027072 0.00026737 -0.00038465
6 4 H 0.00002536 0.00025623 0.00009876
7 3 O -0.00157474 0.00277115 0.00209313
8 4 H 0.00003024 -0.00036715 -0.00022287
9 4 H 0.00046452 -0.00029094 -0.00044461
SUM OF ATOMIC FORCES 0.00000000 0.00000000 0.00000000
0.00000000
ATOMIC FORCES in [a.u.]
# Atom Kind Element X Y Z
1 1 H 0.00649874 0.00414139 0.07816522
2 1 H 0.07417144 0.00048362 -0.02114914
3 2 O -0.08150537 -0.00333247 -0.05532788
SUM OF ATOMIC FORCES -0.00083518 0.00129255 0.00168820
0.00228434
ATOMIC FORCES in [a.u.]
# Atom Kind Element X Y Z
1 1 O -0.07889568 -0.00432811 -0.05465838
2 2 H 0.00649874 0.00414139 0.07816522
3 2 H 0.07417144 0.00048362 -0.02114914
4 3 O -0.00387886 -0.00415042 -0.00762564
5 4 H 0.00170896 0.00130686 0.00098811
6 4 H 0.00012168 0.00069388 0.00175512
7 3 O 0.00041762 0.00200859 0.00649940
8 4 H -0.00064255 0.00003140 -0.00099570
9 4 H -0.00036341 -0.00010049 -0.00293265
SUM OF ATOMIC FORCES -0.00086204 0.00008672 0.00004634
0.00086763
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.
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.
Could anyone in the group help me figuring out what happened?
Thanks in advance!
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20120705/95c69bb6/attachment.htm>
More information about the CP2K-user
mailing list