diff --git a/src/somd2/_utils/_somd1.py b/src/somd2/_utils/_somd1.py index eddd872..826c6b2 100644 --- a/src/somd2/_utils/_somd1.py +++ b/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() diff --git a/tests/_utils/test_somd1.py b/tests/_utils/test_somd1.py index 6cd0f75..740cb41 100644 --- a/tests/_utils/test_somd1.py +++ b/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): """ diff --git a/tests/conftest.py b/tests/conftest.py index a4ee48c..93f8e0b 100644 --- a/tests/conftest.py +++ b/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)