Skip to content

Remove cross-bond terms for ring-making/breaking perturbations #517

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
lohedges merged 1 commit into from
May 7, 2026
Merged

Remove cross-bond terms for ring-making/breaking perturbations #517

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
91 changes: 91 additions & 0 deletions src/BioSimSpace/Align/_merge.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down
157 changes: 157 additions & 0 deletions tests/Align/test_align.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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})"
)
Loading