Remove cross-bond dihedrals spanning ring-making/breaking bonds at absent state.

Author: lohedges

src/ghostly/_ghostly.py

@@ -451,12 +451,15 @@ def modify(
             mol, ghosts0, modifications, skip_ghosts=linearised0, is_lambda1=False
         )
 
-        # Remove angles that span ring-making bonds in the state where those
-        # bonds do not yet exist.
+        # Remove angles and dihedrals that span ring-making bonds in the state
+        # where those bonds do not yet exist.
         if _ring_making_bonds:
             mol = _remove_cross_bond_angles(
                 mol, _ring_making_bonds, modifications, is_lambda1=False
             )
+            mol = _remove_cross_bond_dihedrals(
+                mol, _ring_making_bonds, modifications, is_lambda1=False
+            )
 
         # Soften any surviving mixed ghost/physical dihedrals.
         mol = _soften_mixed_dihedrals(
@@ -572,12 +575,15 @@ def modify(
             mol, ghosts1, modifications, skip_ghosts=linearised1, is_lambda1=True
         )
 
-        # Remove angles that span ring-breaking bonds in the state where those
-        # bonds no longer exist.
+        # Remove angles and dihedrals that span ring-breaking bonds in the state
+        # where those bonds no longer exist.
         if _ring_breaking_bonds:
             mol = _remove_cross_bond_angles(
                 mol, _ring_breaking_bonds, modifications, is_lambda1=True
             )
+            mol = _remove_cross_bond_dihedrals(
+                mol, _ring_breaking_bonds, modifications, is_lambda1=True
+            )
 
         # Soften any surviving mixed ghost/physical dihedrals.
         mol = _soften_mixed_dihedrals(
@@ -2508,6 +2514,88 @@ def _remove_cross_bond_angles(mol, changing_bonds, modifications, is_lambda1=Fal
     return mol
 
 
+def _remove_cross_bond_dihedrals(mol, changing_bonds, modifications, is_lambda1=False):
+    r"""
+    Remove dihedral terms whose central bond spans a ring-making or
+    ring-breaking bond in the end state where that bond does not exist.
+    Such dihedrals encode the bonded geometry and couple atoms across the
+    missing bond, causing poor overlap between adjacent lambda windows.
+
+        A - B ~~~ C - D
+              (missing bond)
+
+    If the B-C bond is absent in this end state, the dihedral A-B-C-D
+    is removed.
+
+    Parameters
+    ----------
+
+    mol : sire.mol.Molecule
+        The perturbable molecule.
+
+    changing_bonds : set of (int, int)
+        Pairs of atom index values for bonds that change between end states.
+        Pass ring-making bonds for is_lambda1=False, ring-breaking bonds for
+        is_lambda1=True.
+
+    modifications : dict
+        A dictionary to store details of the modifications made.
+
+    is_lambda1 : bool, optional
+        Whether to modify dihedrals at lambda = 1.
+
+    Returns
+    -------
+
+    mol : sire.mol.Molecule
+        The updated molecule.
+    """
+
+    if not changing_bonds:
+        return mol
+
+    info = mol.info()
+
+    if is_lambda1:
+        mod_key = "lambda_1"
+        suffix = "1"
+    else:
+        mod_key = "lambda_0"
+        suffix = "0"
+
+    dihedrals = mol.property("dihedral" + suffix)
+    new_dihedrals = _SireMM.FourAtomFunctions(mol.info())
+    modified = False
+
+    for p in dihedrals.potentials():
+        idx0 = info.atom_idx(p.atom0())
+        idx1 = info.atom_idx(p.atom1())
+        idx2 = info.atom_idx(p.atom2())
+        idx3 = info.atom_idx(p.atom3())
+
+        i, j, k, l = idx0.value(), idx1.value(), idx2.value(), idx3.value()
+        central = (min(j, k), max(j, k))
+
+        if central in changing_bonds:
+            _logger.debug(
+                f"  Removing cross-bond dihedral: [{i}-{j}-{k}-{l}], {p.function()}"
+            )
+            modifications[mod_key]["removed_dihedrals"].append(f"{i},{j},{k},{l}")
+            modified = True
+        else:
+            new_dihedrals.set(idx0, idx1, idx2, idx3, p.function())
+
+    if modified:
+        mol = (
+            mol.edit()
+            .set_property("dihedral" + suffix, new_dihedrals)
+            .molecule()
+            .commit()
+        )
+
+    return mol
+
+
 def _soften_mixed_dihedrals(
     mol, ghosts, modifications, soften_anchors=1.0, is_lambda1=False
 ):