<jupyter_start><jupyter_text>Chempy example notebook Credit goes to CHEMPY authors MOST OF THIS NOTEBOOK IS ADAPTED FROM THEIR EXAMPLES  First we are going to import the chempy package, see the install instructions on the chempy github, and use the "chempy._ _version_" as a sanity check for installing the correct version. Make sure to install the optional dependencies for the chempy package because other useful packages such as thermo,sympy, etc are useful for ode solving, thermodynamic constants, etc. __This is to remind myself that I am in the virtual env kittens where I installed chempy__  MOST OF THE RATE CONSTANTS VALUES ARE ARBITRARY. will need to ask elizabeth for more clarification, but this notebook is a good start to figuring out the kinetic modelling<jupyter_code>import chempy
import numpy as np
import sympy as sp
sp.init_printing()
import matplotlib.pyplot as plt
%matplotlib inline
sp.__version__, chempy.__version__<jupyter_output><empty_output><jupyter_text>We will first go through an example of balancing stoichiometry and mass fraction calculations that comes with the chempy package<jupyter_code># import classes
from chempy import balance_stoichiometry
from pprint import pprint
from chempy import mass_fractions
reactants, products = balance_stoichiometry({'HCOOH'},{'H2','CO2'})#,'CO','H2O'})
# Creates an ordered dictionary object with the reactants 
reactants
products
# Creates a dictionary object with the stoichiometric constants as the items and the reactants as the keys
pprint(dict(reactants))

pprint(dict(products))
mass_fractions
# Calculates the weight percents based on the mass 
for fractions in map(mass_fractions, [reactants, products]):
    pprint({k: '{0:.3g} wt%'.format(v*100) for k, v in fractions.items()})
    # .3g is for the number of sig figs 
reaction, product = balance_stoichiometry({'HCOOH'},{'CO','H2O'})
for fractions in map(mass_fractions, [reaction, product]):
    pprint({k: '{0:.3g} wt%'.format(v*100) for k, v in fractions.items()})
    # .4g is for the number of sig figs<jupyter_output>{'HCOOH': '100 wt%'}
{'CO': '60.9 wt%', 'H2O': '39.1 wt%'}<jupyter_text>Reaction Chemical Kinetics<jupyter_code># Import packages
from chempy import ReactionSystem 
from chempy.kinetics.ode import get_odesys
from collections import defaultdict
import numpy as np
import matplotlib.pyplot as plt
# Reaction System from string 
# The rate constants below are arbitrary
rsys = ReactionSystem.from_string("""
HCOOH -> H2 + CO2; 42
HCOOH -> H2O + CO; 17
""")  # "[H2O]" = 1.0 (actually 55.4 at RT)

# Ordinary Differential Systems of Equations 
odesys, extra = get_odesys(rsys)
# Time axis 
tout = sorted(np.concatenate((np.linspace(0, 23), np.logspace(-8, 1))))
# Initial Concentrations
c0 = defaultdict(float, {'HCOOH': 1.00, 'H2': 0.00, 'CO2': 0.000, 'H2O': 0.00, 'CO': 0.00})
# Integrated ode system 
result = odesys.integrate(tout, c0, atol=1e-12, rtol=1e-14)
# Plotting the results 
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
for ax in axes:
    result.plot(names=[k for k in rsys.substances if k != 'H2O'], ax=ax)
    ax.legend(loc='best', prop={'size': 9})
    ax.set_xlabel('Time')
    ax.set_ylabel('Concentration')
