Discrepancy in Tripeptide Aggregation Behavior Between GROMACS 2020.1 and 2023.3

User discussions
1

GROMACS version:2020.1 and 2023.3

Hi,

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?

Thanks in advance,

Gianluca Dell’Orletta

2

There should not be any changed in the code that affect the physics. Have you tried continuing in the same way using 2020.1?

Note that MD is chaotic. You would need to run at least 3 repeats to be able to draw conclusions with some certainty.

3

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.

area_75-175_shifted
area_75-175_shifted792×612 5.85 KB

4

That looks like a serious artifact. How did you run the continuation simulations? GPUs? How many MPI ranks?

5

It looks like the artifacts appear very quickly. Could you check if the issue is also present with the latest version: 2025.2?

6

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.

Confronto_versioni_gromacs
Confronto_versioni_gromacs792×612 6.3 KB

7

Could you post your mdp parameters?

8

These are the parameters I used in my simulations.
integrator = md
nsteps = 50000000
dt = 0.002

; Output control
nstxout = 250000
nstvout = 250000
nstenergy = 5000
nstlog = 5000
nstxout-compressed = 5000

; Neighbor searching
cutoff-scheme = Verlet
nstlist = 20
ns_type = grid
rcoulomb = 1.2
rvdw = 1.2
vdw-modifier = Force-switch
rvdw-switch = 1.0

; Electrostatics
coulombtype = PME
pme-order = 4
fourierspacing = 0.16

; Temperature coupling
tcoupl = Nose-Hoover
tc-grps = System
tau-t = 1.0
ref-t = 350

; Pressure coupling
pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau-p = 5.0
compressibility = 4.5e-5
ref-p = 1.0

; GENERATE VELOCITIES FOR STARTUP RUN
gen_vel = no

; Constraints
constraints = h-bonds
constraint-algorithm = LINCS
lincs-order = 4
lincs-iter = 1

; Periodic boundary conditions
pbc = xyz
; Center of mass motion removal
comm-mode = Linear
comm-grps = System

9

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?

10

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).

area
area792×612 5.65 KB

Thank you for your help.
Best regards,
Gianluca Dell’Orletta