import matplotlib.pyplot as plt
import numpy as np
from chempy import ReactionSystem
from chempy.units import to_unitless, SI_base_registry as si, default_units as u, default_constants as const
from chempy.kinetics.ode import get_odesys
%matplotlib inline
R = 8.314472
T_K = 300
dH=80e3
dS=10
rsys1 = ReactionSystem.from_string("""
NOBr -> NO + Br; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol)
""".format(dH=dH, dS=dS))
kref = 20836643994.118652*T_K*np.exp(-(dH - T_K*dS)/(R*T_K))
kref
NOBr0_M = 0.7
init_cond = dict(
NOBr=NOBr0_M*u.M,
NO=0*u.M,
Br=0*u.M
)
t = 5*u.second
params = dict(
temperature=T_K*u.K
)
def integrate_and_plot(rsys):
odesys, extra = get_odesys(rsys, unit_registry=si, constants=const)
fig, axes = plt.subplots(1, 3, figsize=(14, 6))
res = odesys.integrate(t, init_cond, params, integrator='cvode')
t_sec = to_unitless(res.xout, u.second)
NOBr_ref = NOBr0_M*np.exp(-kref*t_sec)
cmp = to_unitless(res.yout, u.M)
ref = np.empty_like(cmp)
ref[:, odesys.names.index('NOBr')] = NOBr_ref
ref[:, odesys.names.index('Br')] = NOBr0_M - NOBr_ref
ref[:, odesys.names.index('NO')] = NOBr0_M - NOBr_ref
assert np.allclose(cmp, ref)
axes[0].plot(t_sec, cmp)
axes[1].plot(t_sec, cmp - ref)
res.plot_invariant_violations(ax=axes[2])
integrate_and_plot(rsys1)
rsys2 = ReactionSystem.from_string("""
NOBr -> NO + Br; MassAction(EyringHS([{dH}*J/mol, {dS}*J/K/mol]))
""".format(dH=dH, dS=dS))
integrate_and_plot(rsys2)