"""Elemental analysis and calculated descriptors of a species (#93).
The elemental analysis is exact arithmetic on the formula: each element's
share of the molar mass, with RDKit's standard atomic weights and an isotope
label's own mass. The analysis line is in the ACS style for combustion
(CHNS) analysis, "Anal. Calcd for C9H8O4: C, 60.00; H, 4.48.", ACS Guide to
Scholarly Communication, "Elemental analysis":
https://pubs.acs.org/doi/book/10.1021/acsguide
The descriptors are RDKit's, and they are estimates, not measurements:
- logP and molar refractivity: atom contributions of S. A. Wildman and
G. M. Crippen, J. Chem. Inf. Comput. Sci. 39 (1999) 868-873,
doi:10.1021/ci990307l (``Crippen.MolLogP``, ``Crippen.MolMR``);
- topological polar surface area: P. Ertl, B. Rohde and P. Selzer,
J. Med. Chem. 43 (2000) 3714-3717, doi:10.1021/jm000942e
(``rdMolDescriptors.CalcTPSA``, N and O only);
- hydrogen-bond donors and acceptors in Lipinski's counting, C. A. Lipinski
et al., Adv. Drug Delivery Rev. 23 (1997) 3-25,
doi:10.1016/S0169-409X(96)00423-1 (``CalcNumLipinskiHBD``/``HBA``);
- rotatable bonds, heavy atoms and rings (SSSR) as RDKit counts them.
"""
from __future__ import annotations
from dataclasses import asdict, dataclass
from typing import Any
__all__ = ["Analysis", "Descriptors", "descriptors", "elemental_analysis"]
#: The elements a combustion analysis reports, in the order the line gives them.
COMBUSTION = ("C", "H", "N", "S")
[docs]
@dataclass(frozen=True)
class Analysis:
"""Mass percent per element in Hill order, and the formatted analysis line."""
formula: str
mass_percent: dict[str, float]
line: str
[docs]
def to_dict(self) -> dict[str, Any]:
return asdict(self)
[docs]
@dataclass(frozen=True)
class Descriptors:
"""Calculated estimates: RDKit's values, unrounded."""
logp: float
molar_refractivity: float
tpsa: float
hbond_donors: int
hbond_acceptors: int
rotatable_bonds: int
heavy_atoms: int
rings: int
[docs]
def to_dict(self) -> dict[str, Any]:
return asdict(self)
def _hill(symbols: list[str]) -> list[str]:
"""Hill order: carbon, then hydrogen, then the rest alphabetically; with no carbon, all."""
present = set(symbols)
if "C" in present:
return ["C", *(["H"] if "H" in present else []), *sorted(present - {"C", "H"})]
return sorted(present)
[docs]
def elemental_analysis(mol: Any) -> Analysis:
"""The mass percent of each element of ``mol``, implicit hydrogens included."""
from rdkit import Chem
from rdkit.Chem import rdMolDescriptors
from vibesketch.chem.species import without_hydrogens
table = Chem.GetPeriodicTable()
hydrogen = float(table.GetAtomicWeight(1))
masses: dict[str, float] = {}
for atom in mol.GetAtoms():
symbol = str(atom.GetSymbol())
masses[symbol] = masses.get(symbol, 0.0) + float(atom.GetMass())
hydrogens = int(atom.GetTotalNumHs())
if hydrogens:
masses["H"] = masses.get("H", 0.0) + hydrogens * hydrogen
total = sum(masses.values())
percent = {symbol: 100 * masses[symbol] / total for symbol in _hill(list(masses))}
# An isotope label shows in the formula, as its mass does in the percentages.
formula = str(rdMolDescriptors.CalcMolFormula(without_hydrogens(mol), separateIsotopes=True))
found = [f"{symbol}, {percent[symbol]:.2f}" for symbol in COMBUSTION if symbol in percent]
line = f"Anal. Calcd for {formula}: {'; '.join(found)}." if found else ""
return Analysis(formula=formula, mass_percent=percent, line=line)
[docs]
def descriptors(mol: Any) -> Descriptors:
"""Crippen logP and molar refractivity, TPSA, Lipinski donors and acceptors, counts."""
from rdkit.Chem import Crippen, rdMolDescriptors
from vibesketch.chem._rdkit import untyped
from vibesketch.chem.species import without_hydrogens
heavy = without_hydrogens(mol)
crippen = untyped(Crippen)
return Descriptors(
logp=float(crippen.MolLogP(heavy)),
molar_refractivity=float(crippen.MolMR(heavy)),
tpsa=float(rdMolDescriptors.CalcTPSA(heavy)),
hbond_donors=int(rdMolDescriptors.CalcNumLipinskiHBD(heavy)),
hbond_acceptors=int(rdMolDescriptors.CalcNumLipinskiHBA(heavy)),
rotatable_bonds=int(rdMolDescriptors.CalcNumRotatableBonds(heavy)),
heavy_atoms=int(heavy.GetNumHeavyAtoms()),
rings=int(rdMolDescriptors.CalcNumRings(heavy)),
)