[gmx-users] Gromacs and Lammps: energy conservation very different
Michael Shirts
mrshirts at gmail.com
Tue May 6 14:00:04 CEST 2014
No worries! Its hard to validate these things.
This does confirm the previous suggestions that that cutoffs = bad for
electrostatics.
On Tue, May 6, 2014 at 7:20 AM, James <jamesresearching at gmail.com> wrote:
> Dear Michael and Mark,
>
> As an update, it seems there was a units error in my Lammps input file. Now
> both Lammps and Gromacs have catastrophic energy drift.
>
> Apologies for taking your time with this unconventional "validation" of
> Gromacs against Lammps. Thank you again for the discussion.
>
> With best regards,
>
> James
>
>
> On 6 May 2014 18:42, James <jamesresearching at gmail.com> wrote:
>
> > Dear Michael and Mark,
> >
> > Thanks a lot for your replies.
> > The final 100ps of the 500ps equilibration for Gromacs and Lammps yield
> > the following average energies (all kJ/mol):
> >
> > Total energy per molecule: Gromacs = -37.1, Lammps = -36.8
> > PE per molecule: Gromacs = -44.5, Lammps = -44.3
> > KE per molecule: Gromacs = 7.5, Lammps = 7.5
> >
> > So they match remarkably well (within 1 standard deviation of the
> > fluctuations).
> >
> > The drift rate for Lammps is 20 kJ/mol/ns/molecule. That for Gromacs is
> > exponential, taking about 50ps to raise the initial 20 kJ/mol/molecule!
> >
> > I have tried with the default LINCS and also shake (the latter with the
> > same tolerance as Lammps), but this made no difference, unfortunately.
> >
> > I will ask the Lammps mailing list to see if anything goes on behind the
> > scenes, although my initial intuition is that 20 kJ/mol/ns/molecule is
> > fairly reasonable (ie, it's poor, but it's not exponential). I will let
> you
> > know the response.
> >
> > I'm interested in this, not because I want to simulate water with a
> simple
> > cut-off per se, but because I'm doing other calculations utilising
> > short-term auto-correlations which seem to work with Lammps but give an
> > over-estimate in Gromacs, and I'm wondering if it's connected to this. At
> > least if I can rule out a trivial error in my input files, there can be
> > confidence to look deeper.
> >
> > Thanks again for your help.
> > With best regards,
> >
> > James
> >
> >
> >
> >
> > On 6 May 2014 15:09, Mark Abraham <mark.j.abraham at gmail.com> wrote:
> >
> >> On May 6, 2014 4:16 AM, "Michael Shirts" <mrshirts at gmail.com> wrote:
> >> >
> >> > What are the energies with the same initial configuration? Perhaps
> >> LAMMPS
> >> > is doing something under the hood like shifting the cutoff which
> Gromacs
> >> > won't do unless you ask. Knowing the step=0 energy will help identify
> if
> >> it
> >> > is a difference in dynamics or in the energy function.
> >>
> >> Using the potential-shift modifier is indeed a good idea, but it should
> >> only affect the noisiness, not the drift.
> >>
> >> > Simulating water with a cutoff is generally a bad idea (LJ = maybe OK,
> >> > water = bad). It is interesting that lammps fails less badly, but this
> >> is
> >> > guaranteed to fail in either program. Lammps working almost suggests
> >> that
> >> > it's doing something behind the scenes.
> >>
> >> Also, using the LINCS settings intended for conservation might help
> there.
> >>
> >> The OP reporting the magnitude of the drift would be nice ;-) Any drift
> >> looks horrific on the wrong scale!
> >>
> >> Mark
> >>
> >> > On Mon, May 5, 2014 at 8:21 PM, James <jamesresearching at gmail.com>
> >> wrote:
> >> >
> >> > > Dear all,
> >> > >
> >> > > (resending due to hitting 50kb limit on mailing messages)
> >> > >
> >> > > Considering a simple SPC/E water simulation, I am getting
> catastrohpic
> >> > > energy drift in Gromacs but only minor drift in Lammps for the same
> >> system
> >> > > (see attached graph).
> >> > >
> >> > > I guess I'm making an embarrassingly trivial mistake somehow with
> >> Gromacs,
> >> > > but if anyone might be able to point out where I'm going wrong it
> >> would
> >> be
> >> > > very much appreciated!
> >> > >
> >> > > Details are below, and the md.log output can be found here:
> >> > > https://gist.github.com/jamesresearching/908e4b8804c96a49c61c
> >> > >
> >> > > The systems start from the same configuration, are equilibrated for
> >> 500ps
> >> > > with NVT, and then run under NVE. Both use a simple cut-off. I
> >> understand
> >> > > this will not conserve energy perfectly, but it shouldn't be as bad
> as
> >> > > Gromacs suggests.
> >> > >
> >> > > Both systems
> >> > > - Update neighbour list at every step
> >> > > - Apply COM velocity correction every 10 steps
> >> > > - Same temperature
> >> > > - Same LJ combining rule
> >> > > - Same cut-off
> >> > > - Same time-step
> >> > >
> >> > > Running Gromacs in double-precision doesn't help.
> >> > > Switching from LINCS to Shake (constraint-algorithm = SHAKE,
> >> shake-tol =
> >> > > 0.000001) in Gromacs doesn't help either (Lammps uses shake with the
> >> same
> >> > > tolerance). Example files below are with LINCS (default).
> >> > >
> >> > > The equilibration files are as follows:
> >> > >
> >> > > ___top.top___
> >> > > #include "FF.itp"
> >> > > #include "mol.itp"
> >> > >
> >> > > [ system ]
> >> > > ; Name
> >> > > SPCEWater
> >> > >
> >> > > [ molecules ]
> >> > > ; Compound #mols
> >> > > SPCE 901
> >> > > ------------------------------
> >> > >
> >> > >
> >> > > ___FF.itp___
> >> > > [ defaults ]
> >> > > ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
> >> > > 1 2 no 0.5 0.5
> >> > >
> >> > > [ atomtypes ]
> >> > > ; name bond_type mass charge ptype sigma
> >> epsilon
> >> > > OW OW 15.9994 -0.8476 A 0.3166 0.65
> >> > > HW HW 1.008 0.4238 A 0.0 0.0
> >> > > ------------------------------
> >> > >
> >> > >
> >> > > ___mol.itp___
> >> > > ; molname nrexcl
> >> > > SPCE 2
> >> > >
> >> > > [ atoms ]
> >> > > ; id at type res nr residu name at name cg nr charge
> mass
> >> > > 1 OW 1 SOL OW 1 -0.8476 15.9994
> >> > > 2 HW 1 SOL HW 2 0.4238 1.008
> >> > > 3 HW 1 SOL HW 3 0.4238 1.008
> >> > >
> >> > >
> >> > > [ settles ]
> >> > > ; OW funct doh dhh
> >> > > 1 1 0.1000 0.1632991
> >> > >
> >> > > [ exclusions ]
> >> > > 1 2 3
> >> > > 2 1 3
> >> > > 3 1 2
> >> > > --------------------------------
> >> > >
> >> > >
> >> > > ___parameters file (Equilibration)___
> >> > > cutoff-scheme = group
> >> > > constraints = all-angles
> >> > > integrator = md-vv
> >> > > dt = 0.002
> >> > > nsteps = 250000
> >> > > nstlist = 1 ; neighbour-list update period
> >> > > ns_type = grid ; restrict building of new neighbour
> lists
> >> to
> >> > > those atoms in neighbouring cells
> >> > > coulombtype = cut-off
> >> > >
> >> > > rlist = 1.2
> >> > > rcoulomb = 1.2 ; Maximum distance for short-range coulomb
> >> > > calculation
> >> > > rvdw = 1.2 ; Maximum distance for short-range vdw
> >> > > calculation
> >> > >
> >> > >
> >> > > nstxout = 0 ; period with which to write out
> co-ordinates
> >> to a
> >> > > trajectory file
> >> > > nstvout = 0 ; " velocities
> >> > > nstfout = 0 ; " forces
> >> > > nstlog = 500 ; period with which to write out energies
> to
> >> log
> >> > > file
> >> > > nstenergy = 500 ; period with which to write out energies
> to
> >> > > energy file
> >> > > energygrps = System ; the groups to write to the energy
> file
> >> > > nstxtcout = 500 ; period with which to write out positions
> >> to
> >> an
> >> > > xtc file
> >> > > xtc-precision = 10000 ; xtc precision
> >> > >
> >> > > nstcomm = 10 ; number of steps between COM removal
> >> > > comm_mode = Linear ; remove translational COM (other
> options
> >> > > include angular)
> >> > > comm_grps = System ; remove COM for the whole system
> >> > >
> >> > > Tcoupl = nose-hoover ; Thermostatting
> >> > > tc-grps = System ; can also couple groups separately to
> >> bath
> >> > > tau_t = 0.5 ; nose-hoover time-constant (ps)
> >> > > ref_t = 300.0 ; Temperature
> >> > >
> >> > > DispCorr = EnerPres
> >> > >
> >> > > Pcoupl = no
> >> > >
> >> > > gen_vel = yes
> >> > > gen_temp = 300.0
> >> > > gen_seed = 173529
> >> > > -------------------------------------
> >> > >
> >> > >
> >> > > The production calculation is then performed as:
> >> > >
> >> > > ___parameters file (production)___
> >> > > cutoff-scheme = group
> >> > > constraints = all-angles
> >> > > integrator = md-vv
> >> > > dt = 0.002
> >> > > nsteps = 250000
> >> > > nstlist = 1 ; neighbour-list update period
> >> > > ns_type = grid ; restrict building of new neighbour
> lists
> >> to
> >> > > those atoms in neighbouring cells
> >> > > coulombtype = cut-off
> >> > >
> >> > > rlist = 1.2
> >> > > rcoulomb = 1.2 ; Maximum distance for short-range coulomb
> >> > > calculation
> >> > > rvdw = 1.2 ; Maximum distance for short-range vdw
> >> > > calculation
> >> > >
> >> > >
> >> > > nstxout = 0 ; period with which to write out
> co-ordinates
> >> to a
> >> > > trajectory file
> >> > > nstvout = 0; " velocities
> >> > > nstfout = 0 ; " forces
> >> > > nstlog = 100 ; period with which to write out energies
> to
> >> log
> >> > > file
> >> > > nstenergy = 100 ; period with which to write out energies
> to
> >> > > energy file
> >> > > energygrps = System ; the groups to write to the energy
> file
> >> > > nstxtcout = 0 ; period with which to write out positions
> to
> >> an
> >> > > xtc file
> >> > > xtc-precision = 10000 ; xtc precision
> >> > >
> >> > > nstcomm = 10 ; number of steps between COM removal
> >> > > comm_mode = Linear ; remove translational COM (other
> options
> >> > > include angular)
> >> > > comm_grps = System ; remove COM for the whole system
> >> > >
> >> > > Tcoupl = no
> >> > > DispCorr = EnerPres
> >> > > Pcoupl = no
> >> > > -------------------------------
> >> > >
> >> > > Production is after the equilibration simulation and run with grompp
> >> before
> >> > > execution:
> >> > > grompp -f parameters.mdp -c ../Equil/confout.gro -p ../Equil/top.top
> >> -t
> >> > > ../Equil/state.cpt
> >> > >
> >> > > For reference, the lammps simulation runs with the following input
> >> file:
> >> > >
> >> > > ___Lammps input file___
> >> > > units real
> >> > > atom_style full
> >> > >
> >> > > read_data lammps.data # Same position configuration as Gromacs
> >> > >
> >> > > pair_style lj/cut/coul/cut 12
> >> > >
> >> > > pair_coeff 1 1 0.1554 3.166 # kcal/mol = 0.65 kJ/mol
> >> > > pair_coeff * 2 0.0000 0.0000
> >> > >
> >> > > bond_style harmonic
> >> > > angle_style harmonic
> >> > >
> >> > > bond_coeff 1 1000.00 1.0
> >> > > angle_coeff 1 100.0 109.47
> >> > >
> >> > > neighbor 0.0 bin
> >> > > neigh_modify every 1
> >> > >
> >> > > timestep 0.002
> >> > >
> >> > > fix 1 all shake 0.000001 20 0 b 1 a 1
> >> > > fix NVT all nvt temp 300.0 300.0 0.5
> >> > >
> >> > > velocity all create 300 432567 dist uniform
> >> > > fix linmom all momentum 10 linear 1 1 1
> >> > >
> >> > > thermo_style custom step temp etotal ke pe
> >> > > variable mystep equal step*0.002
> >> > > variable myetotal equal etotal
> >> > > variable mype equal pe
> >> > > variable myke equal ke
> >> > > variable myevdwl equal evdwl
> >> > > variable myecoul equal ecoul
> >> > > variable mytemp equal temp
> >> > >
> >> > > fix printenergy all print 100 "${mystep} ${myetotal} ${mype} ${myke}
> >> > > ${myevdwl} ${myecoul}" file sim.energy.equil
> >> > > fix printtemp all print 100 "${mystep} ${mytemp}" file
> sim.temp.equil
> >> > >
> >> > > # NVT equilibration
> >> > > run 250000
> >> > >
> >> > > # NVE production
> >> > > unfix NVT
> >> > > fix NVE all nve
> >> > > reset_timestep 0
> >> > >
> >> > > fix printenergy all print 100 "${mystep} ${myetotal} ${mype} ${myke}
> >> > > ${myevdwl} ${myecoul}" file sim.energy.prod
> >> > > fix printtemp all print 100 "${mystep} ${mytemp}" file sim.temp.prod
> >> > >
> >> > > run 250000
> >> > > -------------------------------
> >> > >
> >> > > Best regards,
> >> > >
> >> > > James
> >> > >
> >> > > --
> >> > > Gromacs Users mailing list
> >> > >
> >> > > * Please search the archive at
> >> > > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> >> > > posting!
> >> > >
> >> > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >> > >
> >> > > * For (un)subscribe requests visit
> >> > > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-usersor
> >> > > send a mail to gmx-users-request at gromacs.org.
> >> > >
> >> > >
> >> > --
> >> > Gromacs Users mailing list
> >> >
> >> > * Please search the archive at
> >> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> >> posting!
> >> >
> >> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >> >
> >> > * For (un)subscribe requests visit
> >> > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> >> send a mail to gmx-users-request at gromacs.org.
> >> --
> >> Gromacs Users mailing list
> >>
> >> * Please search the archive at
> >> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> >> posting!
> >>
> >> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >>
> >> * For (un)subscribe requests visit
> >> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> >> send a mail to gmx-users-request at gromacs.org.
> >>
> >
> >
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>
More information about the gromacs.org_gmx-users
mailing list