notebook.community

Edit and run

Kinetic Mechanism Pathway Analyzer

This analyzer is designed to provide a better user experience and improved analyzing capacity (e.g., large mechanisms) compared with Chemkin mechanism analyzer. This analyzer takes reaction ROP file (e.g. Chemkin ckcsv file), builds a reaction network, helps answer

  • if the mechanism has certain pathway (from start species to end species)

  • what is the biggest flux from/towards a certain species

When using this script, please pay attention to the sections marked with [user input] which require user to specify certain input parameters.





In [ ]:



    


from rmgpy.rmg.model import CoreEdgeReactionModel
from rmgpy.chemkin import loadChemkinFile
import os

[user input] Prepare files

Please add

  • annotated chemkin files (mechanism and species dictionary) and

  • chemkin simulation ckcsv file

to a folder with name stored in mech variable in the next cell.





In [ ]:



    


mechPath = '/home/mjliu/Documents/Models/Tetralin/acetylene/run1'
chemkinPath= os.path.join(mechPath, 'chem_annotated.inp')
dictionaryPath = os.path.join(mechPath, 'species_dictionary.txt')
ckcsvPath= os.path.join(mechPath, 'CKSoln.ckcsv')

model = CoreEdgeReactionModel()
model.core.species, model.core.reactions = loadChemkinFile(chemkinPath,dictionaryPath)



    











In [ ]:



    


# generate paris for reactions that don't have flux pairs
for rxn in model.core.reactions:
    if not rxn.pairs: rxn.generatePairs()

Create reaction network





In [ ]:



    


import networkx as nx
import matplotlib.pyplot as plt
from rmgpy.tools.extractInfoFromckcsv import getConcentrationDictFromCKCSV, getROPFromCKCSV, getFluxGraphEdgesDict
from rmgpy.chemkin import getSpeciesIdentifier
from IPython.display import display
import numpy as np
from rmgpy.rmg.pdep import PDepReaction
%matplotlib inline

step1: prepare a dict for graph edges





In [ ]:



    


firstColDict, spc_mf_dict = getConcentrationDictFromCKCSV(ckcsvPath)
first_col_dict, spc_total_rop_dict, spc_rop_dict = getROPFromCKCSV(ckcsvPath)

graph_edges_dict = getFluxGraphEdgesDict(spc_rop_dict, model.core.reactions)

graph_edges_dict_simple = {}
for pair in graph_edges_dict:
    node1 = getSpeciesIdentifier(pair[0])
    node2 = getSpeciesIdentifier(pair[1])
    graph_edges_dict_simple[(node1, node2)] = graph_edges_dict[pair]

[user input] step2: pick a time to analyze pathways

Please specify

  • the time point you want to investigate

  • the chemkin name of your initial species/reactant

  • the reaction temperature in K

  • the reaction pressure in Pa if the mechanism includes pressure-dependent reactions

  • showIntegratedFlux: True if users want to use time-integrated fluxes (integrated up to the investigated time point) in mole for analysis. Flase if users want to use instantaneous fluxes in mole/cm3/s at the investigated time point for analysis





In [ ]:



    


########## User Input ####################
time_investigated = 1.0/3600 # hour
reactant_name = 'C10H12(1)'
T = 1200 # K
P = 16*101325 # Pa
showIntegratedFlux = True 
##########################################


timepoint_index = (np.abs(firstColDict['Time_(sec)']-time_investigated*3600)).argmin()
timeList = firstColDict['Time_(sec)']
Plist = firstColDict['Pressure_(atm)']
Tlist = firstColDict['Temperature_(k)']
Vlist = firstColDict['Volume_(cm3)']

print "Investigated time point is {0} secs".format(timeList[timepoint_index])
print "At this moment, the reactant's mole fraction is {}.".format(spc_mf_dict[reactant_name][timepoint_index])

# Create moleFraction_dict at the investigated time point
moleFraction_dict = {}
for i, v in spc_mf_dict.items():
    moleFraction_dict.update({i:v[timepoint_index]})

# Create species_identifier_didct for the core species    
species_identifier_dict = {}
for i, s in enumerate(model.core.species):
    species_identifier_dict.update({getSpeciesIdentifier(model.core.species[i]): s})

# find the total consumption of the reactant in mol (positive value means consumption)
reactantConsumption = 0 # 
for t in range(timepoint_index):
    reactantROP0 = spc_total_rop_dict[reactant_name][1][t] # in mole/cm3*s
    reactantROP1 = spc_total_rop_dict[reactant_name][1][t+1] # in mole/cm3*s
    #integrate the flux using a trapezoidal rule
    reactantConsumption += 0.5 * (reactantROP1 * Vlist[t+1] + reactantROP0 * Vlist[t]) * (timeList[t+1] - timeList[t])
