The Reaction That Gives You the Wrong Product Because You Were Impatient Science Journaling Club -- kinetics-vs-equilibrium.py numpy 2.4.2, python 3.12.3 No random number generator is used anywhere in this script. ============================================================================== VALIDATION 1 of 5 -- RK4 vs the textbook A<=>B closed form ============================================================================== Single reversible reaction only (k2 = k2r = 0), baseline forward/ reverse rate constants at T = 298.15 K, A0 = 1.000 mol/L. k1 = 1.708764e-03 s^-1 k1r = 4.025278e-06 s^-1 (K1 = k1/k1r = 4.245084e+02) t (s) B numeric B analytic abs diff rel diff 0.0583843 0.0000997600 0.0000997600 1.355e-18 1.359e-14 0.116769 0.0001995100 0.0001995100 1.130e-17 5.665e-14 0.175153 0.0002992501 0.0002992501 2.114e-17 7.065e-14 0.350306 0.0005984104 0.0005984104 3.133e-17 5.236e-14 0.583843 0.0009971512 0.0009971512 3.014e-17 3.023e-14 1.05092 0.0017941545 0.0017941545 4.987e-17 2.780e-14 1.92668 0.0032868183 0.0032868183 1.214e-17 3.694e-15 3.38629 0.0057696212 0.0057696212 3.296e-17 5.713e-15 6.13035 0.0104205202 0.0104205202 5.031e-17 4.828e-15 10.9179 0.0184827007 0.0184827007 4.163e-17 2.253e-15 19.5004 0.0327711805 0.0327711805 6.939e-17 2.117e-15 34.9138 0.0579106765 0.0579106765 1.596e-16 2.756e-15 62.4128 0.1011462033 0.1011462033 3.747e-16 3.705e-15 111.631 0.1736236226 0.1736236226 3.331e-16 1.918e-15 199.674 0.2889710760 0.2889710760 2.087e-14 7.223e-14 357.078 0.4564425710 0.4564425710 4.663e-15 1.022e-14 638.549 0.6634623132 0.6634623132 3.542e-14 5.338e-14 1141.82 0.8565169161 0.8565169161 4.841e-14 5.651e-14 2041.93 0.9674450240 0.9674450240 4.663e-14 4.820e-14 3651.47 0.9957320248 0.9957320248 5.551e-15 5.575e-15 6529.76 0.9976360095 0.9976360095 1.921e-14 1.925e-14 11676.9 0.9976498682 0.9976498682 4.097e-14 4.106e-14 Max absolute difference over 22 checkpoints: 4.840572e-14 mol/L Max relative difference: 7.222935e-14 (RK4 step count 200000 over 1.168e+04 s, dt = 0.05838 s) PASS (rel. error < 1e-6): True ============================================================================== VALIDATION 2 of 5 -- RK4 vs the exact linear solution, full model ============================================================================== Now both channels active: A<=>B and A<=>C together, baseline parameters. exact_linear is closed form (eigendecomposition); RK4 is the from-scratch step-by-step integrator. They solve the same equations by unrelated methods, so agreement is a real check. k1 = 1.708764e-03 k1r = 4.025278e-06 K1 = 4.245084e+02 k2 = 3.025227e-05 k2r = 2.972085e-12 K2 = 1.017881e+07 K2/K1 = 2.397786e+04 (thermodynamic limit) k2/k1 = 1.770418e-02 (kinetic limit) t (s) C-B numeric C-B exact abs diff rel diff (of A0) 0.05 -0.0000839219 -0.0000839219 1.084e-19 1.084e-19 0.1 -0.0001678366 -0.0001678366 1.432e-16 1.432e-16 0.25 -0.0004195366 -0.0004195366 9.639e-17 9.639e-17 0.5 -0.0008388904 -0.0008388904 3.191e-16 3.191e-16 1.05 -0.0017608258 -0.0017608258 7.199e-17 7.199e-17 2.2 -0.0036856550 -0.0036856550 8.934e-17 8.934e-17 4.75 -0.0079400161 -0.0079400161 1.214e-17 1.214e-17 10.15 -0.0168870687 -0.0168870687 1.943e-16 1.943e-16 21.65 -0.0356626273 -0.0356626273 5.690e-16 5.690e-16 46.2 -0.0745068971 -0.0745068971 5.274e-16 5.274e-16 98.7 -0.1522010864 -0.1522010864 4.802e-15 4.802e-15 210.8 -0.2961040451 -0.2961040451 1.593e-14 1.593e-14 450.15 -0.5235708971 -0.5235708971 5.307e-14 5.307e-14 961.4 -0.7826933367 -0.7826933367 1.061e-13 1.061e-13 2053.2 -0.9359316288 -0.9359316288 5.895e-14 5.895e-14 4384.95 -0.9620328239 -0.9620328239 1.221e-15 1.221e-15 9364.75 -0.9618118482 -0.9618118482 0.000e+00 0.000e+00 20000 -0.9603547681 -0.9603547681 3.342e-14 3.342e-14 Max absolute difference in (C-B) over 18 checkpoints: 1.061373e-13 mol/L Max difference relative to A0: 1.061373e-13 PASS (rel. error < 1e-6 of A0): True ============================================================================== VALIDATION 3 of 5 -- short-time (kinetic) and long-time (thermodynamic) limits ============================================================================== Analytic kinetic limit, lim t->0 C/B = k2/k1 = 1.7704179195e-02 Analytic thermodynamic limit, lim t->inf C/B = K2/K1 = 2.3977864865e+04 Simulated ratio C(t)/B(t) at decreasing t, via exact_linear: t (s) C/B simulated rel. dev. from k2/k1 1 1.770421483690e-02 2.013215e-06 0.01 1.770417954494e-02 1.978990e-08 0.0001 1.770417846749e-02 4.106881e-08 1e-06 1.770415369051e-02 1.440567e-06 Below about 1e-6 s, exact_linear itself loses precision: it computes B(t) as a small difference of two O(1) numbers (the equilibrium value and an exponential correction that nearly cancels it), and that subtraction eats significant digits as t shrinks. A single small RK4 step has no such cancellation -- it just multiplies an O(1) rate by a small dt -- so it is the more trustworthy probe of the true t->0 limit: dt (s) C/B, one RK4 step rel. dev. from k2/k1 1e-06 1.770417919461e-02 2.012588e-12 1e-08 1.770417919458e-02 1.998870e-14 1e-10 1.770417919458e-02 1.959677e-16 1e-12 1.770417919458e-02 0.000000e+00 1e-14 1.770417919458e-02 0.000000e+00 1e-16 1.770417919458e-02 0.000000e+00 At dt = 1.0e-14 s: single RK4 step gives C/B = 1.7704179195e-02, analytic k2/k1 = 1.7704179195e-02, rel. diff = 0.000e+00 Simulated ratio C(t)/B(t) at increasing t, approaching equilibrium: t (s) C/B simulated rel. dev. from K2/K1 1000 1.774967422067e-02 9.999993e-01 100000 2.481437632895e-02 9.999990e-01 1e+07 1.048881895314e+00 9.999563e-01 1e+09 2.397786486503e+04 7.889564e-15 1e+12 2.397786486503e+04 7.889564e-15 1e+15 2.397786486503e+04 7.889564e-15 1e+20 2.397786486503e+04 7.889564e-15 1e+40 2.397786486503e+04 7.889564e-15 At t = 1.0e+40 s: simulated C/B = 2.3977864865e+04, analytic K2/K1 = 2.3977864865e+04, rel. diff = 7.890e-15 PASS (both limits within 1e-6): True ============================================================================== VALIDATION 4 of 5 -- mass conservation ============================================================================== RK4 trajectory (Validation 2 run, 400000 steps): max |A+B+C-A0|/A0 = 2.010e-14 exact_linear, 40 points from 1e-6 s to 1e30 s: max |A+B+C-A0|/A0 = 1.110e-16 PASS (mass conserved to < 1e-9 relative): True ============================================================================== VALIDATION 5 of 5 -- RK4 convergence order (step-size halving) ============================================================================== Global error of RK4 at fixed t_final, against the exact A<=>B closed form, as the number of steps doubles. A fourth-order method should quarter... no: divide the error by 16 each time the step count doubles (halving h multiplies a 4th-order error by 2^-4). N steps dt (s) |error| observed order 8 3.649018e+02 7.255668e-05 -- 16 1.824509e+02 3.469850e-06 4.386 32 9.122545e+01 1.902101e-07 4.189 64 4.561272e+01 1.113667e-08 4.094 128 2.280636e+01 6.737181e-10 4.047 256 1.140318e+01 4.142708e-11 4.023 512 5.701591e+00 2.568279e-12 4.012 1024 2.850795e+00 1.604272e-13 4.001 2048 1.425398e+00 1.021405e-14 3.973 Mean observed order (excluding the last one or two points, where floating-point noise floors the error): 4.125 (RK4 theory: 4.000) PASS (observed order within 0.3 of theoretical 4): True ============================================================================== VALIDATION SUMMARY ============================================================================== 1. RK4 vs single-reaction analytic solution ......... PASS 2. RK4 vs exact linear solution, full model .......... PASS 3. Kinetic and thermodynamic limits ................. PASS 4. Mass conservation ................................. PASS 5. RK4 convergence order ............................. PASS All five validations pass. Proceeding to the sweeps below using exact_linear as the workhorse (it is closed form and its agreement with RK4 has just been demonstrated to better than one part in a million), and RK4 directly for the figure-1 time course. ============================================================================== RESULT 1 -- the baseline reaction: a time course ============================================================================== Ea1 = 90.0 kJ/mol, Ea2 = 100.0 kJ/mol (thermo channel 10.0 kJ/mol higher) DeltaG_B = -15.0 kJ/mol, DeltaG_C = -40.0 kJ/mol (thermo product 25.0 kJ/mol more stable) T = 298.15 K, A_pref = 1.00e+13 s^-1, A0 = 1.000 mol/L k1 = 1.708764e-03 s^-1 (half-life to first-order decay of A via this channel alone: 405.6 s) k1r = 4.025278e-06 s^-1 K1 = k1/k1r = 4.245084e+02 k2 = 3.025227e-05 s^-1 (half-life to first-order decay of A via this channel alone: 2.291e+04 s) k2r = 2.972085e-12 s^-1 K2 = k2/k2r = 1.017881e+07 Kinetic-limit ratio k2/k1 = 1.770418e-02 (B forms 56.5 x faster than C at t->0) Thermodynamic-limit ratio K2/K1 = 2.397786e+04 (C outnumbers B by this factor at equilibrium) CROSSOVER TIME (C(t) first equals B(t)): status = crossed t* = 9.655e+06 s (111.7 day) at t*: A = 0.001156, B = 0.499422, C = 0.499422 mol/L (B and C should match) B is not monotonic. Its own reverse reaction (k1r) is fast enough that B tracks A in a fast pre-equilibrium (B/A ~ K1) while the slow channel keeps draining A into C; as A falls, B is dragged down with it. B peaks at t = 5809 s (1.614 hr), B_max = 0.979976 mol/L (98.00% of A0) ...before the crossover at t* = 9.655e+06 s (111.7 day), B has already fallen to 0.499422 mol/L (a 49.0% decline from its own peak, even before C catches up to it) Concentrations at selected times along the curve used for Figure 1: t (s) A B C C/B 1 0.998262 0.001707 0.000030 1.770421e-02 10 0.982761 0.016940 0.000300 1.770454e-02 60 0.900930 0.097347 0.001724 1.770635e-02 600 0.352888 0.635839 0.011273 1.772922e-02 3600 0.004148 0.978320 0.017532 1.792035e-02 21600 0.002266 0.978935 0.018799 1.920370e-02 86400 0.002256 0.974513 0.023231 2.383889e-02 604800 0.002176 0.939849 0.057975 6.168563e-02 2.6e+06 0.001892 0.817560 0.180547 2.208368e-01 3.15e+07 0.000251 0.108577 0.891172 8.207761e+00 3.15e+08 0.000000 0.000042 0.999958 2.397771e+04 ============================================================================== RESULT 2 -- temperature sweep: how heat moves the crossover ============================================================================== Same molecular parameters (Ea1=90.0, Ea2=100.0, DeltaG_B=-15.0, DeltaG_C=-40.0 kJ/mol) swept across temperature only. T (K) k1 (s^-1) k2 (s^-1) K2/K1 crossover t* 250.00 1.569990e-06 1.278062e-08 1.672484e+05 7.359e+10 s (2332 yr) 260.00 8.301021e-06 8.131040e-08 1.053085e+05 8.757e+09 s (277.5 yr) 270.00 3.879657e-05 4.510397e-07 6.861930e+04 1.22e+09 s (38.66 yr) 280.00 1.624135e-04 2.213786e-06 4.610161e+04 1.956e+08 s (6.199 yr) 290.00 6.159794e-04 9.736589e-06 3.183452e+04 3.559e+07 s (1.128 yr) 298.15 1.708764e-03 3.025227e-05 2.397786e+04 9.655e+06 s (111.7 day) 310.00 6.845428e-03 1.413986e-04 1.630755e+04 1.637e+06 s (18.95 day) 320.00 2.038423e-02 4.753283e-04 1.204347e+04 4.056e+05 s (4.694 day) 340.00 1.490917e-01 4.336812e-03 6.929587e+03 3.183e+04 s (8.842 hr) 360.00 8.741681e-01 3.095003e-02 4.239696e+03 3313 s (55.22 min) 380.00 4.254803e+00 1.796017e-01 2.731633e+03 437.3 s (7.289 min) 400.00 1.767810e+01 8.741681e-01 1.839083e+03 70.66 s (1.178 min) Linear fit of ln(t*) against 1/T (an Arrhenius plot for the crossover time itself, not for any single rate constant): slope = 1.384228e+04 K intercept = -30.345414 effective activation energy of the crossover time, Ea_eff = slope*R = 115.091 kJ/mol for comparison: Ea2 - Ea1 = 10.0 kJ/mol (the kinetic penalty of the slow channel) worst relative residual of the log-linear fit across the 12 points: 0.0024 Room temperature (298.15 K) crossover: 9.655e+06 s (111.7 day) Cooling from 298.15 K to 250 K multiplies the wait by 7.62e+03 x Heating from 298.15 K to 400 K divides the wait by 1.37e+05 x ============================================================================== RESULT 3 -- the regime map: barrier gap vs stability gap ============================================================================== Fixed: Ea1 = 90.0 kJ/mol, DeltaG_B = -15.0 kJ/mol, T = 298.15 K Swept: dEa = Ea2-Ea1 in [0,50] kJ/mol (26 values), dG = DeltaG_C-DeltaG_B in [-40,10] kJ/mol (26 values). Ea2 = Ea1+dEa, DeltaG_C = DeltaG_B+dG. Out of 676 grid cells (26 x 26): within a minute 3.0% of cells (cumulative) within an hour 3.0% of cells (cumulative) within a day 3.0% of cells (cumulative) within a year 20.7% of cells (cumulative) within a century 38.3% of cells (cumulative) within the age of the universe 77.5% of cells (cumulative) never (thermo product not more stable) 22.5% of cells crosses, but only beyond the age of the universe 0.0% of cells Slice at dG = -26.0 kJ/mol (a fairly typical stability gap), varying dEa: dEa (kJ/mol) crossover t* status 0.0 2.371e-06 s crossed 6.0 1.852e+06 s (21.43 day) crossed 12.0 2.173e+07 s (251.5 day) crossed 18.0 2.453e+08 s (7.772 yr) crossed 24.0 2.76e+09 s (87.46 yr) crossed 30.0 3.105e+10 s (984 yr) crossed 36.0 3.493e+11 s (1.107e+04 yr) crossed 42.0 3.93e+12 s (1.245e+05 yr) crossed 48.0 4.421e+13 s (1.401e+06 yr) crossed Slice at dEa = 14.0 kJ/mol (this article's baseline gap), varying dG: dG (kJ/mol) crossover t* status -40.0 4.879e+07 s (1.546 yr) crossed -34.0 4.879e+07 s (1.546 yr) crossed -28.0 4.879e+07 s (1.546 yr) crossed -22.0 4.879e+07 s (1.546 yr) crossed -16.0 4.882e+07 s (1.547 yr) crossed -10.0 4.918e+07 s (1.558 yr) crossed -4.0 5.375e+07 s (1.703 yr) crossed 2.0 never (thermodynamic product is not actually more stable) never 8.0 never (thermodynamic product is not actually more stable) never Zooming in on dG near zero, at dEa = 14.0 kJ/mol -- what 'more stable' is actually buying you. As dG approaches zero from below (the two products approach equal stability), the crossover time diverges, smoothly, with no discontinuity: it is not that the thermodynamic product suddenly stops winning, it just takes arbitrarily long to. dG (kJ/mol) crossover t* -8.000 4.967e+07 s (1.574 yr) -4.000 5.375e+07 s (1.703 yr) -2.000 6.258e+07 s (1.983 yr) -1.000 7.593e+07 s (2.406 yr) -0.500 9.293e+07 s (2.945 yr) -0.200 1.194e+08 s (3.783 yr) -0.100 1.414e+08 s (4.48 yr) -0.050 1.644e+08 s (5.209 yr) -0.020 1.958e+08 s (6.203 yr) -0.010 2.199e+08 s (6.967 yr) ============================================================================== RESULT 4 -- results table (printed here, reproduced in the article) ============================================================================== Crossover time at T = 298.15 K for a grid of (dEa, dG) combinations, chosen to span the interesting range at coarser resolution than the full sweep. dEa\dG -40 kJ/mol -25 kJ/mol -10 kJ/mol -2 kJ/mol 10 kJ/mol 0 kJ/mol 0.0002698 s 0.0002028 s 0.0001755 s instant never 5 kJ/mol 1.203e+06 s 1.203e+06 s 1.214e+06 s 1.607e+06 s never 10 kJ/mol 9.655e+06 s 9.655e+06 s 9.733e+06 s 1.244e+07 s never 15 kJ/mol 7.307e+07 s 7.307e+07 s 7.365e+07 s 9.368e+07 s never 20 kJ/mol 5.497e+08 s 5.497e+08 s 5.54e+08 s 7.043e+08 s never 30 kJ/mol 3.105e+10 s 3.105e+10 s 3.13e+10 s 3.978e+10 s never 40 kJ/mol 1.754e+12 s 1.754e+12 s 1.768e+12 s 2.247e+12 s never 50 kJ/mol 9.907e+13 s 9.907e+13 s 9.985e+13 s 1.269e+14 s never ============================================================================== Total runtime and reproducibility statement ============================================================================== Total wall-clock time for this script: 38.62 s No random number generator was used; there is no seed to report. Every number above is deterministic given this file and this Python/numpy version. numpy 2.4.2, Python 3.12.3