Fix ENERGY CV gromacs2024/25 + fix ENERGY CV with multi MPI and negative float to double (v2.10, any MD engine)

[invemichele-peptone]

Fix ENERGY with GROMACS 2024, and support it with the native GROMACS 2025 interface

Branch: fix-1205-gromacs2024-energy-v2.10 → base v2.10 · Fixes: #1205 (and the root cause behind #1384)
Companions: the same 2024 fix against v2.9; the GROMACS 2026 counterpart stacked on #1439,
which needs part 4 below and so should be merged after this one

Disclaimer. This PR was entirely generated by Claude, supervised by @invemichele.

Description

Energy-based PLUMED input is broken or unavailable on every GROMACS newer than 2023:

• 2024 — ENERGY silently returns 0.0 except on nstcalcenergy steps (ENERGY bug, gromacs 2024 + plumed 2.9 #1205).
• 2025 — the native PLUMED module never calls setEnergy, so ENERGY stays 0.0 for the
  whole run with no error. The GROMACS manual lists "interaction with GROMACS energy" under the
  interface Limitations.
• Single-precision builds, every version — the core converts the scalars passed by the MD code
  from float to double by rounding them through std::log10, which is NaN for a negative
  argument. The energy is negative, so wherever it did arrive it arrived as NaN.
• Every version, more than one MPI rank — applying a force on ENERGY writes past the end of
  the force buffer of the MD code and corrupts the heap. This one is a v2.10 regression in the
  core and is not specific to GROMACS.

This PR fixes all four, so ENERGY, OPES_EXPANDED with ECV_MULTITHERMAL, the well-tempered
ensemble and energy reweighting behave on 2024 and 2025 as they do on 2023.

Scope. This branch carries three bug fixes that belong on a maintenance branch (the 2024
patch, identical to the v2.9 companion, and the DataPassingObject and domain decomposition
fixes) plus one new feature (ENERGY in the 2025 patch). If needed we can split it into smaller chunks.

Part 1 — GROMACS 2024: ENERGY is zero on most steps. The patch asks for the energy by
OR-ing bits into force_flags after isEnergyNeeded. In 2023 that worked because force_flags
went straight into the do_force() call below; in 2024 do_force() no longer takes flags and
the workload is snapshotted ~120 lines earlier by setupStepWorkload(), so the PLUMED line
became dead code. do_force() zeroes term[F_EPOT] every step and only refills it when
stepWork.computeEnergy is set, so PLUMED received 0.0 on 99 steps out of 100 with the default
nstcalcenergy = 100. The virial was lost the same way, so the force rescaling used for energy
biasing was wrong too. In patches/gromacs-2024.3.diff/.../md.cpp: carry the request forward
before setupStepWorkload(); rebuild runScheduleWork_->stepWork right after isEnergyNeeded
answers for the current step, which makes the fix exact rather than one step late; gmx_fatal if
the energy is still not scheduled; treat plumedNeedsEnergy like bCalcVir in the MD GPU graph
reset and reuse conditions.

Part 2 — the float to double conversion returns NaN for negative scalars.
DataPassingObjectTyped<float>::saveValueAsDouble rounds to 6 significant figures using
std::log10(bvalue), to keep single-precision values from acquiring meaningless digits when
widened. std::log10 is NaN for a negative argument and -inf for zero, which then gives
0/0, so any negative or zero scalar coming from a single-precision engine is turned into NaN.
The potential energy is essentially always negative, so ENERGY arrived as NaN and poisoned
every bias built on it. Fixed by taking the magnitude from std::fabs(bvalue) and leaving exact
zeros alone. Nothing here is specific to GROMACS or to ENERGY, but ENERGY is the only value
affected in practice: array data keeps the pointer of the MD code and is never rounded, timestep
and kBT take the same path but are strictly positive, and EXTRACV values are declared
MUTABLE, which sends them down the pointer branch instead.

Part 3 — GROMACS 2025: teaching the native interface about the energy. From 2025 the
interface is an IForceProvider rather than a source patch, and three things block energy
biasing: nobody can ask GROMACS to compute the energy, since a force provider cannot influence
stepWork; at calculateForces() time the energy does not exist yet, as it is called from
computeSpecialForces() several hundred lines before accumulatePotentialEnergies(); and the
force-rescaling trick is not expressible, because biasing the energy multiplies the entire
force array by 1 - ∂V/∂E (Energy::apply() → rescaleForces) while a provider only owns a
separate additive ForceWithVirial buffer. Two hooks are added to IForceProvider, both
defaulting to no-ops so no other module is affected:

//! Whether this provider needs the MD potential energy on `step`.
virtual bool requestsPotentialEnergy(int64_t /*step*/) { return false; }

//! Called once the potential energy is accumulated and the force buffer is complete.
virtual void applyAfterPotentialEnergy(bool /*energyWasComputed*/, real* /*potentialEnergy*/,
                                       ArrayRef<RVec> /*force*/, tensor /*virial*/) {}

md.cpp calls the first immediately before setupStepWorkload(), the last moment the workload
can change. The answer must be exact for this step rather than inherited from the previous one,
so requestsPotentialEnergy() splits prepareCalc() in two: whether PLUMED needs the energy is
known after prepareDependencies(), and that half needs only the step number
(setStepLong → prepareDependencies → isEnergyNeeded, the sequence src/generic/Plumed.cpp
already uses, so no new PLUMED command is required). calculateForces() runs the second half
(shareData) once positions are available, falling back to the full prepareCalc() when the
early hook did not run for this step — which is what makes energy minimisation and mdrun -rerun
work, since neither goes through the md.cpp loop but both always compute the energy.
sim_util.cpp calls the second hook after accumulatePotentialEnergies() and postProcessForces(),
where term[F_EPOT] is valid and the total force buffer is complete; there PLUMED does
setEnergy / setForces on the total force / setVirial / performCalc, with the same 2 ×
/ 0.5 × virial convention and virial replacement as the legacy 2023 and 2024 patches. When
the energy is not needed the existing additive path is unchanged.

Part 4 — energy biasing overran the force buffer under domain decomposition. Applying a force
on ENERGY makes PLUMED rescale all of the forces of the MD code by 1 - ∂V/∂E
(Energy::apply() → ActionToPutData::rescaleForces). The loop length comes from
getNumberOfForcesToRescale(), which returned copyOutput(0)->getNumberOfValues() — the total
number of atoms — while under domain decomposition the buffer of the MD code only holds the atoms
local to the rank. On 4 ranks PLUMED therefore wrote about four times past the end of that buffer,
which shows up as double free or corruption or a segfault the moment the bias switches on. The
domain-aware branch that returns the local count exists but is unreachable: it is guarded by
getName()!="ENERGY", whereas the values actually being rescaled are the posx/posy/posz
objects created as PUT FROM_DOMAINS, and it would have hit its own
plumed_assert(getDependencies().size()==1) if it had ever been entered. The fix keys the branch
off the existing from_domains flag instead. v2.9 is not affected — it rescales over gatindex,
i.e. the local atoms — so this is a v2.10 regression introduced with the new data-passing
architecture, and it is present in 2.10.1 and in current master as well.

This is the same combination as #1205: OPES multithermal and the well-tempered ensemble are
precisely the inputs that put a force on ENERGY, so anyone running them on more than one rank
was hitting silent memory corruption.

Verification

216 SPC waters (648 atoms), OPLS-AA, 0.8 nm cut-offs, dt = 2 fs, v-rescale at 300 K, LINCS on
h-bonds, nstcalcenergy = 100 unless stated, mixed precision. GROMACS 2023.5 with the legacy
patch is the reference throughout, since it predates the regression and works as shipped.

	2023.5 (reference)	2025.0 + this PR
ENERGY every step, 250 steps	251 / 251	251 / 251
temperatures chosen by OPES	10	10
steps with ENERGY == 0	0 / 1501	0 / 1501
bias switches on at	0.402 ps	0.402 ps
∂V/∂E non-zero after observation	1300 / 1300	1300 / 1300
DeltaF at 350 K after 3 ps (kJ/mol)	1435	1442

Unpatched 2024.3 manages 4 / 251, exactly the nstcalcenergy steps. The temperature count is the
sharpest check, since OPES sizes that grid from the energy fluctuations seen during observation, so
a zero or intermittent energy collapses it to 2 temperatures; a control at nstcalcenergy = 1
picks the same 10. Trajectories are not bit-identical across GROMACS versions, so the residual
DeltaF spread is MD noise. Also checked: PRINT ARG=ene STRIDE=7 against nstcalcenergy = 100,
a stride sharing no factor with the energy period and the case an inherited, one-step-late request
cannot serve; restart via -cpi from step 150, not a multiple of nstcalcenergy, with the energy
continuous across the restart; DISTANCE + RESTRAINT unchanged, confirming the additive path is
untouched; and plumed patch -p -e gromacs-2025.0 applying to a pristine tarball with no fuzz and
no rejects, with plumed patch -r restoring the tree byte-for-byte.

Part 4 needs a bigger box, so that one uses 1728 SPC waters (5184 atoms) with
OpenMPI, -DGMX_MPI=ON, on 1, 2 and 4 ranks (DD grids 2 x 1 x 1 and 4 x 1 x 1). Before the
fix, every run with a force on ENERGY died — double free or corruption on 2023.5 and on 2026.3
alike, in NVT and NPT, always at the step the bias switched on; ENERGY without a force, and a
RESTRAINT on a DISTANCE, both survived, which is what isolates the rescaling. After the fix:

check	result
OPES multithermal, 1 vs 4 ranks (2023.5 NVT / NPT, 2026.3 NPT)	clean; DeltaF at 350 K agrees to 1 kJ/mol
biased ÷ unbiased force, per atom, single force evaluation	uniform over all 5184 atoms, = 1 - ∂V/∂E, spread 2 × 10⁻⁵, same at 1 / 2 / 4 ranks
ENERGY vs .edr Potential, DispCorr = EnerPres	≤ 0.015 kJ/mol, i.e. ≤ 1.5 × 10⁻⁷ relative, at every rank count
regtest/basic and regtest/opes	identical results with and without the fix

The force-ratio row is the direct check: because a force on the energy only ever multiplies the
existing forces by a scalar, every atom on every rank must show the same ratio, and it must equal
1 - ∂V/∂E. It does. PLUMED 2.9 built with the same MPI stack was run through the same test and
is clean at 1, 2 and 4 ranks, confirming the regression is 2.10-only.

Also: the long tail correction warning in the docs is stale (#567). While checking what the
energy PLUMED receives actually contains, it turned out that the \bug note on ENERGY — "does
not include long tail corrections … GROMACS DispCorr Ener" — no longer holds for GROMACS. The
correction is computed inside do_force() and folded into F_EPOT by sum_epot(), both before
PLUMED reads the energy. With DispCorr = EnerPres, ENERGY matches the Potential term in the
GROMACS energy file to 0.005 kJ/mol on 2023.5, 2024.3 and 2026.3, while the correction itself
is −154 kJ/mol. Under NPT, where the volume varies by up to 29% and the correction with it by
40 kJ/mol, ENERGY still tracks GROMACS to 0.005 kJ/mol over 401 frames — so the correction is
not merely present but correctly updated. The note is reworded rather than deleted, and kept short
and engine-agnostic: codes other than GROMACS, such as LAMMPS with pair_modify tail yes, were not
tested here, and where the mismatch is real the advice to reweight with the energy of the engine
still stands.

Remaining gaps. GPU was tested through the companion 2026 patch, which carries the identical
design, since GROMACS 2025 does not build against the CUDA toolkit available here: -nb gpu,
-nb gpu -pme gpu, fully GPU-resident -update gpu and GMX_CUDA_GRAPH=1 all give correct
energies, matching GROMACS' own .edr Potential term to ≤ 5 × 10⁻⁸ relative, and OPES
multithermal reproduces the CPU temperature grid and DeltaF. The modular simulator is the one
unsupported path: it bypasses the md.cpp loop, so
the early hook never runs. It is the default for md-vv, where energy-dependent input now
fails with a NotImplementedError naming it and pointing at GMX_DISABLE_MODULAR_SIMULATOR=1;
with that set, md-vv works. Input that does not use the energy is unaffected. Under MTS the
total force is taken from the combined buffer only on slow steps. The two hooks are the minimal
shape of what would have to go upstream; GROMACS issue
#4939 does not list energy biasing among the
requirements for the PLUMED interface, which is probably why it was never designed in.

Target release

I would like my code to appear in release 2.10.

Type of contribution

• changes to code or doc authored by PLUMED developers, or additions of code in the core or within the default modules
• changes to a module not authored by you
• new module contribution or edit of a module authored by you

Copyright

• I agree to transfer the copyright of the code I have written to the PLUMED developers or to the author of the code I am modifying.

Tests

• I added a new regtest or modified an existing regtest to validate my changes.
• I verified that all regtests are passed successfully on GitHub Actions.

Parts 1 and 3 live entirely under patches/, which PLUMED does not compile, so the test suite
cannot exercise them; the GROMACS builds and MD runs above are the substitute. Part 2 is core
code and is covered by existing regtests, since it is on the path every float-precision engine
uses to pass a scalar. Part 4 is core code too, but no regtest can reach it today: the failure
needs gatindex.size() < natoms, and with no MPI communicator DomainDecomposition::shareAll()
takes the branch that assumes every atom is local, so a serial fake-MD harness cannot set it up —
and the regtest framework has no multi-rank test type. It would need a small MPI harness; happy to
add one, here or separately, if you want the infrastructure. regtest/basic and regtest/opes
were run with and without the fix and give identical results (the same two pre-existing failures,
rt-average and rt-multi-1, in both).