[gmx-users] Gromacs and Lammps: energy conservation very different
Mark Abraham
mark.j.abraham at gmail.com
Tue May 6 08:09:18 CEST 2014
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.
More information about the gromacs.org_gmx-users
mailing list