MentDB WeakMQL Recipes
← innov-ai.com All profiles Portal
Profiles / Chemistry & life
AP

Advanced physical chemist

Kinetics, stiff ODEs, activation energies and equilibria from thermodynamics.

10 recipes on this page · open in the interactive view

Recipes
  1. Activation energy from rate constants
  2. Stiff kinetics and the steady-state approximation
  3. Equilibrium yield from the Gibbs energy
  4. Maxwell-Boltzmann speeds of nitrogen
  5. Huckel theory of butadiene: energies, bond orders and beta
  6. Heat capacity of diamond: Einstein vs Debye
  7. Consecutive reactions A to B to C by eigen-decomposition
  8. Activation enthalpy and entropy by the Eyring equation
  9. BET specific surface area from a nitrogen isotherm
  10. Rotational populations of HCl
executed on MentDB

Activation energy from rate constants

First order rate constants were measured between 290 and 330 K. The script fits the Arrhenius law directly on k (non-linear least squares, centred on Tref = 310 K for a well conditioned fit) and returns Ea with its standard error, the pre-exponential factor A and the half-life at 37 degC.

Pipeline
01 · fitArrhenius k = kref.exp(-Ea/R(1/T - 1/Tref))stat fit_model→02 · extractEa, its error and krefjson select→03 · reportEa rounded with its uncertaintymath round_uncertainty
Script
admin
#1. Rate constants measured from 290 to 330 K, in 1/s;
-> "[t]" "[290, 300, 310, 320, 330]";
-> "[k]" "[0.0005762, 0.001556, 0.00384, 0.009142, 0.02068]";
#2. Arrhenius fit centred on Tref = 310 K, a well conditioned form;
json load "rx_ch_ea" (stat fit_model "kref*exp(-Ea/8.314462618*(1/T - 1/310))" "T" [t] [k] "{\"kref\": 0.004, \"Ea\": 60000}");
-> "[ea]" (json select "rx_ch_ea" /parameters/Ea); -> "[dea]" (json select "rx_ch_ea" /errors/Ea); -> "[kref]" (json select "rx_ch_ea" /parameters/kref);
#3. Pre-exponential factor A and half-life at 37 degC;
-> "[a]" (math calc "[kref]*exp([ea]/(_R*310))");
concat "Ea = " (math round_uncertainty (math calc "[ea]/1000") (math calc "[dea]/1000")) " kJ/mol, A = " (math calc "[a]") " 1/s, t1/2 at 310 K = " (math round (math calc "ln(2)/[kref]") 1) " s"
Result
mentdb
Ea = 71.52 ± 0.15 kJ/mol, A = 4.331304805188175E9 1/s, t1/2 at 310 K = 180.2 s
Go further
  • Fit the Eyring equation instead to get the activation enthalpy and entropy.
  • Weight the points with the measurement errors (6th argument of stat fit_model).
executed on MentDB

Stiff kinetics and the steady-state approximation

Mechanism A <-> I -> P with k1 = 10, k-1 = 1e5 and k2 = 1e4 per s: the intermediate is very fast, so the ODE system is stiff. The script integrates it with the implicit Radau method and checks the steady-state approximation k_eff = k1.k2/(k-1 + k2).

Pipeline
01 · odeStiff ODE system (Radau IIA)num ode→02 · extractP at t = 2 sjson select→03 · approxSteady-state k_eff and P(t)math calc
Script
admin
#1. Pre-equilibrium A <-> I -> P: k1 = 10, k-1 = 1e5, k2 = 1e4 per s, a stiff system;
json load "rx_ch_ss" (num ode "{\"A\": \"-10*A + 1e5*I\", \"I\": \"10*A - 1e5*I - 1e4*I\", \"P\": \"1e4*I\"}" "{\"A\": 1, \"I\": 0, \"P\": 0}" 0 2 "{\"points\": 3, \"method\": \"stiff\"}");
-> "[p]" (json select "rx_ch_ss" /P[2]);
#2. Steady-state approximation: k_eff = k1.k2 / k-1 + k2;
-> "[keff]" (math calc "10*1e4/(1e5+1e4)");
-> "[p_ss]" (math calc "1-exp(-[keff]*2)");
concat "Radau: P(2 s) = " (math round [p] 6) ", steady state: " (math round [p_ss] 6) " with k_eff = " (math round [keff] 5) " 1/s, solver " (json select "rx_ch_ss" /method)
Result
mentdb
Radau: P(2 s) = 0.837654, steady state: 0.837679 with k_eff = 0.90909 1/s, solver radau
Go further
  • Draw the concentrations with plot ode.
  • Add an event option to num ode to find when P reaches 50 %.
