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

[Andrew DeYoung]

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?

(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.

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

Thank you very much for your time.  I truly appreciate it.

Sincerely,

Andrew DeYoung
Carnegie Mellon University