From e116413c1f2074315accf58e420278311ab5083f Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Mon, 2 Feb 2026 13:33:24 +0000 Subject: [PATCH 1/2] Update Morse restraint functionality to support automatic parametrisation when a bond is being alchemically created in addition to annihiliation [ci skip] --- doc/source/tutorial/part06/03_restraints.rst | 4 +-- src/sire/restraints/_restraints.py | 34 +++++++++++++++---- tests/conftest.py | 4 +++ .../test_morse_potential_restraints.py | 18 ++++++++-- 4 files changed, 49 insertions(+), 11 deletions(-) diff --git a/doc/source/tutorial/part06/03_restraints.rst b/doc/source/tutorial/part06/03_restraints.rst index 7a43feeea..da56cd4e9 100644 --- a/doc/source/tutorial/part06/03_restraints.rst +++ b/doc/source/tutorial/part06/03_restraints.rst @@ -360,14 +360,14 @@ Morse Potential Restraints --------------------------- The :func:`sire.restraints.morse_potential` function is used to create Morse potential restraints, -which can be used to carry harmonic bond annihilations alchemical relative binding free energy calculations. +which can be used to carry harmonic bond annihilations or creations in alchemical relative binding free energy calculations. To create a Morse potential restraint, you need to specify the two atoms to be restrained. Like the distance restraints, the atoms can be specified using a search string, passing lists of atom indexes, or molecule views holding the atoms. You have to specify the bond force constants, equilibrium bond distance value and the dissociation energy for the restraints. If not supplied, automatic parametrisation feature can be used, which will detect the bond being alchemically -annihilated and set the parameters accordingly (dissociation energy value still needs to be provided). For example, +annihilated or created and set the parameters accordingly (dissociation energy value still needs to be provided). For example, >>> mols = sr.load_test_files("cyclopentane_cyclohexane.bss") >>> morse_restraints = sr.restraints.morse_potential( diff --git a/src/sire/restraints/_restraints.py b/src/sire/restraints/_restraints.py index cc443facc..c6dab643d 100644 --- a/src/sire/restraints/_restraints.py +++ b/src/sire/restraints/_restraints.py @@ -790,8 +790,8 @@ def morse_potential( If True, will attempt to automatically parametrise the Morse potential from a perturbation that annihilates a bond. This requires that 'mols' contains exactly one molecule that is perturbable, and that this - molecule contains exactly one bond that is annihilated at lambda=1. - The atoms involved in the annihilated bond will be used as 'atoms0' + molecule contains exactly one bond that is annihilated or created. + The atoms involved in this bond will be used as 'atoms0' and 'atoms1', the equilibrium distance r0 will be set to the original bond length, and the force constant k will be set to the force constant of the bond in the unperturbed state. Note that 'de' must still be provided. @@ -838,24 +838,44 @@ def morse_potential( if len(ref_mol) != 1: raise ValueError( - "We need exactly one molecule that is perturbable to automatically " + "Exactly one perturbable molecule is required to automatically " "set up the Morse potential restraints" ) perturbable_mol = ref_mol[0] pert = perturbable_mol.perturbation(map=map) pert_omm = pert.to_openmm() changed_bonds = pert_omm.changed_bonds(to_pandas=False) + changed_bonds_df = pert_omm.changed_bonds(to_pandas=True) - # Attempt to find the bond that is annihilated at lambda=1 + # Check that exactly one bond is being created or annihilated, i.e that k0 column + # or k1 column has exactly one zero value + n_bonds_created = (changed_bonds_df["k0"] == 0).sum() + n_bonds_annihilated = (changed_bonds_df["k1"] == 0).sum() + if n_bonds_created + n_bonds_annihilated != 1: + raise ValueError( + "Exactly one bond must be created or annihilated to automatically " + "set up the Morse potential restraints" + ) + + # Attempt to find this bond now for bond in changed_bonds: bond_name, length0, length1, k0, k1 = bond - if k1 == 0: + if k1 == 0 or k0 == 0: - atom0_idx = [bond_name.atom0().index().value()][0] - atom1_idx = [bond_name.atom1().index().value()][0] + # If the bond is being created (k0 == 0), then we should + # use the parameters from the final state (length1, k1). + # If the bond is being annihilated (k1 == 0), then we don't + # need to do anything as length0 and k0 are already selected. + if k0 == 0: + length0 = length1 + k0 = k1 + length0 = u(f"{length0} nm") + atom0_idx = [bond_name.atom0().index().value()][0] + atom1_idx = [bond_name.atom1().index().value()][0] + # Divide k0 by 2 to convert from force constant to sire half # force constant k if k is None: diff --git a/tests/conftest.py b/tests/conftest.py index 8f6d4fad8..c734adff9 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -155,6 +155,10 @@ def cyclopentane_cyclohexane(): return sr.load_test_files("cyclopentane_cyclohexane.bss") +@pytest.fixture(scope="session") +def propane_cyclopropane(): + return sr.load_test_files("propane_cyclopropane.bss") + @pytest.fixture(scope="session") def pentane_cyclopentane(): return sr.load_test_files("pentane_cyclopentane.bss") diff --git a/tests/restraints/test_morse_potential_restraints.py b/tests/restraints/test_morse_potential_restraints.py index 4addbcfd1..2eb2fcd29 100644 --- a/tests/restraints/test_morse_potential_restraints.py +++ b/tests/restraints/test_morse_potential_restraints.py @@ -21,8 +21,9 @@ def test_morse_potential_restraints_setup(cyclopentane_cyclohexane): assert restraints[0].de().value() == 50.0 -def test_morse_potential_restraint_auto_param(cyclopentane_cyclohexane): - """Tests that morse_potential restraints can be set up correctly with automatic parametrisation.""" +def test_morse_potential_restraint_annihiliation_auto_param(cyclopentane_cyclohexane): + """Tests that morse_potential restraints can be set up correctly with automatic parametrisation + when a bond is to be annihilated.""" mols = cyclopentane_cyclohexane.clone() restraints = sr.restraints.morse_potential( mols, @@ -35,6 +36,19 @@ def test_morse_potential_restraint_auto_param(cyclopentane_cyclohexane): assert restraints[0].k().value() == pytest.approx(228.89, rel=0.1) assert restraints[0].de().value() == 25.0 +def test_morse_potential_restraint_creation_auto_param(propane_cyclopropane): + """Tests that morse_potential restraints can be set up correctly with automatic parametrisation + when a bond is to be created.""" + mols = propane_cyclopropane.clone() + restraints = sr.restraints.morse_potential( + mols, + de="25 kcal mol-1", + auto_parametrise=True, + ) + assert restraints.num_restraints() == 1 + assert restraints[0].atom0() == 0 + assert restraints[0].atom1() == 2 + def test_morse_potential_restraint_auto_param_override(cyclopentane_cyclohexane): """Tests that morse_potential restraints can be set up correctly with automatic parametrisation and some parameters can be overwritten.""" From 3ee47c90b87be0ed322a5c57b8efd6a4ab32868e Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Mon, 2 Feb 2026 13:38:23 +0000 Subject: [PATCH 2/2] Update changelog entry --- doc/source/changelog.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index e6f6f24f6..baeb027a0 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -32,6 +32,9 @@ organisation on `GitHub `__. * Fix :meth:`Dynamics.get_rest2_scale()` method. +* Add automatic parametrisation for Morse potential restraints when bonds are + being created in alchemical simulations. + `2025.3.0 `__ - November 2025 ---------------------------------------------------------------------------------------------