[gmx-users] Settles vs. 3 normal constraints (energy conservation problem)
S. Alireza Bagherzadeh
s.a.bagherzadeh.h at gmail.com
Tue Jun 25 01:34:02 CEST 2013
Dear All,
I have a box of 3073 tip4p water molecules. I do a 250ps nvt, then 250 ps
npt and finally a 1 ns nve (production run).
I used the opls forcefield and I copied the tip4p.itp to my working
directory (in order to be able to make changes).
In one case I used the [ settles ] directive to constraint water molecules.
.top file:
; Topology for 3087 TIPT4P waters
#include
"/global/software/gromacs/gromacs-4.5.5/share/gromacs/top/oplsaa.ff/forcefield.itp"
;SOL
;-----------------------------------------------
; water topology - liquid phase
#include "./tip4p.itp"
;------------------------------------------------
[ system ]
; Name
A box of 216 tip4p for protocol testing
[ molecules ]
; Compound #mols
SOL 3073
.itp file:
; Note the strange order of atoms to make it faster in gromacs.
;
[ moleculetype ]
; molname nrexcl
SOL 2
[ atoms ]
; id at type res nr residu name at name cg nr charge
1 opls_113 1 SOL OW 1 0.0
2 opls_114 1 SOL HW1 1 0.52
3 opls_114 1 SOL HW2 1 0.52
4 opls_115 1 SOL MW 1 -1.04
#ifndef FLEXIBLE
[ settles ]
; OW funct doh dhh
1 1 0.09572 0.15139
#else
[ bonds ]
; i j funct length force.c.
1 2 1 0.09572 502416.0 0.09572 502416.0
1 3 1 0.09572 502416.0 0.09572 502416.0
[ angles ]
; i j k funct angle force.c.
2 1 3 1 104.52 628.02 104.52 628.02
#endif
[ exclusions ]
1 2 3 4
2 1 3 4
3 1 2 4
4 1 2 3
; The position of the virtual site is computed as follows:
;
; O
;
; D
;
; H H
;
; const = distance (OD) / [ cos (angle(DOH)) * distance (OH) ]
; 0.015 nm / [ cos (52.26 deg) * 0.09572 nm ]
; Vsite pos x4 = x1 + a*(x2-x1) + b*(x3-x1)
[ virtual_sites3 ]
; Vsite from funct a b
4 1 2 3 1 0.128012065 0.128012065
.mdp file:
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Run control
integrator = md ; Leap-frog algorithm
tinit = 0 ; starting time [ps]
dt = 0.001 ; time step [ps]
nsteps = 250000 ; number of steps
nstcomm = 100 ; frequency for center of mass motion
removal [steps]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Output control
nstxout = 0 ; frequency to write coordinates to
output trajectory file [steps]
nstvout = 0 ; frequency to write velocities to
output trajectory file [steps]
nstfout = 0 ; frequency to write forces to output
trajectory file [steps]
nstlog = 500 ; frequency to write energies to log
file [steps]
nstenergy = 500 ; frequency to write energies to energy
file [steps]
nstxtcout = 000 ; frequency to write coordinates to xtc
trajectory [steps]
xtc-precision = 1000 ; precision to write to xtc trajectory
[real]
xtc_grps = SOL
energygrps = SOL
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Neighborsearching and short-range nonbonded interactions
nstlist = 10 ; frequency to update the neighbor list
[steps]
ns_type = grid ; (grid / simple) search for
neighboring list
pbc = xyz ; priodic boundary conditions (xyz / no
/ xy)
rlist = 1.7 ; cut-off distance for the short-range
neighbor list [nm]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Electrostatics
coulombtype = PME-Switch
rcoulomb_switch = 1.3 ; where to switch the Coulomb potential
[nm]
rcoulomb = 1.5 ; distance for the Coulomb cut-off [nm]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; van der Waals
vdw-type = shift ; (the LJ normal out at rvdw_switch to
reach zero at rvdw)
rvdw-switch = 1.3 ; where to strat switching the LJ
potential [nm]
rvdw = 1.5 ; cut-off distance for vdw potenrial
[nm]
DispCorr = EnerPres ; (Apply long range dispersion
corrections for Energy and Pressure)
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; EWALD/PME/PPPM parameters
fourierspacing = 0.12 ; maximum grid spacing for the FFT when
using PPPM or PME [nm]
pme_order = 6 ; interpolation order for PME
ewald_rtol = 1e-06 ; relative strength of the
Ewald-shifted direct potential at rcoulomb
epsilon_surface = 0 ; dipole correction to the Ewald
summation in 3d
optimize_fft = no ; optimal FFT plan for the grid at
startup (yes / no)
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Temperature coupling
tcoupl = no
;tcoupl = nose-hoover
nh-chain-length = 1
tc_grps = system
tau_t = 0.2 ; time constatn for coupling [ps], 1
for each group
ref_t = 300 ; refernce temperature [K], 1 for each
group
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Pressure coupling
Pcoupl = no
;Pcoupl = Parrinello-Rahman
tau_p = 0.5 ; time constant for coupling
compressibility = 4.5e-05
ref_p = 1.0 ; reference pressure for coupling [bar]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Velocity generation
gen_vel = no ; generating velocities according to
maxwell distribution?
gen_temp = 300 ; [K]
gen_seed = -1 ; initialize random generator based on
the process ID number [integer]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Bonds
constraints = all-angles
constraint-algorithm = shake
continuation = no ; DO apply constraints to the start
configuration
shake-tol = 1e-5
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
This case I obtain a very good energy conservation.
However, in another case where I change settles to 3 normal constraints,
there is a rather substantial energy drift (~2% per nano second) and also
the T does not remain constant and it decreases (~8 K per ns).
.itp file:
; Note the strange order of atoms to make it faster in gromacs.
;
[ moleculetype ]
; molname nrexcl
SOL 2
[ atoms ]
; id at type res nr residu name at name cg nr charge
1 opls_113 1 SOL OW 1 0.0
2 opls_114 1 SOL HW1 1 0.52
3 opls_114 1 SOL HW2 1 0.52
4 opls_115 1 SOL MW 1 -1.04
#ifndef FLEXIBLE
[ constraints ]
; ai aj funct b0
1 2 1 0.09572
1 3 1 0.09572
2 3 1 0.15139
#else
[ bonds ]
; i j funct length force.c.
1 2 1 0.09572 502416.0 0.09572 502416.0
1 3 1 0.09572 502416.0 0.09572 502416.0
[ angles ]
; i j k funct angle force.c.
2 1 3 1 104.52 628.02 104.52 628.02
#endif
[ exclusions ]
1 2 3 4
2 1 3 4
3 1 2 4
4 1 2 3
; The position of the virtual site is computed as follows:
;
; O
;
; D
;
; H H
;
; const = distance (OD) / [ cos (angle(DOH)) * distance (OH) ]
; 0.015 nm / [ cos (52.26 deg) * 0.09572 nm ]
; Vsite pos x4 = x1 + a*(x2-x1) + b*(x3-x1)
[ virtual_sites3 ]
; Vsite from funct a b
4 1 2 3 1 0.128012065 0.128012065
The .mdp and .top files are the same.
I was wondering if there is any difference between [settles] and 3 normal [
constraints]?
* I am asking this question because in my real simulation, the above is a
simplified version where I can check for energy conservation, I need to use
water in two different molecule types and if I use [ settles ] for both of
them I get this error:
Fatal error:
The [molecules] section of your topology specifies more than one block of
a [moleculetype] with a [settles] block. Only one such is allowed. If you
are trying to partition your solvent into different *groups* (e.g. for
freezing, T-coupling, etc.) then you are using the wrong approach. Index
files specify groups. Otherwise, you may wish to change the least-used
block of molecules with SETTLE constraints into 3 normal constraints.
I appreciate being advised on what I might be doing wrong.
--
*S. Alireza Bagherzadeh, M.A.Sc. <https://circle.ubc.ca/handle/2429/26269>
*
* *
*PhD Candidate
*
* *
*Dept. of Chem. & Bio. Eng. <http://www.chbe.ubc.ca/>*
* *
*University of BC <http://www.ubc.ca/>
*
More information about the gromacs.org_gmx-users
mailing list