reactantConsumption = reactantConsumption * -1.0 # in mole

if showIntegratedFlux:
    print "At this moment, the reactant's total molar consumption is {0:.3E} mole.".format(reactantConsumption)

# create DiGraph for this moment
G = nx.DiGraph()
for pair in graph_edges_dict:
    node1 = getSpeciesIdentifier(pair[0])
    node2 = getSpeciesIdentifier(pair[1])
    e_rawdata = graph_edges_dict[pair]
    total_flux = 0
    for rxn in e_rawdata:
        if showIntegratedFlux:
            for t in range(timepoint_index):
                rxnROP0 = e_rawdata[rxn][t] # in mole/cm3
                rxnROP1 = e_rawdata[rxn][t+1] # in mole/cm3
                #integrate the flux using a trapezoidal rule
                total_flux += 0.5 * (rxnROP1 * Vlist[t+1] + rxnROP0 * Vlist[t]) * (timeList[t+1] - timeList[t]) # in mole
        else:
            total_flux += e_rawdata[rxn][timepoint_index] # in mole/cm3*s
    if total_flux >= 0:
        G.add_edge(node2, node1, {"total_flux":total_flux}) # in G, positive means production of node1
    else: 
        G.add_edge(node1, node2, {"total_flux":-total_flux}) # in G, negative means consumption of node1

[user input] Functionality 1: find most dominant species at the investigated time point

Pleae specify

  • dominant_list: If you want to visualize the first top 10 dominant species, specify the list as [0, 9]. If you want to visualize from the 5th most dominant species to the 10th most dominant species, specify the list as [4, 9].




In [ ]:



    


########## User Input ####################
dominant_list = [0, 9] 
##########################################


# sort out the moleFraction_dict in the descending order
v = list(moleFraction_dict.values())
k = list(moleFraction_dict.keys())
index_list = [i[0] for i in sorted(enumerate(v), key=lambda x:x[1], reverse=True)]
print '=============================================================='
print 'Order \t Species \t Mole Frac. \t MW (g/mol)'
print '=============================================================='

for order in range(dominant_list[0], dominant_list[-1]+1):
    index_number = index_list[order]
    species_dominant = k[index_number]
    molecule_dominant = species_identifier_dict.get(species_dominant)
    mole_fraction = v[index_number]
    molecular_weight = molecule_dominant.molecule[0].getMolecularWeight()*1000
    if len(species_dominant) < 6: # if the species label is too short, put more space to allign it
        print '{0} \t {1} \t\t {2} \t {3:.0f}'.format(order, species_dominant, mole_fraction, molecular_weight)
    else:
        print '{0} \t {1} \t {2} \t {3:.0f}'.format(order, species_dominant, mole_fraction, molecular_weight)
    display(molecule_dominant)
    print '--------------------------------------------------------------'

[user input] Functionality 2: find paths between two species

Pleae specify

  • source species

  • target species





In [ ]:



    


########## User Input ####################
src_species = 'C10H12(1)'
tgt_species = 'CH3(3)'
##########################################

output1: sort out all the paths and most significant one





In [ ]:



    


paths = list(nx.all_simple_paths(G, source=src_species, target=tgt_species, cutoff=6))



    











In [ ]:



    


path_fluxes = []
if showIntegratedFlux:
    print '========================================================================================================'
    print 'Path index \t Flux (mols)    Path steps'
else:
    print '================================================'
    print 'Path index \t Flux (mol/cm3/s)    Path steps'
print ''

path_fluxes = [(path, min([G[path[step]][path[step+1]]['total_flux'] for step in range(len(path)-1)])) for path in paths if 'H(5)' not in path]

path_fluxes.sort(key=lambda x: x[1], reverse=True)

for i, (path, flux) in enumerate(path_fluxes):
    if flux < 1e-10:
        break
    path_fluxes.append(flux)
    if showIntegratedFlux:
        print 'Path #{0}: \t {1:.3E} \t {2}'.format(i, flux, path)
    else:
        print 'Path #{0}: \t {1:.3E} \t {2}'.format(i, flux, path)

[user input] output2: visualize the path of interest

Please specify

  • the path index you want to visualize




In [ ]:



    


########## User Input ####################
path_index_investigate = 1
##########################################

