Fix make_compatible dropping unique non-ghost terms for ring-breaking systems

Author: lohedges

src/somd2/_utils/_somd1.py

@@ -264,6 +264,19 @@ def make_compatible(system, fix_perturbable_zero_sigmas=False):
                 new_bonds0.set(idx0, idx1, p0.function())
                 new_bonds1.set(idx0, idx1, p1.function())
 
+        # Pass through unique terms that have no ghost in the state they exist in.
+        for b_idx in bonds0_unique_idx.values():
+            p = bonds0[b_idx]
+            a0, a1 = p.atom0(), p.atom1()
+            if not _has_ghost(mol, [a0, a1]):
+                new_bonds0.set(a0, a1, p.function())
+
+        for b_idx in bonds1_unique_idx.values():
+            p = bonds1[b_idx]
+            a0, a1 = p.atom0(), p.atom1()
+            if not _has_ghost(mol, [a0, a1], True):
+                new_bonds1.set(a0, a1, p.function())
+
         # Set the new bonded terms.
         edit_mol = edit_mol.set_property("bond0", new_bonds0).molecule()
         edit_mol = edit_mol.set_property("bond1", new_bonds1).molecule()
@@ -366,6 +379,19 @@ def make_compatible(system, fix_perturbable_zero_sigmas=False):
                 new_angles0.set(idx0, idx1, idx2, p0.function())
                 new_angles1.set(idx0, idx1, idx2, p1.function())
 
+        # Pass through unique terms that have no ghost in the state they exist in.
+        for a_idx in angles0_unique_idx.values():
+            p = angles0[a_idx]
+            a0, a1, a2 = p.atom0(), p.atom1(), p.atom2()
+            if not _has_ghost(mol, [a0, a1, a2]):
+                new_angles0.set(a0, a1, a2, p.function())
+
+        for a_idx in angles1_unique_idx.values():
+            p = angles1[a_idx]
+            a0, a1, a2 = p.atom0(), p.atom1(), p.atom2()
+            if not _has_ghost(mol, [a0, a1, a2], True):
+                new_angles1.set(a0, a1, a2, p.function())
+
         # Set the new angle terms.
         edit_mol = edit_mol.set_property("angle0", new_angles0).molecule()
         edit_mol = edit_mol.set_property("angle1", new_angles1).molecule()
@@ -479,6 +505,25 @@ def make_compatible(system, fix_perturbable_zero_sigmas=False):
                 new_dihedrals0.set(idx0, idx1, idx2, idx3, p0.function())
                 new_dihedrals1.set(idx0, idx1, idx2, idx3, p1.function())
 
+        # Pass through unique terms that have no ghost in the state they exist in.
+        for d_idx in dihedrals0_unique_idx.values():
+            p = dihedrals0[d_idx]
+            a0 = info.atom_idx(p.atom0())
+            a1 = info.atom_idx(p.atom1())
+            a2 = info.atom_idx(p.atom2())
+            a3 = info.atom_idx(p.atom3())
+            if not _has_ghost(mol, [a0, a1, a2, a3]):
+                new_dihedrals0.set(a0, a1, a2, a3, p.function())
+
+        for d_idx in dihedrals1_unique_idx.values():
+            p = dihedrals1[d_idx]
+            a0 = info.atom_idx(p.atom0())
+            a1 = info.atom_idx(p.atom1())
+            a2 = info.atom_idx(p.atom2())
+            a3 = info.atom_idx(p.atom3())
+            if not _has_ghost(mol, [a0, a1, a2, a3], True):
+                new_dihedrals1.set(a0, a1, a2, a3, p.function())
+
         # Set the new dihedral terms.
         edit_mol = edit_mol.set_property("dihedral0", new_dihedrals0).molecule()
         edit_mol = edit_mol.set_property("dihedral1", new_dihedrals1).molecule()
@@ -605,6 +650,25 @@ def make_compatible(system, fix_perturbable_zero_sigmas=False):
                 new_impropers0.set(idx0, idx1, idx2, idx3, p0.function())
                 new_impropers1.set(idx0, idx1, idx2, idx3, p1.function())
 
+        # Pass through unique terms that have no ghost in the state they exist in.
+        for i_idx in impropers0_unique_idx.values():
+            p = impropers0[i_idx]
+            a0 = info.atom_idx(p.atom0())
+            a1 = info.atom_idx(p.atom1())
+            a2 = info.atom_idx(p.atom2())
+            a3 = info.atom_idx(p.atom3())
+            if not _has_ghost(mol, [a0, a1, a2, a3]):
+                new_impropers0.set(a0, a1, a2, a3, p.function())
+
+        for i_idx in impropers1_unique_idx.values():
+            p = impropers1[i_idx]
+            a0 = info.atom_idx(p.atom0())
+            a1 = info.atom_idx(p.atom1())
+            a2 = info.atom_idx(p.atom2())
+            a3 = info.atom_idx(p.atom3())
+            if not _has_ghost(mol, [a0, a1, a2, a3], True):
+                new_impropers1.set(a0, a1, a2, a3, p.function())
+
         # Set the new improper terms.
         edit_mol = edit_mol.set_property("improper0", new_impropers0).molecule()
         edit_mol = edit_mol.set_property("improper1", new_impropers1).molecule()

