From 691b0b70a33bb1237df57a433cc666f0617dd5ad Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 17 Sep 2026 09:12:41 +0100 Subject: [PATCH] Backport fix from PR #17. [ci skip] --- src/ghostly/_ghostly.py | 52 ++++++++++--- tests/test_ghostly.py | 157 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 199 insertions(+), 10 deletions(-) diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index 21231fa..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,27 +1127,45 @@ 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 + + # 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()) diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index ad263cf..334d3be 100644 --- a/tests/test_ghostly.py +++ b/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.