[gmx-users] Gromacs and Lammps: energy conservation very different
Michael Shirts
mrshirts at gmail.com
Tue May 6 04:13:17 CEST 2014
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.
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.
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.
>
>
More information about the gromacs.org_gmx-users
mailing list