<jupyter_start><jupyter_text>_Kinetics of Reactions_ ( _Kinexns_ ) 1. Building and solving the chemical kinetic model from chemkin input files<jupyter_code>import sys, os
import matplotlib.pyplot as plt
sys.path.append("/Users/prguser/Documents/kinexns/")

from kinexns.ode_builder import *
from kinexns.ode_solver import *
from kinexns.parse_chemkin import *<jupyter_output>Could not find GLIMDA.<jupyter_text>setting up the paths for the input files. Since these input files are not in chemkin format, they can be directly used with *kinexns*.<jupyter_code>cwd = os.getcwd()

file_reactionlist, file_rateconstantlist, file_free_energy = set_paths(cwd)
file_rateconstantlist<jupyter_output><empty_output><jupyter_text>__What are these files? Lets Look at the reactionlist file first__   __Now the rate constant and free energy files__   __Now, we have all the files in the required format__Its time to create some additional variables required to build the model<jupyter_code># reac_prod_list = A list of reactants and products with their stoichiometric ratios
# speciesindices = a dictionary where the key is the species smiles and 
# value is the number to represent the species
# unique_species = a list of unique species present in the mechanism

reac_prod_list, speciesindices, unique_species = build_species_list(file_reactionlist)
first5species = {k: speciesindices[k] for k in list(speciesindices)[:5]}
first5species

# indices_to_species = dict(zip(speciesindices.values(), speciesindices.keys()))<jupyter_output><empty_output><jupyter_text>__Create the stoichiometric matrix, where the number of rows is equal to number of reactions andnumber of columns is equal to number of unique species__<jupyter_code>stoic_matrix = build_stoic_matrix(reac_prod_list, unique_species, speciesindices)
stoic_matrix<jupyter_output><empty_output><jupyter_text>__No Third body matrix required__  Now build the rate constant values, both forward and reverse<jupyter_code># be sure to use the correct value of the variable convert. This varialble mentions the unit of the activation energy 
# listed in the rate constant file. The values of convert can be 'cal', 'kcal', 'KJ' and 'J' (the unit in the
# rate constant file)
temp = 773
forward_rate_constants_all = build_forward_rates(file_rateconstantlist, temp, convert='kcal')<jupyter_output><empty_output><jupyter_text>__Since build_forward_rates fucntion returns both the forward rate constants and parameters, we need to collect the rate parameters in a seperate variable__<jupyter_code>forward_rate_constants = forward_rate_constants_all[0]<jupyter_output><empty_output><jupyter_text>__To include the reverse rate parameters and calculate them, we first need to calcualte the free energy of each species at the given temp__<jupyter_code>thermo_dict = build_free_energy_dict(file_free_energy, temp)<jupyter_output><empty_output><jupyter_text>__Next, if you want to include the reverse rate constants, calculate them using the function build_reverse_rates__<jupyter_code># The value 1 correspondce the converstion factor of the unit of energy in thermo_dict to Joule/mol.
# Since our free energies are calcualted directly from Gaussian, they are in hartrees. 
# So, we need to set the value of the factor carefully.

factor = 2625.5 * 1000

reverse_rate_constants = build_reverse_rates( thermo_dict, unique_species,stoic_matrix, factor,
                                             forward_rate_constants,temp, chemkin=False)
reverse_rate_constants

# If you do not want to include reverse rates in your calcualation just initiate them with zero
# reverse_rate_constants = np.zeros (len(forward_rate_constants))<jupyter_output>[-3.95538139e+04 -6.67824018e+04 -6.46713160e+04  6.39605931e+03
  1.48020702e+04 -3.93713416e+04 -8.59817381e+04 -5.41311150e+04
 -5.63827963e+04  1.00546936e+04 -4.30840349e+04 -6.21059137e+04
 -7.26189670e+04 -5.02363170e+04 -7.54174087e+04 -7.74649049e+04
 -3.70974748e+04 -1.50210894e+04 -5.64602748e+04 -6.25459738e+04
 -8.66387170e+04 -1.16119301e+04 -8.15002984e+04 -1.33900500e+01
 -6.45466573e+04 -7.15094833e+04 -1.43448840e+05 -6.79424264e+04
 -3.56133059e+04 -8.01613459e+04 -8.20628905e+04 -3.20507387e+04
 -5.39011999e+04  1.19623819e+04 -9.47590997e+04 -1.14485715e+04
 -7.54107137e+04 -7.30823678e+04 -9.49324614e+04 -4.18893537e+04
 -4.48966539e+04 -4.28586095e+04  2.40670266e+06  2.43393125e+06
  2.42387655e+06  2.46696059e+06  2.46696059e+06  2.42529944e+06
  2.42172375e+06  2.47818402e+06  2.45638943e+06  2.41701809e+06
  2.42869521e+06  2.46430852e+06  2.44281568e+06  2.47486642e+06
  2.45551189e+06  2.46696046e+06  2.48820296e+06  2.48416756e+06
  2.56482274e+06  2.49680[...]<jupyter_text>__Setting up the initial Condition__<jupyter_code>y0 = np.zeros((len(unique_species)), dtype = float)
