[gmx-users] Plan for energy minimization, NPT equilibration, and NPT dynamics

[Andrew DeYoung]

Hi,

In the past few weeks, I have run some simulations of 254 water molecules.
Although I got some results, I did not do things in the best way.  For
example, I threw away my equilibrated ensemble before doing my production
run, as was kindly pointed out to me on this mailing list last week.  In
trying to think more carefully about the input parameters and the commands I
use, I have done some more reading, particularly in the very helpful
lysozyme tutorial
(http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lyso
zyme/index.html) by Justin Lemkul.  

In planning another simulation, I would like to understand better reasonable
parameters to use and the correct commands to use.  If you have time, can
you please look at the pdf file at
http://www.andrew.cmu.edu/user/adeyoung/may28.pdf that I have created?  If
you have time, I have some questions about it.  On page 1, I have typed some
possible parameter files (adapted from Justin Lemkul's lysozyme tutorial);
on page 2, I have typed a list of commands (also adapted from the lysozyme
tutorial) I may use in my run.  I plan to do energy minimization, NPT
equilibration, and NPT MD (for NPT, I plan to use Berendsen pressure
coupling and Nose-Hoover temperature coupling) using Gromacs 4.5.4.  If you
have time, here are my questions about the pdf file:  
---
(i) In equilibration, what is a reasonable value of nstlist to use?  The
Gromacs manual on page 18 says that "nstlist is typically 10," but is it
true that, in principle, the smaller that nstlist is (that is, the more
often the pairlist is updated), the more accurate the simulation is, at
greater computational cost?  As a rule of thumb, should I use nstlist = 10
or something smaller?  I have access to a cluster, so I don't have too many
worries about computational cost, to a reasonable extent.  

(ii) In equilibration, what is a reasonable value of the time constant tau_p
to use for Berendsen pressure coupling?  The default is 1 ps, although I
have seen tutorials where people use tau_p as high as 2 ps or as low as 0.1
ps.  The Gromacs manual on page 32 gives a relationship for the change in
pressure with respect to time for Berendsen pressure coupling:  

dP/dt = (P_0 - P)/tau_p

where P_0 is the reference pressure (ref_p in Gromacs).  So it seems that
the smaller the tau_p, the faster the pressure relaxes toward the reference
pressure.  So, naively, it would be best to use a small tau_p.  However, am
I correct in thinking that, on the other hand, if I make tau_p too small for
a fixed step size dt, significant inaccuracies may result?  

(iii) In equilibration, is it reasonable to set DispCorr = no?  I won't be
using my equilibration results for any detailed analysis; I only wish to
generate a (hopefully) reasonable ensemble and relax the pressure to the
reference pressure.

(iv) Similarly, in production, is it reasonable to set DispCorr = no?  I
suppose that it depends on what one ultimately wants to calculate using the
final MD production results.  Since I am still at the stage of getting
familiar with Gromacs, I will be doing only simple analyses (RDF, MSD,
velocity autocorrelation function, and extracting plots of pressure and
density versus time from the final edr file); I will not be doing
complicated free energy calculations, for example.  

(v) What does the continuation parameter do?  Looking under section 7.3.18
(the section called "Bonds") on page 181 of the Gromacs manual, it seems
that this is mainly if I choose to constrain particular atoms or bonds.
Here, since my system consists only of water, I will not use any
constraints, I think.  However, I am getting myself confused, I think.  I
do, in fact, want my (NPT) MD production run to be a "continuation" of my
(NPT) equilibration; does this mean that the continuation parameter should
be set to "yes"?  I am not sure, because I am not sure if the continuation
parameter is exactly relevant to whether my production run is a
"continuation" of my equilibration (that is, I want my production run to use
the ensemble that was created during equilibration).  

(vi) On page 2 of my pdf file, I have highlighted parts of steps (6) and (7)
in magenta.  Could you please look at these parts and help me see whether or
not I am using the correct commands to, firstly, create a checkpoint file at
the end of equilibration and to, secondly, use that checkpoint file as input
into the grompp step of my MD production run?  About a week ago, Justin
pointed out that I should use the -t option in the grompp program to do
this.  But, I want to make sure that I understand correctly...

...In the entry for grompp in the Gromacs manual (for example, see
http://manual.gromacs.org/current/online/grompp.html or page 300 in the
manual), the listing for -t says "Full precision trajectory: trr trj cpt" in
the description column.  Does this mean that I can give grompp EITHER the
previous trr trajectory file OR the cpt checkpoint file that I created using
-cpo in equilibration?  Or must I definitely pass the cpt file to grompp in
order to continue using in production the ensemble that was created during
equilibration?  
---

Thank you very much for your time!

Andrew DeYoung
Carnegie Mellon University