Add reaction free energy and reduction potential helpers

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

ThermoScreening/thermo/reactions.py

@@ -0,0 +1,121 @@
+"""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"] / (PhysicalConstants["cal"] * 1000.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``).
+
+    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

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
 --------------------
 

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
 ------------------
 

tests/thermo/test_reactions.py

@@ -0,0 +1,103 @@
+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
+_H_TO_KJ = 2625.4996394798254  # 1 Hartree in kJ/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
+    )
+    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():
+    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_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
+    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