Intermittent LatticeReduction::sort assertion m[0]<=m[1]*onePlusEpsilon during well-tempered metadynamics (GROMACS 2025.4 runtime kernel, v2.10.0)

[PricklyE]

Summary

PLUMED 2.10.0, loaded as runtime kernel by GROMACS 2025.4, occasionally aborts with an internal assertion in LatticeReduction::sort during a well-tempered 2D metadynamics run on an explicit-solvent polymer-drug system. The crash is intermittent — restarts from the same checkpoint sometimes reproduce the crash within a few ns, sometimes run cleanly for tens of thousands of steps past the last crash point. No equivalent fix commit exists ahead of v2.10.0 on either v2.10 maintenance branch or master (verified via git log v2.10.0..origin/v2.10 -- src/tools/LatticeReduction.cpp → empty).

Environment

• PLUMED: 2.10.0 (tag 0f2be2d617d7b8d577f1856c9172759543f5c9f0 Release v2.10.0)
• GROMACS: 2025.4, mixed precision, CUDA 13.1, linked against PLUMED via runtime kernel loader (PLUMED_KERNEL env var; PLUMED is not statically patched into GROMACS)
• Hardware: NVIDIA RTX PRO 6000 Blackwell Max-Q Workstation Edition, AMD Threadripper 7960X, Ubuntu 24.04
• mdrun flags: -ntmpi 1 -ntomp 16 -nb gpu -pme gpu -bonded gpu -update cpu -pin on -pinoffset 16 -gpu_id 1 -noappend
• PLUMED CV stack: DISTANCE between two COM groups + COORDINATION with NLIST NL_CUTOFF=1.8 + GYRATION, biased by 2D well-tempered METAD
• Integrator: 1 fs timestep, NPT with standard GROMACS pressure coupling (Parrinello-Rahman or C-rescale — system mdp standard)

Assertion message (verbatim)

-------------------------------------------------------
Program:     gmx mdrun, version 2025.4
Source file: src/gromacs/applied_forces/plumed/plumedforceprovider.cpp (line 212)
Function:    virtual void gmx::PlumedForceProvider::calculateForces(
               const gmx::ForceProviderInput&, gmx::ForceProviderOutput*)

Internal error (bug):
An error occurred while PLUMED was calculating the forces
:
(tools/LatticeReduction.cpp:43) static void
PLMD::LatticeReduction::sort(PLMD::Vector*)
+++ assertion failed: m[0]<=m[1]*onePlusEpsilon
-------------------------------------------------------

Observation history

Attempt	Kernel	Start step	End step / outcome	Progress (ns)
1 — Apr 18	2.10.0	17 252 300	crashed at assertion (step not logged in summary)	~4-5 ns
2 — Apr 18 (retry)	2.10.0	20 661 100	crashed at step 22 910 000 at the same assertion	2.25 ns
3 — Apr 20 (15 min diagnostic run on sandbox copy of same checkpoint)	2.10.0	22 902 100	clean exit at -maxh 0.25, reached step 22 936 500	0.034 ns, no crash

Key observation: attempts 1 and 2 crashed; attempt 3, resuming from a checkpoint just 8 ps earlier than attempt 2's crash point, ran past that point cleanly. The bug is not deterministically triggered by a particular trajectory state — GPU GROMACS + PLUMED restarts are not bit-exact (non-associative FP reduction on GPU; PLUMED stochastic elements), so each restart explores a slightly different numerical trajectory.

Why this is probably numerical

The assertion m[0]<=m[1]*onePlusEpsilon is in LatticeReduction::sort(), called during lattice-vector magnitude sorting for minimum-image calculations. Given the box evolves under NPT with moderate anisotropy, two lattice-vector magnitudes occasionally become degenerate within the onePlusEpsilon tolerance (current value ≈ 1 + 1e-6), making the sort's ordering invariant unsatisfied.

Reproducer

