"""Calculated NMR shieldings and shifts of a run's result (#97).
ORCA prints the isotropic shielding of every nucleus in a summary table
(ORCA 6.1 manual, "NMR chemical shifts", whose ethanol example the tests
follow):
https://www.faccts.de/docs/orca/6.1/manual/contents/spectroscopyproperties/nmr.html
The strict reader keeps it on the result by atom index in the job's input.
A chemical shift is a difference of shieldings, δ = σ(reference) − σ, with
the reference computed at the same level of theory: the reference run
(TMS for ¹H and ¹³C, CFCl₃ for ¹⁹F, …) is named and recorded, its job
specification must equal the run's, and its atoms of the element must all
be equivalent. No literature constant is ever used in its place. Nuclei
equivalent by the molecule's symmetry (RDKit's canonical ranking without
tie-breaking, on the graph perceived from the run's geometry) are averaged,
as a spectrum shows one line for a methyl group.
Couplings are not read: vibe-sketch's ORCA writer requests shieldings only.
"""
from __future__ import annotations
import csv
import io
import re
from collections.abc import Sequence
from typing import Any
from vibesketch.document.model import (
Canvas,
Document,
Job,
ResultNucleus,
ResultShieldings,
Run,
)
from vibesketch.document.tables import MEDIA_TYPES, TableFormat, _latex, _markdown
__all__ = [
"MEDIA_TYPES",
"NmrRefused",
"nmr_shifts",
"nmr_table",
"orca_shieldings",
"with_nmr_labels",
]
_NUMBER = r"[+-]?(?:\d+(?:\.\d*)?|\.\d+)(?:[EeDd][+-]?\d+)?"
_HEADER = re.compile(r"^[ \t]*Nucleus[ \t]+Element[ \t]+Isotropic[ \t]+Anisotropy[ \t]*$", re.M)
_ROW = re.compile(rf"^[ \t]*(\d+)[ \t]+([A-Z][a-z]?)[ \t]+({_NUMBER})[ \t]+({_NUMBER})[ \t]*$")
[docs]
class NmrRefused(ValueError):
"""Shieldings that cannot be read, or shifts that cannot be worked out as asked."""
[docs]
def orca_shieldings(text: str, symbols: Sequence[str]) -> ResultShieldings | None:
"""The shielding summary of an ORCA 6 output, checked against the input's atoms."""
headers = list(_HEADER.finditer(text))
if not headers:
return None
if len(headers) > 1:
raise NmrRefused("the ORCA output prints its shielding summary more than once")
lines = text[headers[0].end() :].lstrip("\n").splitlines()
if not lines or not re.fullmatch(r"[ \t]*-+[ \t]+-+[ \t]+-+[ \t]+-+[ \t]*", lines[0]):
raise NmrRefused("the ORCA shielding summary has no column rule under its header")
nuclei: list[ResultNucleus] = []
for line in lines[1:]:
if not line.strip():
break
match = _ROW.match(line)
if match is None:
raise NmrRefused(f"unreadable ORCA shielding row: {line.strip()!r}")
nuclei.append(
ResultNucleus(
atom=int(match[1]),
element=match[2],
isotropic=float(match[3].replace("D", "E")),
anisotropy=float(match[4].replace("D", "E")),
)
)
if [n.atom for n in nuclei] != list(range(len(symbols))):
raise NmrRefused(
f"the ORCA shielding summary lists {len(nuclei)} nuclei, not the input's "
f"{len(symbols)} atoms in order"
)
for nucleus in nuclei:
if nucleus.element != symbols[nucleus.atom]:
raise NmrRefused(
f"ORCA's nucleus {nucleus.atom} is {nucleus.element}, the input's "
f"{symbols[nucleus.atom]}"
)
number = text[: headers[0].start()].count("\n") + 1
return ResultShieldings(source=f"line {number}: shielding summary", nuclei=nuclei)
def _run(document: Document, identifier: str) -> tuple[Job, Run, ResultShieldings]:
for job in document.jobs:
for run in job.runs or []:
if run.id == identifier:
if run.result is None or run.result.shieldings is None:
raise NmrRefused(f"run {identifier!r} has no imported NMR shieldings")
return job, run, run.result.shieldings
raise NmrRefused(f"no run {identifier!r} in this document")
def _perceived(document: Document, job: Job) -> Any:
"""The run's molecule in its input's atom order, hydrogens included."""
from vibesketch.chem.perceive import perceive
from vibesketch.chem.species import mol_from_smiles
species = document.species_by_id(job.species)
if not species.conformers:
raise NmrRefused(f"species {species.id} has no 3D structure to compare atoms with")
conformer = species.conformers[0]
template = mol_from_smiles(species.smiles) if species.smiles else None
found = perceive(
conformer.symbols, conformer.xyz_angstrom, charge=species.charge, template=template
)
if found is None:
raise NmrRefused(f"the bonds of species {species.id} cannot be perceived")
return found.mol
def _classes(mol: Any) -> list[int]:
"""Each atom's symmetry class: equal for atoms the graph cannot tell apart."""
from rdkit import Chem
return [int(rank) for rank in Chem.CanonicalRankAtoms(mol, breakTies=False)]
def _reference(document: Document, job: Job, reference: str) -> dict[str, float]:
"""The reference shielding per element, from a run at the same level of theory."""
ref_job, _, shieldings = _run(document, reference)
if ref_job.spec != job.spec or ref_job.program != job.program:
raise NmrRefused(
f"the reference run {reference!r} was not computed at the run's level of theory"
)
classes = _classes(_perceived(document, ref_job))
by_element: dict[str, list[ResultNucleus]] = {}
for nucleus in shieldings.nuclei:
by_element.setdefault(nucleus.element, []).append(nucleus)
values: dict[str, float] = {}
for element, nuclei in by_element.items():
if len({classes[n.atom] for n in nuclei}) == 1:
values[element] = sum(n.isotropic for n in nuclei) / len(nuclei)
return values
[docs]
def nmr_shifts(
document: Document,
run: str,
*,
reference: str | None = None,
nuclei: Sequence[str] | None = None,
) -> dict[str, Any]:
"""One row per set of equivalent nuclei: its atoms, shielding and, with a reference, shift."""
job, _, shieldings = _run(document, run)
if reference == run:
raise NmrRefused("a run cannot be its own reference")
references = _reference(document, job, reference) if reference else {}
mol = _perceived(document, job)
classes = _classes(mol)
groups: dict[int, list[ResultNucleus]] = {}
for nucleus in shieldings.nuclei:
if nuclei and nucleus.element not in nuclei:
continue
groups.setdefault(classes[nucleus.atom], []).append(nucleus)
rows: list[dict[str, Any]] = []
for members in sorted(groups.values(), key=lambda m: m[0].atom):
element = members[0].element
shielding = sum(n.isotropic for n in members) / len(members)
shift = references[element] - shielding if element in references else None
attached = [
f"{a.GetSymbol()}{a.GetIdx() + 1}"
for a in mol.GetAtomWithIdx(members[0].atom).GetNeighbors()
if a.GetAtomicNum() != 1
]
rows.append(
{
"element": element,
"atoms": [n.atom for n in members],
"labels": [f"{element}{n.atom + 1}" for n in members],
"on": attached[0] if element == "H" and attached else None,
"shielding_ppm": shielding,
"shift_ppm": shift,
}
)
elements = {str(r["element"]) for r in rows}
missing = sorted(elements - set(references)) if reference else []
return {
"run": run,
"reference": reference,
"reference_shielding_ppm": references,
"without_reference": missing,
"source": shieldings.source,
"rows": rows,
}
[docs]
def nmr_table(
document: Document,
run: str,
*,
reference: str | None = None,
nuclei: Sequence[str] | None = None,
format: TableFormat = "csv",
) -> bytes:
"""The per-nucleus line list as CSV (full precision), Markdown or a LaTeX fragment."""
if format not in MEDIA_TYPES:
raise NmrRefused("table format must be csv, markdown or latex")
found = nmr_shifts(document, run, reference=reference, nuclei=nuclei)
if format == "csv":
stream = io.StringIO(newline="")
writer = csv.writer(stream, lineterminator="\n")
writer.writerow(["element", "atoms", "count", "on", "shielding_ppm", "shift_ppm"])
for row in found["rows"]:
writer.writerow(
[
row["element"],
" ".join(row["labels"]),
len(row["atoms"]),
row["on"] or "",
row["shielding_ppm"],
"" if row["shift_ppm"] is None else row["shift_ppm"],
]
)
return stream.getvalue().encode("utf-8")
headers = ["Nucleus", "Atoms", "n", "σ (ppm)", "δ (ppm)"]
values = [headers]
for row in found["rows"]:
atoms = ", ".join(row["labels"]) + (f" (on {row['on']})" if row["on"] else "")
values.append(
[
row["element"],
atoms,
str(len(row["atoms"])),
f"{row['shielding_ppm']:.2f}",
"" if row["shift_ppm"] is None else f"{row['shift_ppm']:.2f}",
]
)
caption = (
f"Run {run}; shifts against run {reference}, computed at the same level."
if reference
else f"Run {run}; absolute shieldings, no reference chosen."
)
if format == "markdown":
lines = ["| " + " | ".join(_markdown(c) for c in line) + " |" for line in values]
lines.insert(1, "| " + " | ".join(["---"] * len(headers)) + " |")
return ("\n".join(lines) + "\n\n" + caption + "\n").encode("utf-8")
lines = [r"% Requires \usepackage{booktabs}.", r"\begin{tabular}{@{}llrrr@{}}", r"\toprule"]
for index, line in enumerate(values):
lines.append(" & ".join(_latex(c) for c in line) + r" \\")
if index == 0:
lines.append(r"\midrule")
lines += [r"\bottomrule", r"\end{tabular}", r"\par\smallskip " + _latex(caption), ""]
return "\n".join(lines).encode("utf-8")
[docs]
def with_nmr_labels(
document: Document, run: str, *, reference: str, nuclei: Sequence[str] = ("H", "C")
) -> Document:
"""``document`` with shift labels beside the run species' first drawn structure.
Each heavy atom gets one text: its own shift (δC 18.1) and those of the
hydrogens on it (δH 1.17). The run's atoms are matched to the drawing by
graph; equivalent atoms have one shift, so any match gives the same labels.
"""
from rdkit import Chem
from vibesketch.chem.species import identity
from vibesketch.document.derive import fragments
from vibesketch.document.ket import text_node
job, _, _ = _run(document, run)
found = nmr_shifts(document, run, reference=reference, nuclei=nuclei)
shift = {atom: row for row in found["rows"] for atom in row["atoms"]}
species = document.species_by_id(job.species)
drawn = next(
(
f
for f in fragments(document)
if f.mol is not None and identity(f.mol).smiles == species.smiles
),
None,
)
if drawn is None:
raise NmrRefused(f"species {species.id} is not drawn on the canvas")
mol = _perceived(document, job)
heavy = Chem.RemoveHs(mol) # the input's heavy atoms come first, in order
picture = Chem.Mol(drawn.mol)
for atom in picture.GetAtoms():
atom.SetIntProp("drawn", atom.GetIdx())
picture = Chem.RemoveHs(picture)
match = picture.GetSubstructMatch(heavy)
if len(match) != heavy.GetNumAtoms():
raise NmrRefused(f"the run's atoms do not match the drawing of species {species.id}")
data: dict[str, Any] = {**document.canvas.data}
root: dict[str, Any] = {**data.get("root", {})}
nodes: list[Any] = list(root.get("nodes", []))
for index in range(heavy.GetNumAtoms()):
lines = []
row = shift.get(index)
if row is not None and row["shift_ppm"] is not None:
lines.append(f"δ{row['element']} {row['shift_ppm']:.1f}")
hydrogens = sorted(
{
round(shift[h.GetIdx()]["shift_ppm"], 2)
for h in mol.GetAtomWithIdx(index).GetNeighbors()
if h.GetAtomicNum() == 1
and h.GetIdx() in shift
and shift[h.GetIdx()]["shift_ppm"] is not None
}
)
if hydrogens:
lines.append("δH " + "/".join(f"{v:.2f}" for v in hydrogens))
if not lines:
continue
drawn_atom = picture.GetAtomWithIdx(match[index])
position = drawn.mol.GetConformer().GetAtomPosition(drawn_atom.GetIntProp("drawn"))
neighbours = [
drawn.mol.GetConformer().GetAtomPosition(n.GetIdx())
for n in drawn.mol.GetAtomWithIdx(drawn_atom.GetIntProp("drawn")).GetNeighbors()
]
dx = position.x - sum(p.x for p in neighbours) / len(neighbours) if neighbours else 0.0
dy = position.y - sum(p.y for p in neighbours) / len(neighbours) if neighbours else 1.0
length = (dx * dx + dy * dy) ** 0.5 or 1.0
away = 0.45 + 0.2 * len(lines)
centre = (position.x + dx / length * away, position.y + dy / length * away)
nodes.append(text_node("\n".join(lines), centre))
root["nodes"] = nodes
data["root"] = root
return document.model_copy(update={"canvas": Canvas(format="ket", data=data)})