Feature force groups

corelib/src/libs/SireIO/amberprm.cpp

@@ -5396,6 +5396,63 @@ System AmberPrm::startSystem(const PropertyMap &map) const
         bar.success();
     }
 
+    if (not in_expected_order)
+    {
+        // Collect all residues in loaded (mol-index) order and check for
+        // residue-number inversions. These arise when covalent bonds cross
+        // AMBER's molecule boundaries (e.g. a metal ion bonded to a small
+        // molecule), forcing Sire to group the bonded atoms into the same
+        // molecule and thereby placing some residues out of their original
+        // residue-number sequence.
+        struct ResEntry
+        {
+            int resnum;
+            QString resname;
+        };
+
+        QVector<ResEntry> all_res;
+
+        for (const auto &mol : mols)
+        {
+            const auto &molinfo = mol.info();
+            const int nres = molinfo.nResidues();
+
+            for (int j = 0; j < nres; ++j)
+            {
+                all_res.append({molinfo.number(ResIdx(j)).value(),
+                                molinfo.name(ResIdx(j)).value()});
+            }
+        }
+
+        QStringList out_of_order;
+        int prev_resnum = 0;
+
+        for (const auto &res : all_res)
+        {
+            if (res.resnum < prev_resnum)
+            {
+                out_of_order.append(
+                    QString("%1(%2)").arg(res.resname).arg(res.resnum));
+            }
+
+            prev_resnum = res.resnum;
+        }
+
+        if (not out_of_order.isEmpty())
+        {
+            qWarning().noquote() << QObject::tr(
+                "WARNING: One or more residues have been reordered relative to the "
+                "original topology when loading this file. This happens when covalent "
+                "bonds cross AMBER's molecule boundaries (e.g. a metal ion bonded to a "
+                "small molecule ligand), which forces Sire to group the bonded atoms "
+                "into the same molecule. The following residues appear out of "
+                "residue-number order in the loaded system: %1. "
+                "To avoid ordering issues, access residues by name rather than by "
+                "position (e.g. mols.residues()[\"resname ML1\"] instead of mols.residues()[1]).")
+                .arg(out_of_order.join(", "));
+        }
+    }
+
     MoleculeGroup molgroup("all");
 
     for (auto mol : mols)

corelib/src/libs/SireMM/bondrestraints.cpp