y0[speciesindices['O[C@H]1[C@H](O)CO[C@@H](O)[C@@H]1O']] = 1<jupyter_output><empty_output><jupyter_text>__Now it is all set to solve the model__<jupyter_code>mod, sim = stiff_ode_solver(stoic_matrix, y0, forward_rate_constants, reverse_rate_constants, third_body = None, sim_time = 10)<jupyter_output>Final Run Statistics: --- 

 Number of steps                                 : 236
 Number of function evaluations                  : 264
 Number of Jacobian evaluations                  : 5
 Number of function eval. due to Jacobian eval.  : 230
 Number of error test failures                   : 2
 Number of nonlinear iterations                  : 260
 Number of nonlinear convergence failures        : 0

Solver options:

 Solver                   : CVode
 Linear multistep method  : BDF
 Nonlinear solver         : Newton
 Linear solver type       : DENSE
 Maximal order            : 5
 Tolerances (absolute)    : 1e-10
 Tolerances (relative)    : 1e-06

Simulation interval    : 0.0 - 10.0 seconds.
Elapsed simulation time: 0.05003202799991868 seconds.<jupyter_text>__You can now calculate the mol/mass fraction of each species with time using the following function__<jupyter_code>from rdkit import Chem
from rdkit.Chem import Descriptors

def generate_mol_mass_fraction(array, species_dict, fraction):
    
    mol_wt = np.zeros((len(species_dict)), dtype = float)

    for i, mol in enumerate(species_dict):
        mol_wt[i] = Descriptors.MolWt(Chem.MolFromSmiles(mol))
    
    if fraction == 'mass':
        mass = array * mol_wt
        frac = mass / mass[1].sum()
    if fraction == 'mol':
        frac = array / array.sum(axis = 1, dtype = float)[:,None]
    
    return frac

mole_frac = generate_mol_mass_fraction(sim, speciesindices, 'mass')

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

plt.plot(mod, mole_frac[:,speciesindices['O[C@H]1[C@H](O)CO[C@@H](O)[C@@H]1O']] * 100, color="b" , label = 'Xylose')
plt.plot(mod, mole_frac[:,speciesindices['O=CCO']] * 100, color="r" , label = 'Glycolaldehyde')
# plt.plot(mod, mole_frac[:,speciesindices['O=O']] * 100, color="r" , label = 'O2')
# plt.plot(mod, mole_frac[:,speciesindices['O']] *100, color="g" , label = 'O2')

#plt.ylim(0,100)
plt.xlim(0,2)
plt.legend()
# plt.xlabel('Time(s)')

plt.xlabel('Time(s)')
plt.ylabel('Mass Fraction (%)')
#plt.yscale('log')<jupyter_output><empty_output><jupyter_text>Sensitivity Analysis - Here, we used SALib to perform SA- The Saltelli sampler generates N*(2D+2) samples, where in is the number of samples, and D the number of model inputs. The keyword argument __calc_second_order=False__ will exclude second-order indices, resulting in a smaller sample matrix with 𝑁 * (𝐷 + 2) rows instead.- We need to supply a list of parameters and their upper and lower bounds to vary in an input file. In the example input file, all the rate parameters are varied one order of magnitude above and below their prescribed value, but the range was written as 1 and -1.- we also need to supply the name of the output files where generated parameter combinations are written. __Setting the path__<jupyter_code>#generate parameters for SA

setpath = '/Users/chowdhury/Documents/kmpy_results/SA_data_parallel/'<jupyter_output><empty_output><jupyter_text>__Importing the module and generating parameter sets__<jupyter_code>from kinexns.sensitivity_analysis import *
gen_params(1000, setpath,'params.txt','param_set.txt')<jupyter_output><empty_output><jupyter_text>__Defining the file names__<jupyter_code>import io
cwd = os.getcwd()
file_path_read = setpath + 'param_set.txt'
file_path_write = setpath +'model_solutions.txt'
file_input = setpath + 'params.txt'<jupyter_output><empty_output><jupyter_text>__Importing the parallel solver__<jupyter_code>from kinexns.parallel_solver import *<jupyter_output><empty_output><jupyter_text>__Solve the the various parameter sets using multiprocessing__<jupyter_code>#multisolve(4, file_path_read, file_path_write, forward_rate_constants, file_rateconstantlist, file_free_energy, stoic_matrix, unique_species, 'kcal', y0, 10, chemkin=False)<jupyter_output><empty_output><jupyter_text>__Analyze sensitivity on a specific model output__<jupyter_code>analyze_sensitivity(setpath, speciesindices['O=CCO'], "\'\t'", 2,'glycolaldehyde', parallel=True)<jupyter_output><empty_output><jupyter_text>__Plot SA results for top n sensitive parameters__<jupyter_code>plot_sensitivity_results(setpath, 'glycolaldehyde', number=5)
draw_top_reactions(setpath, 'glycolaldehyde', reac_prod_list)
# import aggdraw
from rdkit import Chem
from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole

mol = 'C1=CC=CC=C[C+]1'
mol = Chem.MolFromSmiles(mol)
Draw.MolToImage(mol)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0942743288590604, 'fraction_numerical': 0.08273909395973154, 'mean_word_length': 5.324270557029178, 'pattern_counts': {'":': 0, '<': 70, '<?xml version=': 0, '>': 70, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '146626', 'license': 'MIT', 'n_tokens_mistral': 3768, 'n_tokens_neox': 3262, 'n_words': 858, 'path': 'Notebooks/example_xylose/build_solve_model.ipynb', 'repo_name': 'cmashraf/kmpy', 'size': '9536'}