Available on request: metad.tpr (45 MB), plumed.dat, metad.cpt (38 MB), start.gro (112 MB), system.top, index.ndx. System: 6285 atoms polymer (10-chain cationic amphiphile), 47 atoms drug (curcumin), ~19 357 solute atoms total in explicit TIP3P water, triclinic box ~6×6×6 nm, 1 fs timestep.

Crash is intermittent; tends to surface within 2-6 ns of restart. Current MTBF estimate (n=2 crashes over 8 ns progress): ~4 ns between crashes, but sample is small.

Suggested investigation

1. Is onePlusEpsilon safely chosen for single-precision float paths? (GROMACS 2025.4 used mixed precision; PLUMED kernel uses double internally AFAIK, but the vectors may originate from GROMACS single floats.)
2. Would std::stable_sort with a relaxed comparator (tolerance, or nominal ordering on tie) be a safer fallback than the assertion?
3. Is there a precision gap between single-precision GROMACS box vectors passed into the PLUMED kernel vs. what the kernel expects?

Workarounds already ruled out

• NOPBC on DISTANCE / COORDINATION: scientifically incorrect for inter-molecular CVs where the drug crosses box boundaries.
• WHOLEMOLECULES: keeps molecules intact but does not fix inter-molecular minimum-image calculation.
• Downgrade to v2.9: shares the same LatticeReduction.cpp code (78c27bc67 Fix lattice reduction from Oct 2020 is already in both trees). Unlikely to help.
• Pull v2.10 maintenance branch HEAD: verified — no commits touch LatticeReduction.cpp since v2.10.0.

Current workaround in use

Auto-restart wrapper (run_metad_restart.sh) that retries on any non-zero gmx mdrun exit from the same checkpoint. Tolerable for this 100 ns run (expected ~20-40 restarts worst-case), but not a long-term solution for ensemble or production workflows.

[GiovanniBussi]

Thanks!

Yes this is a numerical issue. It would be great to have a simple input to reproduce it. Would it be possible to have a box shape that triggers the bug?

Alternatively: can you rerun with a higher tolerance (e.g. 1 + 1e-5) and see if it works?

I am really puzzled by this. In principle the check if not even necessary, but I suspect that the bug is trigger only if you have crazily tilted boxes.

[PricklyE]

Thanks for the quick reply Giovanni.

I dug into the box geometry and I think the tilted-box hypothesis can be ruled out for this case. A couple of data points:

Box is cubic, not tilted. The .mdp uses pbc = xyz with pcoupl = Parrinello-Rahman and pcoupltype = isotropic. Extracting Box-X / Box-Y / Box-Z from the .edr file (every 10 ps over 22 ns of simulated time) shows Box-X == Box-Y == Box-Z to six decimal places at every recorded frame, with RMS deviation of ~0.003 nm on a 25.23 nm edge. There are no off-diagonal Box-XY/XZ/YZ terms in the .edr (GROMACS only writes them for triclinic boxes). Example window immediately before the crash (step 22 910 000):

    time (ps)    Box-X        Box-Y        Box-Z
    22800.00    25.227020    25.227020    25.227020
    22810.00    25.220860    25.220860    25.220860
    22820.00    25.229870    25.229870    25.229870
    ...
    22890.00    25.227095    25.227095    25.227095
    22900.00    25.224274    25.224274    25.224274
    22910.00    25.233963    25.233963    25.233963   ← crash between this and next frame

So the three lattice vectors are supposed to be geometrically identical. My read is that the assertion is firing on normal GPU single-precision FP rounding in the Box-X/Y/Z values — the sort's strict m[0] <= m[1] * (1 + 1e-6) invariant gets violated by ~1 ppm noise in whichever direction the reduction happens to land that frame. That also fits the intermittency: GROMACS + PLUMED restarts aren't bit-exact on GPU (non-associative FP reduction), so each restart gives a different rounding pattern past the crash point.

This also fits your remark that "the check is not even necessary" — for an isotropic-cubic box the ordering of three equal-magnitude vectors has no geometric meaning, so any permutation should be valid.

