Clarify that calibration is required (not just recommended) for pKa

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)
@@ -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
     ----------

docs/usage.rst

@@ -227,34 +227,47 @@ 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
+   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; 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, ``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
+   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
 -------------------------------------

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