From c995e7c7ecb0e36c2bab10633ae9e3ecd72871db Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 25 Jun 2026 13:28:54 +0100 Subject: [PATCH 1/4] Filter bridge-extension dihedrals. --- CHANGELOG.md | 1 + src/ghostly/_ghostly.py | 39 ++++++++++++++++-- tests/test_ghostly.py | 91 +++++++++++++++++++++++++++++++++++++++-- 3 files changed, 124 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 564ce89..887bde0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,7 @@ Changelog * Please add an item to this CHANGELOG for any new features or bug fixes when creating a PR. * Add linear spacer modification for ring-breaking ghost bridges. * Remove cross-bond angles spanning ring-making/breaking bonds in the state where the bond is absent. +* Fixed missing removal of bridge-extension dihedrals (`real–ghost–ghost–ghost`) that arise when a ghost group contains a ring, e.g. cyclopropyl, where the ring topology creates spurious torsional coupling between the real scaffold and the ghost ring interior. [2025.2.0](https://github.com/openbiosim/loch/compare/2025.1.0...2025.2.0) - Mar 2026 ------------------------------------------------------------------------------------- diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index 2566d9f..ba03534 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -1964,7 +1964,7 @@ def _remove_impropers(mol, ghosts, modifications, is_lambda1=False): def _remove_residual_ghost_dihedrals(mol, ghosts, modifications, is_lambda1=False): r""" Remove dihedral terms that couple ghost and physical regions but were - not caught by the per-bridge junction handlers. This covers two cases: + not caught by the per-bridge junction handlers. This covers three cases: 1. Cross-bridge: both terminal atoms are ghost and both middle atoms are physical. This arises when two ghost groups have adjacent bridge @@ -1991,6 +1991,22 @@ def _remove_residual_ghost_dihedrals(mol, ghosts, modifications, is_lambda1=Fals Removed dihedrals: e.g. R1-X1-DR-X2, X1-DR-X2-R2 + 3. Bridge extension: one terminal is a physical bridge atom and the + remaining three atoms are ghost. This arises when the ghost group + contains a ring (e.g. cyclopropyl): the ring creates additional + bond paths so that dihedrals of the form X-DR1-DR2-DR3 reach from + the bridge into the ring interior. The per-bridge handlers only + remove dihedrals where the ghost terminal is the directly-bonded + ghost neighbour of the bridge; deeper ring atoms are missed. + + DR2---DR3 + / \ + X---DR1 DR4 + \ + DR5 + + Removed dihedrals: e.g. X-DR1-DR2-DR3, DR3-DR2-DR1-X + Parameters ---------- @@ -2060,9 +2076,26 @@ def _remove_residual_ghost_dihedrals(mol, ghosts, modifications, is_lambda1=Fals and (idx1 in ghosts or idx2 in ghosts) ) - if cross_bridge or ghost_middle: + # Case 3: One terminal physical (bridge), three atoms ghost + # (bridge-into-ring extension). The per-bridge handlers only match + # the directly-bonded ghost neighbour as the ghost terminal, so + # dihedrals that continue into a ghost ring are missed. + bridge_extension = ( + idx0 not in ghosts and idx1 in ghosts and idx2 in ghosts and idx3 in ghosts + ) or ( + idx0 in ghosts and idx1 in ghosts and idx2 in ghosts and idx3 not in ghosts + ) + + if cross_bridge or ghost_middle or bridge_extension: + case = ( + "cross-bridge" + if cross_bridge + else "ghost-middle" + if ghost_middle + else "bridge-extension" + ) _logger.debug( - f" Removing residual ghost dihedral: " + f" Removing residual ghost dihedral ({case}): " f"[{idx0.value()}-{idx1.value()}-{idx2.value()}-{idx3.value()}], " f"{p.function()}" ) diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index 4f78b0a..a3d0c6c 100644 --- a/tests/test_ghostly.py +++ b/tests/test_ghostly.py @@ -26,8 +26,12 @@ def test_hexane_to_propane(): # No angles should be removed. assert angles.num_functions() == new_angles.num_functions() - # Six dihedrals should be removed. - assert dihedrals.num_functions() - 6 == new_dihedrals.num_functions() + # Nine dihedrals should be removed: six cross-bridge terms caught by the + # terminal junction handler, plus three bridge-extension terms + # (3-4-5-{17,18,19}) where the real bridge atom (3) is at one terminal + # and three ghost atoms form the rest of the dihedral path into the + # ghost chain. + assert dihedrals.num_functions() - 9 == new_dihedrals.num_functions() # Create dihedral IDs for the missing dihedrals. @@ -40,6 +44,10 @@ def test_hexane_to_propane(): (AtomIdx(11), AtomIdx(2), AtomIdx(3), AtomIdx(14)), (AtomIdx(12), AtomIdx(2), AtomIdx(3), AtomIdx(14)), (AtomIdx(12), AtomIdx(2), AtomIdx(3), AtomIdx(13)), + # Bridge-extension dihedrals into the ghost chain. + (AtomIdx(3), AtomIdx(4), AtomIdx(5), AtomIdx(17)), + (AtomIdx(3), AtomIdx(4), AtomIdx(5), AtomIdx(18)), + (AtomIdx(3), AtomIdx(4), AtomIdx(5), AtomIdx(19)), ] # Store the molecular info. @@ -315,8 +323,11 @@ 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 four fewer at lambda = 1. - assert dihedrals1.num_functions() - 4 == 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}) where the real bridge atom (17) is at one + # terminal and three ghost atoms continue into the ghost group. + 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() @@ -354,6 +365,11 @@ def test_ejm49_to_ejm31(): (AtomIdx(18), AtomIdx(17), AtomIdx(20), AtomIdx(25)), (AtomIdx(20), AtomIdx(17), AtomIdx(16), AtomIdx(33)), (AtomIdx(14), AtomIdx(16), AtomIdx(17), AtomIdx(20)), + # Bridge-extension dihedrals into the ghost group. + (AtomIdx(17), AtomIdx(20), AtomIdx(21), AtomIdx(22)), + (AtomIdx(17), AtomIdx(20), AtomIdx(21), AtomIdx(34)), + (AtomIdx(17), AtomIdx(20), AtomIdx(25), AtomIdx(24)), + (AtomIdx(17), AtomIdx(20), AtomIdx(25), AtomIdx(38)), ] # Check that the missing dihedrals are in the original dihedrals at lambda = 1. @@ -464,6 +480,73 @@ def test_ejm49_to_ejm31(): ) +def test_ejm31_to_jmc28(): + """ + Test ghost atom modifications for the TYK2 ligands EJM31 to JMC28. + + This perturbation involves an appearing methylcyclopropyl group (atoms + 32-42) attached at atom 32 to real bridge atom 17. The cyclopropyl ring + creates dihedral paths of the form 17-32-33-* and 17-32-34-* that extend + from the real bridge into the ghost ring interior. These bridge-extension + dihedrals must be removed to avoid spurious torsional coupling between + the ghost ring and the real scaffold at lambda = 0. + """ + + mols = sr.load_test_files("ejm31_jmc28.s3") + + dihedrals0 = mols[0].property("dihedral0") + dihedrals1 = mols[0].property("dihedral1") + + new_mols, _ = modify(mols) + + new_dihedrals0 = new_mols[0].property("dihedral0") + new_dihedrals1 = new_mols[0].property("dihedral1") + + from sire.legacy.Mol import AtomIdx + + info = mols[0].info() + + # 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}). + assert dihedrals0.num_functions() - 11 == new_dihedrals0.num_functions() + + # These six bridge-extension dihedrals should be absent after modification. + bridge_extension0 = [ + (AtomIdx(17), AtomIdx(32), AtomIdx(33), AtomIdx(34)), + (AtomIdx(17), AtomIdx(32), AtomIdx(33), AtomIdx(35)), + (AtomIdx(17), AtomIdx(32), AtomIdx(33), AtomIdx(38)), + (AtomIdx(17), AtomIdx(32), AtomIdx(34), AtomIdx(33)), + (AtomIdx(17), AtomIdx(32), AtomIdx(34), AtomIdx(36)), + (AtomIdx(17), AtomIdx(32), AtomIdx(34), AtomIdx(37)), + ] + + # Check all bridge-extension dihedrals were present in the original. + assert all( + check_dihedral(info, dihedrals0.potentials(), *d) for d in bridge_extension0 + ) + + # Check all bridge-extension dihedrals are absent after modification. + assert not any( + check_dihedral(info, new_dihedrals0.potentials(), *d) for d in bridge_extension0 + ) + + # The anchor dihedrals (16-17-32-{33,34,42}) must survive: they are + # real-real-ghost-ghost and are intentionally kept to prevent flapping. + 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 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 + # arise since atom 19 is terminal (no further ghost neighbours). + assert dihedrals1.num_functions() - 5 == new_dihedrals1.num_functions() + + def check_angle(info, potentials, idx0, idx1, idx2): """ Check if an angle potential is in a list of potentials. From b925aa44e4fd34e51541a48f499cef75435a8299 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 25 Jun 2026 15:05:58 +0100 Subject: [PATCH 2/4] Auto-zero anchor dihedrals for ring-constrained immediate ghosts. --- CHANGELOG.md | 1 + src/ghostly/_ghostly.py | 139 +++++++++++++++++++++++++++++++++++++--- tests/test_ghostly.py | 28 +++++--- 3 files changed, 149 insertions(+), 19 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 887bde0..541c28d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,7 @@ Changelog * Add linear spacer modification for ring-breaking ghost bridges. * Remove cross-bond angles spanning ring-making/breaking bonds in the state where the bond is absent. * Fixed missing removal of bridge-extension dihedrals (`real–ghost–ghost–ghost`) that arise when a ghost group contains a ring, e.g. cyclopropyl, where the ring topology creates spurious torsional coupling between the real scaffold and the ghost ring interior. +* Auto-zero anchor dihedrals when the immediate ghost atom lies on a ring within the ghost subgraph. The ring topology already constrains the ghost orientation relative to the bridge, making the anchor redundant; retaining it can introduce a free-energy bias. [2025.2.0](https://github.com/openbiosim/loch/compare/2025.1.0...2025.2.0) - Mar 2026 ------------------------------------------------------------------------------------- diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index ba03534..91cf613 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -429,12 +429,20 @@ 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, ) @@ -543,12 +551,19 @@ 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, ) @@ -2429,8 +2444,79 @@ def _remove_ghost_centre_angles( return mol +def _ghost_adjacency(ghosts, connectivity): + """Build an adjacency dict for the ghost-atom subgraph.""" + ghost_set = set(ghosts) + adj = {g: [] for g in ghost_set} + for g in ghost_set: + for c in connectivity.connections_to(g): + if c in ghost_set: + adj[g].append(c) + return adj + + +def _ghost_in_cycle(atom, adj): + """ + Return True if ``atom`` lies on a cycle in the ghost subgraph described + by ``adj``. + + Removes ``atom`` from the graph and does a BFS from the first neighbour. + If any other neighbour of ``atom`` is reachable without passing through + ``atom``, the three form a cycle. + """ + neighbors = adj.get(atom, []) + if len(neighbors) < 2: + return False + + visited = {neighbors[0]} + queue = [neighbors[0]] + while queue: + node = queue.pop(0) + for n in adj.get(node, []): + if n != atom and n not in visited: + visited.add(n) + queue.append(n) + + return any(n in visited for n in neighbors[1:]) + + +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, adj): + ring_ghosts.add(g) + + if ring_ghosts: + _logger.debug( + f"Ring-constrained immediate ghosts (anchors will be zeroed): " + f"[{', '.join(str(g.value()) for g in ring_ghosts)}]" + ) + + return ring_ghosts + + def _soften_mixed_dihedrals( - mol, ghosts, modifications, soften_anchors=1.0, is_lambda1=False + mol, ghosts, modifications, soften_anchors=1.0, ring_ghosts=None, is_lambda1=False ): r""" Soften surviving mixed ghost/physical dihedral terms by scaling their @@ -2443,9 +2529,15 @@ 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), no modifications are made. - When 0.0, all mixed dihedrals are removed. Intermediate values scale - the force constants. + 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. Parameters ---------- @@ -2462,6 +2554,12 @@ def _soften_mixed_dihedrals( soften_anchors : float, optional Scale factor for mixed dihedral force constants (0.0 to 1.0). + ring_ghosts : set, optional + 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. @@ -2472,8 +2570,12 @@ def _soften_mixed_dihedrals( The updated molecule. """ - # Nothing to do if there are no ghost atoms or no softening is requested. - if not ghosts or soften_anchors >= 1.0: + 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): return mol # Store the molecular info. @@ -2509,12 +2611,31 @@ def _soften_mixed_dihedrals( if has_ghost and has_physical: dih_idx_str = ",".join(str(a.value()) for a in atoms) - if soften_anchors > 0.0: - scaled = p.function() * soften_anchors + + # 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 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 new_dihedrals.set(idx0, idx1, idx2, idx3, scaled) _logger.debug( f" Softening mixed dihedral: [{dih_idx_str}], " - f"scale={soften_anchors}" + f"scale={effective_scale}" ) modifications[mod_key]["softened_dihedrals"].append(dih_idx_str) else: diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index a3d0c6c..30b42b0 100644 --- a/tests/test_ghostly.py +++ b/tests/test_ghostly.py @@ -323,11 +323,12 @@ 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 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}) where the real bridge atom (17) is at one - # terminal and three ghost atoms continue into the ghost group. - assert dihedrals1.num_functions() - 8 == new_dihedrals1.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 impropers should be six fewer at lambda = 1. assert improper1.num_functions() - 6 == new_improper1.num_functions() @@ -370,6 +371,9 @@ 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. @@ -508,8 +512,10 @@ 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}). - assert dihedrals0.num_functions() - 11 == new_dihedrals0.num_functions() + # 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() # These six bridge-extension dihedrals should be absent after modification. bridge_extension0 = [ @@ -531,15 +537,17 @@ 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}) must survive: they are - # real-real-ghost-ghost and are intentionally kept to prevent flapping. + # 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. 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 all(check_dihedral(info, new_dihedrals0.potentials(), *d) for d in anchor0) + assert not any( + 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 From 153bb5a2efb226aee4aa5cf6535777076f1ab615 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 29 Jun 2026 10:29:55 +0100 Subject: [PATCH 3/4] Fix intermittent ring-constrained ghost detection by using integer atom indices --- src/ghostly/_ghostly.py | 52 +++++++++++++++++++++++------------------ 1 file changed, 29 insertions(+), 23 deletions(-) diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index 91cf613..b1d832e 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -2445,26 +2445,32 @@ def _remove_ghost_centre_angles( def _ghost_adjacency(ghosts, connectivity): - """Build an adjacency dict for the ghost-atom subgraph.""" - ghost_set = set(ghosts) - adj = {g: [] for g in ghost_set} - for g in ghost_set: + """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): - if c in ghost_set: - adj[g].append(c) + j = c.value() + if j in ghost_indices: + adj[i].append(j) return adj -def _ghost_in_cycle(atom, adj): +def _ghost_in_cycle(atom_idx, adj): """ - Return True if ``atom`` lies on a cycle in the ghost subgraph described - by ``adj``. + Return True if ``atom_idx`` lies on a cycle in the ghost subgraph described + by ``adj`` (an integer-keyed adjacency dict). - Removes ``atom`` from the graph and does a BFS from the first neighbour. - If any other neighbour of ``atom`` is reachable without passing through - ``atom``, the three form a cycle. + Removes the atom from the graph and does a BFS from the first neighbour. + If any other neighbour is reachable without passing through the atom, + they form a cycle. """ - neighbors = adj.get(atom, []) + neighbors = adj.get(atom_idx, []) if len(neighbors) < 2: return False @@ -2473,7 +2479,7 @@ def _ghost_in_cycle(atom, adj): while queue: node = queue.pop(0) for n in adj.get(node, []): - if n != atom and n not in visited: + if n != atom_idx and n not in visited: visited.add(n) queue.append(n) @@ -2503,13 +2509,13 @@ def _ring_constrained_ghosts(bridges, ghosts, connectivity): ring_ghosts = set() for ghost_list in bridges.values(): for g in ghost_list: - if _ghost_in_cycle(g, adj): - ring_ghosts.add(g) + 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(g.value()) for g in ring_ghosts)}]" + f"[{', '.join(str(i) for i in sorted(ring_ghosts))}]" ) return ring_ghosts @@ -2554,11 +2560,11 @@ def _soften_mixed_dihedrals( soften_anchors : float, optional Scale factor for mixed dihedral force constants (0.0 to 1.0). - ring_ghosts : set, optional - 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``. + 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. @@ -2620,7 +2626,7 @@ def _soften_mixed_dihedrals( effective_scale = soften_anchors if ring_ghosts: for i, a in enumerate(atoms): - if a in ring_ghosts: + if a.value() in ring_ghosts: neighbors_in_dih = [] if i > 0: neighbors_in_dih.append(atoms[i - 1]) From b8738410a951b9a02cb25d962e8e99d47ee3238c Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 29 Jun 2026 10:58:20 +0100 Subject: [PATCH 4/4] Fix _ghost_in_cycle to handle pendant ghost neighbours on cyclopropyl atoms --- src/ghostly/_ghostly.py | 34 +++++++++++++++++++++------------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index b1d832e..abd1fa2 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -2466,24 +2466,32 @@ 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 the first neighbour. - If any other neighbour is reachable without passing through the atom, - they form a cycle. + 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 - visited = {neighbors[0]} - queue = [neighbors[0]] - 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) - - return any(n in visited for n in neighbors[1:]) + 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):