[gmx-users] Recommended parameters for NVE simulation of SPCE water

[Andrew DeYoung]

Hi,

I am a Gromacs novice and am trying to simulate 1000 SPC/E water molecules
(with the usual charges) in the NVE ensemble. I am trying to obtain good
energy conservation. 

My configuration contains a 10-by-10-by-10 cubic crystal of SPC/E waters,
with a lattice constant of approximately 4.7 Angstroms.  I first run energy
minimization, using integrator = steep, emtol = 1000 kJ/mol/nm, emstep =
0.01 nm, and coulombtype = PME.  This energy minimization converges in only
1 step -- I guess because my value of emtol is so large and because the
water molecules are relatively far apart from each other.

Then, when I run my initial equilibration/dynamics (I am calling my run both
equilibration and dynamics because, for the time being, I am just doing a
single run) using coulombtype = PME, I get the following two notes messages:
------
NOTE 1 [file nve-eq.mdp]:
You are using a cut-off for VdW interactions with NVE, for good energy
conservation use vdwtype = Shift (possibly with DispCorr)

NOTE 2 [file nve-eq.mdp]:
You are using a cut-off for electrostatics with NVE, for good energy
conservation use coulombtype = PME-Switch or Reaction-Field-zero
------

This warning is also discussed here:
http://www.gromacs.org/Documentation/Terminology/NVE

And, indeed, just as Gromacs says, I get terrible energy conservation using
coulombtype = PME (I use g_energy to extract the "Total energy" from the
.edr file).  I have posted plots of my results using coulombtype = PME here:
http://www.andrew.cmu.edu/user/adeyoung/july13/july13.pdf .  As you can see
in the light red plot on the right-hand side of that PDF file, energy
conservation is terrible using coulombtype = PME, even when I do a 2 ns run;
total energy appears to increase nearly linearly as a function of time.

If you have time, I am wondering if you have any advice about the following
two questions:

1.  NOTE 1 above suggests that I use vdwtype = Shift.  When I do this, do
you recommend that I apply long range dispersion corrections for both energy
and pressure, using DispCorr = EnerPres, or for only energy, using DispCorr
= Ener?  Typically, for various (non-NVE) calculations, I have been using
DispCorr = no, but I am not sure if this is a good idea.  Pages 97-98 of the
Gromacs 4.5.4 manual seem to suggest that the energy correction due to
DispCorr is small and usually only significant for free energy calculations
(which I will not be doing here).  As a rule of thumb, do you typically turn
dispersion corrections off?

2.  NOTE 2 above suggests that I use either coulombtype = PME-Switch or
coulombtype = Reaction-Field-zero.  Do you have any advice or
recommendation?  My concern (and I am very much a novice in MD) about using
coulombtype = Reaction-Field-zero is that the manual says
(http://manual.gromacs.org/current/online/mdp_opt.html#el) that it can only
be used with an infinite dielectric constant.  I know that water has a
relatively high dielectric constant (~80 at ~298 K), but strictly speaking
it is not infinite.  On the other hand, about coulombtype = PME-Switch or
Switch, the manual says
(http://manual.gromacs.org/current/online/mdp_opt.html#el) that "Switching
the Coulomb potential can lead to serious artifacts, advice: use
Reaction-Field-zero instead."  It doesn't explicitly specify what kinds of
serious artifacts I can expect, but it does sound like a dangerous choice.
Do you have any experience with NVE simulations?

Thank you very much for your time!

Andrew DeYoung
Carnegie Mellon University