Finds all instances of homolytic fission (i.e. radical recombination, in reverse) reactions in a chemkin file, and adds the reactions you'd get with roaming.
import sys
import os
import re
import numpy
import csv
import cPickle as pickle
import itertools
import logging
from collections import Counter, defaultdict, OrderedDict
import IPython
from IPython.display import display, Markdown, HTML
def mprint(s):
"A convenience to format things in Jupyter notebooks via MarkDown syntax"
display(Markdown(s))
#import ck2cti # the customized one stored alongside this notebook
from cantera import ck2cti # we'll monkey-patch it if needed
# Add $RMGpy to front of your $PYTHONPATH in case you didn't already
sys.path.insert(1,os.getenv('RMGpy',os.path.expanduser('~/Code/RMG-Py')))
from rmgpy.molecule import Molecule
import rmgpy.kinetics
from rmgpy.data.rmg import RMGDatabase
from rmgpy.data.kinetics.library import KineticsLibrary
First, read in the chemkin file to the Cantera ck2cti
parser.
model_name = 'butanol'
model_path = 'butanol-sarathy2012'
reactions_filepath = os.path.join(model_path, 'model.txt')
thermo_filepath = os.path.join(model_path, 'thermo.txt')
transport_filepath = os.path.join(model_path, 'transport.txt')
outdir = 'CanteraFiles'
os.path.exists(outdir) or os.mkdir(outdir)
for f in os.listdir(outdir):
if f.startswith(model_name):
os.unlink(os.path.join(outdir, f))
parser = ck2cti.Parser()
surfaces = parser.convertMech(inputFile=reactions_filepath,
thermoFile=thermo_filepath,
transportFile=transport_filepath,
surfaceFile=None,
phaseName='gas',
outName=os.path.join(outdir, model_name+'.original.cti'),
permissive=True)
parser.reactions[0], str(parser.reactions[0])
Now read in the species dictionary file to get RMG Species objects for each species
get_species_dict = KineticsLibrary().getSpecies
species_dict = get_species_dict(os.path.join(model_path, 'RMG-Py-kinetics-library', 'dictionary.txt'))
species_dict
speciesList = species_dict.values()
len(speciesList)
print "These species have not been identified:"
rs = [str(r).split() for r in parser.reactions]
for s in parser.speciesList:
if s.label not in species_dict:
print s.label,
print "which participates in {} reactions".format(sum([s.label in r for r in rs]))
def is_unimolecular_dissociation(cti_reaction):
"One reactant and two products"
return ( sum([nu for nu,spec in cti_reaction.reactants]) == 1
and sum([nu for nu,spec in cti_reaction.products]) == 2)
def is_bimolecular_recombination(cti_reaction):
"Two reactants and one product"
return ( sum([nu for nu,spec in cti_reaction.reactants]) == 2
and sum([nu for nu,spec in cti_reaction.products]) == 1)
def get_rmg_species(cti_species):
"get the RMG species corresponding to the cti parser species"
return species_dict[cti_species.label]
def is_radical_recombination(r):
if not is_bimolecular_recombination(r): return False
if r.reactants[0][0] == 2:
if not get_rmg_species(r.reactants[0][1]).isRadical(): return False
def is_homolytic_fission(r):
"reverse is radical recombination"
raise NotImplementedError
# Test it out
for r in parser.reactions:
print r
print r.reactants
print r.products
print r.kinetics
if is_unimolecular_dissociation(r): break
if is_radical_recombination(r): break
rxnFamilies = ['R_Recombination'] # Only looking at R_Recombination
rmgDatabase = RMGDatabase()
databasePath = os.path.abspath(os.path.join(os.getenv('RMGpy', '..'), '..', 'RMG-database', 'input'))
print(databasePath)
rmgDatabase.load(databasePath,
kineticsFamilies=rxnFamilies,
transportLibraries=[],
reactionLibraries=[],
seedMechanisms=[],
thermoLibraries=['primaryThermoLibrary', 'BurkeH2O2', 'thermo_DFT_CCSDTF12_BAC', 'CBS_QB3_1dHR' ],
solvation=False,
)
str(r)
def expand_reagent_list(cti_reagent_list):
"""
Expand a list of tuples like
[(2,'A'),(1,'B'),(1,'B')]
(which is how cantera reactions have their reactants and products)
into a straightforward list like
['A', 'A', 'B', 'B']
"""
# cti_reagent_list = [(2,'A'),(1,'B'),(1,'B')]
reagents=[]
for nu,species in cti_reagent_list:
reagents.extend([species]*nu)
return reagents
# Generate resonance isomers
print "These species have more than one resonance isomer:"
for species in speciesList:
species.generateResonanceIsomers()
if len(species.molecule)>1:
print(species.label, species.molecule)
print r
expand_reagent_list(r.products)
def make_rmg_reaction(cti_reaction):
"""
Make an RMG reaction from the cantera reaction.
Reactant and species lists are Species, not Molecules.
"""
return rmgpy.reaction.Reaction(reactants=[get_rmg_species(s) for s in expand_reagent_list(r.reactants)],
products=[get_rmg_species(s) for s in expand_reagent_list(r.products)],
reversible=True)
make_rmg_reaction(r)
RMG normally doesn't store the atom labels once it has generated the reaction (and determined what nodes in the kinetics tree it will use), so it needs modifying to make it do so. Furthermore, because we're interested in the labels on the "product" structure (of a radical recombination reaction), these were not stored at all, previously. So RMG needs to be edited.
Rather than require users of this notebook to get and compile a separate fork of RMG, we will instead monkey-patch it. Unfortunately this involves replacing a rather long and private method rmgpy.data.kinetics.family.KineticsFamily.__generateReactions
. Fortunately in Python, "private" just means "harder to modify by mistake", and we can still modify it with a bit of name mangling. Most of the following cell is copied from RMG source code. To see the changes, refer to the git history of this repository.
def storeLabeledProductAtoms(reaction):
"""
Store the labeled product atoms on the reaction so we can recover them later
"""
labeledProductAtoms = []
for product in reaction.products:
for label, atom in product.getLabeledAtoms().items():
labeledProductAtoms.append((label, atom))
reaction.labeledProductAtoms = labeledProductAtoms
def __generateReactions(self, reactants, products=None, forward=True):
from rmgpy.data.base import ForbiddenStructureException
"""
Generate a list of all of the possible reactions of this family between
the list of `reactants`. The number of reactants provided must match
the number of reactants expected by the template, or this function
will return an empty list. Each item in the list of reactants should
be a list of :class:`Molecule` objects, each representing a resonance
isomer of the species of interest.
"""
rxnList = []; speciesList = []
# Wrap each reactant in a list if not already done (this is done to
# allow for passing multiple resonance structures for each molecule)
# This also makes a copy of the reactants list so we don't modify the
# original
reactants = [reactant if isinstance(reactant, list) else [reactant] for reactant in reactants]
if forward:
template = self.forwardTemplate
elif self.reverseTemplate is None:
return []
else:
template = self.reverseTemplate
# Unimolecular reactants: A --> products
if len(reactants) == 1 and len(template.reactants) == 1:
# Iterate over all resonance isomers of the reactant
for molecule in reactants[0]:
mappings = self._KineticsFamily__matchReactantToTemplate(molecule, template.reactants[0])
for map in mappings:
reactantStructures = [molecule]
try:
productStructures = self._KineticsFamily__generateProductStructures(reactantStructures, [map], forward)
except ForbiddenStructureException:
pass
else:
if productStructures is not None:
rxn = self._KineticsFamily__createReaction(reactantStructures, productStructures, forward)
if rxn:
storeLabeledProductAtoms(rxn)
rxnList.append(rxn)
# Bimolecular reactants: A + B --> products
elif len(reactants) == 2 and len(template.reactants) == 2:
moleculesA = reactants[0]
moleculesB = reactants[1]
# Iterate over all resonance isomers of the reactant
for moleculeA in moleculesA:
for moleculeB in moleculesB:
# Reactants stored as A + B
mappingsA = self._KineticsFamily__matchReactantToTemplate(moleculeA, template.reactants[0])
mappingsB = self._KineticsFamily__matchReactantToTemplate(moleculeB, template.reactants[1])
# Iterate over each pair of matches (A, B)
for mapA in mappingsA:
for mapB in mappingsB:
reactantStructures = [moleculeA, moleculeB]
try:
productStructures = self._KineticsFamily__generateProductStructures(reactantStructures, [mapA, mapB], forward)
except ForbiddenStructureException:
pass
else:
if productStructures is not None:
rxn = self._KineticsFamily__createReaction(reactantStructures, productStructures, forward)
if rxn:
storeLabeledProductAtoms(rxn)
rxnList.append(rxn)
# Only check for swapped reactants if they are different
if reactants[0] is not reactants[1]:
# Reactants stored as B + A
mappingsA = self._KineticsFamily__matchReactantToTemplate(moleculeA, template.reactants[1])
mappingsB = self._KineticsFamily__matchReactantToTemplate(moleculeB, template.reactants[0])
# Iterate over each pair of matches (A, B)
for mapA in mappingsA:
for mapB in mappingsB:
reactantStructures = [moleculeA, moleculeB]
try:
productStructures = self._KineticsFamily__generateProductStructures(reactantStructures, [mapA, mapB], forward)
except ForbiddenStructureException:
pass
else:
if productStructures is not None:
rxn = self._KineticsFamily__createReaction(reactantStructures, productStructures, forward)
if rxn:
storeLabeledProductAtoms(rxn)
rxnList.append(rxn)
# If products is given, remove reactions from the reaction list that
# don't generate the given products
if products is not None:
products = [product.generateResonanceIsomers() for product in products]
rxnList0 = rxnList[:]
rxnList = []
index = 0
for reaction in rxnList0:
products0 = reaction.products if forward else reaction.reactants
# Skip reactions that don't match the given products
match = False
if len(products) == len(products0) == 1:
for product in products[0]:
if products0[0].isIsomorphic(product):
match = True
break
elif len(products) == len(products0) == 2:
for productA in products[0]:
for productB in products[1]:
if products0[0].isIsomorphic(productA) and products0[1].isIsomorphic(productB):
match = True
break
elif products0[0].isIsomorphic(productB) and products0[1].isIsomorphic(productA):
match = True
break
elif len(products) == len(products0):
raise NotImplementedError("Can't yet filter reactions with {} products".format(len(products)))
if match:
rxnList.append(reaction)
# The reaction list may contain duplicates of the same reaction
# These duplicates should be combined (by increasing the degeneracy of
# one of the copies and removing the others)
index0 = 0
while index0 < len(rxnList):
reaction0 = rxnList[index0]
products0 = reaction0.products if forward else reaction0.reactants
products0 = [product.generateResonanceIsomers() for product in products0]
# Remove duplicates from the reaction list
index = index0 + 1
while index < len(rxnList):
reaction = rxnList[index]
products = reaction.products if forward else reaction.reactants
# We know the reactants are the same, so we only need to compare the products
match = False
if len(products) == len(products0) == 1:
for product in products0[0]: # for each resonance isomer of the only product0
if products[0].isIsomorphic(product):
match = True
break
elif len(products) == len(products0) == 2:
for productA in products0[0]:
for productB in products0[1]:
if products[0].isIsomorphic(productA) and products[1].isIsomorphic(productB):
match = True
break
elif products[0].isIsomorphic(productB) and products[1].isIsomorphic(productA):
match = True
break
elif len(products) == len(products0) == 3:
for productA, productB, productC in itertools.product(products0[0], products0[1], products0[2]):
# This is equivalent to three nested for loops,
# but allows us to break out of them all at once
# with a single break statement.
if ( products[0].isIsomorphic(productA) and
products[1].isIsomorphic(productB) and
products[2].isIsomorphic(productC) ):
match = True
break
elif ( products[0].isIsomorphic(productA) and
products[1].isIsomorphic(productC) and
products[2].isIsomorphic(productB) ):
match = True
break
elif ( products[0].isIsomorphic(productB) and
products[1].isIsomorphic(productA) and
products[2].isIsomorphic(productC) ):
match = True
break
elif ( products[0].isIsomorphic(productC) and
products[1].isIsomorphic(productA) and
products[2].isIsomorphic(productB) ):
match = True
break
elif ( products[0].isIsomorphic(productB) and
products[1].isIsomorphic(productC) and
products[2].isIsomorphic(productA) ):
match = True
break
elif ( products[0].isIsomorphic(productC) and
products[1].isIsomorphic(productB) and
products[2].isIsomorphic(productA) ):
match = True
break
elif len(products) == len(products0):
raise NotImplementedError(
"Can't yet check degeneracy of reactions with {0} products".format(len(products))
)
# If we found a match, remove it from the list
# Also increment the reaction path degeneracy of the remaining reaction
if match:
rxnList.remove(reaction)
reaction0.degeneracy += 1
else:
index += 1
index0 += 1
# For R_Recombination reactions, the degeneracy is twice what it should
# be, so divide those by two
# This is hardcoding of reaction families!
if self.label.lower().startswith('r_recombination'):
for rxn in rxnList:
assert(rxn.degeneracy % 2 == 0)
rxn.degeneracy /= 2
# Determine the reactant-product pairs to use for flux analysis
# Also store the reaction template (useful so we can easily get the kinetics later)
for reaction in rxnList:
# Restore the labeled atoms long enough to generate some metadata
for reactant in reaction.reactants:
reactant.clearLabeledAtoms()
for label, atom in reaction.labeledAtoms:
atom.label = label
# Generate metadata about the reaction that we will need later
reaction.pairs = self.getReactionPairs(reaction)
reaction.template = self.getReactionTemplateLabels(reaction)
if not forward:
reaction.degeneracy = self.calculateDegeneracy(reaction)
# Unlabel the atoms
for label, atom in reaction.labeledAtoms:
atom.label = ''
# We're done with the labeled atoms, so delete the attribute
# MODIFICATION (commented the following line)
# del reaction.labeledAtoms
# END MOD
# This reaction list has only checked for duplicates within itself, not
# with the global list of reactions
return rxnList
rmgpy.data.kinetics.family.KineticsFamily._KineticsFamily__generateReactions = __generateReactions
class NoMatchingReactionsException(Exception):
"Couldn't find a matching radical recombination reaction."
pass
class MultipleMatchingReactionsException(Exception):
"Found more than one radical recombination reaction that matches."
def __init__(self, message, reactions=None):
super(Exception, self).__init__(self, message)
self.reactions = reactions or []
class CombinedMoleculeIsRadicalException(Exception):
"The product of the matching radical recombination is itself a radical."
def __init__(self, message, species=None):
super(Exception, self).__init__(self, message)
self.species = species
def make_labeled_rmg_reaction(cantera_reaction):
unlabeled_rmg_reaction = make_rmg_reaction(cantera_reaction)
reactant_molecules = [species.molecule for species in unlabeled_rmg_reaction.reactants]
# reactant_molecules is a list of lists of resonance isomers,
# eg. a bimolecular reaction where the second reactant has 2 isomers is: [[r1],[r2i1,r2i2]]
products = [species.molecule[0] for species in unlabeled_rmg_reaction.products]
# products is a list of molecule objects (only one resonance form of each product), eg [p1, p2]
matching_reactions = []
for reactants in itertools.product(*reactant_molecules):
# reactants is now a tuple of molecules, one for each reactant,
# eg. (r1, r2i1), and on the next iteration (r1, r2i2)
matching_reactions.extend(rmgDatabase.kinetics.generateReactionsFromFamilies(reactants, products, only_families=rxnFamilies))
logging.debug("Generated these reactions:")
for reaction in matching_reactions:
logging.debug(reaction)
if len(matching_reactions) == 0:
raise NoMatchingReactionsException(
"Couldn't generate any reactions matching {} in families {}".format(unlabeled_rmg_reaction, rxnFamilies)
)
if len(matching_reactions) > 1:
raise MultipleMatchingReactionsException(
"Generated {} reactions matching {} in family {}".format(len(matching_reactions),unlabeled_rmg_reaction, rxnFamilies),
matching_reactions,
)
labeled_rmg_reaction = matching_reactions[0]
logging.debug("The reaction of interest is as follows: ")
logging.debug(labeled_rmg_reaction)
logging.debug("asserting that the matching reaction is Isomorphic")
assert unlabeled_rmg_reaction.isIsomorphic(labeled_rmg_reaction)
assert len(labeled_rmg_reaction.products) == 1
product = labeled_rmg_reaction.products[0]
if product.isRadical():
raise CombinedMoleculeIsRadicalException(
'The product of the radical recombination reaction {} is itself a radical'.format(product.toSMILES()),
product)
if not unlabeled_rmg_reaction.isIsomorphic(labeled_rmg_reaction,eitherDirection=False):
# Need to reverse the direction
labeled_rmg_reaction.reactants, labeled_rmg_reaction.products = labeled_rmg_reaction.products, labeled_rmg_reaction.reactants
logging.debug("reaction: {!r}".format(labeled_rmg_reaction))
return labeled_rmg_reaction
make_labeled_rmg_reaction(r)
r2 =make_labeled_rmg_reaction(r)
print r2.labeledAtoms
print r2.labeledProductAtoms
no_matching_reactions = OrderedDict()
multiple_matching_reactions = OrderedDict()
parent_is_radical = OrderedDict()
matched_reactions = OrderedDict()
for r in parser.reactions:
if is_radical_recombination(r) or is_unimolecular_dissociation(r):
print '{:50s}'.format(r) ,
try:
labeled_rmg_reaction = make_labeled_rmg_reaction(r)
except NoMatchingReactionsException as e:
no_matching_reactions[r] = make_rmg_reaction(r)
print("❗️ No match")
except MultipleMatchingReactionsException as e:
multiple_matching_reactions[r] = [make_rmg_reaction(r2) for r2 in e.reactions]
print("‼️ Multiple matches")
except CombinedMoleculeIsRadicalException as e:
parent_is_radical[r] = e.species
print("⚠️ Parent is radical")
else:
matched_reactions[r] = labeled_rmg_reaction
print("✅ Unique match")
mprint('# Unmatched reactions ({})\n These reactions are unimolecular dissociation (or bimolecular recombination) but did *not* match any radical recombination in RMG'.format(len(no_matching_reactions)))
for r, r2 in no_matching_reactions.iteritems():
display(r2)
print str(r)
print '-'*40
mprint('# Multiple matching reactions ({})\n These reactions are unimolecular dissociation (or bimolecular recombination) but matched *multiple* radical recombination in RMG'.format(len(multiple_matching_reactions)))
if not multiple_matching_reactions: print "None!"
for r,rlist in multiple_matching_reactions.iteritems():
display(make_rmg_reaction(r))
print str(r)
for rx in rlist:
display(rx)
print '-'*40
mprint(('# Parent molecule is a radical ({})\n'
'These reactions are unimolecular dissociation (or bimolecular recombination) '
'but the unimolecular parent (or product of recombination) is itself a radical '
'so we are excluding them').format(len(parent_is_radical)))
if not parent_is_radical: print "None!"
for r,species in parent_is_radical.iteritems():
display(make_rmg_reaction(r))
print str(r)
print '-'*40
mprint('# Matched reactions ({})\n These reactions each matched a unique radical recombinations in RMG'.format(len(matched_reactions)))
for r,r2 in matched_reactions.iteritems():
display(r2)
print str(r)
print '-'*40
# Testing
display(r2)
print r2.labeledAtoms
print r2.labeledAtoms[0][1].bonds
print r2.reactants[0].getLabeledAtoms()
print r2.labeledProductAtoms
r2.labeledProductAtoms[0][1].bonds
# preserve the previously loaded family(s) but add/replace any new ones
old_families = rmgDatabase.kinetics.families
rmgDatabase.kinetics.loadFamilies(path='.', families=['R_Roaming', 'R_Roaming_alpha'])
old_families.update(rmgDatabase.kinetics.families)
rmgDatabase.kinetics.families = old_families
rmgDatabase.kinetics.families
We have to perform another monkey-patch on RMG to make it use atom labels (if present) to restrict the subgraph isomorphism when matching a species to a reaction template, so that we only generate the roaming reaction that corresponds to our identified homolytic fission reaction (otherwise we'd get all the roamings associated with all the homolytic fissions from the parent molecule).
# Another monkey-patch
def __matchReactantToTemplate(self, reactant, templateReactant):
"""
Return a complete list of the mappings if the provided reactant
matches the provided template reactant, or an empty list if not.
"""
from rmgpy.data.base import LogicNode
from rmgpy.molecule import Group
if isinstance(templateReactant, list): templateReactant = templateReactant[0]
struct = templateReactant.item
if isinstance(struct, LogicNode):
mappings = []
for child_structure in struct.getPossibleStructures(self.groups.entries):
mappings.extend(reactant.findSubgraphIsomorphisms(child_structure))
return mappings
elif isinstance(struct, Group):
## MOD
template_labels = struct.getLabeledAtoms()
initial_mapping = {}
for label, atom in reactant.getLabeledAtoms().items():
initial_mapping[atom] = template_labels[label]
## END MOD
return reactant.findSubgraphIsomorphisms(struct, initial_mapping)
else:
raise NotImplementedError("Not expecting template of type {}".format(type(struct)))
rmgpy.data.kinetics.family.KineticsFamily._KineticsFamily__matchReactantToTemplate = __matchReactantToTemplate
mprint("# Some roaming reactions!")
all_roaming_reactions_parents = {}
all_roaming_reactions_by_parent = defaultdict(list)
for r1,r2 in matched_reactions.iteritems():
roaming = rmgDatabase.kinetics.families['R_Roaming']
roaming_alpha = rmgDatabase.kinetics.families['R_Roaming_alpha']
print "The homolytic fission reaction"
display(r2)
if len(r2.reactants) == 1:
reactant = r2.reactants[0]
reactant.clearLabeledAtoms()
for label, atom in r2.labeledProductAtoms:
atom.label = label
reverse = False
elif len(r2.products) == 1:
raise NotImplementedError, "Not got the label passing fixed in this direction"
reactant = r2.products[0]
for label,atom in r2.labeledAtoms:
atom.label = label
reverse = True
else:
raise Exception
#print reactant.toAdjacencyList()
roaming_reactions = roaming.generateReactions([reactant])
roaming_reactions.extend(roaming_alpha.generateReactions([reactant]))
# Restore the labels in case they got erased
reactant.clearLabeledAtoms()
for label, atom in r2.labeledProductAtoms:
atom.label = label
# now swap the labels
old_labels = reactant.getLabeledAtoms()
old_labels['*1'].label='*2'
old_labels['*2'].label='*1'
roaming_reactions.extend(roaming.generateReactions([reactant]))
roaming_reactions.extend(roaming_alpha.generateReactions([reactant]))
print "Has {} roaming reaction alternative pathways:".format(len(roaming_reactions))
for r in roaming_reactions:
display(r)
if r.family != 'R_Roaming':
print "Generated by reaction family {}".format(r.family)
all_roaming_reactions_parents[r] = r1
all_roaming_reactions_by_parent[r1].append(r)
print '-'*40
n = len(all_roaming_reactions_parents)
m = len(all_roaming_reactions_by_parent)
print "In total there are {} roaming reactions generated from {} parent reactions".format(n,m)
# What types of kinetics do we need to deal with?
kinetics_types = Counter()
for r in all_roaming_reactions_by_parent.keys():
kinetics_types[type(r.kinetics)] += 1
kinetics_types
species_lists_by_formula = defaultdict(list)
for species in speciesList:
species_lists_by_formula[species.molecule[0].getFormula()].append(species)
all_new_species = []
for r, parent in all_roaming_reactions_parents.iteritems():
# convert the reactants from Molecules into Species, although these won't be new
assert len(r.reactants) == 1, "Going in unexpected direction"
molecule = r.reactants[0]
if type(molecule) == rmgpy.species.Species:
# already converted (run this cell before?). Convert back!
molecule = molecule.molecule[0]
formula = molecule.getFormula()
possibles = species_lists_by_formula[formula]
for possible in possibles:
if possible.isIsomorphic(molecule):
species = possible
break
else:
raise Exception("Wasn't expecting a new species as a reactant!")
r.reactants[0] = species
# now convert the products
new_species = []
for i, molecule in enumerate(r.products):
if type(molecule) == rmgpy.species.Species:
# already converted (run this cell before?). Convert back!
molecule = molecule.molecule[0]
formula = molecule.getFormula()
possibles = species_lists_by_formula[formula]
for possible in possibles:
if possible.isIsomorphic(molecule):
species = possible
break
else: # didn't break, so didn't find existing species
species = rmgpy.species.Species(molecule=[molecule])
species.generateResonanceIsomers()
species_lists_by_formula[formula].append(species)
new_species.append(species)
r.products[i] = species
if new_species:
print "The roaming reaction"
display(r)
print "generated {} new species".format(len(new_species))
for s in new_species:
display(s)
print s.molecule[0].getFormula(), s.molecule[0].toSMILES()
print '-'*40
all_new_species.extend(new_species)
del(r.pairs) # no longer points to reagents, and not useful anyway
mprint("## In total there are {} new species".format(len(all_new_species)))
mprint("We will need to get names, thermo, etc. for these")
for s in all_new_species:
display(s)
print s.molecule[0].getFormula(), s.molecule[0].toSMILES()
# Save/add these in a csv file, for use in the other notebook.
csv_filename = 'new-species-{}.csv'.format(model_name)
import pandas as pd
if not os.path.exists(csv_filename):
with open(csv_filename,'w') as f:
print "Making new file {}".format(csv_filename)
f.write('SMILES, Name\n')
df = pd.read_csv(csv_filename, index_col=0)
for s in all_new_species:
smiles = s.molecule[0].toSMILES()
if smiles not in df.index:
print "Adding {} to csv file".format(smiles)
df.loc[smiles] = None
df.to_csv(csv_filename)
df
Now go off and run the second workbook to generate data for the new species, then come back here to read it in.
# Naming new species
df = pd.read_csv(csv_filename, index_col=0)
for s in all_new_species:
smiles = s.molecule[0].toSMILES()
s.label = df.loc[smiles][0]
display(s)
print s.label
new_parser = ck2cti.Parser()
new_thermo_filename = 'new-species-{}-thermo.txt'.format(model_name)
new_transport_filename = 'new-species-{}-transport.txt'.format(model_name)
surfaces = new_parser.convertMech(inputFile=new_thermo_filename,
thermoFile=None,
transportFile=new_transport_filename,
surfaceFile=None,
phaseName='gas',
outName=os.path.join(outdir, model_name+'.new-species.cti'),
permissive=True)
# Copy the new species definitions into the previous parser
parser.speciesDict.update(new_parser.speciesDict)
for s in new_parser.speciesList:
if s not in parser.speciesList:
parser.speciesList.append(s)
parser.species_tokens.update(new_parser.species_tokens)
Now how about finally making some alterative models?
"""
Remove reactions that make products we don't have in the model
because we don't have names, thermo, transport data yet.
For this check, we assume that anything with a label, we have data for,
and anything we don't have data for, does not have a label.
Now that we have generated data for the new species, this cell
should do nothing and report nothing.
"""
to_remove = []
for roaming,parent in all_roaming_reactions_parents.iteritems():
if all([s.label for s in roaming.products]):
# all species are identified and in model
continue # to next reaction
to_remove.append((roaming, parent))
for roaming,parent in to_remove:
print("Removing roaming reaction {} because a product is not in the model".format(roaming))
del(all_roaming_reactions_parents[roaming])
all_roaming_reactions_by_parent[parent].remove(roaming)
def to_cantera(rmg_reaction, cti_parent):
"Make a cantera reaction, and copy across the kinetics"
if not all([s.label for s in rmg_reaction.products]):
raise NotImplementedError("Generated species {} with no data".format(s.molecule[0].toSMILES()))
reactants = [(1, parser.speciesDict[s.label]) for s in rmg_reaction.reactants]
products = [(1, parser.speciesDict[s.label]) for s in rmg_reaction.products]
cti_reaction = ck2cti.Reaction(reactants=reactants, products=products)
cti_reaction.comment = "ROAMING: from parent reaction {} {!s}".format(cti_parent.index, cti_parent)
k = cti_parent.kinetics
if isinstance(k, ck2cti.Arrhenius):
cti_reaction.kinetics = ck2cti.Arrhenius(A=k.A, T0=k.T0, b=k.b, Ea=k.Ea)
elif isinstance(k, ck2cti.Falloff):
klo = k.arrheniusLow
khi = k.arrheniusHigh
cti_reaction.kinetics = ck2cti.Falloff(
arrheniusLow = ck2cti.Arrhenius(A=klo.A, T0=klo.T0, b=klo.b, Ea=klo.Ea, parser=parser),
arrheniusHigh = ck2cti.Arrhenius(A=khi.A, T0=khi.T0, b=khi.b, Ea=khi.Ea, parser=parser),
efficiencies = k.efficiencies,
F = k.F,
)
cti_reaction.thirdBody = cti_parent.thirdBody
elif isinstance(k, ck2cti.PDepArrhenius):
cti_reaction.kinetics = ck2cti.PDepArrhenius(
pressures = k.pressures,
arrhenius = [ck2cti.Arrhenius(A=k0.A, T0=k0.T0, b=k0.b, Ea=k0.Ea, parser=parser) for k0 in k.arrhenius],
parser = parser,
)
cti_reaction.thirdBody = cti_parent.thirdBody
else:
cti_reaction.kinetics = ck2cti.Arrhenius(A=(0, '1/s'), b=0, Ea=(0, 'cal/mol'))
cti_reaction.comment += "\n Parent kinetics of type {} couldn't be copied.".format(type(cti_parent.kinetics))
cti_reaction.kinetics.parser = cti_parent.kinetics.parser
cti_reaction.line_number = cti_parent.line_number
return cti_reaction
We want to store the number of roaming reactions as an attribute number_of_roaming_reactions
on the reaction object, that it knows how to print itself later.
# More monkeypatching, so that it works for FallOff and PLOG reactions
# When you try to set the number of roaming reactions on a Falloff kinetics type
# store that on BOTH of the Arrhenius objects (high and low),
# so they can both report the ALPHA thing
def falloff_roaming_setter(self, value):
self.arrheniusHigh.number_of_roaming_reactions = value
self.arrheniusLow.number_of_roaming_reactions = value
def falloff_roaming_getter(self):
return self.arrheniusHigh.number_of_roaming_reactions
ck2cti.Falloff.number_of_roaming_reactions = property(falloff_roaming_getter, falloff_roaming_setter)
# When you try to set the number of roaming reactions on a PDepArrhenius kinetics type
# store that on EACH of the Arrhenius objects, so they can each report the ALPHA thing
def plog_roaming_setter(self, value):
for k in self.arrhenius:
k.number_of_roaming_reactions = value
def plog_roaming_getter(self):
return self.arrhenius[0].number_of_roaming_reactions
ck2cti.PDepArrhenius.number_of_roaming_reactions = property(plog_roaming_getter, plog_roaming_setter)
# Monkey patch the ck2cti translator
def arrhenius_rateStr(self):
number_of_roaming_reactions = getattr(self, 'number_of_roaming_reactions', 0)
if number_of_roaming_reactions > 0:
Anumber = '(ALPHA/{1:d}.)*{0:e}'.format(self.A[0], number_of_roaming_reactions)
elif number_of_roaming_reactions < 0:
Anumber = '(1-ALPHA)*{0:e}'.format(self.A[0])
else:
Anumber = '{0:e}'.format(self.A[0])
if ck2cti.compatible_quantities(self.parser.output_quantity_units, self.A[1]):
A = '{0}'.format(Anumber)
else:
A = "({0}, '{1}')".format(Anumber, self.A[1])
if 'k' + self.Ea[1] == self.parser.output_energy_units:
# If stored in kcal but output is in cal, just convert
self.Ea[1] = self.parser.output_energy_units
self.Ea[0] *= 1000.0
if self.Ea[1] == self.parser.output_energy_units:
Ea = str(self.Ea[0])
else:
Ea = "({0}, '{1}')".format(*self.Ea)
return '[{0}, {1}, {2}]'.format(A, self.b, Ea)
ck2cti.Arrhenius.rateStr = arrhenius_rateStr
# Another monkey-patch
def pdeparrhenius_to_cti(self, reactantstr, arrow, productstr, indent=0):
"""
This is actually the method in the dev version of cantera,
as of 2017-02-20 (7b7aea2) but not yet in the 2.3.0a3 release.
"""
rxnstring = reactantstr + arrow + productstr
lines = ['pdep_arrhenius({0!r},'.format(rxnstring)]
prefix = ' '*(indent+15)
template = '[({0}, {1!r}), {2}],'
for pressure, arrhenius in zip(self.pressures[0], self.arrhenius):
lines.append(prefix + template.format(pressure,
self.pressures[1],
arrhenius.rateStr()[1:-1]))
lines[-1] = lines[-1][:-1] + ')'
return '\n'.join(lines)
ck2cti.PDepArrhenius.to_cti = pdeparrhenius_to_cti
def print_comment(r):
if r.comment.strip():
for line in r.comment.split('\n'):
print "# "+line
def reactions_are_duplicate(r1,r2):
"""
Compare two cantera-style reactions and
return True iff they're equivalent (ignoring
pressure-depndence, and direction)
"""
if ((Counter(expand_reagent_list(r1.reactants)) == Counter(expand_reagent_list(r2.reactants))) and
(Counter(expand_reagent_list(r1.products)) == Counter(expand_reagent_list(r2.products)))
): return True
if ((Counter(expand_reagent_list(r1.reactants)) == Counter(expand_reagent_list(r2.products))) and
(Counter(expand_reagent_list(r1.products)) == Counter(expand_reagent_list(r2.reactants)))
): return True
return False
def find_duplicate(reaction, list_of_reactions=parser.reactions):
"""
If the reaction has a duplicate already in the list_of_reactions,
returns that duplicate, else None.
"""
# no attempt at efficiency here
for r in list_of_reactions:
if reactions_are_duplicate(r,reaction):
return r
else:
return None
# convert to cti
rmg_to_cti_reaction_dict = {}
for rmg_reaction, cti_parent in all_roaming_reactions_parents.items():
cantera_reaction = to_cantera(rmg_reaction, cti_parent)
rmg_to_cti_reaction_dict[rmg_reaction] = cantera_reaction
print "Converted {} reactions into cantera objects".format(len(rmg_to_cti_reaction_dict))
# Detect duplicates
mprint("## Duplicates\nthat you should manually remove from the CTI file if not intended.")
roaming_duplicates = OrderedDict()
for rmg_reaction, cti_parent in all_roaming_reactions_parents.items():
cantera_reaction = rmg_to_cti_reaction_dict[rmg_reaction]
dup = find_duplicate(cantera_reaction, list_of_reactions=parser.reactions)
if dup:
cantera_reaction.comment += '\n This is a DUPLICATE of original rxn {} {!s}'.format(dup.index, dup)
dup.comment += '\n This DUPLICATED by new roaming rxn {!s}'.format(cantera_reaction)
roaming_duplicates[cantera_reaction] = dup
display(rmg_reaction)
cantera_reaction.duplicate = True
dup.duplicate = True
print "{!s} is a duplicate of {}:".format(cantera_reaction, dup.index, dup)
print dup.to_cti()
rmg_to_cti_reaction_dict[rmg_reaction] = cantera_reaction
# print cantera_reaction.comment
# print cantera_reaction.to_cti()
print "\nOf the {} roaming reactions, {} are duplicates of things already there".format(len(all_roaming_reactions_parents),len(roaming_duplicates))
mprint("## Roaming pathways, from different parents, that end up duplicates.")
new_roaming_reactions = rmg_to_cti_reaction_dict.values()
for i, cantera_reaction in enumerate(new_roaming_reactions):
dup = find_duplicate(cantera_reaction, list_of_reactions=new_roaming_reactions[i+1:])
if dup:
cantera_reaction.comment += '\n This is a DUPLICATE of another roaming pathway {!s}'.format(dup)
cantera_reaction.duplicate = True
dup.comment += '\n This is a DUPLICATE of another roaming pathway {!s}'.format(dup)
dup.duplicate = True
print "{!s} is a duplicate:".format(cantera_reaction)
print_comment(cantera_reaction)
print cantera_reaction.to_cti()
print_comment(dup)
print dup.to_cti()
print
for index, cti_reaction in enumerate(parser.reactions):
roaming = all_roaming_reactions_by_parent[cti_reaction]
if roaming:
cti_reaction.comment += "\nThe following reaction has generated these roaming pathways:"
for r in roaming:
cti_reaction.comment +="\n {0!s}".format(r)
cti_reaction.kinetics.number_of_roaming_reactions = -1*len(roaming)
print_comment(cti_reaction)
print cti_reaction.to_cti()
is_arrhenius = isinstance(cti_reaction.kinetics, ck2cti.Arrhenius)
for rmg_reaction in roaming:
roaming_reaction = rmg_to_cti_reaction_dict[rmg_reaction]
roaming_reaction.kinetics.number_of_roaming_reactions = len(roaming)
print_comment(roaming_reaction)
print roaming_reaction.to_cti()
parser.reactions.append(roaming_reaction)
header = ["# For roaming:","ALPHA = 0.1",""]
parser.writeCTI(outName=os.path.join(outdir,model_name+'.roaming.cti'),
header = header,
)
parser.checkDuplicateReactions()