executed on MentDB

Equilibrium yield from the Gibbs energy

Haber-Bosch at 700 K and 200 bar: from dH0 = -92.2 kJ/mol and dS0 = -198.7 J/mol.K the script computes dG0 = dH0 - T.dS0, the equilibrium constant Kp = exp(-dG0/RT), then solves the extent of reaction from a 1:3 N2/H2 feed (ideal gases) and gives the ammonia mole fraction.

Pipeline
01 · gibbsdG0 = dH0 - T.dS0 with unitsunit calc→02 · kKp = exp(-dG0/RT)math calc→03 · extentSolve the equilibrium for the extentnum solve
Script
admin
#1. N2 + 3 H2 = 2 NH3 at 700 K: dH0 = -92.2 kJ/mol, dS0 = -198.7 J/mol.K;
-> "[dg]" (unit calc "-92.2 kJ/mol - 700 K * (-198.7 J/(mol*K))" "kJ/mol");
-> "[kp]" (math calc "exp(-[dg]*1000/(_R*700))");
#2. Extent x from 1 mol N2 + 3 mol H2 at 200 bar, P0 = 1 bar, total moles 4 - 2x;
-> "[x]" (num solve (concat "(2*x)^2*(4-2*x)^2 - " [kp] "*200^2*(1-x)*(3-3*x)^3") "x" 0 0.999);
#3. Ammonia mole fraction at equilibrium;
concat "dG0 = " (math round [dg] 2) " kJ/mol, Kp = " (math calc "[kp]") ", extent = " (math round [x] 4) ", y_NH3 = " (math round (math calc "100*2*[x]/(4-2*[x])") 1) " %"
Result
mentdb
dG0 = 46.89 kJ/mol, Kp = 3.1702509874736713E-4, extent = 0.5784, y_NH3 = 40.7 %
Go further
executed on MentDB

Maxwell-Boltzmann speeds of nitrogen

For N2 at 300 K, the script gets the mean speed from an exact integral of the CAS, the most probable speed by maximisation, and the fraction of molecules faster than 1000 m/s with the error function.

Pipeline
01 · aa = m/2kTmath calc→02 · meanExact moment integralcas integrate→03 · modeMost probable speednum maximize→04 · tailFraction above 1000 m/serf(x)
Script
admin
#1. Maxwell-Boltzmann speeds of N2 at 300 K: a = m / 2kT;
-> "[m]" (math calc "28.014e-3/_NA");
-> "[a]" (math calc "[m]/(2*_k*300)");
#2. Mean speed with the CAS: with u = v sqrt a, mean = integral of 4/sqrt pi u^3 exp -u^2 / sqrt a;
-> "[I]" (cas integrate "4/sqrt(pi)*u^3*exp(-u^2)" "u" 0 "Infinity");
#3. Most probable speed: maximum of the distribution;
json load "rx_mb" (num maximize (concat "v^2*exp(-" [a] "*v^2)") "v" 1 2000);
#4. Fraction faster than 1000 m/s: 1 - erf of x + 2x/sqrt pi * exp of -x^2, x = v sqrt a;
-> "[x]" (math calc "1000*sqrt([a])");
-> "[frac]" (math calc "1-erf([x])+2*[x]/sqrt(pi)*exp(-[x]^2)");
concat "{\"integral\": \"" [I] "\", \"mean_speed_m_s\": " (math round (math calc (concat "(" [I] ")/sqrt([a])")) 1) ", \"most_probable_m_s\": " (math round (json select "rx_mb" /x) 1) ", \"fraction_above_1000_m_s\": " (math round [frac] 4) "}"
Result
mentdb
{"integral": "2/sqrt(pi)", "mean_speed_m_s": 476.2, "most_probable_m_s": 422, "fraction_above_1000_m_s": 0.0105}
Go further
  • Compare with the rms speed sqrt(3kT/m) given by math calc.
  • Change the gas: hydrogen at 300 K has 2.5 % of molecules faster than 3000 m/s.
