[gmx-users] 4xRe: Energy conservation of rigid water models

Janne Hirvi janne.hirvi at joensuu.fi
Mon Feb 20 11:11:53 CET 2006

Thanks for all help!

Now I understand at least something about the factors which have an effect on
energy conservation but however I still have one question. I need to calculate
the intrinsic level of energy conservation in a water model with employed time
step of 2fs. So should I use same parameters as in NVT/NPT simulations or
should I eliminate all other errors in NVE ensemble like vdw cut-off and single
precision etc. to be able to see the real effect of time step of 2fs? 


Original messages:


I have read several messages about energy conservation problems and solutions
but however I am not sure if I have done something stupid. I have used time
step of 2.0fs for rigid SPC and SPC/E water models in simulations of bulk water
(NPT) and water droplets on frozen polymer surfaces (NVT). Everything seemed to
be just fine but when I tested energy conservation of bulk model in NVE
ensemble I was terrified.

NVE bulk simulation of 1372 rigid SPC water molecules with 2fs time step (total
time 500ps) resulted in huge drift in total energy. I used same parameters as
for other simulations and even those aren't optimal I expected at least
something better. Energy conservation seemed to be adequate only when I used as
short time step as 0.5fs.

Time step (fs)      Total Energy (kJ/mol)      Drift

2.0    			-45596.6 		1235.5
1.0     			-46147.5 		330.974
0.5        		-46327.6 		22.7932


comm-mode		 = None or Linear
nstlist                         	 = 5 or 10
rlist                   	 = 1.2
coulombtype		 = PME
rcoulomb                	 = 1.2
vdw-type                	 = Cut-off
rvdw                     	 = 1.2
fourierspacing         	 = 0.12
pme_order                    = 4
ewald_rtol                    = 1e-05
ewald_geometry          = 3d
tcoupl                           = No
Pcoupl                          = No

Bad energy conservation indicates that employed time step of 2.0fs is too long
but however other results from NPT bulk simulations are reasonable. Moreover I
think that its quite normal to use this long time step for rigid water
molecules. So I am willing to know if this large energy conservation problems
are acceptable and are they at least "partially compensated" in NPT or NVT
ensemble when temperature coupling prevents heating of the system?

Thanks for any help or comments!


Have you compiled gromacs in single precision? I think double precision
compilation is critical for energy conservation in NVE.

Milton Sonoda

Yes, I have used only single precision compilation because it has been enough
accurate until now (NPT/NVT). Actually I just want to know if employed time
step of 2fs is applicable and so I think that I should use same compilation
also for energy conservation simulation in NVE ensemble. Or may I prove
suitability of this long time step with double precision compilation even I
have done other simulations with single precision compilation?


It is important to realize that you need this kind of precision only in special
cases, for instance when you want to derive properties from the fluctuations of
T, P, E.

The algorithms in GROMACS are not perfect, that is, there are better algroithms
for integrators etc. We are working on these things.

However, what is implemented, i.e. leap-frog, when used carefully, does conserve
energy. Things that spoil energy conservation are:

- T coupling
- P coupling
- Cut-offs (both LJ and Coulomb)
- Constraints (but not SETTLE as used for water)
- Long time steps (in relation to the fastest vibrations)
- Too long intervals between neighborsearching
- To some extent, single precision (due to the fact that you sometimes add small
forces to large one, and by this lose precision).


I strongly sugest you to try double precision compilation for energy
conservation in NVE.
Single precision is enough for NPT/NVT due to convenient algorithms (velocity
and volume reescaling) that keeps the system at the desired termodynamic state
(constante pressure and constante temperature) masking the truncation errors of
the single precision simulation.

In NVE no algorithm is used, so the double precision is necessary to
prevent/postpone energy drift due to this error acumulation.

Milton Sonoda


i'm certainly not a gromacs wizard, but i found that i had better luck for
things like this when i reduced the ft spacing to 0.08 and increased pme order
to 6.  (but i think you need to make sure that you have the patched version of
pme.c for gromacs-3.3.)  and of course i would agree that running double
precision will help.  i hope that's helpful.


will noid

I'd suggest compiling in double precision, and running with the following

dt                                   = 0.002
nstlist                             = 10
rlist                                = 1.2
coulombtype                 = PME
rcoulomb                       = 1.1
vdw-type                       = Switch
rvdw                              = 1.1
rvdw-switch                  = 1.05
tcoupl                            = No
pcoupl                           = No
constraints                     = all-bonds
constraint-algorithm     = shake
unconstrained-start       = no
Shake-SOR                   = no
shake-tol                        = 1e-12

This should give fairly good results for energy conservation.

You may possibly need to used improved Ewald parameters; I had decent success

fourierspacing            = 0.10
pme_order                  = 6
ewald_rtol                  = 1e-06

Michael Shirts

See Chiu, S.-W., Clark, M., Subramaniam, S., and E. Jakobsson. 2000.  Collective
motion artefacts arising in long-duration molecular dynamics simulations.   J.
Comp. Chem. 21:121-131

Eric Jakobsson

Janne Hirvi, MSc(Physical Chemistry), Researcher
University of Joensuu, Department of Chemistry, P.O.Box 111 80101 Joensuu, FI
Tel: +358 13 2514544 & +358 50 3474223
E-mail: Janne.Hirvi at joensuu.fi & hirvi at cc.joensuu.fi

More information about the gromacs.org_gmx-users mailing list