"""3D structures for drawn species: ETKDG, a force field, the stereochemistry kept.
Coordinates come from RDKit's distance geometry with experimental torsion
preferences, ETKDG version 3 (S. Riniker, G. A. Landrum, *J. Chem. Inf.
Model.* 55, 2562-2574 (2015); S. Wang, J. Witek, G. A. Landrum, S. Riniker,
*J. Chem. Inf. Model.* 60, 2044-2058 (2020)). Several embeddings are refined
with MMFF94 (T. A. Halgren, *J. Comput. Chem.* 17, 490-519 (1996)), or with
UFF (A. K. Rappé et al., *J. Am. Chem. Soc.* 114, 10024-10035 (1992)) where
MMFF94 has no parameters, and ranked by the force field's energy. A species
neither covers, such as most metal complexes, keeps ETKDG's geometry with a
warning: distance geometry knows nothing of square-planar or octahedral
centres.
The drawing's stereochemistry is the contract. Every conformer is read back
(``AssignStereochemistryFrom3D``) and must give the R/S and E/Z the drawing
specifies; one that does not is dropped, and if none is left, embedding
fails with that reason instead of returning the wrong isomer (roadmap,
section 6). Stereocentres the drawing leaves open are counted, because the
3D structure has to pick one of them.
Embedding is deterministic: a fixed seed and one thread.
"""
from __future__ import annotations
from dataclasses import dataclass, field
from typing import Any
from vibesketch.document.model import Conformer, ConformerSource
__all__ = ["EmbedFailed", "Embedded", "embed_3d"]
#: The seed every embedding starts from, so the same drawing gives the same structure.
SEED = 0x5EED
[docs]
class EmbedFailed(ValueError):
"""No 3D structure with the drawing's stereochemistry could be made."""
[docs]
@dataclass(frozen=True)
class Embedded:
"""3D structures for one species, lowest force-field energy first."""
conformers: list[Conformer]
#: Force-field energies in kcal/mol, in the same order; ``None`` without a force field.
energies: list[float | None]
#: Stereocentres and stereogenic double bonds the drawing leaves open.
open_stereo: int = 0
warnings: list[str] = field(default_factory=list)
def _target(mol: Any) -> tuple[str, set[int], set[int]]:
"""The canonical isomeric SMILES, and the atoms and bonds with specified stereo."""
from rdkit import Chem
atoms = {
atom.GetIdx()
for atom in mol.GetAtoms()
if atom.GetChiralTag() != Chem.ChiralType.CHI_UNSPECIFIED
}
bonds = {
bond.GetIdx()
for bond in mol.GetBonds()
if bond.GetStereo() not in (Chem.BondStereo.STEREONONE, Chem.BondStereo.STEREOANY)
}
return str(Chem.MolToSmiles(mol)), atoms, bonds
def _stereo_kept(
with_h: Any, conformer_id: int, heavy: Any, target: tuple[str, set[int], set[int]]
) -> bool:
"""Whether a conformer gives back the stereochemistry the drawing specifies."""
from rdkit import Chem
smiles, atoms, bonds = target
read = Chem.Mol(with_h, confId=conformer_id)
Chem.AssignStereochemistryFrom3D(read)
# Hydrogens were appended, so the drawing's atoms and bonds keep their indices.
for atom in read.GetAtoms():
if atom.GetIdx() < heavy.GetNumAtoms() and atom.GetIdx() not in atoms:
atom.SetChiralTag(Chem.ChiralType.CHI_UNSPECIFIED)
for bond in read.GetBonds():
if bond.GetIdx() < heavy.GetNumBonds() and bond.GetIdx() not in bonds:
bond.SetStereo(Chem.BondStereo.STEREONONE)
return bool(Chem.MolToSmiles(Chem.RemoveHs(read)) == smiles)
def _open_stereo(mol: Any) -> int:
from rdkit import Chem
return sum(
1
for info in Chem.FindPotentialStereo(mol)
if info.specified == Chem.StereoSpecified.Unspecified
)
def _conformer(with_h: Any, conformer_id: int, method: str) -> Conformer:
positions = with_h.GetConformer(conformer_id).GetPositions()
return Conformer(
method=method,
symbols=[atom.GetSymbol() for atom in with_h.GetAtoms()],
xyz_angstrom=[[round(float(v), 6) for v in row] for row in positions],
)
[docs]
def embed_3d(smiles: str, *, count: int = 1, attempts: int = 10) -> Embedded:
"""3D conformers of the species ``smiles`` describes, best first.
``attempts`` embeddings are made and refined; the ``count`` lowest in
energy that keep the stereochemistry are returned. Raises
:class:`EmbedFailed` with the reason when there is none.
"""
from rdkit import Chem, rdBase
from rdkit.Chem import rdDistGeom
from rdkit.Chem import rdForceFieldHelpers as force_fields
from vibesketch.chem._rdkit import untyped
from vibesketch.chem.species import mol_from_smiles
if count < 1 or attempts < count:
raise ValueError("ask for at least one conformer, and no more than the attempts")
heavy = mol_from_smiles(smiles)
with_h = Chem.AddHs(heavy)
open_stereo = _open_stereo(heavy)
warnings: list[str] = []
if open_stereo:
warnings.append(
f"{smiles}: the drawing leaves {open_stereo} stereocentre"
f"{'s' if open_stereo > 1 else ''} open; the 3D structure picks one"
)
if with_h.GetNumAtoms() == 1:
atom = Conformer(method="atom", symbols=[with_h.GetAtomWithIdx(0).GetSymbol()],
xyz_angstrom=[[0.0, 0.0, 0.0]]) # fmt: skip
return Embedded([atom], [None], open_stereo, warnings)
distance_geometry = untyped(rdDistGeom)
fields = untyped(force_fields)
params = distance_geometry.ETKDGv3()
params.randomSeed = SEED
params.numThreads = 1
params.pruneRmsThresh = 0.5
blocker = rdBase.BlockLogs()
try:
ids = list(distance_geometry.EmbedMultipleConfs(with_h, attempts, params))
if not ids: # hard cases embed from random coordinates
params.useRandomCoords = True
ids = list(distance_geometry.EmbedMultipleConfs(with_h, attempts, params))
if not ids:
raise EmbedFailed(f"{smiles}: distance geometry found no 3D structure")
energies: list[float | None] = [None] * len(ids)
termination: list[int | None] = [None] * len(ids)
if fields.MMFFHasAllMoleculeParams(with_h):
method = "etkdg+mmff94"
results = fields.MMFFOptimizeMoleculeConfs(with_h, maxIters=2000, numThreads=1)
energies = [float(energy) for _, energy in results]
termination = [int(code) for code, _ in results]
elif fields.UFFHasAllMoleculeParams(with_h):
method = "etkdg+uff"
results = fields.UFFOptimizeMoleculeConfs(with_h, maxIters=2000, numThreads=1)
energies = [float(energy) for _, energy in results]
termination = [int(code) for code, _ in results]
else:
method = "etkdg"
warnings.append(
f"{smiles}: no force field covers it, so the geometry is distance geometry's "
"alone; check metal centres, and optimize before trusting it"
)
for conformer_id, code in zip(ids, termination, strict=True):
if code not in (None, 0):
warnings.append(
f"{smiles}: conformer {conformer_id} {method} termination code {code}; "
"force-field convergence was not achieved"
)
target = _target(heavy)
kept = [
(energy, conformer_id)
for conformer_id, energy in zip(ids, energies, strict=True)
if _stereo_kept(with_h, conformer_id, heavy, target)
]
finally:
del blocker
if not kept:
raise EmbedFailed(
f"{smiles}: no 3D structure kept the drawn stereochemistry "
f"({len(ids)} tried); check the wedges, or draw the isomer explicitly"
)
kept.sort(key=lambda pair: (pair[0] is None, pair[0] or 0.0))
chosen = kept[:count]
records = [
{"conformer_id": int(cid), "termination_code": code, "energy_kcal_mol": energy}
for cid, code, energy in zip(ids, termination, energies, strict=True)
]
conformers = []
for rank, (_, conformer_id) in enumerate(chosen):
record = next(row for row in records if row["conformer_id"] == conformer_id)
preparation = {
"kind": "preparation",
"algorithm": "ETKDGv3",
"seed": SEED,
"requested_attempts": attempts,
"embedded_count": len(ids),
"threads": 1,
"prune_rms_threshold_angstrom": 0.5,
"random_coordinate_fallback": bool(params.useRandomCoords),
"force_field": method.removeprefix("etkdg+") if method != "etkdg" else None,
"max_iterations": 2000,
"selected_conformer_id": int(conformer_id),
"selected_rank": rank,
"termination_code": record["termination_code"],
"force_field_converged": record["termination_code"] == 0,
"energy_kcal_mol": record["energy_kcal_mol"],
"attempts": records,
"warnings": list(warnings),
"verified_quantum_minimum": False,
"input_smiles": smiles,
"atom_order": "input SMILES atoms followed by RDKit explicit hydrogens",
"atom_map_numbers": [atom.GetAtomMapNum() for atom in with_h.GetAtoms()],
}
conformers.append(
_conformer(with_h, conformer_id, method).model_copy(
update={
"source": ConformerSource(
program="RDKit", version=rdBase.rdkitVersion, relaxation=preparation
)
}
)
)
return Embedded(
conformers,
[energy for energy, _ in chosen],
open_stereo,
warnings,
)