Add a pKa helper for proton-coupled redox

ThermoScreening/thermo/__init__.py

@@ -11,3 +11,4 @@
 from .reactions import reaction_free_energy, reduction_potential
 from .ensemble import boltzmann_weights, ensemble_free_energy, lowest_gibbs
 from .kinetics import eyring_rate_constant, wigner_tunneling_correction
+from .pka import pKa, calibrate_proton_reference, PROTON_AQUEOUS_FREE_ENERGY_KCAL

ThermoScreening/thermo/pka.py

@@ -0,0 +1,121 @@
+"""pKa (acid dissociation) from computed ``Thermo`` objects.
+
+A pure post-processing helper combining the absolute Gibbs free energies
+(``Thermo.total_EeGtot()``) of an acid and its conjugate base into a pKa, via
+the "direct method" thermodynamic cycle ``HA(soln) -> A-(soln) + H+(soln)``.
+Compute both species with the same engine and conditions -- in solution (e.g.
+``solvent="water"``) -- then combine them here.
+"""
+
+import math
+
+from ..utils.physicalConstants import PhysicalConstants
+from ._units import HARTREE_TO_KCAL_PER_MOL
+
+# Absolute free energy of the proton in aqueous solution at 298.15 K, in
+# kcal/mol, referenced to the same 1 atm gas-phase-style standard state
+# ThermoScreening's own RRHO thermochemistry uses (so it adds directly to
+# Thermo.total_EeGtot() with no further standard-state correction):
+#
+#   G(H+, gas, 1 atm)      = -6.28 kcal/mol   (ideal monatomic gas, Sackur-
+#                                               Tetrode; Bartmess, J. Phys.
+#                                               Chem. 1994, 98, 6420)
+# + dG_solv(H+, aq, 1 atm) = -264.0 kcal/mol  (Tissandier et al., J. Phys.
+#                                               Chem. A 1998, 102, 7787, and
+#                                               Kelly, Cramer & Truhlar,
+#                                               J. Phys. Chem. B 2006, 110,
+#                                               16066, recommend -265.9
+#                                               kcal/mol at a 1 mol/L
+#                                               gas-phase reference state;
+#                                               converting to the 1 atm
+#                                               reference used throughout here
+#                                               subtracts RT ln(24.46) = 1.89
+#                                               kcal/mol at 298.15 K -- see
+#                                               their Section 2 and Eq. 1-3)
+# = G(H+, aq, 1 atm)       = -270.28 kcal/mol
+#
+# Convention-dependent like SHE_ABSOLUTE_POTENTIAL in reactions.py; override
+# reference_free_energy for a value calibrated to your method (see
+# calibrate_proton_reference), which the literature recommends for
+# quantitative accuracy.
+PROTON_AQUEOUS_FREE_ENERGY_KCAL = -270.28
+
+_R_KCAL_PER_MOL_K = PhysicalConstants["R"] / (PhysicalConstants["cal"] * 1000.0)
+
+
+def _delta_g_kcal(acid, base):
+    """G(base) - G(acid), in kcal/mol (Hartree per particle -> kcal/mol)."""
+    return (base.total_EeGtot() - acid.total_EeGtot()) * HARTREE_TO_KCAL_PER_MOL
+
+
+def pKa(acid, base, temperature=298.15, reference_free_energy=PROTON_AQUEOUS_FREE_ENERGY_KCAL):
+    """
+    pKa for ``HA(soln) -> A-(soln) + H+(soln)`` (the "direct method").
+
+    ``pKa = dG / (R T ln 10)``, with ``dG = G(base) - G(acid) +
+    reference_free_energy`` the deprotonation free energy including the
+    proton's aqueous reference free energy (the proton has no electronic
+    structure, so no engine can compute it directly).
+
+    Parameters
+    ----------
+    acid, base : Thermo
+        The acid (HA) and its conjugate base (A-), computed consistently (same
+        engine, same conditions -- ideally in solution, e.g. ``solvent="water"``
+        -- and open-shell/charged as appropriate).
+    temperature : float
+        Temperature in K. Should match the temperature ``acid``/``base`` were
+        computed at. Default 298.15.
+    reference_free_energy : float
+        The proton's absolute aqueous free energy in kcal/mol. Defaults to
+        :data:`PROTON_AQUEOUS_FREE_ENERGY_KCAL`; pass a value calibrated to
+        your method (see :func:`calibrate_proton_reference`) for quantitative
+        accuracy.
+
+    Returns
+    -------
+    float
+        The (dimensionless) pKa.
+
+    Notes
+    -----
+    The raw "direct method" with a literature proton reference is known to
+    have several-pKa-unit systematic error even at DFT+continuum-solvent
+    levels (e.g. Ho & Coote, Theor. Chem. Acc. 2010, 125, 3); GFN-xTB/DFTB
+    absolute pKa is expected to be considerably less accurate still. Use for
+    **relative** comparisons among structurally similar acids, or calibrate
+    ``reference_free_energy`` against one experimental pKa of a similar
+    reference acid.
+    """
+    delta_g_kcal = _delta_g_kcal(acid, base) + reference_free_energy
+    return delta_g_kcal / (_R_KCAL_PER_MOL_K * temperature * math.log(10))
+
+
+def calibrate_proton_reference(acid, base, experimental_pKa, temperature=298.15):
+    """
+    The ``reference_free_energy`` that reproduces ``experimental_pKa`` for a
+    reference acid/base pair.
+
+    Solves :func:`pKa` for ``reference_free_energy`` instead of ``pKa``. Use
+    the result as ``pKa(..., reference_free_energy=...)`` for other acids
+    computed the same way (same engine/conditions), ideally structurally
+    similar to the calibration pair -- the literature-recommended approach for
+    quantitative pKa accuracy, since the raw literature proton reference alone
+    is not quantitatively reliable (see :func:`pKa`'s Notes).
+
+    Parameters
+    ----------
+    acid, base : Thermo
+        A reference acid/base pair with a known experimental pKa.
+    experimental_pKa : float
+        The reference acid's experimental pKa.
+    temperature : float
+        Temperature in K, matching ``acid``/``base``. Default 298.15.
+
+    Returns
+    -------
+    float
+        The calibrated ``reference_free_energy`` (kcal/mol).
+    """
+    target_delta_g_kcal = experimental_pKa * _R_KCAL_PER_MOL_K * temperature * math.log(10)
+    return target_delta_g_kcal - _delta_g_kcal(acid, base)