path = path_fluxes[path_index_investigate][0]
path_steps = len(path) - 1
print ''
print '\t Pathway Report '
print '======================================'
print 'The pathway you are intested in has {0} steps.'.format(len(path))
for step in range(path_steps):
    step_pair = (path[step], path[step+1])
    h_abs_rxns = []
    disp_rxns = []
    Pdep_rxns = []
    
    print ""
    if showIntegratedFlux:
        print "Step{0} \t\t\t\t Integrated Flux (mole)".format(step+1)
        print '***********************************************************************'
        print "{1} --> {2} \t\t {3:.3E}".\
        format(step+1, step_pair[0], step_pair[1], G[step_pair[0]][step_pair[1]]['total_flux'])
        print '***********************************************************************'
    else:
        print "Step{0} \t\t\t\t Flux (mole/cm3/s)".format(step+1)
        print '*****************************************************'
        print "{1} --> {2} \t\t {3:.3E}".\
        format(step+1, step_pair[0], step_pair[1], G[step_pair[0]][step_pair[1]]['total_flux'])
        print '*****************************************************'
    
    if step_pair not in graph_edges_dict_simple:
        step_pair = (step_pair[1], step_pair[0])
    
    print ''
    for rxn in graph_edges_dict_simple[step_pair]:
        if isinstance(rxn, PDepReaction):
            Pdep_rxns.append(rxn)
        elif rxn.family == "H_Abstraction":
            h_abs_rxns.append(rxn)
        elif rxn.family == "Disproportionation":
            disp_rxns.append(rxn)
        else:
            print "Example reaction: rxn#{0}".format(rxn.index)
            print ''
            print str(rxn)
            display(rxn)
            print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, rxn.getRateCoefficient(T))
            reverseRate = rxn.generateReverseRateCoefficient()
            print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))
            print ''
    if len(h_abs_rxns) > 0: 
        
        print "Example H_Abstraction: rxn#{0} (1/{1} H_Abs): ".format(h_abs_rxns[0].index, len(h_abs_rxns)) 
        print ''
        print str(h_abs_rxns[0])
        display(h_abs_rxns[0])
        print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, h_abs_rxns[0].getRateCoefficient(T))
        reverseRate = h_abs_rxns[0].generateReverseRateCoefficient()
        print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))
        print ''
    if len(disp_rxns) > 0: 
        
        print "Example Disproportionation: rxn#{0} (1/{1} Disp): ".format(disp_rxns[0].index, len(disp_rxns)) 
        print ''
        print str(disp_rxns[0])
        display(disp_rxns[0])
        print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, disp_rxns[0].getRateCoefficient(T))
        reverseRate = disp_rxns[0].generateReverseRateCoefficient()
        print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))
        print ''
    if len(Pdep_rxns) > 0: 
        
        print "Example Pressure-dependent Rxn: rxn#{0} (1/{1} Pdep_rxns): ".format(Pdep_rxns[0].index, len(Pdep_rxns)) 
        print ''
        print str(Pdep_rxns[0])
        display(Pdep_rxns[0])
        print 'Forward rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, Pdep_rxns[0].getRateCoefficient(T, P))
        reverseRate = Pdep_rxns[0].generateReverseRateCoefficient()
        print 'Reverse rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, reverseRate.getRateCoefficient(T, P))
        print ''

[user input] Functionality 3: find paths from a particular species

Please specify

  • source species

  • depth of pathway you want search, e.g., depth = 2 means search pathways 2 steps down from source species

  • path_top_list, e.g., [0, 3] means 1st step of the pathway takes the most significant branch and 2nd step would be the third most siginificant branch





In [ ]:



    


########## User Input ####################
source = "S(16320)"
depth = 2
path_top_list = [0, 0]
##########################################


current_node = source

