From 771c8cd918c77ce4b8b05c8c23a82b07b30440b3 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 24 Aug 2026 12:11:11 +0100 Subject: [PATCH 1/2] Add support for auto-generating Morse restraints for ring-breaking. --- CHANGELOG.md | 1 + src/somd2/config/_config.py | 99 +++++++++ src/somd2/runner/_base.py | 88 ++++++++ tests/runner/test_config.py | 40 ++++ tests/runner/test_ring_break_restraints.py | 232 +++++++++++++++++++++ 5 files changed, 460 insertions(+) create mode 100644 tests/runner/test_ring_break_restraints.py diff --git a/CHANGELOG.md b/CHANGELOG.md index d9dbe38..ddfcb67 100644 --- a/CHANGELOG.md +++ b/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 -------------------------------------------------------------------------------------- diff --git a/src/somd2/config/_config.py b/src/somd2/config/_config.py index adaa0e6..28f2470 100644 --- a/src/somd2/config/_config.py +++ b/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. diff --git a/src/somd2/runner/_base.py b/src/somd2/runner/_base.py index 044e3a9..2f4eea8 100644 --- a/src/somd2/runner/_base.py +++ b/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 diff --git a/tests/runner/test_config.py b/tests/runner/test_config.py index 39bd948..27d8e4a 100644 --- a/tests/runner/test_config.py +++ b/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") diff --git a/tests/runner/test_ring_break_restraints.py b/tests/runner/test_ring_break_restraints.py new file mode 100644 index 0000000..255931a --- /dev/null +++ b/tests/runner/test_ring_break_restraints.py @@ -0,0 +1,232 @@ +import tempfile + +import pytest + +from somd2.config import Config +from somd2.runner import Runner + + +def _config(tmpdir, **kwargs): + """Return a minimal ring-breaking config rooted at 'tmpdir'.""" + options = { + "output_directory": tmpdir, + "lambda_schedule": "ring_break_morph", + "num_lambda": 3, + "runtime": "12fs", + "energy_frequency": "4fs", + "frame_frequency": "4fs", + "checkpoint_frequency": "4fs", + "equilibration_time": "0fs", + "minimise": False, + "platform": "CPU", + "max_threads": 1, + } + options.update(kwargs) + return Config(**options) + + +def _restraint(restraints, name): + """Return the single restraint from the set called 'name'.""" + for restraint_set in restraints: + if str(restraint_set.name()) == name: + assert len(restraint_set) == 1 + return restraint_set[0] + raise AssertionError(f"No restraint set named {name!r} in {restraints}") + + +def test_restraints_are_generated(syk_ring_break_mols): + """ + Ensure that a pair of Morse restraints is automatically generated for a + ring-breaking perturbation when no restraint is supplied, and that they act + on the same pair of atoms, at the same equilibrium distance. + """ + with tempfile.TemporaryDirectory() as tmpdir: + runner = Runner(syk_ring_break_mols.clone(), _config(tmpdir)) + + restraints = runner._config.restraints + assert restraints is not None + assert len(restraints) == 2 + + hard = _restraint(restraints, "morse_hard") + soft = _restraint(restraints, "morse_soft") + + # Both restraints act on the bond that is broken. + assert hard.atom0() == soft.atom0() + assert hard.atom1() == soft.atom1() + assert hard.r0() == soft.r0() + + # The restraints are passed through to the dynamics. + assert runner._dynamics_kwargs["restraints"] is restraints + + +def test_restraints_match_config(syk_ring_break_mols): + """ + Ensure that the generated restraints use the well depths and force constant + from the config, and that the hard restraint inherits the force constant of + the bond that it replaces. + """ + import sire as sr + + with tempfile.TemporaryDirectory() as tmpdir: + config = _config( + tmpdir, + morse_hard_well_depth="123 kcal mol-1", + morse_soft_well_depth="45 kcal mol-1", + morse_soft_force_constant="67 kcal mol-1 A-2", + ) + runner = Runner(syk_ring_break_mols.clone(), config) + + hard = _restraint(runner._config.restraints, "morse_hard") + soft = _restraint(runner._config.restraints, "morse_soft") + + assert hard.de() == sr.u("123 kcal mol-1") + assert soft.de() == sr.u("45 kcal mol-1") + assert soft.k() == sr.u("67 kcal mol-1 A-2") + + # The hard restraint is auto-parametrised from the broken bond, so its + # force constant comes from the bond, not the config. + assert hard.k() != soft.k() + assert hard.k().value() > 0 + + +def test_broken_bond_is_replaced(syk_ring_break_mols): + """ + Ensure that the hard restraint replaces the harmonic bond that is broken by + the perturbation, i.e. that the bond is removed from the runner's system. + Leaving both in place would double count the interaction. + """ + with tempfile.TemporaryDirectory() as tmpdir: + mols = syk_ring_break_mols.clone() + runner = Runner(mols, _config(tmpdir)) + + hard = _restraint(runner._config.restraints, "morse_hard") + + def num_bonds(system, idx0, idx1): + """Count the bond potentials between a pair of atom indices.""" + atoms = system.atoms() + atom0 = atoms[idx0] + atom1 = atoms[idx1] + mol = system[atom0.molecule().number()] + info = mol.info() + expected = {atom0.index().value(), atom1.index().value()} + + count = 0 + for bond_prop in ("bond0", "bond1"): + for potential in mol.property(bond_prop).potentials(): + idxs = { + info.atom_idx(potential.atom0()).value(), + info.atom_idx(potential.atom1()).value(), + } + if idxs == expected: + count += 1 + return count + + # The unmodified input still has the bond, in the reference end state + # only, since it is broken by the perturbation. + assert num_bonds(syk_ring_break_mols, hard.atom0(), hard.atom1()) == 1 + + # The runner's system has it removed, replaced by the Morse restraint. + assert num_bonds(runner._system, hard.atom0(), hard.atom1()) == 0 + + +def test_restraints_are_deterministic(syk_ring_break_mols): + """ + Ensure that generating the restraints twice from the same input gives + identical restraints. A restart regenerates them from the input system + rather than reloading them, so they must not drift between runs, otherwise + the accumulated free energy would be invalidated. + """ + with tempfile.TemporaryDirectory() as tmpdir0: + runner0 = Runner(syk_ring_break_mols.clone(), _config(tmpdir0)) + + with tempfile.TemporaryDirectory() as tmpdir1: + runner1 = Runner(syk_ring_break_mols.clone(), _config(tmpdir1)) + + for name in ("morse_hard", "morse_soft"): + assert _restraint(runner0._config.restraints, name) == _restraint( + runner1._config.restraints, name + ) + + +def test_reverse_schedule_generates_restraints(syk_ring_break_mols): + """ + Ensure that restraints are also generated for the ring-making direction, + which uses the reversed schedule. + """ + with tempfile.TemporaryDirectory() as tmpdir: + config = _config(tmpdir, lambda_schedule="reverse_ring_break_morph") + runner = Runner(syk_ring_break_mols.clone(), config) + + assert runner._config.restraints is not None + assert len(runner._config.restraints) == 2 + + +def test_user_restraints_are_not_overridden(syk_ring_break_mols): + """ + Ensure that a user-supplied restraint is left alone, and that the system is + not modified behind their back. + """ + import sire as sr + + mols = syk_ring_break_mols.clone() + + restraints = sr.restraints.distance( + mols, + atoms0=0, + atoms1=1, + k="10 kcal mol-1 A-2", + ) + + with tempfile.TemporaryDirectory() as tmpdir: + runner = Runner(mols, _config(tmpdir, restraints=restraints)) + + assert len(runner._config.restraints) == 1 + assert runner._config.restraints[0] == restraints + + +@pytest.mark.parametrize("schedule", ["standard_morph", "charge_scaled_morph"]) +def test_no_restraints_for_other_schedules(schedule, ethane_methanol): + """ + Ensure that Morse restraints are only generated for ring-breaking + schedules. + """ + with tempfile.TemporaryDirectory() as tmpdir: + config = _config(tmpdir, lambda_schedule=schedule) + runner = Runner(ethane_methanol.clone(), config) + + assert runner._config.restraints is None + + +def test_no_broken_bond_raises(ethane_methanol): + """ + Ensure that a clear error is raised when a ring-breaking schedule is used + for a perturbation that doesn't break (or form) a bond. + """ + with tempfile.TemporaryDirectory() as tmpdir: + with pytest.raises(RuntimeError, match="Unable to generate Morse restraints"): + Runner(ethane_methanol.clone(), _config(tmpdir)) + + +def test_already_applied_raises(syk_ring_break_mols): + """ + Ensure that a helpful error is raised if the Morse potential has already + been applied to the input system, but the corresponding restraints were not + passed via the config. The replacement must not be applied twice. + """ + import sire as sr + + mols = syk_ring_break_mols.clone() + + # Apply the Morse replacement, as a user following the existing workflow + # would, but don't pass the restraints to the config. + _, mols = sr.restraints.morse_potential( + mols, + de="150 kcal mol-1", + auto_parametrise=True, + direct_morse_replacement=True, + name="morse_hard", + ) + + with tempfile.TemporaryDirectory() as tmpdir: + with pytest.raises(RuntimeError, match="already been applied"): + Runner(mols, _config(tmpdir)) From 57f3693ef5e8be61a9c89ccd265f0928dc4b3356 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 24 Aug 2026 12:34:23 +0100 Subject: [PATCH 2/2] Document schedules, ABFE, and ring-breaking. [ci skip] --- README.md | 81 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) diff --git a/README.md b/README.md index dd22fec..26966d0 100644 --- a/README.md +++ b/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