Parrinello-Rahman oscillates and blows up, even after Berendsen equilibration

User discussions
1

GROMACS version: 5.0.7
GROMACS modification: No

Hi all,

Recently I started doing MD in GROMACS. I’d like to simulate a box of oxygen molecules in an NpT ensemble. In order to do so, I use the TraPPE force field from The Siepmann Group (http://trappe.oit.umn.edu/). The implementation can be found in the attachements, see topol.top.

Before starting the actual NpT simulation, I start with energy minimization followed by an NvT simulation to relax the system. After 0.5ns of NvT simulation I switch on the Berendsen barostat to pressure the system to 1.0 bar. Velocity rescaling is used as thermostat during these steps. After Berendsen (1ns) I switch to the Parrinello-Rahman barostat and the Nosé-hoover thermostat with respective typical time scales of 1ps and 0.1ps. Doing so makes the box eventually blow up. I get

“Warning: Pressure scaling more than 1%. This may mean your system is not yet equilibrated. Use of Parrinello-Rahman pressure coupling during equilibration can lead to simulation instability, and is discouraged.”

and eventually a fatal error:

“Fatal error: 3 particles communicated to PME rank 0 are more than 2/3 times the cut-off out of the domain decomposition cell of their charge group in dimension x. This usually means that your system is not well equilibrated.”

I can mitigate the blow up by increasing the time constant for the PR barostat, but it still shows fluctuations. Besides these fluctuations, the density is off by a factor of 2. Can this discrepancy be explained by incorrect time constants?

Another ‘fix’ is to use the Berendsen barostat, but this is not recommended as it does not truly simulate an NpT ensemble. It however works, although the density is again not correct. I am not sure why this is.

So to summarize, my questions are:

  • What can I do to make sure that PR works with a time constant of 1ps, even after equilibration with Berendsen?
  • Increasing the time constants seem to work, but what effect does it have on the ‘physical validity’ of the simulation?
  • Is there an obvious reason the density in the NpT simulation with Berendsen is incorrect?

Help is very much appreciated, thanks in advance!

Files can be found on: GitHub - Stanvk/gromacs_dioxygen: Oxygen simulation in GROMACS

2

You coupling times are much too short for second order coupling algorithms. I would suggest to use a time of 2 ps for NH and a time of at least 5 ps for PR.

But why are you using a decade old version of GROMACS? In the newest releases you could use the v-rescale thermostat and c-rescale barostat instead.

3

Thanks @hess for your answer. Indeed I notice that increasing the time coupling is beneficial to the stability. Unfortunately I have no access to the latest version of GROMACS yet.

Many do use the PR-barostat which seems to work well. However, I still get a wrong density when using PR instead of Berendsen. What could be the reason for this behaviour?

Thanks in advance,

Stan

4

The only reason I can the that using PR would result in the density is that the coupling time, which should gives oscillations with that period if the compressibility is set correctly, causes resonances with either the thermostat or some mode in the system. What is the autocorrelation time of the box volume? Does increasing tau_p solve the issue?

5

In my case, Nosé-hoover thermostat created this problem, I changed it to V-rescale during the MD run and it solved the pressure scaling issue.