print ''
print '\t Pathway Report '
print '======================================'
print 'The pathway you are intested in has {0} steps.'.format(depth)
for step in range(depth):
    print "\n"    
    nextNode_flux_list = [(next_node, G[current_node][next_node]['total_flux']) for next_node in G[current_node]]
    sorted_nextNode_flux_list = sorted(nextNode_flux_list, key=lambda tup: -tup[1])
    
    # choose the top one as next node
    tup = sorted_nextNode_flux_list[path_top_list[step]]
    next_node = tup[0]
    step_flux = tup[1]
    
    print ""
    if showIntegratedFlux:
        print "Step{0} \t\t\t\t Integrated Flux (mole)".format(step+1)
        print '************************************************************************'
        print "{0} --> {1} \t\t {2:.3E}".\
        format(current_node, next_node, step_flux)
        print '************************************************************************'
    else:
        print "Step{0} \t\t\t\t Flux (mole/cm3/s)".format(step+1)
        print '****************************************************'
        print "{0} --> {1} \t\t {2:.3E}".\
        format(current_node, next_node, step_flux)
        print '****************************************************'

    
    step_pair = (current_node, next_node)
    if step_pair not in graph_edges_dict_simple:
        step_pair = (next_node, current_node)
    
    h_abs_rxns = []
    disp_rxns = []
    Pdep_rxns = []
    for rxn in graph_edges_dict_simple[step_pair]:
        if isinstance(rxn, PDepReaction):
            Pdep_rxns.append(rxn)
        elif rxn.family == "H_Abstraction":
            h_abs_rxns.append(rxn)
        elif rxn.family == "Disproportionation":
            disp_rxns.append(rxn)
        else:
            print ''
            print "Example reaction: rxn#{0}".format(rxn.index)
            print ''
            if showIntegratedFlux:
            #integrate the flux using a trapezoidal rule
                integrated_flux = 0
                for t in range(timepoint_index):
                    rxnROP0 = graph_edges_dict_simple[step_pair][rxn][t]
                    rxnROP1 = graph_edges_dict_simple[step_pair][rxn][t+1]
                    integrated_flux += 0.5 * (rxnROP1 * Vlist[t+1] + rxnROP0 * Vlist[t]) * (timeList[t+1] - timeList[t]) # in mole
                print "{0}: Integrated Flux = {1:.3E} mole".format(str(rxn), integrated_flux)
            else:
                print "{0}: Flux = {1:.3E} mole/cm3/s".format(str(rxn), graph_edges_dict_simple[step_pair][rxn][timepoint_index])
            display(rxn)
            print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, rxn.getRateCoefficient(T))
            reverseRate = rxn.generateReverseRateCoefficient()
            print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))
            print ''
    if len(h_abs_rxns) > 0:
        print ''
        print "Example H_Abstraction: rxn#{0}(1/{1} H_Abs)".format(h_abs_rxns[0].index, len(h_abs_rxns))
        print ''
        print str(h_abs_rxns[0])
        display(h_abs_rxns[0])
        print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, h_abs_rxns[0].getRateCoefficient(T))
        reverseRate = h_abs_rxns[0].generateReverseRateCoefficient()
        print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))
        print ''
    if len(disp_rxns) > 0: 
        print ''
        print "Example Disproportionation: rxn#{0}(1/{1} Disp)".format(disp_rxns[0].index, len(disp_rxns))
        print ''
        print str(disp_rxns[0])
        display(disp_rxns[0])
        print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, disp_rxns[0].getRateCoefficient(T))
        reverseRate = disp_rxns[0].generateReverseRateCoefficient()
        print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))
        print ''
    if len(Pdep_rxns) > 0: 
        print ''
        print "Example Pressure-dependent Rxn: rxn#{0} (1/{1} Pdep_rxns): ".format(Pdep_rxns[0].index, len(Pdep_rxns)) 
        print ''
        print str(Pdep_rxns[0])
        display(Pdep_rxns[0])
        print 'Forward rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, Pdep_rxns[0].getRateCoefficient(T, P))
        reverseRate = Pdep_rxns[0].generateReverseRateCoefficient()
        print 'Reverse rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, reverseRate.getRateCoefficient(T, P))
        print ''
    
    current_node = next_node

More detail for the H_Abstration and Disproportionation in last step of the pathway





In [ ]:



    


print ''
print 'The step you see currently is between the two species {0}'.format(step_pair)
print ''
print 'Positive flux in below means positive production of 1st node of the pair.'
print ''
print 'This step have {0} H_Abstration and Disproportionation in total.'.format(len(h_abs_rxns + disp_rxns))

flux_threshold = 1e-9

print ''
print 'TO avoid too much printout, the reactions shown below have fluxes above {0}.'.format(flux_threshold)

total_flux = 0
rxn_flux_tups = []
for rxn in h_abs_rxns + disp_rxns:
    if showIntegratedFlux:
    #integrate the flux using a trapezoidal rule
        flux = 0
        for t in range(timepoint_index):
            rxnROP0 = graph_edges_dict_simple[step_pair][rxn][t]
            rxnROP1 = graph_edges_dict_simple[step_pair][rxn][t+1]
            flux += 0.5 * (rxnROP1 * Vlist[t+1] + rxnROP0 * Vlist[t]) * (timeList[t+1] - timeList[t]) # in mole
    else:
        flux = graph_edges_dict_simple[step_pair][rxn][timepoint_index]
    rxn_flux_tups.append((rxn, flux))
    

