From b4882a10e9c8665d37331be8ef1fb935a8bd4f31 Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Wed, 22 Jul 2026 13:42:36 +0200 Subject: [PATCH] Add reference-calibrated redox potentials --- README.md | 1 + ThermoScreening/thermo/__init__.py | 6 +- ThermoScreening/thermo/reactions.py | 46 ++++++++++- docs/api.rst | 1 + docs/usage.rst | 88 ++++++++++++++++++++- tests/thermo/test_reactions.py | 114 +++++++++++++++++++++++++++- 6 files changed, 245 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index e6d7a52..82a6af4 100644 --- a/README.md +++ b/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 diff --git a/ThermoScreening/thermo/__init__.py b/ThermoScreening/thermo/__init__.py index dd5ef4f..1130f92 100644 --- a/ThermoScreening/thermo/__init__.py +++ b/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 diff --git a/ThermoScreening/thermo/reactions.py b/ThermoScreening/thermo/reactions.py index 7987c72..929090c 100644 --- a/ThermoScreening/thermo/reactions.py +++ b/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 diff --git a/docs/api.rst b/docs/api.rst index 202cf44..905226f 100644 --- a/docs/api.rst +++ b/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) ------------------------ diff --git a/docs/usage.rst b/docs/usage.rst index 5cdf15e..eb58215 100644 --- a/docs/usage.rst +++ b/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,6 +251,84 @@ 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 @@ -257,7 +336,10 @@ conditions), combine them into reaction free energies and reduction potentials: 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) ------------------------------------------- diff --git a/tests/thermo/test_reactions.py b/tests/thermo/test_reactions.py index d8b168f..403ba9d 100644 --- a/tests/thermo/test_reactions.py +++ b/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