I'll queue up the higher-tolerance test (1 + 1e-5) against the crashing checkpoint and report back. I'm also now confident a minimal reproducer is feasible — a small water box (~10 nm cube) under pcoupltype = isotropic + mixed-precision CUDA GROMACS + any 2D METAD hook should reproduce the assertion in a few ns. I can put that together and share it if useful.

Environment detail that may be relevant: GROMACS 2025.4 mixed precision, CUDA 13.1, RTX PRO 6000 Blackwell Max-Q (sm_120). PLUMED 2.10.0 loaded via PLUMED_KERNEL, not statically patched.

[GiovanniBussi]

I am cross checking the code now. With the current implementation, epsilon is 1e-14, not 1e-5 (unless you changed it manually).

In the past (see #524 ) I found an issue with this part of the code using a buggy compiler (see more in the discussion following #524). On the other hand, the current setting (which is very strict) seems to pass all tests and I do not recall recent reports similar to yours.

I think it would be useful to replace that specific assert with an if followed by some information about the exact vectors. I will do that now

[GiovanniBussi]

In the end, if you have a setting such that epsilon = 1e-10 or so will work, I would be happy to replace it. However, I still think that the requirement to use a very large epsilon might be just a symptom of a compiler bug or another serious numerical issue, so I would be a bit skeptical about changing the tolerance.

Anyway, you could try your example with the modified code I just pushed on v2.9/v2.10/master branches and see what's the message now.

[PricklyE]

Thanks Giovanni — and you're right about the 1e-14, I misread the source.

I've pulled your patch (v2.10 FETCH_HEAD including 1d8519067) and rebuilt the kernel natively on the production machine. The pmf_ccm_2d metadynamics has resumed from step 29 907 100 with the patched PLUMED loaded, and the auto-restart wrapper will keep it going until the next would-be-crash. Given the ~2-4 ns MTBF on this system (~9.5 ns/day throughput), the first diagnostic dump should arrive within 5-10 hours. I'll post the vectors here as soon as they land.

On the tolerance question: my system is mixed-precision GROMACS (single on GPU, double on some CPU paths). Box-X/Y/Z arrive at PLUMED with ~1e-7 single-precision noise, and pcoupltype = isotropic forces the three magnitudes to be supposed to be identical — so the comparison is purely against FP rounding, not geometric content. I'll wait for the actual values before speculating further on whether 1e-10 is adequate or there's something upstream worth looking at.

Separately — a minimal 10 nm pure-water reproducer (same mdp, amber99sb-ildn + TIP3P, 98 k atoms, 2D METAD on DISTANCE+COORDINATION) ran cleanly for 9.6 ns without hitting the assertion, so the trigger isn't purely "cubic box + isotropic NPT + GPU FP". Size (production is ~500 k atoms) or chemistry (production has a cationic polymer + drug) is apparently also a factor. Happy to scale the reproducer up to ~500 k atoms if that's useful, but I suspect the diagnostic from the patched real run will point more directly at the issue.

[PricklyE]

Quick update — the patched kernel has been running cleanly for about 6 hours now, through step 32,540,000 on the production system. That's ~2.3 ns past the resumed checkpoint, which is at the optimistic end of the 2–4 ns MTBF window the pre-patch runs showed. No assertion has fired yet and the new diagnostic log line has not been emitted.

One interpretation is that the trajectory region around the original crash points (steps 22–30 M) was numerically pathological in some sense — we have now "walked past" it and the current box-vector rounding pattern may simply not trigger the invariant. The run is holding benchmark throughput (≈ 9.8 ns/day on an RTX PRO 6000 Max-Q), so the patched code path is not costing anything when the if-branch doesn't trip.

An auto-poll is watching the logs every 5 min; I'll post the vectors the moment the diagnostic fires. If the run reaches its 100 ns target (~8 more days at current rate) without the assertion ever triggering, that itself is a data point worth recording and I'll follow up then.