[gmx-users] Parrinello-Rahman and Pairtype1 causing pressure instabilities

[Wesley Michaels]

Hello,

I hope this message finds you well! I'm running a pretty simple simulation
of a polyethylene melt (N = 100, 50 chains in the simulation, T = 500K, P =
1atm, NPT ensemble, dt = 1 fs, OPLS-AA FF). My simulations are exhibiting
large, unstable volume fluctuations when I use the combination of the
Parrinello-Rahman barostat and pairtype = 1 for 1-4 pairs. I don't see this
behavior with {PR barostat and pairtype = 2} or {Berendsen barostat and
pairtype = 1}.

I've specified the 1-4 pairs explicitly in the "[ pairs ]" section and set
"gen_pairs" to "yes" in the defaults section. I'm sure I'm missing
something obvious, but here are the troubleshooting steps I've tried:

1) ensure the simulation is equilibrated (up to 100 nanoseconds with
Berendsen barostat beforehand, which is more than plenty I believe);
2) use a smaller time step (down to 0.01 fs);
3) increase the pressure coupling time constant (up to 10 ps);
4) using both GROMACS 2018 and GROMACS 2019.

As far as energies go, I don't think there are worryingly large terms that
would blow the simulation up. Here's an example prepared with `gmx energy`:
Time (ps) Bond G96Angle RB LJ-14 Coulomb-14 LJ (SR) Coulomb (SR) Total
Energy Conserved En. Volume Density
0 3.10E+04 5.40E+04 2.29E+04 8.76E+03 -1.66E+03 -4.71E+03 5.27E+03 2.17E+05
2.17E+05 291.0 400.8
2 5.11E+04 5.71E+04 8.24E+03 5.48E+02 -3.00E+03 -3.68E+03 9.01E+03 2.12E+05
1.02E+06 3720.8 31.3
4 2.56E+03 9.53E+03 8.18E+03 4.80E+03 -4.16E+03 -4.27E+03 9.40E+03 1.19E+05
1.25E+06 21855.1 5.3
6 4.76E+03 1.48E+04 1.15E+04 5.36E+03 -3.71E+03 -4.02E+03 8.68E+03 1.31E+05
1.52E+06 88422.7 1.3
8 6.58E+03 2.07E+04 1.48E+04 6.13E+03 -3.27E+03 -4.02E+03 7.95E+03 1.44E+05
1.80E+06 246764.8 0.5
10 1.05E+04 2.80E+04 1.64E+04 6.46E+03 -3.03E+03 -3.98E+03 7.61E+03 1.57E+05
2.08E+06 538599.2 0.2


Eventually, the volume becomes unphysically large and the density extremely
low. I'm sure I'm missing something obvious here (sorry for yet another "my
simulation is blowing up" question), so any insight into other steps I
could try would be greatly appreciated. If I can provide any additional
information (e.g. topology/mdp files), just let me know.

Thanks,
Wes