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

Sun May 29 12:42:26 CEST 2011

Hi Andrew,

looking at your protocol, perhaps I can help thermo- and barostating.
This are very important issues as they decide if your final trajectory
samples a canonical phase space that should correspond to the desired
ensemble( Tuckerman and Martyna, J. Phys. Chem. B 2000, 104, 159-178,
DOI:10.1021/jp992433y is a very good review for this topic ).

You are using Nose-Hoover(NH) for both, equilibration and production.
However to drive a system into equilibrium you should perhaps switch to
the v-rescale thermostat. It is more stable and samples the phase
canonically as well as NH. The shorter you make the coupling time, the
faster you reach equilibrium, but the stronger you disturb the
correlations in the system. For the production run the NC scheme allows
to control the time scales were correlations are affected with the
coupling parameter tau_t, while for v-rescale tau_t is just a measure
how strong the thermostat is coupled to the system.

As you want to simulate a canoncial NPT ensemble, the only possibility
for pressure coupling is parrinello-rahman or mttk, while the latter one
is just an extension of the first one. However to equilibrate the system
berendsen should be the method of choice, though it does not sample
canonically. The rules for the coupling parameter are similar as in the
case of temperature coupling. For mttk and pr, tau_p is a measure for
the time scale where correlations are conserved, while for berendsen
tau_p measures the strength of the coupling. For the pressure coupling,
short coupling times introduce very large fluctuations of the box size,
so that it can take also a very long time to reach equilibrium. 

Hope this helps a little bit.

/Flo     

On Sat, 2011-05-28 at 21:40 -0400, Andrew DeYoung wrote: 
> 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
> 

-- 
Florian Dommert
Dipl. - Phys.

Institute for Computational Physics
University Stuttgart

Pfaffenwaldring 27
70569 Stuttgart

EMail: dommert at icp.uni-stuttgart.de
Homepage: http://www.icp.uni-stuttgart.de/~icp/Florian_Dommert

Tel.: +49 - (0)711 - 68563613
Fax.: +49 - (0)711 - 68563658
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 198 bytes
Desc: This is a digitally signed message part
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110529/68ce3b3e/attachment.sig>