[gmx-users] Reference pressure and pressure fluctuations in an NPT simulation

[Justin A. Lemkul]

Andrew DeYoung wrote:
> Greetings,
> 
> I have been running simulations of 254 SPC/E water molecules using the OPLS
> force field.  As Nilesh mentioned earlier today, I am using the NPT
> ensemble.  I am setting the reference pressure at 1 bar, but when I do a
> sequence of minimization, equilibration, and dynamics steps, and then
> finally use g_energy to determine the average pressure, I find an average
> pressure of ~2 bar, rather than 1 bar.  
> 
> Earlier today, Justin pointed out that this probably means that the system
> has not been equilibrated long enough, so I did another run with longer
> equilibration.  However, I get a similar, puzzling result of ~2 bar.
> Finally, I tried changing the temperature and pressure coupling methods;
> this time, I get a seemingly more reasonable result for average pressure,
> but, as I will describe below and in my PDF file, still somewhat puzzling to
> me (I am new to the field of molecular dynamics).
> 
> If you have time, could you please look at my PDF file at
> http://www.andrew.cmu.edu/user/adeyoung/may13.pdf that summarizes what I
> have tried?  In case it is helpful, also, here is a text description of what
> I have tried:
> 
> (i) T coupling = v-rescale; P coupling = parrinello-rahman.  Reference
> pressure = 1 bar.  1 ns equilibration; 2 ns dynamics.  Gromacs tells me that
> the average pressure is ~2.77 bar, but when I use g_energy to extract
> pressure as a function of time to an xvg list and then use software such as
> Mathematica or Matlab to compute the mean of the list, I find that the
> average pressure is ~1.24 bar.
> 
> (ii) Use the same parameters as before, except equilibrate for 5 ns instead
> of 1 ns.  Gromacs says that the average pressure is ~2.82 bar, but when I
> extract the pressure data and compute the mean of the list, I find that the
> average pressure is ~6.81 bar.
> 
> (iii) Use the same parameters as in (ii), except use berendsen temperature
> coupling and berendsen pressure coupling.  Gromacs says that the average
> pressure is ~1.00 bar, but when I extract the pressure data and compute the
> mean of the list, I find that the average pressure is ~3.19 bar.
> 
> In my PDF file at http://www.andrew.cmu.edu/user/adeyoung/may13.pdf, I have
> plotted pressure vs time for the dynamics runs of (i), (ii), and (iii), and
> these plots show VERY large magnitude oscillations.
> 
> If you have time, I have three questions about these results:
> 
> (1) Does it seem reasonable that I obtain an average pressure closer to my
> reference pressure (1 bar) when I use berendsen temperature coupling and
> berendsen pressure coupling -- instead of v-rescale temperature coupling and
> parrinello-rahman pressure coupling?
> 

Yes.  Berendsen coupling methods produce artificially narrow distributions of 
temperature and/or pressure values, and hence using them produces an incorrect 
ensemble.  There are more algorithmically detailed reasons for this, but that's 
the net effect.  It is for this reason that Berendsen methods can be extremely 
useful for initial equilibration, to relax the T and P values to the target 
values relatively quickly, after which a more rigorously correct ensemble can be 
employed with other thermostats and barostats.

> (2) When I use g_energy to extract pressure vs time, and then compute the
> mean of the list of pressures, why is my answer so different from Gromacs'
> "black box" calculation (that is, using g_energy to have Gromacs simply
> print the average pressure over the dynamics run) of the average pressure?  
> 
> For all of my runs, I am using a step size of 1 fs (=0.001 ps).  However,
> when I extract the pressures (using the default settings), I get a pressure
> value for every 0.1 ps, not 0.001 ps.  Could it be that the difference in
> average pressures is due to the fact that "under the hood" Gromacs is using
> the pressure data at every step (0.001 ps), instead of every 0.1 ps?  But,
> even if this is true, I still wouldn't necessarily expect such a big
> disagreement between the two calculations.

The output from g_energy is correct over all MD steps.  The resulting .xvg file 
contains output only every nstenergy steps.  So unless you've got nstenergy = 1, 
the results will be somewhat different.  It looks like you're setting nstenergy 
= 100.

For more well-behaved quantities (i.e. things that do not fluctuate insanely 
like pressure), the difference is very small.  For pressure, I do see the same 
sort of results.  For all of my protein-in-water type systems, I have average 
pressures (from g_energy) around 1.0 bar that fluctuate on the order of 10^2. 
These systems contain tens of thousands of atoms, so I expect the fluctuations 
to be considerably less than your small system.  You can test these effects by 
expanding your system.  The fluctuations should decrease a lot, as explained here:

http://www.gromacs.org/Documentation/Terminology/Pressure

The fluctuations you're observing are a bit larger than what might be expected, 
but still not completely unreasonable.

> 
> (3) When I plot pressure vs time (as I have done in my PDF file), why is
> there such a large magnitude of pressure fluctuation?  For example, in my
> runs, the maximum pressure is on the order of 3000 bar, whereas the minimum
> pressure is on the order of -3000 bar.  These pressures are unreasonably
> large in magnitude (despite the fact that the average pressure nevertheless
> turns out to be of the correct order of magnitude in the long run).  Is this
> true?  Also, is "negative pressure" unphysical?  Or, does "negative
> pressure" correspond to "compression" and "positive pressure" corresponds to
> "expansion," or something like this?
> 

Something like that.  See the link above.

Are you using Gromacs 4.0.7?  There have been improvements to the P-R barostat 
since then, and version 4.5.4 contains at least one potentially relevant bug fix 
(which may only affect versions in the 4.5.x series, but I am not sure).

-Justin

> Thank you very much for your time.  I truly appreciate it.
> 
> Sincerely,
> 
> 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

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