Fix 1205 gromacs2026 energy

[invemichele-peptone]

Support ENERGY in the GROMACS 2026 patch

Branch: fix-1205-gromacs2026-energy → base v2.10 · Stacked on: #1439, please merge that first
Related: #1205, #1384

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

Description

#1439 brings patches/gromacs-2026.0 into the tree, which is what makes a PLUMED-enabled GROMACS
2026 buildable at all. This PR adds one thing on top: the ENERGY collective variable and energy
biasing, which the native GROMACS PLUMED module does not support on any version from 2025 onwards.

On stock GROMACS 2025/2026 an input containing ENERGY runs to completion and produces 0.0 for
the whole trajectory with no warning: plumedforceprovider.cpp never issues setEnergy, and on
the PLUMED side ActionToPutData::wait() returns early when wasset == false. So OPES_EXPANDED
with ECV_MULTITHERMAL, the well-tempered ensemble, PT-WTE and energy reweighting are silently
wrong rather than unavailable. The GROMACS 2026 manual lists "interaction with GROMACS energy"
under the interface Limitations.

This branch starts from #1439 rather than v2.10 because it edits
patches/gromacs-2026.0.diff/.../plumedforceprovider.{cpp,h}, which #1439 introduces. It does not
re-do anything #1439 covers — the setMPIComm address-of fix, the GREX/multisim block,
setNumOMPthreads, PLUMED_LOG_FILE and the replex guard are taken as-is. The replica-exchange
guard is factored into a small checkReplicaExchangeBias() helper so the new deferred path
performs the same check; its behaviour is unchanged. The companion PR against v2.10 does the
same for patches/gromacs-2025.0 and also fixes the GROMACS 2024 regression (#1205), which has
the same visible symptom.

What blocks energy biasing in the ForceProvider model. Nobody can ask GROMACS to compute the
energy: stepWork.computeEnergy is only true on nstcalcenergy steps and a force provider cannot
influence the step workload. At force-provider time the energy does not exist yet:
calculateForces() is called from computeSpecialForces(), several hundred lines before
accumulatePotentialEnergies() fills the potential-energy term. And the force-rescaling trick is
not expressible: biasing the energy multiplies the entire force array by 1 - ∂V/∂E (PLUMED does
this in Energy::apply() via rescaleForces), while a provider only owns a separate additive
ForceWithVirial buffer and the total force is not assembled at that point.

Design. Two additions to IForceProvider, both defaulting to no-ops so no other module is
affected, forwarded by ForceProviders to every registered provider:

//! 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
still change, and ORs in GMX_FORCE_ENERGY | GMX_FORCE_VIRIAL. The answer has to be exact for
this step rather than inherited from the previous one, so
PlumedForceProvider::requestsPotentialEnergy() splits prepareCalc() in two: whether PLUMED
needs the energy depends on which actions are active, which is known after prepareDependencies(),
and that half needs only the step number (setStepLong → prepareDependencies →
isEnergyNeeded). This is the sequence the PLUMED action already uses in
src/generic/Plumed.cpp, so it needs no new PLUMED command and works with any 2.10 kernel;
prepareDependencies() only activates actions, with all atom and domain-decomposition handling
staying inside shareData(), so splitting there is safe. calculateForces() runs the second half
(shareData) once positions are available and falls back to the full prepareCalc() when the
early hook did not run for this step — that fallback 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 late hook after accumulatePotentialEnergies() and
postProcessForces(), where the energy 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 untouched. If a step that
never asked for the energy reaches the late hook needing it, PLUMED throws a NotImplementedError
pointing at the modular simulator and GMX_DISABLE_MODULAR_SIMULATOR=1 rather than biasing on zero.

Four files are new to patches/gromacs-2026.0.diff (iforceprovider.{h,cpp}, md.cpp,
sim_util.cpp), two are edits to files #1439 already carries (plumedforceprovider.{cpp,h}), and
the .config blurb gains a line. As #1439 explains, every patched file is stored twice in full, so
GitHub reports ~10.7k added lines for what is 76 changed lines of GROMACS source.

Verification

GROMACS 2026.3, mixed precision, 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.
GROMACS 2023.5 with the legacy patch is the reference, since it predates the whole problem.
plumed patch -p -e gromacs-2026.0 applies with no fuzz and no rejects, plumed patch -r restores
the tree byte-for-byte, and the patched tree builds with no errors or warnings in any of the six
changed translation units. ENERGY every step over 250 steps gives 251 / 251 finite non-zero
values, where unpatched the same input yields 0.0 on every step.

OPES multithermal, 300–350 K, PACE=25, OBSERVATION_STEPS=8, 1500 steps	2023.5 (ref)	2026.3 + this PR
temperatures chosen by OPES	10	11
steps with ENERGY == 0	0 / 1501	0 / 1501
DeltaF at 350 K after 3 ps (kJ/mol)	1435	1441

The temperature count is the sharpest check, because OPES derives the grid from the energy
fluctuations observed during the first 0.4 ps: with the energy stuck at zero it collapses to 2
temperatures and every DeltaF degenerates. Trajectories are not bit-identical across GROMACS
versions, so the residual spread is MD noise. Also checked: PRINT ARG=ene STRIDE=7 against
nstcalcenergy = 100, a stride sharing no factor with the energy period and precisely 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 it; and DISTANCE + RESTRAINT unchanged,
confirming the additive ForceWithVirial path is untouched.

GPU. Rebuilt with -DGMX_GPU=CUDA (NVIDIA GB10, one rank) and re-run with PME so the GPU has
real work. ENERGY every step is correct
under -nb cpu (reference), -nb gpu, -nb gpu -pme gpu, fully GPU-resident
-nb gpu -pme gpu -update gpu, and the last of those with GMX_CUDA_GRAPH=1: 121 / 121 in every
case, no zeros, no NaN. OPES multithermal gives 11–12 temperatures and DeltaF 1437–1438 across
CPU, -nb gpu and GPU-resident, matching the CPU baseline. The strongest check is against GROMACS
itself rather than another run: on every step written to the .edr the value PLUMED receives
matches GROMACS' own Potential term to ≤ 5 × 10⁻⁸ relative in all five configurations,
including fully GPU-resident stepping. Restarting across a non-nstcalcenergy step with
-update gpu is clean. Note that GROMACS only takes the CUDA-graph path when nothing computes
forces on the CPU, which PLUMED always does, so GMX_CUDA_GRAPH=1 exercises the surrounding
bookkeeping rather than graph capture itself.

MPI and domain decomposition. Rebuilt with -DGMX_MPI=ON (OpenMPI) and run on 1728 SPC waters
(5184 atoms) on 1, 2 and 4 ranks, DD grids 2 x 1 x 1 and 4 x 1 x 1. The deferred hook behaves:
ENERGY is finite and non-zero on all 1001 frames at every rank count, OPES multithermal runs
clean and DeltaF at 350 K agrees between 1 and 4 ranks to 1 kJ/mol, and the value PLUMED receives
matches GROMACS' own Potential term to ≤ 1.5 × 10⁻⁷ relative. This also exercises the
setMPIComm address-of fix that #1439 brings in, without which the run aborts before it starts.

Depends on the core fix in the companion v2.10 PR. Energy biasing under domain
decomposition also needs a one-function fix to ActionToPutData::getNumberOfForcesToRescale(),
which rescales the forces of the MD code over the global atom count instead of the atoms local
to the rank and so writes past the end of the buffer. It is a v2.10 core regression, present in
2.10.1 and in master, unrelated to this patch but hit by exactly the same inputs — without it
every OPES multithermal run above dies with double free or corruption on more than one rank,
on 2023.5 just as much as on 2026.3. It is not duplicated here, to keep this PR to the
patch; the MPI numbers above were obtained with it applied. Please merge the v2.10 PR first,
or say the word and I will carry the commit here instead.

Remaining gaps. The modular simulator is the one
unsupported path: it bypasses the md.cpp loop, so the early hook never runs and the energy is
never scheduled. It is the default for the md-vv integrator, where energy-dependent input now
fails with the NotImplementedError above; with GMX_DISABLE_MODULAR_SIMULATOR=1, md-vv works.
PLUMED input that does not use the energy is unaffected and still runs under the modular simulator.
Energy minimisation and mdrun -rerun do work. 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 need 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, alongside #1439 and
patches/gromacs-2025.0.diff.

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.

No regtest is added, for the reason #1439 gives: nothing under patches/ is compiled by PLUMED, so
the test suite cannot exercise it. It is consumed only when a user runs plumed patch against a
GROMACS source tree, at which point GROMACS compiles it. The GROMACS builds and MD runs above are
the substitute. This branch touches no core code, so the existing suite is unaffected.