rxn_flux_tups = sorted(rxn_flux_tups, key=lambda tup: tup[1], reverse=False)
for tup in rxn_flux_tups:
    rxn = tup[0]
    flux = tup[1]
    if flux > flux_threshold:
        total_flux += flux
        print ''
        print "**********************************************************************************"
        if showIntegratedFlux:
            print "rxn#{0}: {1}: Integrated Flux = {2:.3E} mole ".format(rxn.index, str(rxn), flux) 
        else:
            print "rxn#{0}: {1}: Flux = {2:.3E} mole/cm3/s".format(rxn.index, str(rxn), flux)
        display(rxn) 
print "***************************************"
print ''
if showIntegratedFlux:
    print "TOTAL integrated flux from h_abs and disp is {0:.3E} mole.".format(total_flux)
else:
    print "TOTAL flux from h_abs and disp is {0:.3E} mole/cm3/s.".format(total_flux)

[user input] Functionality 4: find paths towards a particular species

Similar to functionality 2, please specify

Please specify

  • target species

  • depth of pathway you want search, e.g., depth = 2 means search pathways 2 steps up beyond target species

  • path_top_list, e.g., [0, 3] means 1st step of the pathway takes the most significant branch and 2nd step would be the third most siginificant branch





In [ ]:



    


########## User Input ####################
target = "CH3(3)"
depth = 5
path_top_list = [0, 0, 0, 0, 0]
##########################################


current_node = target
print ''
print '\t Pathway Report '
print '======================================'
print 'The pathway you are intested in has {0} steps.'.format(depth)
print ''
if showIntegratedFlux:
    product_flux = 0
    for t in range(timepoint_index):
        productROP0 = spc_total_rop_dict[target][1][t]
        productROP1 = spc_total_rop_dict[target][1][t+1]
        product_flux += 0.5 * (productROP1 * Vlist[t+1] + productROP0 * Vlist[t]) * (timeList[t+1] - timeList[t]) # in mole
    print "total integrated product flux for {0} is {1:.3E} mole.".format(target, product_flux)
else:
    print "total product flux for {0} is {1:.3E} mole/cm3/s.".format(target, spc_total_rop_dict[target][1][timepoint_index])