axes[1].set_ylim([1e-13, 1e-1])
axes[1].set_xscale('log')
axes[1].set_yscale('log')
fig.tight_layout()
#fig.savefig('examples/kinetics.png', dpi=72)
rsys
odesys
extra<jupyter_output><empty_output><jupyter_text>Example kinetic model with 2 sets of products<jupyter_code># EXAMPLE KINETIC MODEL
from __future__ import absolute_import, division, print_function
from collections import defaultdict
from ipywidgets import interact
import matplotlib.pyplot as plt
from chempy import Reaction, Substance, ReactionSystem
from chempy.kinetics.ode import get_odesys
from chempy.kinetics.analysis import plot_reaction_contributions
from chempy.printing.tables import UnimolecularTable, BimolecularTable
from chempy.util.graph import rsys2graph
import sympy
sympy.init_printing()
%matplotlib inline
# Substance ??
A, B, C, D, E, F= map(Substance, 'ABCDEF')
# Sympy ??
One = sympy.S.One
# Reactions object 
reactions = r0, r1 = [
    Reaction({'A'}, {'B', 'C'}, 4*One/100, name='R1: A cons.'),
    Reaction({'A'}, {'D', 'E'}, 3*One/100, name='R2: A cons.'),
]
# Reaction system 
rsys = ReactionSystem(reactions, (A, B, C, D, E))
# UnimolecularTable and BimolecularTable from Reaction system
uni, not_uni = UnimolecularTable.from_ReactionSystem(rsys)
bi, not_bi = BimolecularTable.from_ReactionSystem(rsys)
assert not (not_uni & not_bi), "There are only uni- & bi-molecular reactions in this set"
# Ordinary Differential Systems of Equations 
odesys, extra = get_odesys(rsys, include_params=True)
odesys.exprs
# Jacobian matrix of Ordinary Differential Systems of Equations 
odesys.get_jac()
reactions
rsys
uni
bi
not_uni
not_bi
# Initial Concentration
c0 = defaultdict(float, {'A': 1})
# Integrated ode system ??
result = odesys.integrate(1e10, c0, integrator='cvode', nsteps=2000)
# Dictionary object ??
{k: v for k, v in result.info.items() if not k.startswith('internal')}
odesys
extra
# Experimental rates ??
extra['rate_exprs_cb'](result.xout, result.yout)
result.plot(xscale='log', yscale='log')
# PLOTTING THE REACTION CONTRIBUTIONS DIRECTLY FROM CHEMPY ROBERTSON
fig, axes = plt.subplots(2, 2, figsize=(14, 6))
# normal y scale
plot_reaction_contributions(result, rsys, extra['rate_exprs_cb'], 'AB', axes=axes[0, :])
plot_reaction_contributions(result, rsys, extra['rate_exprs_cb'], 'AD', axes=axes[1, :])


plt.tight_layout()


fig, axes = plt.subplots(2, 2, figsize=(14, 6))
# Linear y scale
plot_reaction_contributions(result, rsys, extra['rate_exprs_cb'], 'AC', axes=axes[0, :],
                            relative=True, yscale='linear')
plot_reaction_contributions(result, rsys, extra['rate_exprs_cb'], 'AE', axes=axes[1, :],
                            relative=True, yscale='linear')

plt.tight_layout()<jupyter_output>c:\users\user1\anaconda3\envs\kittens\lib\site-packages\chempy\kinetics\analysis.py:86: RuntimeWarning: divide by zero encountered in true_divide
  factor = 1/xyp[1][:, rsys.as_substance_index(sk)] if relative else 1<jupyter_text>Now we will incorporate mass action into the kinetic model<jupyter_code>str_massaction = """
A -> B + C; 'k1'
A -> D + E; 'k2'
"""
rsys3 = ReactionSystem.from_string(str_massaction, substance_factory=lambda formula: Substance(formula))
str_massaction
rsys3
odesys3, extra3 = get_odesys(rsys3, include_params=False, lower_bounds=[0, 0, 0,0,0])
extra3['param_keys'], extra3['unique']
odesys3.exprs, odesys3.params, odesys3.names, odesys3.param_names
def integrate_and_plot(A0=1.0, B0=0.0, C0=0.0, D0=0.0, E0=0.0, lg_k1=-2, lg_k2=4, lg_tend=9):
    plt.figure(figsize=(14, 4))
    tout, yout, info = odesys3.integrate(
        10**lg_tend, {'A': A0, 'B': B0, 'C': C0,'D': D0, 'E': E0},
        {'k1': 10**lg_k1, 'k2': 10**lg_k2},
        integrator='cvode', nsteps=3000)
    plt.subplot(1, 2, 1)
    odesys3.plot_result(xscale='log', yscale='log')
    plt.legend(loc='best')
    plt.subplot(1, 2, 2)
    plt.plot(tout[tout<.05], yout[tout<.05, odesys3.names.index('B')])
    _ = plt.legend('best')
