[gmx-users] Gromacs and Lammps: energy conservation very different
James
jamesresearching at gmail.com
Tue May 6 02:22:00 CEST 2014
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
More information about the gromacs.org_gmx-users
mailing list