"""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
from typing import Any
import numpy as np
from numpy.typing import NDArray
from vibesketch.chem.mapping import AtomMapping, Side
__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 vibesketch.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 vibesketch.chem.embed import _stereo_kept as kept
from vibesketch.chem.embed import _target as target
from vibesketch.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 vibesketch.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]:
"""Hydrogens on each heavy atom pair, and the moving ones, paired by position."""
r_mol, p_mol = mapping.reactants.mol, mapping.products.mol
pairs = list(pairs)
groups: list[tuple[list[int], list[int]]] = []
moving_r: list[int] = []
moving_p_taken: set[int] = set()
for atom in r_mol.GetAtoms():
if atom.GetAtomicNum() == 1:
continue
a = atom.GetIdx()
on_a = [n.GetIdx() for n in atom.GetNeighbors() if n.GetAtomicNum() == 1]
on_b = [
n.GetIdx()
for n in p_mol.GetAtomWithIdx(pairs[a]).GetNeighbors()
if n.GetAtomicNum() == 1
]
stay = [h for h in on_a if pairs[h] in on_b]
groups.append((stay, [pairs[h] for h in stay]))
moving_r.extend(h for h in on_a if pairs[h] not in on_b)
moving_p_taken.update(pairs[h] for h in stay)
moving_r.extend(
a.GetIdx() for a in r_mol.GetAtoms() if a.GetAtomicNum() == 1 and not a.GetNeighbors()
)
groups.append((moving_r, [pairs[h] for h in moving_r]))
for hydrogens, partners in groups:
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)
for h, p in zip(hydrogens, best, strict=True):
pairs[h] = p
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()
)
[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")
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
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
)