[gmx-users] Settles vs. 3 normal constraints (energy conservation problem)
Mark Abraham
mark.j.abraham at gmail.com
Tue Jun 25 13:55:08 CEST 2013
Sure, worth trying.
Mark
On Tue, Jun 25, 2013 at 10:07 AM, Broadbent, Richard
<richard.broadbent09 at imperial.ac.uk> wrote:
> 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
>
> --
> 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