Store ring-breaking and ring-making bond pairs as molecule properties.

Author: lohedges

src/BioSimSpace/Align/_merge.py

@@ -1341,6 +1341,24 @@ def merge(
         if _ring_making or _ring_breaking:
             _mol_info = edit_mol.info()
 
+            # Store the changing bond pairs as molecule properties so the
+            # simulation engine can apply a softcore force to the pair in the
+            # end state where the bond is absent, preventing large repulsion at
+            # the bonded/nonbonded lambda boundary. Ring-breaking (absent at
+            # λ=1) and ring-making (absent at λ=0) are stored separately so
+            # the engine can apply different alpha schedules to each.
+            for _pairs, _prop in [
+                (_ring_breaking, "ring_breaking_bonds"),
+                (_ring_making, "ring_making_bonds"),
+            ]:
+                if _pairs:
+                    edit_mol.set_property(
+                        _prop,
+                        _SireBase.IntegerArrayProperty(
+                            [_idx for _pair in sorted(_pairs) for _idx in _pair]
+                        ),
+                    )
+
             for _changing, _suffix in [(_ring_making, "0"), (_ring_breaking, "1")]:
                 if not _changing:
                     continue

tests/Align/test_align.py

@@ -1225,3 +1225,19 @@ def test_ring_breaking_cross_bond_cleanup():
                 assert not (a in atoms and b in atoms), (
                     f"improper{suffix} spans absent bond ({a},{b})"
                 )
+
+    # Check that the ring-breaking and ring-making bond properties are set.
+    def _read_pairs(prop_name):
+        if not sire_mol.has_property(prop_name):
+            return set()
+        flat = list(sire_mol.property(prop_name).to_list())
+        return {(flat[i], flat[i + 1]) for i in range(0, len(flat), 2)}
+
+    stored_breaking = _read_pairs("ring_breaking_bonds")
+    stored_making = _read_pairs("ring_making_bonds")
+    assert stored_breaking == ring_breaking, (
+        f"ring_breaking_bonds property mismatch: {stored_breaking} != {ring_breaking}"
+    )
+    assert stored_making == ring_making, (
+        f"ring_making_bonds property mismatch: {stored_making} != {ring_making}"
+    )