interact(integrate_and_plot) #, **kw) 
# kw for additional optional arguments 
# We could also have used SymPy to construct symbolic rates:
import sympy
rsys_sym = ReactionSystem.from_string("""
A -> B + C; sp.Symbol('k1')
A -> D + E; sp.Symbol('k2')
""", rxn_parse_kwargs=dict(globals_={'sp': sympy}), substance_factory=lambda formula: Substance(formula))
odesys_sym, _ = get_odesys(rsys_sym, params=True)
for attr in 'exprs params names param_names'.split():
    print(getattr(odesys_sym, attr))<jupyter_output>(-k1*y_0 - k2*y_0, k1*y_0, k1*y_0, k2*y_0, k2*y_0)
(k1, k2)
('A', 'B', 'C', 'D', 'E')
()<jupyter_text>For larger systems it is easy to loose track of what substances are actually playing a part, here the html tables can help (note the yellow background color):<jupyter_code>rsys.substances['F'] = F
uni, not_uni = UnimolecularTable.from_ReactionSystem(rsys)
uni
bi, not_bi = BimolecularTable.from_ReactionSystem(rsys)
bi
# # Saving a Picture 
# rsys2graph(rsys, 'pciture.png', save='.')
# from IPython.display import Image; Image('picture.png')<jupyter_output><empty_output><jupyter_text>Arrhenius Kinetic model with 2 sets of products<jupyter_code>from collections import defaultdict
from itertools import chain
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
from ipywidgets import interact
from chempy import Substance, Reaction, ReactionSystem
from chempy.kinetics.rates import Arrhenius, MassAction
from chempy.kinetics.ode import get_odesys
from chempy.printing.numbers import number_to_scientific_latex
from chempy.units import SI_base_registry, default_units as u
sp.init_printing()
%matplotlib inline
A, B, C, D, E, F= map(Substance, 'ABCDEF')
r1 = Reaction({'A'}, {'B','C'}, MassAction(Arrhenius(unique_keys=('A1', 'Ea_R_1'))))
r2 = Reaction({'A'}, {'D','E'}, MassAction(Arrhenius(unique_keys=('A2', 'Ea_R_2'))))
rsys = ReactionSystem([r1, r2])
rsys
# Temperature is at 15 C just like the calibration
# Activation energies and rates are arbitrary
params = {'A1': 1e11/u.s, 'A2': 2e11/u.s, 'Ea_R_1': 8e3*u.K, 'Ea_R_2': 8.5e3*u.K, 'temperature': 288*u.K}
c0 = defaultdict(lambda: 0*u.molar, {'A': 1*u.molar})
variables = c0.copy()
variables.update(params)
rsys.rates(variables)
odesys, extra = get_odesys(rsys, include_params=False, lower_bounds=0)
print(dict(zip(odesys.dep, odesys.names)))
print(dict(zip(odesys.params, odesys.param_names)))
odesys.exprs
params2 = params.copy()
pk = 'Ea_R_2'
params2[pk] = [(320+273), (340+273), (360+273)]*u.K
params2
params2[pk]<jupyter_output><empty_output><jupyter_text>Run time errors happen for 1 second or over for arrhenius model Will need to fix<jupyter_code>res2 = odesys.integrate(.001*u.s, c0, params2, integrator='cvode')
# plotting at the 3 temperatures 
fig, axes = plt.subplots(1, len(res2), figsize=(14, 4))
for r, ax in zip(res2, axes):
    r.plot(ax=ax)
    ax.set_title('$%s = %s$' % (pk.replace('_', '\\_'), number_to_scientific_latex(r.named_param('Ea_R_2'))))<jupyter_output><empty_output><jupyter_text>Using ipywidgets to get interactive controls:<jupyter_code>def integrate_and_plot(T_C=25):
    res = odesys.integrate(7*u.s, c0, dict(params, temperature=(T_C+273.15)*u.K), integrator='cvode')
    res.plot()
interact(integrate_and_plot)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13323079880803615, 'fraction_numerical': 0.029222339709699126, 'mean_word_length': 5.853754940711463, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 40, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '219643', 'license': 'MIT', 'n_tokens_mistral': 4030, 'n_tokens_neox': 3765, 'n_words': 1016, 'path': 'dev_2/chemical_kinetics_dev.ipynb', 'repo_name': 'raman-noodles/Raman-noodles', 'size': '10403'}