Add reference-calibrated redox potentials

README.md

@@ -23,6 +23,7 @@ python -m sphinx -b html docs docs/_build/html   # open docs/_build/html/index.h
 
 - Thermochemistry calculations for molecular systems
 - DFTB+ geometry optimization, Hessian, and normal-mode integration
+- Stepwise and overall reference-calibrated reduction potentials
 - Readers for DFTB+ `.gen`, XYZ, and vibrational frequency files
 - Runtime type checking for public API calls
 - Test coverage for parsing, thermochemistry, and optional DFTB+ execution paths

ThermoScreening/thermo/__init__.py

@@ -8,7 +8,11 @@
 from .thermo import Thermo
 from .screening import screen, rank_by_gibbs
 from .conformers import generate as generate_conformers, write_conformers, generate_thermo_ensemble
-from .reactions import reaction_free_energy, reduction_potential
+from .reactions import (
+    calibrate_reduction_reference,
+    reaction_free_energy,
+    reduction_potential,
+)
 from .ensemble import boltzmann_weights, ensemble_free_energy, lowest_gibbs, EnsembleThermo
 from .kinetics import eyring_rate_constant, wigner_tunneling_correction
 from .pka import pKa, calibrate_proton_reference, PROTON_AQUEOUS_FREE_ENERGY_KCAL

ThermoScreening/thermo/reactions.py