executed on MentDB

Huckel theory of butadiene: energies, bond orders and beta

Diagonalize the Huckel matrix of butadiene: orbital energies alpha + x.beta, pi energy and delocalization energy, Coulson bond orders of the C1-C2 and C2-C3 bonds from the MO coefficients, and a spectroscopic beta from the 217 nm absorption.

Pipeline
01 · MOsEigenvalues and eigenvectorsmatrix eigen→02 · energyE_pi and delocalizationmath calc→03 · bond ordersp_rs = 2 sum c_r.c_sjson select→04 · betahc/lambda in eV over the gapunit convert
Script
admin
#1. Huckel matrix of butadiene in units of beta, E = alpha + x.beta;
-> "[H]" "[[0, 1, 0, 0], [1, 0, 1, 0], [0, 1, 0, 1], [0, 0, 1, 0]]";
json load "rx_tu1_hk" (matrix eigen [H]);
-> "[x1]" (json select "rx_tu1_hk" "/values/[0]"); -> "[x2]" (json select "rx_tu1_hk" "/values/[1]"); -> "[x3]" (json select "rx_tu1_hk" "/values/[2]");
#2. pi energy of 4 electrons in the two bonding MOs, delocalization energy against two ethylenes 4 alpha + 4 beta;
-> "[Epi]" (math calc "2*[x1]+2*[x2]");
#3. Coulson bond orders p12 and p23 from the MO coefficients, sum over the occupied MOs of 2.c_r.c_s;
-> "[c1]" (json select "rx_tu1_hk" "/vectors/[0]"); -> "[c2]" (json select "rx_tu1_hk" "/vectors/[1]");
json load "rx_tu1_c1" [c1]; json load "rx_tu1_c2" [c2];
-> "[p12]" (math calc (concat "2*(" (json select "rx_tu1_c1" "/[0]") "*" (json select "rx_tu1_c1" "/[1]") "+" (json select "rx_tu1_c2" "/[0]") "*" (json select "rx_tu1_c2" "/[1]") ")"));
-> "[p23]" (math calc (concat "2*(" (json select "rx_tu1_c1" "/[1]") "*" (json select "rx_tu1_c1" "/[2]") "+" (json select "rx_tu1_c2" "/[1]") "*" (json select "rx_tu1_c2" "/[2]") ")"));
#4. Spectroscopic beta from the observed lambda max of butadiene, 217 nm: HOMO-LUMO gap = x2 - x3 times beta;
-> "[Eobs]" (unit convert (math calc (concat (const get "h") "*" (const get "c") "/217e-9")) "J" "eV");
-> "[beta]" (math calc "[Eobs]/([x2]-[x3])");
concat "x = " (math round [x1] 4) ", " (math round [x2] 4) ", " (math round [x3] 4) ", " (math round (json select "rx_tu1_hk" "/values/[3]") 4) "; E_pi = 4 alpha + " (math round [Epi] 4) " beta; delocalization = " (math round (- [Epi] 4) 4) " beta = " (math round (math calc "([Epi]-4)*75") 1) " kJ/mol with beta = -75 kJ/mol; p12 = " (math round [p12] 3) ", p23 = " (math round [p23] 3) "; gap 217 nm = " (math round [Eobs] 3) " eV, spectroscopic |beta| = " (math round [beta] 3) " eV"
Result
mentdb
x = 1.618, 0.618, -0.618, -1.618; E_pi = 4 alpha + 4.4721 beta; delocalization = 0.4721 beta = 35.4 kJ/mol with beta = -75 kJ/mol; p12 = 0.894, p23 = 0.447; gap 217 nm = 5.714 eV, spectroscopic |beta| = 4.622 eV
Go further
  • Do benzene (6-ring matrix): delocalization energy 2 beta, the aromatic stabilization
  • Get the pi charges q_r = 2 sum c_r^2 and check they are all 1 for an alternant hydrocarbon
executed on MentDB

Heat capacity of diamond: Einstein vs Debye

