Preserve intraghost angle when a physical neighbour is skipped for stiffening.

Author: lohedges

CHANGELOG.md

@@ -5,7 +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 [#XX](https://github.com/OpenBioSim/ghostly/pull/XX).
+* 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,25 +1106,36 @@ 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
@@ -1133,11 +1144,15 @@ def _is_heavy(atom_idx):
                 # 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.
+                # 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
 
-                expression = _SireMM.AmberAngle(k_hard, pi).to_expression(
+                k_intraghost = 0.25 * k_hard
+
+                expression = _SireMM.AmberAngle(k_intraghost, pi).to_expression(
                     Symbol("theta")
                 )
                 new_angles.set(idx0, idx1, idx2, expression)
@@ -1148,7 +1163,7 @@ def _is_heavy(atom_idx):
                 ang_idx = (idx0.value(), idx1.value(), idx2.value())
                 ang_idx = ",".join([str(i) for i in ang_idx])
                 modifications[mod_key]["stiffened_angles"][ang_idx] = {
-                    "k": k_hard,
+                    "k": k_intraghost,
                     "theta0": 180.0,
                 }
             else:

tests/test_ghostly.py

@@ -595,13 +595,19 @@ def test_pfkfb3_48_to_47():
     assert modifications["lambda_0"]["removed_angles"] == []
 
     # The four physical-bridge-ghost angles are stiffened to 90 degrees and the
-    # intraghost angle to 180 degrees.
+    # 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_hard} [theta - 3.14159]^2",
+        (
+            AtomIdx(47),
+            AtomIdx(27),
+            AtomIdx(48),
+        ): f"{k_intraghost:g} [theta - 3.14159]^2",
     }
 
     found = {}
@@ -617,7 +623,7 @@ def test_pfkfb3_48_to_47():
     assert found == expected
 
     assert modifications["lambda_0"]["stiffened_angles"]["47,27,48"] == {
-        "k": k_hard,
+        "k": k_intraghost,
         "theta0": 180.0,
     }
 
@@ -652,6 +658,65 @@ def pos(i):
     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.