"""Explicit aqueous reduction and finite-cluster binding recipes (#146).
Ho, Coote, Cramer and Truhlar, Organic Electrochemistry, 5th ed., ch. 4
(2015), eq. 33 and Table 3: https://comp.chem.umn.edu/truhlar/docs/C86preprint.pdf.
Reduction energies include the chosen electron chemical potential. Solute
RRHO free energies at 1 atm receive RT ln(cRT/p) for a 1 M standard state;
proton solvation must use equal gas/solution concentrations (Kelly et al.).
"""
from __future__ import annotations
import math
from collections import Counter
from typing import Any
from pydantic import TypeAdapter
from vibesketch.document.electrochem import FARADAY, balance
from vibesketch.document.model import (
BindingAnalysis,
CalculationStage,
Document,
GeometryRef,
Job,
JobSpec,
Recipe,
RecipeAnalysis,
RedoxAnalysis,
Stage,
)
from vibesketch.document.recipe import set_recipe
from vibesketch.document.results import ResultRefused, _run, species_digest
from vibesketch.units import HARTREE_TO_KCAL_MOL
# CODATA 2018, exact R = N_A k_B, J/(mol K); 1 atm = 101325 Pa.
R = 8.31446261815324
HARTREE_TO_KJ_MOL = HARTREE_TO_KCAL_MOL * 4.184
# Choices are returned for review, never selected implicitly. The different SHE
# conventions are not interchangeable physical definitions. Isse/Gennaro use a
# zero-temperature electron and chemical (not real interfacial) proton potential.
CONSTANTS: dict[str, dict[str, Any]] = {
"she-trasatti-1986": {
"value": 4.44,
"units": "V",
"source": "Trasatti (1986), eq. 19, https://doi.org/10.1351/pac198658070955",
},
"she-isse-gennaro-2010": {
"value": 4.281,
"units": "V",
"source": "Isse and Gennaro (2010), https://doi.org/10.1021/jp100402x",
},
"proton-kelly-2006": {
"value": -265.9 * 4.184,
"units": "kJ/mol",
"source": "Kelly, Cramer and Truhlar (2006), equal-concentration gas/solution transfer, https://doi.org/10.1021/jp063552y",
},
}
# Manufacturer's table, 25 C. No extrapolation or liquid-junction correction.
_AGCL_SOURCE = "Metrohm, Reference electrodes, Table 1 (2021), https://www.metrohm.cn/content/dam/metrohm/shared/documents/fact-sheets/81098050EN.pdf"
# SCE offset from Isse/Gennaro's aqueous reference convention (2010).
_OFFSETS = {
"SHE": (0.0, "SHE definition"),
"Ag/AgCl (sat. KCl)": (0.197, _AGCL_SOURCE),
"Ag/AgCl (3 M KCl)": (0.207, _AGCL_SOURCE),
"SCE": (0.241, "Isse and Gennaro (2010), https://doi.org/10.1021/jp100402x"),
}
[docs]
def constants() -> dict[str, Any]:
"""Cited values to choose explicitly, including aqueous reference offsets."""
import copy
return {
"choices": copy.deepcopy(CONSTANTS),
"offsets_vs_she_v": {
name: {"value": value, "source": source} for name, (value, source) in _OFFSETS.items()
},
"temperature_k": 298.15,
}
[docs]
def reduction_potential(delta_g_kj_mol: float, electrons: int, absolute_she_v: float) -> float:
"""E vs SHE = -delta G/(nF) - absolute SHE; delta G includes electrons."""
if isinstance(electrons, bool) or not isinstance(electrons, int) or electrons < 1:
raise ResultRefused("a reduction consumes a positive integer electron count")
if not all(math.isfinite(v) for v in (delta_g_kj_mol, absolute_she_v)):
raise ResultRefused("redox energies and reference potentials must be finite")
value = -delta_g_kj_mol * 1000 / (electrons * FARADAY) - absolute_she_v
if not math.isfinite(value):
raise ResultRefused("redox potential overflow")
return value
def _redox_terms(
document: Document, analysis: RedoxAnalysis
) -> tuple[dict[str, float], float, int]:
reaction = next((r for r in document.reactions if r.id == analysis.reaction), None)
if reaction is None or reaction.electrons is None:
raise ResultRefused("choose a drawn half-reaction with electrons")
if reaction.arrow in ("resonance", "retrosynthetic") or reaction.electrons.side != "reactants":
raise ResultRefused("redox recipes require a forward reduction with electrons consumed")
if not balance(document, reaction.id)["balanced"]:
raise ResultRefused("the half-reaction must balance atoms and charge including electrons")
if any(state.state != "aq" for state in reaction.states):
raise ResultRefused("aqueous redox recipes accept only aqueous participant states")
terms: dict[str, float] = {}
protons = 0.0
for side, sign in ((reaction.reactants, -1), (reaction.products, 1)):
for sid in side:
species = document.species_by_id(sid)
factor = sign * reaction.coefficient(sid)
if species.smiles == "[H+]":
protons += factor
continue
if not species.smiles or species.solid or species.nanostructure:
raise ResultRefused("aqueous redox recipes require concrete molecular species")
if species.conformers:
from rdkit import Chem
from vibesketch.chem.species import mol_from_smiles
graph = Chem.AddHs(mol_from_smiles(species.smiles))
if Counter(a.GetSymbol() for a in graph.GetAtoms()) != Counter(
species.conformers[0].symbols
):
raise ResultRefused(f"{sid}: conformer atoms do not match the half-reaction")
terms[sid] = terms.get(sid, 0.0) + factor
terms = {sid: factor for sid, factor in terms.items() if factor}
if not terms:
raise ResultRefused("the half-reaction has no molecular energy terms")
if protons and analysis.proton_solvation_kj_mol is None:
raise ResultRefused(
"proton-coupled reduction needs explicit cited gas and solvation energies"
)
return terms, protons, reaction.electrons.count
def _binding_terms(document: Document, analysis: BindingAnalysis) -> dict[str, float]:
ids = [analysis.complex, analysis.cluster, analysis.adsorbate]
if len(set(ids)) != 3:
raise ResultRefused("complex, cluster and adsorbate must be distinct species")
if analysis.surface and analysis.surface not in {s.id for s in document.surfaces}:
raise ResultRefused("the binding recipe names an unknown schematic surface")
atoms: Counter[str] = Counter()
charge = 0
for sid, sign in zip(ids, [1, -1, -1], strict=True):
try:
species = document.species_by_id(sid)
except KeyError:
raise ResultRefused(f"unknown binding species {sid}") from None
if not species.smiles:
raise ResultRefused(f"{sid}: a finite molecular identity is required by result import")
if not species.conformers or species.solid:
raise ResultRefused(
f"{sid}: supply an explicit finite-cluster geometry; "
"a surface picture is not a model"
)
for element in species.conformers[0].symbols:
atoms[element] += sign
charge += sign * species.charge
if any(atoms.values()) or charge:
raise ResultRefused("complex atoms and charge must equal cluster plus adsorbate")
return dict(zip(ids, [1.0, -1.0, -1.0], strict=True))
def _terms(document: Document, analysis: RecipeAnalysis) -> dict[str, float]:
return (
_redox_terms(document, analysis)[0]
if isinstance(analysis, RedoxAnalysis)
else _binding_terms(document, analysis)
)
[docs]
def with_energy_recipe(
document: Document,
identifier: str,
analysis: RecipeAnalysis | dict[str, Any],
spec: JobSpec | dict[str, Any],
program: str,
) -> Document:
"""Save an ordinary calculation stage plus explicit reduction/binding analysis."""
chosen = TypeAdapter(RecipeAnalysis).validate_python(analysis)
job_spec = JobSpec.model_validate(spec)
terms = _terms(document, chosen)
if isinstance(chosen, RedoxAnalysis):
if (
job_spec.task != "opt+freq"
or job_spec.solvent is None
or job_spec.solvent.name.lower() != "water"
):
raise ResultRefused("redox recipes require opt+freq in explicit CPCM/SMD water")
elif job_spec.task != "sp":
raise ResultRefused("binding energies require sp on explicitly supplied geometries")
calculation = CalculationStage.model_validate({"program": program, "spec": job_spec})
recipe = Recipe(
id=identifier,
species=list(terms),
analysis=chosen,
stages=[
Stage(
id="energy",
kind="calculation",
geometry=GeometryRef(conformer=0),
calculation=calculation,
)
],
)
return set_recipe(document, recipe)
def _chosen(document: Document, identifier: str) -> tuple[Recipe, CalculationStage]:
recipe = next((r for r in document.recipes if r.id == identifier), None)
if recipe is None or recipe.analysis is None:
raise ResultRefused("choose a saved redox or binding recipe")
if (
len(recipe.stages) != 1
or recipe.stages[0].calculation is None
or recipe.stages[0].geometry != GeometryRef(conformer=0)
):
raise ResultRefused("energy recipe requires one calculation stage from conformer zero")
calculation = recipe.stages[0].calculation
# Revalidate on every use: edited reactions and handcrafted recipes cannot
# bypass the builder's scientific admission checks.
rebuilt = with_energy_recipe(
document, identifier, recipe.analysis, calculation.spec, calculation.program
)
wanted = next(r for r in rebuilt.recipes if r.id == identifier)
if recipe.species != wanted.species:
raise ResultRefused("recipe species no longer match its reaction or binding partners")
return recipe, calculation
[docs]
def compile_energy_recipe(document: Document, identifier: str) -> tuple[Document, dict[str, Any]]:
"""Compile ordinary jobs and return input files, or refuse the entire recipe."""
from vibesketch.inputs import get_writer
recipe, calculation = _chosen(document, identifier)
if isinstance(recipe.analysis, BindingAnalysis) and recipe.analysis.counterpoise:
# Boys and Bernardi, Mol. Phys. 19, 553 (1970), doi:10.1080/00268977000101561.
# E_AB^AB - E_A^AB - E_B^AB needs ghost-basis fragment jobs at complex geometry.
raise ResultRefused(
"counterpoise needs ghost-basis fragment jobs at the complex geometry; "
"neither writer exposes that contract, so no uncorrected substitute is compiled"
)
writer = get_writer(calculation.program)
jobs = []
files: list[dict[str, str]] = []
refusals: list[str] = []
taken = (
{s.id for s in document.species}
| {r.id for r in document.reactions}
| {j.id for j in document.jobs}
| {r.id for r in document.recipes}
| {s.id for s in document.surfaces}
| {g.id for g in document.generics}
)
for sid in recipe.species:
species = document.species_by_id(sid)
reasons = writer.refusals(species, calculation.spec)
from vibesketch.elements import ATOMIC_NUMBER
if species.conformers and any(ATOMIC_NUMBER[e] > 36 for e in species.conformers[0].symbols):
reasons.append(
"energy recipes do not yet record an explicit heavy-element "
"ECP/relativistic contract; refusing this model rather than assuming one"
)
refusals.extend(f"{sid}: {reason}" for reason in reasons)
job_id = f"{identifier}-{sid}"
if job_id in taken:
refusals.append(f"job id {job_id} already exists; use a new recipe id")
if not reasons:
bundle = writer.write(species, calculation.spec)
files.extend(
{"job": job_id, "name": f.name, "content": f.content, "media_type": f.media_type}
for f in bundle.files
)
jobs.append(Job(id=job_id, species=sid, program=calculation.program, spec=calculation.spec))
if refusals:
raise ResultRefused("; ".join(refusals))
data = document.to_json_data()
data["jobs"] = [j.model_dump(mode="json", exclude_none=True) for j in [*document.jobs, *jobs]]
return Document.model_validate(data), {
"recipe": identifier,
"files": files,
"jobs": [j.id for j in jobs],
}
[docs]
def energy_recipe_result(document: Document, identifier: str, runs: list[str]) -> dict[str, Any]:
"""Evaluate explicitly selected, imported result runs; never select the latest run."""
recipe, calculation = _chosen(document, identifier)
analysis = recipe.analysis
assert analysis is not None
if isinstance(analysis, BindingAnalysis) and analysis.counterpoise:
raise ResultRefused(
"counterpoise results require verified ghost-basis jobs, which are unsupported"
)
terms = _terms(document, analysis)
if len(runs) != len(set(runs)) or len(runs) != len(terms):
raise ResultRefused("select exactly one distinct imported run per recipe species")
values: dict[str, float] = {}
sources = []
versions = set()
for run_id in runs:
job, run = _run(document, run_id)
result = run.result
species = document.species_by_id(job.species)
if job.species not in terms or job.species in values:
raise ResultRefused("selected runs must cover exactly the recipe species")
if (
result is None
or result.specification != job.spec
or result.species_sha256 != species_digest(species)
or result.parser != {"orca": "orca-text/1", "vibeqc": "qvf/1"}.get(job.program)
):
raise ResultRefused(f"{run_id}: missing or stale imported result")
if job.program != calculation.program or job.spec.model_dump(
exclude={"resources"}
) != calculation.spec.model_dump(exclude={"resources"}):
raise ResultRefused(f"{run_id}: result does not match recipe program/specification")
versions.add(result.program_version)
value = result.energy.value * HARTREE_TO_KJ_MOL
if isinstance(analysis, RedoxAnalysis):
thermal = result.thermochemistry
if (
thermal is None
or abs(thermal.temperature_k - analysis.temperature_k) > 1e-6
or abs(thermal.pressure_atm - 1) > 1e-6
):
raise ResultRefused("redox results require Gibbs energies at 298.15 K and 1 atm")
if result.vibrations is None or result.vibrations.imaginary:
raise ResultRefused("redox results require frequencies confirming a minimum")
value = thermal.gibbs.value * HARTREE_TO_KJ_MOL
values[job.species] = value
sources.append(
{
"species": job.species,
"job": job.id,
"run": run.id,
"coefficient": terms[job.species],
"result": result.model_dump(mode="json", exclude_none=True),
}
)
if len(versions) != 1:
raise ResultRefused("selected results use different program versions")
delta = math.fsum(terms[sid] * value for sid, value in values.items())
report: dict[str, Any] = {
"recipe": identifier,
"analysis": analysis.model_dump(mode="json"),
"specification": calculation.spec.model_dump(mode="json", exclude_none=True),
"program": calculation.program,
"sources": sources,
"temperature_k": None,
"pressure_atm": None,
}
if isinstance(analysis, RedoxAnalysis):
_, protons, electrons = _redox_terms(document, analysis)
correction = (
R * analysis.temperature_k * math.log(1000 * R * analysis.temperature_k / 101325) / 1000
)
delta += sum(terms.values()) * correction
if protons:
assert (
analysis.proton_gas_gibbs_kj_mol is not None
and analysis.proton_solvation_kj_mol is not None
)
delta += protons * (
analysis.proton_gas_gibbs_kj_mol.value
+ correction
+ analysis.proton_solvation_kj_mol.value
)
delta -= electrons * analysis.electron_gibbs_kj_mol.value
standard = reduction_potential(delta, electrons, analysis.absolute_she_v.value)
# Nernst: proton activity=10^-pH; all other activities unity.
nernst = R * analysis.temperature_k * math.log(10) / FARADAY
at_ph = standard + protons * nernst * analysis.ph / electrons
if analysis.reference == "RHE":
offset, source = (
-nernst * analysis.ph,
"Nernst equation, RHE at unit H2 pressure; Ho et al. (2015), ch. 4",
)
else:
offset, source = _OFFSETS[analysis.reference]
report.update(
reaction=analysis.reaction,
electrons=electrons,
proton_coefficient=protons,
standard_state_correction_kj_mol=correction,
standard_potential_she_v=standard,
potential_v=at_ph - offset,
reference_offset_v=offset,
reference_source=source,
temperature_k=analysis.temperature_k,
pressure_atm=1.0,
)
assert calculation.spec.solvent is not None
report["annotation"] = [
f"E = {at_ph - offset:.3f} V vs {analysis.reference}, pH {analysis.ph:g}",
f"{calculation.spec.solvent.model} water",
]
energy_kind = "gibbs"
else:
report["annotation"] = ["Finite-cluster binding energy (uncorrected)"]
energy_kind = "energy"
if not math.isfinite(delta):
raise ResultRefused("energy evaluation overflow")
report["energies"] = {
energy_kind: {
"value": delta / HARTREE_TO_KJ_MOL,
"units": "hartree",
"kcal_mol": delta / 4.184,
"kj_mol": delta,
}
}
return report