Compare the molar heat capacity of diamond at 300 K given by the Einstein model (sinh form), by the Debye integral with the low-temperature thetaD of 2230 K, and the Dulong-Petit limit, and compute the vibrational partition function of an oscillator.

Pipeline
01 · einsteinCv = 3R.y^2/sinh^2 ymath sinh→02 · debyeDebye integral up to thetaD/Tnum integrate→03 · partitionq_vib = -1/expm1 of -theta/Tmath expm1→04 · compareCv of both models vs 3Rconcat
Script
admin
#1. Diamond at 300 K: Einstein and Debye temperatures in K;
-> "[t]" 300; -> "[te]" 1320; -> "[td]" 2230;
#2. Einstein model: Cv = 3R.y^2 / sinh^2 y with y = thetaE/2T;
-> "[y]" (math calc "[te]/(2*[t])");
-> "[cv_e]" (math calc (concat "3*_R*[y]^2/" (math sinh [y]) "^2"));
#3. Debye model: Cv = 9R.T/thetaD^3 . integral of x^4.e^x / e^x-1 ^2 from 0 to thetaD/T;
-> "[int]" (num integrate "x^4*exp(x)/(exp(x)-1)^2" "x" 0.000001 (math calc "[td]/[t]"));
-> "[cv_d]" (math calc "9*_R*([t]/[td])^3*[int]");
#4. Vibrational partition function of one Einstein oscillator, q = 1/1-e^-x = -1/expm1 of -x;
-> "[q]" (math calc (concat "-1/" (math expm1 (math calc "-[te]/[t]"))));
concat "Cv Einstein = " (math round [cv_e] 2) " J/mol/K, Cv Debye = " (math round [cv_d] 2) " J/mol/K, Dulong-Petit = " (math round (math calc "3*_R") 2) " J/mol/K, measured about 6.1; q_vib = " (math round [q] 5)
Result
mentdb
Cv Einstein = 6.08 J/mol/K, Cv Debye = 4.13 J/mol/K, Dulong-Petit = 24.94 J/mol/K, measured about 6.1; q_vib = 1.01243
Go further
executed on MentDB

Consecutive reactions A to B to C by eigen-decomposition

Solve the linear kinetic system of two consecutive first-order reactions exactly with the eigenvalues and eigenvectors of the rate matrix, at the time where the intermediate B is maximum, and check against the analytic formula.

Pipeline
01 · eigenEigenvalues and eigenvectors of Kmatrix eigen→02 · exponentialdiag of exp of lambda.tmatrix diag→03 · solutionc = V.e^Lt.V^-1.c0matrix mul→04 · checkAnalytic B max = k1/k2 ^ k2/k2-k1math power
Script
admin
#1. Consecutive first-order reactions A to B to C, k1 = 0.3 and k2 = 0.1 per min: dc/dt = K.c;
-> "[k]" "[[-0.3, 0, 0], [0.3, -0.1, 0], [0, 0.1, 0]]";
json load "rx_p4c_eig" (matrix eigen [k]);
#2. Eigenvectors as columns V, time of the maximum of B = ln k1/k2 / k1-k2;
-> "[v]" (matrix transpose (json select "rx_p4c_eig" /vectors));
-> "[t]" (math calc "ln(0.3/0.1)/(0.3-0.1)");
#3. Diagonal matrix of exp of lambda.t;
-> "[e]" (matrix diag (concat "[" (math exp (math calc (concat (json select "rx_p4c_eig" /values/[0]) "*[t]"))) ", " (math exp (math calc (concat (json select "rx_p4c_eig" /values/[1]) "*[t]"))) ", " (math exp (math calc (concat (json select "rx_p4c_eig" /values/[2]) "*[t]"))) "]"));
#4. Concentrations c of t = V.exp of Lambda.t.V^-1.c0, with c0 = 1 mol/L of A;
-> "[c]" (matrix mul (matrix mul (matrix mul [v] [e]) (matrix inverse [v])) "[[1], [0], [0]]");
concat "t max of B = " (math round [t] 3) " min; A = " (math round (matrix get [c] 1 1) 4) ", B = " (math round (matrix get [c] 2 1) 4) " (analytic B max = " (math round (math power 3 -0.5) 4) "), C = " (math round (matrix get [c] 3 1) 4) " mol/L"
Result
mentdb
t max of B = 5.493 min; A = 0.1925, B = 0.5774 (analytic B max = 0.5774), C = 0.2302 mol/L
Go further
  • Add a reversible step: K stays linear and the same code works.
  • Compare with a numerical integration by num ode.
