Feature automated Morse restraints during bond creation

doc/source/changelog.rst

@@ -32,6 +32,9 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 
 * 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 <https://github.com/openbiosim/sire/compare/2025.2.0...2025.3.0>`__ - November 2025
 ---------------------------------------------------------------------------------------------
 

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(

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:

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")

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."""