Add a pKa helper for proton-coupled redox - #115
Merged
Merged
Conversation
Add ThermoScreening.thermo.pka: a pure post-processing helper combining the absolute Gibbs free energies of an acid and its conjugate base into a pKa, via the "direct method" thermodynamic cycle HA(soln) -> A-(soln) + H+(soln), pKa = dG / (R T ln 10). Mirrors reaction_free_energy/ reduction_potential's design: takes already-computed Thermo objects, raises no exceptions of its own (pure arithmetic). The proton has no electronic structure, so its aqueous free energy is a literature reference constant. Verified the default (-270.28 kcal/mol) against primary sources rather than memory: - G(H+, gas, 1 atm) = -6.28 kcal/mol (Sackur-Tetrode ideal monatomic gas; Bartmess, J. Phys. Chem. 1994, 98, 6420) - dG_solv(H+, aq) = -264.0 kcal/mol (Tissandier et al., J. Phys. Chem. A 1998, 102, 7787; Kelly, Cramer & Truhlar, J. Phys. Chem. B 2006, 110, 16066 recommend -265.9 kcal/mol at a 1 mol/L gas-phase reference state; converted here to the 1 atm reference ThermoScreening's own RRHO thermochemistry uses, subtracting RT ln(24.46) = 1.89 kcal/mol) Also add calibrate_proton_reference(acid, base, experimental_pKa), which solves for the reference_free_energy that reproduces a known experimental pKa -- the literature-recommended way to get quantitative accuracy, since the raw direct method has known several-pKa-unit systematic error even at DFT+continuum-solvent levels (documented prominently, matching the existing GFN2-redox-accuracy caveat in reactions.py). Docs: a new "Acid dissociation (pKa)" usage section (hydroquinone/HQ- worked example, matching the existing AQ-domain examples) and API reference entries. Verified by building the docs locally and visually checking the rendered warning admonitions/code blocks.
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #115 +/- ##
==========================================
+ Coverage 97.29% 97.31% +0.01%
==========================================
Files 31 32 +1
Lines 1921 1935 +14
==========================================
+ Hits 1869 1883 +14
Misses 52 52
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
…tures Review nit (and a chemistry bug beyond the reported undefined variable): the example built 'hq_anion' from the SAME hydroquinone structure at charge=-1, which represents reduction (adding an electron), not deprotonation. The acid and conjugate base differ by one H atom -- build them as distinct structures via SMILES (matching the existing generate(...) pattern used elsewhere in this doc), so HQ- genuinely has one fewer H than HQ.
galjos
added a commit
that referenced
this pull request
Jul 9, 2026
…117) **Found by running real xtb-cli calculations** on the pKa helper from #115, the same "test with real data" pattern used for the QM-import engines this session. ## The finding The docs said the raw direct-method pKa has "several-pKa-unit systematic error" (citing DFT+continuum-solvent literature) and suggested calibration for "quantitative accuracy." Running it for real on GFN2-xTB (phenol/phenolate, water) tells a very different story: the raw pKa comes out at **≈-100** against the experimental **9.99** — an error of **~110 pKa units**, not "several." This isn't a bug — I verified the arithmetic is internally self-consistent (the ~150 kcal/mol calibration shift needed exactly accounts for the gap). It's a real, well-known property of semiempirical tight-binding methods: they don't preserve an absolute ab-initio/experimental energy scale, so mixing their electronic energies directly with a literature constant derived assuming DFT/ab-initio-quality absolute energies is fundamentally invalid — not just imprecise. Calibrating against **one** reference compound completely fixes this (it's a molecule-independent additive offset): hydroquinone calibrated against phenol (pKa_exp = 9.99) came out at **11.00** vs. the experimental **10.35** — a good result for a cheap semiempirical method. ## The fix - `pka.py`: strengthen the `Notes` sections — calibration is **required**, not merely recommended, for GFN-xTB/DFTB, with the real verified numbers. - `docs/usage.rst`: replace the old placeholder calibration example (`calibrate_proton_reference(ref_acid, ref_base, experimental_pKa=4.20)`, using undefined variables — not even runnable) with this real, working, verified phenol/hydroquinone example. - Tests: a skippable real-xtb-cli integration test reproducing both numbers end to end, gated on the native `xtb` binary (matching `test_reactions.py`'s existing pattern). ## Verification Ran locally with the real `xtb` binary: 10/10 tests pass including the new real-engine test. Full suite: 391 passed, 13 skipped (correctly gated when `xtb` is absent). 100% patch coverage. Docs build clean. Closes #116
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Adds
ThermoScreening.thermo.pka, mirroringreaction_free_energy/reduction_potential's design (pure post-processing over already-computedThermoobjects), for proton-coupled equilibria like the hydroquinone/semiquinone/quinone protonation states relevant to AQ redox.The method
The "direct method" thermodynamic cycle:
HA(soln) -> A-(soln) + H+(soln),pKa = dG / (R T ln 10). The proton has no electronic structure (no engine can compute it), so its aqueous free energy is a literature reference constant.The default constant — verified against primary sources, not memory
PROTON_AQUEOUS_FREE_ENERGY_KCAL = -270.28kcal/mol, derived as:G(H+, gas, 1 atm) = -6.28 kcal/mol— the standard Sackur-Tetrode ideal-monatomic-gas value (Bartmess, J. Phys. Chem. 1994, 98, 6420)dG_solv(H+, aq) = -264.0 kcal/mol— Tissandier et al. (J. Phys. Chem. A 1998, 102, 7787) / Kelly, Cramer & Truhlar (J. Phys. Chem. B 2006, 110, 16066) recommend -265.9 kcal/mol at a 1 mol/L gas-phase reference state; I converted this to the 1 atm reference ThermoScreening's own RRHO thermochemistry actually uses (subtractingRT ln(24.46) = 1.89 kcal/mol), since mixing the two standard-state conventions is a well-documented, easy-to-make ~1.2-pKa-unit error in this literature (the papers explicitly warn about it).I verified every step numerically (
RT ln10 = 1.364 kcal/molat 298.15 K; the acetic-acid pKa=4.756 sign/magnitude check; the standard-state conversion arithmetic) before writing any code — see the module's inline derivation comment for the full chain.Calibration helper
calibrate_proton_reference(acid, base, experimental_pKa)solves for thereference_free_energythat reproduces one known experimental pKa — the literature-recommended path to quantitative accuracy, since the raw direct method has known several-pKa-unit systematic error even at DFT+continuum-solvent levels (documented prominently in the docstring/docs, matching the existing GFN2-redox-accuracy caveat pattern already inreactions.py).Verification
9 tests: the default-constant regression check, sign/scale sanity (weak acid → positive pKa, cross-checked against real acetic acid pKa), reference-energy linearity, the
calibrate_proton_reference/pKaround-trip, and public API export. Full suite passes (391 passed, 12 pre-existing skips), 100% patch coverage. Docs: new "Acid dissociation (pKa)" usage section + API reference entries, built locally and visually verified in the browser (rendered warning admonitions, code blocks, and the Sphinx search index correctly resolving both new functions).Closes #114