Skip to content

Commit

Permalink
Find multiple transition states for reactions that match >1 reaction …
Browse files Browse the repository at this point in the history 
…templates.

Preliminary fix for #142

Instead of calculating degeneracies first, find reaction templates after
generating reactions, then combine degenerate reactions only if reaction
templates match. These isomorphic reactions remain separate
(unlike in RMG-Java) in case someday we want to calculate the separate
transition states.  The reactions are added to the core and printed
to chemkin as duplicates.

Issues still remaining:
- How to identify the correct reverse reaction when having multiple TS's.
 My solution is to match the correct reaction pairs- but this only works in
the case of a bimolecular A+B=C+D reaction.  If the reaction is unimolecular,
the reaction pairs will be the same and the correct path cannot be distinguished.
- If kinetics are identical, they should be combined through degeneracy, in this
commit they remain as duplicates.
  • Loading branch information
@connie
connie committed Feb 3, 2015
1 parent 160ad1b commit 9ea67ff
Show file tree
Hide file tree
Showing 2 changed files with 76 additions and 37 deletions.
106 changes: 71 additions & 35 deletions rmgpy/data/kinetics/family.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1221,16 +1221,27 @@ def generateReactions(self, reactants, failsSpeciesConstraints=None):
# for each reaction, make its reverse reaction and store in a 'reverse' attribute
for rxn in reactionList:
reactions = self.__generateReactions(rxn.products, products=rxn.reactants, forward=True, failsSpeciesConstraints=failsSpeciesConstraints)
#assert len(reactions) == 1, "Expecting one matching reverse reaction, not {0}. Forward reaction {1!s} : {1!r}".format(len(reactions), rxn)
if len(reactions) != 1:
logging.error("Expecting one matching reverse reaction, not {0} in reaction family {1} for forward reaction {2}.\n".format(len(reactions), self.label, str(rxn)))
for reactant in rxn.reactants:
logging.info("Reactant")
logging.info(reactant.toAdjacencyList())
for product in rxn.products:
logging.info("Product")
logging.info(product.toAdjacencyList())
raise KineticsError("Did not find reverse reaction in reaction family {0} for reaction {1}.".format(self.label, str(rxn)))
rxn.reverse = reactions[0]
for possibleReaction in reactions:
if possibleReaction.pairs[0][0].isIsomorphic(rxn.pairs[0][0]) and possibleReaction.pairs[0][1].isIsomorphic(rxn.pairs[0][1]):
rxn.reverse = possibleReaction
break
if possibleReaction.pairs[0][0].isIsomorphic(rxn.pairs[0][1]) and possibleReaction.pairs[0][1].isIsomorphic(rxn.pairs[0][0]):
rxn.reverse = possibleReaction
break
if possibleReaction.pairs[0][0].isIsomorphic(rxn.pairs[1][0]) and possibleReaction.pairs[0][1].isIsomorphic(rxn.pairs[1][1]):
rxn.reverse = possibleReaction
break
if possibleReaction.pairs[0][0].isIsomorphic(rxn.pairs[1][1]) and possibleReaction.pairs[0][1].isIsomorphic(rxn.pairs[1][0]):
rxn.reverse = possibleReaction
break
else:
raise Exception("Could not match to reverse reaction due to no matching pairs. Forward reaction {1!s} : {1!r}".format(len(reactions), rxn))

else:
# Only one matching reverse reaction
rxn.reverse = reactions[0]

else: # family is not ownReverse
# Reverse direction (the direction in which kinetics is not defined)
Expand All @@ -1256,12 +1267,23 @@ def calculateDegeneracy(self, reaction):
"""
reactions = self.__generateReactions(reaction.reactants, products=reaction.products, forward=True)
if len(reactions) != 1:
for reactant in reaction.reactants:
logging.error(reactant)
for product in reaction.products:
logging.error(product)
if not reaction.pairs:
reaction.pairs = self.getReactionPairs(reaction)
for possibleReaction in reactions:
# Match the correct reaction if there are multiple reaction pathways
if possibleReaction.pairs[0][0].isIsomorphic(reaction.pairs[0][0]) and possibleReaction.pairs[0][1].isIsomorphic(reaction.pairs[0][1]):
return possibleReaction.degeneracy
if possibleReaction.pairs[0][0].isIsomorphic(reaction.pairs[0][1]) and possibleReaction.pairs[0][1].isIsomorphic(reaction.pairs[0][0]):
return possibleReaction.degeneracy
if possibleReaction.pairs[0][0].isIsomorphic(reaction.pairs[1][0]) and possibleReaction.pairs[0][1].isIsomorphic(reaction.pairs[1][1]):
return possibleReaction.degeneracy
if possibleReaction.pairs[0][0].isIsomorphic(reaction.pairs[1][1]) and possibleReaction.pairs[0][1].isIsomorphic(reaction.pairs[1][0]):
return possibleReaction.degeneracy

raise Exception('Unable to calculate degeneracy for reaction {0} in reaction family {1}.'.format(reaction, self.label))
return reactions[0].degeneracy

else:
return reactions[0].degeneracy

def __generateReactions(self, reactants, products=None, forward=True, failsSpeciesConstraints=None):
"""
Expand Down Expand Up @@ -1396,6 +1418,29 @@ def __generateReactions(self, reactants, products=None, forward=True, failsSpeci
if match:
rxnList.append(reaction)



# 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.getReactionTemplate(reaction)

# Unlabel the atoms
for label, atom in reaction.labeledAtoms:
atom.label = ''

# We're done with the labeled atoms, so delete the attribute
del reaction.labeledAtoms

# 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)
Expand Down Expand Up @@ -1433,45 +1478,36 @@ def __generateReactions(self, reactants, products=None, forward=True, failsSpeci
# 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
if reaction0.template == reaction.template:
# template must match, otherwise we may have found an alternate transition state
reaction0.degeneracy += 1
rxnList.remove(reaction)
else:
index += 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!
# For reactions of the form A + A -> products, the degeneracy is twice
# what it should be, so divide those by two
if sameReactants or self.label.lower().startswith('r_recombination'):
for rxn in rxnList:
print rxn
print rxn.template
print rxn.degeneracy
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.getReactionTemplate(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
del reaction.labeledAtoms


# This reaction list has only checked for duplicates within itself, not
# with the global list of reactions
Expand Down
7 changes: 5 additions & 2 deletions rmgpy/rmg/model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -483,11 +483,14 @@ def checkForExistingReaction(self, rxn):
if not rxn.duplicate:
return True, rxn0
else:
return True, rxn0
if rxn0.template == rxn.template:
print 'here'
return True, rxn0

if isinstance(family,KineticsFamily) and family.ownReverse:
if (rxn0.reactants == rxn.products and rxn0.products == rxn.reactants):
return True, rxn0
if rxn0.template == rxn.template:
return True, rxn0

# Now check seed mechanisms
# We want to check for duplicates in *other* seed mechanisms, but allow
Expand Down

0 comments on commit 9ea67ff

Please sign in to comment.