[gmx-users] Re: Settles vs. 3 normal constraints (energy conservation problem)
S. Alireza Bagherzadeh
s.a.bagherzadeh.h at gmail.com
Wed Jun 26 21:07:53 CEST 2013
Dear Mark and Richard,
Thank you so much for your help.
I used shake tolerance of 1e-10 and it solved the issue.
I am getting very good energy conservation.
I would have never guessed that shake tolerance has such a huge impact!
I also experimented LINC with the following options:
constraint-algorithm = lincs
continuation = no ; DO apply constraints to the start
configuration
lincs-order = 6
lincs-iter = 2
lincs-warnangle = 30
I was able to stop the temperature from decreasing however the there was
again substantial total energy drift.
Many thanks,
Alireza
> Message: 5
> Date: Tue, 25 Jun 2013 13:55:08 +0200
> From: Mark Abraham <mark.j.abraham at gmail.com>
> Subject: Re: [gmx-users] Settles vs. 3 normal constraints (energy
> conservation problem)
> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> Message-ID:
> <
> CAMNuMARO6HV6+u+p6pVmCkbfQpiHb3iwTt_1HdAjErm_HYEFoA at mail.gmail.com>
> Content-Type: text/plain; charset=ISO-8859-1
>
> 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
>
>
> ------------------------------
>
> --
> 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!
>
> End of gmx-users Digest, Vol 110, Issue 128
> *******************************************
>
More information about the gromacs.org_gmx-users
mailing list