model.py

#!/usr/bin/env python3

###############################################################################
#                                                                             #
# RMG - Reaction Mechanism Generator                                          #
#                                                                             #
# Copyright (c) 2002-2023 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 gc
import itertools
import logging
import os

import numpy as np

import rmgpy.data.rmg
from rmgpy import settings
from rmgpy.constraints import fails_species_constraints, pass_cutting_threshold
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.data.vaporLiquidMassTransfer import vapor_liquid_mass_transfer
from rmgpy.molecule.group import Group
from rmgpy.data.rmg import get_db
from rmgpy.display import display
from rmgpy.exceptions import ForbiddenStructureException
from rmgpy.kinetics import KineticsData, Arrhenius
from rmgpy.quantity import Quantity
from rmgpy.reaction import Reaction
from rmgpy.rmg.pdep import PDepReaction, PDepNetwork
from rmgpy.rmg.react import react_all
from rmgpy.species import Species
from rmgpy.thermo.thermoengine import submit
from rmgpy.rmg.decay import decay_species
from rmgpy.rmg.reactors import PhaseSystem, Phase, Interface, Reactor
from rmgpy.molecule.fragment import Fragment

################################################################################


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, phases=None, interfaces={}):
        self.species = species or []
        self.reactions = reactions or []
        if phases is None:
            phases = {"Default": Phase(), "Surface": Phase()}
            interfaces = {frozenset({"Default", "Surface"}): Interface(list(phases.values()))}
        self.phase_system = PhaseSystem(phases, interfaces)

    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
        final_model = ReactionModel()

        # Put the current model into the merged model as-is
        final_model.species.extend(self.species)
        final_model.reactions.extend(self.reactions)

        # Determine which species in other are already in self
        common_species = {}
        unique_species = []
        for spec in other.species:
            for spec0 in final_model.species:
                if spec.is_isomorphic(spec0):
                    common_species[spec] = spec0
                    if spec0.label not in ["Ar", "N2", "Ne", "He"]:
                        if not spec0.thermo.is_identical_to(spec.thermo):
                            print("Species {0} thermo from model 1 did not match that of model 2.".format(spec.label))

                    break
            else:
                unique_species.append(spec)

        # Determine which reactions in other are already in self
        common_reactions = {}
        unique_reactions = []
        for rxn in other.reactions:
            for rxn0 in final_model.reactions:
                if rxn.is_isomorphic(rxn0, either_direction=True):
                    common_reactions[rxn] = rxn0
                    if not rxn0.kinetics.is_identical_to(rxn.kinetics):
                        print("Reaction {0} kinetics from model 1 did not match that of model 2.".format(str(rxn0)))
                    break
            else:
                unique_reactions.append(rxn)

        # Add the unique species from other to the final model
        final_model.species.extend(unique_species)

        # Make sure unique reactions only refer to species in the final model
        for rxn in unique_reactions:
            for i, reactant in enumerate(rxn.reactants):
                try:
                    rxn.reactants[i] = common_species[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] = common_species[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
        final_model.reactions.extend(unique_reactions)

        # Return the merged model
        return final_model


################################################################################


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
    `network_dict`             A dictionary of pressure-dependent reaction networks (:class:`Network` objects) indexed by source.
    `network_list`             A list of pressure-dependent reaction networks (:class:`Network` objects)
    `network_count`            A counter for the number of pressure-dependent networks created
    `index_species_dict`       A dictionary with a unique index pointing to the species objects
    `solvent_name`             String describing solvent name for liquid reactions. Empty for non-liquid estimation
    `surface_site_density`     The surface site density (a SurfaceConcentration quantity) or None if no heterogeneous catalyst.
    =========================  ==============================================================


    """

    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.network_dict = {}
        self.network_list = []
        self.network_count = 0
        self.species_dict = {}
        self.reaction_dict = {}
        self.species_cache = [None for i in range(4)]
        self.species_counter = 0
        self.reaction_counter = 0
        self.new_species_list = []
        self.new_reaction_list = []
        self.output_species_list = []
        self.output_reaction_list = []
        self.pressure_dependence = None
        self.quantum_mechanics = None
        self.verbose_comments = False
        self.kinetics_estimator = "rate rules"
        self.index_species_dict = {}
        self.save_edge_species = False
        self.iteration_num = 0
        self.thermo_tol_keep_spc_in_edge = np.inf
        self.Gfmax = np.inf
        self.Gmax = np.inf
        self.Gmin = -np.inf
        self.min_core_size_for_prune = 50
        self.maximum_edge_species = 100000
        self.Tmax = 0
        self.reaction_systems = []
        self.new_surface_spcs_add = set()
        self.new_surface_rxns_add = set()
        self.new_surface_spcs_loss = set()
        self.new_surface_rxns_loss = set()
        self.solvent_name = ""
        self.surface_site_density = None
        self.unrealgroups = [
            Group().from_adjacency_list(
                """
            1  O u0 p2 c0 {2,S} {4,S}
            2  O u0 p2 c0 {1,S} {3,S}
            3  R!H u1 px c0 {2,S}
            4  H u0 p0  c0 {1,S}
            """
            )
        ]

    def check_for_existing_species(self, molecule):
        """
        Check to see if an existing species contains the same
        :class:`molecule.Molecule` as `molecule`. Comparison is done using
        isomorphism without consideration of electrons. Therefore, resonance
        structures of a species will all match each other.

        Returns the matched species if found and `None` otherwise.
        """

        # First check cache and return if species is found
        for i, spec in enumerate(self.species_cache):
            if spec is not None and spec.is_isomorphic(molecule, strict=False):
                self.species_cache.pop(i)
                self.species_cache.insert(0, spec)
                return spec

        # If not found in cache, check all species with matching formula
        formula = molecule.get_formula()
        try:
            species_list = self.species_dict[formula]
        except KeyError:
            pass
        else:
            for spec in species_list:
                if spec.is_isomorphic(molecule, strict=False):
                    self.species_cache.pop()
                    self.species_cache.insert(0, spec)
                    return spec

        # At this point we can conclude that the species is new
        return None

    def make_new_species(self, object, label="", reactive=True, check_existing=True, generate_thermo=True, check_decay=False, check_cut=False):
        """
        Formally create a new species from the specified `object`, which can be
        either a :class:`Molecule` object or an :class:`rmgpy.species.Species`
        object. It is emphasized that `reactive` relates to the :Class:`Species` attribute, while `reactive_structure`
        relates to the :Class:`Molecule` attribute.
        """

        if isinstance(object, rmgpy.species.Species):
            molecule = object.molecule[0]
            label = label if label != "" else object.label
            reactive = object.reactive
        else:
            molecule = object

        molecule.clear_labeled_atoms()

        # If desired, check to ensure that the species is new; return the
        # existing species if not new
        if check_existing:
            spec = self.check_for_existing_species(molecule)
            if spec is not None:
                return spec, False

        # If we're here then we're ready to make the new species
        if check_cut:
            try:
                mols = molecule.cut_molecule(cut_through=False)
            except AttributeError:
                # it's Molecule object, change it to Fragment and then cut
                molecule = Fragment().from_adjacency_list(molecule.to_adjacency_list())
                mols = molecule.cut_molecule(cut_through=False)
            if len(mols) == 1:
                molecule = mols[0]
            else:
                return [self.make_new_species(mol, check_decay=check_decay) for mol in mols]

        try:
            spec = Species(label=label, molecule=[molecule], reactive=reactive, thermo=object.thermo, transport_data=object.transport_data)
        except AttributeError:
            spec = Species(label=label, molecule=[molecule], reactive=reactive)

        spec.generate_resonance_structures()

        if check_decay:
            spcs = decay_species(spec)
            if len(spcs) == 1:
                spec = spcs[0]
            else:
                return [self.make_new_species(spc) for spc in spcs]

        if reactive:
            self.species_counter += 1  # count only reactive species
            spec.index = self.species_counter
        else:
            spec.index = -1

        spec.creation_iteration = self.iteration_num
        spec.molecular_weight = Quantity(spec.molecule[0].get_molecular_weight() * 1000.0, "amu")

        if generate_thermo:
            self.generate_thermo(spec)

        # If the species still does not have a label, set initial label as the SMILES
        # This may change later after getting thermo in self.generate_thermo()
        if not spec.label:
            spec.label = spec.smiles

        # ensure species labels are unique
        orilabel = spec.label
        label = orilabel
        i = 2
        while any([label in phase.names for phase in self.edge.phase_system.phases.values()]):
            label = orilabel + "-" + str(i)
            i += 1
        spec.label = label

        logging.debug("Creating new species %s", spec.label)

        formula = molecule.get_formula()
        if formula in self.species_dict:
            self.species_dict[formula].append(spec)
        else:
            self.species_dict[formula] = [spec]

        # Since the species is new, add it to the list of new species
        self.new_species_list.append(spec)

        if spec.reactive:
            self.index_species_dict[spec.index] = spec

        return spec, True

    def check_for_existing_reaction(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

        family_obj = get_family_library_object(rxn.family)
        shortlist = self.search_retrieve_reactions(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 = generate_reaction_id(rxn)

        for rxn0 in shortlist:
            rxn_id0 = generate_reaction_id(rxn0)

            if rxn_id == rxn_id0 and are_identical_species_references(rxn, rxn0):
                if isinstance(family_obj, KineticsLibrary) or isinstance(family_obj, KineticsFamily):
                    if not rxn.duplicate:
                        return True, rxn0
                else:
                    return True, rxn0
            elif isinstance(family_obj, KineticsFamily) and rxn_id == rxn_id0[::-1] and are_identical_species_references(rxn, rxn0):
                if not rxn.duplicate:
                    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 = generate_reaction_key(rxn)
        _, r1_rev, r2_rev = generate_reaction_key(rxn, useProducts=True)

        for library in self.reaction_dict:
            lib_obj = get_family_library_object(library)
            if isinstance(lib_obj, 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 = generate_reaction_id(rxn0)
                    if (rxn_id == rxn_id0) or (rxn_id == rxn_id0[::-1]):
                        if are_identical_species_references(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 are_identical_species_references(rxn, rxn0):
                        return True, rxn0

        return False, None

    def make_new_reaction(self, forward, check_existing=True, generate_thermo=True, generate_kinetics=True, perform_cut=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.new_reaction_list if it is new.
        """

        # Determine the proper species objects for all reactants and products
        if forward.family and forward.is_forward:
            reactants = [self.make_new_species(reactant, generate_thermo=generate_thermo)[0] for reactant in forward.reactants]
            products = []
            if perform_cut:
                # check if the product is too large so that we can cut
                # maybe species do not contain element C
                try:
                    if len(forward.reactants) > len(forward.products) and pass_cutting_threshold(forward.products[0]):
                        # need to cut
                        check_cut = True
                    else:
                        check_cut = False
                except KeyError:
                    check_cut = False
            else:
                check_cut = False
            for product in forward.products:
                spcs = self.make_new_species(product, generate_thermo=generate_thermo, check_decay=True, check_cut=check_cut)
                if type(spcs) == tuple:
                    products.append(spcs[0])
                elif type(spcs) == list:
                    products.extend([spc[0] for spc in spcs])
            # change the reaction to irreversible if we cut the product into fragments
            if len(products) != 1 and check_cut:
                forward.reversible = False
        else:
            try:
                reactants = [self.make_new_species(reactant, generate_thermo=generate_thermo)[0] for reactant in forward.reactants]
                products = [self.make_new_species(product, generate_thermo=generate_thermo)[0] for product in forward.products]
            except:
                logging.error(f"Error when making species in reaction {forward:s} from {forward.family:s}")
                raise

        if forward.specific_collider is not None:
            forward.specific_collider = self.make_new_species(forward.specific_collider)[0]

        if forward.pairs is not None:
            for pairIndex in range(len(forward.pairs)):
                reactant_index = forward.reactants.index(forward.pairs[pairIndex][0])
                product_index = forward.products.index(forward.pairs[pairIndex][1])
                forward.pairs[pairIndex] = (reactants[reactant_index], products[product_index])
                if hasattr(forward, "reverse"):
                    if forward.reverse:
                        forward.reverse.pairs[pairIndex] = (products[product_index], reactants[reactant_index])
        forward.reactants = reactants
        forward.products = products

        if check_existing:
            found, rxn = self.check_for_existing_reaction(forward)
            if found:
                return rxn, False

        # Generate the reaction pairs if not yet defined
        if forward.pairs is None or len(forward.pairs) != max(len(forward.reactants), len(forward.products)):
            forward.generate_pairs()
            if hasattr(forward, "reverse"):
                if forward.reverse:
                    forward.reverse.generate_pairs()

        # Note in the log
        if isinstance(forward, TemplateReaction):
            logging.debug("Creating new %s template reaction %s", forward.family, forward)
        elif isinstance(forward, DepositoryReaction):
            logging.debug("Creating new %s reaction %s", forward.get_source(), forward)
        elif isinstance(forward, LibraryReaction):
            logging.debug("Creating new library reaction %s", forward)
        else:
            raise Exception("Unrecognized reaction type {0!s}".format(forward.__class__))

        self.register_reaction(forward)

        forward.index = self.reaction_counter + 1
        self.reaction_counter += 1

        if generate_kinetics:
            if forward.kinetics is None:
                self.apply_kinetics_to_reaction(forward)

            if isinstance(forward.kinetics, KineticsData):
                forward.kinetics = forward.kinetics.to_arrhenius()
            #  correct barrier heights of estimated kinetics
            if isinstance(forward, (TemplateReaction, DepositoryReaction)):  # i.e. not LibraryReaction
                forward.fix_barrier_height()  # also converts ArrheniusEP to Arrhenius.

            if self.pressure_dependence and forward.is_unimolecular():
                # If this is going to be run through pressure dependence code,
                # we need to make sure the barrier is positive.
                forward.fix_barrier_height(force_positive=True)

        # Since the reaction is new, add it to the list of new reactions
        self.new_reaction_list.append(forward)

        # Return newly created reaction
        return forward, True

    def make_new_pdep_reaction(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 reaction_counter 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.generate_pairs()

        # Set reaction index and increment the counter
        forward.index = self.reaction_counter + 1
        self.reaction_counter += 1

        return forward

    def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bimolecular_react=None, trimolecular_react=None):
        """
        Enlarge a reaction model by processing the objects in the list `new_object`.
        If `new_object` 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 `new_object`
        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 `react_edge` flag is `True`, then no new_object is needed,
        and instead the algorithm proceeds to react the core species together
        to form edge reactions.
        """

        num_old_core_species = len(self.core.species)
        num_old_core_reactions = len(self.core.reactions)
        num_old_edge_species = len(self.edge.species)
        num_old_edge_reactions = len(self.edge.reactions)
        reactions_moved_from_edge = []
        self.new_reaction_list = []
        self.new_species_list = []

        # Determine number of parallel processes.
        from rmgpy.rmg.main import determine_procnum_from_ram

        procnum = determine_procnum_from_ram()

        if react_edge is False:
            # We are adding core species
            new_reactions = []
            pdep_network = None
            object_was_in_edge = False

            if isinstance(new_object, Species):
                new_species = new_object

                object_was_in_edge = new_species in self.edge.species

                if not new_species.reactive:
                    logging.info("NOT generating reactions for unreactive species {0}".format(new_species))
                else:
                    logging.info("Adding species {0} to model core".format(new_species))
                    display(new_species)  # if running in IPython --pylab mode, draws the picture!

                # Add new species
                if new_species not in self.core.species:
                    reactions_moved_from_edge = self.add_species_to_core(new_species)
                else:
                    reactions_moved_from_edge = []

            elif isinstance(new_object, tuple) and isinstance(new_object[0], PDepNetwork) and self.pressure_dependence:
                pdep_network, new_species = new_object
                new_reactions.extend(pdep_network.explore_isomer(new_species))

                self.process_new_reactions(new_reactions, new_species, pdep_network)

            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(new_object, new_object.__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.network_list:
                network.update_configurations(self)
                index = 0
                isomers = [isomer.species[0] for isomer in network.isomers]
                while index < len(self.core.species):
                    species = self.core.species[index]
                    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:
                            new_reactions = network.explore_isomer(species)

                            self.process_new_reactions(new_reactions, species, network)
                            network.update_configurations(self)
                            index = 0
                            break
                    else:
                        index += 1

            if isinstance(new_object, Species) and object_was_in_edge:
                # moved one species from edge to core
                num_old_edge_species -= 1
                # moved these reactions from edge to core
                num_old_edge_reactions -= len(reactions_moved_from_edge)

        else:
            # Generate reactions between all core species which have not been
            # reacted yet and exceed the reaction filter thresholds
            rxn_lists, spcs_tuples = react_all(
                self.core.species, num_old_core_species, unimolecular_react, bimolecular_react, trimolecular_react=trimolecular_react, procnum=procnum
            )

            for rxnList, spcTuple in zip(rxn_lists, spcs_tuples):
                if rxnList:
                    # Identify a core species which was used to generate the reaction
                    # This is only used to determine the reaction direction for processing
                    spc = spcTuple[0]
                    self.process_new_reactions(rxnList, spc)

        ################################################################
        # Begin processing the new species and reactions

        # Generate thermo for new species
        if self.new_species_list:
            logging.info("Generating thermo for new species...")
            self.apply_thermo_to_species(procnum)

        # Do thermodynamic filtering
        if not np.isinf(self.thermo_tol_keep_spc_in_edge) and self.new_species_list != []:
            self.thermo_filter_species(self.new_species_list)

        # Update unimolecular (pressure dependent) reaction networks
        if self.pressure_dependence:
            # Recalculate k(T,P) values for modified networks
            self.update_unimolecular_reaction_networks()
            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.
        new_core_reactions = self.core.reactions[num_old_core_reactions:]
        new_edge_reactions = self.edge.reactions[num_old_edge_reactions:]
        checked_reactions = self.core.reactions[:num_old_core_reactions] + self.edge.reactions[:num_old_edge_reactions]
        from rmgpy.chemkin import mark_duplicate_reaction

        for rxn in new_core_reactions:
            mark_duplicate_reaction(rxn, checked_reactions)
            checked_reactions.append(rxn)
        if self.save_edge_species:
            for rxn in new_edge_reactions:
                mark_duplicate_reaction(rxn, checked_reactions)
                checked_reactions.append(rxn)
        self.log_enlarge_summary(
            new_core_species=self.core.species[num_old_core_species:],
            new_core_reactions=self.core.reactions[num_old_core_reactions:],
            reactions_moved_from_edge=reactions_moved_from_edge,
            new_edge_species=self.edge.species[num_old_edge_species:],
            new_edge_reactions=self.edge.reactions[num_old_edge_reactions:],
            react_edge=react_edge,
        )

        logging.info("")

    def add_new_surface_objects(self, obj, new_surface_species, new_surface_reactions, reaction_system):
        """
        obj is the list of objects for enlargement coming from simulate
        new_surface_species and new_surface_reactions are the current lists of surface species and surface reactions
        following simulation
        reaction_system 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
        """
        surf_spcs = set(self.surface.species)
        surf_rxns = set(self.surface.reactions)

        new_surface_species = set(new_surface_species)
        new_surface_reactions = set(new_surface_reactions)

        added_rxns = {k for k in obj if isinstance(k, Reaction)}
        added_surface_rxns = new_surface_reactions - surf_rxns

        added_bulk_rxns = added_rxns - added_surface_rxns
        lost_surface_rxns = (surf_rxns - new_surface_reactions) | added_bulk_rxns

        added_spcs = {k for k in obj if isinstance(k, Species)} | {
            k.get_maximum_leak_species(reaction_system.T.value_si, reaction_system.P.value_si) for k in obj if isinstance(k, PDepNetwork)
        }
        lost_surface_spcs = (surf_spcs - new_surface_species) | added_spcs
        added_surface_spcs = new_surface_species - surf_spcs

        self.new_surface_spcs_add = self.new_surface_spcs_add | added_surface_spcs
        self.new_surface_rxns_add = self.new_surface_rxns_add | added_surface_rxns
        self.new_surface_spcs_loss = self.new_surface_spcs_loss | lost_surface_spcs
        self.new_surface_rxns_loss = self.new_surface_rxns_loss | lost_surface_rxns

        return not (
            self.new_surface_rxns_add != set()
            or self.new_surface_rxns_loss != set()
            or self.new_surface_spcs_loss != set()
            or self.new_surface_spcs_add != set()
        )

    def adjust_surface(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.new_surface_spcs_add) - self.new_surface_spcs_loss) & set(self.core.species))
        self.surface.reactions = list(
            ((set(self.surface.reactions) | self.new_surface_rxns_add) - self.new_surface_rxns_loss) & set(self.core.reactions)
        )
        self.clear_surface_adjustments()

    def clear_surface_adjustments(self):
        """
        empties surface tracking varaibles
        """
        self.new_surface_spcs_add = set()
        self.new_surface_rxns_add = set()
        self.new_surface_spcs_loss = set()
        self.new_surface_rxns_loss = set()

    def process_new_reactions(self, new_reactions, new_species, pdep_network=None, generate_thermo=True, generate_kinetics=True):
        """
        Process a list of newly-generated reactions involving the new core
        species or explored isomer `new_species` in network `pdep_network`.

        Makes a reaction and decides where to put it: core, edge, or PDepNetwork.
        """
        for rxn in new_reactions:
            rxn, is_new = self.make_new_reaction(rxn, generate_thermo=generate_thermo, generate_kinetics=generate_kinetics)
            if rxn is None:
                # Skip this reaction because there was something wrong with it
                continue
            if is_new:
                # We've made a new reaction, so make sure the species involved
                # are in the core or edge
                all_species_in_core = 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
                for spec in rxn.reactants:
                    if spec not in self.core.species:
                        all_species_in_core = False
                        if spec not in self.edge.species:
                            self.add_species_to_edge(spec)
                for spec in rxn.products:
                    if spec not in self.core.species:
                        all_species_in_core = False
                        if spec not in self.edge.species:
                            self.add_species_to_edge(spec)

            isomer_atoms = 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.pressure_dependence:
                # The pressure dependence option is turned off entirely
                pdep = False
            elif self.pressure_dependence.maximum_atoms is not None and self.pressure_dependence.maximum_atoms < isomer_atoms:
                # The reaction involves so many atoms that pressure-dependent effects are assumed to be negligible
                pdep = False
            elif not (rxn.is_isomerization() or rxn.is_dissociation() or rxn.is_association()):
                # The reaction is not unimolecular in either direction, so it cannot be pressure-dependent
                pdep = False
            elif isinstance(rxn, LibraryReaction):
                # Try generating the high pressure limit kinetics. If successful, set pdep to ``True``, and vice versa.
                pdep = rxn.generate_high_p_limit_kinetics()
            elif any([any([x.is_subgraph_isomorphic(q) for q in self.unrealgroups]) for y in rxn.reactants + rxn.products for x in y.molecule]):
                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 is_new:
                    # The reaction is not new, so it should already be in the core or edge
                    continue
                if all_species_in_core:
                    self.add_reaction_to_core(rxn)
                else:
                    self.add_reaction_to_edge(rxn)
            else:
                # Add the reaction to the appropriate unimolecular reaction network
                # If pdep_network 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.add_reaction_to_unimolecular_networks(rxn, new_species=new_species, network=pdep_network)
                if isinstance(rxn, LibraryReaction) and not rxn.kinetics.is_pressure_dependent():
                    # If the reaction came from a library, and it does not have PDep kinetics,
                    # 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)

    def apply_thermo_to_species(self, procnum):
        """
        Generate thermo for species. QM calculations are parallelized if requested.
        """
        from rmgpy.rmg.input import get_input

        quantum_mechanics = get_input("quantum_mechanics")

        if quantum_mechanics:
            quantum_mechanics.run_jobs(self.new_species_list, procnum=procnum)

        # Serial thermo calculation for other methods
        for spc in self.new_species_list:
            self.generate_thermo(spc, rename=True)

    def generate_thermo(self, spc, rename=False):
        """
        Generate thermo for species.
        """
        if not spc.thermo:
            submit(spc, self.solvent_name)

            if rename and spc.thermo and spc.thermo.label != "":  # check if thermo libraries have a name for it
                if isinstance(spc.molecule[0], Fragment):
                    logging.info("Species {0} NOT renamed {1} but get thermo based on thermo library".format(spc.label, spc.thermo.label))
                    spc.label = spc.smiles
                else:
                    logging.info("Species {0} renamed {1} based on thermo library name".format(spc.label, spc.thermo.label))
                    spc.label = spc.thermo.label

        if vapor_liquid_mass_transfer.enabled:
            spc.get_liquid_volumetric_mass_transfer_coefficient_data()
            spc.get_henry_law_constant_data()

        spc.generate_energy_transfer_model()

    def process_coverage_dependence(self, kinetics):
        """Process the coverage dependence kinetics.

        This ensures that species that are used in coverage-dependent kinetic
        expressions exist in the model. (Before this is called, they may have
        only existed in a reaction libary instance).

        If <CoreEdgeReactionModel>self.coverage_dependence is False then
        it instead removes any coverage_dependence from the kinetics.
        """
        if not self.coverage_dependence:
            kinetics.coverage_dependence = None
            return
        cov_dep = {}
        for species, values in kinetics.coverage_dependence.items():
            species_in_model, is_new = self.make_new_species(species)
            cov_dep[species_in_model] = values
        kinetics.coverage_dependence = cov_dep

    def apply_kinetics_to_reaction(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 get_db

        # Find the reaction kinetics
        kinetics, source, entry, is_forward = self.generate_kinetics(reaction)
        # Flip the reaction direction if the kinetics are defined in the reverse direction
        if not is_forward:
            family = get_db("kinetics").families[reaction.family]
            reaction.reactants, reaction.products = reaction.products, reaction.reactants
            reaction.pairs = [(p, r) for r, p in reaction.pairs]
            if family.own_reverse 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 generate_kinetics(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 = get_family_library_object(reaction.family)

        # Get the kinetics for the reaction
        kinetics, source, entry, is_forward = family.get_kinetics(
            reaction, template_labels=reaction.template, degeneracy=reaction.degeneracy, estimator=self.kinetics_estimator, return_all_kinetics=False
        )
        # Get the gibbs free energy of reaction at 298 K
        G298 = reaction.get_free_energy_of_reaction(298)
        gibbs_is_positive = G298 > -1e-8

        if family.own_reverse and len(reaction.products) == len(reaction.reactants) 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_is_forward = family.get_kinetics(
                    reaction.reverse,
                    template_labels=reaction.reverse.template,
                    degeneracy=reaction.reverse.degeneracy,
                    estimator=self.kinetics_estimator,
                    return_all_kinetics=False,
                )
                # Now decide which direction's kinetics to keep
                keep_reverse = 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.
                    keep_reverse = 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
                    keep_reverse = gibbs_is_positive and is_forward and rev_is_forward
                    reason = "Both directions matched the same entry in {0}, but this direction is exergonic.".format(reaction.family)
                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:
                        keep_reverse = 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:
                        keep_reverse = 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:
                        keep_reverse = gibbs_is_positive and is_forward and rev_is_forward
                        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
                    keep_reverse = gibbs_is_positive and is_forward and rev_is_forward
                    reason = "Both directions are estimates, but this direction is exergonic."

                if keep_reverse:
                    kinetics = rev_kinetics
                    source = rev_source
                    entry = rev_entry
                    is_forward = not rev_is_forward
                    G298 = -G298

                if self.verbose_comments:
                    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.0)

        # 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.verbose_comments:
            # 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, is_forward

    def log_enlarge_summary(
        self, new_core_species, new_core_reactions, new_edge_species, new_edge_reactions, reactions_moved_from_edge=None, react_edge=False
    ):
        """
        Output a summary of a model enlargement step to the log. The details of
        the enlargement are passed in the `new_core_species`, `new_core_reactions`,
        `new_edge_species`, and `new_edge_reactions` objects.
        """

        logging.info("")
        if react_edge:
            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(new_core_species)))
        for spec in new_core_species:
            display(spec)
            logging.info("    {0}".format(spec))

        logging.info("Created {0:d} new edge species".format(len(new_edge_species)))
        for spec in new_edge_species:
            display(spec)
            logging.info("    {0}".format(spec))

        if reactions_moved_from_edge:
            logging.info("Moved {0:d} reactions from edge to core".format(len(reactions_moved_from_edge)))
            for rxn in reactions_moved_from_edge:
                for r in new_core_reactions:
                    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))
                        new_core_reactions.remove(r)
                        break
        logging.info("Added {0:d} new core reactions".format(len(new_core_reactions)))
        for rxn in new_core_reactions:
            logging.info("    {0}".format(rxn))

        logging.info("Created {0:d} new edge reactions".format(len(new_edge_reactions)))
        for rxn in new_edge_reactions:
            logging.info("    {0}".format(rxn))

        core_species_count, core_reaction_count, edge_species_count, edge_reaction_count = self.get_model_size()

        # 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(core_species_count, core_reaction_count))
        logging.info("    The model edge has {0:d} species and {1:d} reactions".format(edge_species_count, edge_reaction_count))
        logging.info("")

    def add_species_to_core(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)

        forbidden_structures = get_db("forbidden")

        # check RMG globally forbidden structures
        if not spec.explicitly_allowed and forbidden_structures.is_molecule_forbidden(spec.molecule[0]):
            rxn_list = []
            if spec in self.edge.species:
                # remove forbidden species from edge
                logging.info("Species {0} was Forbidden and not added to Core...Removing from Edge.".format(spec))
                self.remove_species_from_edge(self.reaction_systems, spec)
                # remove any empty pdep networks as a result of species removal
                if self.pressure_dependence:
                    self.remove_empty_pdep_networks()

                return []

        # Add the species to the core
        self.core.species.append(spec)

        rxn_list = []
        if spec in self.edge.species:
            self.edge.phase_system.pass_species(spec.label, self.core.phase_system)
            # If species was in edge, remove it
            logging.debug("Removing species %s from edge.", 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:
                all_core = True
                for reactant in rxn.reactants:
                    if reactant not in self.core.species:
                        all_core = False
                for product in rxn.products:
                    if product not in self.core.species:
                        all_core = False
                if all_core:
                    rxn_list.append(rxn)

            # Move any identified reactions to the core
            for rxn in rxn_list:
                self.add_reaction_to_core(rxn)
                logging.debug("Moving reaction from edge to core: %s", rxn)
        else:
            if spec.molecule[0].contains_surface_site():
                self.core.phase_system.phases["Surface"].add_species(spec, edge_phase=self.edge.phase_system.phases["Surface"])
                self.edge.phase_system.species_dict[spec.label] = spec
                self.core.phase_system.species_dict[spec.label] = spec
            else:
                self.core.phase_system.phases["Default"].add_species(spec, edge_phase=self.edge.phase_system.phases["Default"])
                self.edge.phase_system.species_dict[spec.label] = spec
                self.core.phase_system.species_dict[spec.label] = spec

        return rxn_list

    def add_species_to_edge(self, spec):
        """
        Add a species `spec` to the reaction model edge.
        """
        self.edge.species.append(spec)
        if spec.molecule[0].contains_surface_site():
            self.edge.phase_system.phases["Surface"].add_species(spec)
            self.edge.phase_system.species_dict[spec.label] = spec
        else:
            self.edge.phase_system.phases["Default"].add_species(spec)
            self.edge.phase_system.species_dict[spec.label] = spec

    def set_thermodynamic_filtering_parameters(
        self, Tmax, thermo_tol_keep_spc_in_edge, min_core_size_for_prune, maximum_edge_species, reaction_systems
    ):
        """
        sets parameters for thermodynamic filtering based on the current core
        Tmax is the maximum reactor temperature in K
        thermo_tol_keep_spc_in_edge is the Gibbs number above which species will be filtered
        min_core_size_for_prune is the core size at which thermodynamic filtering will start
        maximum_edge_species is the maximum allowed number of edge species
        reaction_systems is a list of reaction_system objects
        """
        self.Tmax = Tmax
        Gs = [spc.thermo.get_free_energy(Tmax) for spc in self.core.species]
        self.Gmax = max(Gs)
        self.Gmin = min(Gs)

        self.Gfmax = thermo_tol_keep_spc_in_edge * (self.Gmax - self.Gmin) + self.Gmax
        self.thermo_tol_keep_spc_in_edge = thermo_tol_keep_spc_in_edge
        self.min_core_size_for_prune = min_core_size_for_prune
        self.reaction_systems = reaction_systems
        self.maximum_edge_species = maximum_edge_species

    def thermo_filter_species(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.get_free_energy(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 thermo_tol_keep_spc_in_edge of "
                    "{3} ".format(spc, G, Gn, self.thermo_tol_keep_spc_in_edge)
                )
                self.remove_species_from_edge(self.reaction_systems, spc)

        # Delete any networks that became empty as a result of pruning
        if self.pressure_dependence:
            self.remove_empty_pdep_networks()

    def thermo_filter_down(self, maximum_edge_species, min_species_exist_iterations_for_prune=0):
        """
        removes species from the edge based on their Gibbs energy until maximum_edge_species
        is reached under the constraint that all removed species are older than
        min_species_exist_iterations_for_prune iterations
        maximum_edge_species is the maximum allowed number of edge species
        min_species_exist_iterations_for_prune is the number of iterations a species must be in the edge
        before it is eligible for thermo filtering
        """
        Tmax = self.Tmax
        num_to_remove = len(self.edge.species) - maximum_edge_species
        logging.debug("Planning to remove %d species", num_to_remove)
        iteration = self.iteration_num

        if num_to_remove > 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 = np.array([spc.thermo.get_free_energy(Tmax) for spc in spcs])
            Gns = (Gfs - self.Gmax) / (self.Gmax - self.Gmin)
            inds = np.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
            remove_spcs = []

            rInds = []
            while ind < len(inds) and num_to_remove > 0:  # find the species we can remove and collect indices for removal
                i = inds[ind]
                spc = spcs[i]
                if iteration - spc.creation_iteration >= min_species_exist_iterations_for_prune:
                    remove_spcs.append(spc)
                    rInds.append(i)
                    num_to_remove -= 1
                ind += 1

            logging.debug("Found %d eligible species for filtering", len(remove_spcs))

            for i, spc in enumerate(remove_spcs):
                logging.info(
                    "Removing species {0} from edge to meet maximum number of edge species, Gibbs " "number is {1}".format(spc, Gns[rInds[i]])
                )
                self.remove_species_from_edge(self.reaction_systems, spc)

            # Delete any networks that became empty as a result of pruning
            if self.pressure_dependence:
                self.remove_empty_pdep_networks()

            # call garbage collection
            collected = gc.collect()
            logging.info("Garbage collector: collected %d objects." % (collected))

    def remove_empty_pdep_networks(self):
        """
        searches for and deletes any empty pdep networks
        """
        networks_to_delete = []
        for network in self.network_list:
            if len(network.path_reactions) == 0 and len(network.net_reactions) == 0:
                networks_to_delete.append(network)

        if len(networks_to_delete) > 0:
            logging.info("Deleting {0:d} empty pressure-dependent reaction networks".format(len(networks_to_delete)))
            for network in networks_to_delete:
                logging.debug("    Deleting empty pressure dependent reaction network #%d", network.index)
                source = tuple(network.source)
                nets_with_this_source = self.network_dict[source]
                nets_with_this_source.remove(network)
                if not nets_with_this_source:
                    del self.network_dict[source]
                self.network_list.remove(network)

    def prune(self, reaction_systems, tol_keep_in_edge, tol_move_to_core, maximum_edge_species, min_species_exist_iterations_for_prune):
        """
        Remove species from the model edge based on the simulation results from
        the list of `reaction_systems`.
        """

        ineligible_species = []  # A list of the species which are not eligible for pruning, for any reason
        prunable_species = reaction_systems[0].prunable_species
        prunable_networks = reaction_systems[0].prunable_networks

        num_prunable_species = len(prunable_species)
        iteration = self.iteration_num
        # All edge species that have not existed for more than two enlarge
        # iterations are ineligible for pruning
        for spec in prunable_species:
            if iteration - spec.creation_iteration <= min_species_exist_iterations_for_prune:
                ineligible_species.append(spec)

        # Get the maximum species rates (and network leak rates)
        # across all reaction systems
        max_edge_species_rate_ratios = np.zeros((num_prunable_species), float)
        for reaction_system in reaction_systems:
            for i in range(num_prunable_species):
                rate_ratio = reaction_system.max_edge_species_rate_ratios[i]
                if max_edge_species_rate_ratios[i] < rate_ratio:
                    max_edge_species_rate_ratios[i] = rate_ratio

            for i, network in enumerate(prunable_networks):
                rate_ratio = reaction_system.max_network_leak_rate_ratios[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.get_leak_branching_ratios(reaction_system.T.value_si, reaction_system.P.value_si)
                for spec, frac in ratios.items():
                    if spec in prunable_species:
                        index = prunable_species.index(spec)
                        max_edge_species_rate_ratios[index] += frac * rate_ratio
                # Mark any species that is explored in any partial network as ineligible for pruning
                for spec in network.explored:
                    if spec not in ineligible_species:
                        ineligible_species.append(spec)

        # Sort the edge species rates by index
        indices = np.argsort(max_edge_species_rate_ratios)
        # Determine which species to prune
        species_to_prune = []
        prune_due_to_rate_counter = 0
        for index in indices:
            spec = prunable_species[index]
            if spec in ineligible_species or not spec in self.edge.species:
                continue
            # Remove the species with rates below the pruning tolerance from the model edge
            if max_edge_species_rate_ratios[index] < tol_keep_in_edge:
                species_to_prune.append((index, spec))
                prune_due_to_rate_counter += 1
            # Keep removing species with the lowest rates until we are below the maximum edge species size
            elif num_prunable_species - len(species_to_prune) > maximum_edge_species:
                if max_edge_species_rate_ratios[index] < tol_move_to_core:
                    logging.info("Pruning species {0} to make num_edge_species smaller than maximum_edge_species".format(spec))
                    species_to_prune.append((index, spec))
                else:
                    logging.warning(
                        "Attempted to prune a species that exceeded tol_move_to_core, pruning settings "
                        "for this run are likely bad, either maximum_edge_species needs to be set higher "
                        "(~100000) or min_species_exist_iterations_for_prune should be reduced (~2)"
                    )
                    break
            else:
                break

        # Actually do the pruning
        if prune_due_to_rate_counter > 0:
            logging.info(
                "Pruning %d species whose rate ratios against characteristic rate did not exceed the " "minimum threshold of %g",
                prune_due_to_rate_counter,
                tol_keep_in_edge,
            )
            for index, spec in species_to_prune[0:prune_due_to_rate_counter]:
                logging.info("Pruning species %s", spec)
                logging.debug("    %-56s    %10.4e", spec, max_edge_species_rate_ratios[index])
                self.remove_species_from_edge(reaction_systems, spec)
        if len(species_to_prune) - prune_due_to_rate_counter > 0:
            logging.info(
                "Pruning %d species to obtain an edge size of %d species", len(species_to_prune) - prune_due_to_rate_counter, maximum_edge_species
            )
            for index, spec in species_to_prune[prune_due_to_rate_counter:]:
                logging.info("Pruning species %s", spec)
                logging.debug("    %-56s    %10.4e", spec, max_edge_species_rate_ratios[index])
                self.remove_species_from_edge(reaction_systems, spec)

        # Delete any networks that became empty as a result of pruning
        if self.pressure_dependence:
            self.remove_empty_pdep_networks()

        logging.info("")

    def remove_species_from_edge(self, reaction_systems, spec):
        """
        Remove species `spec` from the reaction model edge.
        """

        # remove the species
        self.edge.species.remove(spec)
        self.index_species_dict.pop(spec.index)
        self.edge.phase_system.remove_species(spec)

        # clean up species references in reaction_systems
        for reaction_system in reaction_systems:
            if not isinstance(reaction_system, Reactor):
                try:
                    reaction_system.species_index.pop(spec)
                except KeyError:
                    pass

                # identify any reactions it's involved in
                rxn_list = []
                for rxn in reaction_system.reaction_index:
                    if spec in rxn.reactants or spec in rxn.products:
                        rxn_list.append(rxn)

                for rxn in rxn_list:
                    reaction_system.reaction_index.pop(rxn)

        # identify any reactions it's involved in
        rxn_list = []
        for rxn in self.edge.reactions:
            if spec in rxn.reactants or spec in rxn.products:
                rxn_list.append(rxn)
        # remove those reactions
        for rxn in rxn_list:
            self.edge.reactions.remove(rxn)

        # Remove the species from any unirxn networks it is in
        if self.pressure_dependence:
            for network in self.network_list:
                # Delete all path reactions involving the species
                rxn_list = []
                for rxn in network.path_reactions:
                    if spec in rxn.reactants or spec in rxn.products:
                        rxn_list.append(rxn)
                for rxn in rxn_list:
                    network.path_reactions.remove(rxn)
                # Delete all net reactions involving the species
                rxn_list = []
                for rxn in network.net_reactions:
                    if spec in rxn.reactants or spec in rxn.products:
                        rxn_list.append(rxn)
                for rxn in rxn_list:
                    network.net_reactions.remove(rxn)
                # Remove the species from list of explored
                if spec in network.explored:
                    network.explored.remove(spec)

                # Recompute the isomers, reactants, and products for this network
                network.update_configurations(self)

        # Remove from the global list of reactions
        # also remove it from the global list of reactions
        for family in self.reaction_dict:
            if spec in self.reaction_dict[family]:
                del self.reaction_dict[family][spec]
            for reactant1 in self.reaction_dict[family]:
                if spec in self.reaction_dict[family][reactant1]:
                    del self.reaction_dict[family][reactant1][spec]
            for reactant1 in self.reaction_dict[family]:
                for reactant2 in self.reaction_dict[family][reactant1]:
                    temp_rxn_delete_list = []
                    for templateReaction in self.reaction_dict[family][reactant1][reactant2]:
                        if spec in templateReaction.reactants or spec in templateReaction.products:
                            temp_rxn_delete_list.append(templateReaction)
                    for tempRxnToBeDeleted in temp_rxn_delete_list:
                        self.reaction_dict[family][reactant1][reactant2].remove(tempRxnToBeDeleted)

        # remove from the global list of species, to free memory
        formula = spec.molecule[0].get_formula()
        self.species_dict[formula].remove(spec)
        if spec in self.species_cache:
            self.species_cache.remove(spec)
            self.species_cache.append(None)

    def add_reaction_to_core(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 not in self.edge.reactions:
                # If a reaction is not in edge but is going to add to core, it is either a seed mechanism or a newly generated reaction where all reactants and products are already in core
                # If the reaction is in edge, then the corresponding rms_rxn was moved from edge phase to core phase in pass_species already.
                rms_species_list = self.core.phase_system.get_rms_species_list()
                species_names = self.core.phase_system.get_species_names()
                bits = np.array([spc.molecule[0].contains_surface_site() for spc in rxn.reactants + rxn.products])
                if all(bits):
                    self.core.phase_system.phases["Surface"].add_reaction(rxn, self.edge.phase_system.phases["Surface"])
                elif all(bits == False):
                    self.core.phase_system.phases["Default"].add_reaction(rxn, self.edge.phase_system.phases["Default"])
                else:
                    self.core.phase_system.interfaces[frozenset({"Default", "Surface"})].add_reaction(rxn, species_names, rms_species_list)
                    self.edge.phase_system.interfaces[frozenset({"Default", "Surface"})].add_reaction(rxn, species_names, rms_species_list)

        if rxn in self.edge.reactions:
            self.edge.reactions.remove(rxn)

    def add_reaction_to_edge(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)
        rms_species_list = self.edge.phase_system.get_rms_species_list()
        species_names = self.edge.phase_system.get_species_names()
        bits = np.array([spc.molecule[0].contains_surface_site() for spc in rxn.reactants + rxn.products])
        if all(bits):
            self.edge.phase_system.phases["Surface"].add_reaction(rxn)
        elif all(bits == False):
            self.edge.phase_system.phases["Default"].add_reaction(rxn)
        else:
            self.edge.phase_system.interfaces[frozenset(["Default", "Surface"])].add_reaction(rxn, species_names, rms_species_list)

    def get_model_size(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.
        """
        core_species_count = len(self.core.species)
        core_reactions_count = len(self.core.reactions)
        edge_species_count = len(self.edge.species)
        edge_reactions_count = len(self.edge.reactions)
        return core_species_count, core_reactions_count, edge_species_count, edge_reactions_count

    def get_species_reaction_lists(self):
        """
        Return lists of all of the species and reactions in the core and the
        edge.
        """
        species_list = []
        species_list.extend(self.core.species)
        species_list.extend(self.edge.species)
        reaction_list = []
        reaction_list.extend(self.core.reactions)
        reaction_list.extend(self.edge.reactions)
        return species_list, reaction_list

    def get_stoichiometry_matrix(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.
        """
        species_list, reaction_list = self.get_species_reaction_lists()
        from scipy import sparse

        stoichiometry = sparse.dok_matrix((self.species_counter, self.reaction_counter), float)
        for rxn in reaction_list:
            j = rxn.index - 1
            spec_list = rxn.reactants[:]
            spec_list.extend(rxn.products)
            for spec in spec_list:
                i = spec.index - 1
                nu = rxn.get_stoichiometric_coefficient(spec)
                if nu != 0:
                    stoichiometry[i, j] = nu
        return stoichiometry.tocsr()

    def add_seed_mechanism_to_core(self, seed_mechanism, react=False):
        """
        Add all species and reactions from `seed_mechanism`, 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

        library_names = list(database.kinetics.libraries.keys())
        family_names = list(database.kinetics.families.keys())

        path = os.path.join(settings["database.directory"], "kinetics", "libraries")
        from rmgpy.rmg.input import rmg

        self.new_reaction_list = []
        self.new_species_list = []
        edge_species_to_move = []
        
        num_old_core_species = len(self.core.species)
        num_old_core_reactions = len(self.core.reactions)

        logging.info("Adding seed mechanism {0} to model core...".format(seed_mechanism))

        seed_mechanism = database.kinetics.libraries[seed_mechanism]

        rxns = seed_mechanism.get_library_reactions()

        for rxn in rxns:
            if (
                isinstance(rxn, LibraryReaction) and not (rxn.library in library_names) and not (rxn.library == "kineticsjobs")
            ):  # if one of the reactions in the library is from another library load that library
                database.kinetics.library_order.append((rxn.library, "Internal"))
                database.kinetics.load_libraries(path=path, libraries=[rxn.library])
                library_names = list(database.kinetics.libraries.keys())
            if isinstance(rxn, TemplateReaction) and not (rxn.family in family_names):
                logging.warning("loading reaction {0} originally from family {1} as a library reaction".format(str(rxn), rxn.family))
                rxn = LibraryReaction(
                    reactants=rxn.reactants[:],
                    products=rxn.products[:],
                    library=seed_mechanism.name,
                    specific_collider=rxn.specific_collider,
                    kinetics=rxn.kinetics,
                    duplicate=rxn.duplicate,
                    reversible=rxn.reversible,
                )
            r, isNew = self.make_new_reaction(rxn)  # updates self.new_species_list and self.new_reaction_list
            for s in rxn.reactants+rxn.products:
                if s in self.edge.species and s not in edge_species_to_move:
                    edge_species_to_move.append(s)
            if getattr(r.kinetics, "coverage_dependence", None):
                self.process_coverage_dependence(r.kinetics)
            if not isNew:
                logging.info("This library reaction was not new: {0}".format(rxn))
            elif (
                self.pressure_dependence
                and rxn.elementary_high_p
                and rxn.is_unimolecular()
                and isinstance(rxn, LibraryReaction)
                and isinstance(rxn.kinetics, Arrhenius)
                and (
                    self.pressure_dependence.maximum_atoms is None
                    or self.pressure_dependence.maximum_atoms >= sum([len(spec.molecule[0].atoms) for spec in r.reactants])
                )
            ):
                # This unimolecular library reaction is flagged as `elementary_high_p` and has Arrhenius type kinetics.
                # We should calculate a pressure-dependent rate for it
                if len(rxn.reactants) == 1:
                    self.process_new_reactions(new_reactions=[rxn], new_species=rxn.reactants[0])
                else:
                    self.process_new_reactions(new_reactions=[rxn], new_species=rxn.products[0])

        # Perform species constraints and forbidden species checks

        for spec in self.new_species_list:
            if database.forbidden_structures.is_molecule_forbidden(spec.molecule[0]):
                if "allowed" in rmg.species_constraints and "seed mechanisms" in rmg.species_constraints["allowed"]:
                    spec.explicitly_allowed = True
                    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, seed_mechanism.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, seed_mechanism.label)
                    )
            if fails_species_constraints(spec):
                if "allowed" in rmg.species_constraints and "seed mechanisms" in rmg.species_constraints["allowed"]:
                    rmg.species_constraints["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, seed_mechanism.label)
                    )

        for spec in edge_species_to_move+self.new_species_list:
            if spec.reactive:
                submit(spec, self.solvent_name)
            if vapor_liquid_mass_transfer.enabled:
                spec.get_liquid_volumetric_mass_transfer_coefficient_data()
                spec.get_henry_law_constant_data()

            self.add_species_to_core(spec)

        for rxn in self.new_reaction_list:
            if self.pressure_dependence and rxn.is_unimolecular():
                # 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.solvent_name)

                rxn.fix_barrier_height(force_positive=True)
            self.add_reaction_to_core(rxn)

        # Check we didn't introduce unmarked duplicates
        self.mark_chemkin_duplicates()

        self.log_enlarge_summary(
            new_core_species=self.core.species[num_old_core_species:],
            new_core_reactions=self.core.reactions[num_old_core_reactions:],
            new_edge_species=[],
            new_edge_reactions=[],
        )

    def add_reaction_library_to_edge(self, reaction_library):
        """
        Add all species and reactions from `reaction_library`, a
        :class:`KineticsPrimaryDatabase` object, to the model edge.
        """

        database = rmgpy.data.rmg.database
        library_names = list(database.kinetics.libraries.keys())
        family_names = list(database.kinetics.families.keys())
        path = os.path.join(settings["database.directory"], "kinetics", "libraries")

        from rmgpy.rmg.input import rmg

        self.new_reaction_list = []
        self.new_species_list = []

        num_old_edge_species = len(self.edge.species)
        num_old_edge_reactions = len(self.edge.reactions)

        logging.info("Adding reaction library {0} to model edge...".format(reaction_library))
        if reaction_library in database.kinetics.libraries:
            reaction_library = database.kinetics.libraries[reaction_library]
        elif reaction_library in database.kinetics.external_library_labels:
            reaction_library = database.kinetics.libraries[database.kinetics.external_library_labels[reaction_library]]
        else:
            raise ValueError(f'Library {reaction_library} not found.')

        rxns = reaction_library.get_library_reactions()
        for rxn in rxns:
            if isinstance(rxn, LibraryReaction) and not (
                rxn.library in library_names
            ):  # if one of the reactions in the library is from another library load that library
                database.kinetics.library_order.append((rxn.library, "Internal"))
                database.kinetics.load_libraries(path=path, libraries=[rxn.library])
                library_names = list(database.kinetics.libraries.keys())
            if isinstance(rxn, TemplateReaction) and not (rxn.family in family_names):
                logging.warning("loading reaction {0} originally from family {1} as a library reaction".format(str(rxn), rxn.family))
                rxn = LibraryReaction(
                    reactants=rxn.reactants[:],
                    products=rxn.products[:],
                    library=reaction_library.name,
                    specific_collider=rxn.specific_collider,
                    kinetics=rxn.kinetics,
                    duplicate=rxn.duplicate,
                    reversible=rxn.reversible,
                )
            r, isNew = self.make_new_reaction(rxn)  # updates self.new_species_list and self.new_reaction_list
            if r is not None and getattr(rxn.kinetics, "coverage_dependence", None):
                self.process_coverage_dependence(r.kinetics)
            if not isNew:
                logging.info("This library reaction was not new: {0}".format(rxn))
            elif (
                self.pressure_dependence
                and rxn.elementary_high_p
                and rxn.is_unimolecular()
                and isinstance(rxn, LibraryReaction)
                and isinstance(rxn.kinetics, Arrhenius)
                and (
                    self.pressure_dependence.maximum_atoms is None
                    or self.pressure_dependence.maximum_atoms >= sum([len(spec.molecule[0].atoms) for spec in r.reactants])
                )
            ):
                # This unimolecular library reaction is flagged as `elementary_high_p` and has Arrhenius type kinetics.
                # We should calculate a pressure-dependent rate for it
                if len(rxn.reactants) == 1:
                    self.process_new_reactions(new_reactions=[rxn], new_species=rxn.reactants[0])
                else:
                    self.process_new_reactions(new_reactions=[rxn], new_species=rxn.products[0])

        # Perform species constraints and forbidden species checks
        for spec in self.new_species_list:
            if not reaction_library.auto_generated:  # No need to check for forbidden species otherwise
                if database.forbidden_structures.is_molecule_forbidden(spec.molecule[0]):
                    if "allowed" in rmg.species_constraints and "reaction libraries" in rmg.species_constraints["allowed"]:
                        spec.explicitly_allowed = True
                        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, reaction_library.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, reaction_library.label)
                        )
            if fails_species_constraints(spec):
                if "allowed" in rmg.species_constraints and "reaction libraries" in rmg.species_constraints["allowed"]:
                    rmg.species_constraints["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, reaction_library.label)
                    )

        for spec in self.new_species_list:
            if spec.reactive:
                submit(spec, self.solvent_name)
            if vapor_liquid_mass_transfer.enabled:
                spec.get_liquid_volumetric_mass_transfer_coefficient_data()
                spec.get_henry_law_constant_data()

            self.add_species_to_edge(spec)

        for rxn in self.new_reaction_list:
            if not (
                self.pressure_dependence
                and rxn.elementary_high_p
                and rxn.is_unimolecular()
                and isinstance(rxn, LibraryReaction)
                and isinstance(rxn.kinetics, Arrhenius)
                and (
                    self.pressure_dependence.maximum_atoms is None
                    or self.pressure_dependence.maximum_atoms >= sum([len(spec.molecule[0].atoms) for spec in r.reactants])
                )
            ):
                # Don't add to the edge library reactions that were already processed
                self.add_reaction_to_edge(rxn)

        if self.save_edge_species:
            from rmgpy.chemkin import mark_duplicate_reaction

            new_edge_reactions = self.edge.reactions[num_old_edge_reactions:]
            checked_reactions = self.core.reactions + self.edge.reactions[:num_old_edge_reactions]
            for rxn in new_edge_reactions:
                mark_duplicate_reaction(rxn, checked_reactions)
                checked_reactions.append(rxn)

        self.log_enlarge_summary(
            new_core_species=[],
            new_core_reactions=[],
            new_edge_species=self.edge.species[num_old_edge_species:],
            new_edge_reactions=self.edge.reactions[num_old_edge_reactions:],
        )

    def add_reaction_library_to_output(self, reaction_library):
        """
        Add all species and reactions from `reaction_library`, 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(reaction_library))

        # 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 == reaction_library:
                    self.output_reaction_list.append(rxn)

                    for species in rxn.reactants + rxn.products:
                        if species not in self.core.species and species not in self.output_species_list:
                            self.output_species_list.append(species)

    def add_reaction_to_unimolecular_networks(self, newReaction, new_species, 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(new_species, Species)

        # Put the reaction in the direction in which the new species is in the reactants
        if new_species 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.network_dict[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.network_dict[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.network_count += 1
                network = PDepNetwork(index=self.network_count, source=reactants[:])
                # should the source passed to PDepNetwork constuctor be a tuple not a list? that's what is used in network_dict
                try:
                    self.network_dict[source].append(network)
                except KeyError:
                    self.network_dict[source] = [network]
                self.network_list.append(network)

        # Add the path reaction to that network
        network.add_path_reaction(newReaction)

    def update_unimolecular_reaction_networks(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.network_dict.values():
            network_count = len(networks)
            for index0, network0 in enumerate(networks):
                index = index0 + 1
                while index < network_count:
                    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.network_list.remove(network)
                        network_count -= 1
                    else:
                        index += 1

        count = sum([1 for network in self.network_list 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.network_list)))

        # Iterate over all the networks, updating the invalid ones as necessary
        # self = reaction_model object
        updated_networks = []
        for network in self.network_list:
            if not network.valid:
                network.update(self, self.pressure_dependence)
                updated_networks.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
        core_reaction_count = len(self.core.reactions)
        while index < core_reaction_count:
            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.get_free_energy_of_reaction(300.0)
                        kf = reaction.get_rate_coefficient(1000, 1e5)
                        kr = reaction.get_rate_coefficient(1000, 1e5) / reaction.get_equilibrium_constant(1000)
                        kf2 = reaction2.get_rate_coefficient(1000, 1e5) / reaction2.get_equilibrium_constant(1000)
                        kr2 = reaction2.get_rate_coefficient(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
                        keep_first = dGrxn < 0
                        # Delete the PDepReaction that we aren't keeping
                        if keep_first:
                            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
                        core_reaction_count -= 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 mark_chemkin_duplicates(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 mark_duplicate_reactions

        rxn_list = self.core.reactions + self.output_reaction_list
        mark_duplicate_reactions(rxn_list)

    def register_reaction(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 = generate_reaction_key(rxn)

        # make dictionary entries if necessary
        if key_family not in self.reaction_dict:
            self.reaction_dict[key_family] = {}

        if key1 not in self.reaction_dict[key_family]:
            self.reaction_dict[key_family][key1] = {}

        if key2 not in self.reaction_dict[key_family][key1]:
            self.reaction_dict[key_family][key1][key2] = []

        # store this reaction at the top of the relevant short-list
        self.reaction_dict[key_family][key1][key2].insert(0, rxn)

    def search_retrieve_reactions(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 = generate_reaction_key(rxn)

        my_reaction_list = []

        rxns = self.retrieve(family_label, r1_fwd, r2_fwd)
        my_reaction_list.extend(rxns)

        family = get_family_library_object(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 = generate_reaction_key(rxn, useProducts=True)

            rxns = self.retrieve(family_label, r1_rev, r2_rev)
            my_reaction_list.extend(rxns)

        return my_reaction_list

    def initialize_index_species_dict(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.index_species_dict[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.reaction_dict[family_label][key1][key2][:]
        except KeyError:  # no such short-list: must be new, unless in seed.
            return []


def generate_reaction_key(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

    if len(spc_list) == 1:
        key1, key2 = get_key(spc_list[0]), None
    else:
        # Use 2 keys, even for trimolecular reactions?
        key1, key2 = sorted([get_key(spc_list[0]), get_key(spc_list[1])])

    return key_family, key1, key2


def generate_reaction_id(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([get_key(reactant) for reactant in rxn.reactants])
    products = sorted([get_key(product) for product in rxn.products])

    return reactants, products


def get_family_library_object(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:
        pass

    try:
        lib = kinetics.libraries[label]
        return lib
    except KeyError:
        pass

    raise Exception(
        f"Could not retrieve the family/library: {label} "
        f"from families {list(kinetics.families.keys())} "
        f"and libraries {list(kinetics.libraries.keys())}."
    )


def get_key(spc):
    """
    Returns a string of the species that can serve as a key in a dictionary.
    """

    return spc.label


def are_identical_species_references(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.specific_collider == rxn2.specific_collider

    return (identical_same_direction or identical_opposite_directions) and identical_collider