executed on MentDB

Activation enthalpy and entropy by the Eyring equation

Fit the Eyring equation directly on rate constants measured at four temperatures to get the activation enthalpy and entropy, then the free energy of activation and the half-life at 25 degC.

Pipeline
01 · fitk = kB.T/h.exp dS/R.exp -dH/RTstat fit_model→02 · free energydG = dH - T.dSmath calc→03 · half-lifet1/2 = ln 2 / kmath log→04 · reportdH, dS, dG, Eaconcat
Script
admin
#1. Rate constants of a hydrolysis in 1/s measured from 15 to 45 degC;
-> "[t]" "[288.15, 298.15, 308.15, 318.15]";
-> "[k]" "[0.000412, 0.00129, 0.00375, 0.0102]";
#2. Eyring equation k = kB.T/h . exp of dS/R . exp of -dH/RT, dH in kJ/mol, dS in J/mol/K;
json load "rx_p4c_ey" (stat fit_model "1.380649e-23*T/6.62607015e-34*exp(dS/8.314462618)*exp(-dH*1000/(8.314462618*T))" "T" [t] [k] "{\"dH\": 80, \"dS\": -30}");
-> "[dh]" (json select "rx_p4c_ey" /parameters/dH); -> "[ds]" (json select "rx_p4c_ey" /parameters/dS);
#3. Free energy of activation and half-life at 25 degC;
-> "[dg]" (math calc "[dh]-298.15*[ds]/1000");
-> "[k25]" (math calc "_k*298.15/_h*exp(-[dg]*1000/(_R*298.15))");
-> "[t12]" (math calc (concat (math log 2) "/[k25]/60"));
concat "dH = " (math round [dh] 1) " kJ/mol, dS = " (math round [ds] 1) " J/mol/K, dG at 25 degC = " (math round [dg] 1) " kJ/mol, half-life = " (math round [t12] 1) " min, Ea = dH + RT = " (math round (math calc "[dh]+_R*298.15/1000") 1) " kJ/mol"
Result
mentdb
dH = 79 kJ/mol, dS = -35.4 J/mol/K, dG at 25 degC = 89.5 kJ/mol, half-life = 9 min, Ea = dH + RT = 81.5 kJ/mol
Go further
  • A negative dS points to an associative transition state; compare solvents.
  • Propagate the fit errors of dH and dS with num uncertainty.
executed on MentDB

BET specific surface area from a nitrogen isotherm

Fit the BET equation on a N2 adsorption isotherm at 77 K in the 0.05-0.30 P/P0 range, get the monolayer volume and C constant, the specific surface area and the relative pressure of monolayer completion.

Pipeline
01 · fitBET isotherm, Vm and Cstat fit_model→02 · surfaceVm.NA.sigma per gramunit calc→03 · monolayerP/P0 where V = Vmnum solve
Script
admin
#1. N2 isotherm at 77 K on an activated alumina: relative pressure and volume adsorbed in cm3 STP/g;
-> "[x]" "[0.05, 0.08, 0.11, 0.14, 0.17, 0.20, 0.25, 0.30]";
-> "[v]" "[52.1, 57.6, 61.9, 65.7, 69.2, 72.6, 78.3, 84.4]";
#2. BET isotherm V = Vm.C.x / 1-x . 1-x+C.x;
json load "rx_p4c_bet" (stat fit_model "Vm*C*x/((1-x)*(1-x+C*x))" "x" [x] [v] "{\"Vm\": 50, \"C\": 100}");
-> "[vm]" (json select "rx_p4c_bet" /parameters/Vm); -> "[c]" (json select "rx_p4c_bet" /parameters/C);
#3. Specific surface: Vm / 22414 . NA . 0.162 nm2 per N2 molecule;
-> "[s]" (unit calc (concat [vm] " cm3/g / (22414 cm3/mol) * _NA * 0.162 nm^2") "m^2/g");
#4. Relative pressure where the monolayer is complete, from the BET equation V = Vm;
-> "[xm]" (num solve (concat [c] "*x/((1-x)*(1-x+" [c] "*x)) - 1") "x" 0.001 0.5);
concat "Vm = " (math round [vm] 2) " cm3 STP/g, C = " (math round [c] 0) ", S_BET = " (math round [s] 0) " m2/g, monolayer at P/P0 = " (math round [xm] 3) ", R2 = " (math round (json select "rx_p4c_bet" /r2) 5)
Result
mentdb
Vm = 60.83 cm3 STP/g, C = 81, S_BET = 265 m2/g, monolayer at P/P0 = 0.1, R2 = 0.9997
Go further
  • Check the Rouquerol criteria: V.(1-x) must increase over the chosen range.
  • Use the t-plot to separate the micropore volume.