tests/_utils/test_somd1.py

@@ -2,11 +2,88 @@
 import sire.legacy.Mol as _SireMol
 
 
+def _unique_nonghost_terms(mol, term_type, n_atoms, final=False):
+    """
+    Return the set of atom-index tuples for terms in `term_type{0 or 1}` that
+    are absent from the other end state and involve no ghost atoms in the state
+    they exist in.
+    """
+    from somd2._utils import _has_ghost
+
+    suffix_own = "1" if final else "0"
+    suffix_other = "0" if final else "1"
+
+    info = mol.info()
+
+    def potentials(suffix):
+        return mol.property(f"{term_type}{suffix}").potentials()
+
+    def key(p):
+        return tuple(
+            info.atom_idx(getattr(p, f"atom{k}")()).value() for k in range(n_atoms)
+        )
+
+    own_keys = {key(p): p for p in potentials(suffix_own)}
+    other_keys = {key(p) for p in potentials(suffix_other)}
+    # also consider reversed keys for symmetric terms
+    other_keys |= {k[::-1] for k in other_keys}
+
+    unique = {}
+    for k, p in own_keys.items():
+        if k not in other_keys:
+            atoms = [info.atom_idx(getattr(p, f"atom{i}")()) for i in range(n_atoms)]
+            if not _has_ghost(mol, atoms, final):
+                unique[k] = p.function()
+    return unique
+
+
 @pytest.fixture
 def mols(request):
     return request.getfixturevalue(request.param)
 
 
+def test_make_compatible_ring_break(ring_break_mols):
+    """
+    Verify that make_compatible preserves non-ghost bonded terms that are
+    unique to one end state, rather than silently dropping them.
+
+    The 6YNGD→intgd perturbation breaks an N-C ring bond. The cross-bond
+    angles, dihedrals, and impropers that span this bond exist only in
+    state0 (the ring is intact there) and involve no ghost atoms, so they
+    must survive make_compatible unchanged.
+    """
+    from somd2._utils._somd1 import make_compatible
+
+    mol_before = ring_break_mols.molecules("property is_perturbable")[0]
+
+    # Collect unique non-ghost terms in state0 before the call.
+    before = {
+        term: _unique_nonghost_terms(mol_before, term, n)
+        for term, n in [("angle", 3), ("dihedral", 4), ("improper", 4)]
+    }
+
+    # Require that there are actually unique non-ghost terms to test against.
+    assert any(before[t] for t in before), (
+        "No unique non-ghost terms found in state0 — test input may be wrong"
+    )
+
+    system_after = make_compatible(ring_break_mols)
+    mol_after = system_after.molecules("property is_perturbable")[0]
+
+    info = mol_after.info()
+
+    for term, n in [("angle", 3), ("dihedral", 4), ("improper", 4)]:
+        after_keys = {
+            tuple(info.atom_idx(getattr(p, f"atom{k}")()).value() for k in range(n))
+            for p in mol_after.property(f"{term}0").potentials()
+        }
+        for atom_key in before[term]:
+            assert atom_key in after_keys or atom_key[::-1] in after_keys, (
+                f"Unique non-ghost {term}0 term {atom_key} was incorrectly "
+                f"removed by make_compatible"
+            )
+
+
 @pytest.mark.parametrize("mols", ["pert_fwd_mols", "pert_rev_mols"], indirect=True)
 def test_reconstruct_intrascale(mols):
     """

tests/conftest.py

@@ -98,3 +98,18 @@ def syk_ring_break_mols():
     """
     mols = sr.load_test_files("syk_5035_5033.s3")
     return sr.morph.link_to_reference(mols)
+
+
+@pytest.fixture(scope="session")
+def ring_break_mols():
+    """
+    Load the 6YNGD→intgd ring-breaking perturbation system.
+
+    Reference state (λ=0): 6YNGD ligand with an intact N-C ring bond.
+    Perturbed state (λ=1): open-chain analogue (intgd) where that bond
+    is absent.  The cross-bond angles, dihedrals, and impropers spanning
+    the breaking bond are non-ghost unique-to-state0 terms and must be
+    preserved by make_compatible.
+    """
+    mols = sr.load_test_files("6yngd_to_intgd.s3")
+    return sr.morph.link_to_reference(mols)