Merge pull request #17 from OpenBioSim/fix_16

Author: lohedges

CHANGELOG.md

@@ -5,6 +5,7 @@ Changelog
 ------------------------------------------------------------------------------------
 
 * Please add an item to this CHANGELOG for any new features or bug fixes when creating a PR.
+* Stiffen the intraghost angle to 180 degrees for dual junctions with two ghost branches instead of removing it. With both ghosts confined to the normal of the physical plane, removal left the same-side arrangement degenerate at the end state, giving large free energy errors for e.g. methoxy and hydroxyl perturbations when the input geometry had both ghosts on the same side [#17](https://github.com/OpenBioSim/ghostly/pull/17).
 
 [2026.2.0](https://github.com/openbiosim/ghostly/compare/2026.1.0...2026.2.0) - Sep 2026
 ----------------------------------------------------------------------------------------

src/ghostly/_ghostly.py

@@ -1106,34 +1106,66 @@ def _is_heavy(atom_idx):
         new_angles = _SireMM.ThreeAtomFunctions(mol.info())
         new_dihedrals = _SireMM.FourAtomFunctions(mol.info())
 
+        # The ghosts only sit on the normal of the physical plane if every
+        # physical-bridge-ghost angle is stiffened. The single-branch handler
+        # skips poorly-scoring physical atoms, so mirror that check here.
+        if bridge_in_ring and not stiffen_ring_bridges:
+            preserve_reason = "is in a ring"
+        elif bridge_is_sp2 and not stiffen_sp2_bridges:
+            preserve_reason = "is sp2"
+        elif phys_scores is not None and any(
+            phys_scores[p] > best_phys_score and (heavy_phys - {p}) for p in physical
+        ):
+            preserve_reason = "has a poorly-scoring physical neighbour"
+        else:
+            preserve_reason = None
+
         # Angles.
         for p in angles.potentials():
             idx0 = info.atom_idx(p.atom0())
             idx1 = info.atom_idx(p.atom1())
             idx2 = info.atom_idx(p.atom2())
 
             if idx0 in ghosts and idx2 in ghosts:
-                # When stiffening is skipped for ring/sp2 bridges, the individual
-                # ghost angles are left at their original values, so the intraghost
-                # angle (e.g. H-bridge-H) should also be preserved.
-                if (bridge_in_ring and not stiffen_ring_bridges) or (
-                    bridge_is_sp2 and not stiffen_sp2_bridges
-                ):
+                # When stiffening is skipped, the individual ghost angles are
+                # left at their original values, so the intraghost angle
+                # (e.g. H-bridge-H) should also be preserved.
+                if preserve_reason is not None:
                     new_angles.set(idx0, idx1, idx2, p.function())
                     _logger.debug(
                         f"  Preserving intraghost angle "
                         f"[{idx0.value()}-{idx1.value()}-{idx2.value()}]: "
-                        f"bridge atom {bridge.value()} "
-                        f"{'is in a ring' if bridge_in_ring else 'is sp2'}, "
+                        f"bridge atom {bridge.value()} {preserve_reason}, "
                         f"not stiffening."
                     )
                     continue
+
+                # Stiffening confines each ghost to the normal of the physical
+                # plane, so set the intraghost angle to 180 degrees to keep the
+                # branches on opposite sides. Removing it would leave the
+                # same-side arrangement degenerate at this end state, with a
+                # large barrier to escape. A quarter of k_hard is enough to
+                # remove that barrier without stiffening the bend mode of the
+                # light ghost atoms to the point of timestep instability.
+                from math import pi
+                from sire.legacy.CAS import Symbol
+
+                k_intraghost = 0.25 * k_hard
+
+                expression = _SireMM.AmberAngle(k_intraghost, pi).to_expression(
+                    Symbol("theta")
+                )
+                new_angles.set(idx0, idx1, idx2, expression)
                 _logger.debug(
-                    f"  Removing angle: [{idx0.value()}-{idx1.value()}-{idx2.value()}], {p.function()}"
+                    f"  Stiffening intraghost angle: [{idx0.value()}-{idx1.value()}-{idx2.value()}], "
+                    f"{p.function()} --> {expression}"
                 )
                 ang_idx = (idx0.value(), idx1.value(), idx2.value())
                 ang_idx = ",".join([str(i) for i in ang_idx])
-                modifications[mod_key]["removed_angles"].append(ang_idx)
+                modifications[mod_key]["stiffened_angles"][ang_idx] = {
+                    "k": k_intraghost,
+                    "theta0": 180.0,
+                }
             else:
                 new_angles.set(idx0, idx1, idx2, p.function())
 

tests/test_ghostly.py

@@ -560,6 +560,163 @@ def test_ejm31_to_jmc28():
     assert dihedrals1.num_functions() - 5 == new_dihedrals1.num_functions()
 
 
+def test_pfkfb3_48_to_47():
+    """
+    Test ghost atom modifications for the PFKFB3 ligands 48 to 47 (methoxy to
+    ethyl). This has a dual junction with two ghost branches at lambda = 0:
+    ghost hydrogens 47 and 48 on the ether oxygen (atom 27), whose physical
+    neighbours are atoms 13 and 28.
+
+    Stiffening the physical-bridge-ghost angles to 90 degrees confines each
+    ghost to the normal of the physical plane. The intraghost angle must be
+    stiffened to 180 degrees so that the two ghosts sit on opposite sides of
+    that plane, otherwise the same-side arrangement is degenerate at lambda = 0
+    and cannot escape once the intraghost angle grows in. The input geometry
+    has both ghosts on the same side of the plane.
+    """
+
+    import numpy as np
+    from sire.legacy.Mol import AtomIdx
+
+    mols = sr.load_test_files("pfkfb3_48_47.s3")
+
+    angles0 = mols[0].property("angle0")
+
+    k_hard = 100
+
+    new_mols, modifications = modify(mols, k_hard=k_hard)
+
+    new_angles0 = new_mols[0].property("angle0")
+
+    info = mols[0].info()
+
+    # No angles should be removed.
+    assert angles0.num_functions() == new_angles0.num_functions()
+    assert modifications["lambda_0"]["removed_angles"] == []
+
+    # The four physical-bridge-ghost angles are stiffened to 90 degrees and the
+    # intraghost angle to 180 degrees with a quarter of the force constant.
+    k_intraghost = 0.25 * k_hard
+
+    expected = {
+        (AtomIdx(13), AtomIdx(27), AtomIdx(47)): f"{k_hard} [theta - 1.5708]^2",
+        (AtomIdx(13), AtomIdx(27), AtomIdx(48)): f"{k_hard} [theta - 1.5708]^2",
+        (AtomIdx(28), AtomIdx(27), AtomIdx(47)): f"{k_hard} [theta - 1.5708]^2",
+        (AtomIdx(28), AtomIdx(27), AtomIdx(48)): f"{k_hard} [theta - 1.5708]^2",
+        (
+            AtomIdx(47),
+            AtomIdx(27),
+            AtomIdx(48),
+        ): f"{k_intraghost:g} [theta - 3.14159]^2",
+    }
+
+    found = {}
+    for p in new_angles0.potentials():
+        idx0 = info.atom_idx(p.atom0())
+        idx1 = info.atom_idx(p.atom1())
+        idx2 = info.atom_idx(p.atom2())
+        if (idx0, idx1, idx2) in expected:
+            found[(idx0, idx1, idx2)] = str(p.function())
+        elif (idx2, idx1, idx0) in expected:
+            found[(idx2, idx1, idx0)] = str(p.function())
+
+    assert found == expected
+
+    assert modifications["lambda_0"]["stiffened_angles"]["47,27,48"] == {
+        "k": k_intraghost,
+        "theta0": 180.0,
+    }
+
+    # Minimising at lambda = 0 must move the ghosts to opposite sides of the
+    # physical plane.
+    def ghost_sides(mol):
+        coords = mol.property("coordinates")
+
+        def pos(i):
+            c = coords[i]
+            return np.array([c.x().value(), c.y().value(), c.z().value()])
+
+        bridge = pos(27)
+        normal = np.cross(pos(13) - bridge, pos(28) - bridge)
+        return (
+            np.sign((pos(47) - bridge) @ normal),
+            np.sign((pos(48) - bridge) @ normal),
+        )
+
+    new_mols = sr.morph.link_to_reference(new_mols)
+
+    s0, s1 = ghost_sides(new_mols[0])
+    assert s0 == s1
+
+    minimised = (
+        new_mols.minimisation(lambda_value=0.0, platform="CPU", cutoff_type="rf")
+        .run()
+        .commit()
+    )
+
+    s0, s1 = ghost_sides(minimised[0])
+    assert s0 == -s1
+
+
+def test_pfkfb3_48_to_47_skipped_neighbour():
+    """
+    Test that the intraghost angle is preserved, not stiffened to 180
+    degrees, when stiffening is skipped through a poorly-scoring physical
+    neighbour of the bridge. Atom 28 is made transmuting so that it scores
+    worse than atom 13. The ghosts then no longer sit on the normal of the
+    physical plane, so a 180 degree intraghost angle would be frustrated.
+    """
+
+    from sire.legacy.Mol import AtomIdx, Element
+
+    mols = sr.load_test_files("pfkfb3_48_47.s3")
+
+    mol = mols[0].edit().atom(28).set_property("element1", Element("N")).molecule()
+    mols.update(mol.commit())
+
+    angles0 = mols[0].property("angle0")
+
+    k_hard = 100
+
+    new_mols, modifications = modify(mols, k_hard=k_hard)
+
+    new_angles0 = new_mols[0].property("angle0")
+
+    info = mols[0].info()
+
+    assert angles0.num_functions() == new_angles0.num_functions()
+    assert modifications["lambda_0"]["removed_angles"] == []
+    assert set(modifications["lambda_0"]["stiffened_angles"]) == {
+        "13,27,47",
+        "13,27,48",
+    }
+
+    # Angles through atom 28 and the intraghost angle keep their original form.
+    original = {}
+    for p in angles0.potentials():
+        idx = (
+            info.atom_idx(p.atom0()),
+            info.atom_idx(p.atom1()),
+            info.atom_idx(p.atom2()),
+        )
+        original[idx] = str(p.function())
+
+    preserved = [
+        (AtomIdx(28), AtomIdx(27), AtomIdx(47)),
+        (AtomIdx(28), AtomIdx(27), AtomIdx(48)),
+        (AtomIdx(47), AtomIdx(27), AtomIdx(48)),
+    ]
+
+    for p in new_angles0.potentials():
+        idx = (
+            info.atom_idx(p.atom0()),
+            info.atom_idx(p.atom1()),
+            info.atom_idx(p.atom2()),
+        )
+        if idx in preserved or idx[::-1] in preserved:
+            assert str(p.function()) == original[idx]
+
+
 def check_angle(info, potentials, idx0, idx1, idx2):
     """
     Check if an angle potential is in a list of potentials.