Implement the DMR approach and make it default for using when using Morse restraints

Author: akalpokas

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 we retaining the harmonic bond, then set the harmonic bond force constant
+                        # to zero. Otherwise, the harmonic bond entirely removed 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()
+    mols_atms = 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 = mols_atms.find(atom0)
+        idxs1 = mols_atms.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):