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
========================================
--
gmx-users mailing list [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [email protected].
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
