[gmx-users] Gromacs and Lammps: energy conservation very different
James
jamesresearching at gmail.com
Tue May 6 11:42:29 CEST 2014
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-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