diff --git a/ThermoScreening/thermo/__init__.py b/ThermoScreening/thermo/__init__.py index e6a6f89..a615c61 100644 --- a/ThermoScreening/thermo/__init__.py +++ b/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 diff --git a/ThermoScreening/thermo/pka.py b/ThermoScreening/thermo/pka.py new file mode 100644 index 0000000..33a371d --- /dev/null +++ b/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) diff --git a/docs/api.rst b/docs/api.rst index 8140c1b..200a1ac 100644 --- a/docs/api.rst +++ b/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 ------------------------------- diff --git a/docs/usage.rst b/docs/usage.rst index 31e4b26..6a2b800 100644 --- a/docs/usage.rst +++ b/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 ------------------------------------- diff --git a/tests/thermo/test_pka.py b/tests/thermo/test_pka.py new file mode 100644 index 0000000..8c99ab5 --- /dev/null +++ b/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