Forecast: written 2026-10-10 20:50 CEST, before any IUPAC number was compared with the model
Data: IUPAC Digitized pKa Dataset v2.4a (Zheng & Lafontant-Joseph; digitized from Serjeant & Dempsey 1979 and Perrin 1965/1972). Reproduced by permission of IUPAC; CC BY-NC 4.0. The data are not shipped with the app; fit_site_model.py reads the CSV from a path you give it.
Question: how good are the pH curves of /pka-density for an arbitrary molecule, and can measured values make them better? Today each ionisable site gets a group-average pKa from a textbook table (qo_chem/protonation.py) and titrates on its own (Henderson–Hasselbalch per site).
Two models are compared on the same molecules:
- M0, today's model: the table pKa, independent sites.
- M1, the new model: one intrinsic pKa per site kind, plus a charge–charge interaction between sites. The interaction falls off with the number of bonds d between the two groups: W = A·ρ^(d−1) pKa units. All 2^n protonation microstates are summed.
- This is Bjerrum's 1923 picture of dicarboxylic acids, with the distance taken through bonds. The fall-off per bond is the factor Perrin, Dempsey and Serjeant use for substituent effects (1981).
- The fit is cross-validated by molecule, five folds.
The comparison uses the sorted macroscopic pKa values: each step at which the molecule as a whole loses a proton. This is what a titration measures and what the dataset reports. A molecule counts only when its number of measured values equals the number of steps the model puts between −1 and 15.
Forecasts
F1. Monoprotic molecules, today's table. RMSE about 1.4 pKa units, median absolute error about 0.8. The table is right on average; the error is the spread inside each kind. Phenols run from 7 to 11, pyridines from 1 to 7, and amines are cut by fluorine, hydroxyl and ether neighbours.
F2. Refitting the table on IUPAC medians lowers the monoprotic RMSE by at most 10%, except for kinds whose textbook value is not the median of the substituted compounds actually measured. Pyridine-type N and aniline should both come out lower than 5.2 and 4.6, because the dataset is full of electron-poor heterocycles and substituted anilines.
F3. Polyprotic molecules (two or more measured steps). M0 makes equal groups titrate on top of each other: succinic acid gets 4.5 and 5.1, measured 4.2 and 5.6. Today's RMSE on the second and later steps is above 1.0. M1 lowers it by 30–50%.
F4. Interaction parameters.
- Fall-off per bond ρ between 0.4 and 0.6; A between 2 and 4 pKa units.
- Ammonium pairs interact more strongly than carboxylate pairs at the same d, because the charge sits on one atom instead of being spread over two oxygens. The residuals will show this as a systematic sign between pair types. A single A will not remove it.
F5. The intrinsic pKa of the α-amino-acid COOH (both neighbours neutral) comes out near 3.5–4. That is between acetic acid (4.76) and the measured 2.3. The difference is the ammonium's charge, which M1 now puts in the interaction instead of in the table value.
F6. What M1 cannot fix: substituent effects inside a kind (nitrophenols, chloroacetic acids, pyrimidines). Those need molecule-specific descriptors, which is the job of the π forecast lines. They are tested separately, on the server's own conformers, by pi_line_benchmark.py, because the descriptors need 3D coordinates and this experiment was run without RDKit.
Acceptance rule, fixed now: the page uses M1 only if its cross-validated RMSE on held-out molecules beats M0 on both the monoprotic and the polyprotic set. A parameter is reported with its fold-to-fold spread, not only its best value.