pKa of every site, tested against IUPAC data

Run on 2026-10-10. The forecast was written before the first number was compared. Models M3–M5p were added after the first test showed what M1 missed. They are labelled that way below and judged only on molecules they were not fitted to.

#The question

On /pka-density a molecule you type gets a pH slider. Until now each ionisable site titrated on its own, at a textbook group value: 4.8 for every carboxylic acid, 10.6 for every amine. Two questions:

• How wrong is that against measured values?
• Can measured values make the pH curves better without hiding what the model is?

#The data

IUPAC Digitized pKa Dataset v2.4a (J. W. Zheng and O. Lafontant-Joseph). It digitizes three IUPAC compilations:

• Serjeant and Dempsey, Ionisation Constants of Organic Acids in Aqueous Solution (1979);
• Perrin, Dissociation Constants of Organic Bases in Aqueous Solution (1965, supplement 1972).

DOI 10.5281/zenodo.7236452. Reproduced by permission of International Union of Pure and Applied Chemistry. Licence CC BY-NC 4.0.

The app ships only the fitted parameters, not the data. fit_site_model.py reads the CSV from a path you give it.

Selection.

• Values: water only, 18–26 °C, IUPAC assessment "Reliable" or "Approximate", A or AH dissociations.
• Molecules: read with the app's own SMILES reader, sites found by the app's own site typing.
• A molecule is used only when the number of measured steps equals the number of steps the model puts in water's range. That leaves 2832 molecules: 2163 with one step and 669 with several, 1417 steps in all.
• Amphoteric molecules split across the two books are merged, with values closer than 0.3 counted once. Glycine, for example, has its acid step in Serjeant and its base step in Perrin.

Test. Five-fold cross-validation by molecule. Every number below is for molecules the model was not fitted to.

#The model

A titration measures macroscopic steps: the moments the molecule as a whole loses a proton. Behind them are the 2ⁿ protonation microstates of n sites. Each microstate h (hᵢ = 1 when site i holds its proton) gets the weight

• qᵢ is the site's charge: −1 for a deprotonated acid, +1 for a protonated base, otherwise 0.
• One site alone gives back Henderson–Hasselbalch.
• The pair term is Bjerrum's 1923 explanation of why the second proton of a dicarboxylic acid comes off later than the first: the first negative charge holds the second proton back.
• Two equal sites therefore differ by W + log 4. The log 4 is pure counting: two ways to lose the first proton, two ways to regain it.

The pair term. W falls off with the number of bonds d between the two groups:

The site pKa. Each site starts from its parent compound's value (acetic acid 4.8, phenol 10.0, methylamine 10.6, pyridine 5.2 and so on). It is then shifted by its substituents:

This is the Perrin–Dempsey–Serjeant recipe (pKa Prediction for Organic Acids and Bases, 1981): a parent value plus substituent effects that fade by a constant factor per bond. Here the substituent strength is Pauling's electronegativity difference instead of a looked-up σ* constant. An sp² or sp carbon counts too, because it is more electronegative than sp³ (Walsh, Bent). The α-amino-acid kinds have no single parent. Their textbook values (2.3, 9.6) already contain the zwitterion's own pull, so their site pKa is fitted.

Nothing in this model is quantum mechanical. It is the empirical layer under the π maps. The π descriptors are tested separately (below).

#Forecast and outcome

#Results

All numbers are held-out RMSE in pKa units. M1c (pair term as B/d instead of A·ρ^(d−1)) is indistinguishable from M1: the data cannot tell the two shapes apart.

On the left, every textbook kind is a horizontal stripe: one value for hundreds of different molecules. On the right the stripes have tilted toward the diagonal. What is left is the spread the model cannot see.

Why M5p and not M5. The rule written before the second iteration picks the lowest pooled error, which is M5, 0.035 better. But M5 lets each kind's value drift toward the average compound of the dataset: acetic acid 4.39 instead of 4.75, methylamine 10.0 instead of 10.67. Those parents are the molecules a visitor types first. M5p fixes them to their measured values and loses almost nothing held out. We chose it after seeing these numbers; that is a judgement, and it is stated here.

Textbook molecules, held out (measured / today / installed):

Over 30 textbook values the RMSE is 0.95 today and 0.46 installed. The largest remaining miss is 1.2 units (the second step of the dicarboxylic acids).

The pair term, seen directly. For molecules with two identical groups, pK₂ − pK₁ − log 4 is W.

