diff --git a/CHANGELOG.md b/CHANGELOG.md index f29cd59..ef18198 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,13 @@ 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 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 [#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. 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/README.md b/README.md index 214002f..cabc0a7 100644 --- a/README.md +++ b/README.md @@ -28,20 +28,24 @@ 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 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, 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 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 λ. -- **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 7d4cdff..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. """, @@ -202,34 +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", - 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, @@ -434,8 +406,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 300be26..86e397b 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -45,6 +45,12 @@ del _sys +# 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, @@ -57,8 +63,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, ): @@ -139,7 +143,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 @@ -151,21 +155,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 - ~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 @@ -239,8 +228,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, } @@ -395,6 +382,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. @@ -423,8 +424,7 @@ def modify( modifications, k_hard=k_hard, bridge_indices=bridge_indices0, - stiffen_ring_bridges=stiffen_ring_bridges, - stiffen_sp2_bridges=stiffen_sp2_bridges, + physical_connectivity=physical_connectivity0, ) # Triple junction. @@ -441,8 +441,7 @@ 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, + physical_connectivity=physical_connectivity0, ) # Higher order junction. @@ -458,11 +457,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, 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( @@ -488,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, @@ -495,6 +497,8 @@ def modify( physical0, ghosts0, ghosts1, + connectivity0, + connectivity1, modifications, is_lambda1=False, soften_anchors=soften_anchors, @@ -505,6 +509,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]) @@ -531,8 +544,7 @@ 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, + physical_connectivity=physical_connectivity1, ) elif junction == 3: @@ -549,8 +561,7 @@ 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, + physical_connectivity=physical_connectivity1, ) # Higher order junction. @@ -567,11 +578,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, 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( @@ -597,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, @@ -604,6 +618,8 @@ def modify( physical1, ghosts1, ghosts0, + connectivity1, + connectivity0, modifications, is_lambda1=True, soften_anchors=soften_anchors, @@ -831,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(): @@ -878,8 +909,7 @@ def _dual( k_hard=100, is_lambda1=False, bridge_indices=None, - stiffen_ring_bridges=False, - stiffen_sp2_bridges=False, + physical_connectivity=None, ): r""" Apply modifications to a dual junction. @@ -954,43 +984,17 @@ def _dual( mod_key = "lambda_0" 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) + bridge_in_ring = physical_connectivity.in_ring(bridge) - from sire.convert import to_rdkit - - try: - rdmol = to_rdkit(end_state_mol) - except Exception: - rdmol = None - - bridge_in_ring = ( - rdmol is not None and rdmol.GetAtomWithIdx(bridge.value()).IsInRing() - ) - - # Check sp2 hybridisation for non-ring bridges. - bridge_is_sp2 = False - if rdmol is not None and not bridge_in_ring: - from rdkit.Chem.rdchem import HybridizationType + # 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 - hyb = rdmol.GetAtomWithIdx(bridge.value()).GetHybridization() - bridge_is_sp2 = hyb == HybridizationType.SP2 - - # 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." - ) + terms = list(_angle_terms(mol, bridge, physical, is_lambda1).values()) + bridge_linear = len(terms) == 1 and degrees(terms[0].theta0()) >= _LINEAR_ANGLE # 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()) @@ -1073,11 +1077,8 @@ 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 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( f" Skipping stiffening for angle " @@ -1086,18 +1087,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. @@ -1169,10 +1158,8 @@ 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" elif phys_scores is not None and any( phys_scores[p] > best_phys_score and (heavy_phys - {p}) for p in physical ): @@ -1272,8 +1259,7 @@ 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, + physical_connectivity=physical_connectivity, ) # Return the updated molecule. @@ -1293,8 +1279,7 @@ def _triple( num_optimise=10, is_lambda1=False, bridge_indices=None, - stiffen_ring_bridges=False, - stiffen_sp2_bridges=False, + physical_connectivity=None, ): r""" Apply modifications to a triple junction. @@ -1375,41 +1360,31 @@ 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 + bridge_in_ring = physical_connectivity.in_ring(bridge) - # Convert the molecule to RDKit to determine the hybridisation of the bridge atom. - try: - rdmol = to_rdkit(end_state_mol) - 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()) - hybridisation = bridge_rdatom.GetHybridization() - bridge_in_ring = bridge_rdatom.IsInRing() - - # 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." - ) + # 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 {'planar' if is_planar else 'non-planar'} junction." + ) + else: + _logger.debug(f" Flattening energy: {flattening:.3f} kcal/mol") + is_planar = flattening < _MAX_PLANAR_FLATTENING # Score the physical neighbours. if bridge_indices is not None: @@ -1429,11 +1404,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 @@ -1710,18 +1681,13 @@ 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 and not stiffen_ring_bridges) or ( - bridge_is_sp2 and not stiffen_sp2_bridges - ): + # Ring bridges aren't stiffened, so keep the sacrificial angle. + 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 @@ -1775,8 +1741,7 @@ 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, + physical_connectivity=physical_connectivity, ) # Return the updated molecule. @@ -1795,6 +1760,7 @@ def _higher( optimise_angles=False, num_optimise=10, is_lambda1=False, + physical_connectivity=None, ): r""" Apply modifications to higher order junctions. @@ -1858,6 +1824,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() @@ -1980,12 +1949,17 @@ def _higher( optimise_angles=optimise_angles, num_optimise=num_optimise, is_lambda1=is_lambda1, + physical_connectivity=physical_connectivity, ) -def _remove_impropers(mol, ghosts, 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. + 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 ---------- @@ -1996,6 +1970,12 @@ 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. + + 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. @@ -2029,11 +2009,27 @@ def _remove_impropers(mol, ghosts, modifications, is_lambda1=False): 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 physical_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()}" ) @@ -2053,6 +2049,90 @@ def _remove_impropers(mol, ghosts, modifications, is_lambda1=False): 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, since they are undefined, + and their forces diverge, when the angle's atoms are collinear. + + 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: + _logger.debug( + f" Skipping non-harmonic angle in linear check: {p.function()}" + ) + 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 @@ -2270,35 +2350,48 @@ 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) - linearised = set() - 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(ghost.value()) - chiral = rdatom.GetChiralTag() - if chiral in ( + + # 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"Unable to check chirality at {_lam_sym} = {int(is_lambda1)}, so " + f"ring-break ghosts that could be stereocentres aren't linearised: {e}" + ) + 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 @@ -2525,6 +2618,8 @@ def _modify_mixed_dihedrals( physical, ghosts, ref_ghosts, + connectivity, + ref_connectivity, modifications, is_lambda1=False, soften_anchors=1.0, @@ -2543,7 +2638,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. @@ -2557,7 +2652,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). @@ -2622,7 +2717,15 @@ 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, + mol.property("dihedral1" if is_lambda1 else "dihedral0"), + ghosts, + is_lambda1=is_lambda1, + warn=not stiffen_rotamers, ) # Nothing to do if no modifications are requested. @@ -2733,16 +2836,28 @@ 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, + dihedrals, + ghosts, + 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, 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 \ / \ / - X ---DR2 X is sp3, P-X is rotatable + X ---DR2 X is non-planar, P-X is rotatable / \ / \ R2 DR3 @@ -2760,6 +2875,19 @@ 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. + + 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. @@ -2779,43 +2907,161 @@ def _find_rotatable_bonds(ref, bridges, physical, is_lambda1=False, warn=True): if not bridges: return rotatable_bonds - from rdkit.Chem import HybridizationType - from sire.convert import to_rdkit - - try: - rdmol = to_rdkit(ref) - except Exception as e: - _logger.warning(f"Failed to convert molecule to RDKit for rotamer check: {e}") - 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: - 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, ref_connectivity, not is_lambda1): + continue - rd_bond = rdmol.GetBondBetweenAtoms(p_idx, b_idx) - if rd_bond is None or rd_bond.IsInRing(): + 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) + or frozenset((p, bridge)) not in anchored + ): continue - hybridisation = rdmol.GetAtomWithIdx(b_idx).GetHybridization() - - 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 +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_connectivity(mol, connectivity, ghosts): + """ + Return the connectivity without the bonds to the given ghost atoms. + """ + info = mol.info() + 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 idx0 not in ghosts and idx1 not in ghosts: + physical.connect(idx0, idx1) + return physical.commit() + + +def _angle_terms(mol, bridge, atoms, is_lambda1=False): + """ + 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 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) + 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 + + # 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 _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 + its angle terms, nor linear. + """ + from math import degrees + + elements = mol.property("element1" if is_lambda1 else "element0") + neighbours = [ + n + for n in connectivity.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): """ Return the molecule linked to the end state where the ghost atoms of the @@ -2827,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 55e40bf..ec506b8 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 @@ -124,10 +125,8 @@ 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 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)), @@ -154,6 +153,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(): """ @@ -175,9 +183,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") @@ -393,9 +400,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. 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, 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)), @@ -405,46 +411,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 = [ @@ -717,6 +698,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 @@ -954,6 +987,205 @@ 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. + """ + + from sire.legacy.Mol import AtomIdx + + connectivity = sr.smiles("C1CCCCCCCCCCCCC1").property("connectivity") + assert connectivity.in_ring(AtomIdx(0)) + assert connectivity.in_ring(AtomIdx(0), AtomIdx(1)) + + mol = sr.smiles("C1CCCCC1CC") + connectivity = mol.property("connectivity") + 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(): + """ + 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, _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)) + 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 physical.in_ring(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_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 _physical_connectivity + + mol = sr.smiles("C1CCCCC1") + connectivity = mol.property("connectivity") + physical = _physical_connectivity(mol, connectivity, [AtomIdx(2), AtomIdx(3)]) + + 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): + """ + 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): + raise RuntimeError("no bond orders") + + monkeypatch.setattr(_ghostly, "_to_rdkit_end_state", fail) + + 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} + + # 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): """ Return a dictionary mapping dihedral indices to their functions. @@ -1018,6 +1250,136 @@ def check_angle(info, potentials, idx0, idx1, idx2): return False +@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 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(file) + info = mols[0].info() + ghosts = { + i + for i, e in enumerate(mols[0].property(f"element{state}")) + if e.num_protons() == 0 + } + + 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) + 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 + if term is not None: + break + assert term is not None + + # Make that angle linear at the end state. + theta = Symbol("theta") + 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(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(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[f"lambda_{state}"]["removed_dihedrals"] + 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.