From aab4f2aa9620b64b6db30201baef0275bdd6f7ce Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 16 Sep 2026 20:31:43 +0100 Subject: [PATCH 1/2] Stiffen intraghost angle to 180 degrees for two-branch dual junctions. --- CHANGELOG.md | 1 + src/ghostly/_ghostly.py | 21 +++++++++- tests/test_ghostly.py | 92 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 112 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1bee816..143aaeb 100644 --- a/CHANGELOG.md +++ b/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 [#XX](https://github.com/OpenBioSim/ghostly/pull/XX). [2026.2.0](https://github.com/openbiosim/ghostly/compare/2026.1.0...2026.2.0) - Sep 2026 ---------------------------------------------------------------------------------------- diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index 21231fa..72ae3b9 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -1128,12 +1128,29 @@ def _is_heavy(atom_idx): 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. + from math import pi + from sire.legacy.CAS import Symbol + + expression = _SireMM.AmberAngle(k_hard, 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_hard, + "theta0": 180.0, + } else: new_angles.set(idx0, idx1, idx2, p.function()) diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index ad263cf..73ebdd9 100644 --- a/tests/test_ghostly.py +++ b/tests/test_ghostly.py @@ -560,6 +560,98 @@ 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. + 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", + } + + 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_hard, + "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 check_angle(info, potentials, idx0, idx1, idx2): """ Check if an angle potential is in a list of potentials. From a96aa4d57cc0ca6329d98b90ec3eb05884389940 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 16 Sep 2026 20:38:20 +0100 Subject: [PATCH 2/2] Preserve intraghost angle when a physical neighbour is skipped for stiffening. --- CHANGELOG.md | 2 +- src/ghostly/_ghostly.py | 37 ++++++++++++++------- tests/test_ghostly.py | 71 +++++++++++++++++++++++++++++++++++++++-- 3 files changed, 95 insertions(+), 15 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 143aaeb..54639e5 100644 --- a/CHANGELOG.md +++ b/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 ---------------------------------------------------------------------------------------- diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index 72ae3b9..f7a0e6f 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -1106,6 +1106,20 @@ 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()) @@ -1113,18 +1127,15 @@ def _is_heavy(atom_idx): 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: diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index 73ebdd9..334d3be 100644 --- a/tests/test_ghostly.py +++ b/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.