vibe-sketcherDocumentation
Source code for vibesketcher.chem.endpoints
"""The two ends of a reaction path, built from a sketch and its atom mapping (#29).
A nudged elastic band needs a reactant complex and a product complex with
the same atoms in the same order and as little needless motion between them
as possible, or the band spends its images rotating molecules. Both
complexes are built the same way:
* each side is embedded as one molecule with ETKDG version 3 (S. Riniker,
G. A. Landrum, *J. Chem. Inf. Model.* 55, 2562-2574 (2015)), the bonds the
reaction will form (on the reactant side) or has broken (on the product
side) held at 1.85 times the sum of the two atoms' covalent radii, a
loosely bound complex rather than a bond;
* an atom that gains one partner and loses another keeps the two opposite
each other (the incoming and the leaving group of a substitution);
* the complex is refined with MMFF94 (T. A. Halgren, *J. Comput. Chem.* 17,
490-519 (1996)), or UFF where MMFF94 has no parameters, with those
distances restrained, and must keep the drawn stereochemistry;
* the product complex is turned onto the reactant complex by a proper
rotation (W. Kabsch, *Acta Cryst.* A32, 922-923 (1976)); its mirror image
is tried too when no stereocentre is specified, and hydrogens that the
mapping leaves interchangeable (the three of a methyl group) are paired
by position, so the band does not twist them past each other.
The geometries are starting points for a calculation, not optima: ORCA's
NEB-TS optimises the ends first, and the vibe-qc script does the same.
"""
from __future__ import annotations
import itertools
from dataclasses import dataclass, field, replace
from typing import Any
import numpy as np
from numpy.typing import NDArray
from vibesketcher.chem.mapping import (
AtomMapping,
MappingFailed,
Side,
_hydrogen_groups,
_select_hydrogen_pins,
validate_correspondence,
)
__all__ = ["EndpointsFailed", "Endpoints", "build_endpoints"]
_Array = NDArray[np.float64]
#: Forming and breaking bonds, in units of the two covalent radii.
_STRETCH = 1.85
#: How far the restrained distances may stray, in angstrom.
_SLACK = 0.15
#: How far a conformation hint from the other end may be missed, in angstrom.
_HINT = 0.25
#: Embeddings tried per side.
_ATTEMPTS = 8
#: The seed every embedding starts from (as vibesketcher.chem.embed).
_SEED = 0x5EED
[docs]
class EndpointsFailed(ValueError):
"""The endpoint geometries could not be built; the message says why."""
[docs]
@dataclass(frozen=True)
class Endpoints:
"""Both ends of the path, atoms in the reactants' order, in angstrom."""
symbols: tuple[str, ...]
reactant: tuple[tuple[float, float, float], ...]
product: tuple[tuple[float, float, float], ...]
#: The mapping with hydrogens paired by position: product atom per reactant atom.
pairs: tuple[int, ...]
#: Root-mean-square displacement between the two ends, in angstrom.
rmsd: float
method: str
warnings: list[str] = field(default_factory=list)
# ── Embedding one side ──────────────────────────────────────────────────
def _target(mol: Any, a: int, b: int) -> float:
from rdkit import Chem
table = Chem.GetPeriodicTable()
radii = table.GetRcovalent(mol.GetAtomWithIdx(a).GetAtomicNum()) + table.GetRcovalent(
mol.GetAtomWithIdx(b).GetAtomicNum()
)
return float(_STRETCH * radii)
def _restraints(
mol: Any, stretched: list[tuple[int, int]], leaving: list[tuple[int, int]]
) -> tuple[dict[tuple[int, int], tuple[float, float]], list[tuple[int, int, float, float]]]:
"""Distance bounds for the stretched pairs and for opposite groups, and force-field restraints.
``leaving`` are bonds of this complex the reaction breaks or has formed:
an atom stretched towards one partner keeps its leaving partner opposite.
"""
from rdkit import Chem
table = Chem.GetPeriodicTable()
bounds: dict[tuple[int, int], tuple[float, float]] = {}
for a, b in stretched:
distance = _target(mol, a, b)
bounds[(min(a, b), max(a, b))] = (distance - _SLACK, distance + _SLACK)
# An atom with a stretched partner and a bonded one that leaves keeps them opposite.
bonded = {frozenset((bond.GetBeginAtomIdx(), bond.GetEndAtomIdx())) for bond in mol.GetBonds()}
leaves = {(min(a, b), max(a, b)) for a, b in leaving}
partners: dict[int, list[int]] = {}
for a, b in stretched:
partners.setdefault(a, []).append(b)
partners.setdefault(b, []).append(a)
for centre, incoming in partners.items():
for x in incoming:
for y in incoming:
if x >= y:
continue
if frozenset((centre, x)) in bonded or frozenset((centre, y)) in bonded:
continue
# Two stretched partners of one atom: keep them apart, not stacked.
far = 0.9 * (_target(mol, centre, x) + _target(mol, centre, y))
bounds[(min(x, y), max(x, y))] = (far, far + 3.0)
for x in incoming:
for neighbour in mol.GetAtomWithIdx(centre).GetNeighbors():
y = neighbour.GetIdx()
if y == x or neighbour.GetAtomicNum() == 1:
continue
if (min(centre, y), max(centre, y)) not in leaves:
continue
bond = table.GetRcovalent(mol.GetAtomWithIdx(centre).GetAtomicNum()) + (
table.GetRcovalent(neighbour.GetAtomicNum())
)
far = _target(mol, centre, x) + bond - 0.1 # within about 30 degrees of a line
bounds[(min(x, y), max(x, y))] = (far, far + 1.0)
restraints = [(a, b, low, high) for (a, b), (low, high) in sorted(bounds.items())]
return bounds, restraints
def _stereo_kept(side: Side, mol: Any, conformer: int) -> bool:
"""Whether every species of the side keeps its drawn stereochemistry in 3D."""
from rdkit import Chem
from vibesketcher.chem.embed import _stereo_kept as kept
from vibesketcher.chem.embed import _target as target
from vibesketcher.chem.species import mol_from_smiles
positions = mol.GetConformer(conformer).GetPositions()
for _, smiles, first in side.species:
heavy = mol_from_smiles(smiles)
with_h = Chem.AddHs(heavy)
count = with_h.GetNumAtoms()
conformer_copy = Chem.Conformer(count)
for i in range(count):
conformer_copy.SetAtomPosition(i, positions[first + i].tolist())
with_h.RemoveAllConformers()
identifier = with_h.AddConformer(conformer_copy, assignId=True)
if not kept(with_h, identifier, heavy, target(heavy)):
return False
return True
def _embed(
side: Side,
stretched: list[tuple[int, int]],
leaving: list[tuple[int, int]],
hints: dict[tuple[int, int], float] | None = None,
) -> tuple[list[_Array], str, list[str]]:
"""Coordinates of the side's complex, best first, and how they were made.
``hints`` are distances between atoms of one molecule the other end
already has; where the molecule's own bounds allow them, they steer its
conformation (a diene s-cis, as in the ring it closes).
"""
from rdkit import Chem, DistanceGeometry, rdBase
from rdkit.Chem import rdDistGeom
from rdkit.Chem import rdForceFieldHelpers as force_fields
from vibesketcher.chem._rdkit import untyped
geometry = untyped(rdDistGeom)
fields = untyped(force_fields)
mol = Chem.Mol(side.mol)
warnings: list[str] = []
if mol.GetNumAtoms() == 1:
return [np.zeros((1, 3))], "atom", warnings
bounds, restraints = _restraints(mol, stretched, leaving)
matrix = geometry.GetMoleculeBoundsMatrix(mol)
for (a, b), distance in (hints or {}).items():
low, high = matrix[b][a], matrix[a][b]
if (a, b) not in bounds and low <= distance <= high:
hint = (max(low, distance - _HINT), min(high, distance + _HINT))
bounds[(a, b)] = hint
restraints.append((a, b, *hint))
for (a, b), (low, high) in bounds.items(): # a < b: upper bound above, lower below
matrix[a][b] = high
matrix[b][a] = low
blocker = rdBase.BlockLogs()
try:
if not untyped(DistanceGeometry).DoTriangleSmoothing(matrix):
raise EndpointsFailed(
"the forming and breaking bonds cannot all be placed at once; pin atoms to "
"choose another mapping"
)
params = geometry.ETKDGv3()
params.randomSeed = _SEED
params.numThreads = 1
params.SetBoundsMat(matrix)
params.useRandomCoords = True
ids = list(geometry.EmbedMultipleConfs(mol, _ATTEMPTS, params))
if not ids:
raise EndpointsFailed("distance geometry found no 3D structure for the complex")
method = "etkdg"
energies: dict[int, float] = {}
if fields.MMFFHasAllMoleculeParams(mol):
method = "etkdg+mmff94"
properties = fields.MMFFGetMoleculeProperties(mol)
for conformer in ids:
field_ = fields.MMFFGetMoleculeForceField(mol, properties, confId=conformer)
for a, b, low, high in restraints:
field_.MMFFAddDistanceConstraint(a, b, False, low, high, 500.0)
field_.Minimize(maxIts=2000)
energies[conformer] = float(field_.CalcEnergy())
elif fields.UFFHasAllMoleculeParams(mol):
method = "etkdg+uff"
for conformer in ids:
field_ = fields.UFFGetMoleculeForceField(mol, confId=conformer)
for a, b, low, high in restraints:
field_.UFFAddDistanceConstraint(a, b, False, low, high, 500.0)
field_.Minimize(maxIts=2000)
energies[conformer] = float(field_.CalcEnergy())
else:
warnings.append(
"no force field covers the complex, so its geometry is distance geometry's alone"
)
kept = [c for c in ids if _stereo_kept(side, mol, c)]
finally:
del blocker
if not kept:
raise EndpointsFailed(
"no 3D structure of the complex kept the drawn stereochemistry; check the wedges"
)
kept.sort(key=lambda c: energies.get(c, 0.0))
coordinates = [np.asarray(mol.GetConformer(c).GetPositions(), dtype=np.float64) for c in kept]
return coordinates, method, warnings
# ── Putting the ends together ───────────────────────────────────────────
def _kabsch(moving: _Array, target: _Array) -> tuple[_Array, _Array, _Array]:
"""The proper rotation and the two centres that put ``moving`` onto ``target``."""
centre_a, centre_b = moving.mean(axis=0), target.mean(axis=0)
if len(moving) < 2:
return np.eye(3), centre_a, centre_b
u, _, vt = np.linalg.svd((moving - centre_a).T @ (target - centre_b))
d = float(np.sign(np.linalg.det(vt.T @ u.T))) or 1.0
rotation: _Array = vt.T @ np.diag([1.0, 1.0, d]) @ u.T
return rotation, centre_a, centre_b
def _apply(points: _Array, transform: tuple[_Array, _Array, _Array]) -> _Array:
rotation, centre_a, centre_b = transform
moved: _Array = (points - centre_a) @ rotation.T + centre_b
return moved
def _rmsd(a: _Array, b: _Array) -> float:
return float(np.sqrt(((a - b) ** 2).sum(axis=1).mean()))
def _repair_hydrogens(
mapping: AtomMapping, pairs: list[int], reactant: _Array, product: _Array
) -> list[int]:
"""Equivalent spectator hydrogens paired by position; declared transfers stay fixed."""
pairs = list(pairs)
pinned = dict(mapping.pinned)
for group in _hydrogen_groups(mapping):
hydrogens = [h for h in group if h not in pinned]
partners = [pairs[h] for h in hydrogens]
if len(hydrogens) < 2:
continue
if len(hydrogens) <= 6:
best = min(
itertools.permutations(partners),
key=lambda order: sum(
float(np.linalg.norm(reactant[h] - product[p]))
for h, p in zip(hydrogens, order, strict=True)
),
)
else: # greedy, nearest first
remaining = list(partners)
chosen = []
for h in hydrogens:
nearest = min(
remaining, key=lambda p: float(np.linalg.norm(reactant[h] - product[p]))
)
remaining.remove(nearest)
chosen.append(nearest)
best = tuple(chosen)
candidate = list(pairs)
for h, p in zip(hydrogens, best, strict=True):
candidate[h] = p
try:
validate_correspondence(mapping, candidate)
except MappingFailed:
continue
pairs = candidate
return pairs
def _changes(mapping: AtomMapping) -> tuple[list[tuple[int, int]], list[tuple[int, int]]]:
"""Bonds the reaction forms and breaks, in reactant atom numbers."""
r_bonds = {
frozenset((b.GetBeginAtomIdx(), b.GetEndAtomIdx()))
for b in mapping.reactants.mol.GetBonds()
}
back = {p: r for r, p in enumerate(mapping.pairs)}
p_bonds = {
frozenset((back[b.GetBeginAtomIdx()], back[b.GetEndAtomIdx()]))
for b in mapping.products.mol.GetBonds()
}
formed = [tuple(sorted(bond)) for bond in p_bonds - r_bonds]
broken = [tuple(sorted(bond)) for bond in r_bonds - p_bonds]
return [(a, b) for a, b in sorted(formed)], [(a, b) for a, b in sorted(broken)]
def _hints(mapping: AtomMapping, pairs: list[int], product: _Array) -> dict[tuple[int, int], float]:
"""Distances in the product between heavy atoms three or more bonds apart in one reactant."""
from rdkit import Chem
mol = mapping.reactants.mol
owners = mapping.reactants.owners
topology = Chem.GetDistanceMatrix(mol)
heavy = [a.GetIdx() for a in mol.GetAtoms() if a.GetAtomicNum() > 1]
found: dict[tuple[int, int], float] = {}
for i, a in enumerate(heavy):
for b in heavy[i + 1 :]:
if owners[a][0] == owners[b][0] and topology[a][b] >= 3:
found[(a, b)] = float(np.linalg.norm(product[pairs[a]] - product[pairs[b]]))
return found
def _specified_centres(side: Side) -> bool:
from rdkit import Chem
return any(
atom.GetChiralTag() != Chem.ChiralType.CHI_UNSPECIFIED for atom in side.mol.GetAtoms()
)
def _validate_stereo_geometry(
mapping: AtomMapping, pairs: list[int], reactant: _Array, product: _Array
) -> None:
"""Read both drawn endpoint stereochemistries back from the actual ordered rows.
Each side is checked independently: an intended SN2 inversion is allowed.
No reflection, rotation or coordinate repair happens during this check.
"""
from rdkit import Chem
back = {b: a for a, b in enumerate(pairs)}
native_product = product[[back[b] for b in range(len(pairs))]]
for side, points in ((mapping.reactants, reactant), (mapping.products, native_product)):
mol = Chem.Mol(side.mol)
conformer = Chem.Conformer(len(pairs))
for i, point in enumerate(points):
conformer.SetAtomPosition(i, point.tolist())
mol.RemoveAllConformers()
identifier = mol.AddConformer(conformer, assignId=True)
if not _stereo_kept(side, mol, identifier):
raise EndpointsFailed(
"the actual endpoint correspondence does not keep drawn stereochemistry"
)
[docs]
def build_endpoints(mapping: AtomMapping) -> Endpoints:
"""The reactant and product complexes for ``mapping``, aligned, in the reactants' order."""
if mapping.ambiguous:
raise EndpointsFailed("the atom mapping is ambiguous; pin atoms first")
requested = mapping
hydrogen_pins = any(
mapping.reactants.mol.GetAtomWithIdx(a).GetAtomicNum() == 1 for a, _ in mapping.pinned
)
if hydrogen_pins:
# Select exactly the same geometries/alignment as without H pins, then
# relabel their fixed rows. A correspondence choice is not a new embed,
# conformer score, rotation or physical-path optimization.
mapping = replace(
mapping,
pairs=mapping.unconstrained_pairs or mapping.pairs,
pinned=tuple(
(a, b)
for a, b in mapping.pinned
if mapping.reactants.mol.GetAtomWithIdx(a).GetAtomicNum() > 1
),
unconstrained_pairs=(),
)
formed, broken = _changes(mapping)
pairs = list(mapping.pairs)
# Product side numbers: what the reaction formed is bonded there, what it broke is stretched.
to_p = {r: p for r, p in enumerate(pairs)}
products, method_p, warnings = _embed(
mapping.products,
[(to_p[a], to_p[b]) for a, b in broken],
[(to_p[a], to_p[b]) for a, b in formed],
)
# The reactants take the product's conformation where their own bonds allow it.
reactants, method_r, more = _embed(
mapping.reactants, formed, broken, _hints(mapping, pairs, products[0])
)
warnings += more
reactant = reactants[0]
heavy = [
i for i, atom in enumerate(mapping.reactants.mol.GetAtoms()) if atom.GetAtomicNum() > 1
]
mirror = not _specified_centres(mapping.products)
best: tuple[float, _Array, list[int]] | None = None
for candidate in products:
for flip in (False, True) if mirror else (False,):
coordinates = candidate * np.array([1.0, 1.0, -1.0]) if flip else candidate
# Heavy atoms first: they fix the orientation, the hydrogens then pair by position.
ordered = coordinates[pairs]
atoms = heavy if len(heavy) >= 3 else list(range(len(pairs)))
aligned = _apply(coordinates, _kabsch(ordered[atoms], reactant[atoms]))
repaired = _repair_hydrogens(mapping, pairs, reactant, aligned)
moved = coordinates[repaired]
final = _apply(moved, _kabsch(moved, reactant))
score = _rmsd(final, reactant)
if best is None or score < best[0] - 1e-9:
best = (score, final, repaired)
assert best is not None
score, product, repaired = best
try:
validate_correspondence(mapping, repaired)
if hydrogen_pins:
selected = _select_hydrogen_pins(requested, repaired)
back = {b: a for a, b in enumerate(repaired)}
product = product[[back[b] for b in selected]]
repaired = list(selected)
score = _rmsd(product, reactant)
validate_correspondence(requested, repaired)
except MappingFailed as exc:
raise EndpointsFailed(str(exc)) from exc
_validate_stereo_geometry(requested, repaired, reactant, product)
symbols = tuple(atom.GetSymbol() for atom in mapping.reactants.mol.GetAtoms())
method = method_r if method_r == method_p else f"{method_r}/{method_p}"
def rows(points: _Array) -> tuple[tuple[float, float, float], ...]:
return tuple(
(round(float(x), 6), round(float(y), 6), round(float(z), 6)) for x, y, z in points
)
return Endpoints(
symbols, rows(reactant), rows(product), tuple(repaired), score, method, warnings
)