I am experiencing a reproducibility issue between two versions of GROMACS and would appreciate your input.
We performed molecular dynamics simulations of tripeptides using GROMACS 2020.1 with the CHARMM36m force field. Over time, we observed the formation of small peptide clusters, as expected.
However, when restarting the simulation from the last frame of the trajectory and running it with GROMACS 2023.3 (using exactly the same parameter files, topology, and input settings), the clusters disaggregate very quickly.
Do you have any suggestions of what could be the cause of this issue?
Thank you for your answer.
No, I haven’t yet continued the simulation with GROMACS 2020.1, but I did try using version 2021.2, and the results look consistent — the clusters remain stable throughout the simulation (up to 750 ns).
On the other hand, with version 2023.3, we observe an unexpected behavior where the clusters rapidly disaggregate.
I understand the need for running replicates. However, the results so far are quite unusual — the disaggregation effect occurs within just ~100 ps, which seems physically unreasonable.
For clarity, I’m attaching the SASA profiles from the simulations run with GROMACS 2020.1 and 2023.3.
I ran the simulation with GROMACS 2020.1 on CPU using 1 MPI rank, and with GROMACS 2023.3 on GPU using 4 MPI ranks. I also performed tests with other versions: 2023.3 on CPU with 8 MPI ranks, 2025.1 on CPU with 16 MPI ranks, and 2025.2 on CPU with 16 MPI ranks.
Unfortunately, I observe very similar SASA behavior across all versions.
Looks all reasonable. The only difference between versions I see is that nsttcouple and nstpcouple might be set differently in new versions. Could you run a check with 2025 with nsttcouple and nstpcouple set in the mdp options to the value reported in the log file for 2020.1?
I ran the test you suggested using the same default values for nsttcouple and nstpcouple as in GROMACS 2020.1. The results are now consistent with those from the original simulation (I’m attaching the SASA as confirmation).
Thank for you this test. Now we know what causes this. Could you also check nsttcouple and nstpcouple independently to see which of the two is the issue?
I don’t think the temperature coupling should be problematic. But a tau_p of 5 ps results in 50 integration steps per period with nstpcouple=100. I don’t see how this could cause issues, unless the actual compressiblity is much lower, for which I see no reason.
I ran two tests: one using the default nsttcouple from GROMACS 2020.1 and the default nstpcouple from GROMACS 2025.2 (labeled “nstt only” in the figure), and another using the default nsttcouple from GROMACS 2025.2 and the default nstpcouple from GROMACS 2020.1 (labeled “nstp only”).
It turns out that the nsttcouple value from version 2025.2 does not introduce any artifacts, whereas the nstpcouple value from the same version does cause the previously discussed artifact. Therefore, the issue can be attributed specifically to the nstpcouple
Can you run an autocorrelation analysis on the box volume from a correct simulation? I’d like to see what the period of the box fluctuations is compared to tau_p.
From the autocorrelation analysis, I found a fluctuation period of approximately 2 ps.
For reference, the pressure coupling constant (tau_p) was set to 5.0 ps.
Ah, very good! Then that explains things and there is no bug in GROMACS. But this must mean that the actual compressibility is a factor (5/2)^2 lower than what you supplied in the mdp file. That is a large factor. Do you only have peptides and water in the system?
I have 538 peptides (500 in the protonated state and 38 in the neutral state), 38407 water molecules, and 500 counterions (Cl⁻). Yes, peptides interact forming small clusters across the unit cell.
So my guess is that that decreases the compressibility. But the effect is stronger than I would expect.
It would be nice if we could detect incorrectly set compressibility. We could analyze the autocorrelation of the pressure and volume, but it could be difficult to draw reliable conclusions from that.
The best solution is using the C-rescale barostat instead.
