Merge pull request #416 from akalpokas/feature_DMR

Feature DMR

Author: lohedges

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

src/sire/restraints/_restraints.py

@@ -732,7 +732,9 @@ def morse_potential(
     de=None,
     use_pbc=None,
     name=None,
-    auto_parametrise=False,
+    auto_parametrise=True,
+    direct_morse_replacement=True,
+    retain_harmonic_bond=False,
     map=None,
 ):
     """
@@ -795,6 +797,16 @@ def morse_potential(
         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.
+        Default is True.
+
+    direct_morse_replacement : bool, optional
+        If True, and if auto_parametrise is True, then the function will attempt to directly
+        replace an existing bond with a Morse potential. Default is True.
+
+    retain_harmonic_bond : bool, optional
+        If True, and if auto_parametrise is True, then the function will only nullify the force
+        constant of the existing harmonic bond, rather than removing the bond potential entirely.
+        If False, then the existing harmonic bond will be removed entirely.
         Default is False.
 
     Returns
@@ -803,6 +815,9 @@ def morse_potential(
         A container of Morse restraints, where the first restraint is
         the MorsePotentialRestraint created. The Morse restraint created can be
         extracted with MorsePotentialRestraints[0].
+
+    mols : sire.system._system.System
+        The system containing the atoms, which will have been modified if auto_parametrise is True and direct_morse_replacement is True.
     """
 
     from .. import u
@@ -875,10 +890,7 @@ def morse_potential(
                 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:
-                    k0 = k0 / 2.0
                     k0 = u(f"{k0} kJ mol-1 nm-2")
                     k = [k0]
 
@@ -897,6 +909,47 @@ def morse_potential(
                     f"molecule property is_perturbable and atomidx {atom1_idx}"
                 ]
                 break
+        if direct_morse_replacement:
+            from ..legacy import MM as _MM
+            import re as _re
+            from ..legacy.CAS import Symbol as _Symbol
+
+            search_pattern = r"r - (\d+\.\d+)"
+            mol = mol[0]
+            info = mol.info()
+
+            # We need to loop through both bond0 and bond1 properties, as we don't know
+            # which one the bond of interest will be in (bond forming or bond breaking)
+            for bond_prop in ("bond0", "bond1"):
+                bonds = mol.property(bond_prop)
+                new_bonds = _MM.TwoAtomFunctions(info)
+
+                for p in bonds.potentials():
+                    idx0 = info.atom_idx(p.atom0())
+                    idx1 = info.atom_idx(p.atom1())
+
+                    # Attempt to match the bond of interest using previously identified atom indices
+                    if idx0.value() == atom0_idx and idx1.value() == atom1_idx:
+                        bond_potential_string = p.function().to_string()
+                        match = _re.search(search_pattern, bond_potential_string)
+
+                        if not match:
+                            raise ValueError(
+                                f"No match found in the string: {bond_potential_string}"
+                            )
+
+                        # If retaining the harmonic bond, then set the harmonic bond force constant
+                        # to zero. Otherwise, remove the harmonic bond from the force list.
+                        if retain_harmonic_bond:
+                            r = float(match.group(1))
+                            amber_bond = _MM.AmberBond(0, r)
+                            expression = amber_bond.to_expression(_Symbol("r"))
+                            new_bonds.set(idx0, idx1, expression)
+                    else:
+                        new_bonds.set(idx0, idx1, p.function())
+
+                mol = mol.edit().set_property(bond_prop, new_bonds).molecule().commit()
+                mols.update(mol)
 
     try:
         atoms0 = _to_atoms(mols, atoms0)
@@ -943,16 +996,16 @@ def morse_potential(
         except:
             raise ValueError(f"Unable to parse 'de' as a Sire GeneralUnit: {de}")
 
-    mols = mols.atoms()
+    atoms = mols.atoms()
 
     if name is None:
         restraints = MorsePotentialRestraints()
     else:
         restraints = MorsePotentialRestraints(name=name)
 
     for i, (atom0, atom1) in enumerate(zip(atoms0, atoms1)):
-        idxs0 = mols.find(atom0)
-        idxs1 = mols.find(atom1)
+        idxs0 = atoms.find(atom0)
+        idxs1 = atoms.find(atom1)
 
         if type(idxs0) is int:
             idxs0 = [idxs0]
@@ -989,7 +1042,7 @@ def morse_potential(
     # Set the use_pbc flag.
     restraints.set_uses_pbc(use_pbc)
 
-    return restraints
+    return restraints, mols
 
 
 def bond(*args, use_pbc=False, **kwargs):

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