Merge pull request #14 from OpenBioSim/backport_13

Author: lohedges

CHANGELOG.md

@@ -1,15 +1,20 @@
 Changelog
 =========
 
-[2026.1.0](https://github.com/openbiosim/loch/compare/2025.2.0...2026.1.0) - Jun 2026
+[2026.2.0](https://github.com/openbiosim/ghostly/compare/2026.1.0...2026.2.0) - ********
+-------------------------------------------------------------------------------------
+
+* 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
 -------------------------------------------------------------------------------------
 
 * 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.
 
-[2025.2.0](https://github.com/openbiosim/loch/compare/2025.1.0...2025.2.0) - Mar 2026
+[2025.2.0](https://github.com/openbiosim/ghostly/compare/2025.1.0...2025.2.0) - Mar 2026
 -------------------------------------------------------------------------------------
 
 * Ensure that there are independent, per-state, physical neighbours. This removes the spurious cross-state ghost filter that resulted in a downgrading of junction types when multiple junctions were present in the molecule.
@@ -23,7 +28,7 @@ Changelog
 * Fixed angle removal logic for ring and sp2 bridge atoms: when stiffening is skipped (``stiffen_ring_bridges=False`` or ``stiffen_sp2_bridges=False``), all ghost-containing angles are now preserved at their original force field values rather than removing the intraghost angles (dual junctions) or the sacrificial physical-bridge-ghost angle (triple junctions).
 * Changed the default for ``--optimise-angles`` to ``False``. Angle optimisation is conformer-dependent: different input geometries can yield different equilibrium angles for the same perturbation, adding variability to the resulting force field. The original force field theta0 is sufficient given the small ``k_soft`` value. Optimisation can still be enabled explicitly to follow Boresch et al. (JCTC 2021) strictly.
 
-[2025.1.0](https://github.com/OpenBioSim/loch/releases/tag/2025.1.0) - Nov 2025
+[2025.1.0](https://github.com/OpenBioSim/ghostly/releases/tag/2025.1.0) - Nov 2025
 -------------------------------------------------------------------------------
 
 * Initial public release.

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:

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