From a5487980920a450fbaff1ffcd5bdbf0dde16a6ae Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Mon, 6 Jul 2026 22:37:11 +0200 Subject: [PATCH 1/2] Add reaction free energy and reduction potential helpers Introduce ThermoScreening.thermo.reactions with two pure post-processing helpers that combine per-species absolute Gibbs free energies (Thermo.total_EeGtot) computed by any engine: - reaction_free_energy(reactants, products, unit): dG = sum(products) - sum(reactants), with (coefficient, Thermo) stoichiometry tuples and H / eV / kcal / kJ output units. - reduction_potential(oxidized, reduced, n_electrons, reference_potential): E = -dG_red / n, referenced to the SHE by default (override for a calibrated reference or 0.0 for the absolute potential). Both are exported from ThermoScreening.thermo. Adds CI-safe exact-math tests plus a skippable native-xtb benzoquinone end-to-end test, and docs (API reference + a usage section) with a prominent caveat that GFN-xTB/DFTB absolute redox potentials are quantitatively poor and best used for relative trends or with a calibrated reference. --- ThermoScreening/thermo/__init__.py | 1 + ThermoScreening/thermo/reactions.py | 111 ++++++++++++++++++++++++++++ docs/api.rst | 6 ++ docs/usage.rst | 35 +++++++++ tests/thermo/test_reactions.py | 93 +++++++++++++++++++++++ 5 files changed, 246 insertions(+) create mode 100644 ThermoScreening/thermo/reactions.py create mode 100644 tests/thermo/test_reactions.py diff --git a/ThermoScreening/thermo/__init__.py b/ThermoScreening/thermo/__init__.py index 4eb2cd4..0ae923e 100644 --- a/ThermoScreening/thermo/__init__.py +++ b/ThermoScreening/thermo/__init__.py @@ -8,3 +8,4 @@ from .thermo import Thermo from .screening import screen from .conformers import generate as generate_conformers, write_conformers +from .reactions import reaction_free_energy, reduction_potential diff --git a/ThermoScreening/thermo/reactions.py b/ThermoScreening/thermo/reactions.py new file mode 100644 index 0000000..eeb9f1e --- /dev/null +++ b/ThermoScreening/thermo/reactions.py @@ -0,0 +1,111 @@ +"""Reaction and redox thermochemistry from computed ``Thermo`` objects. + +These are pure post-processing helpers that combine the absolute Gibbs free +energies (``Thermo.total_EeGtot()`` = electronic + thermal) of individual species +into reaction free energies and reduction potentials. Compute each species with +any engine (``dftbplus_thermo`` / ``xtb_thermo`` / ``xtb_cli_thermo``) under the +same conditions -- ideally the same solvent, and open-shell/charged as +appropriate -- then combine them here. +""" + +from ..utils.physicalConstants import PhysicalConstants + +# Absolute potential of the standard hydrogen electrode (V). Convention-dependent +# (values from ~4.28 to ~4.44 V are used); this is the IUPAC-recommended value. +# Override ``reference_potential`` for a calibrated reference. +SHE_ABSOLUTE_POTENTIAL = 4.44 + +_HARTREE_TO_EV = PhysicalConstants["H"] / PhysicalConstants["eV"] +_HARTREE_TO_KCAL_PER_MOL = PhysicalConstants["H"] * PhysicalConstants["N_A"] / 4184.0 +_HARTREE_TO_KJ_PER_MOL = PhysicalConstants["H"] * PhysicalConstants["N_A"] / 1000.0 + +_UNIT_FACTORS = { + "H": 1.0, + "eV": _HARTREE_TO_EV, + "kcal": _HARTREE_TO_KCAL_PER_MOL, + "kJ": _HARTREE_TO_KJ_PER_MOL, +} + + +def _total_gibbs(species): + """ + Absolute Gibbs free energy (Hartree) of a species entry. + + ``species`` is a ``Thermo`` (coefficient 1) or a ``(coefficient, Thermo)`` + tuple. + """ + if isinstance(species, tuple): + coefficient, thermo = species + else: + coefficient, thermo = 1.0, species + return coefficient * thermo.total_EeGtot() + + +def reaction_free_energy(reactants, products, unit="kcal"): + """ + Reaction free energy dG = sum(products) - sum(reactants). + + Parameters + ---------- + reactants, products : iterable + Each entry is a ``Thermo`` (stoichiometric coefficient 1) or a + ``(coefficient, Thermo)`` tuple. All species should be computed with the + same engine and conditions for the difference to be meaningful. + unit : str + Output unit: ``"kcal"`` (kcal/mol, default), ``"kJ"`` (kJ/mol), ``"eV"`` + (per particle), or ``"H"`` (Hartree per particle). + + Returns + ------- + float + The reaction free energy in the requested unit. + + Raises + ------ + ValueError + If ``unit`` is not supported. + """ + try: + factor = _UNIT_FACTORS[unit] + except KeyError: + known = ", ".join(_UNIT_FACTORS) + raise ValueError(f"Unknown unit {unit!r}; choose one of: {known}.") + + delta = ( + sum(_total_gibbs(product) for product in products) + - sum(_total_gibbs(reactant) for reactant in reactants) + ) + return delta * factor + + +def reduction_potential( + oxidized, reduced, n_electrons=1, reference_potential=SHE_ABSOLUTE_POTENTIAL +): + """ + Reduction potential (V) for ``Ox + n e- -> Red``. + + Uses ``E = -dG_red / (n F) - E_reference`` with ``dG_red = G(Red) - G(Ox)`` + and the free energy of the electron taken as zero. Because one electron-volt + per electron is one volt, the absolute potential simplifies to + ``-dG_red[eV] / n``, which is then referenced to ``reference_potential``. + + Parameters + ---------- + oxidized, reduced : Thermo + The oxidised and reduced species, computed consistently (same engine, + same solvent, and open-shell/charged as appropriate). + n_electrons : int + Number of electrons transferred. Default 1. + reference_potential : float + Absolute potential (V) of the reference electrode to report against. + Defaults to the SHE (:data:`SHE_ABSOLUTE_POTENTIAL`); pass ``0.0`` for the + absolute reduction potential. + + Returns + ------- + float + The reduction potential in volts (versus ``reference_potential``). + """ + delta_g_hartree = reduced.total_EeGtot() - oxidized.total_EeGtot() + absolute = -delta_g_hartree * _HARTREE_TO_EV / n_electrons + return absolute - reference_potential diff --git a/docs/api.rst b/docs/api.rst index a5e3e6e..21ca726 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -23,6 +23,12 @@ Conformer generation .. autofunction:: ThermoScreening.thermo.conformers.generate .. autofunction:: ThermoScreening.thermo.conformers.write_conformers +Reactions and redox +------------------- + +.. autofunction:: ThermoScreening.thermo.reactions.reaction_free_energy +.. autofunction:: ThermoScreening.thermo.reactions.reduction_potential + Thermochemistry core -------------------- diff --git a/docs/usage.rst b/docs/usage.rst index 2ac9aea..423a03b 100644 --- a/docs/usage.rst +++ b/docs/usage.rst @@ -72,6 +72,41 @@ Conformers, then screen the ensemble: write_conformers(conformers, "butanol_confs") results = screen("butanol_confs", engine="xtb-cli") +Reactions and redox +------------------- + +Once you have a ``Thermo`` object per species (from any engine, under the same +conditions), combine them into reaction free energies and reduction potentials: + +.. code-block:: python + + from ase.build import molecule + from ThermoScreening.thermo.api import xtb_cli_thermo + from ThermoScreening.thermo.conformers import generate + from ThermoScreening.thermo import reaction_free_energy, reduction_potential + + # reaction free energy (stoichiometry via (coefficient, Thermo) tuples) + h2 = xtb_cli_thermo(molecule("H2")) + o2 = xtb_cli_thermo(molecule("O2"), spin=1.0) # S = 1, triplet + h2o = xtb_cli_thermo(molecule("H2O")) + dG = reaction_free_energy([(2, h2), o2], [(2, h2o)], unit="kcal") # 2 H2 + O2 -> 2 H2O + + # one-electron reduction potential (Ox + e- -> Red), both in solvent + mol = generate("O=C1C=CC(=O)C=C1", max_conformers=1)[0] # benzoquinone + neutral = xtb_cli_thermo(mol, charge=0, solvent="water") + anion = xtb_cli_thermo(mol, charge=-1, solvent="water") # radical anion, auto open-shell + E = reduction_potential(neutral, anion) # vs SHE (default reference 4.44 V) + E_abs = reduction_potential(neutral, anion, reference_potential=0.0) + +.. warning:: + + ``reaction_free_energy`` / ``reduction_potential`` are exact given the input + energies, but the *accuracy* is set by the underlying method. GFN-xTB and DFTB + give poor **absolute** electron affinities and redox potentials (benzoquinone's + GFN2 EA is ~7 eV vs ~1.9 eV experimental). Use them for **relative** trends + across similar species, with a higher-accuracy method, or with a + ``reference_potential`` calibrated against experiment. + End-to-end example ------------------ diff --git a/tests/thermo/test_reactions.py b/tests/thermo/test_reactions.py new file mode 100644 index 0000000..a66a08d --- /dev/null +++ b/tests/thermo/test_reactions.py @@ -0,0 +1,93 @@ +import math +import os +import shutil + +import pytest + +from ThermoScreening.thermo.reactions import ( + SHE_ABSOLUTE_POTENTIAL, + reaction_free_energy, + reduction_potential, +) + +_H_TO_EV = 27.211386245988034 +_H_TO_KCAL = 627.5094740631 # 1 Hartree in kcal/mol + + +class _FakeThermo: + """A stand-in exposing only what the reactions helpers use.""" + + def __init__(self, eegtot): + self._eegtot = eegtot + + def total_EeGtot(self): + return self._eegtot + + +def test_reaction_free_energy_units(): + a, b, c = _FakeThermo(-1.0), _FakeThermo(-2.0), _FakeThermo(-3.5) + # A + B -> C : dG = -3.5 - (-3.0) = -0.5 Hartree + assert reaction_free_energy([a, b], [c], unit="H") == pytest.approx(-0.5) + assert reaction_free_energy([a, b], [c], unit="eV") == pytest.approx(-0.5 * _H_TO_EV) + assert reaction_free_energy([a, b], [c], unit="kcal") == pytest.approx( + -0.5 * _H_TO_KCAL, rel=1e-4 + ) + + +def test_reaction_free_energy_stoichiometry(): + a, b = _FakeThermo(-1.0), _FakeThermo(-2.5) + # 2 A -> B : dG = -2.5 - 2*(-1.0) = -0.5 Hartree + assert reaction_free_energy([(2.0, a)], [b], unit="H") == pytest.approx(-0.5) + + +def test_reaction_free_energy_rejects_unknown_unit(): + with pytest.raises(ValueError, match="Unknown unit"): + reaction_free_energy([_FakeThermo(0.0)], [_FakeThermo(0.0)], unit="furlong") + + +def test_reduction_potential_absolute_and_vs_reference(): + ox, red = _FakeThermo(0.0), _FakeThermo(-0.5) + # dG_red = -0.5 Hartree -> E_abs = 0.5 * 27.2114 V + assert reduction_potential(ox, red, reference_potential=0.0) == pytest.approx( + 0.5 * _H_TO_EV + ) + assert reduction_potential(ox, red) == pytest.approx( + 0.5 * _H_TO_EV - SHE_ABSOLUTE_POTENTIAL + ) + + +def test_reduction_potential_n_electrons(): + ox, red = _FakeThermo(0.0), _FakeThermo(-1.0) + # n = 2 : E_abs = 1.0 * 27.2114 / 2 + assert reduction_potential(ox, red, n_electrons=2, reference_potential=0.0) == pytest.approx( + _H_TO_EV / 2 + ) + + +def test_public_api_exported(): + from ThermoScreening.thermo import reaction_free_energy as rfe + from ThermoScreening.thermo import reduction_potential as rp + from ThermoScreening.thermo import reactions + + assert rfe is reactions.reaction_free_energy + assert rp is reactions.reduction_potential + + +xtb_available = shutil.which("xtb") is not None or "XTB_COMMAND" in os.environ + + +@pytest.mark.skipif(not xtb_available, reason="the native xtb binary is not available.") +def test_reduction_potential_end_to_end(tmp_path): + # benzoquinone + e- -> radical anion, in water, via xtb-cli. GFN2 absolute + # redox potentials are quantitatively poor, so this only checks that the + # pipeline runs and the reduction is (correctly) favorable for a good acceptor. + from ThermoScreening.thermo.api import xtb_cli_thermo + from ThermoScreening.thermo.conformers import generate + + bq = generate("O=C1C=CC(=O)C=C1", max_conformers=3)[0] + neutral = xtb_cli_thermo(bq, charge=0, solvent="water", directory=str(tmp_path / "n")) + anion = xtb_cli_thermo(bq, charge=-1, solvent="water", directory=str(tmp_path / "a")) + + e_abs = reduction_potential(neutral, anion, reference_potential=0.0) + assert math.isfinite(e_abs) + assert e_abs > 0 # reduction of a quinone is favorable From 47b117d22e8d08157e5e74671a856baa8cec7f6e Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Mon, 6 Jul 2026 22:41:55 +0200 Subject: [PATCH 2/2] Guard zero electrons, test kJ unit, reuse cal constant Address review feedback: raise ValueError for n_electrons == 0 rather than an opaque ZeroDivisionError, add a kJ/mol unit assertion, and derive the kcal factor from PhysicalConstants['cal'] for codebase consistency. --- ThermoScreening/thermo/reactions.py | 12 +++++++++++- tests/thermo/test_reactions.py | 10 ++++++++++ 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/ThermoScreening/thermo/reactions.py b/ThermoScreening/thermo/reactions.py index eeb9f1e..ce6052c 100644 --- a/ThermoScreening/thermo/reactions.py +++ b/ThermoScreening/thermo/reactions.py @@ -16,7 +16,9 @@ SHE_ABSOLUTE_POTENTIAL = 4.44 _HARTREE_TO_EV = PhysicalConstants["H"] / PhysicalConstants["eV"] -_HARTREE_TO_KCAL_PER_MOL = PhysicalConstants["H"] * PhysicalConstants["N_A"] / 4184.0 +_HARTREE_TO_KCAL_PER_MOL = ( + PhysicalConstants["H"] * PhysicalConstants["N_A"] / (PhysicalConstants["cal"] * 1000.0) +) _HARTREE_TO_KJ_PER_MOL = PhysicalConstants["H"] * PhysicalConstants["N_A"] / 1000.0 _UNIT_FACTORS = { @@ -105,7 +107,15 @@ def reduction_potential( ------- float The reduction potential in volts (versus ``reference_potential``). + + Raises + ------ + ValueError + If ``n_electrons`` is zero. """ + if n_electrons == 0: + raise ValueError("n_electrons must be non-zero.") + delta_g_hartree = reduced.total_EeGtot() - oxidized.total_EeGtot() absolute = -delta_g_hartree * _HARTREE_TO_EV / n_electrons return absolute - reference_potential diff --git a/tests/thermo/test_reactions.py b/tests/thermo/test_reactions.py index a66a08d..d8b168f 100644 --- a/tests/thermo/test_reactions.py +++ b/tests/thermo/test_reactions.py @@ -12,6 +12,7 @@ _H_TO_EV = 27.211386245988034 _H_TO_KCAL = 627.5094740631 # 1 Hartree in kcal/mol +_H_TO_KJ = 2625.4996394798254 # 1 Hartree in kJ/mol class _FakeThermo: @@ -32,6 +33,9 @@ def test_reaction_free_energy_units(): assert reaction_free_energy([a, b], [c], unit="kcal") == pytest.approx( -0.5 * _H_TO_KCAL, rel=1e-4 ) + assert reaction_free_energy([a, b], [c], unit="kJ") == pytest.approx( + -0.5 * _H_TO_KJ, rel=1e-4 + ) def test_reaction_free_energy_stoichiometry(): @@ -64,6 +68,12 @@ def test_reduction_potential_n_electrons(): ) +def test_reduction_potential_rejects_zero_electrons(): + ox, red = _FakeThermo(0.0), _FakeThermo(-0.5) + with pytest.raises(ValueError, match="non-zero"): + reduction_potential(ox, red, n_electrons=0) + + def test_public_api_exported(): from ThermoScreening.thermo import reaction_free_energy as rfe from ThermoScreening.thermo import reduction_potential as rp