Lammps and Gromacs: energy conservation very different

LAMMPS LAMMPS Mailing List Mirror
1

Dear Lammps developers and users,

I am getting significantly better energy conservation in Lammps compared to Gromacs for SPC/E water under an NVE ensemble with a simple cut-off. I enquired on the Gromacs mailing list, and in turn it was asked whether there is something going on behind the scenes in Lammps which is helping to achieve this?

The total energy drift is about 20 kJ/mol/ns/molecule in Lammps, while it is catastrophically exponential in Gromacs (raising 20 kJ/mol/ns/molecule after the initial 50ps). I have done my best to match the conditions of simulation for both software.

Both software:

  • 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

  • Same shake tolerance (also tested LINCS in Gromacs)

  • Same equilibration time (500ps)

  • Same atomic starting positions

  • Same energies (total, PE, KE) after equilibration (easily within 1 standard deviation of fluctuations measured during the final 100ps of equilibration)

I realise even 20 kJ/mol/ns/molecule is poor but it’s enough for my current short calculations, and so if anyone could provide any insight as to why Lammps should be so much more stable I would very much appreciate learning more. I’d like to rule out any possibility that I’m inadvertently doing something in the Lammps input file which helps energy conservation, or that Lammps automatically tries to help somehow.

I include all details regarding input files for Lammps (and Gromacs, for your reference) below. If you need any further information, I would be happy to provide it.

Thank you.

With best regards,

James

The Lammps input file is as follows:

Lammps input

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

2

Dear Lammps developers and users,

I am getting significantly better energy conservation in Lammps compared to
Gromacs for SPC/E water under an NVE ensemble with a simple cut-off. I
enquired on the Gromacs mailing list, and in turn it was asked whether there
is something going on behind the scenes in Lammps which is helping to
achieve this?

The total energy drift is about 20 kJ/mol/ns/molecule in Lammps, while it is
catastrophically exponential in Gromacs (raising 20 kJ/mol/ns/molecule after
the initial 50ps). I have done my best to match the conditions of simulation
for both software.

did you compile gromacs for single or double precision math?

how "good" is your initial configuration?

3

Something else to consider—your time units are not consistent. LAMMPS uses femtoseconds for the time scale in real units, so you’re asking LAMMPS for a 0.002 fs timestep, and using a NVT damping constant of 0.500 fs.

This also has a major effect on the SHAKE algorithm, which does not work well with short timesteps.

—AEI

4

Something else to consider—your time units are not consistent. LAMMPS uses
femtoseconds for the time scale in real units, so you're asking LAMMPS for a
0.002 fs timestep, and using a NVT damping constant of 0.500 fs.

good call. 3 orders of magnitude in the time step should make a huge difference.
perhaps we should have a unit set "gromacs" to make
comparison/conversion easier?

5

Dear Ahmed,

Could you expand on the claim that Shake does not work well for small time step? I would think the opposite. The smaller the time step the smaller the correction necessary.

-s-

6

I should have been more careful in my statement.

The issue with SHAKE arises at the very beginning of the simulation. If you create the molecules far away from their intended configuration, the corrections are done right away and usually all at once. Combined with the small time step, this will lead to exaggeratedly large accelerations and correction forces.

—AEI

7

> Something else to consider—your time units are not consistent. LAMMPS
uses
> femtoseconds for the time scale in real units, so you're asking LAMMPS
for a
> 0.002 fs timestep, and using a NVT damping constant of 0.500 fs.

good call. 3 orders of magnitude in the time step should make a huge
difference.
perhaps we should have a unit set "gromacs" to make
comparison/conversion easier?

I could see how this would be useful, partly for users coming from Gromacs,
but also for the ability to work in SI-like units that are around unity.