docs/api.rst

@@ -33,6 +33,12 @@ Reactions and redox
 .. autofunction:: ThermoScreening.thermo.reactions.reaction_free_energy
 .. autofunction:: ThermoScreening.thermo.reactions.reduction_potential
 
+Acid dissociation (pKa)
+------------------------
+
+.. autofunction:: ThermoScreening.thermo.pka.pKa
+.. autofunction:: ThermoScreening.thermo.pka.calibrate_proton_reference
+
 Transition states and kinetics
 -------------------------------
 

docs/usage.rst

@@ -214,6 +214,48 @@ conditions), combine them into reaction free energies and reduction potentials:
    across similar species, with a higher-accuracy method, or with a
    ``reference_potential`` calibrated against experiment.
 
+Acid dissociation (pKa)
+------------------------
+
+Proton-coupled redox (e.g. hydroquinone/semiquinone/quinone protonation
+states) needs a pKa alongside the reduction potential. ``pKa`` uses the
+"direct method" thermodynamic cycle -- ``HA(soln) -> A-(soln) + H+(soln)`` --
+combining the acid and conjugate base's computed free energies with a
+literature reference free energy for the (uncomputable) aqueous proton:
+
+.. code-block:: python
+
+   from ThermoScreening.thermo.api import xtb_cli_thermo
+   from ThermoScreening.thermo.conformers import generate
+   from ThermoScreening.thermo import pKa
+
+   # the acid and its conjugate base are different structures (one fewer H),
+   # not the same structure at a different charge (that would be reduction)
+   hq       = generate("Oc1ccc(O)cc1", max_conformers=1)[0]     # hydroquinone, HQ
+   hq_anion = generate("[O-]c1ccc(O)cc1", max_conformers=1)[0]  # phenolate, HQ-
+
+   acid = xtb_cli_thermo(hq, charge=0, solvent="water")
+   base = xtb_cli_thermo(hq_anion, charge=-1, solvent="water")
+   p_ka = pKa(acid, base)
+
+.. warning::
+
+   The default proton reference (``PROTON_AQUEOUS_FREE_ENERGY_KCAL``) is a
+   literature constant; the raw direct method is known to have several-pKa-unit
+   systematic error even with DFT and an explicit continuum solvent model (Ho
+   & Coote, *Theor. Chem. Acc.* **2010**, *125*, 3), and GFN-xTB/DFTB absolute
+   pKa is expected to be considerably less accurate still. Use for **relative**
+   comparisons among structurally similar acids, or calibrate against one
+   known experimental pKa:
+
+   .. code-block:: python
+
+      from ThermoScreening.thermo import calibrate_proton_reference
+
+      # a reference acid/base pair with a known experimental pKa
+      ref_g = calibrate_proton_reference(ref_acid, ref_base, experimental_pKa=4.20)
+      p_ka = pKa(acid, base, reference_free_energy=ref_g)  # more accurate
+
 Transition states and rate constants
 -------------------------------------
 

tests/thermo/test_pka.py

