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
========================================
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
