Merge pull request #517 from OpenBioSim/fix_ring_break_merge

Author: lohedges

src/BioSimSpace/Align/_merge.py

@@ -1314,6 +1314,97 @@ def merge(
     else:
         edit_mol.set_property("connectivity0", conn0)
         edit_mol.set_property("connectivity1", conn1)
+
+        # Remove bonded terms that span ring-making or ring-breaking bonds in
+        # the end state where those bonds are absent. Such terms constrain atoms
+        # toward a bonded geometry that doesn't exist at that end state, causing
+        # large repulsion and poor overlap at the nonbonded/bonded lambda boundary.
+        _bonds0 = {
+            (
+                min(b.atom0().value(), b.atom1().value()),
+                max(b.atom0().value(), b.atom1().value()),
+            )
+            for b in conn0.get_bonds()
+        }
+        _bonds1 = {
+            (
+                min(b.atom0().value(), b.atom1().value()),
+                max(b.atom0().value(), b.atom1().value()),
+            )
+            for b in conn1.get_bonds()
+        }
+        # Bonds present at λ=1 only: remove spanning terms from the λ=0 properties.
+        _ring_making = _bonds1 - _bonds0
+        # Bonds present at λ=0 only: remove spanning terms from the λ=1 properties.
+        _ring_breaking = _bonds0 - _bonds1
+
+        if _ring_making or _ring_breaking:
+            _mol_info = edit_mol.info()
+
+            for _changing, _suffix in [(_ring_making, "0"), (_ring_breaking, "1")]:
+                if not _changing:
+                    continue
+
+                # Angles: remove if the i-j or j-k pair is a changing bond.
+                if "angle" in shared_props:
+                    _angles = edit_mol.property("angle" + _suffix)
+                    _new_angles = _SireMM.ThreeAtomFunctions(_mol_info)
+                    for _p in _angles.potentials():
+                        _i = _mol_info.atom_idx(_p.atom0()).value()
+                        _j = _mol_info.atom_idx(_p.atom1()).value()
+                        _k = _mol_info.atom_idx(_p.atom2()).value()
+                        if (min(_i, _j), max(_i, _j)) not in _changing and (
+                            min(_j, _k),
+                            max(_j, _k),
+                        ) not in _changing:
+                            _new_angles.set(
+                                _mol_info.atom_idx(_p.atom0()),
+                                _mol_info.atom_idx(_p.atom1()),
+                                _mol_info.atom_idx(_p.atom2()),
+                                _p.function(),
+                            )
+                    edit_mol.set_property("angle" + _suffix, _new_angles)
+
+                # Dihedrals: remove if the central j-k pair is a changing bond.
+                if "dihedral" in shared_props:
+                    _dihedrals = edit_mol.property("dihedral" + _suffix)
+                    _new_dihedrals = _SireMM.FourAtomFunctions(_mol_info)
+                    for _p in _dihedrals.potentials():
+                        _j = _mol_info.atom_idx(_p.atom1()).value()
+                        _k = _mol_info.atom_idx(_p.atom2()).value()
+                        if (min(_j, _k), max(_j, _k)) not in _changing:
+                            _new_dihedrals.set(
+                                _mol_info.atom_idx(_p.atom0()),
+                                _mol_info.atom_idx(_p.atom1()),
+                                _mol_info.atom_idx(_p.atom2()),
+                                _mol_info.atom_idx(_p.atom3()),
+                                _p.function(),
+                            )
+                    edit_mol.set_property("dihedral" + _suffix, _new_dihedrals)
+
+                # Impropers: remove if both atoms of any changing bond appear.
+                if "improper" in shared_props:
+                    _impropers = edit_mol.property("improper" + _suffix)
+                    _new_impropers = _SireMM.FourAtomFunctions(_mol_info)
+                    for _p in _impropers.potentials():
+                        _atoms = {
+                            _mol_info.atom_idx(_p.atom0()).value(),
+                            _mol_info.atom_idx(_p.atom1()).value(),
+                            _mol_info.atom_idx(_p.atom2()).value(),
+                            _mol_info.atom_idx(_p.atom3()).value(),
+                        }
+                        if not any(
+                            _a in _atoms and _b in _atoms for _a, _b in _changing
+                        ):
+                            _new_impropers.set(
+                                _mol_info.atom_idx(_p.atom0()),
+                                _mol_info.atom_idx(_p.atom1()),
+                                _mol_info.atom_idx(_p.atom2()),
+                                _mol_info.atom_idx(_p.atom3()),
+                                _p.function(),
+                            )
+                    edit_mol.set_property("improper" + _suffix, _new_impropers)
+
     # Build the intrascale matrices from the per-state connectivity.
     ff = molecule0.property(ff0)
     sf14 = _SireMM.CLJScaleFactor(

tests/Align/test_align.py

@@ -1068,3 +1068,160 @@ def test_ring_breaking_intrascale_m338():
             assert intra1.get(idx_i, idx_j).lj() == pytest.approx(
                 ref_intra1.get(idx_i, idx_j).lj()
             ), f"intra1 lj mismatch at ({i},{j})"
+
+
+@pytest.mark.skipif(
+    not has_antechamber or not has_tleap,
+    reason="Requires antechamber and tLEaP to be installed.",
+)
+@pytest.mark.skipif(
+    not has_openff,
+    reason="Requires OpenFF to be installed.",
+)
+def test_ring_breaking_cross_bond_cleanup():
+    """
+    Test that bonded terms spanning a ring-breaking bond are removed from the
+    end-state properties where that bond is absent (SYK 5035→5033).
+
+    In this perturbation a ring opens, leaving one bond present at λ=0 but
+    absent at λ=1. Any angle, dihedral or improper whose geometry depends on
+    that bond must be removed from the λ=1 properties; retaining them would
+    constrain atoms toward a bonded geometry that no longer exists and cause
+    large repulsion at the nonbonded/bonded lambda boundary.
+    """
+
+    # MCS mapping: {5033_idx: 5035_idx} — mol0=5033 (ring present), mol1=5035
+    # (ring absent), so the ring bond appears in connectivity0 but not
+    # connectivity1, giving ring_breaking = {(1, 7)} in the merged molecule.
+    mapping = {
+        6: 0,
+        5: 1,
+        4: 2,
+        3: 3,
+        34: 4,
+        8: 5,
+        32: 6,
+        31: 7,
+        11: 8,
+        12: 9,
+        13: 10,
+        14: 11,
+        15: 12,
+        16: 13,
+        17: 14,
+        18: 15,
+        19: 16,
+        20: 17,
+        21: 18,
+        22: 19,
+        23: 20,
+        24: 21,
+        25: 22,
+        26: 23,
+        27: 24,
+        28: 25,
+        29: 26,
+        30: 27,
+        10: 28,
+        9: 29,
+        7: 30,
+        38: 31,
+        37: 32,
+        35: 33,
+        36: 34,
+        1: 35,
+        33: 36,
+        53: 38,
+        52: 39,
+        43: 40,
+        44: 41,
+        45: 42,
+        46: 43,
+        47: 44,
+        48: 45,
+        49: 46,
+        50: 47,
+        51: 48,
+        42: 49,
+        41: 50,
+    }
+
+    mol0 = BSS.Parameters.openff_unconstrained_2_2_1(
+        BSS.IO.readMolecules(f"{url}/5033.sdf")[0]
+    ).getMolecule()
+    mol1 = BSS.Parameters.openff_unconstrained_2_2_1(
+        BSS.IO.readMolecules(f"{url}/5035.sdf")[0]
+    ).getMolecule()
+
+    mol0_aligned = BSS.Align.rmsdAlign(mol0, mol1, mapping)
+    merged = BSS.Align.merge(
+        mol0_aligned,
+        mol1,
+        mapping,
+        allow_ring_breaking=True,
+    )
+
+    sire_mol = merged._sire_object
+    mol_info = sire_mol.info()
+
+    conn0 = sire_mol.property("connectivity0")
+    conn1 = sire_mol.property("connectivity1")
+
+    bonds0 = {
+        (
+            min(b.atom0().value(), b.atom1().value()),
+            max(b.atom0().value(), b.atom1().value()),
+        )
+        for b in conn0.get_bonds()
+    }
+    bonds1 = {
+        (
+            min(b.atom0().value(), b.atom1().value()),
+            max(b.atom0().value(), b.atom1().value()),
+        )
+        for b in conn1.get_bonds()
+    }
+
+    # Bonds present only at λ=1 must not appear in angle0/dihedral0/improper0.
+    ring_making = bonds1 - bonds0
+    # Bonds present only at λ=0 must not appear in angle1/dihedral1/improper1.
+    ring_breaking = bonds0 - bonds1
+
+    # This perturbation must have at least one ring-breaking bond.
+    assert ring_breaking, "Expected ring-breaking bonds in SYK 5035→5033"
+
+    for changing, suffix in [(ring_making, "0"), (ring_breaking, "1")]:
+        if not changing:
+            continue
+
+        for p in sire_mol.property(f"angle{suffix}").potentials():
+            i = mol_info.atom_idx(p.atom0()).value()
+            j = mol_info.atom_idx(p.atom1()).value()
+            k = mol_info.atom_idx(p.atom2()).value()
+            assert (min(i, j), max(i, j)) not in changing, (
+                f"angle{suffix} ({i},{j},{k}) spans absent bond "
+                f"({min(i, j)},{max(i, j)})"
+            )
+            assert (min(j, k), max(j, k)) not in changing, (
+                f"angle{suffix} ({i},{j},{k}) spans absent bond "
+                f"({min(j, k)},{max(j, k)})"
+            )
+
+        for p in sire_mol.property(f"dihedral{suffix}").potentials():
+            j = mol_info.atom_idx(p.atom1()).value()
+            k = mol_info.atom_idx(p.atom2()).value()
+            assert (min(j, k), max(j, k)) not in changing, (
+                f"dihedral{suffix} central bond ({j},{k}) spans absent bond"
+            )
+
+        for p in sire_mol.property(f"improper{suffix}").potentials():
+            atoms = {
+                mol_info.atom_idx(p.atom0()).value(),
+                mol_info.atom_idx(p.atom1()).value(),
+                mol_info.atom_idx(p.atom2()).value(),
+                mol_info.atom_idx(p.atom3()).value(),
+            }
+            for a, b in changing:
+                assert not (a in atoms and b in atoms), (
+                    f"improper{suffix} spans absent bond ({a},{b})"
+                )