From d8ce416bd9ade637c00cbdebaf4ab830f24f30f6 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 9 Jul 2026 10:41:00 +0100 Subject: [PATCH] Revert auto-zeroing of anchor dihedrals for ring-constrained ghosts. --- CHANGELOG.md | 1 + src/ghostly/_ghostly.py | 153 +++------------------------------------- tests/test_ghostly.py | 41 ++++++----- 3 files changed, 33 insertions(+), 162 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 011861e..19c3201 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. +* Revert auto-zeroing of anchor dihedrals for ring-constrained immediate ghosts: a singly-attached ring still has one unconstrained rotational degree of freedom about the bridge-ghost bond, so removing the anchor reintroduces flapping risk. [2026.1.0](https://github.com/openbiosim/ghostly/compare/2025.2.0...2026.1.0) - Jun 2026 ------------------------------------------------------------------------------------- diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index abd1fa2..ba03534 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -429,20 +429,12 @@ def modify( mol, ghosts0, modifications, skip_ghosts=linearised0, is_lambda1=False ) - # Identify immediate ghost atoms that lie on a ring in the ghost - # subgraph: anchor dihedrals through these are redundant and will be - # zeroed automatically. Appearing ghosts (ghost at λ=0) form their - # ring at λ=1, so connectivity1 is used since that is where the ring - # bonds actually exist. - ring_ghosts0 = _ring_constrained_ghosts(bridges0, ghosts0, connectivity1) - # Soften any surviving mixed ghost/physical dihedrals. mol = _soften_mixed_dihedrals( mol, ghosts0, modifications, soften_anchors=soften_anchors, - ring_ghosts=ring_ghosts0, is_lambda1=False, ) @@ -551,19 +543,12 @@ def modify( mol, ghosts1, modifications, skip_ghosts=linearised1, is_lambda1=True ) - # Identify immediate ghost atoms that lie on a ring in the ghost - # subgraph: anchor dihedrals through these are redundant and will be - # zeroed automatically. Disappearing ghosts (ghost at λ=1) have their - # ring bonds in connectivity0 since that is where they are real. - ring_ghosts1 = _ring_constrained_ghosts(bridges1, ghosts1, connectivity0) - # Soften any surviving mixed ghost/physical dihedrals. mol = _soften_mixed_dihedrals( mol, ghosts1, modifications, soften_anchors=soften_anchors, - ring_ghosts=ring_ghosts1, is_lambda1=True, ) @@ -2444,93 +2429,8 @@ def _remove_ghost_centre_angles( return mol -def _ghost_adjacency(ghosts, connectivity): - """Build an adjacency dict for the ghost-atom subgraph, keyed by integer atom index. - - Using integer indices avoids any AtomIdx Python-wrapper hashing - inconsistencies when looking up atoms returned by connections_to(). - """ - ghost_indices = {g.value() for g in ghosts} - adj = {i: [] for i in ghost_indices} - for g in ghosts: - i = g.value() - for c in connectivity.connections_to(g): - j = c.value() - if j in ghost_indices: - adj[i].append(j) - return adj - - -def _ghost_in_cycle(atom_idx, adj): - """ - Return True if ``atom_idx`` lies on a cycle in the ghost subgraph described - by ``adj`` (an integer-keyed adjacency dict). - - Removes the atom from the graph and does a BFS from each neighbour in turn. - If any neighbour can reach another neighbour of the atom without passing - through the atom itself, they lie on a cycle together. - - We must try each neighbour as a start rather than only ``neighbors[0]`` - because pendant neighbours (bonded only to ``atom_idx`` within the ghost - subgraph) cannot reach the cycle members — e.g. a methyl substituent on a - cyclopropyl carbon would cause a false-negative if picked first. - """ - neighbors = adj.get(atom_idx, []) - if len(neighbors) < 2: - return False - - for start in neighbors: - visited = {start} - queue = [start] - while queue: - node = queue.pop(0) - for n in adj.get(node, []): - if n != atom_idx and n not in visited: - visited.add(n) - queue.append(n) - if any(n in visited for n in neighbors if n != start): - return True - - return False - - -def _ring_constrained_ghosts(bridges, ghosts, connectivity): - """ - Return the set of immediate ghost atoms (directly bonded to a bridge) - that lie on a cycle within the ghost subgraph. - - For these atoms the ring topology already constrains their orientation - relative to the bridge, making anchor dihedrals through them redundant. - The ring bonds prevent free rotation around the bridge-ghost bond just as - well as an anchor dihedral would, so zeroing the anchor introduces no - flapping risk while removing a potential free-energy bias. - - Pass the connectivity of the end state where the ghost atoms are real, - since that is where their ring bonds exist: connectivity1 for appearing - ghosts (ghost at lambda=0), connectivity0 for disappearing ghosts - (ghost at lambda=1). - """ - if not ghosts: - return set() - - adj = _ghost_adjacency(ghosts, connectivity) - ring_ghosts = set() - for ghost_list in bridges.values(): - for g in ghost_list: - if _ghost_in_cycle(g.value(), adj): - ring_ghosts.add(g.value()) - - if ring_ghosts: - _logger.debug( - f"Ring-constrained immediate ghosts (anchors will be zeroed): " - f"[{', '.join(str(i) for i in sorted(ring_ghosts))}]" - ) - - return ring_ghosts - - def _soften_mixed_dihedrals( - mol, ghosts, modifications, soften_anchors=1.0, ring_ghosts=None, is_lambda1=False + mol, ghosts, modifications, soften_anchors=1.0, is_lambda1=False ): r""" Soften surviving mixed ghost/physical dihedral terms by scaling their @@ -2543,15 +2443,9 @@ def _soften_mixed_dihedrals( atoms start gaining softcore nonbonded interactions but are constrained too tightly by bonded terms. - When ``soften_anchors`` is 1.0 (default) and ``ring_ghosts`` is empty, - no modifications are made. When ``soften_anchors`` is 0.0, all mixed - dihedrals are removed. Intermediate values scale the force constants. - - Anchor dihedrals whose bridge-adjacent ghost atom lies on a ring within - the ghost subgraph (supplied via ``ring_ghosts``) are always zeroed, - regardless of ``soften_anchors``: the ring topology already constrains - the ghost orientation, so the anchor is redundant and may introduce a - free-energy bias. + When ``soften_anchors`` is 1.0 (default), no modifications are made. + When 0.0, all mixed dihedrals are removed. Intermediate values scale + the force constants. Parameters ---------- @@ -2568,12 +2462,6 @@ def _soften_mixed_dihedrals( soften_anchors : float, optional Scale factor for mixed dihedral force constants (0.0 to 1.0). - ring_ghosts : set of int, optional - Integer atom indices of ghost atoms that are both directly bonded to a - bridge and part of a ring in the ghost subgraph. Any surviving mixed - dihedral whose bridge-adjacent ghost atom is in this set is removed - regardless of ``soften_anchors``. - is_lambda1 : bool, optional Whether to modify dihedrals at lambda = 1. @@ -2584,12 +2472,8 @@ def _soften_mixed_dihedrals( The updated molecule. """ - if ring_ghosts is None: - ring_ghosts = set() - - # Nothing to do if there are no ghost atoms, no softening is requested, - # and no ring-constrained ghosts require automatic zeroing. - if not ghosts or (soften_anchors >= 1.0 and not ring_ghosts): + # Nothing to do if there are no ghost atoms or no softening is requested. + if not ghosts or soften_anchors >= 1.0: return mol # Store the molecular info. @@ -2625,31 +2509,12 @@ def _soften_mixed_dihedrals( if has_ghost and has_physical: dih_idx_str = ",".join(str(a.value()) for a in atoms) - - # Determine the effective scale for this dihedral. If the - # bridge-adjacent ghost (the ghost atom bonded to a physical atom - # within the four-atom sequence) lies on a ring in the ghost - # subgraph, zero it unconditionally: the ring constrains the ghost - # orientation and the anchor is redundant. - effective_scale = soften_anchors - if ring_ghosts: - for i, a in enumerate(atoms): - if a.value() in ring_ghosts: - neighbors_in_dih = [] - if i > 0: - neighbors_in_dih.append(atoms[i - 1]) - if i < 3: - neighbors_in_dih.append(atoms[i + 1]) - if any(n not in ghosts for n in neighbors_in_dih): - effective_scale = 0.0 - break - - if effective_scale > 0.0: - scaled = p.function() * effective_scale + if soften_anchors > 0.0: + scaled = p.function() * soften_anchors new_dihedrals.set(idx0, idx1, idx2, idx3, scaled) _logger.debug( f" Softening mixed dihedral: [{dih_idx_str}], " - f"scale={effective_scale}" + f"scale={soften_anchors}" ) modifications[mod_key]["softened_dihedrals"].append(dih_idx_str) else: diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index 30b42b0..ad263cf 100644 --- a/tests/test_ghostly.py +++ b/tests/test_ghostly.py @@ -323,12 +323,10 @@ def test_ejm49_to_ejm31(): # The number of angles should remain the same at lambda = 1. assert angles1.num_functions() == new_angles1.num_functions() - # The number of dihedrals should be ten fewer at lambda = 1: four caught - # by the triple junction handler, four bridge-extension terms - # (17-20-{21,25}-{22,24,34,38}), plus two anchor dihedrals - # (16-17-20-{21,25}) that are auto-zeroed because atom 20 (the immediate - # ghost) lies on a ring within the ghost subgraph. - assert dihedrals1.num_functions() - 10 == new_dihedrals1.num_functions() + # The number of dihedrals should be eight fewer at lambda = 1: four caught + # by the triple junction handler, plus four bridge-extension terms + # (17-20-{21,25}-{22,24,34,38}). + assert dihedrals1.num_functions() - 8 == new_dihedrals1.num_functions() # The number of impropers should be six fewer at lambda = 1. assert improper1.num_functions() - 6 == new_improper1.num_functions() @@ -371,9 +369,6 @@ def test_ejm49_to_ejm31(): (AtomIdx(17), AtomIdx(20), AtomIdx(21), AtomIdx(34)), (AtomIdx(17), AtomIdx(20), AtomIdx(25), AtomIdx(24)), (AtomIdx(17), AtomIdx(20), AtomIdx(25), AtomIdx(38)), - # Anchor dihedrals auto-zeroed: atom 20 is ring-constrained. - (AtomIdx(16), AtomIdx(17), AtomIdx(20), AtomIdx(21)), - (AtomIdx(16), AtomIdx(17), AtomIdx(20), AtomIdx(25)), ] # Check that the missing dihedrals are in the original dihedrals at lambda = 1. @@ -388,6 +383,16 @@ def test_ejm49_to_ejm31(): for dihedral in missing_dihedrals1 ) + # The anchor dihedrals (16-17-20-{21,25}) must survive: they are + # real-real-ghost-ghost and are kept to prevent flapping, even though + # atom 20 lies on a ring within the ghost subgraph. + anchor1 = [ + (AtomIdx(16), AtomIdx(17), AtomIdx(20), AtomIdx(21)), + (AtomIdx(16), AtomIdx(17), AtomIdx(20), AtomIdx(25)), + ] + + assert all(check_dihedral(info, new_dihedrals1.potentials(), *d) for d in anchor1) + # Bridge atom 17 is sp2 and not in a ring. With stiffen_sp2_bridges=False # (the default), stiffening is skipped to avoid ~30° strain. Check that # these angles are unchanged at both end states. @@ -512,10 +517,8 @@ def test_ejm31_to_jmc28(): # At lambda = 0, the cyclopropyl group (atoms 32-42) is appearing (ghost). # The per-bridge handlers remove five dihedrals; the bridge-extension pass - # removes six more (17-32-33-{34,35,38} and 17-32-34-{33,36,37}); and the - # three anchor dihedrals (16-17-32-{33,34,42}) are auto-zeroed because - # atom 32 lies on a ring in the ghost subgraph. - assert dihedrals0.num_functions() - 14 == new_dihedrals0.num_functions() + # removes six more (17-32-33-{34,35,38} and 17-32-34-{33,36,37}). + assert dihedrals0.num_functions() - 11 == new_dihedrals0.num_functions() # These six bridge-extension dihedrals should be absent after modification. bridge_extension0 = [ @@ -537,17 +540,19 @@ def test_ejm31_to_jmc28(): check_dihedral(info, new_dihedrals0.potentials(), *d) for d in bridge_extension0 ) - # The anchor dihedrals (16-17-32-{33,34,42}) should also be absent: atom 32 - # is ring-constrained so the ring topology already prevents flapping. + # The anchor dihedrals (16-17-32-{33,34,42}) must survive: they are + # real-real-ghost-ghost and are intentionally kept to prevent flapping, + # even though atom 32 lies on a ring (the cyclopropyl) in the ghost + # subgraph. A single-bond attachment to the bridge leaves one rotational + # degree of freedom (rigid rotation of the ring about the bridge-ghost + # bond) that the ring's internal bonds do not constrain. anchor0 = [ (AtomIdx(16), AtomIdx(17), AtomIdx(32), AtomIdx(33)), (AtomIdx(16), AtomIdx(17), AtomIdx(32), AtomIdx(34)), (AtomIdx(16), AtomIdx(17), AtomIdx(32), AtomIdx(42)), ] - assert not any( - check_dihedral(info, new_dihedrals0.potentials(), *d) for d in anchor0 - ) + assert all(check_dihedral(info, new_dihedrals0.potentials(), *d) for d in anchor0) # At lambda = 1, the single-carbon group (atom 19) is disappearing (ghost). # The per-bridge handlers remove five dihedrals; no bridge-extension terms