@@ -188,6 +188,7 @@ BondRestraint &BondRestraint::operator=(const BondRestraint &other)
 {
     if (this != &other)
     {
+        Property::operator=(other);
         atms0 = other.atms0;
         atms1 = other.atms1;
         _k = other._k;

corelib/src/libs/SireMM/inversebondrestraints.cpp

@@ -188,6 +188,7 @@ InverseBondRestraint &InverseBondRestraint::operator=(const InverseBondRestraint
 {
     if (this != &other)
     {
+        Property::operator=(other);
         atms0 = other.atms0;
         atms1 = other.atms1;
         _k = other._k;

corelib/src/libs/SireMM/morsepotentialrestraints.cpp

@@ -190,6 +190,7 @@ MorsePotentialRestraint &MorsePotentialRestraint::operator=(const MorsePotential
 {
     if (this != &other)
     {
+        Property::operator=(other);
         atms0 = other.atms0;
         atms1 = other.atms1;
         _k = other._k;

corelib/src/libs/SireMM/positionalrestraints.cpp

@@ -146,6 +146,7 @@ PositionalRestraint &PositionalRestraint::operator=(const PositionalRestraint &o
 {
     if (this != &other)
     {
+        Property::operator=(other);
         atms = other.atms;
         pos = other.pos;
         _k = other._k;

doc/source/changelog.rst

@@ -36,6 +36,15 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 
 * Add support for 4- and 5-point water models in the OpenMM conversion layer.
 
+* Add per-force-group energy caching to the OpenMM integration layer. Named
+  forces (``clj``, ``bond``, ``angle``, ``torsion``, ``cmap``, ghost forces,
+  restraints, etc.) are each assigned a unique OpenMM force group when the
+  system is built. :meth:`~sire.mol.Dynamics.get_potential_energy` now
+  re-evaluates only the groups whose parameters changed since the last lambda
+  update, using cached values for all others. Call
+  :meth:`~sire.mol.Dynamics.clear_energy_cache` to force a full re-evaluation
+  (e.g. after a replica-exchange position swap).
+
 * Add functionality for coupling one lambda lever to another.
 
 * Added support for Direct Morse Replacement (DMR) feature in ``sire.restraints.morse_potential``
@@ -45,6 +54,8 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 
 * Store OpenMM state at start of a dynamics run to use for crash recovery.
 
+* Print warning when ``sire.legacy.IO.AmberPrm`` parser re-orders residues due to covalent bonds between molecules.
+
 `2025.4.0 <https://github.com/openbiosim/sire/compare/2025.3.0...2025.4.0>`__ - February 2026
 ---------------------------------------------------------------------------------------------
 

src/sire/mol/_dynamics.py

@@ -483,6 +483,10 @@ def _exit_dynamics_block(
                 nrg_sim_lambda_value = nrg
 
                 if lambda_windows is not None:
+                    # Positions have just changed (dynamics completed), so
+                    # invalidate all cached per-group energies before the scan.
+                    self._omm_mols.clear_energy_cache()
+
                     # get the index of the simulation lambda value in the
                     # lambda windows list
                     try:
@@ -542,6 +546,10 @@ def _exit_dynamics_block(
             self._nrgs = nrgs
             self._nrgs_array = nrgs_array
 
+            # Repex synchronisation point: a peer replica may push new
+            # positions into this context, so the cache must be invalidated.
+            self._omm_mols.clear_energy_cache()
+
             # update the interpolation lambda value
             if self._is_interpolate:
                 if delta_lambda:
@@ -853,14 +861,11 @@ def current_potential_energy(self):
         if self.is_null():
             return 0
         else:
-            from openmm.unit import kilocalorie_per_mole as _omm_kcal_mol
-            from ..units import kcal_per_mol as _sire_kcal_mol
+            return self._omm_mols.get_potential_energy(to_sire_units=True)
 
-            state = self._get_current_state()
-
-            nrg = state.getPotentialEnergy()
-
-            return nrg.value_in_unit(_omm_kcal_mol) * _sire_kcal_mol
+    def clear_energy_cache(self):
+        if not self.is_null():
+            self._omm_mols.clear_energy_cache()
 
     def current_kinetic_energy(self):
         if self.is_null():
@@ -989,6 +994,9 @@ def run_minimisation(
             timeout=timeout.to(second),
         )
 
+        # Positions changed during minimisation; invalidate the energy cache.
+        self._omm_mols.clear_energy_cache()
+
     def _rebuild_and_minimise(self):
         if self.is_null():
             return
@@ -1030,38 +1038,28 @@ def _rebuild_and_minimise(self):
         if self._save_crash_report:
             import openmm
             import numpy as np
-            from copy import deepcopy
             from uuid import uuid4
 
             # Create a unique identifier for this crash report.
             crash_id = str(uuid4())[:8]
 
-            # Get the current context and system.
             context = self._omm_mols
-            system = deepcopy(context.getSystem())
-
-            # Add each force to a unique group.
-            for i, f in enumerate(system.getForces()):
-                f.setForceGroup(i)
-
-            # Create a new context.
-            new_context = openmm.Context(system, deepcopy(context.getIntegrator()))
-            new_context.setPositions(context.getState(getPositions=True).getPositions())
 
-            # Write the  energies for each force group.
+            # Write per-force-group energies using the groups already assigned
+            # by sire_to_openmm_system.
             with open(f"crash_{crash_id}.log", "w") as f:
                 f.write(f"Current lambda: {str(self.get_lambda())}\n")
-                for i, force in enumerate(system.getForces()):
-                    state = new_context.getState(getEnergy=True, groups={i})
-                    f.write(f"{force.getName()}, {state.getPotentialEnergy()}\n")
+                for name, grp in context._force_group_map.items():
+                    state = context.getState(getEnergy=True, groups=(1 << grp))
+                    f.write(f"{name}, {state.getPotentialEnergy()}\n")
 
             # Save the serialised system.
             with open(f"system_{crash_id}.xml", "w") as f:
-                f.write(openmm.XmlSerializer.serialize(system))
+                f.write(openmm.XmlSerializer.serialize(context.getSystem()))
 
             # Save the positions.
             positions = (
-                new_context.getState(getPositions=True).getPositions(asNumpy=True)
+                context.getState(getPositions=True).getPositions(asNumpy=True)
                 / openmm.unit.nanometer
             )
             np.savetxt(f"positions_{crash_id}.txt", positions)
@@ -2252,6 +2250,14 @@ def _current_energy_array(self):
         """
         return self._d._current_energy_array()
 
+    def clear_energy_cache(self):
+        """
+        Invalidate the per-force-group energy cache. Call this whenever
+        positions have been changed externally (e.g. after a replica-exchange
+        swap) so that the next energy evaluation fully re-computes all groups.
+        """
+        self._d.clear_energy_cache()
+
     def to_xml(self, f=None):
         """
         Save the current state of the dynamics to XML.