@@ -78,7 +78,7 @@ def reduction_potential(
         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.
+        Positive 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
@@ -92,11 +92,49 @@ def reduction_potential(
     Raises
     ------
     ValueError
-        If ``n_electrons`` is zero.
+        If ``n_electrons`` is not positive.
     """
-    if n_electrons == 0:
-        raise ValueError("n_electrons must be non-zero.")
+    if n_electrons <= 0:
+        raise ValueError("n_electrons must be positive.")
 
     delta_g_hartree = reduced.total_EeGtot() - oxidized.total_EeGtot()
     absolute = -delta_g_hartree * _HARTREE_TO_EV / n_electrons
     return absolute - reference_potential
+
+
+def calibrate_reduction_reference(
+    oxidized, reduced, experimental_potential, n_electrons=1
+):
+    """
+    Calibrate a reduction-potential reference against one known redox pair.
+
+    The returned value can be passed as ``reference_potential`` to
+    :func:`reduction_potential` for other compounds computed with the same
+    method and conditions. This assumes that the systematic offset between the
+    computed absolute energy scale and the experimental reference electrode is
+    transferable from the reference pair to the target compounds.
+
+    Parameters
+    ----------
+    oxidized, reduced : Thermo
+        The reference redox pair, computed consistently with the target
+        compounds.
+    experimental_potential : float
+        Experimental reduction potential of the reference pair, in volts.
+    n_electrons : int
+        Positive number of electrons transferred. Default 1.
+
+    Returns
+    -------
+    float
+        The calibrated value for ``reference_potential``, in volts.
+
+    Raises
+    ------
+    ValueError
+        If ``n_electrons`` is not positive.
+    """
+    absolute = reduction_potential(
+        oxidized, reduced, n_electrons=n_electrons, reference_potential=0.0
+    )
+    return absolute - experimental_potential

docs/api.rst

@@ -34,6 +34,7 @@ Reactions and redox
 
 .. autofunction:: ThermoScreening.thermo.reactions.reaction_free_energy
 .. autofunction:: ThermoScreening.thermo.reactions.reduction_potential
+.. autofunction:: ThermoScreening.thermo.reactions.calibrate_reduction_reference
 
 Acid dissociation (pKa)
 ------------------------

docs/usage.rst

@@ -101,8 +101,9 @@ The ensemble free energy lies at or below the lowest conformer's by the mixing
 (conformational) entropy; use ``lowest_gibbs`` when you instead want the single
 dominant structure to carry forward (e.g. into :func:`reaction_free_energy`).
 
-``pKa``, ``calibrate_proton_reference``, ``reduction_potential``, and
-``reaction_free_energy`` all take whatever they're given and call
+``pKa``, ``calibrate_proton_reference``, ``reduction_potential``,
+``calibrate_reduction_reference``, and ``reaction_free_energy`` all take
+whatever they're given and call
 ``.total_EeGtot()`` on it -- they don't care whether that's a single ``Thermo``
 or something else with the same method. ``EnsembleThermo`` wraps a conformer
 ensemble's Boltzmann free energy behind that same method, so a flexible acid
@@ -250,14 +251,95 @@ conditions), combine them into reaction free energies and reduction potentials:
    E = reduction_potential(neutral, anion)          # vs SHE (default reference 4.44 V)
    E_abs = reduction_potential(neutral, anion, reference_potential=0.0)
 
+For a screening series, calibrate the computed energy scale against a
+structurally similar reference compound measured on the same potential scale
+and under the same conditions. Anthraquinone has two successive one-electron
+reductions, so each step gets its own calibration:
+
+.. code-block:: python
+
+   from ThermoScreening.thermo import (
+       calibrate_reduction_reference,
+       reduction_potential,
+   )
+
+   e1_reference = calibrate_reduction_reference(
+       aq_neutral,
+       aq_radical_anion,
+       experimental_potential=measured_e1,
+   )
+   e2_reference = calibrate_reduction_reference(
+       aq_radical_anion,
+       aq_dianion,
+       experimental_potential=measured_e2,
+   )
+
+   candidate_e1 = reduction_potential(
+       candidate_neutral,
+       candidate_radical_anion,
+       reference_potential=e1_reference,
+   )
+   candidate_e2 = reduction_potential(
+       candidate_radical_anion,
+       candidate_dianion,
+       reference_potential=e2_reference,
+   )
+
+Here ``measured_e1`` and ``measured_e2`` are the experimental anthraquinone
+formal potentials in volts versus the chosen reference electrode. They describe
+different reactions and must not be interchanged:
+
+.. math::
+
+   \mathrm{AQ + e^- \rightarrow AQ^-} \qquad E_1
+
+.. math::
+
+   \mathrm{AQ^- + e^- \rightarrow AQ^{2-}} \qquad E_2
+
+The overall two-electron potential describes a third reaction. For two
+successive one-electron reductions under identical conditions, it is the
+arithmetic mean of the stepwise formal potentials:
+
+.. math::
+
+   \mathrm{AQ + 2e^- \rightarrow AQ^{2-}} \qquad
+   E_{2e} = \frac{E_1 + E_2}{2}
+
+Calculate and calibrate it from the neutral species directly to the dianion:
+
+.. code-block:: python
+
+   measured_e2e = (measured_e1 + measured_e2) / 2
+   e2e_reference = calibrate_reduction_reference(
+       aq_neutral,
+       aq_dianion,
+       experimental_potential=measured_e2e,
+       n_electrons=2,
+   )
+   candidate_e2e = reduction_potential(
+       candidate_neutral,
+       candidate_dianion,
+       n_electrons=2,
+       reference_potential=e2e_reference,
+   )
+
+Do not compare ``candidate_e2e`` with the experimental second reduction
+potential ``measured_e2``. The averaging relationship requires equilibrium or
+formal potentials; it should not be applied directly to irreversible or
+kinetically distorted peak potentials.
+
 .. 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.
+   ``reference_potential`` calibrated against experiment. Calibration assumes
+   that the reference compound's systematic offset transfers to the target
+   compounds; keep the engine, solvent model, temperature, and charge-state
+   treatment identical across the calibration and target set.
 
 Reactivity site prediction (Fukui indices)
 -------------------------------------------

tests/thermo/test_reactions.py

@@ -6,6 +6,7 @@
 
 from ThermoScreening.thermo.reactions import (
     SHE_ABSOLUTE_POTENTIAL,
+    calibrate_reduction_reference,
     reaction_free_energy,
     reduction_potential,
 )
@@ -68,17 +69,124 @@ def test_reduction_potential_n_electrons():
     )
 
 
-def test_reduction_potential_rejects_zero_electrons():
+@pytest.mark.parametrize("n_electrons", [0, -1])
+def test_reduction_potential_rejects_non_positive_electrons(n_electrons):
     ox, red = _FakeThermo(0.0), _FakeThermo(-0.5)
-    with pytest.raises(ValueError, match="non-zero"):
-        reduction_potential(ox, red, n_electrons=0)
+    with pytest.raises(ValueError, match="positive"):
+        reduction_potential(ox, red, n_electrons=n_electrons)
+
+
+@pytest.mark.parametrize("n_electrons", [1, 2])
+def test_calibrate_reduction_reference_reproduces_experiment(n_electrons):
+    ox, red = _FakeThermo(0.0), _FakeThermo(-0.2)
+    experimental = -0.75
+
+    reference = calibrate_reduction_reference(
+        ox, red, experimental, n_electrons=n_electrons
+    )
+
+    assert reduction_potential(
+        ox,
+        red,
+        n_electrons=n_electrons,
+        reference_potential=reference,
+    ) == pytest.approx(experimental)
+
+
+def test_calibrated_reference_transfers_relative_potential():
+    reference_ox, reference_red = _FakeThermo(0.0), _FakeThermo(-0.2)
+    target_ox, target_red = _FakeThermo(-1.0), _FakeThermo(-1.3)
+    experimental = -0.75
+
+    reference = calibrate_reduction_reference(
+        reference_ox, reference_red, experimental
+    )
+    target = reduction_potential(
+        target_ox, target_red, reference_potential=reference
+    )
+
+    assert target == pytest.approx(experimental + 0.1 * _H_TO_EV)
+
+
+def test_two_step_reduction_uses_one_calibration_per_step():
+    ref_neutral = _FakeThermo(-100.0)
+    ref_anion = _FakeThermo(-100.13)
+    ref_dianion = _FakeThermo(-100.23)
+    target_neutral = _FakeThermo(-200.0)
+    target_anion = _FakeThermo(-200.14)
+    target_dianion = _FakeThermo(-200.25)
+    experimental_e1 = -0.75
+    experimental_e2 = -1.4
+
+    e1_reference = calibrate_reduction_reference(
+        ref_neutral, ref_anion, experimental_e1
+    )
+    e2_reference = calibrate_reduction_reference(
+        ref_anion, ref_dianion, experimental_e2
+    )
+    e2e_reference = calibrate_reduction_reference(
+        ref_neutral,
+        ref_dianion,
+        (experimental_e1 + experimental_e2) / 2,
+        n_electrons=2,
+    )
+
+    target_e1 = reduction_potential(
+        target_neutral, target_anion, reference_potential=e1_reference
+    )
+    target_e2 = reduction_potential(
+        target_anion, target_dianion, reference_potential=e2_reference
+    )
+    target_e2e = reduction_potential(
+        target_neutral,
+        target_dianion,
+        n_electrons=2,
+        reference_potential=e2e_reference,
+    )
+
+    assert target_e1 == pytest.approx(experimental_e1 + 0.01 * _H_TO_EV)
+    assert target_e2 == pytest.approx(experimental_e2 + 0.01 * _H_TO_EV)
+    assert e2e_reference == pytest.approx((e1_reference + e2_reference) / 2)
+    assert target_e2e == pytest.approx((target_e1 + target_e2) / 2)
+
+
+def test_overall_two_electron_potential_is_mean_of_stepwise_potentials():
+    neutral = _FakeThermo(-100.0)
+    radical_anion = _FakeThermo(-100.13)
+    dianion = _FakeThermo(-100.23)
+    reference = 4.44
+
+    first = reduction_potential(
+        neutral, radical_anion, reference_potential=reference
+    )
+    second = reduction_potential(
+        radical_anion, dianion, reference_potential=reference
+    )
+    overall = reduction_potential(
+        neutral, dianion, n_electrons=2, reference_potential=reference
+    )
+
+    assert overall == pytest.approx((first + second) / 2)
+
+
+@pytest.mark.parametrize("n_electrons", [0, -1])
+def test_calibrate_reduction_reference_rejects_non_positive_electrons(
+    n_electrons,
+):
+    ox, red = _FakeThermo(0.0), _FakeThermo(-0.5)
+    with pytest.raises(ValueError, match="positive"):
+        calibrate_reduction_reference(
+            ox, red, -0.75, n_electrons=n_electrons
+        )
 
 
 def test_public_api_exported():
+    from ThermoScreening.thermo import calibrate_reduction_reference as crr
     from ThermoScreening.thermo import reaction_free_energy as rfe
     from ThermoScreening.thermo import reduction_potential as rp
     from ThermoScreening.thermo import reactions
 
+    assert crr is reactions.calibrate_reduction_reference
     assert rfe is reactions.reaction_free_energy
     assert rp is reactions.reduction_potential