executed on MentDB

Rotational populations of HCl

From the rotational constant of HCl, compute the rotational temperature and partition function at 298 K, the most populated J level and its population, and where the strongest R-branch line appears.

Pipeline
01 · constantstheta_rot = hcB/kmath calc→02 · levelJ max rounded to an integermath rint→03 · populationBoltzmann fraction of J maxmath calc→04 · lineWavenumber of the R lineconcat
Script
admin
#1. HCl at 298 K, rotational constant B = 10.59 cm-1;
-> "[b]" 10.59; -> "[t]" 298.15;
#2. Rotational temperature and high-temperature partition function q = kT / hcB;
-> "[trot]" (math calc "_h*_c*[b]*100/_k");
-> "[q]" (math calc "[t]/[trot]");
#3. Most populated level J max = sqrt of kT/2hcB - 1/2, rounded to an integer;
-> "[jm]" (math rint (math calc "sqrt([t]/(2*[trot]))-0.5"));
#4. Fraction of the molecules in J max: 2J+1.exp of -J.J+1.theta/T / q;
-> "[p]" (math calc "(2*[jm]+1)*exp(-[jm]*([jm]+1)*[trot]/[t])/[q]");
#5. Wavenumber of the R branch line from J max: 2B.J+1;
concat "theta_rot = " (math round [trot] 2) " K, q_rot = " (math round [q] 2) ", J max = " [jm] " with " (math round (math calc "100*[p]") 1) " % of the molecules, strongest R line near " (math calc "2*[b]*([jm]+1)") " cm-1 from the band origin"
Result
mentdb
theta_rot = 15.24 K, q_rot = 19.57, J max = 3 with 19.4 % of the molecules, strongest R line near 84.72 cm-1 from the band origin
Go further
  • Tabulate all populations with num table and draw them with plot chart bars.
  • Add the centrifugal distortion constant D for high J.

Also in Chemistry & life

CH
Chemist
Stoichiometry, concentrations, dilutions, pH and gases with real units.
10 recipes
CR
Corrosion researcher
Mass loss, Faraday, Tafel / Stern-Geary, Nernst and Arrhenius lifetimes.
10 recipes
EB
Electrochemist & battery engineer
Nernst, cell energy, Peukert and state of charge.
10 recipes
MS
Materials scientist
Tensile curves, thermal expansion and hardness statistics.
10 recipes
PR
Process & reactor engineer
Batch and continuous reactors, conversion, residence time and heat of reaction.
10 recipes
PP
Pharmacist & pharmacokinetics
Compartment models, half-life, AUC and dosing regimens.
10 recipes
BM
Biologist & microbiologist
Growth curves, enzyme kinetics and dilution series.
10 recipes
EP
Epidemiologist
Epidemic models, reproduction numbers and incidence statistics: ODE models, growth fits and exact Poisson limits.
10 recipes
GB
Geneticist & bioinformatician
Sequences, allele frequencies and Mendelian risks: regex on DNA, Hardy-Weinberg tests, exact fractions.
10 recipes
MP
Medical physicist
Dose, attenuation, decay correction and image quality: MIRD dose, HVL, dual energy, MTF.
10 recipes
AG
Agronomist
Crop development, fertilization and irrigation: degree days, response curves and rainfall statistics.
10 recipes