Skip to content

Backport fix from #13 #14

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
lohedges merged 1 commit into from
Jul 9, 2026
Merged

Backport fix from #13 #14

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 9 additions & 4 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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.
Expand All @@ -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.
153 changes: 9 additions & 144 deletions src/ghostly/_ghostly.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
)

Expand Down Expand Up @@ -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,
)

Expand Down Expand Up @@ -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
Expand All @@ -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
----------
Expand All @@ -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.

Expand All @@ -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.
Expand Down Expand Up @@ -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:
Expand Down
41 changes: 23 additions & 18 deletions tests/test_ghostly.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down Expand Up @@ -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.
Expand All @@ -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.
Expand Down Expand Up @@ -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 = [
Expand All @@ -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
Expand Down