From 746d0f45ae5529e8661d7d49fb0d90e052e15d3c Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Mon, 30 Mar 2026 16:37:13 +0100 Subject: [PATCH 1/3] Implement the DMR approach and make it default for using when using Morse restraints --- src/sire/restraints/_restraints.py | 69 ++++++++++++++++++++++++++---- 1 file changed, 61 insertions(+), 8 deletions(-) diff --git a/src/sire/restraints/_restraints.py b/src/sire/restraints/_restraints.py index c7b23ea83..b8be0dd17 100644 --- a/src/sire/restraints/_restraints.py +++ b/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,7 +996,7 @@ 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() @@ -951,8 +1004,8 @@ def morse_potential( 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): From 94e89802b8061e6d5394542b3f98d8a0ed090e5d Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Tue, 31 Mar 2026 12:11:54 +0100 Subject: [PATCH 2/3] Update morse potential tests to support DMR, update changelog and update tutorial describing Morse potentials [ci skip] --- doc/source/changelog.rst | 3 + doc/source/tutorial/part06/03_restraints.rst | 16 +++-- .../test_morse_potential_restraints.py | 67 +++++++++++++++++-- 3 files changed, 75 insertions(+), 11 deletions(-) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 9257e602a..65156f4d2 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -38,6 +38,9 @@ organisation on `GitHub `__. * 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 `__ - February 2026 --------------------------------------------------------------------------------------------- diff --git a/doc/source/tutorial/part06/03_restraints.rst b/doc/source/tutorial/part06/03_restraints.rst index da56cd4e9..e262a326e 100644 --- a/doc/source/tutorial/part06/03_restraints.rst +++ b/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 --------------------------- diff --git a/tests/restraints/test_morse_potential_restraints.py b/tests/restraints/test_morse_potential_restraints.py index f994614a5..e25dfd075 100644 --- a/tests/restraints/test_morse_potential_restraints.py +++ b/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,7 +25,7 @@ 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, @@ -33,7 +33,7 @@ def test_morse_potential_restraint_annihiliation_auto_param(cyclopentane_cyclohe 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 @@ -41,7 +41,7 @@ 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 From b364f8bfba0bf481e3dcefafdf551f972c0ede51 Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Tue, 31 Mar 2026 14:24:03 +0100 Subject: [PATCH 3/3] Update formatting in the Morse restraints --- src/sire/restraints/_restraints.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/sire/restraints/_restraints.py b/src/sire/restraints/_restraints.py index b8be0dd17..37bc78723 100644 --- a/src/sire/restraints/_restraints.py +++ b/src/sire/restraints/_restraints.py @@ -938,8 +938,8 @@ def morse_potential( 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 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) @@ -996,7 +996,7 @@ def morse_potential( except: raise ValueError(f"Unable to parse 'de' as a Sire GeneralUnit: {de}") - mols_atms = mols.atoms() + atoms = mols.atoms() if name is None: restraints = MorsePotentialRestraints() @@ -1004,8 +1004,8 @@ def morse_potential( restraints = MorsePotentialRestraints(name=name) for i, (atom0, atom1) in enumerate(zip(atoms0, atoms1)): - idxs0 = mols_atms.find(atom0) - idxs1 = mols_atms.find(atom1) + idxs0 = atoms.find(atom0) + idxs1 = atoms.find(atom1) if type(idxs0) is int: idxs0 = [idxs0]