From 365584b816e211855babc3e5730bee220937d356 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 10 Nov 2026 15:03:53 +0000 Subject: [PATCH 01/20] Perceive bond orders on ghost-free end states rather than with ghosts attached. --- src/ghostly/_ghostly.py | 91 +++++++++++++++++++++++++++-------------- 1 file changed, 61 insertions(+), 30 deletions(-) diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index 300be26..2c11f4f 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -955,28 +955,24 @@ def _dual( suffix = "0" # Check if the bridge atom is in a ring using RDKit. - if is_lambda1: - end_state_mol = _morph.link_to_perturbed(mol) - else: - end_state_mol = _morph.link_to_reference(mol) - - from sire.convert import to_rdkit - try: - rdmol = to_rdkit(end_state_mol) - except Exception: - rdmol = None + rdmol, index = _to_rdkit_end_state(mol, is_lambda1) + bridge_rdatom = rdmol.GetAtomWithIdx(index[bridge.value()]) + except Exception as e: + _logger.warning( + f"Unable to check the ring membership and hybridisation of bridge " + f"atom {bridge.value()}: {e}" + ) + bridge_rdatom = None - bridge_in_ring = ( - rdmol is not None and rdmol.GetAtomWithIdx(bridge.value()).IsInRing() - ) + bridge_in_ring = bridge_rdatom is not None and bridge_rdatom.IsInRing() # Check sp2 hybridisation for non-ring bridges. bridge_is_sp2 = False - if rdmol is not None and not bridge_in_ring: + if bridge_rdatom is not None and not bridge_in_ring: from rdkit.Chem.rdchem import HybridizationType - hyb = rdmol.GetAtomWithIdx(bridge.value()).GetHybridization() + hyb = bridge_rdatom.GetHybridization() bridge_is_sp2 = hyb == HybridizationType.SP2 # Warn if sp2 stiffening is enabled: the ~30° strain from 90° vs natural @@ -1375,29 +1371,25 @@ def _triple( # Store the molecular info. info = mol.info() - # Link the molecule to the desired end state. if is_lambda1: - end_state_mol = _morph.link_to_perturbed(mol) mod_key = "lambda_1" suffix = "1" else: - end_state_mol = _morph.link_to_reference(mol) mod_key = "lambda_0" suffix = "0" from rdkit.Chem import HybridizationType - from sire.convert import to_rdkit # Convert the molecule to RDKit to determine the hybridisation of the bridge atom. try: - rdmol = to_rdkit(end_state_mol) + rdmol, index = _to_rdkit_end_state(mol, is_lambda1) except Exception as e: msg = "Failed to convert molecule to RDKit for hybridisation check: " + str(e) _logger.error(msg) raise RuntimeError(msg) # Get the hybridisation and ring membership of the bridge atom. - bridge_rdatom = rdmol.GetAtomWithIdx(bridge.value()) + bridge_rdatom = rdmol.GetAtomWithIdx(index[bridge.value()]) hybridisation = bridge_rdatom.GetHybridization() bridge_in_ring = bridge_rdatom.IsInRing() @@ -2270,11 +2262,9 @@ def _linearise_ring_break( # Build an RDKit molecule for the physical end state so we can query # chirality directly rather than inferring it from substituent counts. - from sire.convert import to_rdkit as _to_rdkit from rdkit.Chem import ChiralType as _ChiralType - phys_mol = _physical_state(mol, is_lambda1) - rdmol = _to_rdkit(phys_mol) + rdmol, index = _to_rdkit_end_state(mol, not is_lambda1) linearised = set() for ghost in ghosts: @@ -2286,7 +2276,7 @@ def _linearise_ring_break( # Check whether this atom is a chiral centre in the physical end # state using RDKit. If so, linearising to 180° would destroy the # geometry, so skip it and warn. - rdatom = rdmol.GetAtomWithIdx(ghost.value()) + rdatom = rdmol.GetAtomWithIdx(index[ghost.value()]) chiral = rdatom.GetChiralTag() if chiral in ( _ChiralType.CHI_TETRAHEDRAL_CW, @@ -2780,24 +2770,29 @@ def _find_rotatable_bonds(ref, bridges, physical, is_lambda1=False, warn=True): return rotatable_bonds from rdkit.Chem import HybridizationType - from sire.convert import to_rdkit + # The ghosts are physical at the opposite end state. try: - rdmol = to_rdkit(ref) + rdmol, index = _to_rdkit_end_state(ref, not is_lambda1) except Exception as e: _logger.warning(f"Failed to convert molecule to RDKit for rotamer check: {e}") return rotatable_bonds + connectivity = ref.property("connectivity") + for bridge in bridges: for p in physical[bridge]: b_idx = bridge.value() p_idx = p.value() - rd_bond = rdmol.GetBondBetweenAtoms(p_idx, b_idx) - if rd_bond is None or rd_bond.IsInRing(): + # Ring membership from the full connectivity, since the physical + # neighbour may be a ghost at the opposite end state. + if not connectivity.are_bonded(p, bridge) or connectivity.in_ring( + p, bridge + ): continue - hybridisation = rdmol.GetAtomWithIdx(b_idx).GetHybridization() + hybridisation = rdmol.GetAtomWithIdx(index[b_idx]).GetHybridization() if hybridisation in ( HybridizationType.SP3, @@ -2816,6 +2811,42 @@ def _find_rotatable_bonds(ref, bridges, physical, is_lambda1=False, warn=True): return rotatable_bonds +def _to_rdkit_end_state(mol, is_lambda1=False): + """ + Return an RDKit molecule for one end state of a perturbable molecule, without + that end state's ghost atoms, and a dictionary mapping atom index values in + 'mol' to atom indices in the RDKit molecule. Ghost atoms are removed since + their bonds can exceed the valence of the atoms they attach to, which stops + RDKit perceiving bond orders, aromaticity and hybridisation. + """ + from sire.convert import to_rdkit + + suffix = "1" if is_lambda1 else "0" + physical = [ + i + for i, element in enumerate(mol.property(f"element{suffix}")) + if element.num_protons() > 0 + ] + + if is_lambda1: + state = _morph.extract_perturbed(mol) + else: + state = _morph.extract_reference(mol) + + rdmol = to_rdkit(state) + + # Atoms are left with unpaired electrons if bond orders couldn't be inferred. + if rdmol.GetNumAtoms() != len(physical) or any( + atom.GetNumRadicalElectrons() > 0 for atom in rdmol.GetAtoms() + ): + raise RuntimeError( + f"RDKit couldn't perceive the bond orders of the molecule at " + f"{_lam_sym} = {int(is_lambda1)}" + ) + + return rdmol, {index: i for i, index in enumerate(physical)} + + def _physical_state(mol, is_lambda1=False): """ Return the molecule linked to the end state where the ghost atoms of the From faa36f22ca42c998678f4205d4ea01af72f62293 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 10 Nov 2026 15:05:04 +0000 Subject: [PATCH 02/20] Treat sp2 bridges without a double or aromatic bond as sp3. --- src/ghostly/_ghostly.py | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index 2c11f4f..58dd6a5 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -972,7 +972,7 @@ def _dual( if bridge_rdatom is not None and not bridge_in_ring: from rdkit.Chem.rdchem import HybridizationType - hyb = bridge_rdatom.GetHybridization() + hyb = _hybridisation(bridge_rdatom) bridge_is_sp2 = hyb == HybridizationType.SP2 # Warn if sp2 stiffening is enabled: the ~30° strain from 90° vs natural @@ -1390,7 +1390,7 @@ def _triple( # Get the hybridisation and ring membership of the bridge atom. bridge_rdatom = rdmol.GetAtomWithIdx(index[bridge.value()]) - hybridisation = bridge_rdatom.GetHybridization() + hybridisation = _hybridisation(bridge_rdatom) bridge_in_ring = bridge_rdatom.IsInRing() # sp2 status for non-ring bridges (ring bridges use bridge_in_ring). @@ -2792,7 +2792,7 @@ def _find_rotatable_bonds(ref, bridges, physical, is_lambda1=False, warn=True): ): continue - hybridisation = rdmol.GetAtomWithIdx(index[b_idx]).GetHybridization() + hybridisation = _hybridisation(rdmol.GetAtomWithIdx(index[b_idx])) if hybridisation in ( HybridizationType.SP3, @@ -2847,6 +2847,23 @@ def _to_rdkit_end_state(mol, is_lambda1=False): return rdmol, {index: i for i, index in enumerate(physical)} +def _hybridisation(rdatom): + """ + Return the hybridisation of an RDKit atom, treating sp2 atoms without a + double or aromatic bond, i.e. conjugated lone pairs, as sp3. + """ + from rdkit.Chem import BondType, HybridizationType + + hybridisation = rdatom.GetHybridization() + if hybridisation == HybridizationType.SP2 and not any( + bond.GetIsAromatic() + or bond.GetBondType() in (BondType.DOUBLE, BondType.AROMATIC) + for bond in rdatom.GetBonds() + ): + return HybridizationType.SP3 + return hybridisation + + def _physical_state(mol, is_lambda1=False): """ Return the molecule linked to the end state where the ghost atoms of the From 252e2f5050def79d574df5c3e02aacbf7e9f9607 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 10 Nov 2026 15:49:29 +0000 Subject: [PATCH 03/20] Use angle terms for triple junction planarity. --- src/ghostly/_ghostly.py | 127 ++++++++++++++++++++++++++-------------- 1 file changed, 83 insertions(+), 44 deletions(-) diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index 58dd6a5..b557304 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -45,6 +45,9 @@ del _sys +# Triple junctions needing more than this energy (kcal/mol) to flatten are non-planar. +_MAX_PLANAR_FLATTENING = 2.0 + def modify( system, @@ -954,26 +957,22 @@ def _dual( mod_key = "lambda_0" suffix = "0" - # Check if the bridge atom is in a ring using RDKit. - try: - rdmol, index = _to_rdkit_end_state(mol, is_lambda1) - bridge_rdatom = rdmol.GetAtomWithIdx(index[bridge.value()]) - except Exception as e: - _logger.warning( - f"Unable to check the ring membership and hybridisation of bridge " - f"atom {bridge.value()}: {e}" - ) - bridge_rdatom = None - - bridge_in_ring = bridge_rdatom is not None and bridge_rdatom.IsInRing() + bridge_in_ring = connectivity.in_ring(bridge) - # Check sp2 hybridisation for non-ring bridges. + # Check sp2 hybridisation for non-ring bridges using RDKit. bridge_is_sp2 = False - if bridge_rdatom is not None and not bridge_in_ring: + if not bridge_in_ring: from rdkit.Chem.rdchem import HybridizationType - hyb = _hybridisation(bridge_rdatom) - bridge_is_sp2 = hyb == HybridizationType.SP2 + try: + rdmol, index = _to_rdkit_end_state(mol, is_lambda1) + hyb = _hybridisation(rdmol.GetAtomWithIdx(index[bridge.value()])) + bridge_is_sp2 = hyb == HybridizationType.SP2 + except Exception as e: + _logger.warning( + f"Unable to check the hybridisation of bridge atom " + f"{bridge.value()}: {e}" + ) # Warn if sp2 stiffening is enabled: the ~30° strain from 90° vs natural # ~120° may cause instabilities or free energy bias. @@ -1378,30 +1377,28 @@ def _triple( mod_key = "lambda_0" suffix = "0" - from rdkit.Chem import HybridizationType - - # Convert the molecule to RDKit to determine the hybridisation of the bridge atom. - try: - rdmol, index = _to_rdkit_end_state(mol, is_lambda1) - except Exception as e: - msg = "Failed to convert molecule to RDKit for hybridisation check: " + str(e) - _logger.error(msg) - raise RuntimeError(msg) + # The junction is planar if the angle terms between the physical neighbours + # of the bridge atom allow them to flatten at thermal energies. Higher order + # junctions pass four or more physical neighbours, so are never planar. + if len(physical) > 3: + is_planar = False + else: + flattening = _flattening_energy(mol, bridge, physical, is_lambda1) + if flattening is None: + _logger.warning( + f"Missing angle terms between the physical neighbours of bridge " + f"atom {bridge.value()} at {_lam_sym} = {int(is_lambda1)}. " + f"Defaulting to planar junction." + ) + is_planar = True + else: + _logger.debug(f" Flattening energy: {flattening:.3f} kcal/mol") + is_planar = flattening < _MAX_PLANAR_FLATTENING - # Get the hybridisation and ring membership of the bridge atom. - bridge_rdatom = rdmol.GetAtomWithIdx(index[bridge.value()]) - hybridisation = _hybridisation(bridge_rdatom) - bridge_in_ring = bridge_rdatom.IsInRing() + bridge_in_ring = connectivity.in_ring(bridge) # sp2 status for non-ring bridges (ring bridges use bridge_in_ring). - bridge_is_sp2 = not bridge_in_ring and hybridisation == HybridizationType.SP2 - - # Warn if the hybridisation is unspecified. - if hybridisation in (HybridizationType.UNSPECIFIED, HybridizationType.OTHER): - _logger.warning( - f"Unspecified hybridisation for bridge atom {bridge.value()} " - f"at {_lam_sym} = {int(is_lambda1)}. Defaulting to planar junction." - ) + bridge_is_sp2 = is_planar and not bridge_in_ring # Score the physical neighbours. if bridge_indices is not None: @@ -1421,11 +1418,7 @@ def _is_heavy(atom_idx): heavy_phys = set(physical) # Non-planar junction. - if hybridisation in ( - HybridizationType.SP3, - HybridizationType.SP3D, - HybridizationType.SP3D2, - ): + if not is_planar: _logger.debug(" Non-planar junction.") # First, modify the force constants of the angle terms between the ghost @@ -2778,14 +2771,14 @@ def _find_rotatable_bonds(ref, bridges, physical, is_lambda1=False, warn=True): _logger.warning(f"Failed to convert molecule to RDKit for rotamer check: {e}") return rotatable_bonds - connectivity = ref.property("connectivity") + connectivity = _end_state_connectivity(ref, is_lambda1) for bridge in bridges: for p in physical[bridge]: b_idx = bridge.value() p_idx = p.value() - # Ring membership from the full connectivity, since the physical + # Ring membership from the connectivity, since the physical # neighbour may be a ghost at the opposite end state. if not connectivity.are_bonded(p, bridge) or connectivity.in_ring( p, bridge @@ -2847,6 +2840,52 @@ def _to_rdkit_end_state(mol, is_lambda1=False): return rdmol, {index: i for i, index in enumerate(physical)} +def _end_state_connectivity(mol, is_lambda1=False): + """ + Return the connectivity of a perturbable molecule at an end state, including + the bonds to ghost atoms. Only ring-breaking molecules have end state + connectivities that differ. + """ + key = "connectivity1" if is_lambda1 else "connectivity0" + if mol.has_property(key): + return mol.property(key) + return mol.property("connectivity") + + +def _flattening_energy(mol, bridge, atoms, is_lambda1=False): + """ + Return the energy in kcal/mol needed to make the bridge atom and three of + its neighbours planar, starting from the minimum of the angle terms between + them, or None if any of those angle terms is missing. + """ + from math import pi + + from sire.legacy.CAS import Symbol + + theta = Symbol("theta") + info = mol.info() + suffix = "1" if is_lambda1 else "0" + neighbours = {atom.value() for atom in atoms} + + terms = {} + for p in mol.property("angle" + suffix).potentials(): + idx0 = info.atom_idx(p.atom0()).value() + idx1 = info.atom_idx(p.atom1()).value() + idx2 = info.atom_idx(p.atom2()).value() + if idx1 == bridge.value() and idx0 != idx2 and {idx0, idx2} <= neighbours: + terms[frozenset((idx0, idx2))] = _SireMM.AmberAngle(p.function(), theta) + + if len(terms) != 3: + return None + + # A planar centre needs angles summing to 2 pi. Minimising the harmonic + # energy subject to that constraint gives shortfall^2 / sum(1 / k). + if any(term.k() <= 0 for term in terms.values()): + return 0.0 + shortfall = max(0.0, 2 * pi - sum(term.theta0() for term in terms.values())) + return shortfall**2 / sum(1 / term.k() for term in terms.values()) + + def _hybridisation(rdatom): """ Return the hybridisation of an RDKit atom, treating sp2 atoms without a From 6d99d1a3e6b4b56f13b3ba507f58b336603667e2 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 10 Nov 2026 16:29:01 +0000 Subject: [PATCH 04/20] Keep ring bridge impropers. --- CHANGELOG.md | 2 ++ README.md | 11 ++++---- src/ghostly/_cli.py | 15 ----------- src/ghostly/_ghostly.py | 60 ++++++++++++++++++++--------------------- tests/test_ghostly.py | 17 +++++++++--- 5 files changed, 51 insertions(+), 54 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f29cd59..67677c4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,8 @@ Changelog * Added ``--single-well-anchors`` and ``--k-single-well`` to replace surviving anchor dihedrals with a single-well cosine potential whose minimum is the geometry of the end state where the ghost is physical. Multi-well anchors, e.g. the n=2 term of a dummy hydrogen on a carbonyl oxygen, can trap ghosts in an alternative well at intermediate lambda, giving hysteresis between independently equilibrated windows. Anchors with zero force constants or an ill-defined reference dihedral are left unchanged. (Default: False, 5 kcal/mol) [#19](https://github.com/OpenBioSim/ghostly/pull/19). * Record the options used for each modification in a ``ghostly_options`` property on each modified molecule, and raise an error when modifying a system that has already been modified, unless ``--force`` is used. Forced modifications are applied on top of the existing ones [#20](https://github.com/OpenBioSim/ghostly/pull/20). * Added support for Python 3.13 and 3.14, and dropped support for Python 3.10 and 3.11. +* Fixed bond order perception for end states with ghost atoms attached, which could fail and leave every heavy atom typed as sp3, and decide triple junction planarity and ring membership from the end state angle terms and connectivity [#23](https://github.com/OpenBioSim/ghostly/pull/23). +* Keep improper dihedrals that hold a single ghost in the plane of a ring bridge, since nothing else stops the ghost crossing the plane when the angles aren't stiffened, and remove the ``--stiffen-ring-bridges`` option [#23](https://github.com/OpenBioSim/ghostly/pull/23). [2026.2.0](https://github.com/openbiosim/ghostly/compare/2026.1.0...2026.2.0) - Sep 2026 ---------------------------------------------------------------------------------------- diff --git a/README.md b/README.md index 214002f..943ddba 100644 --- a/README.md +++ b/README.md @@ -28,12 +28,13 @@ handle the diversity of perturbations encountered in practice: - **Anchor selection scoring:** physical anchor atoms are scored to avoid transmuting or bridge atoms, preventing geometrically inconsistent constraints. -- **Ring and sp2 bridge handling:** angle stiffening is skipped by default for - ring and sp2 bridges, where local geometry already constrains the ghost and - 90° stiffening would introduce significant strain. It can be re-enabled via - `--stiffen-ring-bridges` and `--stiffen-sp2-bridges`. +- **Ring and sp2 bridge handling:** ring bridges are not stiffened. Instead, + their improper dihedrals with a single ghost are kept, holding the ghost in + the plane of the ring, where the rigid ring keeps the coupling acceptable. + Angle stiffening is skipped by default for sp2 bridges, and can be enabled + via `--stiffen-sp2-bridges`. - **Residual term cleanup:** a post-processing pass removes mixed improper - dihedrals and cross-bridge dihedrals missed by the per-bridge junction + dihedrals (other than those at ring bridges) and cross-bridge dihedrals missed by the per-bridge junction handlers, as well as angles where a ghost atom is the central atom and both terminal atoms are physical. - **Mixed dihedral softening:** surviving mixed ghost/physical dihedrals can diff --git a/src/ghostly/_cli.py b/src/ghostly/_cli.py index 7d4cdff..a733686 100644 --- a/src/ghostly/_cli.py +++ b/src/ghostly/_cli.py @@ -202,20 +202,6 @@ def run(): required=False, ) - parser.add_argument( - "--stiffen-ring-bridges", - action="store_true", - help=""" - Apply 90 degree angle stiffening when the bridge atom is in a - ring. By default this is skipped since the ring geometry already - constrains the ghost position and stiffening introduces strain. - A warning is logged when this is enabled and a ring bridge is - detected. - """, - default=False, - required=False, - ) - parser.add_argument( "--stiffen-sp2-bridges", action="store_true", @@ -434,7 +420,6 @@ def run(): k_single_well=k_single_well.value(), stiffen_rotamers=args.stiffen_rotamers, k_rotamer=k_rotamer.value(), - stiffen_ring_bridges=args.stiffen_ring_bridges, stiffen_sp2_bridges=args.stiffen_sp2_bridges, linearise_ring_break=args.linearise_ring_break, force=args.force, diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index b557304..38225e2 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -60,7 +60,6 @@ def modify( k_single_well=5, stiffen_rotamers=False, k_rotamer=50, - stiffen_ring_bridges=False, stiffen_sp2_bridges=False, linearise_ring_break=False, force=False, @@ -154,13 +153,6 @@ def modify( rotamer anchor dihedrals (in kcal/mol). The resulting barrier height is ``2 * k_rotamer``. Only used when ``stiffen_rotamers`` is True. - stiffen_ring_bridges : bool, optional - Whether to apply 90° angle stiffening when the bridge atom is in a - ring. Boresch notes that for rigid ring systems the partition function - coupling is small and acceptable, so stiffening is skipped by default - (False). Enabling this restores the original behaviour but introduces - strain where the 90° target fights the natural ring geometry. - stiffen_sp2_bridges : bool, optional Whether to apply 90° angle stiffening when the bridge atom is sp2 and not in a ring (e.g. carbonyl or alkene carbons). The natural angle is @@ -242,7 +234,6 @@ def modify( "k_single_well": k_single_well, "stiffen_rotamers": stiffen_rotamers, "k_rotamer": k_rotamer, - "stiffen_ring_bridges": stiffen_ring_bridges, "stiffen_sp2_bridges": stiffen_sp2_bridges, "linearise_ring_break": linearise_ring_break, } @@ -426,7 +417,6 @@ def modify( modifications, k_hard=k_hard, bridge_indices=bridge_indices0, - stiffen_ring_bridges=stiffen_ring_bridges, stiffen_sp2_bridges=stiffen_sp2_bridges, ) @@ -444,7 +434,6 @@ def modify( optimise_angles=optimise_angles, num_optimise=num_optimise, bridge_indices=bridge_indices0, - stiffen_ring_bridges=stiffen_ring_bridges, stiffen_sp2_bridges=stiffen_sp2_bridges, ) @@ -534,7 +523,6 @@ def modify( k_hard=k_hard, is_lambda1=True, bridge_indices=bridge_indices1, - stiffen_ring_bridges=stiffen_ring_bridges, stiffen_sp2_bridges=stiffen_sp2_bridges, ) @@ -552,7 +540,6 @@ def modify( num_optimise=num_optimise, is_lambda1=True, bridge_indices=bridge_indices1, - stiffen_ring_bridges=stiffen_ring_bridges, stiffen_sp2_bridges=stiffen_sp2_bridges, ) @@ -881,7 +868,6 @@ def _dual( k_hard=100, is_lambda1=False, bridge_indices=None, - stiffen_ring_bridges=False, stiffen_sp2_bridges=False, ): r""" @@ -1068,11 +1054,10 @@ def _is_heavy(atom_idx): # Identify the physical atom in this angle. phys_atom = idx2 if idx0 in ghosts else idx0 - # Ring bridges: the ring geometry already constrains the ghost - # position and prevents flapping (Boresch et al. note this is - # "acceptable if the physical molecule is a rigid ring system"). - # Stiffening to 90° would introduce strain without benefit. - if bridge_in_ring and not stiffen_ring_bridges: + # Ring bridges keep their ghost impropers, which hold the ghost + # in the plane. The rigid ring keeps the redundant coupling + # acceptable (Boresch et al.), so the angles are not stiffened. + if bridge_in_ring: new_angles.set(idx0, idx1, idx2, p.function()) _logger.debug( f" Skipping stiffening for angle " @@ -1164,7 +1149,7 @@ def _is_heavy(atom_idx): # 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: + if bridge_in_ring: preserve_reason = "is in a ring" elif bridge_is_sp2 and not stiffen_sp2_bridges: preserve_reason = "is sp2" @@ -1267,7 +1252,6 @@ def _is_heavy(atom_idx): k_hard=k_hard, is_lambda1=is_lambda1, bridge_indices=bridge_indices, - stiffen_ring_bridges=stiffen_ring_bridges, stiffen_sp2_bridges=stiffen_sp2_bridges, ) @@ -1288,7 +1272,6 @@ def _triple( num_optimise=10, is_lambda1=False, bridge_indices=None, - stiffen_ring_bridges=False, stiffen_sp2_bridges=False, ): r""" @@ -1698,9 +1681,7 @@ def _is_heavy(atom_idx): # When stiffening is skipped for ring/sp2 bridges, preserve all # ghost-containing angles at their original values rather than # removing the sacrificial one. - if (bridge_in_ring and not stiffen_ring_bridges) or ( - bridge_is_sp2 and not stiffen_sp2_bridges - ): + if bridge_in_ring or (bridge_is_sp2 and not stiffen_sp2_bridges): new_angles.set(idx0, idx1, idx2, p.function()) _logger.debug( f" Preserving angle " @@ -1760,7 +1741,6 @@ def _is_heavy(atom_idx): k_hard=k_hard, is_lambda1=is_lambda1, bridge_indices=bridge_indices, - stiffen_ring_bridges=stiffen_ring_bridges, stiffen_sp2_bridges=stiffen_sp2_bridges, ) @@ -1970,7 +1950,9 @@ def _higher( def _remove_impropers(mol, ghosts, modifications, is_lambda1=False): """ - Remove improper dihedral terms that bridge the ghost and physical systems. + Remove improper dihedral terms that bridge the ghost and physical systems, + other than those centred on a physical ring atom with a single ghost. These + hold the ghost in the plane of a ring bridge, which is not stiffened. Parameters ---------- @@ -2008,17 +1990,35 @@ def _remove_impropers(mol, ghosts, modifications, is_lambda1=False): # Initialise a container to store the updated bonded terms. new_impropers = _SireMM.FourAtomFunctions(mol.info()) + connectivity = _end_state_connectivity(mol, is_lambda1) + # Loop over the improper dihedrals. for p in impropers.potentials(): idx0 = info.atom_idx(p.atom0()) idx1 = info.atom_idx(p.atom1()) idx2 = info.atom_idx(p.atom2()) idx3 = info.atom_idx(p.atom3()) + idxs = (idx0, idx1, idx2, idx3) + num_ghosts = sum(idx in ghosts for idx in idxs) + + # The central atom is bonded to the other three. + centre = next( + ( + idx + for idx in idxs + if all(connectivity.are_bonded(idx, x) for x in idxs if x != idx) + ), + None, + ) + ring_bridge = ( + num_ghosts == 1 + and centre is not None + and centre not in ghosts + and connectivity.in_ring(centre) + ) # Remove any improper dihedrals that bridge the ghost and physical systems. - if not all(idx in ghosts for idx in (idx0, idx1, idx2, idx3)) and any( - idx in ghosts for idx in (idx0, idx1, idx2, idx3) - ): + if 0 < num_ghosts < len(idxs) and not ring_bridge: _logger.debug( f" Removing improper dihedral: [{idx0.value()}-{idx1.value()}-{idx2.value()}-{idx3.value()}], {p.function()}" ) diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index 55e40bf..5a505ca 100644 --- a/tests/test_ghostly.py +++ b/tests/test_ghostly.py @@ -124,10 +124,10 @@ def test_toluene_to_pyridine(): for dihedral in missing_dihedrals ) - # Bridge atom 1 is in an aromatic ring. Boresch notes that for rigid ring - # systems the coupling is acceptable, and the ring geometry already - # constrains the ghost position, so no stiffening is applied. The angles - # involving the ring-bridge ghost should be unchanged. + # Bridge atom 1 is in an aromatic ring, so no stiffening is applied. Its + # improper holds the ghost in the plane of the ring, and the rigid ring + # keeps the coupling acceptable (Boresch et al.). The angles involving the + # ring-bridge ghost should be unchanged. ring_bridge_angles = [ (AtomIdx(0), AtomIdx(1), AtomIdx(2)), (AtomIdx(0), AtomIdx(1), AtomIdx(6)), @@ -154,6 +154,15 @@ def test_toluene_to_pyridine(): if angle_idx in orig_angle_map and angle_idx in new_angle_map: assert orig_angle_map[angle_idx] == new_angle_map[angle_idx] + # The improper centred on the ring bridge with the ghost should be kept. + ring_bridge_improper = (AtomIdx(0), AtomIdx(1), AtomIdx(6), AtomIdx(2)) + assert check_improper( + info, mols[0].property("improper1").potentials(), *ring_bridge_improper + ) + assert check_improper( + info, new_mols[0].property("improper1").potentials(), *ring_bridge_improper + ) + def test_acetone_to_propenol(): """ From 7fc906dacbded24292f08e6112a09397c4563668 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 10 Nov 2026 16:30:32 +0000 Subject: [PATCH 05/20] Stiffen all non-ring bridges. --- CHANGELOG.md | 1 + README.md | 17 +++++----- src/ghostly/_cli.py | 15 --------- src/ghostly/_ghostly.py | 71 +++-------------------------------------- tests/test_ghostly.py | 66 ++++++++++++-------------------------- 5 files changed, 35 insertions(+), 135 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 67677c4..beb4de0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ Changelog * Added support for Python 3.13 and 3.14, and dropped support for Python 3.10 and 3.11. * Fixed bond order perception for end states with ghost atoms attached, which could fail and leave every heavy atom typed as sp3, and decide triple junction planarity and ring membership from the end state angle terms and connectivity [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Keep improper dihedrals that hold a single ghost in the plane of a ring bridge, since nothing else stops the ghost crossing the plane when the angles aren't stiffened, and remove the ``--stiffen-ring-bridges`` option [#23](https://github.com/OpenBioSim/ghostly/pull/23). +* Always stiffen ghost angles at bridges that aren't in a ring, and remove the ``--stiffen-sp2-bridges`` option. Without impropers, the angle terms of an sp2 bridge don't hold a ghost in the plane, so skipping stiffening let it cross freely [#23](https://github.com/OpenBioSim/ghostly/pull/23). [2026.2.0](https://github.com/openbiosim/ghostly/compare/2026.1.0...2026.2.0) - Sep 2026 ---------------------------------------------------------------------------------------- diff --git a/README.md b/README.md index 943ddba..e32bec7 100644 --- a/README.md +++ b/README.md @@ -28,15 +28,16 @@ handle the diversity of perturbations encountered in practice: - **Anchor selection scoring:** physical anchor atoms are scored to avoid transmuting or bridge atoms, preventing geometrically inconsistent constraints. -- **Ring and sp2 bridge handling:** ring bridges are not stiffened. Instead, - their improper dihedrals with a single ghost are kept, holding the ghost in - the plane of the ring, where the rigid ring keeps the coupling acceptable. - Angle stiffening is skipped by default for sp2 bridges, and can be enabled - via `--stiffen-sp2-bridges`. +- **Ring bridge handling:** ring bridges are not stiffened. Instead, their + improper dihedrals with a single ghost are kept, holding the ghost in the + plane of the ring, where the rigid ring keeps the coupling acceptable. All + other bridges are stiffened. +- **Junction planarity:** whether a triple junction is planar is decided from + the energy needed to flatten the angle terms between its physical neighbours. - **Residual term cleanup:** a post-processing pass removes mixed improper - dihedrals (other than those at ring bridges) and cross-bridge dihedrals missed by the per-bridge junction - handlers, as well as angles where a ghost atom is the central atom and both - terminal atoms are physical. + dihedrals (other than those at ring bridges) and cross-bridge dihedrals + missed by the per-bridge junction handlers, as well as angles where a ghost + atom is the central atom and both terminal atoms are physical. - **Mixed dihedral softening:** surviving mixed ghost/physical dihedrals can be softened via `--soften-anchors` to allow ghost groups to reorient and avoid steric clashes at small λ. diff --git a/src/ghostly/_cli.py b/src/ghostly/_cli.py index a733686..08920e2 100644 --- a/src/ghostly/_cli.py +++ b/src/ghostly/_cli.py @@ -202,20 +202,6 @@ def run(): required=False, ) - parser.add_argument( - "--stiffen-sp2-bridges", - action="store_true", - help=""" - Apply 90 degree angle stiffening when the bridge atom is sp2 and - not in a ring. By default this is skipped since sp2 hybridisation - already constrains the ghost to the molecular plane and stiffening - introduces ~30 degrees of strain. A warning is logged when this - is enabled and an sp2 bridge is detected. - """, - default=False, - required=False, - ) - parser.add_argument( "--linearise-ring-break", action=argparse.BooleanOptionalAction, @@ -420,7 +406,6 @@ def run(): k_single_well=k_single_well.value(), stiffen_rotamers=args.stiffen_rotamers, k_rotamer=k_rotamer.value(), - stiffen_sp2_bridges=args.stiffen_sp2_bridges, linearise_ring_break=args.linearise_ring_break, force=args.force, ) diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index 38225e2..e90503f 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -60,7 +60,6 @@ def modify( k_single_well=5, stiffen_rotamers=False, k_rotamer=50, - stiffen_sp2_bridges=False, linearise_ring_break=False, force=False, ): @@ -153,14 +152,6 @@ def modify( rotamer anchor dihedrals (in kcal/mol). The resulting barrier height is ``2 * k_rotamer``. Only used when ``stiffen_rotamers`` is True. - stiffen_sp2_bridges : bool, optional - Whether to apply 90° angle stiffening when the bridge atom is sp2 and - not in a ring (e.g. carbonyl or alkene carbons). The natural angle is - ~120°, so stiffening to 90° introduces ~30° strain and may cause - instabilities or free energy bias. Stiffening is skipped by default - (False). Enabling this restores the original behaviour but a warning - will be logged to flag the potential strain issue. - linearise_ring_break : bool, optional Whether to apply a linear spacer modification to ghost atoms that bridge two physical atoms (the ring-breaking topology). Instead of @@ -234,7 +225,6 @@ def modify( "k_single_well": k_single_well, "stiffen_rotamers": stiffen_rotamers, "k_rotamer": k_rotamer, - "stiffen_sp2_bridges": stiffen_sp2_bridges, "linearise_ring_break": linearise_ring_break, } @@ -417,7 +407,6 @@ def modify( modifications, k_hard=k_hard, bridge_indices=bridge_indices0, - stiffen_sp2_bridges=stiffen_sp2_bridges, ) # Triple junction. @@ -434,7 +423,6 @@ def modify( optimise_angles=optimise_angles, num_optimise=num_optimise, bridge_indices=bridge_indices0, - stiffen_sp2_bridges=stiffen_sp2_bridges, ) # Higher order junction. @@ -523,7 +511,6 @@ def modify( k_hard=k_hard, is_lambda1=True, bridge_indices=bridge_indices1, - stiffen_sp2_bridges=stiffen_sp2_bridges, ) elif junction == 3: @@ -540,7 +527,6 @@ def modify( num_optimise=num_optimise, is_lambda1=True, bridge_indices=bridge_indices1, - stiffen_sp2_bridges=stiffen_sp2_bridges, ) # Higher order junction. @@ -868,7 +854,6 @@ def _dual( k_hard=100, is_lambda1=False, bridge_indices=None, - stiffen_sp2_bridges=False, ): r""" Apply modifications to a dual junction. @@ -945,31 +930,6 @@ def _dual( bridge_in_ring = connectivity.in_ring(bridge) - # Check sp2 hybridisation for non-ring bridges using RDKit. - bridge_is_sp2 = False - if not bridge_in_ring: - from rdkit.Chem.rdchem import HybridizationType - - try: - rdmol, index = _to_rdkit_end_state(mol, is_lambda1) - hyb = _hybridisation(rdmol.GetAtomWithIdx(index[bridge.value()])) - bridge_is_sp2 = hyb == HybridizationType.SP2 - except Exception as e: - _logger.warning( - f"Unable to check the hybridisation of bridge atom " - f"{bridge.value()}: {e}" - ) - - # Warn if sp2 stiffening is enabled: the ~30° strain from 90° vs natural - # ~120° may cause instabilities or free energy bias. - if bridge_is_sp2 and stiffen_sp2_bridges: - _logger.warning( - f"Bridge atom {bridge.value()} is sp2 and not in a ring. " - f"Stiffening ghost angles to 90° introduces ~30° strain " - f"from the natural ~120° sp2 geometry. This may cause " - f"instabilities or free energy bias." - ) - # Score the physical neighbours. if bridge_indices is not None: phys_scores = {p: _score_atom(mol, p, bridge_indices) for p in physical} @@ -1066,18 +1026,6 @@ def _is_heavy(atom_idx): ) continue - # sp2 non-ring bridges: sp2 hybridisation constrains the ghost - # to the molecular plane, preventing flapping. Stiffening to - # 90° from the natural ~120° introduces ~30° strain. - if bridge_is_sp2 and not stiffen_sp2_bridges: - new_angles.set(idx0, idx1, idx2, p.function()) - _logger.debug( - f" Skipping stiffening for angle " - f"[{idx0.value()}-{idx1.value()}-{idx2.value()}]: " - f"bridge atom {bridge.value()} is sp2." - ) - continue - # Skip stiffening through poorly-scoring atoms if better ones # exist AND at least one heavy atom would remain as anchor. # Without a heavy anchor, the ghost is under-constrained. @@ -1151,8 +1099,6 @@ def _is_heavy(atom_idx): # skips poorly-scoring physical atoms, so mirror that check here. if bridge_in_ring: 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 ): @@ -1252,7 +1198,6 @@ def _is_heavy(atom_idx): k_hard=k_hard, is_lambda1=is_lambda1, bridge_indices=bridge_indices, - stiffen_sp2_bridges=stiffen_sp2_bridges, ) # Return the updated molecule. @@ -1272,7 +1217,6 @@ def _triple( num_optimise=10, is_lambda1=False, bridge_indices=None, - stiffen_sp2_bridges=False, ): r""" Apply modifications to a triple junction. @@ -1380,9 +1324,6 @@ def _triple( bridge_in_ring = connectivity.in_ring(bridge) - # sp2 status for non-ring bridges (ring bridges use bridge_in_ring). - bridge_is_sp2 = is_planar and not bridge_in_ring - # Score the physical neighbours. if bridge_indices is not None: phys_scores = {p: _score_atom(mol, p, bridge_indices) for p in physical} @@ -1678,16 +1619,15 @@ def _is_heavy(atom_idx): idxs = [idx0, idx1, idx2] if idx in idxs and any([x in ghosts for x in idxs]): - # When stiffening is skipped for ring/sp2 bridges, preserve all - # ghost-containing angles at their original values rather than - # removing the sacrificial one. - if bridge_in_ring or (bridge_is_sp2 and not stiffen_sp2_bridges): + # Ring bridges aren't stiffened, so preserve all ghost-containing + # angles at their original values rather than removing the + # sacrificial one. + if bridge_in_ring: new_angles.set(idx0, idx1, idx2, p.function()) _logger.debug( f" Preserving 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()} is in a ring, " f"not stiffening." ) continue @@ -1741,7 +1681,6 @@ def _is_heavy(atom_idx): k_hard=k_hard, is_lambda1=is_lambda1, bridge_indices=bridge_indices, - stiffen_sp2_bridges=stiffen_sp2_bridges, ) # Return the updated molecule. diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index 5a505ca..78cf8ce 100644 --- a/tests/test_ghostly.py +++ b/tests/test_ghostly.py @@ -184,9 +184,8 @@ def test_acetone_to_propenol(): k_hard = 100 # Apply the ghost atom modifications. Bridge atom 3 is sp2 (carbonyl - # carbon): explicitly enable sp2 stiffening to test the canonical Boresch - # planar triple junction behaviour. - new_mols, _ = modify(mols, k_hard=k_hard, stiffen_sp2_bridges=True) + # carbon), giving the canonical Boresch planar triple junction. + new_mols, _ = modify(mols, k_hard=k_hard) # Get the new angles and dihedrals. new_angles0 = new_mols[0].property("angle0") @@ -402,9 +401,9 @@ def test_ejm49_to_ejm31(): 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. + # Bridge atom 17 is sp2 and not in a ring. Its angle terms don't hold the + # ghost in the plane once the impropers are removed, so the + # physical-bridge-ghost angles are stiffened to 90 degrees at both end states. sp2_bridge_angles0 = [ (AtomIdx(16), AtomIdx(17), AtomIdx(39)), (AtomIdx(18), AtomIdx(17), AtomIdx(39)), @@ -414,46 +413,21 @@ def test_ejm49_to_ejm31(): (AtomIdx(18), AtomIdx(17), AtomIdx(20)), ] - orig_angle_map0 = { - ( - info.atom_idx(p.atom0()), - info.atom_idx(p.atom1()), - info.atom_idx(p.atom2()), - ): str(p.function()) - for p in angles0.potentials() - } - new_angle_map0 = { - ( - info.atom_idx(p.atom0()), - info.atom_idx(p.atom1()), - info.atom_idx(p.atom2()), - ): str(p.function()) - for p in new_angles0.potentials() - } - orig_angle_map1 = { - ( - info.atom_idx(p.atom0()), - info.atom_idx(p.atom1()), - info.atom_idx(p.atom2()), - ): str(p.function()) - for p in angles1.potentials() - } - new_angle_map1 = { - ( - info.atom_idx(p.atom0()), - info.atom_idx(p.atom1()), - info.atom_idx(p.atom2()), - ): str(p.function()) - for p in new_angles1.potentials() - } - - for angle_idx in sp2_bridge_angles0: - if angle_idx in orig_angle_map0 and angle_idx in new_angle_map0: - assert orig_angle_map0[angle_idx] == new_angle_map0[angle_idx] - - for angle_idx in sp2_bridge_angles1: - if angle_idx in orig_angle_map1 and angle_idx in new_angle_map1: - assert orig_angle_map1[angle_idx] == new_angle_map1[angle_idx] + stiffened = "100 [theta - 1.5708]^2" + for new_angles, sp2_bridge_angles in ( + (new_angles0, sp2_bridge_angles0), + (new_angles1, sp2_bridge_angles1), + ): + new_angle_map = { + ( + info.atom_idx(p.atom0()), + info.atom_idx(p.atom1()), + info.atom_idx(p.atom2()), + ): str(p.function()) + for p in new_angles.potentials() + } + for angle_idx in sp2_bridge_angles: + assert new_angle_map[angle_idx] == stiffened # Create improper IDs for the missing impropers at lambda = 0. missing_impropers0 = [ From 74dc93ea69ee838ff9df943a9ae03e28f60a4eb7 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 10 Nov 2026 16:33:01 +0000 Subject: [PATCH 06/20] Detect rotatable bonds from angle terms. --- CHANGELOG.md | 1 + README.md | 12 +++-- src/ghostly/_cli.py | 2 +- src/ghostly/_ghostly.py | 105 +++++++++++++++++++++------------------- 4 files changed, 63 insertions(+), 57 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index beb4de0..60a6af3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,7 @@ Changelog * Fixed bond order perception for end states with ghost atoms attached, which could fail and leave every heavy atom typed as sp3, and decide triple junction planarity and ring membership from the end state angle terms and connectivity [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Keep improper dihedrals that hold a single ghost in the plane of a ring bridge, since nothing else stops the ghost crossing the plane when the angles aren't stiffened, and remove the ``--stiffen-ring-bridges`` option [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Always stiffen ghost angles at bridges that aren't in a ring, and remove the ``--stiffen-sp2-bridges`` option. Without impropers, the angle terms of an sp2 bridge don't hold a ghost in the plane, so skipping stiffening let it cross freely [#23](https://github.com/OpenBioSim/ghostly/pull/23). +* Detect rotatable bridge bonds from the angle terms of the bridge, rather than its RDKit hybridisation, so that rotamer detection is consistent with junction planarity [#23](https://github.com/OpenBioSim/ghostly/pull/23). [2026.2.0](https://github.com/openbiosim/ghostly/compare/2026.1.0...2026.2.0) - Sep 2026 ---------------------------------------------------------------------------------------- diff --git a/README.md b/README.md index e32bec7..cabc0a7 100644 --- a/README.md +++ b/README.md @@ -32,8 +32,9 @@ handle the diversity of perturbations encountered in practice: improper dihedrals with a single ghost are kept, holding the ghost in the plane of the ring, where the rigid ring keeps the coupling acceptable. All other bridges are stiffened. -- **Junction planarity:** whether a triple junction is planar is decided from - the energy needed to flatten the angle terms between its physical neighbours. +- **Junction planarity:** whether a triple junction, or a bridge for rotamer + detection, is planar is decided from the energy needed to flatten the angle + terms between its neighbours. - **Residual term cleanup:** a post-processing pass removes mixed improper dihedrals (other than those at ring bridges) and cross-bridge dihedrals missed by the per-bridge junction handlers, as well as angles where a ghost @@ -41,9 +42,10 @@ handle the diversity of perturbations encountered in practice: - **Mixed dihedral softening:** surviving mixed ghost/physical dihedrals can be softened via `--soften-anchors` to allow ghost groups to reorient and avoid steric clashes at small λ. -- **Rotamer stiffening:** `--stiffen-rotamers` replaces rotatable sp3 anchor - dihedrals with a stiff single-well cosine to control ghost orientation - through flexible bonds. +- **Rotamer stiffening:** `--stiffen-rotamers` replaces anchor dihedrals + across rotatable bonds, i.e. those not in a ring whose bridge atom is neither + planar nor linear, with a stiff single-well cosine to control ghost + orientation through flexible bonds. - **Ring-breaking perturbations:** adjacent bridges with independent ghost groups retain each other as physical neighbours; angles with a ghost central atom spanning two physical neighbours are replaced by a linear spacer diff --git a/src/ghostly/_cli.py b/src/ghostly/_cli.py index 08920e2..01ba431 100644 --- a/src/ghostly/_cli.py +++ b/src/ghostly/_cli.py @@ -182,7 +182,7 @@ def run(): help=""" Whether to replace rotamer anchor dihedrals with a stiff single-well cosine potential. When a bridge-physical bond is - rotatable (not in a ring, sp3 bridge), surviving anchor + rotatable (not in a ring, non-planar bridge), surviving anchor dihedrals can allow rotameric transitions of ghost atoms at intermediate lambda. """, diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index e90503f..657a64d 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -140,7 +140,7 @@ def modify( stiffen_rotamers : bool, optional Whether to replace rotamer anchor dihedrals with a stiff single-well cosine potential. When a bridge--physical bond is rotatable (not in a - ring, sp3 bridge), surviving anchor dihedrals can allow rotameric + ring, non-planar bridge), surviving anchor dihedrals can allow rotameric transitions of ghost atoms at intermediate lambda. Enabling this replaces those dihedrals with a single cosine well of depth ``2 * k_rotamer``, with its minimum at the dihedral angle of the end @@ -2465,7 +2465,7 @@ def _modify_mixed_dihedrals( 1. Removed, when ``soften_anchors`` is 0.0. 2. Stiffened, when ``stiffen_rotamers`` is True and the dihedral spans a - rotatable sp3 bridge--physical bond with a ghost terminal atom. + rotatable bridge--physical bond with a ghost terminal atom. 3. Replaced with a single well of ``k_single_well``, when ``single_well_anchors`` is True and the dihedral is a non-zero anchor, i.e. a single ghost terminal atom and three physical atoms. @@ -2479,7 +2479,7 @@ def _modify_mixed_dihedrals( state where the ghost is physical, and a barrier of 2k. When phi_ref is ill-defined the term falls through to the next action. - Rotatable sp3 bridge--physical bonds are always logged as warnings when + Rotatable bridge--physical bonds are always logged as warnings when ``stiffen_rotamers`` is False, since rotameric transitions of ghost atoms at intermediate lambda can degrade convergence (see Boresch et al., JCTC 2021). @@ -2658,13 +2658,14 @@ def _modify_mixed_dihedrals( def _find_rotatable_bonds(ref, bridges, physical, is_lambda1=False, warn=True): r""" Find the rotatable bridge--physical bonds, i.e. those that are not in a - ring and have an sp3 bridge atom, so that surviving anchor dihedrals - across them can allow rotameric transitions of ghost atoms. + ring and whose bridge atom is neither planar nor linear in the end state + where the ghosts are physical, so that surviving anchor dihedrals across + them can allow rotameric transitions of ghost atoms. R1 DR1 \ / \ / - X ---DR2 X is sp3, P-X is rotatable + X ---DR2 X is non-planar, P-X is rotatable / \ / \ R2 DR3 @@ -2701,22 +2702,14 @@ def _find_rotatable_bonds(ref, bridges, physical, is_lambda1=False, warn=True): if not bridges: return rotatable_bonds - from rdkit.Chem import HybridizationType - - # The ghosts are physical at the opposite end state. - try: - rdmol, index = _to_rdkit_end_state(ref, not is_lambda1) - except Exception as e: - _logger.warning(f"Failed to convert molecule to RDKit for rotamer check: {e}") - return rotatable_bonds - connectivity = _end_state_connectivity(ref, is_lambda1) for bridge in bridges: - for p in physical[bridge]: - b_idx = bridge.value() - p_idx = p.value() + # The ghosts are physical at the opposite end state. + if not _is_rotor(ref, bridge, not is_lambda1): + continue + for p in physical[bridge]: # Ring membership from the connectivity, since the physical # neighbour may be a ghost at the opposite end state. if not connectivity.are_bonded(p, bridge) or connectivity.in_ring( @@ -2724,21 +2717,14 @@ def _find_rotatable_bonds(ref, bridges, physical, is_lambda1=False, warn=True): ): continue - hybridisation = _hybridisation(rdmol.GetAtomWithIdx(index[b_idx])) - - if hybridisation in ( - HybridizationType.SP3, - HybridizationType.SP3D, - HybridizationType.SP3D2, - ): - rotatable_bonds.add((p, bridge)) - if warn: - _logger.warning( - f"Potential rotamer anchor at {_lam_sym} = {int(is_lambda1)}: " - f"bond {p_idx}-{b_idx} is a rotatable sp3 bond. " - f"Surviving anchor dihedrals may allow rotameric " - f"transitions of ghost atoms." - ) + rotatable_bonds.add((p, bridge)) + if warn: + _logger.warning( + f"Potential rotamer anchor at {_lam_sym} = {int(is_lambda1)}: " + f"bond {p.value()}-{bridge.value()} is a rotatable bond. " + f"Surviving anchor dihedrals may allow rotameric " + f"transitions of ghost atoms." + ) return rotatable_bonds @@ -2791,14 +2777,11 @@ def _end_state_connectivity(mol, is_lambda1=False): return mol.property("connectivity") -def _flattening_energy(mol, bridge, atoms, is_lambda1=False): +def _angle_terms(mol, bridge, atoms, is_lambda1=False): """ - Return the energy in kcal/mol needed to make the bridge atom and three of - its neighbours planar, starting from the minimum of the angle terms between - them, or None if any of those angle terms is missing. + Return the angle terms centred on the bridge atom between pairs of the + given atoms at an end state, keyed by the frozenset of the pair. """ - from math import pi - from sire.legacy.CAS import Symbol theta = Symbol("theta") @@ -2813,7 +2796,18 @@ def _flattening_energy(mol, bridge, atoms, is_lambda1=False): idx2 = info.atom_idx(p.atom2()).value() if idx1 == bridge.value() and idx0 != idx2 and {idx0, idx2} <= neighbours: terms[frozenset((idx0, idx2))] = _SireMM.AmberAngle(p.function(), theta) + return terms + + +def _flattening_energy(mol, bridge, atoms, is_lambda1=False): + """ + Return the energy in kcal/mol needed to make the bridge atom and three of + its neighbours planar, starting from the minimum of the angle terms between + them, or None if any of those angle terms is missing. + """ + from math import pi + terms = _angle_terms(mol, bridge, atoms, is_lambda1) if len(terms) != 3: return None @@ -2825,21 +2819,30 @@ def _flattening_energy(mol, bridge, atoms, is_lambda1=False): return shortfall**2 / sum(1 / term.k() for term in terms.values()) -def _hybridisation(rdatom): +def _is_rotor(mol, bridge, is_lambda1=False, max_angle=170): """ - Return the hybridisation of an RDKit atom, treating sp2 atoms without a - double or aromatic bond, i.e. conjugated lone pairs, as sp3. + Return whether rotation about the bonds of the bridge atom can move its + other substituents at an end state, i.e. the bridge is neither planar, by + its angle terms, nor linear. """ - from rdkit.Chem import BondType, HybridizationType + from math import degrees - hybridisation = rdatom.GetHybridization() - if hybridisation == HybridizationType.SP2 and not any( - bond.GetIsAromatic() - or bond.GetBondType() in (BondType.DOUBLE, BondType.AROMATIC) - for bond in rdatom.GetBonds() - ): - return HybridizationType.SP3 - return hybridisation + elements = mol.property("element1" if is_lambda1 else "element0") + neighbours = [ + n + for n in _end_state_connectivity(mol, is_lambda1).connections_to(bridge) + if elements[n.value()].num_protons() > 0 + ] + + if len(neighbours) < 2: + return False + if len(neighbours) == 2: + terms = list(_angle_terms(mol, bridge, neighbours, is_lambda1).values()) + return bool(terms) and degrees(terms[0].theta0()) < max_angle + if len(neighbours) == 3: + flattening = _flattening_energy(mol, bridge, neighbours, is_lambda1) + return flattening is not None and flattening >= _MAX_PLANAR_FLATTENING + return True def _physical_state(mol, is_lambda1=False): From a42ec06f3c998807d88a1bed8e08b57d7a17942e Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 10 Nov 2026 16:34:06 +0000 Subject: [PATCH 07/20] Add CHANGELOG entry for the perception fix. --- CHANGELOG.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 60a6af3..06d35fe 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,7 +10,8 @@ Changelog * Added ``--single-well-anchors`` and ``--k-single-well`` to replace surviving anchor dihedrals with a single-well cosine potential whose minimum is the geometry of the end state where the ghost is physical. Multi-well anchors, e.g. the n=2 term of a dummy hydrogen on a carbonyl oxygen, can trap ghosts in an alternative well at intermediate lambda, giving hysteresis between independently equilibrated windows. Anchors with zero force constants or an ill-defined reference dihedral are left unchanged. (Default: False, 5 kcal/mol) [#19](https://github.com/OpenBioSim/ghostly/pull/19). * Record the options used for each modification in a ``ghostly_options`` property on each modified molecule, and raise an error when modifying a system that has already been modified, unless ``--force`` is used. Forced modifications are applied on top of the existing ones [#20](https://github.com/OpenBioSim/ghostly/pull/20). * Added support for Python 3.13 and 3.14, and dropped support for Python 3.10 and 3.11. -* Fixed bond order perception for end states with ghost atoms attached, which could fail and leave every heavy atom typed as sp3, and decide triple junction planarity and ring membership from the end state angle terms and connectivity [#23](https://github.com/OpenBioSim/ghostly/pull/23). +* Fixed RDKit conversion of end states with ghost atoms attached. The extra bonds to ghosts could exceed the valence of the bridge atoms, so bond orders couldn't be perceived and every heavy atom was silently typed as sp3, with no aromaticity. Ghost atoms are now removed before conversion, and a failed perception raises an error rather than returning a molecule without bond orders [#23](https://github.com/OpenBioSim/ghostly/pull/23). +* Decide triple junction planarity from the energy needed to flatten the angle terms between the physical neighbours of the bridge, rather than its RDKit hybridisation, and ring membership from the end state connectivity, including bonds to ghosts [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Keep improper dihedrals that hold a single ghost in the plane of a ring bridge, since nothing else stops the ghost crossing the plane when the angles aren't stiffened, and remove the ``--stiffen-ring-bridges`` option [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Always stiffen ghost angles at bridges that aren't in a ring, and remove the ``--stiffen-sp2-bridges`` option. Without impropers, the angle terms of an sp2 bridge don't hold a ghost in the plane, so skipping stiffening let it cross freely [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Detect rotatable bridge bonds from the angle terms of the bridge, rather than its RDKit hybridisation, so that rotamer detection is consistent with junction planarity [#23](https://github.com/OpenBioSim/ghostly/pull/23). From 2cc4cedc5b35c1c012baa7b99d06725c27a37ce0 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 10 Nov 2026 16:52:46 +0000 Subject: [PATCH 08/20] Find macrocycle rings and ignore rings closed by ghosts. --- src/ghostly/_ghostly.py | 100 +++++++++++++++++++++++++++++----------- 1 file changed, 74 insertions(+), 26 deletions(-) diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index 657a64d..f5edb9b 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -48,6 +48,9 @@ # Triple junctions needing more than this energy (kcal/mol) to flatten are non-planar. _MAX_PLANAR_FLATTENING = 2.0 +# The longest path searched when checking for rings, which covers macrocycles. +_MAX_RING_PATH = 50 + def modify( system, @@ -440,8 +443,10 @@ def modify( num_optimise=num_optimise, ) - # Remove any improper dihedrals connecting ghosts to the physical region. - mol = _remove_impropers(mol, ghosts0, modifications, is_lambda1=False) + # Remove any improper dihedrals connecting ghosts to the physical region. + mol = _remove_impropers( + mol, ghosts0, connectivity0, modifications, is_lambda1=False + ) # Remove any residual ghost dihedrals not caught by the per-bridge # junction handlers (cross-bridge and ghost-middle patterns). @@ -475,6 +480,8 @@ def modify( physical0, ghosts0, ghosts1, + connectivity0, + connectivity1, modifications, is_lambda1=False, soften_anchors=soften_anchors, @@ -545,8 +552,10 @@ def modify( is_lambda1=True, ) - # Remove any improper dihedrals connecting ghosts to the physical region. - mol = _remove_impropers(mol, ghosts1, modifications, is_lambda1=True) + # Remove any improper dihedrals connecting ghosts to the physical region. + mol = _remove_impropers( + mol, ghosts1, connectivity1, modifications, is_lambda1=True + ) # Remove any residual ghost dihedrals not caught by the per-bridge # junction handlers (cross-bridge and ghost-middle patterns). @@ -580,6 +589,8 @@ def modify( physical1, ghosts1, ghosts0, + connectivity1, + connectivity0, modifications, is_lambda1=True, soften_anchors=soften_anchors, @@ -928,7 +939,9 @@ def _dual( mod_key = "lambda_0" suffix = "0" - bridge_in_ring = connectivity.in_ring(bridge) + bridge_in_ring = _in_ring( + _physical_connectivity(mol, connectivity, is_lambda1), bridge + ) # Score the physical neighbours. if bridge_indices is not None: @@ -1322,7 +1335,9 @@ def _triple( _logger.debug(f" Flattening energy: {flattening:.3f} kcal/mol") is_planar = flattening < _MAX_PLANAR_FLATTENING - bridge_in_ring = connectivity.in_ring(bridge) + bridge_in_ring = _in_ring( + _physical_connectivity(mol, connectivity, is_lambda1), bridge + ) # Score the physical neighbours. if bridge_indices is not None: @@ -1887,7 +1902,7 @@ def _higher( ) -def _remove_impropers(mol, ghosts, modifications, is_lambda1=False): +def _remove_impropers(mol, ghosts, connectivity, modifications, is_lambda1=False): """ Remove improper dihedral terms that bridge the ghost and physical systems, other than those centred on a physical ring atom with a single ghost. These @@ -1902,6 +1917,9 @@ def _remove_impropers(mol, ghosts, modifications, is_lambda1=False): ghosts : List[sire.legacy.Mol.AtomIdx] The list of ghost atoms. + connectivity : sire.legacy.Mol.Connectivity + The connectivity of the molecule at the end state. + modifications : dict A dictionary to store details of the modifications made. @@ -1929,7 +1947,7 @@ def _remove_impropers(mol, ghosts, modifications, is_lambda1=False): # Initialise a container to store the updated bonded terms. new_impropers = _SireMM.FourAtomFunctions(mol.info()) - connectivity = _end_state_connectivity(mol, is_lambda1) + physical_connectivity = _physical_connectivity(mol, connectivity, is_lambda1) # Loop over the improper dihedrals. for p in impropers.potentials(): @@ -1953,7 +1971,7 @@ def _remove_impropers(mol, ghosts, modifications, is_lambda1=False): num_ghosts == 1 and centre is not None and centre not in ghosts - and connectivity.in_ring(centre) + and _in_ring(physical_connectivity, centre) ) # Remove any improper dihedrals that bridge the ghost and physical systems. @@ -2447,6 +2465,8 @@ def _modify_mixed_dihedrals( physical, ghosts, ref_ghosts, + connectivity, + ref_connectivity, modifications, is_lambda1=False, soften_anchors=1.0, @@ -2544,7 +2564,13 @@ def _modify_mixed_dihedrals( ref = _physical_state(mol, is_lambda1) rotatable_bonds = _find_rotatable_bonds( - ref, bridges, physical, is_lambda1=is_lambda1, warn=not stiffen_rotamers + ref, + bridges, + physical, + connectivity, + ref_connectivity, + is_lambda1=is_lambda1, + warn=not stiffen_rotamers, ) # Nothing to do if no modifications are requested. @@ -2655,7 +2681,9 @@ def _modify_mixed_dihedrals( return mol -def _find_rotatable_bonds(ref, bridges, physical, is_lambda1=False, warn=True): +def _find_rotatable_bonds( + ref, bridges, physical, connectivity, ref_connectivity, is_lambda1=False, warn=True +): r""" Find the rotatable bridge--physical bonds, i.e. those that are not in a ring and whose bridge atom is neither planar nor linear in the end state @@ -2683,6 +2711,12 @@ def _find_rotatable_bonds(ref, bridges, physical, is_lambda1=False, warn=True): physical : dict A dictionary mapping bridge atoms to their physical neighbours. + connectivity : sire.legacy.Mol.Connectivity + The connectivity of the molecule at the end state of the bridges. + + ref_connectivity : sire.legacy.Mol.Connectivity + The connectivity of the molecule at the opposite end state. + is_lambda1 : bool, optional Whether the bridges are at lambda = 1. @@ -2702,18 +2736,16 @@ def _find_rotatable_bonds(ref, bridges, physical, is_lambda1=False, warn=True): if not bridges: return rotatable_bonds - connectivity = _end_state_connectivity(ref, is_lambda1) - for bridge in bridges: # The ghosts are physical at the opposite end state. - if not _is_rotor(ref, bridge, not is_lambda1): + if not _is_rotor(ref, bridge, ref_connectivity, not is_lambda1): continue for p in physical[bridge]: # Ring membership from the connectivity, since the physical # neighbour may be a ghost at the opposite end state. - if not connectivity.are_bonded(p, bridge) or connectivity.in_ring( - p, bridge + if not connectivity.are_bonded(p, bridge) or _in_ring( + connectivity, p, bridge ): continue @@ -2765,16 +2797,32 @@ def _to_rdkit_end_state(mol, is_lambda1=False): return rdmol, {index: i for i, index in enumerate(physical)} -def _end_state_connectivity(mol, is_lambda1=False): +def _physical_connectivity(mol, connectivity, is_lambda1=False): """ - Return the connectivity of a perturbable molecule at an end state, including - the bonds to ghost atoms. Only ring-breaking molecules have end state - connectivities that differ. + Return the connectivity without the bonds to the ghost atoms of an end state. """ - key = "connectivity1" if is_lambda1 else "connectivity0" - if mol.has_property(key): - return mol.property(key) - return mol.property("connectivity") + info = mol.info() + elements = mol.property("element1" if is_lambda1 else "element0") + physical = _SireMol.Connectivity(info).edit() + for bond in connectivity.get_bonds(): + idx0 = info.atom_idx(bond.atom0()) + idx1 = info.atom_idx(bond.atom1()) + if all(elements[idx.value()].num_protons() > 0 for idx in (idx0, idx1)): + physical.connect(idx0, idx1) + return physical.commit() + + +def _in_ring(connectivity, atom0, atom1=None): + """ + Return whether an atom, or the bond between two atoms, is in a ring. Sire's + in_ring only finds rings of up to 12 atoms, which misses macrocycles. + """ + if atom1 is not None: + return len(connectivity.find_paths(atom0, atom1, _MAX_RING_PATH)) > 1 + return any( + _in_ring(connectivity, atom0, other) + for other in connectivity.connections_to(atom0) + ) def _angle_terms(mol, bridge, atoms, is_lambda1=False): @@ -2819,7 +2867,7 @@ def _flattening_energy(mol, bridge, atoms, is_lambda1=False): return shortfall**2 / sum(1 / term.k() for term in terms.values()) -def _is_rotor(mol, bridge, is_lambda1=False, max_angle=170): +def _is_rotor(mol, bridge, connectivity, is_lambda1=False, max_angle=170): """ Return whether rotation about the bonds of the bridge atom can move its other substituents at an end state, i.e. the bridge is neither planar, by @@ -2830,7 +2878,7 @@ def _is_rotor(mol, bridge, is_lambda1=False, max_angle=170): elements = mol.property("element1" if is_lambda1 else "element0") neighbours = [ n - for n in _end_state_connectivity(mol, is_lambda1).connections_to(bridge) + for n in connectivity.connections_to(bridge) if elements[n.value()].num_protons() > 0 ] From 654ee31c6f82b83dd01ec49109fd3b312f5055e4 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 10 Nov 2026 16:55:51 +0000 Subject: [PATCH 09/20] Use safer fallbacks for missing terms and chirality. --- src/ghostly/_ghostly.py | 27 ++++++++++++++------------- tests/test_ghostly.py | 11 ++++------- 2 files changed, 18 insertions(+), 20 deletions(-) diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index f5edb9b..4eef767 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -1027,9 +1027,7 @@ def _is_heavy(atom_idx): # Identify the physical atom in this angle. phys_atom = idx2 if idx0 in ghosts else idx0 - # Ring bridges keep their ghost impropers, which hold the ghost - # in the plane. The rigid ring keeps the redundant coupling - # acceptable (Boresch et al.), so the angles are not stiffened. + # Ring bridge impropers hold the ghost in the plane, so don't stiffen. if bridge_in_ring: new_angles.set(idx0, idx1, idx2, p.function()) _logger.debug( @@ -1317,9 +1315,7 @@ def _triple( mod_key = "lambda_0" suffix = "0" - # The junction is planar if the angle terms between the physical neighbours - # of the bridge atom allow them to flatten at thermal energies. Higher order - # junctions pass four or more physical neighbours, so are never planar. + # A bridge with four or more physical neighbours is never planar. if len(physical) > 3: is_planar = False else: @@ -1328,9 +1324,9 @@ def _triple( _logger.warning( f"Missing angle terms between the physical neighbours of bridge " f"atom {bridge.value()} at {_lam_sym} = {int(is_lambda1)}. " - f"Defaulting to planar junction." + f"Defaulting to non-planar junction." ) - is_planar = True + is_planar = False else: _logger.debug(f" Flattening energy: {flattening:.3f} kcal/mol") is_planar = flattening < _MAX_PLANAR_FLATTENING @@ -1634,9 +1630,7 @@ def _is_heavy(atom_idx): idxs = [idx0, idx1, idx2] if idx in idxs and any([x in ghosts for x in idxs]): - # Ring bridges aren't stiffened, so preserve all ghost-containing - # angles at their original values rather than removing the - # sacrificial one. + # Ring bridges aren't stiffened, so keep the sacrificial angle. if bridge_in_ring: new_angles.set(idx0, idx1, idx2, p.function()) _logger.debug( @@ -2214,9 +2208,16 @@ def _linearise_ring_break( # chirality directly rather than inferring it from substituent counts. from rdkit.Chem import ChiralType as _ChiralType - rdmol, index = _to_rdkit_end_state(mol, not is_lambda1) - linearised = set() + + try: + rdmol, index = _to_rdkit_end_state(mol, not is_lambda1) + except Exception as e: + _logger.warning( + f"Skipping ring-break linearisation at {_lam_sym} = {int(is_lambda1)}, " + f"since chirality can't be checked: {e}" + ) + return mol, linearised for ghost in ghosts: all_neighbors = list(connectivity.connections_to(ghost)) physical_neighbors = [n for n in all_neighbors if n not in ghost_set] diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index 78cf8ce..4887267 100644 --- a/tests/test_ghostly.py +++ b/tests/test_ghostly.py @@ -124,10 +124,8 @@ def test_toluene_to_pyridine(): for dihedral in missing_dihedrals ) - # Bridge atom 1 is in an aromatic ring, so no stiffening is applied. Its - # improper holds the ghost in the plane of the ring, and the rigid ring - # keeps the coupling acceptable (Boresch et al.). The angles involving the - # ring-bridge ghost should be unchanged. + # Bridge atom 1 is in an aromatic ring, so the angles involving its ghost + # are unchanged and its improper is kept. ring_bridge_angles = [ (AtomIdx(0), AtomIdx(1), AtomIdx(2)), (AtomIdx(0), AtomIdx(1), AtomIdx(6)), @@ -401,9 +399,8 @@ def test_ejm49_to_ejm31(): assert all(check_dihedral(info, new_dihedrals1.potentials(), *d) for d in anchor1) - # Bridge atom 17 is sp2 and not in a ring. Its angle terms don't hold the - # ghost in the plane once the impropers are removed, so the - # physical-bridge-ghost angles are stiffened to 90 degrees at both end states. + # Bridge atom 17 is sp2 and not in a ring, so its physical-bridge-ghost + # angles are stiffened to 90 degrees at both end states. sp2_bridge_angles0 = [ (AtomIdx(16), AtomIdx(17), AtomIdx(39)), (AtomIdx(18), AtomIdx(17), AtomIdx(39)), From 0b82bc25ba427d4b63e8a3772a652c2f98b1c7a9 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 10 Nov 2026 16:57:48 +0000 Subject: [PATCH 10/20] Add tests for ring, planarity and RDKit helpers. --- tests/test_ghostly.py | 140 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 140 insertions(+) diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index 4887267..36a9fda 100644 --- a/tests/test_ghostly.py +++ b/tests/test_ghostly.py @@ -934,6 +934,146 @@ def test_single_well_skips_reference_ghost(): assert _single_well(ref, idxs, 5, {AtomIdx(1)}) is None +def test_flattening_energy(): + """ + Test the energy needed to flatten a bridge atom and three of its neighbours. + """ + + import sire.morph as morph + from sire.legacy.Mol import AtomIdx + + from ghostly._ghostly import ( + _MAX_PLANAR_FLATTENING, + _create_connectivity, + _flattening_energy, + ) + + # An aromatic carbon in toluene is planar. + mol = sr.load_test_files("tol2pyr.s3")[0] + atoms = [AtomIdx(0), AtomIdx(2), AtomIdx(6)] + assert _flattening_energy(mol, AtomIdx(1), atoms) < 0.01 + + # An sp3 carbon in hexane is not. + mol = sr.load_test_files("hex2prp.s3")[0] + connectivity = _create_connectivity(morph.link_to_reference(mol)) + neighbours = list(connectivity.connections_to(AtomIdx(1))) + assert _flattening_energy(mol, AtomIdx(1), neighbours[:3]) > _MAX_PLANAR_FLATTENING + + # There is no energy without three angle terms. + assert _flattening_energy(mol, AtomIdx(1), neighbours[:2]) is None + + +def test_is_rotor(): + """ + Test whether bridge atoms are rotors, i.e. neither planar nor linear. + """ + + import sire.morph as morph + from sire.legacy.Mol import AtomIdx + + from ghostly._ghostly import _create_connectivity, _is_rotor + + # Four neighbours: an sp3 carbon in hexane. + mol = sr.load_test_files("hex2prp.s3")[0] + connectivity = _create_connectivity(morph.link_to_reference(mol)) + assert _is_rotor(mol, AtomIdx(1), connectivity) + + # Three planar neighbours: the carbonyl carbon in acetone. + mol = sr.load_test_files("acepol.s3")[0] + connectivity0 = _create_connectivity(morph.link_to_reference(mol)) + assert not _is_rotor(mol, AtomIdx(1), connectivity0) + + # Two neighbours: the hydroxyl oxygen in propenol, which is bent at about + # 108 degrees, so only counts as linear with a lower threshold. + connectivity1 = _create_connectivity(morph.link_to_perturbed(mol)) + assert _is_rotor(mol, AtomIdx(9), connectivity1, is_lambda1=True) + assert not _is_rotor(mol, AtomIdx(9), connectivity1, is_lambda1=True, max_angle=100) + + +def test_in_ring(): + """ + Test ring detection, including macrocycles that Sire's in_ring misses. + """ + + from sire.legacy.Mol import AtomIdx + + from ghostly._ghostly import _in_ring + + connectivity = sr.smiles("C1CCCCCCCCCCCCC1").property("connectivity") + assert not connectivity.in_ring(AtomIdx(0)) + assert _in_ring(connectivity, AtomIdx(0)) + assert _in_ring(connectivity, AtomIdx(0), AtomIdx(1)) + + mol = sr.smiles("C1CCCCC1CC") + connectivity = mol.property("connectivity") + assert _in_ring(connectivity, AtomIdx(0)) + assert not _in_ring(connectivity, AtomIdx(mol.num_atoms() - 1)) + assert not _in_ring(connectivity, AtomIdx(5), AtomIdx(6)) + + +def test_physical_connectivity_ring(): + """ + Test that a ring bridge is found from the physical atoms of an end state. + """ + + import sire.morph as morph + from sire.legacy.Mol import AtomIdx + + from ghostly._ghostly import _create_connectivity, _in_ring, _physical_connectivity + + # Toluene to pyridine: at lambda = 1 the methyl carbon is a ghost bonded to + # the ring, which stays a ring once the bonds to ghosts are removed. + mol = sr.load_test_files("tol2pyr.s3")[0] + connectivity = _create_connectivity(morph.link_to_perturbed(mol)) + physical = _physical_connectivity(mol, connectivity, is_lambda1=True) + assert not physical.are_bonded(AtomIdx(0), AtomIdx(1)) + assert _in_ring(physical, AtomIdx(1)) + + +def test_to_rdkit_end_state(): + """ + Test the conversion of an end state to RDKit without its ghost atoms. + """ + + from ghostly._ghostly import _to_rdkit_end_state + + mol = sr.load_test_files("tol2pyr.s3")[0] + + # At lambda = 1 the methyl group is a ghost, so atom 1, the pyridine + # nitrogen, comes first. + rdmol, index = _to_rdkit_end_state(mol, is_lambda1=True) + num_physical = sum(e.num_protons() > 0 for e in mol.property("element1")) + assert rdmol.GetNumAtoms() == num_physical + assert 0 not in index + assert index[1] == 0 + atom = rdmol.GetAtomWithIdx(index[1]) + assert atom.GetSymbol() == "N" and atom.GetIsAromatic() + + # At lambda = 0 every atom is physical and the ring is aromatic. + rdmol, index = _to_rdkit_end_state(mol, is_lambda1=False) + assert rdmol.GetNumAtoms() == mol.num_atoms() + assert rdmol.GetAtomWithIdx(index[1]).GetIsAromatic() + + +def test_linearise_ring_break_without_chirality(monkeypatch): + """ + Test that a failed chirality check skips linearisation rather than raising. + """ + + import ghostly._ghostly as _ghostly + + def fail(*args, **kwargs): + raise RuntimeError("no bond orders") + + monkeypatch.setattr(_ghostly, "_to_rdkit_end_state", fail) + + mols = sr.load_test_files("hex2prp.s3") + _, modifications = modify(mols, linearise_ring_break=True) + + for state in ("lambda_0", "lambda_1"): + assert modifications[state]["linearised_ring_break"] == [] + + def dihedral_functions(mol, prop): """ Return a dictionary mapping dihedral indices to their functions. From 6caf58163a96769b85d777575d4f050f9d422f12 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 10 Nov 2026 16:58:30 +0000 Subject: [PATCH 11/20] Update CHANGELOG. --- CHANGELOG.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 06d35fe..2aa4cf7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,8 +10,9 @@ Changelog * Added ``--single-well-anchors`` and ``--k-single-well`` to replace surviving anchor dihedrals with a single-well cosine potential whose minimum is the geometry of the end state where the ghost is physical. Multi-well anchors, e.g. the n=2 term of a dummy hydrogen on a carbonyl oxygen, can trap ghosts in an alternative well at intermediate lambda, giving hysteresis between independently equilibrated windows. Anchors with zero force constants or an ill-defined reference dihedral are left unchanged. (Default: False, 5 kcal/mol) [#19](https://github.com/OpenBioSim/ghostly/pull/19). * Record the options used for each modification in a ``ghostly_options`` property on each modified molecule, and raise an error when modifying a system that has already been modified, unless ``--force`` is used. Forced modifications are applied on top of the existing ones [#20](https://github.com/OpenBioSim/ghostly/pull/20). * Added support for Python 3.13 and 3.14, and dropped support for Python 3.10 and 3.11. -* Fixed RDKit conversion of end states with ghost atoms attached. The extra bonds to ghosts could exceed the valence of the bridge atoms, so bond orders couldn't be perceived and every heavy atom was silently typed as sp3, with no aromaticity. Ghost atoms are now removed before conversion, and a failed perception raises an error rather than returning a molecule without bond orders [#23](https://github.com/OpenBioSim/ghostly/pull/23). -* Decide triple junction planarity from the energy needed to flatten the angle terms between the physical neighbours of the bridge, rather than its RDKit hybridisation, and ring membership from the end state connectivity, including bonds to ghosts [#23](https://github.com/OpenBioSim/ghostly/pull/23). +* Fixed bridge atoms being typed by RDKit on end states with ghost atoms attached, where bond orders often couldn't be inferred, leaving every heavy atom typed as sp3 with no aromaticity. Bridge typing no longer uses RDKit. The chirality check used by ``--linearise-ring-break`` now converts the ghost-free end state, and skips linearisation with a warning if bond orders can't be inferred [#23](https://github.com/OpenBioSim/ghostly/pull/23). +* Decide triple junction planarity from the energy needed to flatten the angle terms between the physical neighbours of the bridge, rather than its RDKit hybridisation. Junctions with a missing angle term are treated as non-planar [#23](https://github.com/OpenBioSim/ghostly/pull/23). +* Decide whether a bridge is in a ring from the physical atoms of the end state, so rings closed only by ghost atoms don't count, and detect rings of up to 50 atoms, since Sire's ``in_ring`` misses rings of more than 12 [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Keep improper dihedrals that hold a single ghost in the plane of a ring bridge, since nothing else stops the ghost crossing the plane when the angles aren't stiffened, and remove the ``--stiffen-ring-bridges`` option [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Always stiffen ghost angles at bridges that aren't in a ring, and remove the ``--stiffen-sp2-bridges`` option. Without impropers, the angle terms of an sp2 bridge don't hold a ghost in the plane, so skipping stiffening let it cross freely [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Detect rotatable bridge bonds from the angle terms of the bridge, rather than its RDKit hybridisation, so that rotamer detection is consistent with junction planarity [#23](https://github.com/OpenBioSim/ghostly/pull/23). From 3ae12856c9ad0bc0f3adfb382e3b1e4e3d7ef7ee Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 10 Nov 2026 20:08:49 +0000 Subject: [PATCH 12/20] Simplify ring checks and connectivity handling. --- src/ghostly/_ghostly.py | 105 +++++++++++++++++++++++++++------------- tests/test_ghostly.py | 7 ++- 2 files changed, 77 insertions(+), 35 deletions(-) diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index 4eef767..3f5aae1 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -48,9 +48,6 @@ # Triple junctions needing more than this energy (kcal/mol) to flatten are non-planar. _MAX_PLANAR_FLATTENING = 2.0 -# The longest path searched when checking for rings, which covers macrocycles. -_MAX_RING_PATH = 50 - def modify( system, @@ -382,6 +379,20 @@ def modify( bridge_indices0 = set(b.value() for b in bridges0) bridge_indices1 = set(b.value() for b in bridges1) + # The end state connectivities without the bonds to ghosts, for rings. + physical_connectivity0 = _physical_connectivity(mol, connectivity0, ghosts0) + physical_connectivity1 = _physical_connectivity(mol, connectivity1, ghosts1) + + # Remove any improper dihedrals connecting ghosts to the physical region. + mol = _remove_impropers( + mol, + ghosts0, + connectivity0, + physical_connectivity0, + modifications, + is_lambda1=False, + ) + # First lambda = 0. for b in bridges0: # Determine the type of junction. @@ -410,6 +421,7 @@ def modify( modifications, k_hard=k_hard, bridge_indices=bridge_indices0, + physical_connectivity=physical_connectivity0, ) # Triple junction. @@ -426,6 +438,7 @@ def modify( optimise_angles=optimise_angles, num_optimise=num_optimise, bridge_indices=bridge_indices0, + physical_connectivity=physical_connectivity0, ) # Higher order junction. @@ -441,13 +454,9 @@ def modify( k_soft=k_soft, optimise_angles=optimise_angles, num_optimise=num_optimise, + physical_connectivity=physical_connectivity0, ) - # Remove any improper dihedrals connecting ghosts to the physical region. - mol = _remove_impropers( - mol, ghosts0, connectivity0, modifications, is_lambda1=False - ) - # Remove any residual ghost dihedrals not caught by the per-bridge # junction handlers (cross-bridge and ghost-middle patterns). mol = _remove_residual_ghost_dihedrals( @@ -492,6 +501,15 @@ def modify( ) # Now lambda = 1. + mol = _remove_impropers( + mol, + ghosts1, + connectivity1, + physical_connectivity1, + modifications, + is_lambda1=True, + ) + for b in bridges1: junction = len(physical1[b]) @@ -518,6 +536,7 @@ def modify( k_hard=k_hard, is_lambda1=True, bridge_indices=bridge_indices1, + physical_connectivity=physical_connectivity1, ) elif junction == 3: @@ -534,6 +553,7 @@ def modify( num_optimise=num_optimise, is_lambda1=True, bridge_indices=bridge_indices1, + physical_connectivity=physical_connectivity1, ) # Higher order junction. @@ -550,13 +570,9 @@ def modify( optimise_angles=optimise_angles, num_optimise=num_optimise, is_lambda1=True, + physical_connectivity=physical_connectivity1, ) - # Remove any improper dihedrals connecting ghosts to the physical region. - mol = _remove_impropers( - mol, ghosts1, connectivity1, modifications, is_lambda1=True - ) - # Remove any residual ghost dihedrals not caught by the per-bridge # junction handlers (cross-bridge and ghost-middle patterns). mol = _remove_residual_ghost_dihedrals( @@ -865,6 +881,7 @@ def _dual( k_hard=100, is_lambda1=False, bridge_indices=None, + physical_connectivity=None, ): r""" Apply modifications to a dual junction. @@ -939,9 +956,7 @@ def _dual( mod_key = "lambda_0" suffix = "0" - bridge_in_ring = _in_ring( - _physical_connectivity(mol, connectivity, is_lambda1), bridge - ) + bridge_in_ring = _in_ring(physical_connectivity, bridge) # Score the physical neighbours. if bridge_indices is not None: @@ -1209,6 +1224,7 @@ def _is_heavy(atom_idx): k_hard=k_hard, is_lambda1=is_lambda1, bridge_indices=bridge_indices, + physical_connectivity=physical_connectivity, ) # Return the updated molecule. @@ -1228,6 +1244,7 @@ def _triple( num_optimise=10, is_lambda1=False, bridge_indices=None, + physical_connectivity=None, ): r""" Apply modifications to a triple junction. @@ -1331,9 +1348,7 @@ def _triple( _logger.debug(f" Flattening energy: {flattening:.3f} kcal/mol") is_planar = flattening < _MAX_PLANAR_FLATTENING - bridge_in_ring = _in_ring( - _physical_connectivity(mol, connectivity, is_lambda1), bridge - ) + bridge_in_ring = _in_ring(physical_connectivity, bridge) # Score the physical neighbours. if bridge_indices is not None: @@ -1690,6 +1705,7 @@ def _is_heavy(atom_idx): k_hard=k_hard, is_lambda1=is_lambda1, bridge_indices=bridge_indices, + physical_connectivity=physical_connectivity, ) # Return the updated molecule. @@ -1708,6 +1724,7 @@ def _higher( optimise_angles=False, num_optimise=10, is_lambda1=False, + physical_connectivity=None, ): r""" Apply modifications to higher order junctions. @@ -1771,6 +1788,9 @@ def _higher( if not ghosts: return mol + # Work on a copy, since the caller's list is also used for rotamer detection. + physical = list(physical) + # Store the molecular info. info = mol.info() @@ -1893,10 +1913,13 @@ def _higher( optimise_angles=optimise_angles, num_optimise=num_optimise, is_lambda1=is_lambda1, + physical_connectivity=physical_connectivity, ) -def _remove_impropers(mol, ghosts, connectivity, modifications, is_lambda1=False): +def _remove_impropers( + mol, ghosts, connectivity, physical_connectivity, modifications, is_lambda1=False +): """ Remove improper dihedral terms that bridge the ghost and physical systems, other than those centred on a physical ring atom with a single ghost. These @@ -1914,6 +1937,9 @@ def _remove_impropers(mol, ghosts, connectivity, modifications, is_lambda1=False connectivity : sire.legacy.Mol.Connectivity The connectivity of the molecule at the end state. + physical_connectivity : sire.legacy.Mol.Connectivity + The connectivity at the end state without the bonds to ghost atoms. + modifications : dict A dictionary to store details of the modifications made. @@ -1941,8 +1967,6 @@ def _remove_impropers(mol, ghosts, connectivity, modifications, is_lambda1=False # Initialise a container to store the updated bonded terms. new_impropers = _SireMM.FourAtomFunctions(mol.info()) - physical_connectivity = _physical_connectivity(mol, connectivity, is_lambda1) - # Loop over the improper dihedrals. for p in impropers.potentials(): idx0 = info.atom_idx(p.atom0()) @@ -2798,32 +2822,45 @@ def _to_rdkit_end_state(mol, is_lambda1=False): return rdmol, {index: i for i, index in enumerate(physical)} -def _physical_connectivity(mol, connectivity, is_lambda1=False): +def _physical_connectivity(mol, connectivity, ghosts): """ - Return the connectivity without the bonds to the ghost atoms of an end state. + Return the connectivity without the bonds to the given ghost atoms. """ info = mol.info() - elements = mol.property("element1" if is_lambda1 else "element0") + ghosts = set(ghosts) physical = _SireMol.Connectivity(info).edit() for bond in connectivity.get_bonds(): idx0 = info.atom_idx(bond.atom0()) idx1 = info.atom_idx(bond.atom1()) - if all(elements[idx.value()].num_protons() > 0 for idx in (idx0, idx1)): + if idx0 not in ghosts and idx1 not in ghosts: physical.connect(idx0, idx1) return physical.commit() def _in_ring(connectivity, atom0, atom1=None): """ - Return whether an atom, or the bond between two atoms, is in a ring. Sire's - in_ring only finds rings of up to 12 atoms, which misses macrocycles. + Return whether an atom, or the bond between two bonded atoms, is in a ring. + A bond is in a ring if its atoms stay connected once it is removed. """ - if atom1 is not None: - return len(connectivity.find_paths(atom0, atom1, _MAX_RING_PATH)) > 1 - return any( - _in_ring(connectivity, atom0, other) - for other in connectivity.connections_to(atom0) - ) + if atom1 is None: + return any( + _in_ring(connectivity, atom0, other) + for other in connectivity.connections_to(atom0) + ) + + seen = {atom0} + queue = [atom0] + while queue: + atom = queue.pop() + for other in connectivity.connections_to(atom): + if atom == atom0 and other == atom1: + continue + if other == atom1: + return True + if other not in seen: + seen.add(other) + queue.append(other) + return False def _angle_terms(mol, bridge, atoms, is_lambda1=False): diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index 36a9fda..991981a 100644 --- a/tests/test_ghostly.py +++ b/tests/test_ghostly.py @@ -1025,7 +1025,12 @@ def test_physical_connectivity_ring(): # the ring, which stays a ring once the bonds to ghosts are removed. mol = sr.load_test_files("tol2pyr.s3")[0] connectivity = _create_connectivity(morph.link_to_perturbed(mol)) - physical = _physical_connectivity(mol, connectivity, is_lambda1=True) + ghosts = [ + AtomIdx(i) + for i, e in enumerate(mol.property("element1")) + if e.num_protons() == 0 + ] + physical = _physical_connectivity(mol, connectivity, ghosts) assert not physical.are_bonded(AtomIdx(0), AtomIdx(1)) assert _in_ring(physical, AtomIdx(1)) From 76cc0bed5faa19f1f5316b52353ac4c411d83879 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 10 Nov 2026 20:10:26 +0000 Subject: [PATCH 13/20] Refine missing angle term and chirality fallbacks. --- src/ghostly/_ghostly.py | 63 +++++++++++++++++++++++------------------ 1 file changed, 36 insertions(+), 27 deletions(-) diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index 3f5aae1..0388dba 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -1332,24 +1332,25 @@ def _triple( mod_key = "lambda_0" suffix = "0" + bridge_in_ring = _in_ring(physical_connectivity, bridge) + # A bridge with four or more physical neighbours is never planar. if len(physical) > 3: is_planar = False else: flattening = _flattening_energy(mol, bridge, physical, is_lambda1) if flattening is None: + # Ring bridges keep their angles, as for any other ring bridge. + is_planar = bridge_in_ring _logger.warning( f"Missing angle terms between the physical neighbours of bridge " f"atom {bridge.value()} at {_lam_sym} = {int(is_lambda1)}. " - f"Defaulting to non-planar junction." + f"Defaulting to {'planar' if is_planar else 'non-planar'} junction." ) - is_planar = False else: _logger.debug(f" Flattening energy: {flattening:.3f} kcal/mol") is_planar = flattening < _MAX_PLANAR_FLATTENING - bridge_in_ring = _in_ring(physical_connectivity, bridge) - # Score the physical neighbours. if bridge_indices is not None: phys_scores = {p: _score_atom(mol, p, bridge_indices) for p in physical} @@ -2234,36 +2235,44 @@ def _linearise_ring_break( linearised = set() + # Find ghost atoms with exactly 2 connections, both physical. + candidates = [ + ghost + for ghost in ghosts + if len([n for n in connectivity.connections_to(ghost) if n not in ghost_set]) + == 2 + ] + if not candidates: + return mol, linearised + try: rdmol, index = _to_rdkit_end_state(mol, not is_lambda1) except Exception as e: _logger.warning( - f"Skipping ring-break linearisation at {_lam_sym} = {int(is_lambda1)}, " - f"since chirality can't be checked: {e}" + f"Unable to check chirality at {_lam_sym} = {int(is_lambda1)}, so " + f"ring-break ghosts that could be stereocentres aren't linearised: {e}" ) - return mol, linearised - for ghost in ghosts: - all_neighbors = list(connectivity.connections_to(ghost)) - physical_neighbors = [n for n in all_neighbors if n not in ghost_set] - - # Find ghost atoms with exactly 2 connections, both physical. - if len(physical_neighbors) == 2: - # Check whether this atom is a chiral centre in the physical end - # state using RDKit. If so, linearising to 180° would destroy the - # geometry, so skip it and warn. - rdatom = rdmol.GetAtomWithIdx(index[ghost.value()]) - chiral = rdatom.GetChiralTag() - if chiral in ( + rdmol = None + physical_state = _create_connectivity(_physical_state(mol, is_lambda1)) + + for ghost in candidates: + # Linearising a chiral centre in the physical end state to 180° would + # destroy its geometry. Without RDKit, skip any possible stereocentre. + if rdmol is None: + chiral = len(physical_state.connections_to(ghost)) >= 4 + else: + chiral = rdmol.GetAtomWithIdx(index[ghost.value()]).GetChiralTag() in ( _ChiralType.CHI_TETRAHEDRAL_CW, _ChiralType.CHI_TETRAHEDRAL_CCW, - ): - _logger.warning( - f" Ring-break ghost atom {ghost.value()} is a chiral " - f"centre in the physical end state. Linearisation skipped " - f"to preserve chirality." - ) - else: - linearised.add(ghost) + ) + if chiral: + _logger.warning( + f" Ring-break ghost atom {ghost.value()} is a possible chiral " + f"centre in the physical end state. Linearisation skipped " + f"to preserve chirality." + ) + else: + linearised.add(ghost) if not linearised: return mol, linearised From 3756ee3d042eedc677e422eea26fc96fb97b97af Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 10 Nov 2026 20:11:55 +0000 Subject: [PATCH 14/20] Test ghost-closed rings and the chirality fallback. --- CHANGELOG.md | 2 +- tests/test_ghostly.py | 67 +++++++++++++++++++++++++++++++++++++++---- 2 files changed, 63 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2aa4cf7..622676d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,7 +12,7 @@ Changelog * Added support for Python 3.13 and 3.14, and dropped support for Python 3.10 and 3.11. * Fixed bridge atoms being typed by RDKit on end states with ghost atoms attached, where bond orders often couldn't be inferred, leaving every heavy atom typed as sp3 with no aromaticity. Bridge typing no longer uses RDKit. The chirality check used by ``--linearise-ring-break`` now converts the ghost-free end state, and skips linearisation with a warning if bond orders can't be inferred [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Decide triple junction planarity from the energy needed to flatten the angle terms between the physical neighbours of the bridge, rather than its RDKit hybridisation. Junctions with a missing angle term are treated as non-planar [#23](https://github.com/OpenBioSim/ghostly/pull/23). -* Decide whether a bridge is in a ring from the physical atoms of the end state, so rings closed only by ghost atoms don't count, and detect rings of up to 50 atoms, since Sire's ``in_ring`` misses rings of more than 12 [#23](https://github.com/OpenBioSim/ghostly/pull/23). +* Decide whether a bridge is in a ring, for stiffening and keeping impropers, from the physical atoms of the end state, so rings closed only by ghost atoms don't count. Rings of any size are detected, since Sire's ``in_ring`` misses rings of more than 12 atoms [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Keep improper dihedrals that hold a single ghost in the plane of a ring bridge, since nothing else stops the ghost crossing the plane when the angles aren't stiffened, and remove the ``--stiffen-ring-bridges`` option [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Always stiffen ghost angles at bridges that aren't in a ring, and remove the ``--stiffen-sp2-bridges`` option. Without impropers, the angle terms of an sp2 bridge don't hold a ghost in the plane, so skipping stiffening let it cross freely [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Detect rotatable bridge bonds from the angle terms of the bridge, rather than its RDKit hybridisation, so that rotamer detection is consistent with junction planarity [#23](https://github.com/OpenBioSim/ghostly/pull/23). diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index 991981a..6bdf006 100644 --- a/tests/test_ghostly.py +++ b/tests/test_ghostly.py @@ -1060,11 +1060,34 @@ def test_to_rdkit_end_state(): assert rdmol.GetAtomWithIdx(index[1]).GetIsAromatic() +def test_ring_closed_by_ghosts(): + """ + Test that a ring closed only through ghost atoms isn't a ring once the bonds + to the ghosts are removed. + """ + + from sire.legacy.Mol import AtomIdx + + from ghostly._ghostly import _in_ring, _physical_connectivity + + mol = sr.smiles("C1CCCCC1") + connectivity = mol.property("connectivity") + physical = _physical_connectivity(mol, connectivity, [AtomIdx(2), AtomIdx(3)]) + + assert _in_ring(connectivity, AtomIdx(0)) + assert not _in_ring(physical, AtomIdx(0)) + assert not _in_ring(physical, AtomIdx(0), AtomIdx(1)) + + def test_linearise_ring_break_without_chirality(monkeypatch): """ - Test that a failed chirality check skips linearisation rather than raising. + Test that, when chirality can't be checked, only ring-break ghosts that + could be stereocentres are left unlinearised. """ + import sire.morph as morph + from sire.legacy.Mol import AtomIdx + import ghostly._ghostly as _ghostly def fail(*args, **kwargs): @@ -1072,11 +1095,45 @@ def fail(*args, **kwargs): monkeypatch.setattr(_ghostly, "_to_rdkit_end_state", fail) - mols = sr.load_test_files("hex2prp.s3") - _, modifications = modify(mols, linearise_ring_break=True) + def modifications(): + def state(): + return { + "removed_angles": [], + "removed_dihedrals": [], + "stiffened_angles": {}, + "softened_angles": {}, + "softened_dihedrals": [], + "stiffened_dihedrals": [], + "single_well_dihedrals": [], + "linearised_ring_break": [], + "softened_bonds": {}, + } + + return {"lambda_0": state(), "lambda_1": state()} + + # Treating the pyridine nitrogen and the methyl carbon as ghosts at lambda = 1 + # leaves the nitrogen bonded to two physical atoms. It has three neighbours + # in toluene, so can't be a stereocentre, and is linearised. + mol = sr.load_test_files("tol2pyr.s3")[0] + connectivity = _ghostly._create_connectivity(morph.link_to_perturbed(mol)) + _, linearised = _ghostly._linearise_ring_break( + mol, [AtomIdx(1), AtomIdx(0)], connectivity, modifications(), is_lambda1=True + ) + assert {ghost.value() for ghost in linearised} == {1} - for state in ("lambda_0", "lambda_1"): - assert modifications[state]["linearised_ring_break"] == [] + # A hexane carbon with its hydrogens as ghosts has four or more neighbours + # in propane, so could be a stereocentre, and is skipped. + mol = sr.load_test_files("hex2prp.s3")[0] + connectivity = _ghostly._create_connectivity(morph.link_to_reference(mol)) + hydrogens = [ + n + for n in connectivity.connections_to(AtomIdx(2)) + if mol.property("element0")[n.value()].symbol() == "H" + ] + _, linearised = _ghostly._linearise_ring_break( + mol, [AtomIdx(2), *hydrogens], connectivity, modifications() + ) + assert linearised == set() def dihedral_functions(mol, prop): From 7d0aab84f3330421f63559ed71cc5b0172e9e8bc Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 10 Nov 2026 21:44:01 +0000 Subject: [PATCH 15/20] Use Sire in_ring now that it finds rings of any size. --- CHANGELOG.md | 2 +- src/ghostly/_ghostly.py | 36 +++++------------------------------- tests/test_ghostly.py | 27 ++++++++++++--------------- 3 files changed, 18 insertions(+), 47 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 622676d..f4b235eff 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,7 +12,7 @@ Changelog * Added support for Python 3.13 and 3.14, and dropped support for Python 3.10 and 3.11. * Fixed bridge atoms being typed by RDKit on end states with ghost atoms attached, where bond orders often couldn't be inferred, leaving every heavy atom typed as sp3 with no aromaticity. Bridge typing no longer uses RDKit. The chirality check used by ``--linearise-ring-break`` now converts the ghost-free end state, and skips linearisation with a warning if bond orders can't be inferred [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Decide triple junction planarity from the energy needed to flatten the angle terms between the physical neighbours of the bridge, rather than its RDKit hybridisation. Junctions with a missing angle term are treated as non-planar [#23](https://github.com/OpenBioSim/ghostly/pull/23). -* Decide whether a bridge is in a ring, for stiffening and keeping impropers, from the physical atoms of the end state, so rings closed only by ghost atoms don't count. Rings of any size are detected, since Sire's ``in_ring`` misses rings of more than 12 atoms [#23](https://github.com/OpenBioSim/ghostly/pull/23). +* Decide whether a bridge is in a ring, for stiffening and keeping impropers, from the physical atoms of the end state, so rings closed only by ghost atoms don't count [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Keep improper dihedrals that hold a single ghost in the plane of a ring bridge, since nothing else stops the ghost crossing the plane when the angles aren't stiffened, and remove the ``--stiffen-ring-bridges`` option [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Always stiffen ghost angles at bridges that aren't in a ring, and remove the ``--stiffen-sp2-bridges`` option. Without impropers, the angle terms of an sp2 bridge don't hold a ghost in the plane, so skipping stiffening let it cross freely [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Detect rotatable bridge bonds from the angle terms of the bridge, rather than its RDKit hybridisation, so that rotamer detection is consistent with junction planarity [#23](https://github.com/OpenBioSim/ghostly/pull/23). diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index 0388dba..68db07f 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -956,7 +956,7 @@ def _dual( mod_key = "lambda_0" suffix = "0" - bridge_in_ring = _in_ring(physical_connectivity, bridge) + bridge_in_ring = physical_connectivity.in_ring(bridge) # Score the physical neighbours. if bridge_indices is not None: @@ -1332,7 +1332,7 @@ def _triple( mod_key = "lambda_0" suffix = "0" - bridge_in_ring = _in_ring(physical_connectivity, bridge) + bridge_in_ring = physical_connectivity.in_ring(bridge) # A bridge with four or more physical neighbours is never planar. if len(physical) > 3: @@ -1990,7 +1990,7 @@ def _remove_impropers( num_ghosts == 1 and centre is not None and centre not in ghosts - and _in_ring(physical_connectivity, centre) + and physical_connectivity.in_ring(centre) ) # Remove any improper dihedrals that bridge the ghost and physical systems. @@ -2778,8 +2778,8 @@ def _find_rotatable_bonds( for p in physical[bridge]: # Ring membership from the connectivity, since the physical # neighbour may be a ghost at the opposite end state. - if not connectivity.are_bonded(p, bridge) or _in_ring( - connectivity, p, bridge + if not connectivity.are_bonded(p, bridge) or connectivity.in_ring( + p, bridge ): continue @@ -2846,32 +2846,6 @@ def _physical_connectivity(mol, connectivity, ghosts): return physical.commit() -def _in_ring(connectivity, atom0, atom1=None): - """ - Return whether an atom, or the bond between two bonded atoms, is in a ring. - A bond is in a ring if its atoms stay connected once it is removed. - """ - if atom1 is None: - return any( - _in_ring(connectivity, atom0, other) - for other in connectivity.connections_to(atom0) - ) - - seen = {atom0} - queue = [atom0] - while queue: - atom = queue.pop() - for other in connectivity.connections_to(atom): - if atom == atom0 and other == atom1: - continue - if other == atom1: - return True - if other not in seen: - seen.add(other) - queue.append(other) - return False - - def _angle_terms(mol, bridge, atoms, is_lambda1=False): """ Return the angle terms centred on the bridge atom between pairs of the diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index 6bdf006..732c570 100644 --- a/tests/test_ghostly.py +++ b/tests/test_ghostly.py @@ -992,23 +992,20 @@ def test_is_rotor(): def test_in_ring(): """ - Test ring detection, including macrocycles that Sire's in_ring misses. + Test ring detection, including macrocycles. """ from sire.legacy.Mol import AtomIdx - from ghostly._ghostly import _in_ring - connectivity = sr.smiles("C1CCCCCCCCCCCCC1").property("connectivity") - assert not connectivity.in_ring(AtomIdx(0)) - assert _in_ring(connectivity, AtomIdx(0)) - assert _in_ring(connectivity, AtomIdx(0), AtomIdx(1)) + assert connectivity.in_ring(AtomIdx(0)) + assert connectivity.in_ring(AtomIdx(0), AtomIdx(1)) mol = sr.smiles("C1CCCCC1CC") connectivity = mol.property("connectivity") - assert _in_ring(connectivity, AtomIdx(0)) - assert not _in_ring(connectivity, AtomIdx(mol.num_atoms() - 1)) - assert not _in_ring(connectivity, AtomIdx(5), AtomIdx(6)) + assert connectivity.in_ring(AtomIdx(0)) + assert not connectivity.in_ring(AtomIdx(mol.num_atoms() - 1)) + assert not connectivity.in_ring(AtomIdx(5), AtomIdx(6)) def test_physical_connectivity_ring(): @@ -1019,7 +1016,7 @@ def test_physical_connectivity_ring(): import sire.morph as morph from sire.legacy.Mol import AtomIdx - from ghostly._ghostly import _create_connectivity, _in_ring, _physical_connectivity + from ghostly._ghostly import _create_connectivity, _physical_connectivity # Toluene to pyridine: at lambda = 1 the methyl carbon is a ghost bonded to # the ring, which stays a ring once the bonds to ghosts are removed. @@ -1032,7 +1029,7 @@ def test_physical_connectivity_ring(): ] physical = _physical_connectivity(mol, connectivity, ghosts) assert not physical.are_bonded(AtomIdx(0), AtomIdx(1)) - assert _in_ring(physical, AtomIdx(1)) + assert physical.in_ring(AtomIdx(1)) def test_to_rdkit_end_state(): @@ -1068,15 +1065,15 @@ def test_ring_closed_by_ghosts(): from sire.legacy.Mol import AtomIdx - from ghostly._ghostly import _in_ring, _physical_connectivity + from ghostly._ghostly import _physical_connectivity mol = sr.smiles("C1CCCCC1") connectivity = mol.property("connectivity") physical = _physical_connectivity(mol, connectivity, [AtomIdx(2), AtomIdx(3)]) - assert _in_ring(connectivity, AtomIdx(0)) - assert not _in_ring(physical, AtomIdx(0)) - assert not _in_ring(physical, AtomIdx(0), AtomIdx(1)) + assert connectivity.in_ring(AtomIdx(0)) + assert not physical.in_ring(AtomIdx(0)) + assert not physical.in_ring(AtomIdx(0), AtomIdx(1)) def test_linearise_ring_break_without_chirality(monkeypatch): From 2544091967a48cfb4f352b5e3c9992f037b0cf58 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 11 Nov 2026 09:28:37 +0000 Subject: [PATCH 16/20] Only flag rotatable bonds crossed by a surviving anchor dihedral. --- src/ghostly/_ghostly.py | 41 +++++++++++++++++++++++++++++++++++------ 1 file changed, 35 insertions(+), 6 deletions(-) diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index 68db07f..c578e59 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -2603,6 +2603,8 @@ def _modify_mixed_dihedrals( physical, connectivity, ref_connectivity, + mol.property("dihedral1" if is_lambda1 else "dihedral0"), + ghosts, is_lambda1=is_lambda1, warn=not stiffen_rotamers, ) @@ -2716,13 +2718,22 @@ def _modify_mixed_dihedrals( def _find_rotatable_bonds( - ref, bridges, physical, connectivity, ref_connectivity, is_lambda1=False, warn=True + ref, + bridges, + physical, + connectivity, + ref_connectivity, + dihedrals, + ghosts, + is_lambda1=False, + warn=True, ): r""" Find the rotatable bridge--physical bonds, i.e. those that are not in a - ring and whose bridge atom is neither planar nor linear in the end state - where the ghosts are physical, so that surviving anchor dihedrals across - them can allow rotameric transitions of ghost atoms. + ring, whose bridge atom is neither planar nor linear in the end state + where the ghosts are physical, and that are the central bond of a + surviving anchor dihedral, which can allow rotameric transitions of ghost + atoms. R1 DR1 \ / @@ -2751,6 +2762,13 @@ def _find_rotatable_bonds( ref_connectivity : sire.legacy.Mol.Connectivity The connectivity of the molecule at the opposite end state. + dihedrals : sire.legacy.MM.FourAtomFunctions + The dihedrals at the end state of the bridges, after the junction + modifications. + + ghosts : set of sire.legacy.Mol.AtomIdx + The ghost atoms at the end state of the bridges. + is_lambda1 : bool, optional Whether the bridges are at lambda = 1. @@ -2770,6 +2788,15 @@ def _find_rotatable_bonds( if not bridges: return rotatable_bonds + info = ref.info() + anchored = set() + for p in dihedrals.potentials(): + if p.function().is_zero(): + continue + idx = [info.atom_idx(a) for a in (p.atom0(), p.atom1(), p.atom2(), p.atom3())] + if idx[0] in ghosts or idx[3] in ghosts: + anchored.add(frozenset((idx[1], idx[2]))) + for bridge in bridges: # The ghosts are physical at the opposite end state. if not _is_rotor(ref, bridge, ref_connectivity, not is_lambda1): @@ -2778,8 +2805,10 @@ def _find_rotatable_bonds( for p in physical[bridge]: # Ring membership from the connectivity, since the physical # neighbour may be a ghost at the opposite end state. - if not connectivity.are_bonded(p, bridge) or connectivity.in_ring( - p, bridge + if ( + not connectivity.are_bonded(p, bridge) + or connectivity.in_ring(p, bridge) + or frozenset((p, bridge)) not in anchored ): continue From 0d456c06370385aa0b9985cdfe18f81daa1edaca Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 11 Nov 2026 09:36:03 +0000 Subject: [PATCH 17/20] Stiffen all ghost angles at linear dual bridges. --- CHANGELOG.md | 2 +- src/ghostly/_ghostly.py | 9 ++++++- tests/test_ghostly.py | 52 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 61 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f4b235eff..79826f9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,7 +14,7 @@ Changelog * Decide triple junction planarity from the energy needed to flatten the angle terms between the physical neighbours of the bridge, rather than its RDKit hybridisation. Junctions with a missing angle term are treated as non-planar [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Decide whether a bridge is in a ring, for stiffening and keeping impropers, from the physical atoms of the end state, so rings closed only by ghost atoms don't count [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Keep improper dihedrals that hold a single ghost in the plane of a ring bridge, since nothing else stops the ghost crossing the plane when the angles aren't stiffened, and remove the ``--stiffen-ring-bridges`` option [#23](https://github.com/OpenBioSim/ghostly/pull/23). -* Always stiffen ghost angles at bridges that aren't in a ring, and remove the ``--stiffen-sp2-bridges`` option. Without impropers, the angle terms of an sp2 bridge don't hold a ghost in the plane, so skipping stiffening let it cross freely [#23](https://github.com/OpenBioSim/ghostly/pull/23). +* Always stiffen ghost angles at bridges that aren't in a ring, and remove the ``--stiffen-sp2-bridges`` option. Without impropers, the angle terms of an sp2 bridge don't hold a ghost in the plane, so skipping stiffening let it cross freely. At linear bridges, angles through poorly-scoring physical neighbours are stiffened too, since the ghost angles can only all be satisfied at 90 degrees [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Detect rotatable bridge bonds from the angle terms of the bridge, rather than its RDKit hybridisation, so that rotamer detection is consistent with junction planarity [#23](https://github.com/OpenBioSim/ghostly/pull/23). [2026.2.0](https://github.com/openbiosim/ghostly/compare/2026.1.0...2026.2.0) - Sep 2026 diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index c578e59..40b5b61 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -958,8 +958,15 @@ def _dual( bridge_in_ring = physical_connectivity.in_ring(bridge) + # A linear bridge can only satisfy its ghost angles if all are 90 degrees, + # so stiffen through every physical neighbour regardless of score. + from math import degrees + + terms = list(_angle_terms(mol, bridge, physical, is_lambda1).values()) + bridge_linear = len(terms) == 1 and degrees(terms[0].theta0()) >= 170 + # Score the physical neighbours. - if bridge_indices is not None: + if bridge_indices is not None and not bridge_linear: phys_scores = {p: _score_atom(mol, p, bridge_indices) for p in physical} best_phys_score = min(phys_scores.values()) diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index 732c570..5f77e05 100644 --- a/tests/test_ghostly.py +++ b/tests/test_ghostly.py @@ -697,6 +697,58 @@ def test_pfkfb3_48_to_47_skipped_neighbour(): assert str(p.function()) == original[idx] +def test_pfkfb3_48_to_47_linear_bridge(): + """ + Test that a linear dual bridge is stiffened through every physical + neighbour, even a poorly-scoring one, since its ghost angles can only all + be satisfied at 90 degrees. + """ + + from math import pi + + from sire.legacy.CAS import Symbol + from sire.legacy.MM import AmberAngle + from sire.legacy.Mol import AtomIdx, Element + + mols = sr.load_test_files("pfkfb3_48_47.s3") + + angles0 = mols[0].property("angle0") + angles0.set( + AtomIdx(13), + AtomIdx(27), + AtomIdx(28), + AmberAngle(100, pi).to_expression(Symbol("theta")), + ) + mol = ( + mols[0] + .edit() + .atom(28) + .set_property("element1", Element("N")) + .molecule() + .set_property("angle0", angles0) + ) + mols.update(mol.commit()) + + _, modifications = modify(mols) + + def key(idx): + i, j, k = idx.split(",") + return ",".join(sorted((i, k)) + [j]) + + stiffened = { + key(idx): params + for idx, params in modifications["lambda_0"]["stiffened_angles"].items() + } + assert set(stiffened) == { + "13,47,27", + "13,48,27", + "28,47,27", + "28,48,27", + "47,48,27", + } + assert stiffened["47,48,27"]["theta0"] == 180.0 + + def test_stiffen_rotamers_reference_geometry(): """ Test that stiffened rotamer anchor dihedrals have their minimum at the From b3adeb6e332763b0c6d253d6d1a98ccd47967440 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 11 Nov 2026 14:15:40 +0000 Subject: [PATCH 18/20] Remove ghost dihedrals and impropers through linear angles. --- CHANGELOG.md | 1 + src/ghostly/_ghostly.py | 104 +++++++++++++++++++++++++++++++++++++++- tests/test_ghostly.py | 60 +++++++++++++++++++++++ 3 files changed, 164 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 79826f9..ef18198 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,7 @@ Changelog * Keep improper dihedrals that hold a single ghost in the plane of a ring bridge, since nothing else stops the ghost crossing the plane when the angles aren't stiffened, and remove the ``--stiffen-ring-bridges`` option [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Always stiffen ghost angles at bridges that aren't in a ring, and remove the ``--stiffen-sp2-bridges`` option. Without impropers, the angle terms of an sp2 bridge don't hold a ghost in the plane, so skipping stiffening let it cross freely. At linear bridges, angles through poorly-scoring physical neighbours are stiffened too, since the ghost angles can only all be satisfied at 90 degrees [#23](https://github.com/OpenBioSim/ghostly/pull/23). * Detect rotatable bridge bonds from the angle terms of the bridge, rather than its RDKit hybridisation, so that rotamer detection is consistent with junction planarity [#23](https://github.com/OpenBioSim/ghostly/pull/23). +* Remove dihedrals and impropers involving ghost atoms that pass through an angle that is linear at the end state, e.g. anchor dihedrals at a terminal junction whose physical neighbour is a nitrile carbon. These terms are undefined for collinear atoms, where their forces diverge [#23](https://github.com/OpenBioSim/ghostly/pull/23). [2026.2.0](https://github.com/openbiosim/ghostly/compare/2026.1.0...2026.2.0) - Sep 2026 ---------------------------------------------------------------------------------------- diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index 40b5b61..c2005c6 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -48,6 +48,9 @@ # Triple junctions needing more than this energy (kcal/mol) to flatten are non-planar. _MAX_PLANAR_FLATTENING = 2.0 +# Angles with an equilibrium value of at least this many degrees are linear. +_LINEAR_ANGLE = 170.0 + def modify( system, @@ -482,6 +485,11 @@ def modify( mol, ghosts0, modifications, skip_ghosts=linearised0, is_lambda1=False ) + # Remove ghost dihedrals and impropers through angles that are linear. + mol = _remove_linear_ghost_dihedrals( + mol, ghosts0, modifications, is_lambda1=False + ) + # Soften, single-well or stiffen surviving mixed ghost/physical dihedrals. mol = _modify_mixed_dihedrals( mol, @@ -598,6 +606,11 @@ def modify( mol, ghosts1, modifications, skip_ghosts=linearised1, is_lambda1=True ) + # Remove ghost dihedrals and impropers through angles that are linear. + mol = _remove_linear_ghost_dihedrals( + mol, ghosts1, modifications, is_lambda1=True + ) + # Soften, single-well or stiffen surviving mixed ghost/physical dihedrals. mol = _modify_mixed_dihedrals( mol, @@ -963,7 +976,7 @@ def _dual( from math import degrees terms = list(_angle_terms(mol, bridge, physical, is_lambda1).values()) - bridge_linear = len(terms) == 1 and degrees(terms[0].theta0()) >= 170 + bridge_linear = len(terms) == 1 and degrees(terms[0].theta0()) >= _LINEAR_ANGLE # Score the physical neighbours. if bridge_indices is not None and not bridge_linear: @@ -2021,6 +2034,95 @@ def _remove_impropers( return mol +def _remove_linear_ghost_dihedrals(mol, ghosts, modifications, is_lambda1=False): + r""" + Remove dihedrals and impropers involving a ghost atom that pass through an + angle that is linear at the current end state, i.e. three consecutive atoms + of the term are the atoms of the angle. These terms are undefined when the + angle's atoms are collinear, where their forces diverge, e.g. an anchor + dihedral at a terminal junction whose physical neighbour is a nitrile + carbon. The ghosts remain held by their angle terms, so they can only + rotate about the linear axis, which doesn't couple to the physical system. + + DR + \ + X---C≡N---R DR-X-C-R removed if X-C-R is linear + + Parameters + ---------- + + mol : sire.mol.Molecule + The perturbable molecule. + + ghosts : List[sire.legacy.Mol.AtomIdx] + The list of ghost atoms at the current end state. + + modifications : dict + A dictionary to store details of the modifications made. + + is_lambda1 : bool, optional + Whether to modify the terms at lambda = 1. + + Returns + ------- + + mol : sire.mol.Molecule + The updated molecule. + """ + from math import degrees + + from sire.legacy.CAS import Symbol + + # Nothing to do if there are no ghost atoms. + if not ghosts: + return mol + + info = mol.info() + suffix = "1" if is_lambda1 else "0" + mod_key = "lambda_1" if is_lambda1 else "lambda_0" + ghosts = set(ghosts) + + theta = Symbol("theta") + linear = set() + for p in mol.property("angle" + suffix).potentials(): + try: + theta0 = degrees(_SireMM.AmberAngle(p.function(), theta).theta0()) + except Exception: + continue + if theta0 >= _LINEAR_ANGLE: + linear.add( + frozenset(info.atom_idx(a) for a in (p.atom0(), p.atom1(), p.atom2())) + ) + + if not linear: + return mol + + for name in ("dihedral", "improper"): + terms = mol.property(name + suffix) + new_terms = _SireMM.FourAtomFunctions(info) + removed = False + for p in terms.potentials(): + idxs = [ + info.atom_idx(a) for a in (p.atom0(), p.atom1(), p.atom2(), p.atom3()) + ] + if any(idx in ghosts for idx in idxs) and ( + frozenset(idxs[:3]) in linear or frozenset(idxs[1:]) in linear + ): + dih_idx = ",".join(str(idx.value()) for idx in idxs) + _logger.debug( + f" Removing {name} through a linear angle: [{dih_idx}], {p.function()}" + ) + modifications[mod_key]["removed_dihedrals"].append(dih_idx) + removed = True + else: + new_terms.set(*idxs, p.function()) + + if removed: + mol = mol.edit().set_property(name + suffix, new_terms).molecule().commit() + + return mol + + def _remove_residual_ghost_dihedrals(mol, ghosts, modifications, is_lambda1=False): r""" Remove dihedral terms that couple ghost and physical regions but were diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index 5f77e05..cbb763c 100644 --- a/tests/test_ghostly.py +++ b/tests/test_ghostly.py @@ -1249,6 +1249,66 @@ def check_angle(info, potentials, idx0, idx1, idx2): return False +def test_hexane_to_propane_linear_anchor(): + """ + Test that an anchor dihedral through an angle that is linear at the end + state where its ghost is a ghost is removed, as for a terminal junction at + a nitrile. + """ + from math import pi + + from sire.legacy.CAS import Symbol + from sire.legacy.MM import AmberAngle + + mols = sr.load_test_files("hex2prp.s3") + info = mols[0].info() + ghosts = { + i for i, e in enumerate(mols[0].property("element1")) if e.num_protons() == 0 + } + + # Find an anchor dihedral that survives at lambda = 1. + new_mols, _ = modify(mols) + anchor = None + for p in new_mols[0].property("dihedral1").potentials(): + idx = [ + info.atom_idx(a).value() + for a in (p.atom0(), p.atom1(), p.atom2(), p.atom3()) + ] + if not p.function().is_zero() and (idx[0] in ghosts) != (idx[3] in ghosts): + if not ghosts & set(idx[1:3]): + anchor = idx + break + assert anchor is not None + angle = anchor[1:] if anchor[0] in ghosts else anchor[:3] + + # Make the physical angle the anchor passes through linear at lambda = 1. + theta = Symbol("theta") + angles = mols[0].property("angle1") + for p in angles.potentials(): + idx = [info.atom_idx(a) for a in (p.atom0(), p.atom1(), p.atom2())] + if [i.value() for i in idx] in (angle, angle[::-1]): + k = AmberAngle(p.function(), theta).k() + angles.set(*idx, AmberAngle(k, pi).to_expression(theta)) + mol = mols[0].edit().set_property("angle1", angles).molecule().commit() + mols.update(mol) + + new_mols, modifications = modify(mols) + + for p in new_mols[0].property("dihedral1").potentials(): + idx = [ + info.atom_idx(a).value() + for a in (p.atom0(), p.atom1(), p.atom2(), p.atom3()) + ] + if ghosts & set(idx): + assert set(angle) not in (set(idx[:3]), set(idx[1:])) + + removed = modifications["lambda_1"]["removed_dihedrals"] + assert ( + ",".join(map(str, anchor)) in removed + or ",".join(map(str, anchor[::-1])) in removed + ) + + def check_dihedral(info, potentials, idx0, idx1, idx2, idx3): """ Check if a dihedral potential is in a list of potentials. From a3053cc3cedd034952e4b4a40aa23e40778054e8 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 11 Nov 2026 14:22:41 +0000 Subject: [PATCH 19/20] Avoid linear anchors and extend the linear ghost term tests. --- src/ghostly/_ghostly.py | 40 ++++++++++++------- tests/test_ghostly.py | 86 +++++++++++++++++++++++++---------------- 2 files changed, 77 insertions(+), 49 deletions(-) diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index c2005c6..86e397b 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -847,13 +847,28 @@ def _terminal( # Initialise a container to store the updated dihedrals. new_dihedrals = _SireMM.FourAtomFunctions(mol.info()) + # Avoid anchors reached through a linear angle where possible, since + # dihedrals through linear angles are removed later. + from math import degrees + + candidates = [ + c + for c in physical2 + if not any( + degrees(term.theta0()) >= _LINEAR_ANGLE + for p in physical + if connectivity.are_bonded(p, c) + for term in _angle_terms(mol, p, [bridge, c], is_lambda1).values() + ) + ] or physical2 + # Select the best anchor atom and remove it from the list. Dihedrals # through the remaining physical2 atoms will be removed; dihedrals through # the anchor are kept to constrain the ghost group's orientation. - if bridge_indices is not None and len(physical2) > 1: - anchor = _select_anchor(mol, physical2, bridge_indices) + if bridge_indices is not None and len(candidates) > 1: + anchor = _select_anchor(mol, candidates, bridge_indices) else: - anchor = physical2[0] + anchor = candidates[0] physical2.remove(anchor) _logger.debug(f" Anchor atom: {anchor.value()}") for p in dihedrals.potentials(): @@ -2037,16 +2052,8 @@ def _remove_impropers( def _remove_linear_ghost_dihedrals(mol, ghosts, modifications, is_lambda1=False): r""" Remove dihedrals and impropers involving a ghost atom that pass through an - angle that is linear at the current end state, i.e. three consecutive atoms - of the term are the atoms of the angle. These terms are undefined when the - angle's atoms are collinear, where their forces diverge, e.g. an anchor - dihedral at a terminal junction whose physical neighbour is a nitrile - carbon. The ghosts remain held by their angle terms, so they can only - rotate about the linear axis, which doesn't couple to the physical system. - - DR - \ - X---C≡N---R DR-X-C-R removed if X-C-R is linear + angle that is linear at the current end state, since they are undefined, + and their forces diverge, when the angle's atoms are collinear. Parameters ---------- @@ -2088,6 +2095,9 @@ def _remove_linear_ghost_dihedrals(mol, ghosts, modifications, is_lambda1=False) try: theta0 = degrees(_SireMM.AmberAngle(p.function(), theta).theta0()) except Exception: + _logger.debug( + f" Skipping non-harmonic angle in linear check: {p.function()}" + ) continue if theta0 >= _LINEAR_ANGLE: linear.add( @@ -3026,7 +3036,7 @@ def _flattening_energy(mol, bridge, atoms, is_lambda1=False): return shortfall**2 / sum(1 / term.k() for term in terms.values()) -def _is_rotor(mol, bridge, connectivity, is_lambda1=False, max_angle=170): +def _is_rotor(mol, bridge, connectivity, is_lambda1=False, max_angle=_LINEAR_ANGLE): """ Return whether rotation about the bonds of the bridge atom can move its other substituents at an end state, i.e. the bridge is neither planar, by @@ -3063,7 +3073,7 @@ def _physical_state(mol, is_lambda1=False): return _morph.link_to_perturbed(mol) -def _single_well(ref, idxs, k, ref_ghosts, max_angle=170): +def _single_well(ref, idxs, k, ref_ghosts, max_angle=_LINEAR_ANGLE): """ Return a single-well n=1 cosine dihedral term, k[1 + cos(phi - phi_ref + pi)], with its minimum at the dihedral angle phi_ref of the atoms ``idxs`` in the diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index cbb763c..dbe84fc 100644 --- a/tests/test_ghostly.py +++ b/tests/test_ghostly.py @@ -1,3 +1,4 @@ +import pytest import sire as sr from ghostly import modify @@ -1249,64 +1250,81 @@ def check_angle(info, potentials, idx0, idx1, idx2): return False -def test_hexane_to_propane_linear_anchor(): +@pytest.mark.parametrize( + "file, state, name", + [ + ("hex2prp.s3", 1, "dihedral"), + ("acepol.s3", 0, "dihedral"), + ("tol2pyr.s3", 1, "improper"), + ], +) +def test_linear_ghost_terms(file, state, name): """ - Test that an anchor dihedral through an angle that is linear at the end - state where its ghost is a ghost is removed, as for a terminal junction at - a nitrile. + Test that a ghost dihedral or improper through an angle that is linear at + the end state where the ghost is a ghost is removed, as for an anchor at a + terminal junction on a nitrile. """ from math import pi from sire.legacy.CAS import Symbol from sire.legacy.MM import AmberAngle - mols = sr.load_test_files("hex2prp.s3") + mols = sr.load_test_files(file) info = mols[0].info() ghosts = { - i for i, e in enumerate(mols[0].property("element1")) if e.num_protons() == 0 + i + for i, e in enumerate(mols[0].property(f"element{state}")) + if e.num_protons() == 0 } - # Find an anchor dihedral that survives at lambda = 1. + def atoms(p, n): + return [info.atom_idx(getattr(p, f"atom{i}")()).value() for i in range(n)] + + # Find a surviving ghost term that passes through an angle term. new_mols, _ = modify(mols) - anchor = None - for p in new_mols[0].property("dihedral1").potentials(): - idx = [ - info.atom_idx(a).value() - for a in (p.atom0(), p.atom1(), p.atom2(), p.atom3()) - ] - if not p.function().is_zero() and (idx[0] in ghosts) != (idx[3] in ghosts): - if not ghosts & set(idx[1:3]): - anchor = idx + angles = { + frozenset(atoms(p, 3)): atoms(p, 3) + for p in new_mols[0].property(f"angle{state}").potentials() + } + term, angle = None, None + for p in new_mols[0].property(f"{name}{state}").potentials(): + idx = atoms(p, 4) + if p.function().is_zero() or not ghosts & set(idx): + continue + for triple in (frozenset(idx[:3]), frozenset(idx[1:])): + if triple in angles: + term, angle = idx, angles[triple] break - assert anchor is not None - angle = anchor[1:] if anchor[0] in ghosts else anchor[:3] + if term is not None: + break + assert term is not None - # Make the physical angle the anchor passes through linear at lambda = 1. + # Make that angle linear at the end state. theta = Symbol("theta") - angles = mols[0].property("angle1") - for p in angles.potentials(): + old = mols[0].property(f"angle{state}") + new = old.__class__(info) + linearised = False + for p in old.potentials(): idx = [info.atom_idx(a) for a in (p.atom0(), p.atom1(), p.atom2())] + function = p.function() if [i.value() for i in idx] in (angle, angle[::-1]): - k = AmberAngle(p.function(), theta).k() - angles.set(*idx, AmberAngle(k, pi).to_expression(theta)) - mol = mols[0].edit().set_property("angle1", angles).molecule().commit() + k = AmberAngle(function, theta).k() + function = AmberAngle(k, pi).to_expression(theta) + linearised = True + new.set(*idx, function) + assert linearised + mol = mols[0].edit().set_property(f"angle{state}", new).molecule().commit() mols.update(mol) new_mols, modifications = modify(mols) - for p in new_mols[0].property("dihedral1").potentials(): - idx = [ - info.atom_idx(a).value() - for a in (p.atom0(), p.atom1(), p.atom2(), p.atom3()) - ] + for p in new_mols[0].property(f"{name}{state}").potentials(): + idx = atoms(p, 4) if ghosts & set(idx): assert set(angle) not in (set(idx[:3]), set(idx[1:])) - removed = modifications["lambda_1"]["removed_dihedrals"] - assert ( - ",".join(map(str, anchor)) in removed - or ",".join(map(str, anchor[::-1])) in removed - ) + removed = modifications[f"lambda_{state}"]["removed_dihedrals"] + assert ",".join(map(str, term)) in removed def check_dihedral(info, potentials, idx0, idx1, idx2, idx3): From ae91e522de6f045ed08bb58dfcab3c628316f6e9 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 11 Nov 2026 14:28:42 +0000 Subject: [PATCH 20/20] Test that terminal junctions avoid linear anchors. --- tests/test_ghostly.py | 53 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index dbe84fc..ec506b8 100644 --- a/tests/test_ghostly.py +++ b/tests/test_ghostly.py @@ -1327,6 +1327,59 @@ def atoms(p, n): assert ",".join(map(str, term)) in removed +def test_hexane_to_propane_avoid_linear_anchor(): + """ + Test that a terminal junction is anchored through a candidate that isn't + reached through a linear angle, when there is one. At lambda = 1, bridge 3 + is anchored through 2 to 1, and 2 is also bonded to physical atoms 11 and 12. + """ + from math import pi + + from sire.legacy.CAS import Symbol + from sire.legacy.MM import AmberAngle + + mols = sr.load_test_files("hex2prp.s3") + info = mols[0].info() + ghosts = { + i for i, e in enumerate(mols[0].property("element1")) if e.num_protons() == 0 + } + + def anchors(mol): + result = set() + for p in mol.property("dihedral1").potentials(): + idx = [info.atom_idx(getattr(p, f"atom{i}")()).value() for i in range(4)] + if p.function().is_zero() or ghosts & set(idx[1:3]): + continue + if (idx[0] in ghosts) != (idx[3] in ghosts): + g, b, p1, p2 = idx if idx[0] in ghosts else idx[::-1] + result.add((b, p1, p2)) + return result + + new_mols, _ = modify(mols) + assert anchors(new_mols[0]) == {(3, 2, 1)} + + # Make the angle through the current anchor linear at lambda = 1. + theta = Symbol("theta") + old = mols[0].property("angle1") + new = old.__class__(info) + linearised = False + for p in old.potentials(): + idx = [info.atom_idx(getattr(p, f"atom{i}")()) for i in range(3)] + function = p.function() + if [i.value() for i in idx] in ([3, 2, 1], [1, 2, 3]): + k = AmberAngle(function, theta).k() + function = AmberAngle(k, pi).to_expression(theta) + linearised = True + new.set(*idx, function) + assert linearised + mols.update(mols[0].edit().set_property("angle1", new).molecule().commit()) + + new_mols, _ = modify(mols) + kept = anchors(new_mols[0]) + assert kept + assert all(b == 3 and p1 == 2 and p2 in (11, 12) for b, p1, p2 in kept) + + def check_dihedral(info, potentials, idx0, idx1, idx2, idx3): """ Check if a dihedral potential is in a list of potentials.