text: #!/usr/bin/python
  # -*- coding: utf-8 -*-
  
  ################################################################################
  #
  #   RMG - Reaction Mechanism Generator
  #
  #   Copyright (c) 2002-2017 Prof. William H. Green (whgreen@mit.edu), 
  #   Prof. Richard H. West (r.west@neu.edu) and the RMG Team (rmg_dev@mit.edu)
  #
  #   Permission is hereby granted, free of charge, to any person obtaining a
  #   copy of this software and associated documentation files (the 'Software'),
  #   to deal in the Software without restriction, including without limitation
  #   the rights to use, copy, modify, merge, publish, distribute, sublicense,
  #   and/or sell copies of the Software, and to permit persons to whom the
  #   Software is furnished to do so, subject to the following conditions:
  #
  #   The above copyright notice and this permission notice shall be included in
  #   all copies or substantial portions of the Software.
  #
  #   THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  #   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  #   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  #   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  #   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  #   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  #   DEALINGS IN THE SOFTWARE.
  #
  ################################################################################
  
  """
  Contains classes for working with the reaction model generated by RMG.
  """
  
  import logging
  import math
  import numpy
  import itertools
  import gc
  
  from rmgpy.display import display
  #import rmgpy.chemkin
  import rmgpy.constants as constants
  from rmgpy.constraints import failsSpeciesConstraints
  from rmgpy.quantity import Quantity
  from rmgpy.species import Species
  from rmgpy.thermo.thermoengine import submit
  from rmgpy.reaction import Reaction
  from rmgpy.exceptions import ForbiddenStructureException
  from rmgpy.data.kinetics.depository import DepositoryReaction
  from rmgpy.data.kinetics.family import KineticsFamily, TemplateReaction
  from rmgpy.data.kinetics.library import KineticsLibrary, LibraryReaction
  
  from rmgpy.kinetics import KineticsData
  import rmgpy.data.rmg
  from .react import reactAll
  
  from pdep import PDepReaction, PDepNetwork
  
  # generateThermoDataFromQM under the Species class imports the qm package
  
  ################################################################################
  
  class ReactionModel:
      """
      Represent a generic reaction model. A reaction model consists of `species`,
      a list of species, and `reactions`, a list of reactions.
      """
  
      def __init__(self, species=None, reactions=None):
          self.species = species or []
          self.reactions = reactions or []
      
      def __reduce__(self):
          """
          A helper function used when pickling an object.
          """
          return (ReactionModel, (self.species, self.reactions))
  
      def merge(self, other):
          """
          Return a new :class:`ReactionModel` object that is the union of this
          model and `other`.
          """
          if not isinstance(other, ReactionModel):
              raise ValueError('Expected type ReactionModel for other parameter, got {0}'.format(other.__class__))
  
          # Initialize the merged model
          finalModel = ReactionModel()
          
          # Put the current model into the merged model as-is
          finalModel.species.extend(self.species)
          finalModel.reactions.extend(self.reactions)
          
          # Determine which species in other are already in self
          commonSpecies = {}; uniqueSpecies = []
          for spec in other.species:
              for spec0 in finalModel.species:
                  if spec.isIsomorphic(spec0):
                      commonSpecies[spec] = spec0
                      if spec0.label not in ['Ar','N2','Ne','He']:
                          if not spec0.thermo.isIdenticalTo(spec.thermo):
                              print 'Species {0} thermo from model 1 did not match that of model 2.'.format(spec.label)
                          
                      break
              else:
                  uniqueSpecies.append(spec)
          
          # Determine which reactions in other are already in self
          commonReactions = {}; uniqueReactions = []
          for rxn in other.reactions:
              for rxn0 in finalModel.reactions:
                  if rxn.isIsomorphic(rxn0, eitherDirection=True):
                      commonReactions[rxn] = rxn0                    
                      if not rxn0.kinetics.isIdenticalTo(rxn.kinetics):
                          print 'Reaction {0} kinetics from model 1 did not match that of model 2.'.format(str(rxn0))
                      break
              else:
                  uniqueReactions.append(rxn)
          
          # Add the unique species from other to the final model
          finalModel.species.extend(uniqueSpecies)
      
          # Renumber the unique species (to avoid name conflicts on save)
          speciesIndex = 0
          for spec in finalModel.species:
              if spec.label not in ['Ar','N2','Ne','He']:
                  spec.index = speciesIndex + 1
                  speciesIndex += 1
          
          # Make sure unique reactions only refer to species in the final model
          for rxn in uniqueReactions:
              for i, reactant in enumerate(rxn.reactants):
                  try:
                      rxn.reactants[i] = commonSpecies[reactant]
                      if rxn.pairs:
                          for j, pair in enumerate(rxn.pairs):
                              if reactant in pair:
                                  rxn.pairs[j] = (rxn.reactants[i],pair[1])
                  except KeyError:
                      pass
              for i, product in enumerate(rxn.products):
                  try:
                      rxn.products[i] = commonSpecies[product]
                      if rxn.pairs:
                          for j, pair in enumerate(rxn.pairs):
                              if product in pair:
                                  rxn.pairs[j] = (pair[0], rxn.products[i])
                  except KeyError:
                      pass
          
          # Add the unique reactions from other to the final model
          finalModel.reactions.extend(uniqueReactions)
      
          # Return the merged model
          return finalModel
  
  ################################################################################
  
  class CoreEdgeReactionModel:
      """
      Represent a reaction model constructed using a rate-based screening
      algorithm. The species and reactions in the model itself are called the
      *core*; the species and reactions identified as candidates for inclusion in
      the model are called the *edge*. The attributes are:
  
      =========================  ==============================================================
      Attribute                  Description
      =========================  ==============================================================
      `core`                     The species and reactions of the current model core
      `edge`                     The species and reactions of the current model edge
      `networkDict`              A dictionary of pressure-dependent reaction networks (:class:`Network` objects) indexed by source.
      `networkList`              A list of pressure-dependent reaction networks (:class:`Network` objects)
      `networkCount`             A counter for the number of pressure-dependent networks created
      `indexSpeciesDict`         A dictionary with a unique index pointing to the species objects
      `solventName`              String describing solvent name for liquid reactions. Empty for non-liquid estimation
      =========================  ==============================================================
  
  
      """
  
      def __init__(self, core=None, edge=None, surface=None):
          if core is None:
              self.core = ReactionModel()
          else:
              self.core = core
          if edge is None:
              self.edge = ReactionModel()
          else:
              self.edge = edge
          if surface is None:
              self.surface = ReactionModel()
          else:
              self.surface = surface
              
          # The default tolerances mimic the original RMG behavior; no edge
          # pruning takes place, and the simulation is interrupted as soon as
          # a species flux higher than the validity
          self.networkDict = {}
          self.networkList = []
          self.networkCount = 0
          self.speciesDict = {}
          self.reactionDict = {}
          self.speciesCache = [None for i in range(4)]
          self.speciesCounter = 0
          self.reactionCounter = 0
          self.newSpeciesList = []
          self.newReactionList = []
          self.outputSpeciesList = []
          self.outputReactionList = []
          self.pressureDependence = None
          self.quantumMechanics = None
          self.verboseComments = False
          self.kineticsEstimator = 'group additivity'
          self.indexSpeciesDict = {}
          self.saveEdgeSpecies = False
          self.iterationNum = 0
          self.toleranceThermoKeepSpeciesInEdge = numpy.inf
          self.Gfmax = numpy.inf
          self.Gmax = numpy.inf
          self.Gmin = -numpy.inf
          self.minCoreSizeForPrune = 50
          self.maximumEdgeSpecies = 100000
          self.Tmax = 0
          self.reactionSystems = []
          self.newSurfaceSpcsAdd = set()
          self.newSurfaceRxnsAdd = set()
          self.newSurfaceSpcsLoss = set()
          self.newSurfaceRxnsLoss = set()
          self.solventName = ''
  
      def checkForExistingSpecies(self, molecule):
          """
          Check to see if an existing species contains the same
          :class:`molecule.Molecule` as `molecule`. Returns ``True`` 
          and the matched species (if found) or
          ``False`` and ``None`` (if not found).
          """
          # Create obj to check against existing species
          # obj can be `Molecule` object or `Species` object
  
          # For non-cyclic molecules, obj is `Molecule` object
          # We expect it to be part of the list of isomers in a species
          # object if it has a match
          obj = molecule
  
          # For cyclic molecules, obj is `Species` object and aromatic resonance
          # isomers are generated.  This is due to the hysteresis of isomer generation
          # for aromatic/polyaromatic compounds: not all kekulized forms can be found
          # within the list of isomers for a species object describing a unique aromatic compound
          if molecule.isCyclic():
              obj = Species(molecule=[molecule])
              from rmgpy.molecule.resonance import generate_aromatic_resonance_structures
              aromaticIsomers = generate_aromatic_resonance_structures(molecule)
              obj.molecule.extend(aromaticIsomers)
  
          # First check cache and return if species is found
          for i, spec in enumerate(self.speciesCache):
              if spec is not None:
                  for mol in spec.molecule:
                      if obj.isIsomorphic(mol):
                          self.speciesCache.pop(i)
                          self.speciesCache.insert(0, spec)
                          return True, spec
  
          # Return an existing species if a match is found
          formula = molecule.getFormula()
          try:
              speciesList = self.speciesDict[formula]
          except KeyError:
              return False, None
          for spec in speciesList:
              if spec.isIsomorphic(obj):
                  self.speciesCache.pop()
                  self.speciesCache.insert(0, spec)
                  return True, spec
          # At this point we can conclude that the structure does not exist
          return False, None
  
      def makeNewSpecies(self, object, label='', reactive=True, checkForExisting=True):
          """
          Formally create a new species from the specified `object`, which can be
          either a :class:`Molecule` object or an :class:`rmgpy.species.Species`
          object.
          """
  
          if isinstance(object, rmgpy.species.Species):
              molecule = object.molecule[0]
              label = label if label != '' else object.label
              reactive = object.reactive
          else:
              molecule = object
              
          molecule.clearLabeledAtoms()
  
          # If desired, check to ensure that the species is new; return the
          # existing species if not new
          if checkForExisting:
              found, spec = self.checkForExistingSpecies(molecule)
              if found: return spec, False
  
          # Check that the structure is not forbidden
  
          # If we're here then we're ready to make the new species
          if reactive:
              self.speciesCounter += 1   # count only reactive species
              speciesIndex = self.speciesCounter
          else:
              speciesIndex = -1
          try:
              spec = Species(index=speciesIndex, label=label, molecule=[molecule], reactive=reactive,
                   thermo=object.thermo, transportData=object.transportData)
          except AttributeError, e:
              spec = Species(index=speciesIndex, label=label, molecule=[molecule], reactive=reactive)
          
          spec.creationIteration = self.iterationNum
          spec.generate_resonance_structures()
          spec.molecularWeight = Quantity(spec.molecule[0].getMolecularWeight()*1000.,"amu")
          
          submit(spec,self.solventName)
          
          if spec.label == '':
              if spec.thermo and spec.thermo.label != '': #check if thermo libraries have a name for it
                  logging.info('Species with SMILES of {0} named {1} based on thermo library name'.format(molecule.toSMILES().replace('/','').replace('\\',''),spec.thermo.label))
                  spec.label = spec.thermo.label
                  label = spec.label
              else:
                  # Use SMILES as default format for label
                  # However, SMILES can contain slashes (to describe the
                  # stereochemistry around double bonds); since RMG doesn't 
                  # distinguish cis and trans isomers, we'll just strip these out
                  # so that we can use the label in file paths
                  label = molecule.toSMILES().replace('/','').replace('\\','')
                  
          logging.debug('Creating new species {0}'.format(label))
          
          spec.generateEnergyTransferModel()
          formula = molecule.getFormula()
          if formula in self.speciesDict:
              self.speciesDict[formula].append(spec)
          else:
              self.speciesDict[formula] = [spec]
  
  
          # Since the species is new, add it to the list of new species
          self.newSpeciesList.append(spec)
  
          if spec.reactive:
              self.indexSpeciesDict[spec.index] = spec
  
          return spec, True
  
      def checkForExistingReaction(self, rxn):
          """
          Check to see if an existing reaction has the same reactants, products, and
          family as `rxn`. Returns :data:`True` or :data:`False` and the matched
          reaction (if found).
  
          First, a shortlist of reaction is retrieved that have the same reaction keys
          as the parameter reaction.
  
          Next, the reaction ID containing an identifier (e.g. label) of the reactants
          and products is compared between the parameter reaction and the each of the
          reactions in the shortlist. If a match is found, the discovered reaction is 
          returned.
  
          If a match is not yet found, the Library (seed mechs, reaction libs)
          in the reaction database are iterated over to check if a reaction was overlooked
          (a reaction with a different "family" key as the parameter reaction).
  
          """
  
          # Make sure the reactant and product lists are sorted before performing the check
          rxn.reactants.sort()
          rxn.products.sort()
  
          # If reactants and products are identical, then something weird happened along
          # the way and we got a symmetrical reaction.
          if rxn.reactants == rxn.products:
              logging.debug("Symmetrical reaction found. Returning no reaction")
              return True, None
          
          familyObj = getFamilyLibraryObject(rxn.family)
          shortlist = self.searchRetrieveReactions(rxn)
  
          # Now use short-list to check for matches. All should be in same forward direction.
  
          # Make sure the reactant and product lists are sorted before performing the check
          rxn_id = generateReactionId(rxn)
  
          for rxn0 in shortlist:
              rxn_id0 = generateReactionId(rxn0)
  
              if (rxn_id == rxn_id0) and isinstance(familyObj, KineticsLibrary):
                  # If the reaction comes from a kinetics library, then we can 
                  # retain duplicates if they are marked
                  if areIdenticalSpeciesReferences(rxn, rxn0) and not rxn.duplicate:
                      return True, rxn0
              elif ((rxn_id == rxn_id0) or (rxn_id == rxn_id0[::-1])) and \
                          isinstance(familyObj, KineticsFamily):
                  # ensure TemplateReactions have the same templates and families in order
                  # to classify this as existing reaction. Also checks for reverse
                  # direction matching. Marks duplicate if identical species and different
                  # templates or families
                  if areIdenticalSpeciesReferences(rxn, rxn0):
                      if rxn.family == rxn0.family:
                          equal_templates = frozenset(rxn.template) == frozenset(rxn0.template)
                          # check reverse template
                          if not equal_templates and familyObj.ownReverse and \
                                      rxn.reverse is not None:
                              equal_templates = frozenset(rxn.reverse.template) == frozenset(rxn0.template)
                          if equal_templates:
                              return True, rxn0
                          else:
                              rxn.duplicate = True
                              rxn0.duplicate = True
                      else:
                          rxn.duplicate = True
                          rxn0.duplicate = True
              elif (rxn_id == rxn_id0):
                  if areIdenticalSpeciesReferences(rxn, rxn0):
                      return True, rxn0
  
          # Now check seed mechanisms
          # We want to check for duplicates in *other* seed mechanisms, but allow
          # duplicated *within* the same seed mechanism
          _, r1_fwd, r2_fwd = generateReactionKey(rxn)
          _, r1_rev, r2_rev = generateReactionKey(rxn, useProducts=True)
  
          for library in self.reactionDict:
              libObj = getFamilyLibraryObject(library)
              if isinstance(libObj, KineticsLibrary) and library != rxn.family:
  
                  # First check seed short-list in forward direction                
                  shortlist = self.retrieve(library, r1_fwd, r2_fwd)
                  
                  for rxn0 in shortlist:
                      rxn_id0 = generateReactionId(rxn0)
                      if (rxn_id == rxn_id0) or (rxn_id == rxn_id0[::-1]):
                          if areIdenticalSpeciesReferences(rxn, rxn0):
                              return True, rxn0
                  
                  # Now get the seed short-list of the reverse reaction
  
                  shortlist = self.retrieve(library, r1_rev, r2_rev)
                  
                  for rxn0 in shortlist:
                      if areIdenticalSpeciesReferences(rxn, rxn0):
                          return True, rxn0
  
          return False, None
  
      def makeNewReaction(self, forward, checkExisting=True):
          """
          Make a new reaction given a :class:`Reaction` object `forward`. 
          The reaction is added to the global list of reactions.
          Returns the reaction in the direction that corresponds to the
          estimated kinetics, along with whether or not the reaction is new to the
          global reaction list.
  
          The forward direction is determined using the "is_reverse" attribute of the
          reaction's family.  If the reaction family is its own reverse, then it is
          made such that the forward reaction is exothermic at 298K.
          
          The forward reaction is appended to self.newReactionList if it is new.
          """
  
          # Determine the proper species objects for all reactants and products
          reactants = [self.makeNewSpecies(reactant)[0] for reactant in forward.reactants]
          products  = [self.makeNewSpecies(product)[0]  for product  in forward.products ]
          if forward.specificCollider is not None:
              forward.specificCollider = self.makeNewSpecies(forward.specificCollider)[0]
  
          if forward.pairs is not None:
              for pairIndex in range(len(forward.pairs)):
                  reactantIndex = forward.reactants.index(forward.pairs[pairIndex][0])
                  productIndex = forward.products.index(forward.pairs[pairIndex][1])
                  forward.pairs[pairIndex] = (reactants[reactantIndex], products[productIndex])
                  if hasattr(forward, 'reverse'):                   
                      if forward.reverse:
                          forward.reverse.pairs[pairIndex] = (products[productIndex], reactants[reactantIndex])
          forward.reactants = reactants
          forward.products  = products
  
          if checkExisting:
              found, rxn = self.checkForExistingReaction(forward)
              if found: return rxn, False
  
          # Generate the reaction pairs if not yet defined
          if forward.pairs is None:
              forward.generatePairs()
              if hasattr(forward, 'reverse'):
                  if forward.reverse:
                      forward.reverse.generatePairs()
              
          # Note in the log
          if isinstance(forward, TemplateReaction):
              logging.debug('Creating new {0} template reaction {1}'.format(forward.family, forward))
          elif isinstance(forward, DepositoryReaction):
              logging.debug('Creating new {0} reaction {1}'.format(forward.getSource(), forward))
          elif isinstance(forward, LibraryReaction):
              logging.debug('Creating new library reaction {0}'.format(forward))
          else:
              raise Exception("Unrecognized reaction type {0!s}".format(forward.__class__))
          
          self.registerReaction(forward)
  
          forward.index = self.reactionCounter + 1
          self.reactionCounter += 1
  
          # Since the reaction is new, add it to the list of new reactions
          self.newReactionList.append(forward)
  
          # Return newly created reaction
          return forward, True
  
      def makeNewPDepReaction(self, forward):
          """
          Make a new pressure-dependent reaction based on a list of `reactants` and a
          list of `products`. The reaction belongs to the specified `network` and
          has pressure-dependent kinetics given by `kinetics`.
  
          No checking for existing reactions is made here. The returned PDepReaction
          object is not added to the global list of reactions, as that is intended
          to represent only the high-pressure-limit set. The reactionCounter is
          incremented, however, since the returned reaction can and will exist in
          the model edge and/or core.
          """
  
          # Don't create reverse reaction: all such reactions are treated as irreversible
          # The reverse direction will come from a different partial network
          # Note that this isn't guaranteed to satisfy thermodynamics (but will probably be close)
          forward.reverse = None
          forward.reversible = False
  
          # Generate the reaction pairs if not yet defined
          if forward.pairs is None:
              forward.generatePairs()
              
          # Set reaction index and increment the counter
          forward.index = self.reactionCounter + 1
          self.reactionCounter += 1
  
          return forward
  
      def enlarge(self, newObject=None, reactEdge=False, unimolecularReact=None, bimolecularReact=None):
          """
          Enlarge a reaction model by processing the objects in the list `newObject`. 
          If `newObject` is a
          :class:`rmg.species.Species` object, then the species is moved from
          the edge to the core and reactions generated for that species, reacting
          with itself and with all other species in the model core. If `newObject`
          is a :class:`rmg.unirxn.network.Network` object, then reactions are
          generated for the species in the network with the largest leak flux.
  
          If the `reactEdge` flag is `True`, then no newObject is needed,
          and instead the algorithm proceeds to react the core species together
          to form edge reactions.
          """
          
          numOldCoreSpecies = len(self.core.species)
          numOldCoreReactions = len(self.core.reactions)
          numOldEdgeSpecies = len(self.edge.species)
          numOldEdgeReactions = len(self.edge.reactions)
          reactionsMovedFromEdge = []
          self.newReactionList = []; self.newSpeciesList = []
  
          if reactEdge is False:
              # We are adding core species 
              newReactions = []
              pdepNetwork = None
              objectWasInEdge = False
          
              if isinstance(newObject, Species):
                  
                  newSpecies = newObject
  
                  objectWasInEdge = newSpecies in self.edge.species
                  
                  if not newSpecies.reactive:
                      logging.info('NOT generating reactions for unreactive species {0}'.format(newSpecies))
                  else:
                      logging.info('Adding species {0} to model core'.format(newSpecies))
                      display(newSpecies) # if running in IPython --pylab mode, draws the picture!
  
                  # Add new species
                  reactionsMovedFromEdge = self.addSpeciesToCore(newSpecies)
  
              elif isinstance(newObject, tuple) and isinstance(newObject[0], PDepNetwork) and self.pressureDependence:
  
                  pdepNetwork, newSpecies = newObject
                  newReactions.extend(pdepNetwork.exploreIsomer(newSpecies))
  
                  for rxn in newReactions:
                      rxn = self.inflate(rxn)
                      try:
                          rxn.reverse = self.inflate(rxn.reverse)
                      except AttributeError, e:
                          pass
                      
                  self.processNewReactions(newReactions, newSpecies, pdepNetwork)
  
              else:
                  raise TypeError('Unable to use object {0} to enlarge reaction model; expecting an object of class rmg.model.Species or rmg.model.PDepNetwork, not {1}'.format(newObject, newObject.__class__))
  
              # If there are any core species among the unimolecular product channels
              # of any existing network, they need to be made included
              for network in self.networkList:
                  network.updateConfigurations(self)
                  index = 0
                  while index < len(self.core.species):
                      species = self.core.species[index]
                      isomers = [isomer.species[0] for isomer in network.isomers]
                      if species in isomers and species not in network.explored:
                          network.explored.append(species)
                          continue
                      for products in network.products:
                          products = products.species
                          if len(products) == 1 and products[0] == species:
                              newReactions = network.exploreIsomer(species)
                              for rxn in newReactions:
                                  rxn = self.inflate(rxn)
                                  try:
                                      rxn.reverse = self.inflate(rxn.reverse)
                                  except AttributeError, e:
                                      pass
  
                              self.processNewReactions(newReactions, species, network)
                              network.updateConfigurations(self)
                              index = 0
                              break
                      else:
                          index += 1
              
              if isinstance(newObject, Species) and objectWasInEdge:
                  # moved one species from edge to core
                  numOldEdgeSpecies -= 1
                  # moved these reactions from edge to core
                  numOldEdgeReactions -= len(reactionsMovedFromEdge)
  
          else:
              # We are reacting the edge
  
              rxns = reactAll(self.core.species, numOldCoreSpecies, unimolecularReact, bimolecularReact)
              spcs = [self.retrieveNewSpecies(rxn) for rxn in rxns]
              
              for rxn, spc in zip(rxns, spcs):
                  rxn = self.inflate(rxn) 
                  try:
                      rxn.reverse = self.inflate(rxn.reverse)
                  except AttributeError, e:
                      pass
                  self.processNewReactions([rxn], spc)
  
          ################################################################
          # Begin processing the new species and reactions
          
          # Generate kinetics of new reactions
          logging.info('Generating kinetics for new reactions...')
          for reaction in self.newReactionList:
              # If the reaction already has kinetics (e.g. from a library),
              # assume the kinetics are satisfactory
              if reaction.kinetics is None:
                  self.applyKineticsToReaction(reaction)
                      
          # For new reactions, convert ArrheniusEP to Arrhenius, and fix barrier heights.
          # self.newReactionList only contains *actually* new reactions, all in the forward direction.
          for reaction in self.newReactionList:
              # convert KineticsData to Arrhenius forms
              if isinstance(reaction.kinetics, KineticsData):
                  reaction.kinetics = reaction.kinetics.toArrhenius()
              #  correct barrier heights of estimated kinetics
              if isinstance(reaction,TemplateReaction) or isinstance(reaction,DepositoryReaction): # i.e. not LibraryReaction
                  reaction.fixBarrierHeight() # also converts ArrheniusEP to Arrhenius.
                  
              if self.pressureDependence and reaction.isUnimolecular():
                  # If this is going to be run through pressure dependence code,
                  # we need to make sure the barrier is positive.
                  reaction.fixBarrierHeight(forcePositive=True)
              
          # Update unimolecular (pressure dependent) reaction networks
          if self.pressureDependence:
              # Recalculate k(T,P) values for modified networks
              self.updateUnimolecularReactionNetworks()
              logging.info('')
              
          # Check new core and edge reactions for Chemkin duplicates
          # The same duplicate reaction gets brought into the core
          # at the same time, so there is no danger in checking all of the edge.
          newCoreReactions = self.core.reactions[numOldCoreReactions:]
          newEdgeReactions = self.edge.reactions[numOldEdgeReactions:]
          checkedReactions = self.core.reactions[:numOldCoreReactions] + self.edge.reactions[:numOldEdgeReactions]
          from rmgpy.chemkin import markDuplicateReaction
          for rxn in newCoreReactions:
              markDuplicateReaction(rxn, checkedReactions)
              checkedReactions.append(rxn)
          if self.saveEdgeSpecies:
              for rxn in newEdgeReactions:
                  markDuplicateReaction(rxn, checkedReactions)
                  checkedReactions.append(rxn)
          self.printEnlargeSummary(
              newCoreSpecies=self.core.species[numOldCoreSpecies:],
              newCoreReactions=self.core.reactions[numOldCoreReactions:],
              reactionsMovedFromEdge=reactionsMovedFromEdge,
              newEdgeSpecies=self.edge.species[numOldEdgeSpecies:],
              newEdgeReactions=self.edge.reactions[numOldEdgeReactions:],
              reactEdge=reactEdge,
          )
  
          logging.info('')
      
      def addNewSurfaceObjects(self,obj,newSurfaceSpecies,newSurfaceReactions,reactionSystem):
          """
          obj is the list of objects for enlargement coming from simulate
          newSurfaceSpecies and newSurfaceReactions are the current lists of surface species and surface reactions
          following simulation
          reactionSystem is the current reactor
          manages surface species and reactions being moved to and from the surface
          moves them to appropriate newSurfaceSpc/RxnsAdd/loss sets
          returns false if the surface has changed
          """
          surfSpcs = set(self.surface.species)
          surfRxns = set(self.surface.reactions)
          
          newSurfaceSpecies = set(newSurfaceSpecies)
          newSurfaceReactions = set(newSurfaceReactions)
                      
          addedRxns = {k for k in obj if isinstance(k,Reaction)}
          addedSurfaceRxns = newSurfaceReactions - surfRxns
                      
          addedBulkRxns = addedRxns-addedSurfaceRxns
          lostSurfaceRxns = (surfRxns - newSurfaceReactions) | addedBulkRxns
                      
          addedSpcs = {k for k in obj if isinstance(k,Species)} | {k.getMaximumLeakSpecies(reactionSystem.T.value_si, reactionSystem.P.value_si) for k in obj if isinstance(k,PDepNetwork)}
          lostSurfaceSpcs = (surfSpcs-newSurfaceSpecies) | addedSpcs
          addedSurfaceSpcs = newSurfaceSpecies - surfSpcs
                      
          self.newSurfaceSpcsAdd = self.newSurfaceSpcsAdd | addedSurfaceSpcs
          self.newSurfaceRxnsAdd = self.newSurfaceRxnsAdd | addedSurfaceRxns
          self.newSurfaceSpcsLoss = self.newSurfaceSpcsLoss | lostSurfaceSpcs
          self.newSurfaceRxnsLoss = self.newSurfaceRxnsLoss | lostSurfaceRxns
          
          return not (self.newSurfaceRxnsAdd != set() or self.newSurfaceRxnsLoss != set() or self.newSurfaceSpcsLoss != set() or self.newSurfaceSpcsAdd != set())
      
      def adjustSurface(self):
          """
          Here we add species intended to be added and remove any species that need to be moved out of the core.  
          For now we remove reactions from the surface that have become part of a PDepNetwork by 
          intersecting the set of surface reactions with the core so that all surface reactions are in the core
          thus the surface algorithm currently (June 2017) is not implemented for pdep networks
          (however it will function fine for non-pdep reactions on a pdep run)
          """
          self.surface.species = list(((set(self.surface.species) | self.newSurfaceSpcsAdd)-self.newSurfaceSpcsLoss) & set(self.core.species))
          self.surface.reactions = list(((set(self.surface.reactions) | self.newSurfaceRxnsAdd)-self.newSurfaceRxnsLoss) & set(self.core.reactions))
          self.clearSurfaceAdjustments()
          
      def clearSurfaceAdjustments(self):
          """
          empties surface tracking varaibles
          """
          self.newSurfaceSpcsAdd = set()
          self.newSurfaceRxnsAdd = set()
          self.newSurfaceSpcsLoss = set()
          self.newSurfaceRxnsLoss = set()
          
      def processNewReactions(self, newReactions, newSpecies, pdepNetwork=None):
          """
          Process a list of newly-generated reactions involving the new core
          species or explored isomer `newSpecies` in network `pdepNetwork`.
          
          Makes a reaction and decides where to put it: core, edge, or PDepNetwork.
          """
          for rxn in newReactions:
              rxn, isNew = self.makeNewReaction(rxn)
              if rxn is None:
                  # Skip this reaction because there was something wrong with it
                  continue
              if isNew:
                  # We've made a new reaction, so make sure the species involved
                  # are in the core or edge
                  allSpeciesInCore = True
                  # Add the reactant and product species to the edge if necessary
                  # At the same time, check if all reactants and products are in the core
                  spcs = []
                  for spec in rxn.reactants:
                      if spec not in self.core.species:
                          allSpeciesInCore = False
                          if spec not in self.edge.species:
                              spcs.append(spec)
                              self.addSpeciesToEdge(spec)
                  for spec in rxn.products:
                      if spec not in self.core.species:
                          allSpeciesInCore = False
                          if spec not in self.edge.species:
                              spcs.append(spec)
                              self.addSpeciesToEdge(spec)
                      
              isomerAtoms = sum([len(spec.molecule[0].atoms) for spec in rxn.reactants])
              
              # Decide whether or not to handle the reaction as a pressure-dependent reaction
              pdep = True
              if not self.pressureDependence:
                  # The pressure dependence option is turned off entirely
                  pdep = False
              elif self.pressureDependence.maximumAtoms is not None and self.pressureDependence.maximumAtoms < isomerAtoms:
                  # The reaction involves so many atoms that pressure-dependent effects are assumed to be negligible
                  pdep = False
              elif not (rxn.isIsomerization() or rxn.isDissociation() or rxn.isAssociation()):
                  # The reaction is not unimolecular in either direction, so it cannot be pressure-dependent
                  pdep = False
              elif rxn.kinetics is not None and rxn.kinetics.isPressureDependent():
                  # The reaction already has pressure-dependent kinetics (e.g. from a reaction library)
                  pdep = False
                  
              # If pressure dependence is on, we only add reactions that are not unimolecular;
              # unimolecular reactions will be added after processing the associated networks
              if not pdep:
                  if not isNew: 
                      # The reaction is not new, so it should already be in the core or edge
                      continue
                  if allSpeciesInCore:
                      self.addReactionToCore(rxn)
                  else:
                      self.addReactionToEdge(rxn)
              else:
                  # Add the reaction to the appropriate unimolecular reaction network
                  # If pdepNetwork is not None then that will be the network the
                  # (path) reactions are added to
                  # Note that this must be done even with reactions that are not new
                  # because of the way partial networks are explored
                  # Since PDepReactions are created as irreversible, not doing so
                  # would cause you to miss the reverse reactions!
                  self.addReactionToUnimolecularNetworks(rxn, newSpecies=newSpecies, network=pdepNetwork)
                  if isinstance(rxn, LibraryReaction):
                      # If reaction came from a reaction library, omit it from the core and edge so that it does 
                      # not get double-counted with the pdep network
                      if rxn in self.core.reactions:
                          self.core.reactions.remove(rxn)
                      if rxn in self.edge.reactions:
                          self.edge.reactions.remove(rxn)
              
              if not numpy.isinf(self.toleranceThermoKeepSpeciesInEdge) and spcs != []: #do thermodynamic filtering
                  self.thermoFilterSpecies(spcs)
                  
      def applyKineticsToReaction(self, reaction):
          """
          retrieve the best kinetics for the reaction and apply it towards the forward 
          or reverse direction (if reverse, flip the direaction).
          """
          from rmgpy.data.rmg import getDB
          # Find the reaction kinetics
          kinetics, source, entry, isForward = self.generateKinetics(reaction)
          # Flip the reaction direction if the kinetics are defined in the reverse direction
          if not isForward:
              family = getDB('kinetics').families[reaction.family]
              reaction.reactants, reaction.products = reaction.products, reaction.reactants
              reaction.pairs = [(p,r) for r,p in reaction.pairs]
              if family.ownReverse and hasattr(reaction,'reverse'):
                  if reaction.reverse:
                      reaction.template = reaction.reverse.template
                      # replace degeneracy
                      reaction.degeneracy = reaction.reverse.degeneracy
                  # We're done with the "reverse" attribute, so delete it to save a bit of memory
                  reaction.reverse = None
          reaction.kinetics = kinetics
  
      def generateKinetics(self, reaction):
          """
          Generate best possible kinetics for the given `reaction` using the kinetics database.
          """
          # Only reactions from families should be missing kinetics
          assert isinstance(reaction, TemplateReaction)
          
          family = getFamilyLibraryObject(reaction.family)
  
          # Get the kinetics for the reaction
          kinetics, source, entry, isForward = family.getKinetics(reaction, templateLabels=reaction.template, degeneracy=reaction.degeneracy, estimator=self.kineticsEstimator, returnAllKinetics=False)
          # Get the gibbs free energy of reaction at 298 K
          G298 = reaction.getFreeEnergyOfReaction(298)
          gibbsIsPositive = G298 > -1e-8
          
          if family.ownReverse and hasattr(reaction,'reverse'):
              if reaction.reverse:
                  # The kinetics family is its own reverse, so we could estimate kinetics in either direction
                  
                  # First get the kinetics for the other direction
                  rev_kinetics, rev_source, rev_entry, rev_isForward = family.getKinetics(reaction.reverse, templateLabels=reaction.reverse.template, degeneracy=reaction.reverse.degeneracy, estimator=self.kineticsEstimator, returnAllKinetics=False)
                  # Now decide which direction's kinetics to keep
                  keepReverse = False
                  if (entry is not None and rev_entry is None):
                      # Only the forward has an entry, meaning an exact match in a depository or template
                      # the reverse must have used an averaged estimated node - so use forward.
                      reason = "This direction matched an entry in {0}, the other was just an estimate.".format(reaction.family)
                  elif (entry is None and rev_entry is not None):
                      # Only the reverse has an entry (see above) - use reverse.
                      keepReverse = True
                      reason = "This direction matched an entry in {0}, the other was just an estimate.".format(reaction.family)
                  elif (entry is not None and rev_entry is not None 
                        and entry is rev_entry):
                      # Both forward and reverse have the same source and entry
                      # Use the one for which the kinetics is the forward kinetics
                      keepReverse = gibbsIsPositive and isForward and rev_isForward
                      reason = "Both directions matched the same entry in {0}, but this direction is exergonic.".format(reaction.family)
                  elif self.kineticsEstimator == 'group additivity' and (kinetics.comment.find("Fitted to 1 rate")>0
                        and not rev_kinetics.comment.find("Fitted to 1 rate")>0) :
                          # forward kinetics were fitted to only 1 rate, but reverse are hopefully better
                          keepReverse = True
                          reason = "Other direction matched a group only fitted to 1 rate."
                  elif self.kineticsEstimator == 'group additivity' and (not kinetics.comment.find("Fitted to 1 rate")>0
                        and rev_kinetics.comment.find("Fitted to 1 rate")>0) :
                          # reverse kinetics were fitted to only 1 rate, but forward are hopefully better
                          keepReverse = False
                          reason = "Other direction matched a group only fitted to 1 rate."
                  elif entry is not None and rev_entry is not None:
                      # Both directions matched explicit rate rules
                      # Keep the direction with the lower (but nonzero) rank
                      if entry.rank < rev_entry.rank and entry.rank != 0:
                          keepReverse = False
                          reason = "Both directions matched explicit rate rules, but this direction has a rule with a lower rank ({0} vs {1}).".format(entry.rank, rev_entry.rank)
                      elif rev_entry.rank < entry.rank and rev_entry.rank != 0:
                          keepReverse = True
                          reason = "Both directions matched explicit rate rules, but this direction has a rule with a lower rank ({0} vs {1}).".format(rev_entry.rank, entry.rank)
                      # Otherwise keep the direction that is exergonic at 298 K
                      else:
                          keepReverse = gibbsIsPositive and isForward and rev_isForward
                          reason = "Both directions matched explicit rate rules, but this direction is exergonic."
                  else:
                      # Keep the direction that is exergonic at 298 K
                      # This must be done after the thermo generation step
                      keepReverse = gibbsIsPositive and isForward and rev_isForward
                      reason = "Both directions are estimates, but this direction is exergonic."
  
                  if keepReverse:
                      kinetics = rev_kinetics
                      source = rev_source
                      entry = rev_entry
                      isForward = not rev_isForward
                      G298 = -G298
                  
                  if self.verboseComments:
                      kinetics.comment += "\nKinetics were estimated in this direction instead of the reverse because:\n{0}".format(reason)
                      kinetics.comment += "\ndGrxn(298 K) = {0:.2f} kJ/mol".format( G298 / 1000.)
              
          # The comments generated by the database for estimated kinetics can
          # be quite long, and therefore not very useful
          # We don't want to waste lots of memory storing these long, 
          # uninformative strings, so here we replace them with much shorter ones
          if not self.verboseComments:
              # Only keep a short comment (to save memory)
              if 'Exact' in kinetics.comment:
                  # Exact match of rate rule
                  pass
              elif 'Matched reaction' in kinetics.comment:
                  # Stems from matching a reaction from a depository
                  pass
              else:
                  # Estimated (averaged) rate rule
                  kinetics.comment =  kinetics.comment[kinetics.comment.find('Estimated'):]
                  
          return kinetics, source, entry, isForward
      
      def printEnlargeSummary(self, newCoreSpecies, newCoreReactions, newEdgeSpecies, newEdgeReactions, reactionsMovedFromEdge=None, reactEdge=False):
          """
          Output a summary of a model enlargement step to the log. The details of
          the enlargement are passed in the `newCoreSpecies`, `newCoreReactions`,
          `newEdgeSpecies`, and `newEdgeReactions` objects. 
          """
  
          logging.info('')
          if reactEdge:
              logging.info('Summary of Secondary Model Edge Enlargement')
          else:
              logging.info('Summary of Model Enlargement')
          logging.info('---------------------------------')
  
          logging.info('Added {0:d} new core species'.format(len(newCoreSpecies)))
          for spec in newCoreSpecies:
              display(spec)
              logging.info('    {0}'.format(spec))
  
          logging.info('Created {0:d} new edge species'.format(len(newEdgeSpecies)))
          for spec in newEdgeSpecies:
              display(spec)
              logging.info('    {0}'.format(spec))
          
          if reactionsMovedFromEdge:
              logging.info('Moved {0:d} reactions from edge to core'.format(len(reactionsMovedFromEdge)))
              for rxn in reactionsMovedFromEdge:
                  for r in newCoreReactions:
                      if ((r.reactants == rxn.reactants and r.products == rxn.products) or
                          (r.products == rxn.reactants and r.reactants == rxn.products)):
                          logging.info('    {0}'.format(r))
                          newCoreReactions.remove(r)
                          break        
          logging.info('Added {0:d} new core reactions'.format(len(newCoreReactions)))
          for rxn in newCoreReactions:
              logging.info('    {0}'.format(rxn))
              
          logging.info('Created {0:d} new edge reactions'.format(len(newEdgeReactions)))
          for rxn in newEdgeReactions:
              logging.info('    {0}'.format(rxn))
  
          coreSpeciesCount, coreReactionCount, edgeSpeciesCount, edgeReactionCount = self.getModelSize()
  
          # Output current model size information after enlargement
          logging.info('')
          logging.info('After model enlargement:')
          logging.info('    The model core has {0:d} species and {1:d} reactions'.format(coreSpeciesCount, coreReactionCount))
          logging.info('    The model edge has {0:d} species and {1:d} reactions'.format(edgeSpeciesCount, edgeReactionCount))
          logging.info('')
  
      def addSpeciesToCore(self, spec):
          """
          Add a species `spec` to the reaction model core (and remove from edge if
          necessary). This function also moves any reactions in the edge that gain
          core status as a result of this change in status to the core.
          If this are any such reactions, they are returned in a list.
          """
  
          assert spec not in self.core.species, "Tried to add species {0} to core, but it's already there".format(spec.label)
  
          # Add the species to the core
          self.core.species.append(spec)
          
          rxnList = []
          if spec in self.edge.species:
  
              # If species was in edge, remove it
              logging.debug("Removing species {0} from edge.".format(spec))
              self.edge.species.remove(spec)
  
              # Search edge for reactions that now contain only core species;
              # these belong in the model core and will be moved there
              for rxn in self.edge.reactions:
                  allCore = True
                  for reactant in rxn.reactants:
                      if reactant not in self.core.species: allCore = False
                  for product in rxn.products:
                      if product not in self.core.species: allCore = False
                  if allCore: rxnList.append(rxn)
  
              # Move any identified reactions to the core
              for rxn in rxnList:
                  self.addReactionToCore(rxn)
                  logging.debug("Moving reaction from edge to core: {0}".format(rxn))
          return rxnList
  
      def addSpeciesToEdge(self, spec):
          """
          Add a species `spec` to the reaction model edge.
          """
          self.edge.species.append(spec)
      
      def setThermodynamicFilteringParameters(self,Tmax, toleranceThermoKeepSpeciesInEdge,minCoreSizeForPrune,maximumEdgeSpecies,reactionSystems):
          """
          sets parameters for thermodynamic filtering based on the current core
          Tmax is the maximum reactor temperature in K
          toleranceThermoKeepSpeciesInEdge is the Gibbs number above which species will be filtered
          minCoreSizeForPrune is the core size at which thermodynamic filtering will start
          maximumEdgeSpecies is the maximum allowed number of edge species
          reactionSystems is a list of reactionSystem objects
          """
          self.Tmax = Tmax
          Gs = [spc.thermo.getFreeEnergy(Tmax) for spc in self.core.species]
          self.Gmax = max(Gs)
          self.Gmin = min(Gs)
          
          self.Gfmax = toleranceThermoKeepSpeciesInEdge*(self.Gmax-self.Gmin)+self.Gmax
          self.toleranceThermoKeepSpeciesInEdge = toleranceThermoKeepSpeciesInEdge
          self.minCoreSizeForPrune = minCoreSizeForPrune
          self.reactionSystems = reactionSystems
          self.maximumEdgeSpecies = maximumEdgeSpecies
      
      def thermoFilterSpecies(self, spcs):
          """
          checks Gibbs energy of the species in species against the
          maximum allowed Gibbs energy
          """
          Tmax = self.Tmax
          for spc in spcs:
              G = spc.thermo.getFreeEnergy(Tmax)
              if G > self.Gfmax:
                  Gn = (G-self.Gmax)/(self.Gmax-self.Gmin)
                  logging.info('Removing species {0} with Gibbs energy {1} from edge because it\'s Gibbs number {2} is greater than the toleranceThermoKeepSpeciesInEdge of {3} '.format(spc,G,Gn,self.toleranceThermoKeepSpeciesInEdge))
                  self.removeSpeciesFromEdge(self.reactionSystems,spc)
                  
          # Delete any networks that became empty as a result of pruning
          if self.pressureDependence:
              self.removeEmptyPdepNetworks()
                          
      def thermoFilterDown(self,maximumEdgeSpecies,minSpeciesExistIterationsForPrune=0):
          """
          removes species from the edge based on their Gibbs energy until maximumEdgeSpecies
          is reached under the constraint that all removed species are older than
          minSpeciesExistIterationsForPrune iterations
          maximumEdgeSpecies is the maximum allowed number of edge species
          minSpeciesExistIterationsForPrune is the number of iterations a species must be in the edge
          before it is eligible for thermo filtering
          """
          Tmax = self.Tmax
          numToRemove = len(self.edge.species) - maximumEdgeSpecies
          logging.debug('Planning to remove {0} species'.format(numToRemove))
          iteration = self.iterationNum
              
          if numToRemove > 0: #implies flux pruning is off or did not trigger
              logging.info('Reached maximum number of edge species')
              logging.info('Attempting to remove excess edge species with Thermodynamic filtering')
              spcs = self.edge.species
              Gfs = numpy.array([spc.thermo.getFreeEnergy(Tmax) for spc in spcs])
              Gns = (Gfs-self.Gmax)/(self.Gmax-self.Gmin) 
              inds = numpy.argsort(Gns) #could actually do this with the Gfs, but want to print the Gn value later
              inds = inds[::-1] #get in order of increasing Gf
  
              ind = 0
              removeSpcs = []
              
              
              while ind < len(inds) and numToRemove > 0: #find the species we can remove and collect indices for removal     
                  i = inds[ind]
                  spc = spcs[i]
                  if iteration - spc.creationIteration >= minSpeciesExistIterationsForPrune:
                      removeSpcs.append(spc)
                      numToRemove -= 1
                  ind += 1
              
              logging.debug('found {0} eligible species for filtering'.format(len(removeSpcs)))
              
              for spc in removeSpcs:
                  logging.info('Removing species {0} from edge to meet maximum number of edge species, Gibbs number is {1}'.format(spc,Gns[i]))
                  self.removeSpeciesFromEdge(self.reactionSystems,spc)
              
              # Delete any networks that became empty as a result of pruning
              if self.pressureDependence:
                  self.removeEmptyPdepNetworks()
              
              #call garbage collection
              collected = gc.collect()
              logging.info('Garbage collector: collected %d objects.' % (collected))
              
      def removeEmptyPdepNetworks(self):
          """
          searches for and deletes any empty pdep networks
          """
          networksToDelete = []
          for network in self.networkList:
              if len(network.pathReactions) == 0 and len(network.netReactions) == 0:
                  networksToDelete.append(network)
                      
          if len(networksToDelete) > 0:
              logging.info('Deleting {0:d} empty pressure-dependent reaction networks'.format(len(networksToDelete)))
              for network in networksToDelete:
                  logging.debug('    Deleting empty pressure dependent reaction network #{0:d}'.format(network.index))
                  source = tuple(network.source)
                  nets_with_this_source = self.networkDict[source]
                  nets_with_this_source.remove(network)
                  if not nets_with_this_source:
                      del(self.networkDict[source])
                  self.networkList.remove(network)
                      
      def prune(self, reactionSystems, toleranceKeepInEdge, maximumEdgeSpecies, minSpeciesExistIterationsForPrune):
          """
          Remove species from the model edge based on the simulation results from
          the list of `reactionSystems`.
          """
  
          ineligibleSpecies = []     # A list of the species which are not eligible for pruning, for any reason
          prunableSpecies = reactionSystems[0].prunableSpecies
          prunableNetworks = reactionSystems[0].prunableNetworks
          
          numPrunableSpecies = len(prunableSpecies)
          iteration = self.iterationNum
          # All edge species that have not existed for more than two enlarge
          # iterations are ineligible for pruning
          for spec in prunableSpecies:
              if iteration - spec.creationIteration <= minSpeciesExistIterationsForPrune:
                  ineligibleSpecies.append(spec)
  
          # Get the maximum species rates (and network leak rates)
          # across all reaction systems
          maxEdgeSpeciesRateRatios = numpy.zeros((numPrunableSpecies), numpy.float64)
          for reactionSystem in reactionSystems:
              for i in range(numPrunableSpecies):
                  rateRatio = reactionSystem.maxEdgeSpeciesRateRatios[i]
                  if maxEdgeSpeciesRateRatios[i] < rateRatio:
                      maxEdgeSpeciesRateRatios[i] = rateRatio
  
              for i in range(len(prunableNetworks)):
                  network = prunableNetworks[i]
                  rateRatio = reactionSystem.maxNetworkLeakRateRatios[i]
                  # Add the fraction of the network leak rate contributed by
                  # each unexplored species to that species' rate
                  # This is to ensure we have an overestimate of that species flux
                  ratios = network.getLeakBranchingRatios(reactionSystem.T.value_si,reactionSystem.P.value_si)
                  for spec, frac in ratios.iteritems():
                      if spec in prunableSpecies:
                          index = prunableSpecies.index(spec)
                          maxEdgeSpeciesRateRatios[index] += frac * rateRatio
                  # Mark any species that is explored in any partial network as ineligible for pruning
                  for spec in network.explored:
                      if spec not in ineligibleSpecies:
                          ineligibleSpecies.append(spec)
  
          # Sort the edge species rates by index
          indices = numpy.argsort(maxEdgeSpeciesRateRatios)
          # Determine which species to prune
          speciesToPrune = []
          pruneDueToRateCounter = 0
          for index in indices:
              spec = prunableSpecies[index]
              if spec in ineligibleSpecies or not spec in self.edge.species:
                  continue
              # Remove the species with rates below the pruning tolerance from the model edge
              if maxEdgeSpeciesRateRatios[index] < toleranceKeepInEdge:
                  speciesToPrune.append((index, spec))
                  pruneDueToRateCounter += 1
              # Keep removing species with the lowest rates until we are below the maximum edge species size
              elif numPrunableSpecies - len(speciesToPrune) > maximumEdgeSpecies:
                  logging.info('Pruning species {0} to make numEdgeSpecies smaller than maximumEdgeSpecies'.format(spec)) # repeated ~15 lines below
                  speciesToPrune.append((index, spec))
              else:
                  break
  
          # Actually do the pruning
          if pruneDueToRateCounter > 0:
              logging.info('Pruning {0:d} species whose rate ratios against characteristic rate did not exceed the minimum threshold of {1:g}'.format(pruneDueToRateCounter, toleranceKeepInEdge))
              for index, spec in speciesToPrune[0:pruneDueToRateCounter]:
                  logging.info('Pruning species {0:<56}'.format(spec))
                  logging.debug('    {0:<56}    {1:10.4e}'.format(spec, maxEdgeSpeciesRateRatios[index]))
                  self.removeSpeciesFromEdge(reactionSystems, spec)
          if len(speciesToPrune) - pruneDueToRateCounter > 0:
              logging.info('Pruning {0:d} species to obtain an edge size of {1:d} species'.format(len(speciesToPrune) - pruneDueToRateCounter, maximumEdgeSpecies))
              for index, spec in speciesToPrune[pruneDueToRateCounter:]:
                  logging.info('Pruning species {0:<56}'.format(spec))
                  logging.debug('    {0:<56}    {1:10.4e}'.format(spec, maxEdgeSpeciesRateRatios[index]))
                  self.removeSpeciesFromEdge(reactionSystems, spec)
  
          # Delete any networks that became empty as a result of pruning
          if self.pressureDependence:
              self.removeEmptyPdepNetworks()
  
          logging.info('')
  
  
      def removeSpeciesFromEdge(self, reactionSystems, spec):
          """
          Remove species `spec` from the reaction model edge.
          """
  
          # remove the species
          self.edge.species.remove(spec)
          self.indexSpeciesDict.pop(spec.index)
  
          # clean up species references in reactionSystems
          for reactionSystem in reactionSystems:
              try:
                  reactionSystem.speciesIndex.pop(spec)
              except KeyError:
                  pass
  
              # identify any reactions it's involved in
              rxnList = []
              for rxn in reactionSystem.reactionIndex:
                  if spec in rxn.reactants or spec in rxn.products:
                      rxnList.append(rxn)
  
              for rxn in rxnList:
                  reactionSystem.reactionIndex.pop(rxn)
  
          # identify any reactions it's involved in
          rxnList = []
          for rxn in self.edge.reactions:
              if spec in rxn.reactants or spec in rxn.products:
                  rxnList.append(rxn)
          # remove those reactions
          for rxn in rxnList:
              self.edge.reactions.remove(rxn)
          
          # Remove the species from any unirxn networks it is in
          if self.pressureDependence:
              for network in self.networkList:
                  # Delete all path reactions involving the species
                  rxnList = []
                  for rxn in network.pathReactions:
                      if spec in rxn.reactants or spec in rxn.products:
                          rxnList.append(rxn)
                  if len(rxnList) > 0:
                      for rxn in rxnList:
                          network.pathReactions.remove(rxn)
                      # Delete all net reactions involving the species
                      rxnList = []
                      for rxn in network.netReactions:
                          if spec in rxn.reactants or spec in rxn.products:
                              rxnList.append(rxn)
                      for rxn in rxnList:
                          network.netReactions.remove(rxn)
                          
                      # Recompute the isomers, reactants, and products for this network
                      network.updateConfigurations(self)
  
          # Remove from the global list of reactions
          # also remove it from the global list of reactions
          for family in self.reactionDict:
              if spec in self.reactionDict[family]:
                  del self.reactionDict[family][spec]
              for reactant1 in self.reactionDict[family]:
                  if spec in self.reactionDict[family][reactant1]:
                      del self.reactionDict[family][reactant1][spec]
              for reactant1 in self.reactionDict[family]:
                  for reactant2 in self.reactionDict[family][reactant1]:
                      tempRxnDeleteList = []
                      for templateReaction in self.reactionDict[family][reactant1][reactant2]:
                          if spec in templateReaction.reactants or spec in templateReaction.products:
                              tempRxnDeleteList.append(templateReaction)
                      for tempRxnToBeDeleted in tempRxnDeleteList:
                          self.reactionDict[family][reactant1][reactant2].remove(tempRxnToBeDeleted)
  
          # remove from the global list of species, to free memory
          formula = spec.molecule[0].getFormula()
          self.speciesDict[formula].remove(spec)
          if spec in self.speciesCache:
              self.speciesCache.remove(spec)
              self.speciesCache.append(None)
  
      def addReactionToCore(self, rxn):
          """
          Add a reaction `rxn` to the reaction model core (and remove from edge if
          necessary). This function assumes `rxn` has already been checked to
          ensure it is supposed to be a core reaction (i.e. all of its reactants
          AND all of its products are in the list of core species).
          """
          if rxn not in self.core.reactions:
              self.core.reactions.append(rxn)
          if rxn in self.edge.reactions:
              self.edge.reactions.remove(rxn)
          
      def addReactionToEdge(self, rxn):
          """
          Add a reaction `rxn` to the reaction model edge. This function assumes
          `rxn` has already been checked to ensure it is supposed to be an edge
          reaction (i.e. all of its reactants OR all of its products are in the
          list of core species, and the others are in either the core or the
          edge).
          """
          self.edge.reactions.append(rxn)
  
      def getModelSize(self):
          """
          Return the numbers of species and reactions in the model core and edge.
          Note that this is not necessarily equal to the lengths of the
          corresponding species and reaction lists.
          """
          coreSpeciesCount = len(self.core.species)
          coreReactionsCount = len(self.core.reactions)
          edgeSpeciesCount = len(self.edge.species)
          edgeReactionsCount = len(self.edge.reactions)
          return (coreSpeciesCount, coreReactionsCount, edgeSpeciesCount, edgeReactionsCount)
  
      def getLists(self):
          """
          Return lists of all of the species and reactions in the core and the
          edge.
          """
          speciesList = []
          speciesList.extend(self.core.species)
          speciesList.extend(self.edge.species)
          reactionList = []
          reactionList.extend(self.core.reactions)
          reactionList.extend(self.edge.reactions)
          return speciesList, reactionList
  
      def getStoichiometryMatrix(self):
          """
          Return the stoichiometry matrix for all generated species and reactions.
          The id of each species and reaction is the corresponding row and column,
          respectively, in the matrix.
          """
          speciesList, reactionList = self.getLists()
          from scipy import sparse
          stoichiometry = sparse.dok_matrix((self.speciesCounter, self.reactionCounter), float)
          for rxn in reactionList:
              j = rxn.index - 1
              specList = rxn.reactants[:]; specList.extend(rxn.products)
              for spec in specList:
                  i = spec.index - 1
                  nu = rxn.getStoichiometricCoefficient(spec)
                  if nu != 0: stoichiometry[i,j] = nu
          return stoichiometry.tocsr()
  
      def addSeedMechanismToCore(self, seedMechanism, react=False):
          """
          Add all species and reactions from `seedMechanism`, a 
          :class:`KineticsPrimaryDatabase` object, to the model core. If `react`
          is ``True``, then reactions will also be generated between the seed
          species. For large seed mechanisms this can be prohibitively expensive,
          so it is not done by default.
          """
  
          if react: raise NotImplementedError("react=True doesn't work yet")
          database = rmgpy.data.rmg.database
  
          from rmgpy.rmg.input import rmg
          
          self.newReactionList = []; self.newSpeciesList = []
  
          numOldCoreSpecies = len(self.core.species)
          numOldCoreReactions = len(self.core.reactions)
  
          logging.info('Adding seed mechanism {0} to model core...'.format(seedMechanism))
  
          seedMechanism = database.kinetics.libraries[seedMechanism]
          
          rxns = seedMechanism.getLibraryReactions()
          for rxn in rxns:
              r, isNew = self.makeNewReaction(rxn) # updates self.newSpeciesList and self.newReactionlist
              if not isNew: logging.info("This library reaction was not new: {0}".format(rxn))
              
          # Perform species constraints and forbidden species checks
          
          for spec in self.newSpeciesList:
              if database.forbiddenStructures.isMoleculeForbidden(spec.molecule[0]):
                  if 'allowed' in rmg.speciesConstraints and 'seed mechanisms' in rmg.speciesConstraints['allowed']:
                      logging.warning("Species {0} from seed mechanism {1} is globally forbidden.  It will behave as an inert unless found in a seed mechanism or reaction library.".format(spec.label, seedMechanism.label))
                  else:
                      raise ForbiddenStructureException("Species {0} from seed mechanism {1} is globally forbidden. You may explicitly allow it, but it will remain inert unless found in a seed mechanism or reaction library.".format(spec.label, seedMechanism.label))
              if failsSpeciesConstraints(spec):
                  if 'allowed' in rmg.speciesConstraints and 'seed mechanisms' in rmg.speciesConstraints['allowed']:
                      rmg.speciesConstraints['explicitlyAllowedMolecules'].extend(spec.molecule)
                  else:
                      raise ForbiddenStructureException("Species constraints forbids species {0} from seed mechanism {1}. Please reformulate constraints, remove the species, or explicitly allow it.".format(spec.label, seedMechanism.label))
  
          for spec in self.newSpeciesList:            
              if spec.reactive:
                  submit(spec,self.solventName)
  
              self.addSpeciesToCore(spec)
  
          for rxn in self.newReactionList:
              if self.pressureDependence and rxn.isUnimolecular():
                  # If this is going to be run through pressure dependence code,
                  # we need to make sure the barrier is positive.
                  # ...but are Seed Mechanisms run through PDep? Perhaps not.
                  for spec in itertools.chain(rxn.reactants, rxn.products):
                      submit(spec,self.solventName)
  
                  rxn.fixBarrierHeight(forcePositive=True)
              self.addReactionToCore(rxn)
          
          # Check we didn't introduce unmarked duplicates
          self.markChemkinDuplicates()
          
          self.printEnlargeSummary(
              newCoreSpecies=self.core.species[numOldCoreSpecies:],
              newCoreReactions=self.core.reactions[numOldCoreReactions:],
              newEdgeSpecies=[],
              newEdgeReactions=[],
          )
  
  
  
      def addReactionLibraryToEdge(self, reactionLibrary):
          """
          Add all species and reactions from `reactionLibrary`, a
          :class:`KineticsPrimaryDatabase` object, to the model edge.
          """
  
          database = rmgpy.data.rmg.database
  
          from rmgpy.rmg.input import rmg
  
          self.newReactionList = []
          self.newSpeciesList = []
  
          numOldEdgeSpecies = len(self.edge.species)
          numOldEdgeReactions = len(self.edge.reactions)
  
          logging.info('Adding reaction library {0} to model edge...'.format(reactionLibrary))
          reactionLibrary = database.kinetics.libraries[reactionLibrary]
  
          rxns = reactionLibrary.getLibraryReactions()
          for rxn in rxns:
              r, isNew = self.makeNewReaction(rxn) # updates self.newSpeciesList and self.newReactionlist
              if not isNew: logging.info("This library reaction was not new: {0}".format(rxn))
  
          # Perform species constraints and forbidden species checks
          for spec in self.newSpeciesList:
              if database.forbiddenStructures.isMoleculeForbidden(spec.molecule[0]):
                  if 'allowed' in rmg.speciesConstraints and 'reaction libraries' in rmg.speciesConstraints['allowed']:
                      logging.warning("Species {0} from reaction library {1} is globally forbidden.  It will behave as an inert unless found in a seed mechanism or reaction library.".format(spec.label, reactionLibrary.label))
                  else:
                      raise ForbiddenStructureException("Species {0} from reaction library {1} is globally forbidden. You may explicitly allow it, but it will remain inert unless found in a seed mechanism or reaction library.".format(spec.label, reactionLibrary.label))
              if failsSpeciesConstraints(spec):
                  if 'allowed' in rmg.speciesConstraints and 'reaction libraries' in rmg.speciesConstraints['allowed']:
                      rmg.speciesConstraints['explicitlyAllowedMolecules'].extend(spec.molecule)
                  else:
                      raise ForbiddenStructureException("Species constraints forbids species {0} from reaction library {1}. Please reformulate constraints, remove the species, or explicitly allow it.".format(spec.label, reactionLibrary.label))
  
          for spec in self.newSpeciesList:
              if spec.reactive: 
                  submit(spec,self.solventName)
  
              self.addSpeciesToEdge(spec)
  
          for rxn in self.newReactionList:
              # Note that we haven't actually evaluated any fluxes at this point
              # Instead, we remove the comment below if the reaction is moved to
              # the core later in the mechanism generation
              self.addReactionToEdge(rxn)
  
          self.printEnlargeSummary(
              newCoreSpecies=[],
              newCoreReactions=[],
              newEdgeSpecies=self.edge.species[numOldEdgeSpecies:],
              newEdgeReactions=self.edge.reactions[numOldEdgeReactions:],
          )
  
      def addReactionLibraryToOutput(self, reactionLib):
          """
          Add all species and reactions from `reactionLibrary`, a
          :class:`KineticsPrimaryDatabase` object, to the output.
          This does not bring any of the reactions or species into the core itself.
          """
  
          logging.info('Adding reaction library {0} to output file...'.format(reactionLib))
          
          # Append the edge reactions that are from the selected reaction library to an output species and output reactions list
          for rxn in self.edge.reactions:
              if isinstance(rxn, LibraryReaction):
                  if rxn.library == reactionLib:
                      self.outputReactionList.append(rxn)
                      
                      for species in rxn.reactants + rxn.products:
                          if species not in self.core.species and species not in self.outputSpeciesList:
                              self.outputSpeciesList.append(species)
                              
  
  
      def addReactionToUnimolecularNetworks(self, newReaction, newSpecies, network=None):
          """
          Given a newly-created :class:`Reaction` object `newReaction`, update the
          corresponding unimolecular reaction network. If no network exists, a new
          one is created. If the new reaction is an isomerization that connects two
          existing networks, the two networks are merged. This function is called
          whenever a new high-pressure limit edge reaction is created. Returns the
          network containing the new reaction.
          """
  
          assert isinstance(newSpecies, Species)
  
          # Put the reaction in the direction in which the new species is in the reactants
          if newSpecies in newReaction.reactants:
              reactants = newReaction.reactants[:]
              products = newReaction.products[:]
          else:
              reactants = newReaction.products[:]
              products = newReaction.reactants[:] 
          reactants.sort()
          products.sort()
          
          source = tuple(reactants)
  
          # Only search for a network if we don't specify it as a parameter
          if network is None:
              if len(reactants) == 1:
                  # Find the network containing the reactant as the source
                  try:
                      networks = self.networkDict[source]
                      assert len(networks) == 1
                      network = networks[0]
                  except KeyError:
                      pass
              elif len(reactants) > 1:
                  # Find the network containing the reactants as the source AND the
                  # product as an explored isomer
                  try:
                      networks = self.networkDict[source]
                      for n in networks:
                          if products[0] in n.explored:
                              assert network is None
                              network = n
                  except KeyError:
                      pass
              else:
                  return
  
              # If no suitable network exists, create a new one
              if network is None:
                  self.networkCount += 1
                  network = PDepNetwork(index=self.networkCount, source=reactants[:])
                          # should the source passed to PDepNetwork constuctor be a tuple not a list? that's what is used in networkDict
                  try:
                      self.networkDict[source].append(network)
                  except KeyError:
                      self.networkDict[source] = [network]
                  self.networkList.append(network)
  
          # Add the path reaction to that network
          network.addPathReaction(newReaction)
  
      def updateUnimolecularReactionNetworks(self):
          """
          Iterate through all of the currently-existing unimolecular reaction
          networks, updating those that have been marked as invalid. In each update,
          the phenomonological rate coefficients :math:`k(T,P)` are computed for
          each net reaction in the network, and the resulting reactions added or
          updated.
          """
  
          # Merge networks if necessary
          # Two partial networks having the same source and containing one or
          # more explored isomers in common must be merged together to avoid
          # double-counting of rates
          for networks in self.networkDict.itervalues():
              networkCount = len(networks)
              for index0, network0 in enumerate(networks):
                  index = index0 + 1
                  while index < networkCount:
                      found = False
                      network = networks[index]
                      if network0.source == network.source:
                          # The networks contain the same source, but do they contain any common included isomers (other than the source)?
                          for isomer in network0.explored:
                              if isomer != network.source and isomer in network.explored:
                                  # The networks contain an included isomer in common, so we need to merge them
                                  found = True
                                  break
                      if found:
                          # The networks contain the same source and one or more common included isomers
                          # Therefore they need to be merged together
                          logging.info('Merging PDepNetwork #{0:d} and PDepNetwork #{1:d}'.format(network0.index, network.index))
                          network0.merge(network)
                          networks.remove(network)
                          self.networkList.remove(network)
                          networkCount -= 1
                      else:
                          index += 1
  
          count = sum([1 for network in self.networkList if not network.valid and not (len(network.explored) == 0 and len(network.source) > 1)])
          logging.info('Updating {0:d} modified unimolecular reaction networks (out of {1:d})...'.format(count, len(self.networkList)))
          
          # Iterate over all the networks, updating the invalid ones as necessary
          # self = reactionModel object
          updatedNetworks = []
          for network in self.networkList:
              if not network.valid:
                  network.update(self, self.pressureDependence)
                  updatedNetworks.append(network)
              
          # PDepReaction objects generated from partial networks are irreversible
          # However, it makes more sense to have reversible reactions in the core
          # Thus we mark PDepReaction objects as reversible and remove the reverse
          # direction from the list of core reactions
          # Note that well-skipping reactions may not have a reverse if the well
          # that they skip over is not itself in the core
          index = 0
          coreReactionCount = len(self.core.reactions)
          while index < coreReactionCount:
              reaction = self.core.reactions[index]
              if isinstance(reaction, PDepReaction):
                  for reaction2 in self.core.reactions[index+1:]:
                      if isinstance(reaction2, PDepReaction) and reaction.reactants == reaction2.products and reaction.products == reaction2.reactants:
                          # We've found the PDepReaction for the reverse direction
                          dGrxn = reaction.getFreeEnergyOfReaction(300.)
                          kf = reaction.getRateCoefficient(1000,1e5)
                          kr = reaction.getRateCoefficient(1000,1e5) / reaction.getEquilibriumConstant(1000)
                          kf2 = reaction2.getRateCoefficient(1000,1e5) / reaction2.getEquilibriumConstant(1000)
                          kr2 = reaction2.getRateCoefficient(1000,1e5)
                          if kf / kf2 < 0.5 or kf / kf2 > 2.0:
                              # Most pairs of reactions should satisfy thermodynamic consistency (or at least be "close")
                              # Warn about the ones that aren't close (but don't abort)
                              logging.warning('Forward and reverse PDepReactions for reaction {0!s} generated from networks {1:d} and {2:d} do not satisfy thermodynamic consistency.'.format(reaction, reaction.network.index, reaction2.network.index))
                              logging.warning('{0!s}:'.format(reaction))
                              logging.warning('{0:.2e} {1:.2e}:'.format(kf, kf2))
                              logging.warning('{0!s}:'.format(reaction2))
                              logging.warning('{0:.2e} {1:.2e}:'.format(kr, kr2))
                          # Keep the exergonic direction
                          keepFirst = dGrxn < 0
                          # Delete the PDepReaction that we aren't keeping
                          if keepFirst:
                              self.core.reactions.remove(reaction2)
                              reaction.reversible = True
                          else:
                              self.core.reactions.remove(reaction)
                              self.core.reactions.remove(reaction2)
                              self.core.reactions.insert(index, reaction2)
                              reaction2.reversible = True
                          coreReactionCount -= 1
                          # There should be only one reverse, so we can stop searching once we've found it
                          break
                  else:
                      reaction.reversible = True
              # Move to the next core reaction
              index += 1
  
  
      def markChemkinDuplicates(self):
          """
          Check that all reactions that will appear the chemkin output have been checked as duplicates.
          
          Call this if you've done something that may have introduced undetected duplicate reactions,
          like add a reaction library or seed mechanism.
          Anything added via the :meth:`expand` method should already be detected.
          """
          from rmgpy.chemkin import markDuplicateReactions
          
          rxnList = self.core.reactions + self.outputReactionList
          markDuplicateReactions(rxnList)
          
      
      def registerReaction(self, rxn):
          """
          Adds the reaction to the reaction database.
  
          The reaction database is structured as a multi-level
          dictionary, for efficient search and retrieval of
          existing reactions.
  
          The database has two types of dictionary keys:
          - reaction family
          - reactant(s) keys
  
          First, the keys are generated for the parameter reaction.
          
          Next, it is checked whether the reaction database already 
          contains similar keys. If not, a new container is created,
          either a dictionary for the family key and first reactant key,
          or a list for the second reactant key.
  
          Finally, the reaction is inserted as the first element in the 
          list.
          """
  
          key_family, key1, key2 = generateReactionKey(rxn)
  
          # make dictionary entries if necessary
          if key_family not in self.reactionDict:
              self.reactionDict[key_family] = {}
  
          if not self.reactionDict[key_family].has_key(key1):
              self.reactionDict[key_family][key1] = {}
  
          if not self.reactionDict[key_family][key1].has_key(key2):
              self.reactionDict[key_family][key1][key2] = []
  
          # store this reaction at the top of the relevant short-list
          self.reactionDict[key_family][key1][key2].insert(0, rxn)
  
  
      def searchRetrieveReactions(self, rxn):
          """
          Searches through the reaction database for 
          reactions with an identical reaction key as the key of the 
          parameter reaction.
  
          Both the reaction key based on the reactants as well as on the products
          is used to search for possible candidate reactions.
          """
  
          # Get the short-list of reactions with the same family, reactant1 and reactant2
          family_label, r1_fwd, r2_fwd = generateReactionKey(rxn)
          
          my_reactionList = []
  
          rxns = self.retrieve(family_label, r1_fwd, r2_fwd)
          my_reactionList.extend(rxns)
              
              
          family = getFamilyLibraryObject(family_label)       
          # if the family is its own reverse (H-Abstraction) then check the other direction
          if isinstance(family,KineticsFamily): 
  
              # Get the short-list of reactions with the same family, product1 and product2
              family_label, r1_rev, r2_rev = generateReactionKey(rxn, useProducts=True)
  
              rxns = self.retrieve(family_label, r1_rev, r2_rev)
              my_reactionList.extend(rxns)
  
          return my_reactionList
  
      def initializeIndexSpeciesDict(self):
          """
          Populates the core species dictionary
          
          integer -> core Species
  
          with the species that are currently in the core.
          """
  
          for spc in itertools.chain(self.core.species, self.edge.species):
              if spc.reactive:
                  self.indexSpeciesDict[spc.index] = spc
  
      def retrieve(self, family_label, key1, key2):
          """
          Returns a list of reactions from the reaction database with the 
          same keys as the parameters.
  
          Returns an empty list when one of the keys could not be found.
          """
          try:
              return self.reactionDict[family_label][key1][key2][:]
          except KeyError: # no such short-list: must be new, unless in seed.
              return []
  
  
  
      def inflate(self, rxn):
          """
          Convert reactions from
          reactants/products that are referring
          to the core species index, to the respective Species objects.
          """
          reactants, products, pairs = [], [], []
  
          for reactant in rxn.reactants:
              reactant = self.getSpecies(reactant)  
              reactants.append(reactant)
  
          for product in rxn.products:
              product = self.getSpecies(product)
              products.append(product)
  
          for reactant, product in rxn.pairs:
              reactant = self.getSpecies(reactant)
              product = self.getSpecies(product)
              pairs.append((reactant, product))
  
          rxn.reactants = reactants  
          rxn.products = products
          rxn.pairs = pairs
  
          return rxn
  
      def getSpecies(self, obj):
          """
          Retrieve species object, by
          polling the index species dictionary.
          """
          if isinstance(obj, int):
              try:
                  spc = self.indexSpeciesDict[obj]
                  return spc
              except KeyError, e:
                  raise e
  
          return obj
  
      def retrieveNewSpecies(self, deflatedRxn):
          """
          Searches for the first reactant or product in the deflated reaction
          that is represented by an integer.
  
          Such an object refers to a core species that was used to generate the
          reaction in the first place. Reactants or products represented by an
          object that is not an integer will be a newly-generated structure.
          """
          for obj in itertools.chain(deflatedRxn.reactants, deflatedRxn.products):
              if isinstance(obj, int):
                  return self.getSpecies(obj)
          raise Exception("No core species were found in either reactants or products of {0}!".format(deflatedRxn))
  
  
  def generateReactionKey(rxn, useProducts=False):
      """
      Returns a tuple with 3 keys:
      - the reaction family (or library) the reaction belongs to
      - the keys of the reactants.
  
      None for the third element in the tuple if there is 
      only 1 reactant.
  
      The keys are sorted alphabetically.
      """
  
      key_family = rxn.family
  
      spc_list = rxn.products if useProducts else rxn.reactants
      key1 = getKey(spc_list[0])
      key2 = None if len(spc_list) == 1 else getKey(spc_list[1])
      key1, key2 = sorted([key1, key2], reverse=True)# ensure None is always at end
  
      return (key_family, key1, key2)
  
  def generateReactionId(rxn):
      """
      Returns a tuple of the reactions reactant and product
      keys.
  
      Both lists are sorted.
  
      The first element in the tuple is the reactants list.
      """
  
  
      reactants = sorted([getKey(reactant) for reactant in rxn.reactants])
      products = sorted([getKey(product) for product in rxn.products])
  
      return (reactants, products)
  
  def getFamilyLibraryObject(label):
      """
      Returns the KineticsFamily or KineticsLibrary object associated with the
      parameter string.
  
      First search through the reaction families, then 
      through the libraries.
      """
  
      kinetics = rmgpy.data.rmg.database.kinetics
  
      try:
          fam = kinetics.families[label]
          return fam
      except KeyError, e:
          pass
  
      try:
          lib = kinetics.libraries[label]
          return lib
      except KeyError, e:
          pass
  
      raise Exception('Could not retrieve the family/library: {}'.format(label))
  
  def getKey(spc):
      """
      Returns a string of the species that can serve as a key in a dictionary.
      """
  
      return spc.label
  
  def areIdenticalSpeciesReferences(rxn1, rxn2):
      """
      Checks if the references of the reactants and products of the two reactions
      are identical, in either direction.
      """
      identical_same_direction = rxn1.reactants == rxn2.reactants and rxn1.products == rxn2.products
      identical_opposite_directions = rxn1.reactants == rxn2.products and rxn1.products == rxn2.reactants
      identical_collider = rxn1.specificCollider == rxn2.specificCollider
      
      return (identical_same_direction or identical_opposite_directions) and identical_collider
  