Source code for vibesketch.document.results

"""Explicit, provenance-bound result imports and reaction energies (#28).

An optimisation's converged geometry comes back with its result, as a new
conformer of the species named ``orca:<task>:<first 12 hex digits of the
source's SHA-256>``, so the conformer and the result it came from name each
other. An ORCA Hessian file (``.hess``) is checked against the run's stationary
point (#85).

ORCA text vocabulary: ORCA 6.1 manual, Hello water and Thermochemistry:
https://www.faccts.de/docs/orca/6.1/manual/contents/quickstartguide/hellowater.html
https://www.faccts.de/docs/orca/6.1/manual/contents/structurereactivity/thermochemistry.html
QVF root energies and thermochemistry: format specification sections 3.3/3.6.
No filesystem discovery, implicit run selection or companion-dependent parsing.
"""

from __future__ import annotations

import hashlib
import io
import json
import math
import re
import zipfile
from collections import Counter
from datetime import UTC, datetime
from typing import Any, Literal

from pydantic import BaseModel, ConfigDict, Field

from vibesketch.document.model import (
    Conformer,
    Document,
    Job,
    ResultEnergy,
    Run,
    RunResult,
    Species,
)
from vibesketch.units import BOHR_TO_ANGSTROM, HARTREE_TO_KCAL_MOL

MAX_RESULT_BYTES = 32 * 1024 * 1024
_NUMBER = r"[+-]?(?:\d+(?:\.\d*)?|\.\d+)(?:[EeDd][+-]?\d+)?"


