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

[Mark Abraham]

On 14/05/2011 11:02 AM, Andrew DeYoung wrote:
> Greetings,
>
> I have been running simulations of 254 SPC/E water molecules using the OPLS

Bear in mind that this is a tiny system, and that temperature and 
pressure are macroscopic quantities. Reliable measurements of them need 
long times and lots of statistics.

> 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,

I think that is partly fortuitous, and partly the nature of Berendsen 
coupling.

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

See second para of g_energy -h.

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

Looks normal. I see that Justin's reply has dealt with the issues below.

Mark

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