"""Identity and properties of one species, from an RDKit molecule."""
from __future__ import annotations
from dataclasses import dataclass
from typing import Any
__all__ = [
"Identity",
"Properties",
"TRANSITION_METALS",
"identity",
"DRAWN_RADICAL",
"inert_pair_radicals",
"mol_from_smiles",
"properties",
"unpaired",
]
#: d-block and f-block elements: their spin state cannot be read off a drawing.
TRANSITION_METALS = frozenset([*range(21, 31), *range(39, 49), *range(57, 81), *range(89, 113)])
[docs]
@dataclass(frozen=True)
class Identity:
"""What identifies a species and what a calculation needs from the drawing."""
smiles: str
inchikey: str | None
charge: int
#: The multiplicity the drawing implies: unpaired radical electrons + 1.
multiplicity: int
electrons: int
has_transition_metal: bool
has_stereo: bool
[docs]
@property
def key(self) -> str:
"""What two drawings of the same species share."""
return self.inchikey or self.smiles
[docs]
@dataclass(frozen=True)
class Properties:
formula: str
average_mass: float
monoisotopic_mass: float
inchi: str | None
def _quiet() -> Any:
from rdkit import rdBase
return rdBase.BlockLogs()
def without_hydrogens(mol: Any) -> Any:
"""A copy without explicit hydrogens, quietly: a lone H or H+ stays, unremarked."""
from rdkit import Chem
blocker = _quiet() # "not removing hydrogen atom without neighbors"
heavy = Chem.RemoveHs(mol)
del blocker
return heavy
#: p-block elements whose ns² pair stays paired (the inert-pair effect), by their
#: valence electron count: Ga, In, Tl; Ge, Sn, Pb; As, Sb, Bi.
_INERT_PAIR = {31: 3, 49: 3, 81: 3, 32: 4, 50: 4, 82: 4, 33: 5, 51: 5, 83: 5}
#: An atom property marking a radical the drawing or file states, which is kept as given.
DRAWN_RADICAL = "vibesketchDrawnRadical"
[docs]
def unpaired(atom: Any) -> int:
"""An atom's unpaired electrons, with the inert pair of Ga to Bi kept paired (#147).
RDKit's valence model gives a bare Pb²⁺ two radical electrons, a triplet,
and re-assigns them whenever a molecule is sanitized; the ion is
[Xe]4f¹⁴5d¹⁰6s², a singlet. For the elements above, the electrons left
after the charge and the bonds fill the ns² pair first and then the np
shell by Hund's rule: Pb²⁺, Sn²⁺, In⁺ and Bi³⁺ have none unpaired, a Pb
atom two (³P), a Bi atom three (⁴S), a Ga atom one (²P). A radical marked
as drawn (:data:`DRAWN_RADICAL`) is kept as given.
"""
radicals = int(atom.GetNumRadicalElectrons())
valence = _INERT_PAIR.get(int(atom.GetAtomicNum()))
if valence is None or not radicals or atom.HasProp(DRAWN_RADICAL):
return radicals
left = valence - int(atom.GetFormalCharge()) - int(atom.GetTotalValence())
if left < 0:
return radicals
shell = left - 2 if left >= 2 else left
return shell if shell <= 3 else 6 - shell
[docs]
def inert_pair_radicals(mol: Any, drawn: Any = None) -> None:
"""Mark the radicals ``drawn(index)`` says were drawn, and set the others to :func:`unpaired`.
The counts set last only until RDKit sanitizes the molecule again;
:func:`unpaired` gives them at any time.
"""
for atom in mol.GetAtoms():
if drawn is not None and atom.GetNumRadicalElectrons() and drawn(int(atom.GetIdx())):
atom.SetBoolProp(DRAWN_RADICAL, True)
continue
count = unpaired(atom)
if count != atom.GetNumRadicalElectrons():
atom.SetNumRadicalElectrons(count)
[docs]
def mol_from_smiles(smiles: str) -> Any:
"""An RDKit molecule from SMILES; raises ``ValueError`` when RDKit refuses it."""
from rdkit import Chem
blocker = _quiet()
mol = Chem.MolFromSmiles(smiles)
del blocker
if mol is None:
raise ValueError(f"RDKit cannot read the SMILES {smiles!r}")
inert_pair_radicals(mol) # SMILES writes no radical on these; RDKit's guess is replaced
return mol
[docs]
def identity(mol: Any) -> Identity:
"""The identity of a sanitized RDKit molecule."""
from rdkit import Chem as chem_module
from rdkit.Chem import inchi as inchi_module
from vibesketch.chem._rdkit import untyped
Chem = untyped(chem_module)
inchi = untyped(inchi_module)
heavy = without_hydrogens(mol)
smiles = str(Chem.MolToSmiles(heavy))
blocker = _quiet()
try:
key = str(inchi.MolToInchiKey(heavy)) or None
except Exception: # InChI refuses some valid structures, for example some metals
key = None
del blocker
charge = int(Chem.GetFormalCharge(heavy))
radicals = sum(unpaired(atom) for atom in heavy.GetAtoms())
with_h = Chem.AddHs(heavy)
electrons = sum(int(atom.GetAtomicNum()) for atom in with_h.GetAtoms()) - charge
numbers = {int(atom.GetAtomicNum()) for atom in heavy.GetAtoms()}
stereo = Chem.FindMolChiralCenters(
heavy, includeUnassigned=False, useLegacyImplementation=False
)
# E/Z or cis/trans: RDKit's newer stereo perception reports the latter.
double = any(
bond.GetStereo() not in (Chem.BondStereo.STEREONONE, Chem.BondStereo.STEREOANY)
for bond in heavy.GetBonds()
)
return Identity(
smiles=smiles,
inchikey=key,
charge=charge,
multiplicity=radicals + 1,
electrons=electrons,
has_transition_metal=bool(numbers & TRANSITION_METALS),
has_stereo=bool(stereo) or double,
)
[docs]
def properties(mol: Any) -> Properties:
"""Formula, masses and InChI of a sanitized RDKit molecule."""
from rdkit.Chem import Descriptors as descriptors_module
from rdkit.Chem import inchi as inchi_module
from rdkit.Chem import rdMolDescriptors
from vibesketch.chem._rdkit import untyped
Descriptors = untyped(descriptors_module)
inchi = untyped(inchi_module)
heavy = without_hydrogens(mol)
blocker = _quiet()
try:
text = str(inchi.MolToInchi(heavy)) or None
except Exception:
text = None
del blocker
return Properties(
formula=str(rdMolDescriptors.CalcMolFormula(heavy)),
average_mass=float(Descriptors.MolWt(heavy)),
monoisotopic_mass=float(Descriptors.ExactMolWt(heavy)),
inchi=text,
)