[gmx-users] Improving Energy Conservation in NVE Simulation of Water

[Andrew DeYoung]

Hi, 

I have been simulating 1000 SPC/E water molecules in the NVE ensemble (by
using tcoupl = no and pcoupl = no).  In these simulations, I am considering
only Lennard-Jones interactions; I have "turned charges off" by setting the
partial charges on oxygen and hydrogen to zero in a local copy of spce.itp
that I include in my topology file.  

If you have time, I have posted a pdf file containing a summary of what I
have tried so far: http://www.andrew.cmu.edu/user/adeyoung/july27/july27.pdf

At first, I tried using vdwtype = Cut-off.  However, Gromacs gave me a note
message: 
---
NOTE:
You are using a cut-off for VdW interactions with NVE, for good energy
conservation use vdwtype = Shift (possibly with DispCorr)
---

I ignored this note message and ran the simulation anyway.  Indeed, as
Gromacs warned me, the (total) energy conservation is very bad; the total
energy decreases seemingly significantly in what appears to be an almost
linear fashion.  I have done energy minimization, equilibration, and
production dynamics, but for simplicity, I have drawn the equilibration and
production dynamics results in the same color in the pdf file that I have
posted.  

(As far as energy minimization, I did a series of two energy minimization
steps.  First, I used integrator = l-bfgs with dt = 0.001 ps, nsteps =
10000, emstep = 0.01 nm, and emtol = 10.0 kJ mol^-1 nm^-1.  I got the
following results:

Low-Memory BFGS Minimizer converged to Fmax < 10 in 0 steps
Potential Energy  = -9.3022546e+02
Maximum force     =  5.5933684e-01 on atom 1123
Norm of force     =  1.9135195e-01

So, it seems that the energy of the configuration is quite low--the water
molecules are relatively far apart from one another (the simulation box is
quite large).  

Next, I used integrator = steep (passing the ending configuration from the
l-bfgs minimization to grompp using -t followed by the trajectory file) with
dt = 0.001 ps, nsteps = 10000, emstep = 0.01 nm, and emtol = 10.0 kJ mol^-1
nm^-1.  I got the following results:

Steepest Descents converged to Fmax < 10 in 1 steps
Potential Energy  = -1.2968502e+03
Maximum force     =  7.0817396e-02 on atom 820
Norm of force     =  2.6016856e-02

This configuration seems probably acceptable.  So much for energy
minimization.)

Secondly, I used vdwtype = Shift as the note message recommended.  I tried
two different sets of values of rlist, rvdw_switch, rvdw, and rcoulomb.
But, in both cases, energy conservation is bad--almost as bad as when I used
vdwtype = Cut-off.  

If you have time, can you please look at the pdf file that I have posted
(http://www.andrew.cmu.edu/user/adeyoung/july27/july27.pdf)?  Do you have
any suggestions of other vdw types or vdw parameters that I should try to
hopefully improve energy conservation?  On page 1, I have displayed the
total energy results.  On pages 2 and 3, I have pasted a typical .mdp file
that I used and also the results of other quantities from the simulations.  

One thing that I notice is, on the bottom of page 2, it appears that the
_potential_ energy conservation is good, whereas the kinetic energy
conservation is quite bad.  Do you know, what sort of clue should this fact
give me as to what is going on with my system?  

Thank you.

Andrew DeYoung
Carnegie Mellon University