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

[Justin A. Lemkul]

Andrew DeYoung wrote:
> 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.
> 

Try minimizing further.  1000 kJ mol^-1 nm^-1 is generally OK for more robust 
systems under "normal" conditions, but it can't hurt to go further, if possible.

> 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?
> 

I always use it.  Any time you truncate van der Waals interactions, you're 
losing a long-range component of the interactions.  The effect may be small, but 
if you can be more accurate, why not do it?  Many of my systems have lipid 
membranes, which are particularly sensitive to such terms, so maybe for pure 
water this isn't strictly required.  I certainly don't think it would hurt.

> 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?
> 

Sorry to say I can't comment here.  It would seem that the documentation 
certainly spins you in circles with this, though.  You can easily do some short 
tests of any of the electrostatics methods, though, to see which one may give 
you the best results.  Somewhat empirical, I know, but at least you can narrow 
it down.

-Justin

> Thank you very much for your time!
> 
> Andrew DeYoung
> Carnegie Mellon University
> 

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================