is use of a restraint for pmf calculations valid?
garold
g.murd... at gmail.com
Wed Mar 14 16:18:28 UTC 2012
Dear all,
I am calculating a PMF curve of a molecule solvated in water, as a
function of depth within the water slab. I am new to this type of
calculation so please forgive me if my questions seem trivial. For a
first attempt, I am wondering if the direct approach of using a
harmonic restraint on the center-of-mass of the solute molecule with
respect to a fixed point in the slab (z_s) is valid. Thus, my
restraining potential (on COLVAR 2) is just K*(z-z_s)**2 and force of
the restraint is 2*K*(z-z_s). I have set K to 1 au (hartree/bohr**2,
I have also tried 0.1 and 0.01). I can easily post-process what the
force is at any frame from the z value printed out (COM of solute
molecule). Then I average the force over the run (after a suitable
equilibration period). The use of an explicit restraint seems simpler
than not using one and instead printing out the lagrange multiplier
coming from the Shake algorithm. The PMF can be obtained as usual
from integrating the average force on the solute (negative of the
restraint force). I have some preliminary results with such
calculations and would like to know if the approach is reasonable
before going further. Please also advise me if my input file (see
below) seems reasonable. For now I using classical, polarizable
potentials.
Thank you,
Garold
---------------------------------------------------------------------------------------------------------------------------------------------------------------------
@SET BASE_NAME run
@SET ID 01
&FORCE_EVAL
METHOD FIST
&MM
&POISSON
&EWALD
EWALD_TYPE EWALD
ALPHA 0.3
GMAX 31 31 143 ! 25
O_SPLINE 6
&MULTIPOLES T
MAX_MULTIPOLE_EXPANSION DIPOLE
POL_SCF SELF_CONSISTENT
EPS_POL 1.0e-6
MAX_IPOL_ITER 100
&END MULTIPOLES
&END EWALD
&END POISSON
&FORCEFIELD
&SPLINE
EMAX_SPLINE 10000000000.0
RCUT_NB 12.00
&END SPLINE
&BEND
ATOMS H O H
THETA0 [deg] 109.500
!K [rad^2kjmol] 627.600
K [rad^2kjmol] 0.000
KIND G87
&END BEND
&BEND
ATOMS H O M
THETA0 [deg] 54.750
!K [rad^2kjmol] 418.400
K [rad^2kjmol] 0.000
KIND G87
&END BEND
&BEND
ATOMS M H O
THETA0 [deg] 109.500
K [rad^2kjmol] 0.000
KIND G87
&END BEND
&BEND
ATOMS M H H
THETA0 [deg] 109.500
K [rad^2kjmol] 0.000
KIND G87
&END BEND
&BEND
ATOMS O H H
THETA0 [deg] 109.500
K [rad^2kjmol] 0.000
KIND G87
&END BEND
&BEND
ATOMS O M H
THETA0 [deg] 109.500
K [rad^2kjmol] 0.000
KIND G87
&END BEND
&BEND
ATOMS H M H
THETA0 [deg] 109.500
K [rad^2kjmol] 0.000
KIND G87
&END BEND
&BEND
ATOMS O2 N2 O2
K 0.
THETA0 1.8
&END BEND
&BOND
ATOMS O H
K [nm^-2kjmol] 502080.0
!K [nm^-2kjmol] 0.0
R0 [nm] 0.09572
KIND G87
&END BOND
&BOND
ATOMS O M
!K [nm^-2kjmol] 753120.0
K [nm^-2kjmol] 00.0
R0 [nm] 0.01750
KIND G87
&END BOND
&BOND
ATOMS H H
K [nm^-2kjmol] 000000.0
R0 [nm] 0.01750
KIND G87
&END BOND
&BOND
ATOMS H M
K [nm^-2kjmol] 000000.0
R0 [nm] 0.01750
KIND G87
&END BOND
&BOND
ATOMS N2 O2
K 0.
R0 1.8
&END BOND
&CHARGE
ATOM O
CHARGE 0.000000
&END CHARGE
&CHARGE
ATOM H
CHARGE 0.5190
&END CHARGE
&CHARGE
ATOM M
CHARGE -1.038
&END CHARGE
&DIPOLE
ATOM M
APOL 1.444
&END DIPOLE
&CHARGE
ATOM N2
CHARGE 0.28
&END CHARGE
&CHARGE
ATOM O2
CHARGE -0.14
&END CHARGE
&NONBONDED
&LENNARD-JONES
ATOMS O O
EPSILON [kcalmol] 0.1825
SIGMA [nm] 3.234E-01
&END LENNARD-JONES
&LENNARD-JONES
ATOMS O M
EPSILON 0.0
SIGMA 0.1
&END LENNARD-JONES
&LENNARD-JONES
ATOMS O H
EPSILON 0.0
SIGMA 0.1
&END LENNARD-JONES
&LENNARD-JONES
ATOMS H H
EPSILON 0.0
SIGMA 0.1
&END LENNARD-JONES
&LENNARD-JONES
ATOMS H M
EPSILON 0.0
SIGMA 0.1
&END LENNARD-JONES
&LENNARD-JONES
ATOMS M M
EPSILON 0.0
SIGMA 0.1
&END LENNARD-JONES
&LENNARD-JONES
ATOMS N2 N2
EPSILON 78.198
SIGMA 3.166
&END LENNARD-JONES
&LENNARD-JONES
ATOMS N2 O2
EPSILON 0.0
SIGMA 3.6705
&END LENNARD-JONES
&LENNARD-JONES
ATOMS O2 O2
EPSILON 0.0
SIGMA 3.30523
&END LENNARD-JONES
&WILLIAMS
ATOMS O N2
A [kcalmol] 56844.355799
B [angstrom^-1] 3.674000
C [angstrom^6kcalmol] 303.547232
&END WILLIAMS
&WILLIAMS
ATOMS O O2
A [kcalmol] 60718.546747
B [angstrom^-1] 3.798000
C [angstrom^6kcalmol] 243.535105
&END WILLIAMS
&WILLIAMS
ATOMS H N2
A [kcalmol] 9600.936470
B [angstrom^-1] 3.420000
C [angstrom^6kcalmol] 96.538689
&END WILLIAMS
&WILLIAMS
ATOMS H O2
A [kcalmol] 10255.282195
B [angstrom^-1] 3.544000
C [angstrom^6kcalmol] 77.452723
&END WILLIAMS
&WILLIAMS
ATOMS M N2
A [kcalmol] 0.0
B [angstrom^-1] 0.0
C [angstrom^6kcalmol] 0.0
&END WILLIAMS
&WILLIAMS
ATOMS M O2
A [kcalmol] 0.0
B [angstrom^-1] 0.0
C [angstrom^6kcalmol] 0.0
&END WILLIAMS
&END NONBONDED
&END FORCEFIELD
&END MM
&SUBSYS
&CELL
ABC 15.0 15.0 71.44
&END CELL
&COORD
O 3.1662481629706754E+01 3.8212607304406419E+01
4.6145495800727794E+01 H2O 1
H 3.2370761924609546E+01 3.7649888127503075E+01
4.6458414083339115E+01 H2O 1
H 3.0947883930321868E+01 3.8065366918359565E+01
4.6765094135574984E+01 H2O 1
M 3.1661322487315456E+01 3.8082341112356445E+01
4.6316597642464885E+01 H2O 1
O -2.2427516707894708E+00 1.8141648260858169E+01
3.2873224002628270E+01 H2O 2
H -2.4001057532262369E+00 1.7199422463306000E+01
3.2812545019755355E+01 H2O 2
H -3.1165745469166026E+00 1.8532055176187448E+01
3.2857484834454446E+01 H2O 2
M -2.4319561026075309E+00 1.8040398347145139E+01
3.2859202497658337E+01 H2O 2
O 1.7436691005838029E-01 3.1843195435862381E+01
4.8711658036883598E+01 H2O 3
H 7.2474912269787284E-01 3.1157291690181463E+01
4.8333706755015449E+01 H2O 3
H -5.8818366249847054E-01 3.1883779640067164E+01
4.8134508795822704E+01 H2O 3
M 1.3543741919403809E-01 3.1724789650926954E+01
4.8536412410738507E+01 H2O 3
.
.
.
[200 molecule coordinates removed]
.
.
.
.
O 2.1548881835647431E+00 2.4479166362552618E+01
2.0991429221689963E+01 H2O 213
H 2.5225875195388769E+00 2.4517118796527051E+01
2.1874372317746786E+01 H2O 213
H 2.6206746315755836E+00 2.3757683165022637E+01
2.0568656611027709E+01 H2O 213
M 2.3078194558140108E+00 2.4353749431273627E+01
2.1075863124625279E+01 H2O 213
O -7.0956688752052228E+00 4.2163194469506976E+01
4.7423581212118272E+01 H2O 214
H -7.9426919155240148E+00 4.2534965312023658E+01
4.7669691166444274E+01 H2O 214
H -6.4551414437399695E+00 4.2810951598408472E+01
4.7717463860964060E+01 H2O 214
M -7.1335575072402806E+00 4.2350261499033550E+01
4.7522661193318939E+01 H2O 214
O 5.5548794114490008E+00 2.6354706038205300E+01
5.0074069011426154E+01 H2O 215
H 5.4713169329547178E+00 2.6619254563897609E+01
4.9157955978260077E+01 H2O 215
H 4.6660240621338263E+00 2.6115214639965210E+01
5.0336388926715927E+01 H2O 215
M 5.3764563376279240E+00 2.6359303619176455E+01
4.9954108461137466E+01 H2O 215
N2 7.64284600146 7.50000000000
30.70968697232 NO2
O2 7.43753762385 6.40143134468
31.12711286950 NO2
O2 7.43740493577 8.59856865532
31.12704758148 NO2
&END COORD
&COLVAR
&DISTANCE
POINTS 1 2
AXIS Z
&POINT
TYPE FIX_POINT
XYZ 15.1065589017861 14.2580626814867 67.7442647138845
&END POINT
&POINT
TYPE GEO_CENTER
ATOMS 1 5 9 13 17 21 25 29 33 37
ATOMS 41 45 49 53 57 61 65 69 73 77
ATOMS 81 85 89 93 97 101 105 109 113 117
ATOMS 121 125 129 133 137 141 145 149 153 157
ATOMS 161 165 169 173 177 181 185 189 193 197
ATOMS 201 205 209 213 217 221 225 229 233 237
ATOMS 241 245 249 253 257 261 265 269 273 277
ATOMS 281 285 289 293 297 301 305 309 313 317
ATOMS 321 325 329 333 337 341 345 349 353 357
ATOMS 361 365 369 373 377 381 385 389 393 397
ATOMS 401 405 409 413 417 421 425 429 433 437
ATOMS 441 445 449 453 457 461 465 469 473 477
ATOMS 481 485 489 493 497 501 505 509 513 517
ATOMS 521 525 529 533 537 541 545 549 553 557
ATOMS 561 565 569 573 577 581 585 589 593 597
ATOMS 601 605 609 613 617 621 625 629 633 637
ATOMS 641 645 649 653 657 661 665 669 673 677
ATOMS 681 685 689 693 697 701 705 709 713 717
ATOMS 721 725 729 733 737 741 745 749 753 757
ATOMS 761 765 769 773 777 781 785 789 793 797
ATOMS 801 805 809 813 817 821 825 829 833 837
ATOMS 841 845 849 853 857
&END POINT
&END DISTANCE
&END COLVAR
&COLVAR
&DISTANCE
POINTS 1 2
AXIS Z
&POINT
TYPE FIX_POINT
XYZ [angstrom] 0.0 0.0 31.0
&END POINT
&POINT
TYPE GEO_CENTER
ATOMS 861 862 863
WEIGHTS 0.30446162366521585741 0.34776918816739207129
0.34776918816739207129
&END POINT
&END DISTANCE
&END COLVAR
&PRINT
&CELL
&END CELL
&END PRINT
&KIND O
ELEMENT O
MASS 15.99940
&END KIND
&KIND H
ELEMENT H
MASS 1.00800
&END KIND
&KIND M
ELEMENT H
MASS 0.00000
&END KIND
&KIND N2
ELEMENT N
&END KIND
&KIND O2
ELEMENT O
&END KIND
&END SUBSYS
&PRINT
&GRID_INFORMATION
&END GRID_INFORMATION
&END PRINT
&END FORCE_EVAL
&GLOBAL
PROJECT ${BASE_NAME}-${ID}
RUN_TYPE MD
! PREFERRED_FFT_LIBRARY FFTW
PRINT_LEVEL MEDIUM
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVT
STEPS 10000
TIMESTEP 1.0
TEMPERATURE 300.0
! TEMP_TOL 50
&THERMOSTAT
TYPE NOSE
REGION MOLECULE
&NOSE
LENGTH 3
YOSHIDA 3
TIMECON 100
MTS 2
&END NOSE
&END THERMOSTAT
&PRINT ON
&ENERGY
&EACH
MD 1
&END EACH
FILENAME =${BASE_NAME}-${ID}.ener
&END ENERGY
&END PRINT
&END MD
&PRINT
&TRAJECTORY ON
&EACH
MD 1
&END EACH
FILENAME =${BASE_NAME}-${ID}.xyz
FORMAT XYZ
&END TRAJECTORY
&VELOCITIES ON
&EACH
MD 1
&END EACH
FILENAME =${BASE_NAME}-${ID}_vel.xyz
FORMAT XYZ
&END VELOCITIES
&FORCES ON
&EACH
MD 1
&END EACH
FILENAME =${BASE_NAME}-${ID}_force.xyz
FORMAT XYZ
&END FORCES
&RESTART ON
&EACH
MD 1
&END EACH
FILENAME =${BASE_NAME}-${ID}.restart
&END RESTART
&END PRINT
&CONSTRAINT
&G3X3
INTERMOLECULAR FALSE
DISTANCES 1.808845716 1.808845716 2.860459332
MOLNAME H2O
ATOMS 1 2 3
&END G3X3
&G3X3
INTERMOLECULAR FALSE
DISTANCES 2.2544431 2.2544431 4.1519869
MOLNAME NO2
ATOMS 1 2 3
&END G3X3
&VIRTUAL_SITE
INTERMOLECULAR FALSE
ATOMS 4 2 1 3
PARAMETERS 0.18348396 0.18348396
MOLNAME H2O
&END VIRTUAL_SITE
CONSTRAINT_INIT T
SHAKE_TOLERANCE 1.0E-12
ROLL_TOLERANCE 1.0E-10
#1st - COM alway set to 0.0
&COLLECTIVE
COLVAR 1
INTERMOLECULAR
TARGET 0.0
&RESTRAINT
K 1.0
&END RESTRAINT
&END COLLECTIVE
&COLLECTIVE
COLVAR 2
INTERMOLECULAR
TARGET 0.0
&RESTRAINT
K 1.0
&END RESTRAINT
&END COLLECTIVE
&END CONSTRAINT
&END MOTION
&EXT_RESTART
RESTART_FILE_NAME ./restart
RESTART_COUNTERS F
RESTART_AVERAGES F
RESTART_VEL F
&END EXT_RESTART
---------------------------------------------------------------------------------------------------------------------------------------------------------------------
More information about the CP2K-user
mailing list