@@ -0,0 +1,118 @@
+import math
+
+import pytest
+
+from ThermoScreening.thermo.pka import (
+    pKa,
+    calibrate_proton_reference,
+    PROTON_AQUEOUS_FREE_ENERGY_KCAL,
+)
+
+_H_TO_KCAL = 627.5094740631  # 1 Hartree in kcal/mol
+_R_KCAL = 8.314462618 / 4184.0  # kcal/(mol K)
+
+
+class _FakeThermo:
+    """A stand-in exposing only what the pKa helpers use."""
+
+    def __init__(self, eegtot):
+        self._eegtot = eegtot
+
+    def total_EeGtot(self):
+        return self._eegtot
+
+
+def test_default_proton_reference_matches_the_documented_derivation():
+    # G(H+, gas, 1 atm) = -6.28 kcal/mol (Bartmess 1994) + dG_solv(H+, aq,
+    # 1 atm) = -264.0 kcal/mol (Tissandier 1998 / Kelly-Cramer-Truhlar 2006,
+    # converted from their 1 mol/L reference state)
+    assert PROTON_AQUEOUS_FREE_ENERGY_KCAL == pytest.approx(-270.28, abs=0.01)
+
+
+def test_pKa_zero_delta_g_gives_pKa_from_reference_alone():
+    acid = _FakeThermo(-10.0)
+    base = _FakeThermo(-10.0)  # G(base) - G(acid) = 0
+
+    result = pKa(acid, base, reference_free_energy=0.0)
+
+    assert result == pytest.approx(0.0)
+
+
+def test_pKa_known_value_acetic_acid_scale():
+    # solve dG (Hartree) so that, with reference_free_energy=0, pKa lands at
+    # the experimental acetic acid value (4.756) -- a numerical sign/scale
+    # sanity check: a weak acid must give a POSITIVE pKa (unfavorable
+    # dissociation), not a negative one
+    temperature = 298.15
+    target_pKa = 4.756
+    dG_kcal = target_pKa * _R_KCAL * temperature * math.log(10)
+    dG_hartree = dG_kcal / _H_TO_KCAL
+
+    acid = _FakeThermo(-10.0)
+    base = _FakeThermo(-10.0 + dG_hartree)
+
+    result = pKa(acid, base, temperature=temperature, reference_free_energy=0.0)
+
+    assert result == pytest.approx(target_pKa)
+    assert result > 0  # weak acid -> positive pKa
+
+
+def test_pKa_stronger_acid_has_lower_pKa():
+    acid = _FakeThermo(-10.0)
+    weak_base = _FakeThermo(-9.98)    # larger dG -> weaker acid -> higher pKa
+    strong_base = _FakeThermo(-10.05)  # smaller (negative) dG -> stronger acid
+
+    weak_pKa = pKa(acid, weak_base, reference_free_energy=0.0)
+    strong_pKa = pKa(acid, strong_base, reference_free_energy=0.0)
+
+    assert strong_pKa < weak_pKa
+
+
+def test_pKa_reference_free_energy_shifts_result_linearly():
+    acid, base = _FakeThermo(-10.0), _FakeThermo(-10.0)
+    rt_ln10 = _R_KCAL * 298.15 * math.log(10)
+
+    p0 = pKa(acid, base, reference_free_energy=0.0)
+    p1 = pKa(acid, base, reference_free_energy=rt_ln10)  # +1 pKa unit worth
+
+    assert p1 - p0 == pytest.approx(1.0)
+
+
+def test_pKa_uses_default_reference_free_energy():
+    acid, base = _FakeThermo(-10.0), _FakeThermo(-10.0)
+
+    default_call = pKa(acid, base)
+    explicit_call = pKa(acid, base, reference_free_energy=PROTON_AQUEOUS_FREE_ENERGY_KCAL)
+
+    assert default_call == pytest.approx(explicit_call)
+
+
+def test_calibrate_proton_reference_round_trips_with_pKa():
+    acid = _FakeThermo(-10.0)
+    base = _FakeThermo(-9.99)
+    target_pKa = 7.2
+
+    calibrated = calibrate_proton_reference(acid, base, target_pKa, temperature=310.0)
+    recovered = pKa(acid, base, temperature=310.0, reference_free_energy=calibrated)
+
+    assert recovered == pytest.approx(target_pKa)
+
+
+def test_calibrate_proton_reference_matches_closed_form():
+    acid, base = _FakeThermo(-10.0), _FakeThermo(-9.995)
+    temperature = 298.15
+
+    calibrated = calibrate_proton_reference(acid, base, 5.0, temperature=temperature)
+
+    dG_kcal = (base.total_EeGtot() - acid.total_EeGtot()) * _H_TO_KCAL
+    expected = 5.0 * _R_KCAL * temperature * math.log(10) - dG_kcal
+    assert calibrated == pytest.approx(expected)
+
+
+def test_public_api_exported():
+    from ThermoScreening.thermo import pKa as pKa_pub
+    from ThermoScreening.thermo import calibrate_proton_reference as cal_pub
+    from ThermoScreening.thermo import pka as pka_module
+
+    assert pKa_pub is pka_module.pKa
+    assert cal_pub is pka_module.calibrate_proton_reference