"""Molecules from 3D coordinates: bonds, bond orders, charges and radicals.
Programs of the vibe family hand structures over as atoms and Cartesian
coordinates: a QVF ``structure`` section, an XYZ file. The sketcher needs a
molecule to draw, so the bonds are perceived, in this order:
1. **A single atom** keeps the stated charge, and its unpaired electrons
from its valence: N is a quartet, Cl a doublet, Cl⁻ a singlet.
2. **A template.** When the file names the species by SMILES, the bonds come
from the coordinates and the bond orders and charges from the SMILES
(RDKit's ``AssignBondOrdersFromTemplate``). The atoms stay in the file's
order, so the coordinates are kept.
3. **Perception** with the stated charge: RDKit's ``DetermineBonds``, the
xyz2mol algorithm of Kim and Kim, *Bull. Korean Chem. Soc.* 36, 1769
(2015), with Hückel connectivity off.
4. **Connectivity alone,** for a neutral structure the bond orders of
which cannot be assigned, such as a radical or a metal hydride: single
bonds, the unpaired electrons where valences are left open, and two
neighbours' unpaired electrons paired into a further bond (N=O). The
caller says so, because the drawing may need checking.
"""
from __future__ import annotations
import math
from dataclasses import dataclass
from typing import Any, Literal
__all__ = ["Perceived", "perceive"]
How = Literal["atom", "template", "perceived", "connectivity"]
[docs]
@dataclass(frozen=True)
class Perceived:
"""A sanitized RDKit molecule, hydrogens explicit, in the coordinates' atom order."""
mol: Any
#: Which of the module's steps gave the bonds.
how: How
def _from_template(mol: Any, template: Any) -> Any | None:
"""Bond orders and charges from ``template`` onto the coordinates' bonds.
The match compares unpaired electrons, which the bare coordinates do not
have yet, so the template's are cleared for it and the result's assigned
from its open valences afterwards.
"""
from rdkit import Chem
from rdkit.Chem import AllChem, rdDetermineBonds, rdmolops
from vibesketch.chem._rdkit import untyped
connected = Chem.Mol(mol)
reference = Chem.AddHs(template)
for atom in reference.GetAtoms():
atom.SetNumRadicalElectrons(0)
try:
untyped(rdDetermineBonds).DetermineConnectivity(connected)
for atom in connected.GetAtoms():
atom.SetNoImplicit(True)
matched = untyped(AllChem).AssignBondOrdersFromTemplate(reference, connected)
for atom in matched.GetAtoms():
atom.SetNoImplicit(True)
Chem.SanitizeMol(matched)
untyped(rdmolops).AssignRadicals(matched)
except (ValueError, RuntimeError):
return None
return matched
def _paired(mol: Any) -> None:
"""Raise the order of each bond between two atoms that both have an unpaired electron.
Single bonds leave nitric oxide with three unpaired electrons; pairing
two of them across the bond gives N=O with one, as it is drawn.
"""
from rdkit import Chem
raised = {
Chem.BondType.SINGLE: Chem.BondType.DOUBLE,
Chem.BondType.DOUBLE: Chem.BondType.TRIPLE,
}
changed = True
while changed:
changed = False
for bond in mol.GetBonds():
ends = (bond.GetBeginAtom(), bond.GetEndAtom())
if bond.GetBondType() in raised and all(a.GetNumRadicalElectrons() for a in ends):
bond.SetBondType(raised[bond.GetBondType()])
for atom in ends:
atom.SetNumRadicalElectrons(atom.GetNumRadicalElectrons() - 1)
changed = True
def _connectivity(mol: Any) -> Any | None:
from rdkit import Chem
from rdkit.Chem import rdDetermineBonds, rdmolops
from vibesketch.chem._rdkit import untyped
connected = Chem.Mol(mol)
try:
untyped(rdDetermineBonds).DetermineConnectivity(connected)
for atom in connected.GetAtoms():
atom.SetNoImplicit(True)
Chem.SanitizeMol(connected)
untyped(rdmolops).AssignRadicals(connected)
_paired(connected)
Chem.SanitizeMol(connected)
except (ValueError, RuntimeError):
return None
return connected
[docs]
def perceive(
symbols: list[str],
positions: list[list[float]],
*,
charge: int = 0,
template: Any = None,
) -> Perceived | None:
"""The molecule the coordinates show, or ``None`` when no step gives one.
``template`` is an RDKit molecule for the species, from the SMILES a file
states; it is used only when it fits the atoms.
"""
from rdkit import Chem, rdBase
from rdkit.Chem import rdDetermineBonds, rdmolops
from vibesketch.chem._rdkit import untyped
blocker = rdBase.BlockLogs()
try:
if not symbols:
return None
# RDKit's XYZ text parser rejects scientific notation. Build the
# conformer directly so near-zero coordinates keep their precision.
try:
atoms = Chem.RWMol()
conformer = Chem.Conformer(len(symbols))
conformer.Set3D(True)
for index, (symbol, position) in enumerate(zip(symbols, positions, strict=True)):
x, y, z = position
if not all(math.isfinite(value) for value in (x, y, z)):
return None
atoms.AddAtom(Chem.Atom(symbol))
conformer.SetAtomPosition(index, (x, y, z))
mol = atoms.GetMol()
mol.AddConformer(conformer, assignId=True)
except (ValueError, RuntimeError):
return None
if mol.GetNumAtoms() == 1:
atom = mol.GetAtomWithIdx(0)
atom.SetFormalCharge(charge)
atom.SetNoImplicit(True)
try:
Chem.SanitizeMol(mol)
except (ValueError, RuntimeError):
return None
untyped(rdmolops).AssignRadicals(mol)
return Perceived(mol, "atom")
if template is not None:
matched = _from_template(mol, template)
if matched is not None:
return Perceived(matched, "template")
perceived = Chem.Mol(mol)
try:
untyped(rdDetermineBonds).DetermineBonds(perceived, charge=charge)
Chem.SanitizeMol(perceived)
except (ValueError, RuntimeError):
pass
else:
return Perceived(perceived, "perceived")
if charge != 0: # single bonds cannot place a charge
return None
connected = _connectivity(mol)
return None if connected is None else Perceived(connected, "connectivity")
finally:
del blocker