Add support for auto-generating Morse restraints for ring-breaking

CHANGELOG.md

@@ -19,6 +19,7 @@ Changelog
 * Skip minimisation on restart [#191](https://github.com/OpenBioSim/somd2/pull/191).
 * Pre-equilibrate the water with GCMC moves before minimising in the regular `Runner`, making it consistent with the `RepexRunner`, which already did so to stop the geometry relaxing into a dry pocket [#191](https://github.com/OpenBioSim/somd2/pull/191).
 * Add a `precision` option for GPU platforms, defaulting to `mixed` [#191](https://github.com/OpenBioSim/somd2/pull/191).
+* Add support for generating Morse restraints for ring-breaking perturbations.
 
 [2026.1.0](https://github.com/openbiosim/somd2/compare/2025.1.0...2026.1.0) - Jun 2026
 --------------------------------------------------------------------------------------

README.md

@@ -289,6 +289,87 @@ geometry. To override this for all groups:
 somd2 perturbable_system.bss --terminal-flip-frequency "1 ps" --terminal-flip-angle "180 degrees"
 ```
 
+## Lambda schedules
+
+The way that the perturbation is applied across the lambda coordinate is
+controlled by the `--lambda-schedule` option, which defaults to
+`standard_morph`. The available schedules are:
+
+| Schedule | Description |
+| --- | --- |
+| `standard_morph` | Linear interpolation between the two end states. |
+| `charge_scaled_morph` | As above, but with charges scaled at intermediate lambda values. |
+| `annihilate` | Absolute binding free energies, removing all non-bonded interactions. |
+| `decouple` | Absolute binding free energies, removing only intermolecular interactions. |
+| `ring_break_morph` | Ring-breaking perturbations. |
+| `reverse_ring_break_morph` | Ring-making perturbations, i.e. the reverse of the above. |
+
+For the `annihilate`, `decouple`, and ring-breaking schedules, appropriate
+restraints can be generated automatically. See the sections below.
+
+## Absolute binding free energies
+
+Absolute binding free energy (ABFE) calculations are supported using the
+`annihilate` and `decouple` lambda schedules. Both first discharge the ligand,
+then remove its Lennard-Jones interactions: `annihilate` removes all non-bonded
+interactions, including those within the ligand, whereas `decouple` retains the
+intramolecular terms.
+
+```
+somd2 perturbable_system.bss --lambda-schedule decouple
+```
+
+The ligand must be restrained within the binding site. If no restraints are
+passed, a Boresch restraint is generated automatically for the bound leg, i.e.
+when the system contains both a protein and water. This is done by minimising
+the system, running a short trajectory at lambda = 0, then choosing the anchor
+atoms and force constants from it. The length of this trajectory and the
+frequency at which frames are saved can be controlled with the
+`--restraint-search-time` and `--restraint-search-frequency` options. By
+default the receptor anchor atoms are chosen from the protein backbone; use
+`--restraint-search-receptor-selection` to pass a `Sire` selection string
+instead.
+
+The restraint is written to `abfe_restraint.s3` in the output directory and is
+reloaded on restart, since the accumulated free energy corresponds to that
+particular restraint. The standard state correction is logged and written to
+the metadata of the energy trajectory, so analysis code can apply it without
+needing to scan the log.
+
+> [!NOTE]
+> The Beutler soft-core form, enabled with `--softcore-form beutler`, is only
+> supported with the ABFE schedules, or a custom schedule.
+
+## Ring-breaking perturbations
+
+Perturbations that break (or form) a ring are supported using the
+`ring_break_morph` schedule, or `reverse_ring_break_morph` for the ring-making
+direction.
+
+```
+somd2 perturbable_system.bss --lambda-schedule ring_break_morph
+```
+
+These perturbations require a pair of Morse restraints on the atoms of the bond
+that is broken. If no restraints are passed, both are generated automatically.
+A "hard" Morse potential replaces the harmonic bond, inheriting its force
+constant and equilibrium length, and is switched off as a weaker "soft" Morse
+restraint holds the fragment in place. Their well depths and the force constant
+of the soft restraint can be controlled with the `--morse-hard-well-depth`,
+`--morse-soft-well-depth`, and `--morse-soft-force-constant` options.
+
+Unlike the ABFE restraints, these are regenerated on each run rather than being
+cached, since they are derived from the bond parameters alone and are therefore
+identical every time.
+
+> [!NOTE]
+> The defaults are a reasonable starting point, but ring-breaking
+> perturbations are demanding. A non-uniform spacing of lambda values, set with
+> `--lambda-values`, is typically needed to obtain good overlap around the point
+> at which the bond is broken. The
+> [alchemate](https://github.com/akalpokas/alchemate) package provides
+> workflows for iteratively optimising the lambda schedule.
+
 ## Debugging with energy components
 
 To help diagnose simulation instabilities, `SOMD2` can record the potential

src/somd2/config/_config.py

@@ -180,6 +180,9 @@ def __init__(
         restraint_search_time="1 ns",
         restraint_search_frequency="10 ps",
         restraint_search_receptor_selection=None,
+        morse_hard_well_depth="150 kcal mol-1",
+        morse_soft_well_depth="50 kcal mol-1",
+        morse_soft_force_constant="125 kcal mol-1 A-2",
     ):
         """
         Constructor.
@@ -622,6 +625,20 @@ def __init__(
             Sire selection string for receptor anchor atom candidates used
             during automatic Boresch restraint generation. If None, the default
             backbone selection is used (CA, C, N atoms in non-water molecules).
+
+        morse_hard_well_depth: str
+            The well depth of the "hard" Morse potential that replaces the
+            broken bond when auto-generating restraints for a ring-breaking
+            simulation.
+
+        morse_soft_well_depth: str
+            The well depth of the "soft" Morse restraint that holds the broken
+            fragment in place when auto-generating restraints for a
+            ring-breaking simulation.
+
+        morse_soft_force_constant: str
+            The force constant of the "soft" Morse restraint used when
+            auto-generating restraints for a ring-breaking simulation.
         """
 
         # Setup logger before doing anything else
@@ -717,6 +734,9 @@ def __init__(
         self.restraint_search_time = restraint_search_time
         self.restraint_search_frequency = restraint_search_frequency
         self.restraint_search_receptor_selection = restraint_search_receptor_selection
+        self.morse_hard_well_depth = morse_hard_well_depth
+        self.morse_soft_well_depth = morse_soft_well_depth
+        self.morse_soft_force_constant = morse_soft_force_constant
         self.write_config = write_config
         self.overwrite = overwrite
 
@@ -2793,6 +2813,85 @@ def restraint_search_receptor_selection(self, restraint_search_receptor_selectio
                 )
         self._restraint_search_receptor_selection = restraint_search_receptor_selection
 
+    @property
+    def morse_hard_well_depth(self):
+        return self._morse_hard_well_depth
+
+    @morse_hard_well_depth.setter
+    def morse_hard_well_depth(self, morse_hard_well_depth):
+        self._morse_hard_well_depth = self._parse_well_depth(
+            morse_hard_well_depth, "morse_hard_well_depth"
+        )
+
+    @property
+    def morse_soft_well_depth(self):
+        return self._morse_soft_well_depth
+
+    @morse_soft_well_depth.setter
+    def morse_soft_well_depth(self, morse_soft_well_depth):
+        self._morse_soft_well_depth = self._parse_well_depth(
+            morse_soft_well_depth, "morse_soft_well_depth"
+        )
+
+    @property
+    def morse_soft_force_constant(self):
+        return self._morse_soft_force_constant
+
+    @morse_soft_force_constant.setter
+    def morse_soft_force_constant(self, morse_soft_force_constant):
+        if not isinstance(morse_soft_force_constant, str):
+            raise TypeError("'morse_soft_force_constant' must be of type 'str'")
+
+        from sire.units import angstrom, kcal_per_mol
+
+        try:
+            k = _sr.u(morse_soft_force_constant)
+        except:
+            raise ValueError(
+                "Unable to parse 'morse_soft_force_constant' as a Sire "
+                f"GeneralUnit: {morse_soft_force_constant}"
+            )
+
+        if not k.has_same_units(kcal_per_mol / (angstrom * angstrom)):
+            raise ValueError("'morse_soft_force_constant' units are invalid.")
+
+        self._morse_soft_force_constant = k
+
+    @staticmethod
+    def _parse_well_depth(value, name):
+        """
+        Internal helper to validate a Morse potential well depth.
+
+        Parameters
+        ----------
+
+        value: str
+            The well depth as a string, e.g. "150 kcal mol-1".
+
+        name: str
+            The name of the option, used in error messages.
+
+        Returns
+        -------
+
+        well_depth: sire.units.GeneralUnit
+            The parsed well depth.
+        """
+        if not isinstance(value, str):
+            raise TypeError(f"'{name}' must be of type 'str'")
+
+        from sire.units import kcal_per_mol
+
+        try:
+            de = _sr.u(value)
+        except:
+            raise ValueError(f"Unable to parse '{name}' as a Sire GeneralUnit: {value}")
+
+        if not de.has_same_units(kcal_per_mol):
+            raise ValueError(f"'{name}' units are invalid.")
+
+        return de
+
     def _reset_logger(self, logger):
         """
         Internal method to reset the logger.

src/somd2/runner/_base.py

@@ -267,6 +267,24 @@ def __init__(self, system, config):
             self._config._extra_args["use_gcmc_lrc"] = True
             self._config._extra_args["num_gcmc_waters"] = self._config.gcmc_num_waters
 
+        # Auto-generate Morse restraints for ring-breaking perturbations with no
+        # user-supplied restraint. This is done before any modification of the
+        # bonded terms, and before the reference system is stored and the restart
+        # checks are performed, since the hard restraint replaces a bond in the
+        # system, which must match the checkpoints.
+        if self._is_ring_break and self._config.restraints is None:
+            try:
+                self._config.restraints = self._generate_morse_restraints()
+            except Exception as e:
+                msg = (
+                    "Unable to generate Morse restraints for ring-breaking "
+                    f"perturbation: {e}. If the Morse potential has already been "
+                    "applied to the input system, then pass the corresponding "
+                    "restraints using the 'restraints' option."
+                )
+                _logger.error(msg)
+                raise RuntimeError(msg)
+
         # We're running in SOMD1 compatibility mode.
         if self._config.somd1_compatibility:
             from .._utils._somd1 import make_compatible
@@ -1002,6 +1020,76 @@ def _is_abfe_bound(self):
             and self._has_water
         )
 
+    @property
+    def _is_ring_break(self):
+        """
+        Whether this is a ring-breaking (or ring-making) simulation, i.e. one
+        using the 'ring_break_morph' lambda schedule, or its reverse.
+        """
+        return self._config._lambda_schedule_name in (
+            "ring_break_morph",
+            "reverse_ring_break_morph",
+        )
+
+    def _generate_morse_restraints(self):
+        """
+        Return the pair of Morse restraints required by the 'ring_break_morph'
+        lambda schedule, or its reverse. Called automatically when running a
+        ring-breaking simulation with no user-supplied restraint.
+
+        The "hard" restraint directly replaces the harmonic bond that is broken
+        (or formed) by the perturbation, inheriting its force constant and
+        equilibrium length. The "soft" restraint acts on the same pair of atoms
+        and holds the broken fragment in place while the hard restraint is
+        switched off.
+
+        Returns
+        -------
+
+        restraints: [sire.mm.MorsePotentialRestraints]
+            The hard and soft Morse restraints, in the order expected by the
+            schedule's 'morse_hard' and 'morse_soft' levers.
+
+        Notes
+        -----
+
+        As a side effect, ``self._system`` is updated with the replacement of
+        the broken bond by the hard Morse potential.
+        """
+        from sire.restraints import morse_potential as _morse_potential
+
+        _logger.info(
+            "No restraints supplied for ring-breaking perturbation. "
+            "Generating default Morse restraints."
+        )
+
+        hard_restraints, self._system = _morse_potential(
+            self._system,
+            de=self._config.morse_hard_well_depth,
+            auto_parametrise=True,
+            direct_morse_replacement=True,
+            name="morse_hard",
+        )
+
+        # Restrain the same pair of atoms as the hard restraint, at the same
+        # equilibrium distance.
+        soft_restraints, _ = _morse_potential(
+            self._system,
+            atoms0=hard_restraints[0].atom0(),
+            atoms1=hard_restraints[0].atom1(),
+            r0=hard_restraints[0].r0(),
+            k=self._config.morse_soft_force_constant,
+            de=self._config.morse_soft_well_depth,
+            auto_parametrise=False,
+            direct_morse_replacement=False,
+            name="morse_soft",
+        )
+
+        _logger.info(f"Hard Morse restraint: {hard_restraints[0]}")
+        _logger.info(f"Soft Morse restraint: {soft_restraints[0]}")
+
+        return [hard_restraints, soft_restraints]
+
     def _generate_boresch_restraint(self, device=None):
         """
         Return a Boresch restraint for the ABFE simulation, either by loading

tests/runner/test_config.py

@@ -83,3 +83,43 @@ def test_logfile_creation():
         assert Path.exists(runner._config.output_directory / runner._config.log_file)
 
         somd2._logger.remove()
+
+
+def test_morse_restraint_options():
+    """Validate that the Morse restraint options are parsed correctly."""
+    import math
+
+    import pytest
+
+    # The defaults are parsed as Sire units.
+    config = Config()
+    assert config.morse_hard_well_depth == sr.u("150 kcal mol-1")
+    assert config.morse_soft_well_depth == sr.u("50 kcal mol-1")
+    assert config.morse_soft_force_constant == sr.u("125 kcal mol-1 A-2")
+
+    # Equivalent units are accepted, and converted.
+    config = Config(morse_hard_well_depth="418.4 kJ mol-1")
+    assert math.isclose(
+        config.morse_hard_well_depth.to(sr.units.kcal_per_mol), 100.0, rel_tol=1e-6
+    )
+
+    # Well depths must be energies.
+    for option in ("morse_hard_well_depth", "morse_soft_well_depth"):
+        with pytest.raises(TypeError):
+            Config(**{option: 150})
+
+        with pytest.raises(ValueError, match="Unable to parse"):
+            Config(**{option: "not a unit"})
+
+        with pytest.raises(ValueError, match="units are invalid"):
+            Config(**{option: "150 kcal mol-1 A-2"})
+
+    # The force constant must be an energy per unit area.
+    with pytest.raises(TypeError):
+        Config(morse_soft_force_constant=125)
+
+    with pytest.raises(ValueError, match="Unable to parse"):
+        Config(morse_soft_force_constant="not a unit")
+
+    with pytest.raises(ValueError, match="units are invalid"):
+        Config(morse_soft_force_constant="125 kcal mol-1")