Update morse potential tests to support DMR, update changelog and update tutorial describing Morse potentials [ci skip]

Author: akalpokas

doc/source/changelog.rst

@@ -38,6 +38,9 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 
 * Add functionality for coupling one lambda lever to another.
 
+* Added support for Direct Morse Replacement (DMR) feature in ``sire.restraints.morse_potential``
+  which is enabled by default.
+
 `2025.4.0 <https://github.com/openbiosim/sire/compare/2025.3.0...2025.4.0>`__ - February 2026
 ---------------------------------------------------------------------------------------------
 

doc/source/tutorial/part06/03_restraints.rst

@@ -360,7 +360,7 @@ 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 or creations in alchemical relative binding free energy calculations.
+which can be used to carry out 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
@@ -370,7 +370,7 @@ If not supplied, automatic parametrisation feature can be used, which will detec
 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(
+>>> morse_restraints, mols = sr.restraints.morse_potential(
 ...     mols,
 ...     atoms0=mols["molecule property is_perturbable and atomidx 0"],
 ...     atoms1=mols["molecule property is_perturbable and atomidx 4"],
@@ -385,14 +385,22 @@ MorsePotentialRestraint( 0 <=> 4, k=100 kcal mol-1 Å-2 : r0=1.5 Å : de=50 kcal
 creates a Morse potential restraint between atoms 0 and 4 using the specified parameters.
 Alternatively, if the molecule contains a bond that is being alchemically annihilated, e.g.
 
->>> morse_restraints = sr.restraints.morse_potential(mols, auto_parametrise=True, de="50 kcal mol-1")
+>>> morse_restraints, mols = sr.restraints.morse_potential(mols, auto_parametrise=True, de="50 kcal mol-1")
 >>> morse_restraint = morse_restraints[0]
 >>> print(morse_restraint)
-MorsePotentialRestraint( 0 <=> 4, k=228.89 kcal mol-1 Å-2 : r0=1.5354 Å : de=50 kcal mol-1 )
+MorsePotentialRestraint( 0 <=> 4, k=457.78 kcal mol-1 Å-2 : r0=1.5354 Å : de=50 kcal mol-1 )
 
 creates a Morse potential restraint between atoms 0 and 4 by attempting to match the
 bond parameters of the bond being alchemically annihilated.
 
+.. note::
+
+   Automatic bond parametrisation is enabled by default and will strip away the bond being
+   alchemically annihilated or created from the system, and use the parameters of that bond
+   to set the parameters of the Morse potential restraint. If you want to keep the original bond
+   in the system, then you can disable automatic parametrisation by setting ``auto_parametrise=False``
+   and explicitly providing the parameters for the Morse potential restraint.
+
 Boresch Restraints
 ---------------------------
 

tests/restraints/test_morse_potential_restraints.py

@@ -5,7 +5,7 @@
 def test_morse_potential_restraints_setup(cyclopentane_cyclohexane):
     """Tests that morse_potential restraints can be set up correctly with custom parameters."""
     mols = cyclopentane_cyclohexane.clone()
-    restraints = sr.restraints.morse_potential(
+    restraints, mols = sr.restraints.morse_potential(
         mols,
         atoms0=mols["molecule property is_perturbable and atomidx 0"],
         atoms1=mols["molecule property is_perturbable and atomidx 4"],
@@ -25,23 +25,23 @@ def test_morse_potential_restraint_annihiliation_auto_param(cyclopentane_cyclohe
     """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(
+    restraints, mols = 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() == 4
-    assert restraints[0].k().value() == pytest.approx(228.89, rel=0.1)
+    assert restraints[0].k().value() == pytest.approx(457.78, 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(
+    restraints, mols = sr.restraints.morse_potential(
         mols,
         de="25 kcal mol-1",
         auto_parametrise=True,
@@ -54,7 +54,7 @@ def test_morse_potential_restraint_creation_auto_param(propane_cyclopropane):
 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."""
     mols = cyclopentane_cyclohexane.clone()
-    restraints = sr.restraints.morse_potential(
+    restraints, mols = sr.restraints.morse_potential(
         mols,
         de="25 kcal mol-1",
         k="100 kcal mol-1 A-2",
@@ -70,18 +70,71 @@ def test_morse_potential_restraint_auto_param_override(cyclopentane_cyclohexane)
 def test_multiple_morse_potential_restraints(cyclopentane_cyclohexane):
     """Tests that multiple morse_potential restraints can be set up correctly."""
     mols = cyclopentane_cyclohexane.clone()
-    restraints = sr.restraints.morse_potential(
+    restraints, mols = sr.restraints.morse_potential(
         mols,
         de="25 kcal mol-1",
         auto_parametrise=True,
+        direct_morse_replacement=False,
     )
-    restraint1 = sr.restraints.morse_potential(
+    restraint1, mols = sr.restraints.morse_potential(
         mols,
         atoms0=mols["molecule property is_perturbable and atomidx 0"],
         atoms1=mols["molecule property is_perturbable and atomidx 1"],
         k="100 kcal mol-1 A-2",
         r0="1.5 A",
         de="50 kcal mol-1",
+        direct_morse_replacement=False,
     )
     restraints.add(restraint1)
     assert restraints.num_restraints() == 2
+
+def test_morse_potential_direct_morse_replacement(cyclopentane_cyclohexane):
+    """Tests that morse_potential restraints by default will remove the annihilated harmonic bond."""
+    mols = cyclopentane_cyclohexane.clone()
+    bonds0_org_mol = mols[0].property("bond0")
+    bonds1_org_mol = mols[0].property("bond1")
+    num_bonds0_org_mol = len(bonds0_org_mol.potentials())
+    num_bonds1_org_mol = len(bonds1_org_mol.potentials())
+
+    restraints, mols = sr.restraints.morse_potential(
+        mols,
+        de="25 kcal mol-1",
+        auto_parametrise=True,
+    )
+   
+    bonds0_mod_mol = mols[0].property("bond0")
+    bonds1_mod_mol = mols[0].property("bond1")
+    num_bonds0_mod_mol = len(bonds0_mod_mol.potentials())
+    num_bonds1_mod_mol = len(bonds1_mod_mol.potentials())
+
+    # New bonds0 should be smaller by 1 than the original bonds0 if the function removed the bond
+    assert num_bonds0_mod_mol == num_bonds0_org_mol - 1
+
+    # New bonds1 should be same as the original bonds1, as the bonds1 will already be smaller
+    # by 1 than the original bonds0 which is introduced during the merge.
+    assert num_bonds1_mod_mol == num_bonds1_org_mol
+    assert num_bonds0_org_mol == num_bonds1_org_mol + 1
+
+def test_morse_potential_direct_morse_replacement_retain_harmonic(cyclopentane_cyclohexane):
+    """Tests that morse_potential restraints can retain the annihilated harmonic bond, if specified."""
+    mols = cyclopentane_cyclohexane.clone()
+    bonds0_org_mol = mols[0].property("bond0")
+    bonds1_org_mol = mols[0].property("bond1")
+    num_bonds0_org_mol = len(bonds0_org_mol.potentials())
+    num_bonds1_org_mol = len(bonds1_org_mol.potentials())
+
+    restraints, mols = sr.restraints.morse_potential(
+        mols,
+        de="25 kcal mol-1",
+        auto_parametrise=True,
+        retain_harmonic_bond=True,
+    )
+    bonds0_mod_mol = mols[0].property("bond0")
+    bonds1_mod_mol = mols[0].property("bond1")
+    num_bonds0_mod_mol = len(bonds0_mod_mol.potentials())
+    num_bonds1_mod_mol = len(bonds1_mod_mol.potentials())
+
+    # If the harmonic bond is retained, the number of bonds in bonds0 and bonds1
+    # should be the same as the original molecule
+    assert num_bonds0_mod_mol == num_bonds0_org_mol
+    assert num_bonds1_mod_mol == num_bonds1_org_mol