#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Where the textbook titration shortcuts actually fail. Science Journaling Club, Volume 1 Issue 3, Spring 2025. THE QUESTION ------------ General chemistry teaches several closed-form shortcuts for acid-base titration curves and buffer pH. Each arrives with a hand-wave about when it stops working ("as long as the acid is weak", "provided the solution is not too dilute"). Nobody says how wrong they get, or where the boundary actually sits. This script measures both. THE MODEL --------- An exact numerical solver for the aqueous acid-base titration of a protic acid H_nA with strong base, built from charge balance and mass balance with no approximations anywhere. Let h = [H+], and let the acid have stepwise dissociation constants Ka_1 ... Ka_n with cumulative products beta_k = prod_{j<=k} Ka_j (beta_0 = 1). The fraction of total acid carrying k removed protons is alpha_k(h) = beta_k h^(n-k) / D(h), D(h) = sum_{k=0..n} beta_k h^(n-k) and the mean number of protons removed per acid molecule is nbar(h) = sum_k k alpha_k(h). With C_A the total (diluted) analytical acid concentration and C_B the concentration of added strong-base cation, charge balance reads C_B + h = Kw/h + C_A nbar(h) (*) Mass balance on the acid is enforced identically by the construction of alpha_k, since sum_k alpha_k = 1 follows from the definition of D(h); the script checks that numerically anyway. The left side of (*) is strictly increasing in h and the right side strictly decreasing, so (*) has exactly one positive root. We find it by bisection in pH on [-3, 17] with 90 halvings, which brackets the root to 2e-26 pH units, far below the resolution of a double. Titrant concentration equals the initial analyte concentration unless stated otherwise, so dilution is real and is carried through C_A and C_B at every point of every curve. THE APPROXIMATIONS UNDER TEST ----------------------------- HH Henderson-Hasselbalch, pH = pKa + log10(phi/(1-phi)), the ratio taken from stoichiometry alone (phi = fraction of the proton titrated). SQRT the initial-pH formula [H+] = sqrt(Ka C), which assumes the acid dissociates negligibly and that water supplies no protons. QUAD the same with dissociation kept, h^2 + Ka h - Ka C = 0, water still ignored. Isolates the two error sources from one another. NOW the full speciation solved with Kw set to zero: water autoionisation dropped, everything else exact. EQV the equivalence-point formula pOH = (pKb - log10 C')/2. AMPH the amphiprotic-salt shortcut pH = (pK1 + pK2)/2. ASSUMPTIONS ----------- * 25 degrees C throughout, Kw = 1.000e-14. * Ideal solution. All activity coefficients are 1, so concentrations stand in for activities everywhere. This is the largest single departure from a real titration and it is not small: at 0.1 M ionic strength the Davies equation puts the activity coefficient of a univalent ion near 0.78, worth roughly 0.11 in pH. Every number below is the pH of an ideal solution, which is the same idealisation the textbook formulas are written in, so the comparison between exact and approximate is fair even though the comparison to a real glass electrode is not. * Strong base is fully dissociated; the acid is the only weak species. * A closed vessel. No atmospheric CO2 dissolves in during the titration. * Stepwise Ka values are independent constants, the usual convention, which hides statistical-factor effects between identical protonation sites. LIMITATIONS ----------- * This is a computation, not an experiment. No solution was mixed, no burette was read, no electrode was calibrated. Every number printed below is output of the equations above, and its agreement with a real titration is limited by the assumptions listed, chiefly activity. * Kw, and the Ka values used for the named real acids, come from published compilations at 25 C. The solver does not derive them. * Junction potentials, electrode drift, carbonate error in the titrant, temperature variation and slow kinetics are all outside the model. * The error of an approximation is reported against this exact solver, not against a measured pH. Those are different quantities. VALIDATION ---------- Eight checks, each printing the club's computed value beside the analytic or accepted value with the difference. Charge-balance and mass-balance residuals are printed at machine precision. Random seed: 20250411, numpy default_rng (PCG64). Fully deterministic. """ import sys import time from decimal import Decimal, getcontext import numpy as np SEED = 20250411 KW = 1.0e-14 PKW = 14.0 LN10 = np.log(10.0) EPS = np.finfo(float).eps T0 = time.time() SOLVE_COUNT = [0] # ---------------------------------------------------------------------------- # core speciation # ---------------------------------------------------------------------------- def nbar(h, Kas): """Mean protons removed per acid molecule at proton concentration h. Kas is a sequence of stepwise dissociation constants, each a scalar or an array broadcastable against h. Evaluated in the h^-k normalisation, which keeps every partial product inside double range over pH -3 to 17. """ h = np.asarray(h, dtype=float) num = np.zeros_like(h) den = np.ones_like(h) beta = np.ones_like(h) for k, Ka in enumerate(Kas, start=1): beta = beta * (Ka / h) den = den + beta num = num + k * beta return num / den def alphas(h, Kas): """All speciation fractions alpha_0 .. alpha_n at proton concentration h.""" h = np.asarray(h, dtype=float) betas = [np.ones_like(h)] for Ka in Kas: betas.append(betas[-1] * (Ka / h)) den = np.sum(betas, axis=0) return np.array([b / den for b in betas]) def charge_residual(h, CA, CB, Kas, Kw=KW): """Signed residual of the charge-balance equation (*), in mol/L.""" return CB + h - Kw / h - CA * nbar(h, Kas) def _bisect_ph(f, shape, iters, lo, hi): """Shared bisection driver. f(pH) must decrease through its single root.""" lo = np.full(shape, float(lo)) hi = np.full(shape, float(hi)) for _ in range(iters): mid = 0.5 * (lo + hi) pos = f(mid) > 0.0 hi = np.where(pos, hi, mid) lo = np.where(pos, mid, lo) return 0.5 * (lo + hi) def solve_ph(CA, CB, Kas, Kw=KW, iters=90, lo=-3.0, hi=17.0): """Exact pH by bisection on the strictly monotone charge balance.""" CA, CB = np.broadcast_arrays(np.asarray(CA, float), np.asarray(CB, float)) SOLVE_COUNT[0] += int(np.asarray(CA).size) def f(ph): return charge_residual(10.0 ** (-ph), CA, CB, Kas, Kw) return _bisect_ph(f, CA.shape, iters, lo, hi) def solve_ph_strong(CA, CB, Kw=KW, iters=90, lo=-3.0, hi=17.0): """Same driver with nbar identically 1: a fully dissociated strong acid.""" CA, CB = np.broadcast_arrays(np.asarray(CA, float), np.asarray(CB, float)) SOLVE_COUNT[0] += int(np.asarray(CA).size) def f(ph): h = 10.0 ** (-ph) return CB + h - Kw / h - CA return _bisect_ph(f, CA.shape, iters, lo, hi) def solve_ph_nowater(CA, CB, Kas, iters=120, lo=-3.0, hi=40.0): """Same solver with water autoionisation removed (Kw = 0 exactly).""" CA, CB = np.broadcast_arrays(np.asarray(CA, float), np.asarray(CB, float)) SOLVE_COUNT[0] += int(np.asarray(CA).size) def f(ph): h = 10.0 ** (-ph) return CB + h - CA * nbar(h, Kas) return _bisect_ph(f, CA.shape, iters, lo, hi) def titration_conc(Ca, phi, ratio=1.0): """Diluted C_A and C_B at titrated fraction phi. Titrant concentration is ratio * Ca. phi counts moles of OH- delivered per mole of the proton being titrated, so equivalence is phi = 1. """ phi = np.asarray(phi, dtype=float) v = phi / ratio # titrant volume per unit initial volume CA = Ca / (1.0 + v) CB = Ca * ratio * v / (1.0 + v) return CA, CB # ---------------------------------------------------------------------------- # the taught approximations # ---------------------------------------------------------------------------- def ph_hh(pKa, phi): """Henderson-Hasselbalch from stoichiometry alone.""" phi = np.asarray(phi, dtype=float) return pKa + np.log10(phi / (1.0 - phi)) def ph_sqrt(pKa, C): """[H+] = sqrt(Ka C): negligible dissociation, no water.""" return 0.5 * (np.asarray(pKa, float) - np.log10(np.asarray(C, float))) def ph_quad(pKa, C): """h^2 + Ka h - Ka C = 0: dissociation kept, water still ignored.""" Ka = 10.0 ** (-np.asarray(pKa, dtype=float)) C = np.asarray(C, dtype=float) h = 0.5 * (-Ka + np.sqrt(Ka * Ka + 4.0 * Ka * C)) return -np.log10(h) def ph_equiv(pKa, Cprime): """Equivalence point: pOH = (pKb - log10 C')/2 with pKb = 14 - pKa.""" pKb = PKW - np.asarray(pKa, dtype=float) pOH = 0.5 * (pKb - np.log10(np.asarray(Cprime, float))) return PKW - pOH def davies_gamma(I, z=1, A=0.5092): """Davies-equation activity coefficient for an ion of charge z. log10 gamma = -A z^2 ( sqrt(I)/(1+sqrt(I)) - 0.3 I ), the form de Levie uses for activity-corrected titration simulation. A = 0.5092 at 25 C. """ I = np.maximum(np.asarray(I, dtype=float), 0.0) sI = np.sqrt(I) return 10.0 ** (-A * z * z * (sI / (1.0 + sI) - 0.3 * I)) def solve_ph_davies(CA, CB, pKa, Kw=KW, outer=25): """Activity-corrected pH of a MONOPROTIC acid titration, Davies equation. HA is neutral and A- is univalent, so the conditional (concentration) constants are Ka' = Ka / gamma^2 and Kw' = Kw / gamma^2. For a monoprotic acid the ionic strength closes exactly on the charge balance: I = C_B + [H+]. Iterated to self-consistency. Returns the activity pH, -log10(gamma [H+]), which is what a calibrated electrode reads, and the activity coefficient itself. """ CA, CB = np.broadcast_arrays(np.asarray(CA, float), np.asarray(CB, float)) Ka = 10.0 ** (-np.asarray(pKa, float)) I = np.maximum(CB, 1e-10) + 1e-10 g = davies_gamma(I) h = np.zeros_like(I) for _ in range(outer): g = davies_gamma(I) pHc = solve_ph(CA, CB, [Ka / (g * g)], Kw=Kw / (g * g)) h = 10.0 ** (-pHc) I = CB + h # pH is minus log of the ACTIVITY of H+, which is gamma times [H+], # so it sits above the concentration pH by -log10(gamma), not below it. return -np.log10(g * h), g, I # ---------------------------------------------------------------------------- # independent cross-checks # ---------------------------------------------------------------------------- def ph_polyroot(CA, CB, Kas, Kw=KW, scale=1e-7): """Independent root of the same charge balance, as a polynomial in h. Multiplying (*) by h D(h) gives (h^2 + C_B h - Kw) D(h) - C_A h N(h) = 0, N(h) = sum_k k beta_k h^(n-k) a polynomial of degree n+2. Solved with numpy.roots, a companion-matrix eigenvalue method that shares no code path with bisection. Substituting h = scale * u first keeps the coefficients within a few orders of one another, which is the whole battle with a polynomial spanning Kw. """ n = len(Kas) betas = [1.0] for Ka in Kas: betas.append(betas[-1] * float(Ka)) D = np.array(betas, dtype=float) N = np.array([k * b for k, b in enumerate(betas)], dtype=float) P1 = np.array([1.0, float(CB), -Kw]) poly = np.polysub(np.convolve(P1, D), float(CA) * np.convolve([1.0, 0.0], N)) deg = len(poly) - 1 poly = poly * (scale ** np.arange(deg, -1, -1)) r = np.roots(poly) cand = [z.real * scale for z in r if abs(z.imag) < 1e-6 * max(1.0, abs(z.real)) and z.real > 1e-12] if not cand: return float("nan") best = min(cand, key=lambda hh: abs(charge_residual(hh, CA, CB, Kas, Kw))) return -np.log10(best) def ph_decimal(CA, CB, Kas, Kw=KW, digits=60, iters=220): """Bisection of the identical equation in 60-digit decimal arithmetic. Nothing here shares a code path with the float solver except the algebra, so agreement near 1e-15 pH is evidence the float answer is limited only by the double itself. """ getcontext().prec = digits CA = Decimal(repr(float(CA))) CB = Decimal(repr(float(CB))) Kw = Decimal(repr(float(Kw))) K = [Decimal(repr(float(k))) for k in Kas] ten = Decimal(10) def f(pH): h = ten ** (-pH) num = Decimal(0) den = Decimal(1) beta = Decimal(1) for i, Ka in enumerate(K, start=1): beta = beta * Ka / h den += beta num += i * beta return CB + h - Kw / h - CA * (num / den) lo, hi = Decimal(-3), Decimal(17) for _ in range(iters): mid = (lo + hi) / 2 if f(mid) > 0: lo = mid else: hi = mid return (lo + hi) / 2 def crossing(y, x, thresh, from_left): """First crossing of thresh scanning from one end, linear between nodes.""" idx = range(len(x) - 1) if from_left else range(len(x) - 2, -1, -1) for i in idx: a, b = y[i], y[i + 1] if (a > thresh) != (b > thresh): t = (thresh - a) / (b - a) return x[i] + t * (x[i + 1] - x[i]) return float("nan") def banner(s): print() print("=" * 78) print(s) print("=" * 78) def sub(s): print() print("-- " + s) # ============================================================================ def main(): rng = np.random.default_rng(SEED) banner("TITRATION APPROXIMATIONS: WHERE THE TEXTBOOK SHORTCUTS FAIL") print("Science Journaling Club, Volume 1 Issue 3, Spring 2025") print("python : %s" % sys.version.split()[0]) print("numpy : %s" % np.__version__) print("master seed : %d (numpy default_rng, PCG64)" % SEED) print("Kw : %.6e (25 C, ideal solution)" % KW) print("solver : bisection on charge balance, 90 halvings, pH in [-3,17]") print("final bracket: 20 / 2**90 = %.3e pH units" % (20.0 / 2.0 ** 90)) print("titrant : same concentration as analyte unless stated; dilution carried") print() print("The computation is the experiment. No solution was mixed and no") print("electrode was read. All activity coefficients are 1.") # ====================================================================== banner("PART 1. VALIDATION") # ====================================================================== sub("V1. Mass balance on the acid: do the speciation fractions sum to 1?") hgrid = 10.0 ** (-np.linspace(-3, 17, 4001)) worst_mass = 0.0 for name, Kas in [ ("monoprotic pKa 4.76", [10 ** -4.76]), ("diprotic 1.25 / 4.27", [10 ** -1.25, 10 ** -4.27]), ("triprotic 2.15/7.20/12.38", [10 ** -2.148, 10 ** -7.198, 10 ** -12.375]), ]: A = alphas(hgrid, Kas) d = float(np.max(np.abs(A.sum(axis=0) - 1.0))) worst_mass = max(worst_mass, d) print(" %-28s max |sum(alpha) - 1| = %.3e over 4001 pH values" % (name, d)) print(" double-precision epsilon = %.3e" % EPS) print(" worst residual is %.1f eps." % (worst_mass / EPS)) sub("V2. Charge balance residual across full titration curves") resid_rel_worst = 0.0 npts = 0 for pKa in [1.0, 2.5, 4.76, 7.0, 9.25, 11.0, 13.0]: for C in [1.0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5]: phi = np.linspace(0.0, 2.0, 401) CA, CB = titration_conc(C, phi) Kas = [10.0 ** -pKa] pH = solve_ph(CA, CB, Kas) h = 10.0 ** (-pH) r = charge_residual(h, CA, CB, Kas) scale = np.maximum.reduce([np.abs(CB), h, KW / h, CA * nbar(h, Kas)]) resid_rel_worst = max(resid_rel_worst, float(np.max(np.abs(r) / scale))) npts += phi.size print(" %d solved points over 42 (pKa, C) combinations, phi 0 to 2" % npts) print(" worst |charge residual| / largest term = %.3e" % resid_rel_worst) print(" that is %.1f eps. Charge balance holds to machine precision." % (resid_rel_worst / EPS)) sub("V3. Strong acid with strong base: closed form, exact arithmetic") print(" For a fully dissociated acid nbar = 1 and (*) collapses to") print(" h - Kw/h = C_A - C_B, whose positive root is") print(" h = (D + sqrt(D^2 + 4Kw))/2 with D = C_A - C_B. The club solver") print(" runs its ordinary bisection driver on the same case.") print() print(" Past equivalence D is negative and that expression subtracts two") print(" nearly equal numbers, so we also evaluate the algebraically") print(" identical h = 2Kw / (sqrt(D^2 + 4Kw) - D), which does not. Both") print(" are printed. The disagreement between them is the closed form's.") print() print(" %-9s %-7s %-18s %-18s %-11s %s" % ("C (M)", "phi", "club solver pH", "closed form pH", "naive-form", "difference")) v3worst = 0.0 v3naive = 0.0 for C in [1e-1, 1e-2, 1e-3, 1e-5]: for phi in [0.0, 0.5, 0.9, 0.999, 1.0, 1.001, 1.5]: CA, CB = titration_conc(C, phi) D = CA - CB root = np.sqrt(D * D + 4.0 * KW) h_stable = (0.5 * (D + root)) if D >= 0 else (2.0 * KW / (root - D)) h_naive = 0.5 * (D + root) ph_exact = -np.log10(h_stable) ph_naive = -np.log10(h_naive) ph_club = float(solve_ph_strong(np.array([CA]), np.array([CB]))[0]) d = ph_club - float(ph_exact) v3worst = max(v3worst, abs(d)) v3naive = max(v3naive, abs(ph_club - float(ph_naive))) print(" %-9.0e %-7.3f %-18.13f %-18.13f %-11.3e %+.3e" % (C, phi, ph_club, ph_exact, ph_club - ph_naive, d)) print(" worst absolute difference over 28 points: %.3e pH units" % v3worst) print(" worst against the cancelling form instead : %.3e pH units" % v3naive) print() print(" The residual left over is at phi exactly 1, where C_B = C_A and the") print(" sum C_B + h - Kw/h - C_A cancels two numbers of order 0.05 to reach") print(" one of order 1e-7. The float residual is identically zero over a") print(" window about 1e-11 pH wide, so nothing finer can be resolved there") print(" by any method that evaluates this expression in doubles. That is a") print(" real limit of the exact solver and it is eleven orders of magnitude") print(" below anything a chemist measures.") print() print(" The bisection agrees with the stable algebra to the last bit and") print(" beats the textbook quadratic by seven orders of magnitude past") print(" equivalence. We found this the wrong way round and had to go") print(" looking for a bug in the solver that was not there.") sub("V4. Pure water, no acid at all: pH must be exactly 7.000000") ph_w = float(solve_ph(np.array([0.0]), np.array([0.0]), [10.0 ** -4.76])[0]) print(" club solver : %.15f" % ph_w) print(" analytic -log10(sqrt(Kw)) : %.15f" % (-np.log10(np.sqrt(KW)))) print(" difference : %+.3e pH units" % (ph_w + np.log10(np.sqrt(KW)))) sub("V5. Bisection against an eigenvalue polynomial root, same equation") print(" %-26s %-7s %-6s %-17s %-17s %s" % ("acid", "C", "phi", "bisection pH", "numpy.roots pH", "difference")) v5worst = 0.0 cases = [ ("monoprotic pKa 4.76", [10 ** -4.76]), ("monoprotic pKa 1.00", [10 ** -1.0]), ("monoprotic pKa 11.0", [10 ** -11.0]), ("oxalic 1.25 / 4.27", [10 ** -1.25, 10 ** -4.27]), ("carbonic 6.35 / 10.33", [10 ** -6.35, 10 ** -10.33]), ("phosphoric 2.15/7.2/12.4", [10 ** -2.148, 10 ** -7.198, 10 ** -12.375]), ] for name, Kas in cases: for C, phi in [(0.1, 0.3), (0.01, 0.85), (1e-3, 0.5)]: CA, CB = titration_conc(C, phi) pb = float(solve_ph(np.array([CA]), np.array([CB]), Kas)[0]) pr = ph_polyroot(float(CA), float(CB), Kas) v5worst = max(v5worst, abs(pb - pr)) print(" %-26s %-7.0e %-6.2f %-17.12f %-17.12f %+.3e" % (name, C, phi, pb, pr, pb - pr)) print(" worst absolute difference over 18 points: %.3e pH units" % v5worst) print(" The companion-matrix root is the limiting factor here, not the") print(" bisection: the polynomial's coefficients span Kw.") sub("V6. Bisection against 60-digit decimal arithmetic, same equation") print(" %-26s %-7s %-6s %-22s %s" % ("acid", "C", "phi", "60-digit pH", "float - decimal")) v6worst = 0.0 for name, Kas in cases: for C, phi in [(0.1, 0.3), (1e-4, 0.5)]: CA, CB = titration_conc(C, phi) pb = float(solve_ph(np.array([CA]), np.array([CB]), Kas)[0]) pd = ph_decimal(float(CA), float(CB), Kas) d = pb - float(pd) v6worst = max(v6worst, abs(d)) print(" %-26s %-7.0e %-6.2f %-22s %+.3e" % (name, C, phi, ("%.18f" % pd), d)) print(" worst absolute difference over 12 points: %.3e pH units" % v6worst) sub("V7. Half equivalence of a weak acid: the textbook says pH = pKa") print(" This is the check that does NOT come out exact, and that is not a") print(" bug in the solver. Charge balance at phi = 0.5 gives") print(" [A-] = C_B + [H+] - [OH-], so [A-]/[HA] equals 1 only when") print(" [H+] - [OH-] is negligible beside C_B. To first order,") print(" pH - pKa = 4([H+] - [OH-]) / (C_A ln10).") print() print(" %-7s %-9s %-17s %-8s %-13s %-13s" % ("pKa", "C (M)", "exact pH", "pKa", "difference", "1st-order est")) for pKa, C in [(4.76, 1.0), (4.76, 0.1), (4.76, 1e-3), (4.76, 1e-5), (2.00, 0.1), (1.00, 0.1), (0.50, 0.1), (10.0, 0.1), (11.0, 1e-3), (12.0, 1e-4), (13.0, 1e-3)]: CA, CB = titration_conc(C, 0.5) pH = float(solve_ph(np.array([CA]), np.array([CB]), [10.0 ** -pKa])[0]) h = 10.0 ** (-pH) est = 4.0 * (h - KW / h) / (CA * LN10) print(" %-7.2f %-9.0e %-17.10f %-8.2f %+-13.6f %+-13.6f" % (pKa, C, pH, pKa, pH - pKa, est)) print() print(" The exact pH sits ABOVE pKa for strong acids, because the acid has") print(" already given up protons on its own, and BELOW pKa for very weak") print(" ones, because the conjugate base is pulling protons off water.") print(" The first-order estimate tracks the measured offset, which is the") print(" point: the deviation belongs to the approximation, not the solver.") sub("V8. The first-order half-equivalence rule against the solver, 675 cases") pk_v8 = np.linspace(0.3, 13.7, 135) for C in [1.0, 1e-1, 1e-2, 1e-3, 1e-4]: CA, CB = titration_conc(C, 0.5) pH = solve_ph(np.full(pk_v8.shape, CA), np.full(pk_v8.shape, CB), [10.0 ** -pk_v8]) h = 10.0 ** (-pH) est = 4.0 * (h - KW / h) / (CA * LN10) act = pH - pk_v8 small = np.abs(act) < 0.05 m1 = float(np.max(np.abs(est[small] - act[small]))) if small.any() else float("nan") print(" C = %-8.0e rule vs exact: max error where |offset| < 0.05" " = %.2e ; overall max = %.3f" % (C, m1, float(np.max(np.abs(est - act))))) print(" A linearisation is excellent where the offset is small and only") print(" indicative where it is large. That is exactly what it does here,") print(" which confirms the algebra behind the whole error analysis.") # ====================================================================== banner("PART 2. HENDERSON-HASSELBALCH ACROSS THE PARAMETER SPACE") # ====================================================================== pk_grid = np.arange(0.0, 14.0001, 0.25) # 57 values lc_grid = np.arange(-5.0, 0.0001, 0.125) # 41 values phi_grid = np.linspace(0.10, 0.90, 81) PK, LC, PHI = np.meshgrid(pk_grid, lc_grid, phi_grid, indexing="ij") Cs = 10.0 ** LC CA3, CB3 = titration_conc(Cs, PHI) pH_exact = solve_ph(CA3, CB3, [10.0 ** -PK]) pH_hh = PK + np.log10(PHI / (1.0 - PHI)) err_hh = pH_hh - pH_exact aerr = np.abs(err_hh) maxerr = aerr.max(axis=2) argphi = phi_grid[aerr.argmax(axis=2)] n_cells = maxerr.size print("grid: %d pKa values (0 to 14 by 0.25) x %d concentrations" % (len(pk_grid), len(lc_grid))) print(" (1e-5 to 1 M, 8 per decade) x %d titrated fractions (0.10 to 0.90)" % len(phi_grid)) print("exact solutions in this map: %d" % aerr.size) print() print("fraction of the %d (pKa, C) cells whose worst buffer-window HH error" % n_cells) print("exceeds:") for t in [0.01, 0.02, 0.05, 0.10, 0.20, 0.50, 1.00]: print(" %.2f pH units : %6.2f %% (%d cells)" % (t, 100.0 * np.mean(maxerr > t), int((maxerr > t).sum()))) iw, jw = np.unravel_index(int(maxerr.argmax()), maxerr.shape) print() print("worst cell on the grid: pKa = %.2f, C = %.3e M, worst error = %.3f pH" % (pk_grid[iw], 10.0 ** lc_grid[jw], maxerr.max())) print(" and the worst point inside that cell is at phi = %.3f" % argphi[iw, jw]) sub("Where in the buffer window the worst error lands") print(" cells whose worst point is at low phi (acid end) : %d" % int((argphi < 0.5).sum())) print(" cells whose worst point is at phi = 0.5 exactly : %d" % int((np.abs(argphi - 0.5) < 1e-9).sum())) print(" cells whose worst point is at high phi (base end) : %d" % int((argphi > 0.5).sum())) print(" HH is worst at the EDGES of the buffer window, never in the middle.") sub("HH error at C = 0.1 M, the concentration of a teaching titration") j01 = int(np.argmin(np.abs(lc_grid + 1.0))) k1 = int(np.argmin(np.abs(phi_grid - 0.1))) k5 = int(np.argmin(np.abs(phi_grid - 0.5))) k9 = int(np.argmin(np.abs(phi_grid - 0.9))) print(" %-7s %-13s %-13s %-13s %-13s" % ("pKa", "err phi=0.1", "err phi=0.5", "err phi=0.9", "worst |err|")) for pkv in [0.0, 1.0, 2.0, 3.0, 3.75, 4.0, 4.75, 6.0, 8.0, 9.75, 10.0, 11.0, 12.0, 13.0, 14.0]: i = int(np.argmin(np.abs(pk_grid - pkv))) print(" %-7.2f %+-13.4f %+-13.4f %+-13.4f %-13.4f" % (pk_grid[i], err_hh[i, j01, k1], err_hh[i, j01, k5], err_hh[i, j01, k9], maxerr[i, j01])) sub("The same sweep at C = 1e-4 M, a thousandfold more dilute") j4 = int(np.argmin(np.abs(lc_grid + 4.0))) print(" %-7s %-13s %-13s %-13s %-13s" % ("pKa", "err phi=0.1", "err phi=0.5", "err phi=0.9", "worst |err|")) for pkv in [0.0, 2.0, 4.0, 5.0, 6.0, 6.75, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0]: i = int(np.argmin(np.abs(pk_grid - pkv))) print(" %-7.2f %+-13.4f %+-13.4f %+-13.4f %-13.4f" % (pk_grid[i], err_hh[i, j4, k1], err_hh[i, j4, k5], err_hh[i, j4, k9], maxerr[i, j4])) # ====================================================================== banner("PART 3. THE BOUNDARY: WHERE HH CROSSES 0.05 AND 0.10 pH UNITS") # ====================================================================== print("For each concentration, HH's worst buffer-window error is large at low") print("pKa (the acid is too strong and dissociates on its own), small in the") print("middle, and large again at high pKa (the conjugate base deprotonates") print("water). The usable band lies between the two.") print() print("%-11s %-11s %-11s %-11s %-11s %-9s %-9s" % ("C (M)", "pKa lo .05", "pKa hi .05", "pKa lo .10", "pKa hi .10", "width .05", "width .10")) band_rows = [] for j, lc in enumerate(lc_grid): col = maxerr[:, j] lo05 = crossing(col, pk_grid, 0.05, True) hi05 = crossing(col, pk_grid, 0.05, False) lo10 = crossing(col, pk_grid, 0.10, True) hi10 = crossing(col, pk_grid, 0.10, False) band_rows.append((10.0 ** lc, lo05, hi05, lo10, hi10)) if abs(lc * 2 - round(lc * 2)) < 1e-9: print("%-11.3e %-11.3f %-11.3f %-11.3f %-11.3f %-9.3f %-9.3f" % (10.0 ** lc, lo05, hi05, lo10, hi10, hi05 - lo05, hi10 - lo10)) sub("A closed-form prediction of the same two boundaries") print(" Expanding both logarithms to first order gives") print(" HH error = -([H+] - [OH-]) / (C_A ln10 phi (1-phi))") print(" and at the acid end [H+] = Ka (1-phi)/phi to zeroth order, so the") print(" error there is Ka / (C_A ln10 phi^2). Setting that to a threshold T") print(" at phi = 0.1 gives pKa_lo = -log10(T C_A ln10 phi^2).") print(" At the base end [OH-] = Kw phi / (Ka (1-phi)) and the same step gives") print(" pKa_hi = 14 + log10(T C_A ln10 (1-phi)^2) at phi = 0.9.") print() print(" %-10s %-6s %-14s %-14s %-11s %-14s %-14s %s" % ("C (M)", "T", "solver lo", "closed form", "diff", "solver hi", "closed form", "diff")) for T in [0.05, 0.10]: for lc in [0.0, -1.0, -2.0, -3.0, -4.0, -5.0]: j = int(np.argmin(np.abs(lc_grid - lc))) C = 10.0 ** lc CA_lo = C / 1.1 # phi = 0.1 CA_hi = C / 1.9 # phi = 0.9 pred_lo = -np.log10(T * CA_lo * LN10 * 0.1 ** 2) pred_hi = PKW + np.log10(T * CA_hi * LN10 * 0.1 ** 2) slo = band_rows[j][1] if T == 0.05 else band_rows[j][3] shi = band_rows[j][2] if T == 0.05 else band_rows[j][4] print(" %-10.0e %-6.2f %-14.3f %-14.3f %+-11.3f %-14.3f %-14.3f %+.3f" % (C, T, slo, pred_lo, slo - pred_lo, shi, pred_hi, shi - pred_hi)) print(" The two agree to a few hundredths of a pKa unit across five decades") print(" of concentration, which is as much as a first-order expansion is") print(" entitled to. The boundary is not empirical; it is algebra.") sub("Where the band closes entirely") print(" Both edges move one pKa unit per decade of concentration, in") print(" opposite directions, so the band has to shut at some dilution.") have = [(np.log10(r[0]), r[1], r[2], r[3], r[4]) for r in band_rows if np.isfinite(r[1]) and np.isfinite(r[2])] for T, ilo, ihi in [(0.05, 1, 2), (0.10, 3, 4)]: xs = np.array([r[0] for r in band_rows if np.isfinite(r[ilo]) and np.isfinite(r[ihi])]) ys_lo = np.array([r[ilo] for r in band_rows if np.isfinite(r[ilo]) and np.isfinite(r[ihi])]) ys_hi = np.array([r[ihi] for r in band_rows if np.isfinite(r[ilo]) and np.isfinite(r[ihi])]) xs = np.log10(xs) slo, blo = np.polyfit(xs, ys_lo, 1) shi, bhi = np.polyfit(xs, ys_hi, 1) xstar = (bhi - blo) / (slo - shi) print(" T = %.2f pH: lo edge = %+0.4f log10C %+0.4f, hi edge = %+0.4f log10C %+0.4f" % (T, slo, blo, shi, bhi)) print(" they meet at log10 C = %.3f, C = %.3e M, pKa = %.3f" % (xstar, 10.0 ** xstar, blo + slo * xstar)) print(" below that concentration NO pKa keeps HH inside %.2f pH" % T) print(" The fitted slopes are %+.4f and %+.4f against an exact prediction of" % (slo, shi)) print(" -1 and +1 from the closed form above.") sub("The band at other tolerances") print(" The tolerance enters inside a logarithm, so changing it slides both") print(" edges by the same amount and moves the closing point. All five") print(" tolerances, located on the same grid:") print() print(" %-8s %-13s %-13s %-13s %-13s %s" % ("T (pH)", "lo at 0.1 M", "hi at 0.1 M", "width at 0.1M", "closes at C", "closing pKa")) for T in [0.01, 0.02, 0.05, 0.10, 0.20]: lo_t, hi_t, xs_t = [], [], [] for j, lc in enumerate(lc_grid): a = crossing(maxerr[:, j], pk_grid, T, True) b = crossing(maxerr[:, j], pk_grid, T, False) if np.isfinite(a) and np.isfinite(b): xs_t.append(lc); lo_t.append(a); hi_t.append(b) xs_t = np.array(xs_t); lo_t = np.array(lo_t); hi_t = np.array(hi_t) sl, bl = np.polyfit(xs_t, lo_t, 1) sh, bh = np.polyfit(xs_t, hi_t, 1) xstar = (bh - bl) / (sl - sh) j0 = int(np.argmin(np.abs(lc_grid + 1.0))) a0 = crossing(maxerr[:, j0], pk_grid, T, True) b0 = crossing(maxerr[:, j0], pk_grid, T, False) print(" %-8.2f %-13.3f %-13.3f %-13.3f %-13.3e %.3f" % (T, a0, b0, b0 - a0, 10.0 ** xstar, bl + sl * xstar)) sub("Real acids at 0.1 M, worst HH error over the buffer window") named = [ ("trichloroacetic", 0.66), ("dichloroacetic", 1.35), ("chloroacetic", 2.865), ("citric (pKa1)", 3.13), ("hydrofluoric", 3.17), ("formic", 3.75), ("benzoic", 4.20), ("acetic", 4.756), ("propanoic", 4.874), ("carbonic (pKa1)", 6.35), ("dihydrogenphosphate", 7.198), ("hypochlorous", 7.53), ("ammonium", 9.25), ("hydrocyanic", 9.21), ("phenol", 9.99), ("bicarbonate", 10.33), ("methylammonium", 10.66), ("hydrogenphosphate", 12.375), ] named.sort(key=lambda t: t[1]) print("%-22s %-8s %-14s %-12s %-12s %s" % ("acid", "pKa", "worst |err|", "at phi", "err phi=0.5", "inside 0.05?")) named_rows = [] for nm, pkv in named: CAn, CBn = titration_conc(0.1, phi_grid) pHn = solve_ph(CAn, CBn, [10.0 ** -pkv]) en = pkv + np.log10(phi_grid / (1 - phi_grid)) - pHn kk = int(np.argmax(np.abs(en))) k5n = int(np.argmin(np.abs(phi_grid - 0.5))) named_rows.append((nm, pkv, abs(en[kk]), phi_grid[kk], en[k5n])) print("%-22s %-8.3f %-14.4f %-12.2f %+-12.5f %s" % (nm, pkv, abs(en[kk]), phi_grid[kk], en[k5n], "yes" if abs(en[kk]) <= 0.05 else "NO")) print(" %d of the %d land inside 0.05 pH at this concentration." % (sum(1 for r in named_rows if r[2] <= 0.05), len(named_rows))) sub("The narrow reading: HH restricted to phi between 0.2 and 0.8") kl = int(np.argmin(np.abs(phi_grid - 0.2))) kh = int(np.argmin(np.abs(phi_grid - 0.8))) maxerr_narrow = aerr[:, :, kl:kh + 1].max(axis=2) print(" %-10s %-13s %-13s %-13s %-13s" % ("C (M)", "lo .05 wide", "lo .05 narrow", "hi .05 wide", "hi .05 narrow")) gains = [] for j in range(len(lc_grid)): a = crossing(maxerr[:, j], pk_grid, 0.05, True) b = crossing(maxerr_narrow[:, j], pk_grid, 0.05, True) gains.append(a - b) for lc in [0.0, -1.0, -2.0, -3.0, -4.0, -5.0]: j = int(np.argmin(np.abs(lc_grid - lc))) print(" %-10.0e %-13.3f %-13.3f %-13.3f %-13.3f" % (10.0 ** lc, crossing(maxerr[:, j], pk_grid, 0.05, True), crossing(maxerr_narrow[:, j], pk_grid, 0.05, True), crossing(maxerr[:, j], pk_grid, 0.05, False), crossing(maxerr_narrow[:, j], pk_grid, 0.05, False))) print(" Mean gain on the strong-acid side: %.3f pKa units." % float(np.nanmean(gains))) print(" Predicted gain, 2 log10(0.2/0.1) = %.3f. Real, and modest." % (2 * np.log10(2.0))) # ====================================================================== banner("PART 4. THE INITIAL-pH FORMULAS") # ====================================================================== print("At phi = 0 the taught answer is [H+] = sqrt(Ka C). Two things get") print("dropped: the loss of HA to its own dissociation, and the protons water") print("supplies. QUAD keeps the first and drops the second, which separates") print("the two sources.") print() PKm, LCm = np.meshgrid(pk_grid, lc_grid, indexing="ij") Cm = 10.0 ** LCm pH0_exact = solve_ph(Cm, np.zeros_like(Cm), [10.0 ** -PKm]) e_sqrt = ph_sqrt(PKm, Cm) - pH0_exact e_quad = ph_quad(PKm, Cm) - pH0_exact print("%-7s %-10s %-14s %-12s %-12s %-12s %-12s" % ("pKa", "C (M)", "exact pH", "sqrt pH", "sqrt err", "quad pH", "quad err")) for pkv, cv in [(4.75, 1.0), (4.75, 0.1), (4.75, 1e-3), (4.75, 1e-5), (2.00, 0.1), (1.00, 0.1), (0.00, 0.1), (0.00, 1e-3), (7.00, 1e-4), (9.00, 1e-4), (10.0, 1e-4), (12.0, 1e-3), (13.0, 1e-2)]: i = int(np.argmin(np.abs(pk_grid - pkv))) j = int(np.argmin(np.abs(lc_grid - np.log10(cv)))) print("%-7.2f %-10.0e %-14.6f %-12.6f %+-12.4f %-12.6f %+-12.4f" % (pk_grid[i], 10.0 ** lc_grid[j], pH0_exact[i, j], ph_sqrt(pk_grid[i], 10.0 ** lc_grid[j]), e_sqrt[i, j], ph_quad(pk_grid[i], 10.0 ** lc_grid[j]), e_quad[i, j])) sub("Fraction of the (pKa, C) grid where each initial-pH formula fails") for label, E in [("sqrt(Ka C)", e_sqrt), ("quadratic", e_quad)]: print(" %-12s >0.05 pH : %6.2f %% >0.10 pH : %6.2f %% worst : %.3f pH" % (label, 100 * np.mean(np.abs(E) > 0.05), 100 * np.mean(np.abs(E) > 0.10), np.max(np.abs(E)))) sub("The five percent rule, tested") print(" Textbooks say sqrt(Ka C) is safe when the dissociated fraction is") print(" under 5 percent, that is sqrt(Ka/C) < 0.05, that is Ka/C < 2.5e-3.") ratio = 10.0 ** (-PKm) / Cm inside = ratio < 2.5e-3 print(" cells satisfying the rule : %d of %d" % (inside.sum(), inside.size)) print(" of those, |sqrt error| > 0.05 pH : %d (%.2f %% of the obedient cells)" % (int((np.abs(e_sqrt)[inside] > 0.05).sum()), 100 * np.mean(np.abs(e_sqrt)[inside] > 0.05))) print(" worst sqrt error inside the rule : %.4f pH" % np.max(np.abs(e_sqrt)[inside])) print(" worst quad error inside the rule : %.4f pH" % np.max(np.abs(e_quad)[inside])) hit = int(np.argmax(np.where(inside, np.abs(e_sqrt), -1.0))) hi_, hj_ = np.unravel_index(hit, e_sqrt.shape) print(" worst offender obeying the rule : pKa %.2f, C %.1e M, error %+.3f pH" % (pk_grid[hi_], 10.0 ** lc_grid[hj_], e_sqrt[hi_, hj_])) print(" The rule controls the dissociation error it was written for and says") print(" nothing at all about the water error, which is what actually bites.") sub("At the equivalence point: pOH = (pKb - log10 C')/2") print(" C' is the conjugate-base concentration after dilution, C/2 here.") print() print("%-7s %-10s %-14s %-14s %-12s" % ("pKa", "C (M)", "exact pH", "formula pH", "difference")) eq_rows = [] for pkv in [2.0, 3.0, 4.0, 4.75, 6.0, 8.0, 10.0, 12.0]: for cv in [0.1, 1e-3]: CAe, CBe = titration_conc(cv, 1.0) pe = float(solve_ph(np.array([CAe]), np.array([CBe]), [10.0 ** -pkv])[0]) pf = float(ph_equiv(pkv, CAe)) eq_rows.append((pkv, cv, pe, pf)) print("%-7.2f %-10.0e %-14.6f %-14.6f %+-12.4f" % (pkv, cv, pe, pf, pf - pe)) print(" The formula treats A- as a weak base in pure water and ignores both") print(" the leftover HA and the water contribution. It fails from below for") print(" strong acids and from above for very weak ones.") # ====================================================================== banner("PART 5. DROPPING WATER AUTOIONISATION") # ====================================================================== print("Kw set to zero, everything else exact. Evaluated over phi 0.10 to 0.90,") print("because with no water term the equation has no solution at all at") print("phi = 1: there is nothing left to balance the charge.") print() pH_now = solve_ph_nowater(CA3, CB3, [10.0 ** -PK]) e_now = np.abs(pH_now - pH_exact).max(axis=2) print("fraction of %d cells where dropping water costs more than:" % e_now.size) for t in [0.01, 0.05, 0.10, 0.50, 1.00]: print(" %.2f pH : %6.2f %%" % (t, 100 * np.mean(e_now > t))) print() print("%-10s %-11s %-11s %-11s %-11s %-11s" % ("C (M)", "pKa 2", "pKa 5", "pKa 8", "pKa 10", "pKa 12")) for lc in [0.0, -1.0, -2.0, -3.0, -4.0, -5.0]: j = int(np.argmin(np.abs(lc_grid - lc))) row = [] for pkv in [2.0, 5.0, 8.0, 10.0, 12.0]: i = int(np.argmin(np.abs(pk_grid - pkv))) row.append(e_now[i, j]) print("%-10.0e %-11.5f %-11.5f %-11.5f %-11.5f %-11.5f" % (10.0 ** lc, row[0], row[1], row[2], row[3], row[4])) sub("The worst cells for dropping water") order = np.argsort(e_now.ravel())[::-1][:6] for o in order: i, j = np.unravel_index(int(o), e_now.shape) print(" pKa %-6.2f C %-10.2e cost of dropping water = %.3f pH" % (pk_grid[i], 10.0 ** lc_grid[j], e_now[i, j])) print(" Water is irrelevant for an ordinary carboxylic acid at 0.1 M and") print(" decisive for a phenol at 1e-5 M. Same equation, same solvent.") # ====================================================================== banner("PART 6. POLYPROTIC ACIDS AND THE OVERLAP PROBLEM") # ====================================================================== poly_acids = [ ("oxalic", [1.25, 4.27], "H2C2O4"), ("maleic", [1.92, 6.27], "C4H4O4"), ("malonic", [2.83, 5.69], "C3H4O4"), ("phthalic", [2.95, 5.41], "C8H6O4"), ("succinic", [4.21, 5.64], "C4H6O4"), ("carbonic", [6.35, 10.33], "H2CO3"), ("citric", [3.13, 4.76, 6.40], "C6H8O7"), ("phosphoric", [2.148, 7.198, 12.375], "H3PO4"), ] print("HH applied to the FIRST buffer region of a polyprotic acid at C = 0.1 M,") print("phi counted against the first proton only, window 0.1 to 0.9. dpK is the") print("gap to the next pKa: the smaller it is, the worse HH does, because the") print("second equilibrium is already running underneath the first.") print() print("%-12s %-9s %-9s %-7s %-13s %-10s %-12s" % ("acid", "pKa1", "formula", "dpK", "worst |err|", "at phi", "err phi=0.5")) poly_rows = [] phi_p = np.linspace(0.1, 0.9, 81) for name, pks, formula in poly_acids: Kas = [10.0 ** -p for p in pks] CAp, CBp = titration_conc(0.1, phi_p) pHp = solve_ph(CAp, CBp, Kas) e = pks[0] + np.log10(phi_p / (1.0 - phi_p)) - pHp k = int(np.argmax(np.abs(e))) k5p = int(np.argmin(np.abs(phi_p - 0.5))) poly_rows.append((name, pks, abs(e[k]), phi_p[k], e[k5p])) print("%-12s %-9.3f %-9s %-7.2f %-13.4f %-10.3f %+-12.4f" % (name, pks[0], formula, pks[1] - pks[0], abs(e[k]), phi_p[k], e[k5p])) sub("Monoprotic control at the same pKa1, for comparison") print("%-12s %-9s %-13s %-13s %s" % ("acid", "pKa1", "poly worst", "mono worst", "extra cost of proton 2")) for (name, pks, formula), row in zip(poly_acids, poly_rows): CAm, CBm = titration_conc(0.1, phi_p) pHm = solve_ph(CAm, CBm, [10.0 ** -pks[0]]) em = float(np.max(np.abs(pks[0] + np.log10(phi_p / (1 - phi_p)) - pHm))) print("%-12s %-9.3f %-13.4f %-13.4f %+.4f" % (name, pks[0], row[2], em, row[2] - em)) sub("The amphiprotic shortcut pH = (pK1 + pK2)/2 for the salt NaHA") print(" A solution of the intermediate salt prepared directly at") print(" concentration C, so C_A = C_B = C, no dilution involved.") print() print("%-12s %-10s %-14s %-14s %-13s %s" % ("salt of", "C (M)", "exact pH", "(pK1+pK2)/2", "difference", "refined err")) amph_rows = [] for name, pks, formula in poly_acids: Kas = [10.0 ** -p for p in pks] for C in [0.1, 1e-3]: pHa = float(solve_ph(np.array([C]), np.array([C]), Kas)[0]) simple = 0.5 * (pks[0] + pks[1]) K1, K2 = Kas[0], Kas[1] hbet = np.sqrt((K1 * K2 * C + K1 * KW) / (K1 + C)) refined = -np.log10(hbet) amph_rows.append((name, C, pHa, simple, refined)) print("%-12s %-10.0e %-14.6f %-14.4f %+-13.4f %+.4f" % (name, C, pHa, simple, simple - pHa, refined - pHa)) print(" The last column is the standard refinement") print(" h = sqrt((K1 K2 C + K1 Kw)/(K1 + C)), printed as its own error.") print(" It is much better and it is not error free either.") # ====================================================================== banner("PART 7. MONTE CARLO OVER THE PARAMETER SPACE") # ====================================================================== NMC = 200000 print("Two regions, %d draws each, from the single seeded stream. A draw is a" % NMC) print("(pKa, C, phi) triple drawn uniformly in pKa, in log10 C and in phi; the") print("statistic is whether HH's error at that point exceeds a threshold.") print("Standard errors are binomial, taken from the runs themselves.") regions = [ ("classroom", (3.0, 10.0), (-2.0, 0.0), (0.2, 0.8)), ("wide", (0.0, 14.0), (-5.0, 0.0), (0.05, 0.95)), ] mc_store = {} for label, pkr, lcr, phr in regions: pk_s = rng.uniform(pkr[0], pkr[1], NMC) lc_s = rng.uniform(lcr[0], lcr[1], NMC) ph_s = rng.uniform(phr[0], phr[1], NMC) CAs, CBs = titration_conc(10.0 ** lc_s, ph_s) pe = solve_ph(CAs, CBs, [10.0 ** -pk_s]) eh = np.abs(pk_s + np.log10(ph_s / (1 - ph_s)) - pe) mc_store[label] = eh sub("region %s: pKa %s, log10 C %s, phi %s" % (label, pkr, lcr, phr)) for t in [0.02, 0.05, 0.10, 0.20]: p = float(np.mean(eh > t)) se = float(np.sqrt(p * (1 - p) / NMC)) print(" P(|HH error| > %.2f pH) = %.5f +/- %.5f (%.3f %% +/- %.3f %%)" % (t, p, se, 100 * p, 100 * se)) print(" mean |HH error| = %.5f +/- %.5f pH" % (float(np.mean(eh)), float(np.std(eh, ddof=1) / np.sqrt(NMC)))) print(" median |HH error| = %.5f pH" % float(np.median(eh))) print(" 90th percentile = %.5f pH" % float(np.percentile(eh, 90))) print(" 99th percentile = %.5f pH" % float(np.percentile(eh, 99))) print(" largest drawn = %.5f pH" % float(eh.max())) sub("Convergence of the classroom-region estimate of P(error > 0.05)") eh = mc_store["classroom"] run = np.cumsum(eh > 0.05) ks = np.unique(np.round(np.logspace(2, np.log10(NMC), 34)).astype(int)) print(" %-10s %-13s %-13s %s" % ("k trials", "running p", "SE", "p +/- 2 SE")) conv = [] for k in ks: p = float(run[k - 1]) / k se = float(np.sqrt(max(p * (1 - p), 1e-12) / k)) conv.append((int(k), p, se)) print(" %-10d %-13.5f %-13.5f %.5f to %.5f" % (k, p, se, p - 2 * se, p + 2 * se)) final_p, final_se = conv[-1][1], conv[-1][2] tail = [c for c in conv if c[0] >= 1000] ok = sum(1 for k, p, se in tail if abs(p - final_p) <= 3 * se) print(" of the %d checkpoints at k >= 1000, %d sit within 3 of their own SE" % (len(tail), ok)) print(" of the final estimate %.5f." % final_p) sub("The same statistic on a deterministic grid, as a cross-check") print(" The Monte Carlo draws uniformly in pKa, log10 C and phi, so the") print(" statistic is a volume fraction and a deterministic quadrature of") print(" the same box has to reproduce it. Two quadratures, printed because") print(" the first one is wrong and the reason is worth having in print.") print() bi = (pk_grid >= 3.0) & (pk_grid <= 10.0) bj = (lc_grid >= -2.0) & (lc_grid <= 0.0) bk = (phi_grid >= 0.2) & (phi_grid <= 0.8) sub_a = aerr[np.ix_(bi, bj, bk)] p_lattice = float(np.mean(sub_a > 0.05)) z_lat = (p_lattice - final_p) / final_se print(" (a) reuse the Part 2 lattice, every node weighted equally.") print(" estimate %.5f over %d nodes, %+.5f from the Monte Carlo," % (p_lattice, sub_a.size, p_lattice - final_p)) print(" which is %.1f Monte Carlo standard errors out. That is a real" % z_lat) print(" disagreement and it is the lattice's fault, not the sampler's.") print(" The set where HH fails inside this box hugs the pKa = 3 and") print(" pKa = 10 faces and the phi = 0.2 and phi = 0.8 faces, and a") print(" closed lattice puts full-weight nodes on all six of them. It") print(" is the trapezoid end-point error, and it inflates the answer.") print() nA, nB, nC = 70, 60, 64 ea = 7.0 / nA eb = 2.0 / nB ec = 0.6 / nC pk_m = 3.0 + (np.arange(nA) + 0.5) * ea lc_m = -2.0 + (np.arange(nB) + 0.5) * eb ph_m = 0.2 + (np.arange(nC) + 0.5) * ec PKq, LCq, PHq = np.meshgrid(pk_m, lc_m, ph_m, indexing="ij") CAq, CBq = titration_conc(10.0 ** LCq, PHq) eq = np.abs(PKq + np.log10(PHq / (1 - PHq)) - solve_ph(CAq, CBq, [10.0 ** -PKq])) p_mid = float(np.mean(eq > 0.05)) z_mid = (p_mid - final_p) / final_se print(" (b) the same box on a midpoint rule, %d x %d x %d cell centres," % (nA, nB, nC)) print(" no node on any face.") print(" estimate %.5f over %d nodes" % (p_mid, eq.size)) print(" Monte Carlo %.5f +/- %.5f" % (final_p, final_se)) print(" difference %+.5f = %.2f Monte Carlo standard errors" % (p_mid - final_p, z_mid)) print() print(" The midpoint quadrature and the sampler agree. The two methods") print(" share no code beyond the solver, so this is the check that the") print(" Monte Carlo is measuring what we think it measures.") # ====================================================================== banner("PART 8. SENSITIVITY TO THE DILUTION CHOICE") # ====================================================================== print("Everything above titrates with base at the same concentration as the") print("analyte, so the solution is diluted by 1 + phi. A more concentrated") print("titrant dilutes less and flatters HH. The whole map rerun at") print("titrant/analyte ratios of 1, 10 and 100:") print() print("%-8s %-16s %-16s %-16s %-16s" % ("ratio", "median max err", "cells >0.05 (%)", "pKa lo at 0.1 M", "pKa hi at 0.1 M")) for r in [1.0, 10.0, 100.0]: CAr, CBr = titration_conc(Cs, PHI, ratio=r) per = np.abs(PK + np.log10(PHI / (1 - PHI)) - solve_ph(CAr, CBr, [10.0 ** -PK])) mr = per.max(axis=2) print("%-8.0f %-16.5f %-16.2f %-16.3f %-16.3f" % (r, float(np.median(mr)), 100 * float(np.mean(mr > 0.05)), crossing(mr[:, j01], pk_grid, 0.05, True), crossing(mr[:, j01], pk_grid, 0.05, False))) print() print("Dilution is a second-order effect on the boundary, worth log10(1+phi)") print("in pKa and no more. The boundary is set by the ratio of [H+] to the") print("acid concentration, and diluting moves both of them together.") # ====================================================================== banner("PART 8b. TWO MODEL INPUTS THAT WOULD MOVE EVERYTHING") # ====================================================================== sub("Kw is an input, not a result") print(" Every number above uses Kw = 1.000e-14 exactly. Kw is measured, it") print(" depends on temperature, and the base-side boundary is set by it.") print(" The whole boundary recomputed at three values of pKw:") print() print(" %-8s %-16s %-16s %-16s" % ("pKw", "pKa lo at 0.1 M", "pKa hi at 0.1 M", "band width")) for pkw in [13.90, 14.00, 14.10]: kw = 10.0 ** -pkw pex = solve_ph(CA3, CB3, [10.0 ** -PK], Kw=kw) mk = np.abs(PK + np.log10(PHI / (1 - PHI)) - pex).max(axis=2) a = crossing(mk[:, j01], pk_grid, 0.05, True) b = crossing(mk[:, j01], pk_grid, 0.05, False) print(" %-8.2f %-16.3f %-16.3f %-16.3f" % (pkw, a, b, b - a)) print(" The acid-side edge does not move at all and the base-side edge") print(" moves one for one with pKw, which is what the closed form says it") print(" should do. A titration at body temperature has a different Kw and") print(" therefore a different band.") sub("Activity coefficients, the assumption that costs the most") print(" Everything above is an ideal solution. A real titration has an") print(" ionic strength that climbs as base goes in, and a glass electrode") print(" reads activity, not concentration. Redone with the Davies equation") print(" on a monoprotic acid: conditional constants Ka/gamma^2 and") print(" Kw/gamma^2, ionic strength I = C_B + [H+] closed self-consistently,") print(" pH reported as -log10(gamma [H+]).") print() print(" %-7s %-9s %-8s %-11s %-11s %-11s %-11s" % ("pKa", "C (M)", "I at 0.5", "gamma", "ideal pH", "Davies pH", "shift")) for pkv, cv in [(4.756, 1.0), (4.756, 0.1), (4.756, 0.01), (4.756, 1e-3), (2.865, 0.1), (7.198, 0.1), (9.99, 0.1), (9.25, 0.1)]: CAd, CBd = titration_conc(cv, 0.5) pid = float(solve_ph(np.array([CAd]), np.array([CBd]), [10.0 ** -pkv])[0]) pda, gg, Id = solve_ph_davies(np.array([CAd]), np.array([CBd]), pkv) print(" %-7.3f %-9.0e %-8.5f %-11.5f %-11.5f %-11.5f %+-11.5f" % (pkv, cv, float(Id[0]), float(gg[0]), pid, float(pda[0]), float(pda[0]) - pid)) print() pk_d = np.arange(0.0, 14.001, 0.5) lc_d = np.arange(-5.0, 0.001, 0.25) PKd, LCd, PHId = np.meshgrid(pk_d, lc_d, phi_grid, indexing="ij") CAd2, CBd2 = titration_conc(10.0 ** LCd, PHId) pHd, gd, Id2 = solve_ph_davies(CAd2, CBd2, PKd) pHi2 = solve_ph(CAd2, CBd2, [10.0 ** -PKd]) e_hh_ideal = np.abs(PKd + np.log10(PHId / (1 - PHId)) - pHi2).max(axis=2) e_hh_dav = np.abs(PKd + np.log10(PHId / (1 - PHId)) - pHd).max(axis=2) e_model = np.abs(pHd - pHi2).max(axis=2) jd = int(np.argmin(np.abs(lc_d + 1.0))) print(" at C = 0.1 M the 0.05 pH band, ideal : pKa %.3f to %.3f" % (crossing(e_hh_ideal[:, jd], pk_d, 0.05, True), crossing(e_hh_ideal[:, jd], pk_d, 0.05, False))) bl = crossing(e_hh_dav[:, jd], pk_d, 0.05, True) bh = crossing(e_hh_dav[:, jd], pk_d, 0.05, False) print(" at C = 0.1 M the 0.05 pH band, Davies : pKa %s to %s" % (("%.3f" % bl) if np.isfinite(bl) else "empty", ("%.3f" % bh) if np.isfinite(bh) else "empty")) print(" cells where HH beats 0.05 pH, ideal : %.2f %%" % (100 * np.mean(e_hh_ideal <= 0.05))) print(" cells where HH beats 0.05 pH, Davies : %.2f %%" % (100 * np.mean(e_hh_dav <= 0.05))) print(" largest ideal-to-Davies gap anywhere : %.4f pH" % e_model.max()) print(" median ideal-to-Davies gap : %.4f pH" % float(np.median(e_model))) print() print(" Switching the solvent model on is worth more than every") print(" approximation error we have measured, at every concentration a") print(" student actually uses. The band we drew is the band for an ideal") print(" solution, and that is the honest name for it.") # ====================================================================== banner("PART 9. FIGURE DATA") # ====================================================================== sub("FIG1: exact titration curves, C = 0.1 M, titrant 0.1 M") phi_f = np.concatenate([np.linspace(0.002, 0.05, 13), np.linspace(0.06, 0.96, 91), np.array([0.97, 0.98, 0.99, 0.995, 0.999, 1.0, 1.001, 1.005, 1.01, 1.02, 1.03]), np.linspace(1.04, 1.30, 27)]) print(" phi: %s" % " ".join("%.4f" % x for x in phi_f)) for pkv in [2.0, 4.75, 10.0]: CAf, CBf = titration_conc(0.1, phi_f) pe = solve_ph(CAf, CBf, [10.0 ** -pkv]) print(" pKa=%.2f exact: %s" % (pkv, " ".join("%.4f" % x for x in pe))) sub("FIG2: HH error against phi at C = 0.1 M") phi_e = np.linspace(0.05, 0.95, 46) print(" phi: %s" % " ".join("%.3f" % x for x in phi_e)) for pkv in [1.0, 2.0, 3.0, 4.75, 10.0, 11.0, 12.0]: CAe2, CBe2 = titration_conc(0.1, phi_e) pe = solve_ph(CAe2, CBe2, [10.0 ** -pkv]) er = pkv + np.log10(phi_e / (1 - phi_e)) - pe print(" pKa=%-6.2f: %s" % (pkv, " ".join("%+.4f" % x for x in er))) sub("FIG3: boundary curves in (log10 C, pKa)") print(" %-10s %-10s %-10s %-10s %-10s" % ("log10C", "lo05", "hi05", "lo10", "hi10")) for j, lc in enumerate(lc_grid): r = band_rows[j] print(" %-10.3f %-10.4f %-10.4f %-10.4f %-10.4f" % (lc, r[1], r[2], r[3], r[4])) sub("FIG4: initial-pH error surface, signed sqrt-formula error in pH") print(" rows pKa 0 to 14 step 0.5, columns log10 C -5 to 0 step 0.5") for pkv in np.arange(0.0, 14.001, 0.5): i = int(np.argmin(np.abs(pk_grid - pkv))) vals = [] for lc in np.arange(-5.0, 0.001, 0.5): j = int(np.argmin(np.abs(lc_grid - lc))) vals.append(e_sqrt[i, j]) print(" pKa %-6.2f %s" % (pk_grid[i], " ".join("%+.4f" % v for v in vals))) sub("FIG5: Monte Carlo convergence trace, classroom region") print(" %-10s %-13s %-13s" % ("k", "p", "se")) for k, p, se in conv: print(" %-10d %-13.6f %-13.6f" % (k, p, se)) # ====================================================================== banner("PART 10. SUMMARY") # ====================================================================== print("exact charge-balance solutions computed : %d" % SOLVE_COUNT[0]) print("charge balance worst relative residual : %.3e (%.1f eps)" % (resid_rel_worst, resid_rel_worst / EPS)) print("mass balance worst |sum alpha - 1| : %.3e (%.1f eps)" % (worst_mass, worst_mass / EPS)) print("strong acid / strong base worst error : %.3e pH over 28 points" % v3worst) print("bisection vs numpy.roots worst difference: %.3e pH over 18 points" % v5worst) print("bisection vs 60-digit decimal worst diff : %.3e pH over 12 points" % v6worst) print("pure water pH : %.12f against 7.000000000000" % ph_w) print() print("HH in the classroom box (pKa 3-10, C 0.01-1 M, phi 0.2-0.8):") print(" P(error > 0.05 pH) = %.5f +/- %.5f" % (final_p, final_se)) print(" mean |error| = %.5f pH" % float(np.mean(mc_store["classroom"]))) pw = float(np.mean(mc_store["wide"] > 0.05)) print("HH over the wide box (pKa 0-14, C 1e-5 to 1 M, phi 0.05-0.95):") print(" P(error > 0.05 pH) = %.5f +/- %.5f" % (pw, np.sqrt(pw * (1 - pw) / NMC))) print() print("at C = 1 M the 0.05 pH band runs pKa %.2f to %.2f" % (band_rows[-1][1], band_rows[-1][2])) print("at C = 0.1 M the 0.05 pH band runs pKa %.2f to %.2f" % (band_rows[j01][1], band_rows[j01][2])) j03 = int(np.argmin(np.abs(lc_grid + 3.0))) print("at C = 1e-3 M it runs pKa %.2f to %.2f" % (band_rows[j03][1], band_rows[j03][2])) print("at C = 1e-4 M and below the 0.05 pH band is empty: no pKa works") print() print("acetic acid, pKa 4.756, at 0.1 M: worst buffer-window HH error = %.4f pH" % float(np.max(np.abs( 4.756 + np.log10(phi_grid / (1 - phi_grid)) - solve_ph(*titration_conc(0.1, phi_grid), [10.0 ** -4.756]))))) print("chloroacetic acid, pKa 2.865, at 0.1 M: worst error = %.4f pH" % float(np.max(np.abs( 2.865 + np.log10(phi_grid / (1 - phi_grid)) - solve_ph(*titration_conc(0.1, phi_grid), [10.0 ** -2.865]))))) print("phenol, pKa 9.99, at 0.1 M: worst error = %.4f pH" % float(np.max(np.abs( 9.99 + np.log10(phi_grid / (1 - phi_grid)) - solve_ph(*titration_conc(0.1, phi_grid), [10.0 ** -9.99]))))) print() print("runtime: %.1f s" % (time.time() - T0)) if __name__ == "__main__": main()