• Amine pairs sit above the single curve and carboxylate pairs below it, as forecast. The ammonium charge sits on one atom; the carboxylate charge is shared by two oxygens.
• The d = 2 points are far above the curve. Two acids on the same carbon (substituted malonic acids) hold their monoanion together with an internal hydrogen bond, which a charge–charge term cannot describe.
• One amplitude per charge pair (M2, M4) captures the split but wins too little held out to justify three amplitudes.

Stability. Across the five folds A = 2.93 ± 0.10, ρ = 0.779 ± 0.011, λ = 0.705 ± 0.011.

Per kind (one-site molecules, held-out RMSE, today → installed): carboxylic acids 1.36 → 0.85, aromatic carboxylic acids 1.09 → 0.70, phenols 3.11 → 1.36, amines 2.08 → 1.18, α-amino groups 2.35 → 0.73, imidazoles 2.22 → 1.97. Anilines (2.37) and pyridine-type N (2.25) did not improve: their ring substituents act through resonance, which a through-bond sum does not see. For these two kinds the substituent slope is switched off by a stated rule (a kind keeps its slope only if it beats today's value held out), so they use their parent value. The page shows each site with its kind's held-out error.

#What it cannot do

• Resonance through rings. Anilines, pyridines and heterocycles depend on where a substituent sits (ortho/para against meta), and on π acceptors and donors. That is what the π descriptors were built for.
• Shape. Bond counts are not distances. Internal hydrogen bonds, folded chains and ring strain are invisible.
• Tautomers. Each site is typed as drawn.
• Conditions. The values mix 18–26 °C and different ionic strengths; amines move about 0.02 per °C. The fit cannot be more precise than this spread.
• Groups the typing does not know. Pyrazoles, triazoles, oxazoles, sulfonamides, hydroxamic acids and others get no site.

#The π lines: the next test

The π forecasts (Q_site, V_probe) need a 3D conformer, and this experiment ran without RDKit. On the server:

nice -n 19 python3 tools/iupac/pi_line_benchmark.py --csv iupac_high-confidence_v2_4.csv --jobs 1

It asks, held out, three questions:

1. How good are today's lines (fitted on 12 + 5 molecules) on hundreds of molecules?
2. Refitted per kind, which descriptor works best?
3. Does a descriptor explain what this site model leaves over? Anilines and pyridines are where it should.

If the answer to 3 is yes for a kind, the π descriptor becomes that kind's substituent term.

#What changed in the app

• /pka-density, any molecule, solution mode:
  • All sites (up to six) titrate together.
  • The maps are the population-weighted mixture over every site.
  • Below the map: the fraction of each charge state from pH 0 to 14, every deprotonation step with the sites that take part, and each site's charge at the chosen pH.
• /api/pi-model returns a speciation block: site pKa, ±, pair terms, macroscopic steps.
• /api/model's protonation uses the same model.
• The site typing no longer calls pyrrole and indole N–H a base. It recognises tetrazoles as acids, and imidazole only as a 1,3-diazole.
• The app's own SMILES reader now reads aromatic SMILES.

#Reproduce

python tools/iupac/fit_site_model.py --csv iupac_high-confidence_v2_4.csv --out iupac_out --write --select M5p
python tools/iupac/plot_results.py iupac_out

Takes about 4 minutes on one core and needs NumPy, SciPy and Matplotlib, not RDKit. Per-molecule predictions go to iupac_out/, outside the package.

#Origins

• L. J. Henderson, Am. J. Physiol. 21, 173 (1908)
• K. A. Hasselbalch, Biochem. Z. 78, 112 (1917)
• N. Bjerrum, Z. Phys. Chem. 106, 219 (1923): the charge between two acid groups
• J. G. Kirkwood and F. H. Westheimer, J. Chem. Phys. 6, 506 (1938): why the effective dielectric between two groups is lower than water's
• L. Pauling, J. Am. Chem. Soc. 54, 3570 (1932): electronegativity
• D. D. Perrin, B. Dempsey and E. P. Serjeant, pKa Prediction for Organic Acids and Bases (Chapman & Hall, 1981): parent value + fall-off per bond
• Microstates and macroscopic constants: the binding-polynomial treatment of Wyman (1964) and of Ullmann for proteins (2003)
• The data: Zheng and Lafontant-Joseph (2026), from Serjeant and Dempsey (1979) and Perrin (1965, 1972), reproduced by permission of IUPAC.