From e8c22366a874c71a326b68bd532d660f71049453 Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Thu, 9 Jul 2026 11:23:38 +0200 Subject: [PATCH 1/2] Clarify that calibration is required (not just recommended) for pKa Verified with real xtb-cli (GFN2-xTB + ALPB water) calculations, not just the closed-form unit tests: phenol's raw pKa with the default literature proton reference comes out at approximately -100 against the experimental 9.99 -- an error of ~110 pKa units, not the "several units" DFT+ continuum-solvent methods show. This isn't a bug (the arithmetic is internally consistent, verified by back-computing the calibration shift); semiempirical tight-binding methods don't preserve an absolute ab-initio/experimental energy scale, so their energies can't be mixed directly with a literature constant derived assuming DFT/ab-initio-quality absolute energies. Calibrating against one reference compound absorbs the offset completely: hydroquinone calibrated against phenol (pKa_exp=9.99) comes out at 11.00 vs. the experimental 10.35. - pka.py: strengthen the Notes to state calibration is required (not just recommended) for GFN-xTB/DFTB, with the verified real numbers. - docs/usage.rst: replace the placeholder calibration example (calibrate_proton_reference(ref_acid, ref_base, experimental_pKa=4.20), which used undefined variables) with this real, working, verified phenol/hydroquinone example. - tests: a skippable real-xtb-cli integration test reproducing both results end to end (gated on the native xtb binary, matching the existing pattern in test_reactions.py). Closes #116 --- ThermoScreening/thermo/pka.py | 33 ++++++++++++++++------- docs/usage.rst | 50 ++++++++++++++++++++++------------- tests/thermo/test_pka.py | 37 ++++++++++++++++++++++++++ 3 files changed, 91 insertions(+), 29 deletions(-) diff --git a/ThermoScreening/thermo/pka.py b/ThermoScreening/thermo/pka.py index 33a371d..ba4cce7 100644 --- a/ThermoScreening/thermo/pka.py +++ b/ThermoScreening/thermo/pka.py @@ -79,13 +79,25 @@ def pKa(acid, base, temperature=298.15, reference_free_energy=PROTON_AQUEOUS_FRE 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. + The raw "direct method" with a literature proton reference has + several-pKa-unit systematic error even at DFT+continuum-solvent levels + (e.g. Ho & Coote, Theor. Chem. Acc. 2010, 125, 3). With a semiempirical + engine (GFN-xTB, DFTB) it is far worse than that -- **not usable at all** + uncalibrated: verified with real xtb-cli (GFN2-xTB + ALPB water) + calculations, phenol's raw pKa comes out at approximately -100 (vs. the + experimental 9.99), an error of about 110 pKa units. This is not a + solvation-model inaccuracy; semiempirical tight-binding methods do not + preserve an absolute ab-initio/experimental energy scale, so their + electronic energies cannot be combined directly with a literature + constant derived assuming DFT/ab-initio-quality absolute energies. + + ``calibrate_proton_reference`` against **one** reference compound absorbs + this offset entirely (it is an additive constant, not molecule-specific) + and is **required**, not merely recommended, for a semiempirical engine: + in the same verification, hydroquinone's pKa calibrated against phenol + (pKa_exp = 9.99) came out at 11.00 vs. the experimental 10.35 -- a good + result. See the "Acid dissociation (pKa)" section of the usage docs for + the full worked example. """ delta_g_kcal = _delta_g_kcal(acid, base) + reference_free_energy return delta_g_kcal / (_R_KCAL_PER_MOL_K * temperature * math.log(10)) @@ -99,9 +111,10 @@ def calibrate_proton_reference(acid, base, experimental_pKa, temperature=298.15) 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). + similar to the calibration pair. With a semiempirical engine (GFN-xTB, + DFTB) this is **required**, not just recommended for extra accuracy: the + raw literature proton reference alone gives an unusable result (errors of + ~100 pKa units, not just a few -- see :func:`pKa`'s Notes). Parameters ---------- diff --git a/docs/usage.rst b/docs/usage.rst index 6a2b800..186faaf 100644 --- a/docs/usage.rst +++ b/docs/usage.rst @@ -227,34 +227,46 @@ literature reference free energy for the (uncomputable) aqueous proton: from ThermoScreening.thermo.api import xtb_cli_thermo from ThermoScreening.thermo.conformers import generate - from ThermoScreening.thermo import pKa + from ThermoScreening.thermo import pKa, calibrate_proton_reference - # 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- + def acid_base_pair(acid_smiles, base_smiles): + # the acid and its conjugate base are different structures (one fewer + # H), not the same structure at a different charge (that would be + # reduction, not deprotonation) + acid = generate(acid_smiles, max_conformers=1)[0] + base = generate(base_smiles, max_conformers=1)[0] + return ( + xtb_cli_thermo(acid, charge=0, solvent="water"), + xtb_cli_thermo(base, charge=-1, solvent="water"), + ) - acid = xtb_cli_thermo(hq, charge=0, solvent="water") - base = xtb_cli_thermo(hq_anion, charge=-1, solvent="water") - p_ka = pKa(acid, base) + phenol, phenolate = acid_base_pair("Oc1ccccc1", "[O-]c1ccccc1") + hq, hq_anion = acid_base_pair("Oc1ccc(O)cc1", "[O-]c1ccc(O)cc1") # hydroquinone + + p_ka = pKa(hq, hq_anion) # -99.45 -- see the warning below .. 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: + literature constant derived assuming DFT/ab-initio-quality absolute + energies. With **GFN-xTB or DFTB it is not usable at all uncalibrated**: + verified with the real xtb-cli calculation above, ``pKa(phenol, + phenolate)`` comes out at **-100.46** against phenol's experimental 9.99 -- + an error of ~110 pKa units, not the "several units" DFT+continuum-solvent + methods show (Ho & Coote, *Theor. Chem. Acc.* **2010**, *125*, 3). + Semiempirical tight-binding methods do not preserve an absolute + ab-initio/experimental energy scale, so their energies cannot be mixed + directly with this constant. + + ``calibrate_proton_reference`` against **one** reference compound absorbs + this offset completely (it is an additive constant, not + molecule-specific) and is **required** for these engines, not just + recommended for extra accuracy: .. 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 + ref_g = calibrate_proton_reference(phenol, phenolate, experimental_pKa=9.99) + p_ka = pKa(hq, hq_anion, reference_free_energy=ref_g) # 11.00 vs. exp 10.35 Transition states and rate constants ------------------------------------- diff --git a/tests/thermo/test_pka.py b/tests/thermo/test_pka.py index 8c99ab5..0be727c 100644 --- a/tests/thermo/test_pka.py +++ b/tests/thermo/test_pka.py @@ -1,4 +1,6 @@ import math +import os +import shutil import pytest @@ -116,3 +118,38 @@ def test_public_api_exported(): assert pKa_pub is pka_module.pKa assert cal_pub is pka_module.calibrate_proton_reference + + +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_pKa_calibration_end_to_end(tmp_path): + # phenol/phenolate and hydroquinone/hydroquinone-monoanion, in water, via + # xtb-cli (GFN2-xTB + ALPB). The raw pKa with the default literature + # proton reference is unusable for a semiempirical engine (~-100, since + # GFN2-xTB energies are not on the ab-initio/experimental absolute energy + # scale that reference assumes) -- calibrating against phenol's + # experimental pKa (9.99) fixes this: hydroquinone comes out within ~1 + # pKa unit of its experimental value (10.35). + from ThermoScreening.thermo.api import xtb_cli_thermo + from ThermoScreening.thermo.conformers import generate + + def acid_base_pair(acid_smiles, base_smiles, tag): + acid = generate(acid_smiles, max_conformers=1)[0] + base = generate(base_smiles, max_conformers=1)[0] + return ( + xtb_cli_thermo(acid, charge=0, solvent="water", directory=str(tmp_path / f"{tag}_a")), + xtb_cli_thermo(base, charge=-1, solvent="water", directory=str(tmp_path / f"{tag}_b")), + ) + + phenol, phenolate = acid_base_pair("Oc1ccccc1", "[O-]c1ccccc1", "phenol") + hq, hq_anion = acid_base_pair("Oc1ccc(O)cc1", "[O-]c1ccc(O)cc1", "hq") + + raw_pKa = pKa(phenol, phenolate) + assert math.isfinite(raw_pKa) + assert raw_pKa < -50 # unusable uncalibrated, as documented + + ref_g = calibrate_proton_reference(phenol, phenolate, experimental_pKa=9.99) + calibrated_pKa = pKa(hq, hq_anion, reference_free_energy=ref_g) + assert calibrated_pKa == pytest.approx(10.35, abs=1.0) # within ~1 unit of experiment From 6eb78104f8dbc6efbd4ae8a8611c1e7e6d254610 Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Thu, 9 Jul 2026 11:30:01 +0200 Subject: [PATCH 2/2] Address review nits: consistent required-not-recommended wording, show phenol pKa in docs example --- ThermoScreening/thermo/pka.py | 4 ++-- docs/usage.rst | 5 +++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/ThermoScreening/thermo/pka.py b/ThermoScreening/thermo/pka.py index ba4cce7..2977cc0 100644 --- a/ThermoScreening/thermo/pka.py +++ b/ThermoScreening/thermo/pka.py @@ -36,8 +36,8 @@ # # 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. +# calibrate_proton_reference) -- required, not just recommended, for +# semiempirical engines (GFN-xTB, DFTB); see pKa's Notes below. PROTON_AQUEOUS_FREE_ENERGY_KCAL = -270.28 _R_KCAL_PER_MOL_K = PhysicalConstants["R"] / (PhysicalConstants["cal"] * 1000.0) diff --git a/docs/usage.rst b/docs/usage.rst index 186faaf..b47fc9e 100644 --- a/docs/usage.rst +++ b/docs/usage.rst @@ -244,14 +244,15 @@ literature reference free energy for the (uncomputable) aqueous proton: hq, hq_anion = acid_base_pair("Oc1ccc(O)cc1", "[O-]c1ccc(O)cc1") # hydroquinone p_ka = pKa(hq, hq_anion) # -99.45 -- see the warning below + phenol_p_ka = pKa(phenol, phenolate) # -100.46 vs. experimental 9.99 .. warning:: The default proton reference (``PROTON_AQUEOUS_FREE_ENERGY_KCAL``) is a literature constant derived assuming DFT/ab-initio-quality absolute energies. With **GFN-xTB or DFTB it is not usable at all uncalibrated**: - verified with the real xtb-cli calculation above, ``pKa(phenol, - phenolate)`` comes out at **-100.46** against phenol's experimental 9.99 -- + verified with the real xtb-cli calculation above, ``phenol_p_ka`` comes out + at **-100.46** against phenol's experimental 9.99 -- an error of ~110 pKa units, not the "several units" DFT+continuum-solvent methods show (Ho & Coote, *Theor. Chem. Acc.* **2010**, *125*, 3). Semiempirical tight-binding methods do not preserve an absolute