Change in Volume Fluctuations / Compressibility for TIP4P/2005 Between GROMACS 2016 and 2024

User discussions
,
1

GROMACS version: 2024
GROMACS modification: No

My research group recently upgraded our cluster from GROMACS 2016 to GROMACS 2024.1, and I’ve been running into a consistent but puzzling difference in the isothermal compressibility of pure water.

System:

  • 1000 TIP4P/2005 waters

  • 298.15 K

  • NVT → NPT (Parrinello–Rahman, isotropic)

  • Nose–Hoover thermostat

What’s happening:
On the old system (2016 build), this exact setup gave a κₜ ≈ 4.62 × 10⁻⁵ bar⁻¹, matching experiment perfectly.
When I moved to 2024.1, using the same workflow and the same force-field files, my κₜ jumped to around 5.0–5.2 × 10⁻⁵ bar⁻¹—a roughly 10 % increase in volume fluctuations. The density is still correct, so the problem is only in σ(V).

What I’ve tried so far:

  • Explicitly set verlet-buffer-tolerance = 0.005 and verlet-buffer-pressure-tolerance = -1 (since defaults changed).

  • Matched nstpcouple, nsttcouple, and other timing parameters to the old run.

  • Verified identical force-field and water model files.

  • Re-ran multiple times in 2024 → results are fully reproducible (so it’s not stochastic noise).

  • Finally, I took the exact TPR that was generated under GROMACS 2016 and ran it unchanged with the 2024 engine.

    • Result: the newer engine reproduces the 2024-style κₜ (≈ 5 × 10⁻⁵), not the old 2016 value.

    • So the discrepancy appears to stem from the engine version itself, not grompp or the MDP.

MDP used (explicit values for changed defaults):

integrator               = md
dt                       = 0.002
nsteps                   = 5000000
nstenergy                = 50
nstlog                   = 500
nstxout-compressed       = 250
constraint-algorithm     = lincs
constraints              = h-bonds
cutoff-scheme            = Verlet
verlet-buffer-tolerance          = 0.005
verlet-buffer-pressure-tolerance = -1
comm-mode                = Linear
nstcomm                  = 100
coulombtype              = PME
rcoulomb                 = 1.0
fourierspacing           = 0.12
pme-order                = 4
ewald-rtol               = 1e-05
vdwtype                  = Cut-off
rvdw                     = 1.0
coulomb-modifier         = Potential-shift
vdw-modifier             = Potential-shift
DispCorr                 = EnerPres
tcoupl                   = Nose-Hoover
tc-grps                  = System
tau-t                    = 2.0
ref-t                    = 298.15
pcoupl                   = Parrinello-Rahman
pcoupltype               = isotropic
tau-p                    = 4.0
compressibility          = 4.46e-5
ref-p                    = 1.0
nstpcouple               = 10

Summary:

  • The same TPR gives slightly larger volume fluctuations when run with the 2024 engine than with 2016.

  • The difference is reproducible across runs and nodes.

  • I’m wondering if any changes were made to the Parrinello–Rahman implementation, LINCS behavior, integration scheme, or precision handling that could alter volume fluctuations at the 1–2 % level.

Has anyone else observed this when upgrading from older releases? Any hints on what part of the integrator/barostat chain might have changed its effective coupling strength would be very helpful.

Thanks in advance for any insight!
— Jakob Schanzer

2

I can’t help but I remember a similar discussion recently. Take a look here and in the related issue here! :)

3

The question is what the correct value is, which is not necessarily the experimental value.

The answer above refers to an issue with the increase of the default value of nstpcouple from 10 to 100. But you set it explicitly to 10. Then I have no idea what could have changed between the two versions. As you can see in here I observe 2-4% too high fluctuations for SPC/E with tau_p=5 and nstpcouple=10. With tau_p=4 this will be slightly higher. Please try nstpcouple=1 and report back what you observed.

4

Thanks so much for this. I tried a very short (2 ns) production with just changing this parameter in V 2024, and got my value down to ~4.5e-5 /bar. Much more in line with what I expect from my previous work. I’m currently doing a 100 ns run to allow for an extended equlibration, but this seems to maybe have done the trick. Now the question is, as you noted, which value best models reality?

5

So you got a change of 10% by changing nstpcouple from 10 to 1 with tau_p=4? That is a bit more than I expected.

You can look up the experimental value. But am MD engine should reproduce the values of the model, not reality.

6

Once I did the longer run of 100 ns, it leveled off to 4.67e-5, which essentially reproduces what I was getting when running this system on version 2016

7

But that’s still a difference of 9%. That’s a factor two larger than the 4% tolerance which I set as a maximum acceptable deviation for getting better performance.

I general I would advise to use the C-rescale barostat, which is less sensitive to the step size.