[docs] class ResultRefused(ValueError): """The supplied evidence cannot support this result or comparison."""
[docs] class Selection(BaseModel): """Select a group or exactly one run per species, never an implicit latest run. Coefficients, when supplied, name every species on both chemical sides of the requested reaction. Unspecified reactions use one of each species; atom/charge balance is checked, so an omitted factor cannot pass silently. """ model_config = ConfigDict(extra="forbid", strict=True) group: str | None = None runs: list[str] | None = None reactions: list[str] | None = None coefficients: dict[str, dict[str, dict[str, float]]] = Field(default_factory=dict)
def _digest(value: Any) -> str: return hashlib.sha256(json.dumps(value, sort_keys=True, allow_nan=False).encode()).hexdigest()
[docs] def species_digest(species: Species) -> str: """Bind identity, charge/spin and the exact input geometry; labels are editable.""" return _digest( { "smiles": species.smiles, "inchikey": species.inchikey, "charge": species.charge, "multiplicity": species.multiplicity, "conformer": species.conformers[0].model_dump(mode="json") if species.conformers else None, } )
def _run(document: Document, identifier: str) -> tuple[Job, Run]: matches = [ (job, run) for job in document.jobs for run in job.runs or [] if run.id == identifier ] if len(matches) != 1: raise ResultRefused(f"run {identifier!r} must identify exactly one recorded submission") return matches[0] def _finite(value: Any, label: str) -> float: if isinstance(value, bool) or not isinstance(value, int | float) or not math.isfinite(value): raise ResultRefused(f"{label} must be a finite number") return float(value) def _quantity(value: Any, label: str) -> ResultEnergy: if not isinstance(value, dict) or str(value.get("units", "")).lower() not in { "hartree", "eh", "ha", }: raise ResultRefused(f"{label} must state hartree units") return ResultEnergy(value=_finite(value.get("value"), label), source=label) #: How far, in angstrom, two printings of one geometry may differ: ORCA prints #: six decimals in angstrom in its output and twelve in bohr in a Hessian file. _SAME_GEOMETRY = 1e-4 _Geometry = tuple[list[str], list[list[float]]] def _coordinate_blocks(text: str) -> list[tuple[int, _Geometry]]: """Every ``CARTESIAN COORDINATES (ANGSTROEM)`` block, with where it starts.""" blocks = [] atom = re.compile(rf"^\s*([A-Z][a-z]?)\s+({_NUMBER})\s+({_NUMBER})\s+({_NUMBER})\s*$") for header in re.finditer(r"^-+\n\s*CARTESIAN COORDINATES \(ANGSTROEM\)\s*\n-+\n", text, re.M): symbols: list[str] = [] xyz: list[list[float]] = [] for line in text[header.end() :].splitlines(): found = atom.match(line) if found is None: break symbols.append(found[1]) xyz.append([float(found[i].replace("D", "E")) for i in (2, 3, 4)]) if symbols: blocks.append((header.start(), (symbols, xyz))) return blocks def _same(first: list[list[float]], second: list[list[float]]) -> bool: return len(first) == len(second) and all( math.dist(a, b) <= _SAME_GEOMETRY for a, b in zip(first, second, strict=True) ) def _same_up_to_translation(first: list[list[float]], second: list[list[float]]) -> bool: """Compare atom-ordered geometries without requiring the same coordinate origin. Real ORCA 6 opt+freq Hessians can translate the output's stationary point. Remove only the centroid: do not rotate, reflect, scale or remap atoms. The existing per-atom tolerance still applies to the remaining differences. """ if not first or len(first) != len(second): return False def centered(points: list[list[float]]) -> list[list[float]]: center = [sum(point[axis] / len(points) for point in points) for axis in range(3)] return [[value - center[axis] for axis, value in enumerate(point)] for point in points] return _same(centered(first), centered(second)) def _checked_atoms(species: Species, symbols: list[str], what: str) -> None: """The atoms are the species' input atoms, in their order: nothing is re-mapped.""" if not species.conformers or symbols != species.conformers[0].symbols: raise ResultRefused(f"{what} does not list this species' atoms in its input order") def _final_geometry(text: str, species: Species) -> _Geometry: """The geometry ORCA printed after the optimisation converged; every printing must agree.""" converged = text.index("THE OPTIMIZATION HAS CONVERGED") end = text.index("ORCA TERMINATED NORMALLY") after = [geometry for start, geometry in _coordinate_blocks(text) if converged < start < end] if not after: raise ResultRefused("the ORCA output prints no geometry after the optimisation converged") symbols, xyz = after[-1] if any(s != symbols or not _same(x, xyz) for s, x in after): raise ResultRefused("the ORCA output prints different geometries after convergence") _checked_atoms(species, symbols, "the optimised geometry") return symbols, xyz def _hessian_geometry(data: bytes) -> _Geometry: """The ``$atoms`` block of an ORCA Hessian file: symbol, mass, x, y, z in bohr.""" try: text = data.decode("utf-8") except UnicodeDecodeError: raise ResultRefused("an ORCA Hessian file must be UTF-8 text") from None if not text.lstrip().startswith("$orca_hessian_file"): raise ResultRefused("not an ORCA Hessian file: it does not start with $orca_hessian_file") blocks = re.findall(r"^\$atoms\s*\n\s*(\d+)\s*\n((?:[^$\n]*\n?)*)", text, re.M) if len(blocks) != 1: raise ResultRefused("expected one $atoms block in the Hessian file") count, body = int(blocks[0][0]), blocks[0][1].splitlines() line = re.compile( rf"^\s*([A-Z][a-z]?)\s+{_NUMBER}\s+({_NUMBER})\s+({_NUMBER})\s+({_NUMBER})\s*$" ) rows = [line.match(row) for row in body[:count]] if count == 0 or len(rows) != count or any(row is None for row in rows): raise ResultRefused("the Hessian file's $atoms block does not list its atoms") symbols = [row[1] for row in rows if row is not None] xyz = [ [float(row[i].replace("D", "E")) * BOHR_TO_ANGSTROM for i in (2, 3, 4)] for row in rows if row is not None ] if not all(math.isfinite(v) for point in xyz for v in point): raise ResultRefused("the Hessian file's geometry is not finite") return symbols, xyz def _method(task: str, digest: str) -> str: """The conformer method naming the result a geometry came from.""" return f"orca:{task}:{digest[:12]}" def _with_hessian(document: Document, run_id: str, data: bytes) -> Document: """``document`` after checking an ORCA Hessian file against its run's stationary point. A frequency job runs at its input geometry, an ``opt+freq`` job at the optimised geometry attached with its output; the Hessian file must hold that geometry. It is attached only where the species does not already hold it. """ job, run = _run(document, run_id) if job.program != "orca" or job.spec.task not in {"freq", "opt+freq"}: raise ResultRefused("a Hessian file belongs to an ORCA freq or opt+freq run") if run.result is None: raise ResultRefused("import the run's ORCA output first; the Hessian file is checked " "against the geometry it verified") # fmt: skip species = document.species_by_id(job.species) symbols, xyz = _hessian_geometry(data) _checked_atoms(species, symbols, "the Hessian file") if job.spec.task == "freq": expected: Conformer | None = species.conformers[0] else: method = _method(job.spec.task, run.result.source_sha256) expected = next((c for c in species.conformers if c.method == method), None) if expected is None or not _same_up_to_translation(expected.xyz_angstrom, xyz): raise ResultRefused("the Hessian file's geometry is not this run's stationary point") return document def _orca(document: Document, job: Job, data: bytes) -> dict[str, Any]: from vibesketch.inputs import write_input if job.program != "orca" or job.spec.task not in {"sp", "opt", "freq", "opt+freq"}: raise ResultRefused("the ORCA reader supports ORCA sp, opt, freq and opt+freq jobs") try: text = data.decode("utf-8") except UnicodeDecodeError: raise ResultRefused("ORCA output must be UTF-8 text") from None versions = re.findall(r"^\s*Program Version\s+(6\.\d+[^\s]*)", text, re.M) if len(versions) != 1 or text.count("ORCA TERMINATED NORMALLY") != 1: raise ResultRefused("expected one normally terminated ORCA 6 calculation") if re.search( r"SCF NOT CONVERGED|SCF CONVERGENCE FAILURE|ORCA finished by error termination", text, re.I ): raise ResultRefused("ORCA reports a failed calculation") if "SCF CONVERGED AFTER" not in text: raise ResultRefused("ORCA output does not confirm SCF convergence") if "opt" in job.spec.task and "THE OPTIMIZATION HAS CONVERGED" not in text: raise ResultRefused("normal termination alone does not establish optimisation convergence") # Match the echoed input, including coordinates, charge/spin and every option. # Reject compound/multi-job outputs instead of selecting an arbitrary energy. echo = re.findall(r"^[ \t]*\|[ \t]*\d+>[ \t]?(.*)$", text, re.M) if echo and echo[-1].strip() == "****END OF INPUT****": echo.pop() # ORCA appends this numbered marker; it is not an input directive. expected = write_input(document, job=job.id).files[0].content.splitlines() def normalize(lines: list[str]) -> list[str]: return [ " ".join(line.split()) for line in lines if line.strip() and not line.lstrip().startswith("#") ] if normalize(echo) != normalize(expected): raise ResultRefused("ORCA's echoed input does not match this job and input geometry") energies = list(re.finditer(rf"^\s*FINAL SINGLE POINT ENERGY\s+({_NUMBER})\s*$", text, re.M)) if not energies: raise ResultRefused("no supported final electronic energy in the ORCA output") if job.spec.task in {"sp", "freq"} and len(energies) != 1: raise ResultRefused("multiple electronic energies are ambiguous for this task") match = energies[-1] energy = ResultEnergy( value=float(match[1].replace("D", "E")), source=f"line {text[: match.start()].count(chr(10)) + 1}: final electronic energy", ) if text.index("ORCA TERMINATED NORMALLY") < match.end(): raise ResultRefused("energy occurs after the normal-termination marker") result: dict[str, Any] = { "parser": "orca-text/1", "program_version": versions[0], "energy": energy, } if "freq" in job.spec.task: patterns = { "temperature_k": rf"^\s*Temperature\s+\.+\s+({_NUMBER})\s+K\s*$", "pressure_atm": rf"^\s*Pressure\s+\.+\s+({_NUMBER})\s+atm\s*$", "enthalpy": rf"^\s*Total [Ee]nthalpy\s+\.+\s+({_NUMBER})\s+Eh\s*$", "gibbs": rf"^\s*Final Gibbs free energy\s+\.+\s+({_NUMBER})\s+Eh\s*$", } thermal: dict[str, Any] = {} for key, pattern in patterns.items(): matches = list(re.finditer(pattern, text, re.M)) # ORCA repeats total enthalpy in its Gibbs summary. Both values # must agree; temperatures and other quantities remain unique. allowed = {1, 2} if key == "enthalpy" else {1} if len(matches) not in allowed: raise ResultRefused( f"expected one thermochemistry {key}; " "multiple temperatures/blocks are unsupported" ) if len(matches) == 2 and matches[0][1] != matches[1][1]: raise ResultRefused("conflicting thermochemistry enthalpy values") m = matches[0] if any( m.start() < match.end() or m.end() > text.index("ORCA TERMINATED NORMALLY") for m in matches ): raise ResultRefused( "thermochemistry must follow the final energy before termination" ) value = float(m[1].replace("D", "E")) thermal[key] = ( ResultEnergy( value=value, source=f"line {text[: m.start()].count(chr(10)) + 1}: {key}" ) if key in {"enthalpy", "gibbs"} else value ) result["thermochemistry"] = thermal if "opt" in job.spec.task: result["geometry"] = _final_geometry(text, document.species_by_id(job.species)) return result def _qvf(document: Document, job: Job, data: bytes) -> dict[str, Any]: from vibesketch.document.qvf import _atoms, _job_spec, _modules if job.program != "vibeqc": raise ResultRefused("QVF results currently require a recorded vibeqc job") _, library = _modules() limits = library.Limits( max_member_bytes=MAX_RESULT_BYTES, max_total_bytes=2 * MAX_RESULT_BYTES, max_json_bytes=MAX_RESULT_BYTES, ) with library.QvfReader(io.BytesIO(data), limits=limits) as reader: validation = library.validate_qvf(io.BytesIO(data), limits=limits) if not validation["ok"]: raise ResultRefused( "QVF validation failed: " + "; ".join(f["code"] for f in validation["findings"] if f["severity"] == "error") ) manifest = reader.manifest provenance = manifest.get("provenance", {}) if provenance.get("run_status") != "converged": raise ResultRefused("QVF does not report a converged result") sections = manifest.get("sections", []) specs = [s for s in sections if s.get("kind") == "job.spec"] records = [s for s in sections if s.get("kind") == "run.record"] if len(specs) != 1 or len(records) != 1 or records[0].get("exit_status") != 0: raise ResultRefused("expected one successful QVF job and run record") if ( records[0].get("program") != "vibe-qc" or records[0].get("run_status") != "converged" or records[0].get("job_ref") != specs[0]["id"] ): raise ResultRefused("QVF run record does not identify this successful vibe-qc job") spec = reader.read_json(specs[0]["members"]["spec"]) if reader.read_json(records[0]["members"]["input"]) != spec: raise ResultRefused("QVF executed input differs from its job specification") for key in ("method", "functional"): if provenance.get(key) is not None and provenance[key] != spec.get(key): raise ResultRefused(f"QVF provenance {key} conflicts with the executed input") species = document.species_by_id(job.species) expected = _job_spec(job, species) if expected is None or not species.conformers: raise ResultRefused("the recorded job has no comparable QVF input") expected = {k: v for k, v in expected.items() if k not in {"label", "structure_ref"}} actual = {k: v for k, v in spec.items() if k not in {"structure_ref"} and v is not None} if actual != expected: raise ResultRefused("QVF job specification does not match the recorded job") ref = specs[0].get("structure_ref") or spec.get("structure_ref") structures = [s for s in sections if s.get("id") == ref and s.get("kind") == "structure"] if len(structures) != 1: raise ResultRefused("QVF has no unique input structure") atoms = reader.read_json(structures[0]["members"]["structure"]).get("atoms") if not isinstance(atoms, list) or [ {k: v for k, v in a.items() if k != "atomic_number"} for a in atoms ] != _atoms(species.conformers[0]): raise ResultRefused("QVF input structure does not match this species and geometry") root_energy = manifest.get("energies", {}).get("total") scf = provenance.get("scf_energy") if root_energy is None and scf is None: raise ResultRefused("QVF has no total energy") energy = _quantity( root_energy if root_energy is not None else scf, "energies.total" if root_energy is not None else "provenance.scf_energy", ) if ( scf is not None and abs(_quantity(scf, "provenance.scf_energy").value - energy.value) > 1e-8 ): raise ResultRefused("QVF total and provenance energies conflict") version = records[0].get("program_version") if not isinstance(version, str) or not version: raise ResultRefused("QVF run record must identify its program version") result: dict[str, Any] = {"parser": "qvf/1", "program_version": version, "energy": energy} thermal = manifest.get("thermochemistry") if thermal: result["thermochemistry"] = { "temperature_k": _finite(thermal.get("temperature_k"), "temperature_k"), "pressure_atm": _finite(thermal.get("pressure_atm"), "pressure_atm"), "enthalpy": ResultEnergy( value=_finite(thermal.get("enthalpy_eh"), "enthalpy_eh"), source="thermochemistry.enthalpy_eh", ), "gibbs": ResultEnergy( value=_finite(thermal.get("gibbs_free_energy_eh"), "gibbs_free_energy_eh"), source="thermochemistry.gibbs_free_energy_eh", ), } if "freq" in job.spec.task and not thermal: raise ResultRefused("the frequency result has no thermochemistry") return result
[docs] def import_result( document: Document, run_id: str, data: bytes, *, format: Literal["orca", "qvf", "orca-hess"], source_name: str, replace: bool = False, ) -> Document: """Return a copy with one verified result. Never read a path supplied over HTTP. An ORCA optimisation's converged geometry is attached to the species as a conformer; ``orca-hess`` checks a Hessian file against a run that already has its verified output. """ if format not in {"orca", "qvf", "orca-hess"}: raise ResultRefused("result format must be orca, orca-hess or qvf") if format == "orca-hess": if not data or len(data) > MAX_RESULT_BYTES: raise ResultRefused("result must be nonempty and at most 32 MiB") return _with_hessian(document, run_id, data) job, run = _run(document, run_id) from vibesketch.chem.species import mol_from_smiles if not data or len(data) > MAX_RESULT_BYTES: raise ResultRefused("result must be nonempty and at most 32 MiB") digest = hashlib.sha256(data).hexdigest() if run.result is not None and run.result.source_sha256 != digest and not replace: raise ResultRefused( "this run already has a different result; explicitly request replacement" ) try: species = document.species_by_id(job.species) if not species.smiles: raise ResultRefused("result import requires a molecular identity") if any(atom.GetIsotope() for atom in mol_from_smiles(species.smiles).GetAtoms()): raise ResultRefused( "isotope-specific results are unsupported: input masses are not verified" ) values = _orca(document, job, data) if format == "orca" else _qvf(document, job, data) geometry: _Geometry | None = values.pop("geometry", None) result = RunResult( retrieved=datetime.now(UTC), source_name=source_name, source_sha256=digest, species_sha256=species_digest(document.species_by_id(job.species)), specification=job.spec.model_copy(deep=True), **values, ) except ResultRefused: raise except (ValueError, KeyError, TypeError, OSError, zipfile.BadZipFile) as error: raise ResultRefused(f"invalid result: {error}") from None updated = document.model_copy(deep=True) _, target = _run(updated, run_id) previous = target.result target.result = result owner = updated.species_by_id(job.species) if previous is not None: # a replaced result takes its geometry with it stale = _method(job.spec.task, previous.source_sha256) owner.conformers = [c for c in owner.conformers if c.method != stale] if geometry is not None: method = _method(job.spec.task, digest) if all(c.method != method for c in owner.conformers): symbols, xyz = geometry owner.conformers.append(Conformer(method=method, symbols=symbols, xyz_angstrom=xyz)) return updated
[docs] def reaction_energies(document: Document, selection: Selection) -> list[dict[str, Any]]: """Products minus reactants, with explicit source runs and balanced coefficients.""" if (selection.group is None) == (selection.runs is None): raise ResultRefused("choose exactly one submission group or an explicit list of runs") if selection.runs is not None and ( not selection.runs or len(selection.runs) != len(set(selection.runs)) ): raise ResultRefused("select distinct run identifiers") candidates = ( [(j, r) for j in document.jobs for r in j.runs or [] if r.group == selection.group] if selection.group is not None else [_run(document, r) for r in selection.runs or []] ) if not candidates: raise ResultRefused("the selected group has no runs") if len({r.id for _, r in candidates}) != len(candidates): raise ResultRefused("selected run identifiers are not unique") chosen: dict[str, tuple[Job, Run]] = {} for job, run in candidates: if job.species in chosen: raise ResultRefused(f"multiple selected runs for {job.species}; select one explicitly") chosen[job.species] = (job, run) reactions = {r.id: r for r in document.reactions} wanted = selection.reactions if selection.reactions is not None else list(reactions) if not wanted or len(set(wanted)) != len(wanted) or not set(wanted) <= reactions.keys(): raise ResultRefused("select known, distinct reactions") if not set(selection.coefficients) <= set(wanted): raise ResultRefused("coefficients name an unselected reaction") rows = [] for identifier in wanted: reaction = reactions[identifier] if reaction.arrow == "resonance": raise ResultRefused(f"{identifier}: a resonance arrow has no reaction energy") reactants, products = reaction.reactants, reaction.products if reaction.arrow == "retrosynthetic": reactants, products = products, reactants coeff = selection.coefficients.get( identifier, {"reactants": dict.fromkeys(reactants, 1.0), "products": dict.fromkeys(products, 1.0)}, ) if ( set(coeff) != {"reactants", "products"} or set(coeff["reactants"]) != set(reactants) or set(coeff["products"]) != set(products) ): raise ResultRefused( f"{identifier}: coefficients must name every species on both chemical sides" ) terms: list[tuple[float, Job, Run, RunResult]] = [] balance: dict[str, float] = {} specifications: set[str] = set() for side, sign in (("reactants", -1), ("products", 1)): for sid, value in coeff[side].items(): factor = _finite(value, "coefficient") if factor <= 0: raise ResultRefused("coefficients must be positive") if sid not in chosen: raise ResultRefused(f"{identifier}: no selected run for {sid}") job, run = chosen[sid] result = run.result species = document.species_by_id(sid) if result is None: raise ResultRefused(f"{identifier}: run {run.id} has no imported result") if ( result.parser != {"orca": "orca-text/1", "vibeqc": "qvf/1"}.get(job.program) or result.specification != job.spec or result.species_sha256 != species_digest(species) ): raise ResultRefused( f"{identifier}: run {run.id} result is stale for this job/species" ) spec = job.spec.model_dump(mode="json", exclude={"resources"}) specifications.add( json.dumps( {"program": job.program, "version": result.program_version, "spec": spec}, sort_keys=True, ) ) if not species.conformers: raise ResultRefused(f"{sid}: no input structure to check atom balance") for element, count in Counter(species.conformers[0].symbols).items(): balance[element] = balance.get(element, 0.0) + sign * factor * count balance["charge"] = balance.get("charge", 0.0) + sign * factor * species.charge terms.append((sign * factor, job, run, result)) if any(not math.isfinite(v) or abs(v) > 1e-8 for v in balance.values()): raise ResultRefused( f"{identifier}: atoms/charge are not balanced; " "supply explicit coefficients or complete the reaction" ) if len(specifications) != 1: raise ResultRefused( f"{identifier}: selected runs have incompatible programs/specifications" ) thermal = [r.thermochemistry for _, _, _, r in terms] if any(t is not None for t in thermal) and not all(t is not None for t in thermal): raise ResultRefused(f"{identifier}: thermochemistry is missing from some selected runs") temperatures = {(t.temperature_k, t.pressure_atm) for t in thermal if t is not None} if len(temperatures) > 1: raise ResultRefused(f"{identifier}: thermochemistry temperature/pressure differs") energies: dict[str, Any] = {} for key in ("energy", "enthalpy", "gibbs"): if key != "energy" and not temperatures: continue try: value = math.fsum( factor * ( result.energy.value if key == "energy" else getattr(result.thermochemistry, key).value ) for factor, _, _, result in terms ) except (OverflowError, ValueError): raise ResultRefused("reaction energy overflow") from None if not math.isfinite(value): raise ResultRefused("reaction energy overflow") converted = value * HARTREE_TO_KCAL_MOL if not math.isfinite(converted): raise ResultRefused("reaction energy conversion overflow") energies[key] = {"value": value, "units": "hartree", "kcal_mol": converted} rows.append( { "reaction": identifier, "group": selection.group, "coefficients": coeff, "specification": terms[0][1].spec.model_dump(mode="json", exclude_none=True), "program": terms[0][1].program, "energies": energies, "temperature_k": next(iter(temperatures))[0] if temperatures else None, "pressure_atm": next(iter(temperatures))[1] if temperatures else None, "sources": [ { "species": job.species, "job": job.id, "run": run.id, "coefficient": factor, "result": result.model_dump(mode="json", exclude_none=True), } for factor, job, run, result in terms ], } ) return rows
[docs] def result_figure(document: Document, selection: Selection) -> bytes: """A derived SVG preview: energies below arrows, without changing the sketch.""" from xml.sax.saxutils import escape from vibesketch.document.derive import _roles, derive from vibesketch.document.ket import read_scene, text_node from vibesketch.document.model import Canvas from vibesketch.render import render_svg from vibesketch.render.fonts import advance from vibesketch.render.styles import resolve rows = reaction_energies(document, selection) derived = derive(document).document if [r.model_dump() for r in derived.reactions] != [r.model_dump() for r in document.reactions]: raise ResultRefused( "the drawn reactions do not match the document; update the drawing first" ) scene = read_scene(document.canvas.data) roles = _roles(scene) if len(roles) != len(document.reactions): raise ResultRefused("each reaction needs a drawn arrow for a result figure") copy = document.model_copy(deep=True) canvas: dict[str, Any] = copy.canvas.data preset = resolve(document.style.preset) unit = preset.size / preset.bond line_step = 1.3 * unit boxes = [fragment.bounds for fragment in scene.fragments] boxes += [ ( t.centre[0] - advance(t.text) * unit / 2, t.centre[1] - unit, t.centre[0] + advance(t.text) * unit / 2, t.centre[1] + unit, ) for t in scene.texts ] labels = {"energy": "ΔE", "enthalpy": "ΔH", "gibbs": "ΔG"} for row in rows: index = next( i for i, reaction in enumerate(document.reactions) if reaction.id == row["reaction"] ) arrow = roles[index][0] x = (arrow.tail[0] + arrow.head[0]) / 2 y = min(arrow.tail[1], arrow.head[1]) - 1.0 nearby = [ t.centre[1] for t in scene.texts if abs(t.centre[0] - x) < 3 and y - 2.5 <= t.centre[1] <= y + 1 ] if nearby: y = min(y, min(nearby) - 0.7) lines = [ f"{labels[key]} = {value['kcal_mol']:.2f} kcal/mol" for key, value in row["energies"].items() ] spec = row["specification"] lines.append(f"{spec['method']}/{spec.get('basis', 'composite')}") if row["temperature_k"] is not None: lines.append(f"{row['temperature_k']:g} K, {row['pressure_atm']:g} atm") width = max(advance(line) for line in lines) * unit + 0.4 height = (len(lines) - 1) * line_step + unit for _ in range(12): box = (x - width / 2, y - height, x + width / 2, y + unit) if all( box[2] < b[0] - 0.2 or box[0] > b[2] + 0.2 or box[3] < b[1] - 0.2 or box[1] > b[3] + 0.2 for b in boxes ): boxes.append(box) break y -= 0.6 else: raise ResultRefused( "no clear space below the reaction arrow; spread out the drawing first" ) for line in lines: canvas["root"]["nodes"].append(text_node(line, (x, y))) y -= line_step copy.canvas = Canvas(format="ket", data=copy.canvas.data) svg = render_svg(copy) metadata = ( "<metadata id='reaction-results'>" + escape(json.dumps(rows, allow_nan=False)) + "</metadata>" ).encode() index = svg.index(b">", svg.index(b"<svg")) + 1 return svg[:index] + metadata + svg[index:]