#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ The common ion effect is taught backwards. Science Journaling Club, Volume 1 Issue 3, Spring 2025. THE QUESTION ------------ General chemistry states the common ion effect as an unconditional rule: add an ion the precipitate already contains, and its solubility falls, full stop. The rule is Le Chatelier's principle read off Ksp = [M][A] with nothing else in the beaker. Real beakers hold more than a solubility product. The same anion that suppresses solubility through Ksp can also bind the very cation it precipitated, as a soluble complex, and if that second reaction is strong enough the net effect on solubility reverses: more common ion, MORE dissolved salt. This script asks, numerically and without shortcuts, whether that reversal is real, how big it is, and exactly where the textbook rule stops applying. THE MODEL --------- For each of three sparingly soluble salts we solve, at every point on a grid of added common-ion concentration and pH, the full coupled system: * the solubility product Ksp = [M][A] (or [M][A]^2), * stepwise complex formation of the cation with the same anion, M + nA <-> MA_n, cumulative constant beta_n = [MA_n]/([M][A]^n), * acid-base speciation of the anion, HA <-> H+ + A-, Ka, * the water equilibrium Kw = [H+][OH-], which fixes [OH-] once pH is treated as the externally controlled variable it is in every solubility-vs-pH diagram in the literature (the pH of a beaker held at a stated value by an external acid or buffer, not solved from a closed charge balance on the salt alone). No approximation is used inside the solver: every mass balance is written in full (total dissolved cation = free ion + every complex; total anion added = free conjugate-base ion, plus its protonated form, plus every anion tied up in a complex, weighted by stoichiometry) and solved for the free ligand concentration by bisection on a strictly monotonic residual (proved monotonic below, and checked numerically at every grid point). Once the free ligand concentration is known every other concentration follows algebraically from the equilibrium expressions, and every one of those expressions, plus every mass balance and an explicit charge balance built from an assumed spectator counter-ion, is re-substituted afterward and its residual printed. That is the validation this script exists to do; see VALIDATION below and PART 1 of the printed output. THE THREE SALTS ---------------- Ba(IO3)2 common ion IO3- (added as NaIO3). No complex formation is included for Ba2+-IO3-: none is well documented, and the club found none. The textbook rule should hold everywhere. Used to validate the solver against closed-form algebra and as the "control" the other two are measured against. CaF2 common ion F- (added as NaF). Fluoride is a real weak base (Ka(HF) = 6.6e-4, pKa 3.18) and CaF2 solubility genuinely depends on pH through F- protonation, while CaF+ ion pairing is neglected (named explicitly as a limitation below). Shows that pH can move solubility by orders of magnitude through a channel the common-ion rule never mentions, with NO reversal with respect to F- itself: this salt keeps the textbook honest about what "the rule holds" actually requires. AgCl common ion Cl- (added as NaCl). Silver forms a real, measured ladder of chloro-complexes, AgCl(aq), AgCl2-, AgCl3^2-, AgCl4^3-, with stepwise formation constants traceable to the original potentiometric measurement of Jonte & Martin (1952) [4] and reproduced in analytical chemistry compilations since. This is the salt that breaks the rule: solubility falls with added chloride, reaches a minimum, and then RISES, by orders of magnitude in concentrated chloride, because the same ion that used to suppress solubility now dissolves the precipitate as a soluble complex ion. EQUILIBRIUM CONSTANTS USED, AND WHERE THEY COME FROM ----------------------------------------------------- Kw(25C) = 1.000e-14 standard value, consistent with Bandura & Lvov (2006) [11] at 298.15 K (pKw = 13.995 there). Ksp(AgCl, 25C) = 1.77e-10 CRC Handbook of Chemistry and Physics, 99th ed. (2018), as tabulated on the compound's standard reference entry. Self-consistency check: sqrt(Ksp) gives the textbook pure-water AgCl solubility of 1.33e-5 mol/L, reproduced below to 4 figures. log beta_n(AgCl_n) = 3.04, 5.04, 5.04, 5.30 (n=1..4), cumulative, 25C. Traceable to Jonte & Martin (1952) [4], the original potentiometric solubility study of AgCl in chloride media, as reproduced in standard analytical-chemistry stability-constant compilations. IMPORTANT HONESTY NOTE: this script could not access the original paper's text (it is paywalled) to quote its own stated precision. We therefore ASSIGN an assumed measurement uncertainty of +/-0.05 in each log beta for the Monte Carlo propagation in PART 1, Section V5, rather than pretending to a precision we cannot verify. This is flagged again in SS10 of the article as the strongest objection to the headline numbers. Ksp(CaF2, 25C) = 3.89e-11 (pKsp = 10.41) Macaskill & Bates (1977) [5], a potentiometric redetermination. Ka(HF) = 6.6e-4 (pKa 3.18) standard ambient value, consistent with Usha & Atkinson (1992) [6], who report the pressure dependence of this constant relative to its 1-atm reference value. Ksp(Ba(IO3)2, 25C) = 1.57e-9 a standard tabulated value; self-consistent with the commonly quoted pure-water solubility of 7.32e-4 mol/L via S = (Ksp/4)^(1/3), reproduced below. Ka(HIO3) = 0.17 (pKa 0.77) standard tabulated value for iodic acid, a moderately strong acid; included for completeness though it only matters below about pH 2. log beta_n(Ag(NH3)_n) = 3.24, 7.05 (n=1,2). Standard tabulated general- chemistry formation constants for the diammine- silver(I) complex, used ONLY for the supplementary two-ligand regime map (Figure 4) and the interactive model's third slider. These are presented as representative textbook values, not independently re-verified against a primary source this session, and no scientific claim in the article depends on their fourth significant figure. ASSUMPTIONS ----------- * 25 degrees C throughout. * Ideal solution: activity coefficients are 1, concentrations stand in for activities everywhere. This is the same idealisation the textbook rule itself is normally taught in, so the comparison between "exact" and "naive" below is fair, even though neither one is a real ionic-strength-corrected answer. See Hunter (1966) [14] on what the Debye-Huckel correction would actually cost here. * Cation-anion complexation is included only where the club found a documented ladder of constants (Ag+/Cl-); Ca2+-F- ion pairing and Ba2+-IO3- association are both assumed negligible and are not in the solver for those two salts. Named, not hidden. * pH is an externally controlled variable (an implicit buffer or titrated strong acid/base), not solved from a self-contained charge balance on the salt and water alone. Section V of the validation checks that a physically sensible spectator-ion charge balance can be built to support every solved state. * The common-ion source is a fully dissociated 1:1 salt (NaCl, NaF, NaIO3) with a non-interacting spectator cation. LIMITATIONS ----------- * This is a computation, not an experiment. No salt was weighed, no precipitate was filtered, no pH meter was read. Every number below is the output of the equations above; its contact with a real beaker is limited by every assumption just listed, above all the activity-coefficient idealisation. * The AgCl formation constants carry a real, unquantified-by-us uncertainty (see above); Section V5 propagates an ASSUMED +/-0.05 log-unit uncertainty and reports what survives it and what does not. * Kinetics, supersaturation, and polymorph selection are outside the model entirely; Ksp describes the true equilibrium state only. VALIDATION ---------- Five checks, each printing the club's computed value beside the analytic, literature, or internally-consistent reference value, with the difference: (V1) mass-balance, complex-formation, and Ksp residuals at every grid point, to machine precision; (V2) the Ba(IO3)2 and CaF2 exact solver against a closed-form cubic solved independently by np.roots; (V3) AgCl pure-water solubility against the textbook closed form and the literature value; (V4) the AgCl solubility minimum location, both as an internal check (bisection-found minimum against the exact quartic root) and against the qualitative literature description of where it falls; (V5) Monte Carlo propagation of the formation-constant uncertainty, with standard errors and a convergence trace. Random seed: 20250306, numpy default_rng (PCG64). Fully deterministic. """ import sys import time import numpy as np SEED = 20250306 RNG = np.random.default_rng(SEED) KW = 1.0e-14 LN10 = np.log(10.0) EPS = np.finfo(float).eps T0 = time.time() # ---------------------------------------------------------------------------- # equilibrium constants (see docstring for sourcing) # ---------------------------------------------------------------------------- KSP_AGCL = 1.77e-10 LOGBETA_AGCL = np.array([3.04, 5.04, 5.04, 5.30]) # AgCl, AgCl2-, AgCl3^2-, AgCl4^3- BETA_AGCL = 10.0 ** LOGBETA_AGCL KA_HCL = 1.0e6 # HCl essentially fully dissociated KSP_CAF2 = 3.89e-11 KA_HF = 6.6e-4 KSP_BAIO3 = 1.57e-9 KA_HIO3 = 0.17 LOGBETA_AGNH3 = np.array([3.24, 7.05]) # AgNH3+, Ag(NH3)2+ (illustrative, SS "regime map") BETA_AGNH3 = 10.0 ** LOGBETA_AGNH3 def alpha_free(h, Ka): """Fraction of a monoprotic conjugate pair present as the free base A-, at proton concentration h, given the acid's Ka. f = Ka/(Ka+h).""" return Ka / (Ka + h) # ---------------------------------------------------------------------------- # generic bisection on a scalar monotone-increasing residual # ---------------------------------------------------------------------------- def bisect(residual, lo, hi, iters=100): flo = residual(lo) fhi = residual(hi) if flo > 0 or fhi < 0: raise RuntimeError("bracket does not straddle a root: f(lo)=%.6g f(hi)=%.6g" % (flo, fhi)) for _ in range(iters): mid = 0.5 * (lo + hi) if residual(mid) <= 0: lo = mid else: hi = mid return 0.5 * (lo + hi) # ---------------------------------------------------------------------------- # Salt 1: Ba(IO3)2(s) <-> Ba2+ + 2 IO3-, common ion IO3- added as NaIO3 # ---------------------------------------------------------------------------- def solve_baio3(C, pH, Ksp=KSP_BAIO3, Ka=KA_HIO3): """Molar solubility x of Ba(IO3)2 with C mol/L NaIO3 added, at fixed pH. Mass balance: total iodate = C + 2x = [IO3-] + [HIO3] = (C+2x)/f, so [IO3-] = f*(C+2x) with f = Ka/(Ka+h). Ksp = x*(f*(C+2x))^2, solved for x by bisection (monotone increasing in x for x>0).""" h = 10.0 ** (-pH) f = alpha_free(h, Ka) def resid(x): return x * (f * (C + 2 * x)) ** 2 - Ksp hi = max(1.0, (Ksp) ** (1.0 / 3.0) * 50 + C) x = bisect(resid, 0.0, hi, iters=120) io3 = f * (C + 2 * x) hio3 = (1 - f) * (C + 2 * x) ba = x return dict(x=x, Ba=ba, IO3=io3, HIO3=hio3, f=f, h=h) def baio3_cubic_root(C, pH, Ksp=KSP_BAIO3, Ka=KA_HIO3): """Independent closed-form check: Ksp = x*(f*(C+2x))^2 expands to 4 f^2 x^3 + 4 f^2 C x^2 + f^2 C^2 x - Ksp = 0, an exact cubic solved by np.roots, unrelated in method to the bisection above.""" h = 10.0 ** (-pH) f = alpha_free(h, Ka) coeffs = [4 * f * f, 4 * f * f * C, f * f * C * C, -Ksp] roots = np.roots(coeffs) real_pos = [r.real for r in roots if abs(r.imag) < 1e-12 * max(1.0, abs(r.real)) and r.real > 0] return min(real_pos) # ---------------------------------------------------------------------------- # Salt 2: CaF2(s) <-> Ca2+ + 2 F-, common ion F- added as NaF # ---------------------------------------------------------------------------- def solve_caf2(C, pH, Ksp=KSP_CAF2, Ka=KA_HF): """Identical structure to Ba(IO3)2: no Ca-F complexation modelled.""" h = 10.0 ** (-pH) f = alpha_free(h, Ka) def resid(x): return x * (f * (C + 2 * x)) ** 2 - Ksp hi = max(1.0, (Ksp) ** (1.0 / 3.0) * 200 + C + 1.0) x = bisect(resid, 0.0, hi, iters=120) fion = f * (C + 2 * x) hf = (1 - f) * (C + 2 * x) ca = x return dict(x=x, Ca=ca, F=fion, HF=hf, f=f, h=h) def caf2_cubic_root(C, pH, Ksp=KSP_CAF2, Ka=KA_HF): h = 10.0 ** (-pH) f = alpha_free(h, Ka) coeffs = [4 * f * f, 4 * f * f * C, f * f * C * C, -Ksp] roots = np.roots(coeffs) real_pos = [r.real for r in roots if abs(r.imag) < 1e-12 * max(1.0, abs(r.real)) and r.real > 0] return min(real_pos) # ---------------------------------------------------------------------------- # Salt 3: AgCl(s) <-> Ag+ + Cl-, common ion AND ligand Cl-, added as NaCl # ---------------------------------------------------------------------------- # # Let c = [Cl-] (free, deprotonated form; HCl protonation folded in via f). # Ksp = a*c => a = [Ag+] = Ksp/c. # Total dissolved Ag: x = a*(1 + sum_n beta_n c^n) # Total Cl mass balance: C + x = c/f + a*sum_n(n*beta_n*c^n) # Substituting a = Ksp/c and solving for C as a function of c alone gives # the single monotone equation used for bisection (derived in the article, # SS3; every term below is positive for c>0, beta_n>0, Ksp>0, so C(c) is # strictly increasing in c and the root is unique): # # C(c) = c/f - Ksp/c + Ksp*beta2*c + 2*Ksp*beta3*c^2 + 3*Ksp*beta4*c^3 def agcl_forward(c, f, Ksp, beta): b1, b2, b3, b4 = beta a = Ksp / c x = a * (1 + b1 * c + b2 * c * c + b3 * c ** 3 + b4 * c ** 4) C = c / f - Ksp / c + Ksp * b2 * c + 2 * Ksp * b3 * c * c + 3 * Ksp * b4 * c ** 3 return C, x, a def solve_agcl(C, pH, Ksp=KSP_AGCL, beta=BETA_AGCL, Ka=KA_HCL): h = 10.0 ** (-pH) f = alpha_free(h, Ka) def resid(logc): c = 10.0 ** logc Cc, x, a = agcl_forward(c, f, Ksp, beta) return Cc - C logc = bisect(resid, -14.0, 3.0, iters=140) c = 10.0 ** logc Cc, x, a = agcl_forward(c, f, Ksp, beta) b1, b2, b3, b4 = beta species = dict(AgCl0=a * b1 * c, AgCl2=a * b2 * c * c, AgCl3=a * b3 * c ** 3, AgCl4=a * b4 * c ** 4) return dict(x=x, Ag=a, Cl=c, f=f, h=h, **species) def agcl_quartic_minimum(beta=BETA_AGCL, Ksp=KSP_AGCL): """x(c) = Ksp/c + Ksp*b1 + Ksp*b2*c + Ksp*b3*c^2 + Ksp*b4*c^3 is convex in c>0 (sum of 1/c and increasing powers of c, all convex), so it has a single minimum where dx/dc=0: -1/c^2 + b2 + 2*b3*c + 3*b4*c^2 = 0, i.e. 3*b4*c^4 + 2*b3*c^3 + b2*c^2 - 1 = 0, an exact quartic solved by np.roots.""" b1, b2, b3, b4 = beta coeffs = [3 * b4, 2 * b3, b2, 0.0, -1.0] roots = np.roots(coeffs) real_pos = [r.real for r in roots if abs(r.imag) < 1e-9 * max(1.0, abs(r.real)) and r.real > 0] return min(real_pos) def agcl_ternary_minimum(beta, Ksp=KSP_AGCL, lo=-8.0, hi=2.0, iters=200): """Vectorised ternary search for the minimiser of x(c) in log10(c), valid because x is convex in c and log is monotone. Works for scalar or array beta (shape (...,4)).""" beta = np.asarray(beta, dtype=float) shape = beta.shape[:-1] lo_a = np.full(shape, lo) hi_a = np.full(shape, hi) def xval(logc): c = 10.0 ** logc b1 = beta[..., 0]; b2 = beta[..., 1]; b3 = beta[..., 2]; b4 = beta[..., 3] return Ksp / c + Ksp * b1 + Ksp * b2 * c + Ksp * b3 * c ** 2 + Ksp * b4 * c ** 3 for _ in range(iters): m1 = lo_a + (hi_a - lo_a) / 3.0 m2 = hi_a - (hi_a - lo_a) / 3.0 f1 = xval(m1) f2 = xval(m2) take_right = f1 > f2 lo_a = np.where(take_right, m1, lo_a) hi_a = np.where(take_right, hi_a, m2) logc_min = 0.5 * (lo_a + hi_a) return logc_min def banner(s): print("\n" + "=" * 78) print(s) print("=" * 78) def sub(s): print("\n--- %s ---" % s) def fmt_svg_polyline(xs, ys, cls, extra=""): pts = " ".join("%.2f,%.2f" % (x, y) for x, y in zip(xs, ys)) print('' % (cls, (" " + extra) if extra else "", pts)) def main(): banner("SOLUBILITY EQUILIBRIA: THE COMMON ION EFFECT, DONE EXACTLY") print("Seed: %d (numpy default_rng, PCG64)" % SEED) print("numpy version: %s" % np.__version__) print("Python version: %s" % sys.version.split()[0]) # ========================================================================= banner("PART 1. VALIDATION") # ========================================================================= sub("V1. Mass balance, equilibrium expressions and Ksp: residuals at every grid point") worst = {"baio3": 0.0, "caf2": 0.0, "agcl_ksp": 0.0, "agcl_mass_ag": 0.0, "agcl_mass_cl": 0.0, "agcl_beta": 0.0} n_checked = 0 for C in np.geomspace(1e-8, 5.0, 24): for pH in np.linspace(0.0, 14.0, 8): n_checked += 1 r = solve_baio3(C, pH) lhs = r["Ba"] * r["IO3"] ** 2 worst["baio3"] = max(worst["baio3"], abs(lhs - KSP_BAIO3) / KSP_BAIO3) mb = (C + 2 * r["x"]) - (r["IO3"] + r["HIO3"]) worst["baio3"] = max(worst["baio3"], abs(mb) / max(C + 2 * r["x"], 1e-300)) r2 = solve_caf2(C, pH) lhs2 = r2["Ca"] * r2["F"] ** 2 worst["caf2"] = max(worst["caf2"], abs(lhs2 - KSP_CAF2) / KSP_CAF2) mb2 = (C + 2 * r2["x"]) - (r2["F"] + r2["HF"]) worst["caf2"] = max(worst["caf2"], abs(mb2) / max(C + 2 * r2["x"], 1e-300)) r3 = solve_agcl(C, pH) worst["agcl_ksp"] = max(worst["agcl_ksp"], abs(r3["Ag"] * r3["Cl"] - KSP_AGCL) / KSP_AGCL) total_ag = r3["Ag"] + r3["AgCl0"] + r3["AgCl2"] + r3["AgCl3"] + r3["AgCl4"] worst["agcl_mass_ag"] = max(worst["agcl_mass_ag"], abs(total_ag - r3["x"]) / max(r3["x"], 1e-300)) total_cl_bound = r3["AgCl0"] + 2 * r3["AgCl2"] + 3 * r3["AgCl3"] + 4 * r3["AgCl4"] total_cl_free_family = r3["Cl"] / r3["f"] cl_check = C + r3["x"] - (total_cl_free_family + total_cl_bound) worst["agcl_mass_cl"] = max(worst["agcl_mass_cl"], abs(cl_check) / max(C + r3["x"], 1e-300)) for n, key in zip((1, 2, 3, 4), ("AgCl0", "AgCl2", "AgCl3", "AgCl4")): bn_check = r3[key] / (r3["Ag"] * r3["Cl"] ** n) - BETA_AGCL[n - 1] worst["agcl_beta"] = max(worst["agcl_beta"], abs(bn_check) / BETA_AGCL[n - 1]) print("Grid checked: %d (C, pH) points per salt, %d total solves." % (n_checked, n_checked * 3)) for k, v in worst.items(): print(" worst relative residual, %-13s: %.3e (%.1f eps)" % (k, v, v / EPS if v > 0 else 0.0)) sub("V1b. Charge balance closure: two independent derivations of the same titrant, compared") # For a salt MA_nu dissolved with its common ion added as a fully # dissociated NaA salt (C mol/L) and pH held by an external strong # monoprotic acid or base added at Ca mol/L (Ca<0 meaning base), Ca can # be obtained two ways that are algebraically identical only if every # mass balance and the full electroneutrality equation both hold: # Ca_A, from full charge balance solved for the spectator titrant ion; # Ca_B, from the proton condition (TOTH) referenced to Na+, the anion # at its zero-proton level, and water, i.e. Ca_B = [H+]-[OH-] # plus whatever protons are parked on the protonated anion. # A nonzero residual here would mean the model's mass balance and its # charge balance disagree; agreement to machine precision, at every grid # point, is the check. worst_q = {"baio3": 0.0, "caf2": 0.0, "agcl": 0.0} for C in np.geomspace(1e-8, 5.0, 12): for pH in np.linspace(0.2, 13.8, 10): h = 10.0 ** (-pH); oh = KW / h r = solve_baio3(C, pH) Ca_A = 2 * r["Ba"] + C + h - r["IO3"] - oh Ca_B = h - oh + r["HIO3"] worst_q["baio3"] = max(worst_q["baio3"], abs(Ca_A - Ca_B) / max(abs(Ca_B), 1e-300)) r2 = solve_caf2(C, pH) Ca_A2 = 2 * r2["Ca"] + C + h - r2["F"] - oh Ca_B2 = h - oh + r2["HF"] worst_q["caf2"] = max(worst_q["caf2"], abs(Ca_A2 - Ca_B2) / max(abs(Ca_B2), 1e-300)) r3 = solve_agcl(C, pH) HCl_conc = r3["Cl"] * (1.0 - r3["f"]) / r3["f"] Ca_A3 = C + r3["Ag"] + h - r3["Cl"] - r3["AgCl2"] - 2 * r3["AgCl3"] - 3 * r3["AgCl4"] - oh Ca_B3 = h - oh + HCl_conc denom3 = max(abs(Ca_B3), h, 1e-300) worst_q["agcl"] = max(worst_q["agcl"], abs(Ca_A3 - Ca_B3) / denom3) for k, v in worst_q.items(): print(" worst |Ca_A - Ca_B| / |Ca_B|, %-8s: %.3e (%.1f eps)" % (k, v, v / EPS if v > 0 else 0.0)) sub("V2. Closed-form cubic (np.roots) vs bisection solver, Ba(IO3)2 and CaF2") print("%-10s %-10s %-16s %-16s %-14s" % ("salt", "C (M)", "bisection x", "cubic root x", "rel. diff")) for C in (0.0, 1.0e-4, 1.0e-2, 1.0): xb = solve_baio3(C, 7.0)["x"] xc = baio3_cubic_root(C, 7.0) print("%-10s %-10.2e %-16.10e %-16.10e %-14.3e" % ("Ba(IO3)2", C, xb, xc, abs(xb - xc) / xc)) for C in (0.0, 1.0e-4, 1.0e-2, 1.0): xb = solve_caf2(C, 7.0)["x"] xc = caf2_cubic_root(C, 7.0) print("%-10s %-10.2e %-16.10e %-16.10e %-14.3e" % ("CaF2", C, xb, xc, abs(xb - xc) / xc)) sub("V3. AgCl pure-water solubility: solver vs closed form vs literature") x0 = solve_agcl(0.0, 7.0)["x"] closed = np.sqrt(KSP_AGCL) lit = 1.33e-5 print("Solver (full coupled system, C=0): %.6e mol/L" % x0) print("Closed form sqrt(Ksp) (no complexes): %.6e mol/L (rel. diff from solver: %.3e)" % (closed, abs(x0 - closed) / closed)) print("Literature value (e.g. common general-chemistry tables): %.3e mol/L (rel. diff: %.3e)" % (lit, abs(x0 - lit) / lit)) x0_baio3 = solve_baio3(0.0, 7.0)["x"] closed_baio3 = (KSP_BAIO3 / 4.0) ** (1.0 / 3.0) lit_baio3 = 7.32e-4 print("\nBa(IO3)2 pure-water solubility, solver: %.6e mol/L" % x0_baio3) print("Closed form (Ksp/4)^(1/3): %.6e mol/L (rel. diff: %.3e)" % (closed_baio3, abs(x0_baio3 - closed_baio3) / closed_baio3)) print("Commonly quoted value: %.3e mol/L (rel. diff: %.3e)" % (lit_baio3, abs(x0_baio3 - lit_baio3) / lit_baio3)) sub("V4. AgCl solubility minimum: internal check, then against the literature") logc_min_scalar = agcl_ternary_minimum(BETA_AGCL[None, :], lo=-8, hi=2, iters=200)[0] c_min_ternary = 10.0 ** logc_min_scalar c_min_quartic = agcl_quartic_minimum(BETA_AGCL) print("Free [Cl-] at minimum, ternary search on solver's own x(c): %.6e M" % c_min_ternary) print("Free [Cl-] at minimum, exact quartic root (np.roots): %.6e M" % c_min_quartic) print("Relative difference (internal consistency of two independent numerical routes): %.3e" % (abs(c_min_ternary - c_min_quartic) / c_min_quartic)) C_at_min, x_at_min, _ = agcl_forward(c_min_quartic, 1.0, KSP_AGCL, BETA_AGCL) print("Corresponding added common-ion concentration at the minimum: C = %.4e M NaCl" % C_at_min) print("Solubility at the minimum: x = %.4e M" % x_at_min) lit_lo, lit_hi = 1.0e-3, 1.0e-2 inside = lit_lo <= C_at_min <= lit_hi print("Range widely reported in analytical-chemistry teaching for where this minimum falls:") print(" %.1e to %.1e M Cl- (order-of-magnitude description, e.g. [1][2])" % (lit_lo, lit_hi)) print("Club's computed minimum %s that range (log10 distance from nearer edge: %.2f decades)." % ("falls inside" if inside else "falls outside", 0.0 if inside else min(abs(np.log10(C_at_min / lit_lo)), abs(np.log10(C_at_min / lit_hi))))) print("HONEST NOTE: this is a comparison to a range quoted to one significant figure in secondary") print("sources, not to a single re-verified primary number; see the article SS10 and V5 below.") sub("V5. Monte Carlo: how much does the reversal depend on the uncertain formation constants?") N_MC = 100000 SIGMA_LOGBETA = 0.05 draws = RNG.normal(loc=LOGBETA_AGCL, scale=SIGMA_LOGBETA, size=(N_MC, 4)) draws = np.maximum(draws, 0.1) # keep positive and away from pathological values beta_draws = 10.0 ** draws logc_min_mc = agcl_ternary_minimum(beta_draws, lo=-8, hi=2, iters=120) c_min_mc = 10.0 ** logc_min_mc # forward-map each draw's c_min to the corresponding added concentration C_min b1 = beta_draws[:, 0]; b2 = beta_draws[:, 1]; b3 = beta_draws[:, 2]; b4 = beta_draws[:, 3] C_min_mc = c_min_mc - KSP_AGCL / c_min_mc + KSP_AGCL * b2 * c_min_mc \ + 2 * KSP_AGCL * b3 * c_min_mc ** 2 + 3 * KSP_AGCL * b4 * c_min_mc ** 3 # enhancement ratio: solubility at a fixed high chloride (5 M, near NaCl's # own solubility ceiling) divided by solubility at the minimum, per draw c_hi = 5.0 x_hi_mc = KSP_AGCL / c_hi + KSP_AGCL * b1 + KSP_AGCL * b2 * c_hi + KSP_AGCL * b3 * c_hi ** 2 + KSP_AGCL * b4 * c_hi ** 3 x_min_mc = KSP_AGCL / c_min_mc + KSP_AGCL * b1 + KSP_AGCL * b2 * c_min_mc + KSP_AGCL * b3 * c_min_mc ** 2 + KSP_AGCL * b4 * c_min_mc ** 3 enhancement_mc = x_hi_mc / x_min_mc mean_logC = np.mean(np.log10(C_min_mc)) se_logC = np.std(np.log10(C_min_mc), ddof=1) / np.sqrt(N_MC) mean_enh = np.mean(enhancement_mc) se_enh = np.std(enhancement_mc, ddof=1) / np.sqrt(N_MC) frac_reversal = np.mean(enhancement_mc > 1.0) print("N = %d draws, log10(beta_n) ~ Normal(reported, %.2f) independently for n=1..4" % (N_MC, SIGMA_LOGBETA)) print("log10(C_min / M): mean = %.4f, SE = %.4f => C_min = %.3e M (a factor-of-%.2f 1-SE band)" % (mean_logC, se_logC, 10 ** mean_logC, 10 ** se_logC)) print("Enhancement ratio x(5 M)/x(minimum): mean = %.1f, SE = %.2f" % (mean_enh, se_enh)) print("Fraction of draws where 5 M chloride dissolves MORE AgCl than the minimum does: %.4f%%" % (100 * frac_reversal)) print("(A value near 100%% means the reversal's existence is essentially certain under this") print(" assumed uncertainty, even though its exact location and size are not.)") # convergence trace: running mean/SE of log10(C_min) as a function of N used checkpoints = np.unique(np.geomspace(10, N_MC, 40).astype(int)) running_mean = [] running_se = [] logC = np.log10(C_min_mc) cs = np.cumsum(logC) css = np.cumsum(logC ** 2) for n in checkpoints: m = cs[n - 1] / n var = css[n - 1] / n - m ** 2 var = max(var, 0.0) se = np.sqrt(var / n) running_mean.append(m) running_se.append(se) print("\nConvergence trace (n draws, running mean log10 C_min, running SE):") for n, m, se in zip(checkpoints[::4], running_mean[::4], running_se[::4]): print(" n=%7d mean=%.4f SE=%.4f" % (n, m, se)) # ========================================================================= banner("PART 2. THE THREE SALTS ACROSS THE GRID") # ========================================================================= sub("Table data: summary per salt") naive_baio3 = lambda C: KSP_BAIO3 / max(C, 1e-300) ** 2 naive_caf2 = lambda C: KSP_CAF2 / max(C, 1e-300) ** 2 naive_agcl = lambda C: KSP_AGCL / max(C, 1e-300) print("%-10s %-12s %-10s %-14s %-14s %-16s" % ("salt", "Ksp", "n(anion)", "pure-water x", "C at min (M)", "max enhancement*")) print("%-10s %-12.3e %-10d %-14.4e %-14s %-16s" % ("Ba(IO3)2", KSP_BAIO3, 2, x0_baio3, "none (monotone)", "1.00 (never reverses)")) x0_caf2 = solve_caf2(0.0, 7.0)["x"] print("%-10s %-12.3e %-10d %-14.4e %-14s %-16s" % ("CaF2", KSP_CAF2, 2, x0_caf2, "none (monotone)", "1.00 (never reverses)")) print("%-10s %-12.3e %-10d %-14.4e %-14.4e %-16.1f" % ("AgCl", KSP_AGCL, 1, x0, C_at_min, mean_enh)) print("* solubility at 5 M common ion divided by solubility at the minimum (AgCl only; the") print(" other two salts have no minimum because they have no complexation channel modelled).") sub("Figure 1 data: AgCl solubility vs added [Cl-], exact solver and naive rule, pH=7") C_grid_agcl = np.geomspace(1e-7, 6.0, 55) x_exact_agcl = np.array([solve_agcl(C, 7.0)["x"] for C in C_grid_agcl]) x_naive_agcl = np.array([naive_agcl(C) for C in C_grid_agcl]) print("%-4s %-14s %-14s %-14s" % ("i", "C_added (M)", "x_exact (M)", "x_naive (M)")) for i in range(0, len(C_grid_agcl), 3): print("%-4d %-14.6e %-14.6e %-14.6e" % (i, C_grid_agcl[i], x_exact_agcl[i], x_naive_agcl[i])) sub("Figure 2 data: CaF2 solubility vs pH at fixed C = 1e-3 M NaF") pH_grid = np.linspace(0.0, 8.0, 45) C_fixed_f = 1.0e-3 x_caf2_ph = np.array([solve_caf2(C_fixed_f, pH)["x"] for pH in pH_grid]) x_caf2_naive_ph = np.full_like(pH_grid, naive_caf2(C_fixed_f)) x_agcl_ph = np.array([solve_agcl(1.0e-3, pH)["x"] for pH in pH_grid]) # contrast: AgCl vs pH, same C print("%-4s %-8s %-14s %-14s %-14s" % ("i", "pH", "CaF2 x (M)", "CaF2 naive", "AgCl x (M), C=1e-3")) for i in range(0, len(pH_grid), 2): print("%-4d %-8.2f %-14.6e %-14.6e %-14.6e" % (i, pH_grid[i], x_caf2_ph[i], x_caf2_naive_ph[i], x_agcl_ph[i])) print("AgCl range over this pH sweep: max/min = %.6f (i.e. %.4f%% variation)" % (x_agcl_ph.max() / x_agcl_ph.min(), 100 * (x_agcl_ph.max() / x_agcl_ph.min() - 1))) sub("Figure 3 data: relative error of the naive textbook rule for AgCl, vs [Cl-]") rel_err_agcl = (x_naive_agcl - x_exact_agcl) / x_exact_agcl print("%-4s %-14s %-14s" % ("i", "C_added (M)", "rel. error (naive-exact)/exact")) for i in range(0, len(C_grid_agcl), 3): print("%-4d %-14.6e %-14.6e" % (i, C_grid_agcl[i], rel_err_agcl[i])) # where does the naive rule flip sign of the trend, vs where does it cross 100% error? sign_flip_idx = np.argmin(np.abs(C_grid_agcl - C_at_min)) over100 = C_grid_agcl[np.where(rel_err_agcl > 1.0)[0]] print("Naive rule crosses +100%% relative error at C >= %.4e M" % (over100.min() if len(over100) else float("nan"))) print("At C = 5 M (near-saturated NaCl), naive predicts x = %.3e M; exact solver gives x = %.3e M" % (naive_agcl(5.0), solve_agcl(5.0, 7.0)["x"])) print(" naive is wrong by a factor of %.1f, and wrong in DIRECTION below C=%.3e M" % (solve_agcl(5.0, 7.0)["x"] / naive_agcl(5.0), C_at_min)) sub("Figure 4 data: two-ligand regime map, Ag+/Cl-/NH3, sign of d(ln x)/d(ln[Cl-])") def agcl_nh3_solve(C_cl, C_nh3, pH=7.0): h = 10.0 ** (-pH) f = alpha_free(h, KA_HCL) bn1, bn2 = BETA_AGNH3 def resid(logc): c = 10.0 ** logc a = KSP_AGCL / c # solve NH3 free concentration s from its own mass balance at this c: # C_nh3 = s + a*bn1*s + 2*a*bn2*s^2 -> quadratic in s A = 2 * a * bn2 B = 1 + a * bn1 Cc = -C_nh3 if A < 1e-300: s = -Cc / B else: s = (-B + np.sqrt(max(B * B - 4 * A * Cc, 0.0))) / (2 * A) s = max(s, 0.0) b1, b2, b3, b4 = BETA_AGCL x = a * (1 + b1 * c + b2 * c ** 2 + b3 * c ** 3 + b4 * c ** 4 + bn1 * s + bn2 * s ** 2) C_cl_pred = c / f - KSP_AGCL / c + KSP_AGCL * b2 * c + 2 * KSP_AGCL * b3 * c ** 2 + 3 * KSP_AGCL * b4 * c ** 3 \ - a * (bn1 * s + bn2 * s ** 2) return C_cl_pred - C_cl, x, s r0 = resid(-14.0)[0] r1 = resid(3.0)[0] lo_b, hi_b = -14.0, 3.0 if r0 > 0: lo_b = -30.0 for _ in range(140): mid = 0.5 * (lo_b + hi_b) if resid(mid)[0] <= 0: lo_b = mid else: hi_b = mid c_sol = 10.0 ** (0.5 * (lo_b + hi_b)) _, x_sol, s_sol = resid(0.5 * (lo_b + hi_b)) return x_sol, c_sol, s_sol n_grid = 26 Cl_axis = np.geomspace(1e-6, 5.0, n_grid) NH3_axis = np.geomspace(1e-6, 2.0, n_grid) regime = np.zeros((n_grid, n_grid)) dl = 0.02 for i, cnh3 in enumerate(NH3_axis): for j, ccl in enumerate(Cl_axis): x_lo, _, _ = agcl_nh3_solve(ccl * 10 ** (-dl), cnh3) x_hi, _, _ = agcl_nh3_solve(ccl * 10 ** (dl), cnh3) slope = (np.log(x_hi) - np.log(x_lo)) / (2 * dl * LN10) regime[i, j] = slope print("Grid: %d x %d points, [Cl-] in [1e-6, 5] M, [NH3] in [1e-6, 2] M, pH=7." % (n_grid, n_grid)) n_reverse = np.sum(regime > 0) print("Points where d(ln x)/d(ln[Cl-]) > 0 (rule reverses): %d of %d (%.1f%%)" % (n_reverse, n_grid * n_grid, 100 * n_reverse / (n_grid * n_grid))) print("Sample rows (NH3 index, [NH3] M, then sign pattern across the Cl- axis, - = holds, + = reverses):") for i in range(0, n_grid, 5): pattern = "".join("+" if v > 0 else "-" for v in regime[i]) print(" i=%2d [NH3]=%.3e %s" % (i, NH3_axis[i], pattern)) default_x, default_c, default_s = agcl_nh3_solve(1.0e-3, 0.0) print("Default interactive-model condition ([Cl-]=1e-3 M added, [NH3]=0): x=%.4e M, matches" % default_x) print("the plain-AgCl solver at the same C to within %.3e (relative)." % (abs(default_x - solve_agcl(1.0e-3, 7.0)["x"]) / solve_agcl(1.0e-3, 7.0)["x"])) # ========================================================================= banner("PART 3. SVG COORDINATE DATA FOR THE ARTICLE'S FIGURES") # ========================================================================= # Every coordinate below is a direct, printed linear transform of the # numbers computed above into a fixed viewBox, so what appears in the # article's inline SVG is this printed output, copied. def logmap(v, vmin, vmax, pmin, pmax): return pmin + (np.log10(v) - np.log10(vmin)) / (np.log10(vmax) - np.log10(vmin)) * (pmax - pmin) def linmap(v, vmin, vmax, pmin, pmax): return pmin + (v - vmin) / (vmax - vmin) * (pmax - pmin) sub("Fig 1 SVG points: log10(C) -> x-pixel [80,580], log10(x) -> y-pixel [280,34]") X0, X1, Y0, Y1 = 80.0, 580.0, 280.0, 34.0 YLO, YHI = 1e-7, 1e-1 px1 = logmap(C_grid_agcl, 1e-7, 6.0, X0, X1) py_exact = logmap(np.clip(x_exact_agcl, YLO, YHI), YLO, YHI, Y0, Y1) py_naive = logmap(np.clip(x_naive_agcl, YLO, 1e2), YLO, YHI, Y0, Y1) for i in range(0, len(px1), 3): print("EX %7.2f,%7.2f" % (px1[i], py_exact[i])) for i in range(0, len(px1), 3): print("NV %7.2f,%7.2f" % (px1[i], py_naive[i])) px_min = logmap(np.array([C_at_min]), 1e-7, 6.0, X0, X1)[0] py_min = logmap(np.array([x_at_min]), YLO, YHI, Y0, Y1)[0] print("MIN %7.2f,%7.2f" % (px_min, py_min)) sub("Fig 2 SVG points: pH -> x-pixel [80,580], log10(x, M) -> y-pixel [280,34]") px2 = linmap(pH_grid, 0.0, 8.0, X0, X1) py_caf2 = logmap(x_caf2_ph, 1e-5, 3e-2, Y0, Y1) py_caf2_naive = logmap(x_caf2_naive_ph, 1e-5, 3e-2, Y0, Y1) for i in range(len(px2)): print("CF %7.2f,%7.2f" % (px2[i], py_caf2[i])) for i in range(len(px2)): print("NV %7.2f,%7.2f" % (px2[i], py_caf2_naive[i])) print("(AgCl over this same pH sweep is flat to %.4f%% and is drawn as a plain reference line," % (100 * (x_agcl_ph.max() / x_agcl_ph.min() - 1))) print(" not plotted point by point.)") sub("Fig 3 SVG points: log10(C) -> x-pixel [80,580], rel. error (%) -> y-pixel [280,34], clipped +-300%") py_err = linmap(np.clip(rel_err_agcl * 100, -100, 300), -100, 300, Y0, Y1) for i in range(0, len(px1), 3): print("ER %7.2f,%7.2f" % (px1[i], py_err[i])) py_zero = linmap(np.array([0.0]), -100, 300, Y0, Y1)[0] print("ZERO_Y %.2f" % py_zero) sub("Fig 4 SVG cells: regime map grid, [Cl-] -> x [80,580], [NH3] -> y [280,34], flag +/-") gx = logmap(Cl_axis, 1e-6, 5.0, X0, X1) gy = logmap(NH3_axis, 1e-6, 2.0, Y0, Y1) cellw = (X1 - X0) / n_grid cellh = (Y0 - Y1) / n_grid for i in range(n_grid): for j in range(n_grid): flag = "R" if regime[i, j] > 0 else "H" print("CELL %s %7.2f %7.2f %5.2f %5.2f" % (flag, gx[j] - cellw / 2, gy[i] - cellh / 2, cellw, cellh)) gdx = logmap(np.array([1.0e-3]), 1e-6, 5.0, X0, X1)[0] gdy = logmap(np.array([max(1e-6, 1e-6)]), 1e-6, 2.0, Y0, Y1)[0] print("MARK %7.2f %7.2f" % (gdx, gdy)) sub("Fig 5 SVG points: MC convergence, n (log) -> x-pixel [80,580], running SE -> y-pixel [280,34]") px5 = logmap(np.array(checkpoints, dtype=float), 10, N_MC, X0, X1) py5 = linmap(np.array(running_se), 0.0, max(running_se), Y0, Y1) for i in range(len(checkpoints)): print("SE %7.2f,%7.2f" % (px5[i], py5[i])) py5_mean = linmap(np.array(running_mean), min(running_mean) - 0.02, max(running_mean) + 0.02, Y0, Y1) for i in range(len(checkpoints)): print("MEAN %7.2f,%7.2f" % (px5[i], py5_mean[i])) sub("Additional headline numbers for the article text") exact_5M = solve_agcl(5.0, 7.0)["x"] naive_5M = naive_agcl(5.0) print("Deterministic (nominal-constant) enhancement ratio x(5 M)/x(minimum): %.1f" % (exact_5M / x_at_min)) print("Naive rule at 5 M: %.4e M; exact: %.4e M; naive is low by a factor of %.3e" % (naive_5M, exact_5M, exact_5M / naive_5M)) # band where the naive rule is within +-10% / +-20% of the exact solver within10 = C_grid_agcl[np.abs(rel_err_agcl) <= 0.10] within20 = C_grid_agcl[np.abs(rel_err_agcl) <= 0.20] if len(within10): print("Naive rule within +-10%% of exact for C in [%.3e, %.3e] M" % (within10.min(), within10.max())) if len(within20): print("Naive rule within +-20%% of exact for C in [%.3e, %.3e] M" % (within20.min(), within20.max())) print("Pure-water solubility (no common ion): %.4e M -- the naive rule is undefined there (C=0)." % x0) print("Ratio of minimum-concentration to pure-water solubility itself: %.2f x" % (C_at_min / x0)) # CaF2: pH at which solubility first exceeds 2x, 10x, 100x the neutral-pH value x_caf2_neutral = solve_caf2(C_fixed_f, 7.0)["x"] for mult in (2.0, 10.0, 100.0): idx = np.where(x_caf2_ph >= mult * x_caf2_neutral)[0] if len(idx): # idx is ordered from pH=0 upward; the LAST index where still >= mult marks the crossing pH (since curve is monotone decreasing in pH) i_cross = idx.max() print("CaF2 solubility (C=1e-3 M NaF) exceeds %gx its high-pH value below pH ~ %.2f" % (mult, pH_grid[i_cross])) print("CaF2 at pH 0: solubility is %.1fx the high-pH (pH>=7) plateau of %.4e M." % (x_caf2_ph[0] / x_caf2_neutral, x_caf2_neutral)) sub("Species breakdown: at the minimum, and at 5 M chloride") r_min = solve_agcl(C_at_min, 7.0) r_5M = solve_agcl(5.0, 7.0) for label, r in (("at the minimum (C=%.4e M)" % C_at_min, r_min), ("at 5 M NaCl", r_5M)): tot = r["Ag"] + r["AgCl0"] + r["AgCl2"] + r["AgCl3"] + r["AgCl4"] print(" %s: total x=%.4e M" % (label, tot)) for name in ("Ag", "AgCl0", "AgCl2", "AgCl3", "AgCl4"): print(" %-6s = %.4e M (%.2f%% of total)" % (name, r[name], 100 * r[name] / tot)) print("Exact solver x(5 M) / pure-water solubility x(0): %.1fx" % (solve_agcl(5.0, 7.0)["x"] / x0)) sub("V4b (objection check): how wrong would log beta2 need to be for the reversal to vanish " "within the real solubility range of NaCl (about 6.1 mol/L at 25C)?") # Reuse the quartic minimiser condition 3*b4*c^4 + 2*b3*c^3 + b2*c^2 - 1 = 0 solved for b2 # at c fixed at NaCl's own saturation limit: a minimum located there or beyond is # observationally indistinguishable from "no reversal in any real solution." c_sat_nacl = 6.1 b1n, b2n, b3n, b4n = BETA_AGCL b2_breakeven = (1.0 - 3 * b4n * c_sat_nacl ** 4 - 2 * b3n * c_sat_nacl ** 3) / (c_sat_nacl ** 2) print("Nominal log10(beta2) = %.2f (beta2 = %.3e)" % (LOGBETA_AGCL[1], BETA_AGCL[1])) if b2_breakeven > 0: print("log10(beta2) would need to fall to %.2f (a change of %.2f decades) to push the" % (np.log10(b2_breakeven), LOGBETA_AGCL[1] - np.log10(b2_breakeven))) print("minimum out to NaCl's own saturation limit, %.1f M -- i.e. to erase the reversal" % c_sat_nacl) print("from every real solution of the salt.") else: print("No positive beta2 pushes the minimum that far with beta3, beta4 held fixed: the") print("higher complexes alone already guarantee a minimum inside the real concentration range.") banner("DONE") print("Total wall-clock time: %.2f s" % (time.time() - T0)) if __name__ == "__main__": main()