[gmx-users] Settles vs. 3 normal constraints (energy conservation problem)
Broadbent, Richard
richard.broadbent09 at imperial.ac.uk
Tue Jun 25 10:07:08 CEST 2013
For my system reducing shake-tol greatly improves the energy conservation
generally 1.0e-7 is the largest I would use. However if you want very good
energy conservation 1.0e-9 or lower might be needed.
This effect might only be for my system but I think it might help here too
Richard
On 25/06/2013 08:45, "Mark Abraham" <mark.j.abraham at gmail.com> wrote:
>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/force
>>field.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
>--
>gmx-users mailing list gmx-users at gromacs.org
>http://lists.gromacs.org/mailman/listinfo/gmx-users
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>* Please don't post (un)subscribe requests to the list. Use the
>www interface or send it to gmx-users-request at gromacs.org.
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
More information about the gromacs.org_gmx-users
mailing list