[gmx-users] Settles vs. 3 normal constraints (energy conservation problem)
Mark Abraham
mark.j.abraham at gmail.com
Tue Jun 25 09:45:29 CEST 2013
On Tue, Jun 25, 2013 at 1:34 AM, S. Alireza Bagherzadeh
<s.a.bagherzadeh.h at gmail.com> wrote:
> 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).
Hmm, that doesn't sound good. Guessing wildly, SHAKE and the v-site
are not properly implemented?
> .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 believe not, but Berk's probably the only expert on the
implementation, here...
> * 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.
I can't see anything wrong. I would suggest trying LINCS, which is
more likely to be better tested with v-sites. Please report back,
either way! :-)
Mark
More information about the gromacs.org_gmx-users
mailing list