for step in range(depth):
    print "\n"
    prev_nodes = []
    for node1 in G:
        if current_node in G[node1]:
            prev_nodes.append(node1)
    prevNode_flux_list = [(prev_node, G[prev_node][current_node]['total_flux']) for prev_node in prev_nodes]
    sorted_prevNode_flux_list = sorted(prevNode_flux_list, key=lambda tup: -tup[1])
    
    # choose the top one as next node
    tup = sorted_prevNode_flux_list[path_top_list[step]]
    prev_node = tup[0]
    step_flux = tup[1]
    
    print ""
    if showIntegratedFlux:
        print "Step{0} \t\t\t\t Integrated Flux (mole)".format(step+1)
        print '************************************************************************'
        print "{0} <-- {1} \t\t {2:.3E}".\
        format(current_node, prev_node, step_flux)
        print '************************************************************************'
    else:
        print "Step{0} \t\t\t\t Flux (mole/cm3/s)".format(step+1)
        print '****************************************************'
        print "{0} <-- {1} \t\t {2:.3E}".\
        format(current_node, prev_node, step_flux)
        print '****************************************************'
    
    step_pair = (prev_node, current_node)
    if step_pair not in graph_edges_dict_simple:
        step_pair = (current_node, prev_node)
    
    h_abs_rxns = []
    disp_rxns = []
    Pdep_rxns = []
    for rxn in graph_edges_dict_simple[step_pair]:
        if isinstance(rxn, PDepReaction):
            Pdep_rxns.append(rxn)
        elif rxn.family == "H_Abstraction":
            h_abs_rxns.append(rxn)
        elif rxn.family == "Disproportionation":
            disp_rxns.append(rxn)
        else:
            print ''
            print "Example reaction: rxn#{0}".format(rxn.index)
            print ''
            if showIntegratedFlux:
            #integrate the flux using a trapezoidal rule
                integrated_flux = 0
                for t in range(timepoint_index):
                    rxnROP0 = graph_edges_dict_simple[step_pair][rxn][t]
                    rxnROP1 = graph_edges_dict_simple[step_pair][rxn][t+1]
                    integrated_flux += 0.5 * (rxnROP1 * Vlist[t+1] + rxnROP0 * Vlist[t]) * (timeList[t+1] - timeList[t]) # in mole
                print "{0}: Integrated Flux = {1:.3E} mole".format(str(rxn), integrated_flux)
            else:
                print "{0}: Flux = {1:.3E} mole/cm3/s".format(str(rxn), graph_edges_dict_simple[step_pair][rxn][timepoint_index])
            display(rxn)
            print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, rxn.getRateCoefficient(T))
            reverseRate = rxn.generateReverseRateCoefficient()
            print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))
            
    if len(h_abs_rxns) > 0:
        print ''
        print "Example H_Abstraction: rxn#{0}(1/{1} H_Abs)".format(h_abs_rxns[0].index, len(h_abs_rxns))
        print ''
        print str(h_abs_rxns[0])
        display(h_abs_rxns[0])
        print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, h_abs_rxns[0].getRateCoefficient(T))
        reverseRate = h_abs_rxns[0].generateReverseRateCoefficient()
        print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))
        print ''
    if len(disp_rxns) > 0: 
        print ''
        print "Example Disproportionation: rxn#{0}(1/{1} Disp)".format(disp_rxns[0].index, len(disp_rxns))
        print ''
        print str(disp_rxns[0])
        display(disp_rxns[0])
        print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, disp_rxns[0].getRateCoefficient(T))
        reverseRate = disp_rxns[0].generateReverseRateCoefficient()
        print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))
        print ''
    if len(Pdep_rxns) > 0: 
        print ''
        print "Example Pressure-dependent Rxn: rxn#{0} (1/{1} Pdep_rxns): ".format(Pdep_rxns[0].index, len(Pdep_rxns)) 
        print ''
        print str(Pdep_rxns[0])
        display(Pdep_rxns[0])
        print 'Forward rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, Pdep_rxns[0].getRateCoefficient(T, P))
        reverseRate = Pdep_rxns[0].generateReverseRateCoefficient()
        print 'Reverse rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, reverseRate.getRateCoefficient(T, P))
        print ''
    
    current_node = prev_node

More detail for the H_Abstration and Disproportionation in last step of the pathway





In [ ]:



    


print ''
print 'The step you see currently is between the two species {0}'.format(step_pair)
print ''
print 'Positive flux in below means positive production of 1st node of the pair.'
print ''
print 'This step have {0} H_Abstration and Disproportionation in total.'.format(len(h_abs_rxns + disp_rxns))

flux_threshold = 1e-9

print ''
print 'TO avoid too much printout, the reactions shown below have fluxes above {0}.'.format(flux_threshold)

total_flux = 0
rxn_flux_tups = []
for rxn in h_abs_rxns + disp_rxns:
    if showIntegratedFlux:
    #integrate the flux using a trapezoidal rule
        flux = 0
        for t in range(timepoint_index):
            rxnROP0 = graph_edges_dict_simple[step_pair][rxn][t]
            rxnROP1 = graph_edges_dict_simple[step_pair][rxn][t+1]
            flux += 0.5 * (rxnROP1 * Vlist[t+1] + rxnROP0 * Vlist[t]) * (timeList[t+1] - timeList[t]) # in mole
    else:
        flux = graph_edges_dict_simple[step_pair][rxn][timepoint_index]
    rxn_flux_tups.append((rxn, flux))

rxn_flux_tups = sorted(rxn_flux_tups, key=lambda tup: tup[1], reverse=False)
for tup in rxn_flux_tups:
    rxn = tup[0]
    flux = tup[1]
    if flux > flux_threshold:
        total_flux += flux
        print ''
        print "**********************************************************************************"
        if showIntegratedFlux:
            print "rxn#{0}: {1}: Integrated Flux = {2:.3E} mole ".format(rxn.index, str(rxn), flux) 
        else:
            print "rxn#{0}: {1}: Flux = {2:.3E} mole/cm3/s".format(rxn.index, str(rxn), flux)
        display(rxn) 
print "***************************************"
print ''
if showIntegratedFlux:
    print "TOTAL integrated flux from h_abs and disp is {0:.3E} mole.".format(total_flux)
else:
    print "TOTAL flux from h_abs and disp is {0:.3E} mole/cm3/s.".format(total_flux)



    











In [ ]:

Content source: mliu49/RMG-stuff

Similar notebooks:

notebook.community | gallery | about