#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ The Reaction That Gives You the Wrong Product Because You Were Impatient Science Journaling Club, Volume 1, Issue 3, Spring 2025 Theme: "Equilibrium Chemistry by Computer" ================================================================================ THE QUESTION ================================================================================ A starting material A can become either of two products through two separate, reversible, one-step reactions: A <=> B (rate constants k1 forward, k1r reverse) "kinetic" product A <=> C (rate constants k2 forward, k2r reverse) "thermodynamic" product B and C do not interconvert directly. This is the standard idealisation behind every "kinetic versus thermodynamic control" example a chemistry course hands out (1,2- vs 1,4-addition to a diene, endo vs exo Diels-Alder adducts, kinetic vs thermodynamic enolates), and it is honest about what such courses usually leave vague: a number. If B forms faster (k1 > k2) but C is the more stable product (its equilibrium constant K2 = k2/k2r exceeds K1 = k1/k1r), the mixture is mostly B at short times and mostly C at long times. Somewhere in between, the two curves cross. How long is that wait, in seconds a person can actually experience, and what two numbers -- a barrier gap and a stability gap -- decide whether the wait is a coffee break or longer than the solar system will exist? ================================================================================ THE MODEL ================================================================================ Each rate constant obeys the Arrhenius law, k(T) = A_pref * exp(-Ea / (R T)) (1) with a single shared pre-exponential factor A_pref for all four steps (a simplification named below). The forward barriers Ea1, Ea2 are free parameters. The two reverse barriers are NOT free: they are fixed by the requested product stability through the thermodynamic identity Keq = kf / kr = exp(-DeltaG / (R T)) (2) which, together with kf = A_pref*exp(-Eaf/RT) and taking the reverse step to share the same A_pref, forces Ea_reverse = Ea_forward - DeltaG (3) DeltaG_B and DeltaG_C are the free energies of B and C relative to A (both taken negative, i.e. both products more stable than the starting material, in the sweeps where that is the intended scenario). Equation (3) is exactly the statement that thermodynamics and kinetics must agree at equilibrium: it is not an extra assumption glued on afterwards, it is what makes k1/k1r equal the Keq that DeltaG_B implies, to machine precision, at every temperature. The rate equations are dA/dt = -(k1+k2) A + k1r B + k2r C dB/dt = k1 A - k1r B dC/dt = k2 A - k2r C (4) with A(0) = A0, B(0) = C(0) = 0, and A+B+C = A0 for all t (mass conservation, checked numerically below rather than assumed). Two solvers are used, deliberately, and cross-checked against each other: rk4_integrate a fixed-step, explicit fourth-order Runge-Kutta integrator, written from scratch (no ODE library), applied to (4) directly. This is "the computation is the experiment": we watch concentrations evolve the way a numerical experiment would, one step at a time. exact_linear Equations (4) are linear with constant coefficients at fixed T. Eliminating C = A0 - A - B leaves a 2x2 linear system dy/dt = M y + b, solved exactly by eigendecomposing M (numpy.linalg.eig) and exponentiating its eigenvalues. This is not a second numerical method in the sense of competing with RK4 -- it is closed form, limited only by floating point and by numpy's eigensolver -- and it lets us evaluate the crossover time for hundreds of parameter combinations, some of which do not cross until 10^15 seconds or more, without ever taking a single small time step out that far. Section VALIDATION below checks rk4_integrate against exact_linear, and both against the textbook closed-form solution of the single reversible reaction A <=> B taken alone, at many time points, before either is trusted for the sweeps in section RESULTS. ================================================================================ ASSUMPTIONS, NAMED ================================================================================ * One shared pre-exponential factor A_pref = 1e13 s^-1 for all four steps (forward and reverse, both channels). Real reactions differ in A_pref by orders of magnitude because entropy of activation differs step to step; fixing it means every quoted stability DeltaG maps to a rate-constant RATIO exactly (equation 3) but the ABSOLUTE rate of any one step still depends on this choice. 1e13 s^-1 is the generic order of magnitude transition-state theory gives for a "normal" unimolecular step (k_B T/h at room temperature is 6e12 s^-1). * B and C never interconvert directly and neither reaction is autocatalytic or catalysed by anything else present. * Elementary, single-step, first-order kinetics in both directions. No intermediate, no pre-equilibrium complex, no diffusion limitation. * Constant temperature and volume; no side reactions, no solvent, no concentration dependence in the rate "constants" beyond concentration itself. * DeltaG here is treated as temperature-independent (no DeltaCp, no entropy term beyond what is folded into the shared A_pref assumption). Equation (3) is therefore an approximation once T is swept far from a single reference temperature; see the temperature-sweep discussion in the article for how much this could bend the numbers. * This is a model of a mechanism, not a measurement of any specific real reaction. No reagent was mixed, no reaction was run, and no yield in this file was measured on a real bench. Every number below is the printed output of the code in this file. Seed: none needed. Every quantity in this file is produced by solving deterministic differential equations, not by sampling; there is no random number generator anywhere in this script. The "convergence" figure means numerical convergence of the integrator with step size, not Monte Carlo convergence with sample size. Author: Science Journaling Club. Written for Python 3.12, using numpy only. """ import sys import time as _time import numpy as np R = 8.314462618 # J / (mol K), CODATA exact value under the SI redefinition # ------------------------------------------------------------------------ # Section 1: rate constants # ------------------------------------------------------------------------ def arrhenius(A_pref, Ea_J, T): """k(T) = A_pref * exp(-Ea/RT), Ea in J/mol, T in kelvin.""" return A_pref * np.exp(-Ea_J / (R * T)) def rate_constants(Ea1_kJ, Ea2_kJ, dGB_kJ, dGC_kJ, T, A_pref=1.0e13): """ Forward barriers Ea1, Ea2 are free parameters (kJ/mol). Reverse barriers are fixed by equation (3) so that k1/k1r = exp(-dGB/RT) exactly and k2/k2r = exp(-dGC/RT) exactly, at the given temperature T (kelvin). Returns (k1, k1r, k2, k2r) in s^-1. """ Ea1 = Ea1_kJ * 1000.0 Ea2 = Ea2_kJ * 1000.0 dGB = dGB_kJ * 1000.0 dGC = dGC_kJ * 1000.0 Ea1r = Ea1 - dGB Ea2r = Ea2 - dGC k1 = arrhenius(A_pref, Ea1, T) k1r = arrhenius(A_pref, Ea1r, T) k2 = arrhenius(A_pref, Ea2, T) k2r = arrhenius(A_pref, Ea2r, T) return k1, k1r, k2, k2r # ------------------------------------------------------------------------ # Section 2: the two solvers # ------------------------------------------------------------------------ def rhs(y, k1, k1r, k2, k2r): A, B, C = y dA = -(k1 + k2) * A + k1r * B + k2r * C dB = k1 * A - k1r * B dC = k2 * A - k2r * C return np.array([dA, dB, dC]) def rk4_integrate(k1, k1r, k2, k2r, A0, t_final, n_steps): """Fixed-step, explicit RK4 on the full 3-species system. Returns arrays t, A, B, C of length n_steps+1, including the initial condition.""" dt = t_final / n_steps y = np.array([A0, 0.0, 0.0]) t = 0.0 ts = np.empty(n_steps + 1) ys = np.empty((n_steps + 1, 3)) ts[0] = 0.0 ys[0] = y for i in range(1, n_steps + 1): k_a = rhs(y, k1, k1r, k2, k2r) k_b = rhs(y + 0.5 * dt * k_a, k1, k1r, k2, k2r) k_c = rhs(y + 0.5 * dt * k_b, k1, k1r, k2, k2r) k_d = rhs(y + dt * k_c, k1, k1r, k2, k2r) y = y + (dt / 6.0) * (k_a + 2 * k_b + 2 * k_c + k_d) t += dt ts[i] = t ys[i] = y return ts, ys[:, 0], ys[:, 1], ys[:, 2] def analytic_single_reaction(k1, k1r, A0, t): """Closed-form solution of dA/dt=-k1 A+k1r B, dB/dt=k1 A-k1r B alone (i.e. the textbook simple reversible first-order reaction, k2=k2r=0).""" Be = A0 * k1 / (k1 + k1r) B = Be * (1.0 - np.exp(-(k1 + k1r) * t)) A = A0 - B return A, B def exact_linear(k1, k1r, k2, k2r, A0, t): """ Exact solution of the full two-channel linear system by eigendecomposing the 2x2 matrix that remains after eliminating C = A0 - A - B: d/dt [A] [ -(k1+k2+k2r) k1r-k2r ] [A] [ k2r*A0 ] [B] = [ k1 -k1r ] [B] + [ 0 ] y(t) = y_ss + V diag(exp(lambda t)) V^-1 (y0 - y_ss), where y_ss is the fixed point of M y + b = 0. t may be a scalar or a numpy array; returns (A(t), B(t), C(t)) with the same shape as t. """ M = np.array([[-(k1 + k2 + k2r), k1r - k2r], [k1, -k1r]], dtype=float) b = np.array([k2r * A0, 0.0]) y0 = np.array([A0, 0.0]) y_ss = np.linalg.solve(M, -b) eigvals, V = np.linalg.eig(M) Vinv = np.linalg.inv(V) coeffs = Vinv @ (y0 - y_ss) # complex in general, real here (tree graph) t_arr = np.atleast_1d(np.asarray(t, dtype=float)) # exponent matrix, shape (len(t), 2) expo = np.exp(np.outer(t_arr, eigvals)) delta = (expo * coeffs[np.newaxis, :]) @ V.T # shape (len(t), 2) A = y_ss[0] + delta[:, 0].real B = y_ss[1] + delta[:, 1].real C = A0 - A - B if np.isscalar(t): return A[0], B[0], C[0] return A, B, C def rk4_single_step(k1, k1r, k2, k2r, A0, dt): """One RK4 step of size dt from (A0,0,0). Used only to probe the t->0 kinetic limit: unlike exact_linear (which subtracts an O(1) equilibrium value from an O(1) exponential term to recover an O(dt) result, and so loses precision to catastrophic cancellation once dt is many orders of magnitude below 1/k), a single small forward step never subtracts two large numbers -- it multiplies an O(1) rate by a small dt directly.""" y = np.array([A0, 0.0, 0.0]) k_a = rhs(y, k1, k1r, k2, k2r) k_b = rhs(y + 0.5 * dt * k_a, k1, k1r, k2, k2r) k_c = rhs(y + 0.5 * dt * k_b, k1, k1r, k2, k2r) k_d = rhs(y + dt * k_c, k1, k1r, k2, k2r) return y + (dt / 6.0) * (k_a + 2 * k_b + 2 * k_c + k_d) def find_peak(k1, k1r, k2, k2r, A0, species=1, t_lo=1e-2, t_hi=1e9, n_scan=4000): """Locate the time and value of the maximum of B(t) (species=1) or C(t) (species=2) by a coarse log-spaced scan followed by golden-section refinement in log(t). Returns (t_peak, value_peak).""" ts = np.geomspace(t_lo, t_hi, n_scan) A_s, B_s, C_s = exact_linear(k1, k1r, k2, k2r, A0, ts) y = B_s if species == 1 else C_s i = int(np.argmax(y)) lo = np.log(ts[max(i - 1, 0)]) hi = np.log(ts[min(i + 1, n_scan - 1)]) gr = (np.sqrt(5.0) - 1.0) / 2.0 a, b = lo, hi def val(logt): A_v, B_v, C_v = exact_linear(k1, k1r, k2, k2r, A0, np.exp(logt)) return B_v if species == 1 else C_v c = b - gr * (b - a) d = a + gr * (b - a) fc, fd = val(c), val(d) for _ in range(80): if fc > fd: b, d, fd = d, c, fc c = b - gr * (b - a) fc = val(c) else: a, c, fc = c, d, fd d = a + gr * (b - a) fd = val(d) t_peak = float(np.exp(0.5 * (a + b))) return t_peak, val(np.log(t_peak)) def crossover_time(k1, k1r, k2, k2r, A0, t_lo=1e-8, t_hi=1e30, tol=1e-10): """ First time t at which C(t) >= B(t), found by bisection on f(t) = C(t) - B(t) in log(t). Returns (t_star, status): status = "instant" if C(t_lo) already >= B(t_lo) (thermo product also forms faster: no meaningful "kinetic" head start) status = "never" if C never catches B even at t_hi (checked at the equilibrium point directly, not just at t_hi) status = "crossed" otherwise, t_star is the bisected crossing time """ _, B_eq, C_eq = exact_linear(k1, k1r, k2, k2r, A0, 1e60) if C_eq <= B_eq * (1 + 1e-12): return np.inf, "never" _, B_lo, C_lo = exact_linear(k1, k1r, k2, k2r, A0, t_lo) if C_lo >= B_lo: return t_lo, "instant" lo, hi = np.log(t_lo), np.log(t_hi) for _ in range(200): mid = 0.5 * (lo + hi) _, Bm, Cm = exact_linear(k1, k1r, k2, k2r, A0, np.exp(mid)) if Cm >= Bm: hi = mid else: lo = mid if hi - lo < tol: break return np.exp(0.5 * (lo + hi)), "crossed" # ------------------------------------------------------------------------ # Section 3: pretty-printing helpers # ------------------------------------------------------------------------ def fmt_time(t): """Render a duration in seconds using whatever unit reads most naturally, all the way out past the age of the universe.""" if not np.isfinite(t): return "never (thermodynamic product is not actually more stable)" units = [ (1.0, "s"), (60.0, "min"), (3600.0, "hr"), (86400.0, "day"), (31557600.0, "yr"), ] AGE_UNIVERSE = 4.35e17 # s, ~13.8 Gyr if t > 1000 * AGE_UNIVERSE: return "%.3e s (%.3e x the age of the universe)" % (t, t / AGE_UNIVERSE) if t > AGE_UNIVERSE: return "%.3e s (%.2f x the age of the universe)" % (t, t / AGE_UNIVERSE) if t < 1.0: return "%.4g s" % t best = units[0] for u in units: if t >= u[0]: best = u val = t / best[0] return "%.4g s (%.4g %s)" % (t, val, best[1]) def hdr(title): print() print("=" * 78) print(title) print("=" * 78) # ========================================================================== # MAIN # ========================================================================== def main(): t_start = _time.time() print("The Reaction That Gives You the Wrong Product Because You Were Impatient") print("Science Journaling Club -- kinetics-vs-equilibrium.py") print("numpy %s, python %s" % (np.__version__, sys.version.split()[0])) print("No random number generator is used anywhere in this script.") A_PREF = 1.0e13 # s^-1, shared pre-exponential factor, section ASSUMPTIONS A0 = 1.0 # mol/L, arbitrary concentration scale (ratios are what matter) # Baseline / "headline" reaction, used for Figures 1, 2, 5 and the # temperature sweep, and reproduced by the interactive model's defaults. EA1_BASE, EA2_BASE = 90.0, 100.0 # kJ/mol forward barriers DGB_BASE, DGC_BASE = -15.0, -40.0 # kJ/mol product stabilities T_REF = 298.15 # K # ====================================================================== hdr("VALIDATION 1 of 5 -- RK4 vs the textbook A<=>B closed form") # ====================================================================== print("Single reversible reaction only (k2 = k2r = 0), baseline forward/") print("reverse rate constants at T = %.2f K, A0 = %.3f mol/L." % (T_REF, A0)) k1, k1r, _, _ = rate_constants(EA1_BASE, EA2_BASE, DGB_BASE, DGC_BASE, T_REF, A_PREF) print("k1 = %.6e s^-1 k1r = %.6e s^-1 (K1 = k1/k1r = %.6e)" % (k1, k1r, k1 / k1r)) t_final_v1 = 20.0 / (k1 + k1r) # about 20 relaxation times: well past equilibrium n_steps_v1 = 200000 ts, A_rk4, B_rk4, _ = rk4_integrate(k1, k1r, 0.0, 0.0, A0, t_final_v1, n_steps_v1) print() print("%14s %16s %16s %16s %16s" % ("t (s)", "B numeric", "B analytic", "abs diff", "rel diff")) check_idx = np.unique(np.round(np.geomspace(1, n_steps_v1, 22)).astype(int)) max_abs, max_rel = 0.0, 0.0 for i in check_idx: t = ts[i] _, B_an = analytic_single_reaction(k1, k1r, A0, t) B_num = B_rk4[i] d = abs(B_num - B_an) rel = d / max(abs(B_an), 1e-300) max_abs = max(max_abs, d) max_rel = max(max_rel, rel) print("%14.6g %16.10f %16.10f %16.3e %16.3e" % (t, B_num, B_an, d, rel)) print() print("Max absolute difference over %d checkpoints: %.6e mol/L" % (len(check_idx), max_abs)) print("Max relative difference: %.6e" % max_rel) print("(RK4 step count %d over %.4g s, dt = %.4g s)" % (n_steps_v1, t_final_v1, t_final_v1 / n_steps_v1)) V1_PASS = max_rel < 1e-6 print("PASS (rel. error < 1e-6): %s" % V1_PASS) # ====================================================================== hdr("VALIDATION 2 of 5 -- RK4 vs the exact linear solution, full model") # ====================================================================== print("Now both channels active: A<=>B and A<=>C together, baseline") print("parameters. exact_linear is closed form (eigendecomposition); RK4") print("is the from-scratch step-by-step integrator. They solve the same") print("equations by unrelated methods, so agreement is a real check.") k1, k1r, k2, k2r = rate_constants(EA1_BASE, EA2_BASE, DGB_BASE, DGC_BASE, T_REF, A_PREF) print() print("k1 = %.6e k1r = %.6e K1 = %.6e" % (k1, k1r, k1 / k1r)) print("k2 = %.6e k2r = %.6e K2 = %.6e" % (k2, k2r, k2 / k2r)) print("K2/K1 = %.6e (thermodynamic limit) k2/k1 = %.6e (kinetic limit)" % (k2 / k2r / (k1 / k1r), k2 / k1)) t_final_v2 = 2.0e4 # s; long enough to see both channels move, short enough for a cheap fixed-step RK4 n_steps_v2 = 400000 ts2, A_rk4b, B_rk4b, C_rk4b = rk4_integrate(k1, k1r, k2, k2r, A0, t_final_v2, n_steps_v2) print() print("%14s %16s %16s %16s %16s" % ("t (s)", "C-B numeric", "C-B exact", "abs diff", "rel diff (of A0)")) idx2 = np.unique(np.round(np.geomspace(1, n_steps_v2, 18)).astype(int)) max_abs2, max_rel2 = 0.0, 0.0 for i in idx2: t = ts2[i] A_e, B_e, C_e = exact_linear(k1, k1r, k2, k2r, A0, t) num = C_rk4b[i] - B_rk4b[i] exact = C_e - B_e d = abs(num - exact) rel = d / A0 max_abs2 = max(max_abs2, d) max_rel2 = max(max_rel2, rel) print("%14.6g %16.10f %16.10f %16.3e %16.3e" % (t, num, exact, d, rel)) print() print("Max absolute difference in (C-B) over %d checkpoints: %.6e mol/L" % (len(idx2), max_abs2)) print("Max difference relative to A0: %.6e" % max_rel2) V2_PASS = max_rel2 < 1e-6 print("PASS (rel. error < 1e-6 of A0): %s" % V2_PASS) # ====================================================================== hdr("VALIDATION 3 of 5 -- short-time (kinetic) and long-time (thermodynamic) limits") # ====================================================================== kin_limit_analytic = k2 / k1 thermo_limit_analytic = (k2 / k2r) / (k1 / k1r) print("Analytic kinetic limit, lim t->0 C/B = k2/k1 = %.10e" % kin_limit_analytic) print("Analytic thermodynamic limit, lim t->inf C/B = K2/K1 = %.10e" % thermo_limit_analytic) print() print("Simulated ratio C(t)/B(t) at decreasing t, via exact_linear:") print("%14s %20s %20s" % ("t (s)", "C/B simulated", "rel. dev. from k2/k1")) for t_probe in [1.0, 1e-2, 1e-4, 1e-6]: A_e, B_e, C_e = exact_linear(k1, k1r, k2, k2r, A0, t_probe) ratio = C_e / B_e rel = abs(ratio - kin_limit_analytic) / kin_limit_analytic print("%14.6g %20.12e %20.6e" % (t_probe, ratio, rel)) print() print("Below about 1e-6 s, exact_linear itself loses precision: it computes") print("B(t) as a small difference of two O(1) numbers (the equilibrium value") print("and an exponential correction that nearly cancels it), and that") print("subtraction eats significant digits as t shrinks. A single small RK4") print("step has no such cancellation -- it just multiplies an O(1) rate by a") print("small dt -- so it is the more trustworthy probe of the true t->0 limit:") print("%14s %20s %20s" % ("dt (s)", "C/B, one RK4 step", "rel. dev. from k2/k1")) for dt_probe in [1e-6, 1e-8, 1e-10, 1e-12, 1e-14, 1e-16]: y = rk4_single_step(k1, k1r, k2, k2r, A0, dt_probe) ratio = y[2] / y[1] rel = abs(ratio - kin_limit_analytic) / kin_limit_analytic print("%14.6g %20.12e %20.6e" % (dt_probe, ratio, rel)) short_t_final = 1e-14 y_short = rk4_single_step(k1, k1r, k2, k2r, A0, short_t_final) short_ratio = y_short[2] / y_short[1] short_rel_err = abs(short_ratio - kin_limit_analytic) / kin_limit_analytic print() print("At dt = %.1e s: single RK4 step gives C/B = %.10e, analytic k2/k1 = %.10e, rel. diff = %.3e" % (short_t_final, short_ratio, kin_limit_analytic, short_rel_err)) print() print("Simulated ratio C(t)/B(t) at increasing t, approaching equilibrium:") print("%14s %20s %20s" % ("t (s)", "C/B simulated", "rel. dev. from K2/K1")) for t_probe in [1e3, 1e5, 1e7, 1e9, 1e12, 1e15, 1e20, 1e40]: A_e, B_e, C_e = exact_linear(k1, k1r, k2, k2r, A0, t_probe) ratio = C_e / B_e rel = abs(ratio - thermo_limit_analytic) / thermo_limit_analytic print("%14.6g %20.12e %20.6e" % (t_probe, ratio, rel)) long_t_final = 1e40 A_e, B_e, C_e = exact_linear(k1, k1r, k2, k2r, A0, long_t_final) long_ratio = C_e / B_e long_rel_err = abs(long_ratio - thermo_limit_analytic) / thermo_limit_analytic print() print("At t = %.1e s: simulated C/B = %.10e, analytic K2/K1 = %.10e, rel. diff = %.3e" % (long_t_final, long_ratio, thermo_limit_analytic, long_rel_err)) V3_PASS = short_rel_err < 1e-6 and long_rel_err < 1e-6 print("PASS (both limits within 1e-6): %s" % V3_PASS) # ====================================================================== hdr("VALIDATION 4 of 5 -- mass conservation") # ====================================================================== total_rk4 = A_rk4b + B_rk4b + C_rk4b max_mass_dev = np.max(np.abs(total_rk4 - A0)) / A0 print("RK4 trajectory (Validation 2 run, %d steps): max |A+B+C-A0|/A0 = %.3e" % (n_steps_v2, max_mass_dev)) t_check = np.geomspace(1e-6, 1e30, 40) A_e, B_e, C_e = exact_linear(k1, k1r, k2, k2r, A0, t_check) max_mass_dev_exact = np.max(np.abs(A_e + B_e + C_e - A0)) / A0 print("exact_linear, %d points from 1e-6 s to 1e30 s: max |A+B+C-A0|/A0 = %.3e" % (len(t_check), max_mass_dev_exact)) V4_PASS = max_mass_dev < 1e-9 and max_mass_dev_exact < 1e-9 print("PASS (mass conserved to < 1e-9 relative): %s" % V4_PASS) # ====================================================================== hdr("VALIDATION 5 of 5 -- RK4 convergence order (step-size halving)") # ====================================================================== print("Global error of RK4 at fixed t_final, against the exact A<=>B") print("closed form, as the number of steps doubles. A fourth-order method") print("should quarter... no: divide the error by 16 each time the step") print("count doubles (halving h multiplies a 4th-order error by 2^-4).") k1v, k1rv, _, _ = rate_constants(EA1_BASE, EA2_BASE, DGB_BASE, DGC_BASE, T_REF, A_PREF) t_final_conv = 5.0 / (k1v + k1rv) # a few relaxation times, still transient (not saturated) _, B_an_conv = analytic_single_reaction(k1v, k1rv, A0, t_final_conv) step_counts = [8, 16, 32, 64, 128, 256, 512, 1024, 2048] errors = [] print() print("%10s %16s %16s %14s" % ("N steps", "dt (s)", "|error|", "observed order")) prev_err = None for N in step_counts: _, A_c, B_c, _ = rk4_integrate(k1v, k1rv, 0.0, 0.0, A0, t_final_conv, N) err = abs(B_c[-1] - B_an_conv) errors.append(err) dt = t_final_conv / N if prev_err is not None and err > 0: order = np.log2(prev_err / err) print("%10d %16.6e %16.6e %14.3f" % (N, dt, err, order)) else: print("%10d %16.6e %16.6e %14s" % (N, dt, err, "--")) prev_err = err orders = [np.log2(errors[i] / errors[i + 1]) for i in range(len(errors) - 1) if errors[i + 1] > 0] mean_order = float(np.mean(orders[:-2])) if len(orders) > 2 else float(np.mean(orders)) print() print("Mean observed order (excluding the last one or two points, where") print("floating-point noise floors the error): %.3f (RK4 theory: 4.000)" % mean_order) V5_PASS = mean_order > 3.7 print("PASS (observed order within 0.3 of theoretical 4): %s" % V5_PASS) ALL_PASS = V1_PASS and V2_PASS and V3_PASS and V4_PASS and V5_PASS hdr("VALIDATION SUMMARY") print("1. RK4 vs single-reaction analytic solution ......... %s" % ("PASS" if V1_PASS else "FAIL")) print("2. RK4 vs exact linear solution, full model .......... %s" % ("PASS" if V2_PASS else "FAIL")) print("3. Kinetic and thermodynamic limits ................. %s" % ("PASS" if V3_PASS else "FAIL")) print("4. Mass conservation ................................. %s" % ("PASS" if V4_PASS else "FAIL")) print("5. RK4 convergence order ............................. %s" % ("PASS" if V5_PASS else "FAIL")) if not ALL_PASS: print("At least one validation failed. Reporting the disagreement plainly") print("rather than silently proceeding; see the numbers printed above.") else: print("All five validations pass. Proceeding to the sweeps below using") print("exact_linear as the workhorse (it is closed form and its agreement") print("with RK4 has just been demonstrated to better than one part in a") print("million), and RK4 directly for the figure-1 time course.") # ====================================================================== hdr("RESULT 1 -- the baseline reaction: a time course") # ====================================================================== k1, k1r, k2, k2r = rate_constants(EA1_BASE, EA2_BASE, DGB_BASE, DGC_BASE, T_REF, A_PREF) K1, K2 = k1 / k1r, k2 / k2r print("Ea1 = %.1f kJ/mol, Ea2 = %.1f kJ/mol (thermo channel %.1f kJ/mol higher)" % (EA1_BASE, EA2_BASE, EA2_BASE - EA1_BASE)) print("DeltaG_B = %.1f kJ/mol, DeltaG_C = %.1f kJ/mol (thermo product %.1f kJ/mol more stable)" % (DGB_BASE, DGC_BASE, DGB_BASE - DGC_BASE)) print("T = %.2f K, A_pref = %.2e s^-1, A0 = %.3f mol/L" % (T_REF, A_PREF, A0)) print() print("k1 = %.6e s^-1 (half-life to first-order decay of A via this channel alone: %.4g s)" % (k1, np.log(2) / k1)) print("k1r = %.6e s^-1 K1 = k1/k1r = %.6e" % (k1r, K1)) print("k2 = %.6e s^-1 (half-life to first-order decay of A via this channel alone: %.4g s)" % (k2, np.log(2) / k2)) print("k2r = %.6e s^-1 K2 = k2/k2r = %.6e" % (k2r, K2)) print() print("Kinetic-limit ratio k2/k1 = %.6e (B forms %.1f x faster than C at t->0)" % (k2 / k1, k1 / k2)) print("Thermodynamic-limit ratio K2/K1 = %.6e (C outnumbers B by this factor at equilibrium)" % (K2 / K1)) t_star, status = crossover_time(k1, k1r, k2, k2r, A0) print() print("CROSSOVER TIME (C(t) first equals B(t)): status = %s" % status) print("t* = %s" % fmt_time(t_star)) A_s, B_s, C_s = exact_linear(k1, k1r, k2, k2r, A0, t_star) print("at t*: A = %.6f, B = %.6f, C = %.6f mol/L (B and C should match)" % (A_s, B_s, C_s)) t_peak, B_peak = find_peak(k1, k1r, k2, k2r, A0, species=1) print() print("B is not monotonic. Its own reverse reaction (k1r) is fast enough that") print("B tracks A in a fast pre-equilibrium (B/A ~ K1) while the slow channel") print("keeps draining A into C; as A falls, B is dragged down with it.") print("B peaks at t = %s, B_max = %.6f mol/L (%.2f%% of A0)" % (fmt_time(t_peak), B_peak, 100 * B_peak / A0)) print("...before the crossover at t* = %s, B has already fallen to %.6f mol/L" % (fmt_time(t_star), B_s)) print("(a %.1f%% decline from its own peak, even before C catches up to it)" % (100 * (1 - B_s / B_peak))) # a log-spaced time grid for Figure 1 and Figure 2, from well before the # fast channel moves to well past the slow channel's equilibrium t_grid = np.geomspace(1e-4, 3e8, 400) A_g, B_g, C_g = exact_linear(k1, k1r, k2, k2r, A0, t_grid) print() print("Concentrations at selected times along the curve used for Figure 1:") print("%14s %12s %12s %12s %14s" % ("t (s)", "A", "B", "C", "C/B")) for tp in [1, 10, 60, 600, 3600, 21600, 86400, 604800, 2.6e6, 3.15e7, 3.15e8]: A_p, B_p, C_p = exact_linear(k1, k1r, k2, k2r, A0, tp) print("%14.6g %12.6f %12.6f %12.6f %14.6e" % (tp, A_p, B_p, C_p, C_p / B_p)) # ====================================================================== hdr("RESULT 2 -- temperature sweep: how heat moves the crossover") # ====================================================================== T_list = np.array([250, 260, 270, 280, 290, 298.15, 310, 320, 340, 360, 380, 400], dtype=float) print("Same molecular parameters (Ea1=%.1f, Ea2=%.1f, DeltaG_B=%.1f, DeltaG_C=%.1f kJ/mol)" % (EA1_BASE, EA2_BASE, DGB_BASE, DGC_BASE)) print("swept across temperature only.") print() print("%8s %14s %14s %14s %20s" % ("T (K)", "k1 (s^-1)", "k2 (s^-1)", "K2/K1", "crossover t*")) t_star_list = [] for T in T_list: k1T, k1rT, k2T, k2rT = rate_constants(EA1_BASE, EA2_BASE, DGB_BASE, DGC_BASE, T, A_PREF) tstarT, statusT = crossover_time(k1T, k1rT, k2T, k2rT, A0) t_star_list.append(tstarT) print("%8.2f %14.6e %14.6e %14.6e %20s" % (T, k1T, k2T, (k2T / k2rT) / (k1T / k1rT), fmt_time(tstarT))) t_star_arr = np.array(t_star_list) # Effective Arrhenius activation energy of the crossover time itself: # fit ln(t*) = Ea_eff/(R T) + const invT = 1.0 / T_list lnT = np.log(t_star_arr) slope, intercept = np.polyfit(invT, lnT, 1) Ea_eff = slope * R / 1000.0 # kJ/mol print() print("Linear fit of ln(t*) against 1/T (an Arrhenius plot for the crossover") print("time itself, not for any single rate constant):") print(" slope = %.6e K intercept = %.6f" % (slope, intercept)) print(" effective activation energy of the crossover time, Ea_eff = slope*R = %.3f kJ/mol" % Ea_eff) print(" for comparison: Ea2 - Ea1 = %.1f kJ/mol (the kinetic penalty of the slow channel)" % (EA2_BASE - EA1_BASE)) fit_pred = np.exp(slope * invT + intercept) fit_resid = np.abs(fit_pred - t_star_arr) / t_star_arr print(" worst relative residual of the log-linear fit across the %d points: %.4f" % (len(T_list), np.max(fit_resid))) print() print("Room temperature (298.15 K) crossover: %s" % fmt_time(t_star_list[T_list.tolist().index(298.15)])) r_cold = t_star_list[0] / t_star_list[T_list.tolist().index(298.15)] r_hot = t_star_list[T_list.tolist().index(298.15)] / t_star_list[-1] print("Cooling from 298.15 K to %.0f K multiplies the wait by %.3g x" % (T_list[0], r_cold)) print("Heating from 298.15 K to %.0f K divides the wait by %.3g x" % (T_list[-1], r_hot)) # ====================================================================== hdr("RESULT 3 -- the regime map: barrier gap vs stability gap") # ====================================================================== print("Fixed: Ea1 = %.1f kJ/mol, DeltaG_B = %.1f kJ/mol, T = %.2f K" % (EA1_BASE, DGB_BASE, T_REF)) print("Swept: dEa = Ea2-Ea1 in [0,50] kJ/mol (26 values), dG = DeltaG_C-DeltaG_B in") print(" [-40,10] kJ/mol (26 values). Ea2 = Ea1+dEa, DeltaG_C = DeltaG_B+dG.") n_grid = 26 dEa_vals = np.linspace(0.0, 50.0, n_grid) dG_vals = np.linspace(-40.0, 10.0, n_grid) tstar_grid = np.full((n_grid, n_grid), np.nan) status_grid = np.empty((n_grid, n_grid), dtype=object) for i, dEa in enumerate(dEa_vals): for j, dG in enumerate(dG_vals): Ea2 = EA1_BASE + dEa dGC = DGB_BASE + dG k1g, k1rg, k2g, k2rg = rate_constants(EA1_BASE, Ea2, DGB_BASE, dGC, T_REF, A_PREF) tstar, status = crossover_time(k1g, k1rg, k2g, k2rg, A0) tstar_grid[i, j] = tstar status_grid[i, j] = status finite = np.isfinite(tstar_grid) n_total = tstar_grid.size n_never = int(np.sum(~finite)) AGE_UNIVERSE = 4.35e17 thresholds = [ ("within a minute", 60.0), ("within an hour", 3600.0), ("within a day", 86400.0), ("within a year", 3.15576e7), ("within a century", 3.15576e9), ("within the age of the universe", AGE_UNIVERSE), ] print() print("Out of %d grid cells (%d x %d):" % (n_total, n_grid, n_grid)) prev_frac = 0.0 for label, thresh in thresholds: frac = np.sum(finite & (tstar_grid <= thresh)) / n_total print(" %-32s %6.1f%% of cells (cumulative)" % (label, 100 * frac)) frac_never_math = n_never / n_total frac_beyond_universe = np.sum(finite & (tstar_grid > AGE_UNIVERSE)) / n_total print(" %-32s %6.1f%% of cells" % ("never (thermo product not more stable)", 100 * frac_never_math)) print(" %-32s %6.1f%% of cells" % ("crosses, but only beyond the age of the universe", 100 * frac_beyond_universe)) # A short readable slice: dG fixed at -25 kJ/mol, dEa varied j_slice = int(np.argmin(np.abs(dG_vals - (-25.0)))) print() print("Slice at dG = %.1f kJ/mol (a fairly typical stability gap), varying dEa:" % dG_vals[j_slice]) print("%10s %20s %14s" % ("dEa (kJ/mol)", "crossover t*", "status")) for i in range(0, n_grid, 3): print("%10.1f %20s %14s" % (dEa_vals[i], fmt_time(tstar_grid[i, j_slice]), status_grid[i, j_slice])) i_slice = int(np.argmin(np.abs(dEa_vals - 15.0))) print() print("Slice at dEa = %.1f kJ/mol (this article's baseline gap), varying dG:" % dEa_vals[i_slice]) print("%10s %20s %14s" % ("dG (kJ/mol)", "crossover t*", "status")) for j in range(0, n_grid, 3): print("%10.1f %20s %14s" % (dG_vals[j], fmt_time(tstar_grid[i_slice, j]), status_grid[i_slice, j])) print() print("Zooming in on dG near zero, at dEa = %.1f kJ/mol -- what 'more stable'" % dEa_vals[i_slice]) print("is actually buying you. As dG approaches zero from below (the two") print("products approach equal stability), the crossover time diverges,") print("smoothly, with no discontinuity: it is not that the thermodynamic") print("product suddenly stops winning, it just takes arbitrarily long to.") print("%12s %20s" % ("dG (kJ/mol)", "crossover t*")) Ea2_fix = EA1_BASE + dEa_vals[i_slice] for dG_probe in [-8, -4, -2, -1, -0.5, -0.2, -0.1, -0.05, -0.02, -0.01]: dGC_probe = DGB_BASE + dG_probe k1p, k1rp, k2p, k2rp = rate_constants(EA1_BASE, Ea2_fix, DGB_BASE, dGC_probe, T_REF, A_PREF) tstar_p, status_p = crossover_time(k1p, k1rp, k2p, k2rp, A0) print("%12.3f %20s" % (dG_probe, fmt_time(tstar_p) if status_p == "crossed" else status_p)) # ====================================================================== hdr("RESULT 4 -- results table (printed here, reproduced in the article)") # ====================================================================== print("Crossover time at T = 298.15 K for a grid of (dEa, dG) combinations,") print("chosen to span the interesting range at coarser resolution than the full sweep.") dEa_table = [0, 5, 10, 15, 20, 30, 40, 50] dG_table = [-40, -25, -10, -2, 10] print() header = "%10s" % "dEa\\dG" for dG in dG_table: header += "%16s" % ("%d kJ/mol" % dG) print(header) table_rows = [] for dEa in dEa_table: row_txt = "%10s" % ("%d kJ/mol" % dEa) row_vals = [] for dG in dG_table: Ea2 = EA1_BASE + dEa dGC = DGB_BASE + dG k1g, k1rg, k2g, k2rg = rate_constants(EA1_BASE, Ea2, DGB_BASE, dGC, T_REF, A_PREF) tstar, status = crossover_time(k1g, k1rg, k2g, k2rg, A0) row_vals.append((tstar, status)) if status == "never" or not np.isfinite(tstar): cell = "never" elif status == "instant": cell = "instant" else: cell = fmt_time(tstar).split(" (")[0] row_txt += "%16s" % cell table_rows.append((dEa, row_vals)) print(row_txt) # ====================================================================== hdr("Total runtime and reproducibility statement") # ====================================================================== elapsed = _time.time() - t_start print("Total wall-clock time for this script: %.2f s" % elapsed) print("No random number generator was used; there is no seed to report.") print("Every number above is deterministic given this file and this Python/numpy version.") print("numpy %s, Python %s" % (np.__version__, sys.version.split()[0])) if __name__ == "__main__": main()