<commit_before>#!/usr/bin/python
# -*- coding: utf-8 -*-

################################################################################
#
#	RMG - Reaction Mechanism Generator
#
#	Copyright (c) 2002-2009 Prof. <NAME> (<EMAIL>) and the
#	RMG Team (<EMAIL>)
#
#	Permission is hereby granted, free of charge, to any person obtaining a
#	copy of this software and associated documentation files (the 'Software'),
#	to deal in the Software without restriction, including without limitation
#	the rights to use, copy, modify, merge, publish, distribute, sublicense,
#	and/or sell copies of the Software, and to permit persons to whom the
#	Software is furnished to do so, subject to the following conditions:
#
#	The above copyright notice and this permission notice shall be included in
#	all copies or substantial portions of the Software.
#
#	THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
#	FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
#	DEALINGS IN THE SOFTWARE.
#
################################################################################

"""
Contains classes for working with the reaction model generated by RMG.
"""

import logging
import math
import numpy
import os
import os.path

from rmgpy.quantity import Quantity, constants
import rmgpy.species
from rmgpy.thermo import Wilhoit, MultiNASA

from rmgpy.data.thermo import *
from rmgpy.data.kinetics import *
from rmgpy.data.states import *

import rmgpy.measure.network

import settings

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

class Species(rmgpy.species.Species):

    def __init__(self, index=-1, label='', thermo=None, states=None, molecule=None, E0=0.0, lennardJones=None, molecularWeight=0.0, reactive=True):
        rmgpy.species.Species.__init__(self, index, label, thermo, states, molecule, E0, lennardJones, molecularWeight, reactive)
        self.coreSizeAtCreation = 0

    def __reduce__(self):
        """
        A helper function used when pickling an object.
        """
        d = {
            'coreSizeAtCreation': self.coreSizeAtCreation,
        }
        return (Species, (self.index, self.label, self.thermo, self.states, self.molecule, self.E0, self.lennardJones, self.molecularWeight, self.reactive), d)

    def __setstate__(self, d):
        """
        A helper function used when unpickling an object.
        """
        self.coreSizeAtCreation = d['coreSizeAtCreation']

    def generateThermoData(self, database, thermoClass=MultiNASA):
        """
        Generate thermodynamic data for the species using the thermo database.

        Generates the thermo data for each structure (resonance isomer),
        picks that with lowest H298 value, and saves it to `self.thermoData`.
        """

        # Ensure molecules are using explicit hydrogens
        implicitH = [mol.implicitHydrogens for mol in self.molecule]
        for molecule in self.molecule:
            molecule.makeHydrogensExplicit()
        
        thermo = []
        for molecule in self.molecule:
            molecule.clearLabeledAtoms()
            molecule.updateAtomTypes()
            tdata = database.thermo.getThermoData(molecule)
            thermo.append(tdata)

        H298 = numpy.array([t.getEnthalpy(298.) for t in thermo])
        indices = H298.argsort()

        # If multiple resonance isomers are present, use the thermo data of
        # the most stable isomer (i.e. one with lowest enthalpy of formation)
        # as the thermo data of the species
        self.thermo = thermo[indices[0]]

        # Sort the structures in order of decreasing stability
        self.molecule = [self.molecule[ind] for ind in indices]
        implicitH = [implicitH[ind] for ind in indices]

        # Convert to desired thermo class
        thermo0 = self.thermo
        linear = self.molecule[0].isLinear()
        nRotors = self.molecule[0].countInternalRotors()
        nFreq = 3 * len(self.molecule[0].atoms) - (5 if linear else 6) - nRotors
        self.thermo = convertThermoModel(self.thermo, Wilhoit, linear=linear, nFreq=nFreq, nRotors=nRotors)
        self.E0 = Quantity(self.thermo.getEnthalpy(1.0)/1000.0,"kJ/mol")
        self.thermo = convertThermoModel(self.thermo, thermoClass, Tmin=100.0, Tmax=5000.0, Tint=1000.0)
        if self.thermo.__class__ != thermo0.__class__:
            # Compute RMS error of overall transformation
            Tlist = numpy.array([300.0, 400.0, 500.0, 600.0, 800.0, 1000.0, 1500.0], numpy.float64)
            err = math.sqrt(numpy.sum((self.thermo.getHeatCapacities(Tlist) - thermo0.getHeatCapacities(Tlist))**2))/constants.R/len(Tlist)
            logging.log(logging.WARNING if err > 0.1 else 0, 'Average RMS error in heat capacity fit to {0} = {1:g}*R'.format(self, err))

        # Restore implicit hydrogens if necessary
        for implicit, molecule in zip(implicitH, self.molecule):
            if implicit: molecule.makeHydrogensImplicit()

        return self.thermo

    def generateStatesData(self, database):
        """
        Generate molecular degree of freedom data for the species. You must
        have already provided a thermodynamics model using e.g.
        :meth:`generateThermoData()`.
        """
        if not self.thermo:
            raise Exception("Unable to determine states model for species {0}: No thermodynamics model found.".format(self))
        molecule = self.molecule[0]
        implicitH = molecule.implicitHydrogens
        molecule.makeHydrogensExplicit()
        self.states = database.states.getStatesData(molecule, self.thermo)
        if implicitH: molecule.makeHydrogensImplicit()

    def generateLennardJonesParameters(self):
        """
        Generate the Lennard-Jones parameters for the species. This "algorithm"
        is *very* much in need of improvement.
        """

        count = sum([1 for atom in self.molecule[0].vertices if atom.isNonHydrogen()])
        self.lennardJones = rmgpy.species.LennardJones()

        if count == 1:
            self.lennardJones.sigma = Quantity(3.758e-10,"m")
            self.lennardJones.epsilon = Quantity(148.6 * constants.kB,"J")
        elif count == 2:
            self.lennardJones.sigma = Quantity(4.443e-10,"m")
            self.lennardJones.epsilon = Quantity(110.7 * constants.kB,"J")
        elif count == 3:
            self.lennardJones.sigma = Quantity(5.118e-10,"m")
            self.lennardJones.epsilon = Quantity(237.1 * constants.kB,"J")
        elif count == 4:
            self.lennardJones.sigma = Quantity(4.687e-10,"m")
            self.lennardJones.epsilon = Quantity(531.4 * constants.kB,"J")
        elif count == 5:
            self.lennardJones.sigma = Quantity(5.784e-10,"m")
            self.lennardJones.epsilon = Quantity(341.1 * constants.kB,"J")
        else:
            self.lennardJones.sigma = Quantity(5.949e-10,"m")
            self.lennardJones.epsilon = Quantity(399.3 * constants.kB,"J")

class PDepReaction(rmgpy.reaction.Reaction):

    def __init__(self, index=-1, reactants=None, products=None, network=None, kinetics=None, reversible=True, transitionState=None, thirdBody=False):
        rmgpy.reaction.Reaction.__init__(self, index, reactants, products, kinetics, reversible, transitionState, thirdBody)
        self.network = network

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

class PressureDependenceError(Exception):
    """
    An exception class to use when an error involving pressure dependence is
    encountered. Pass a string describing the circumstances of the exceptional
    behavior.
    """
    pass

class PDepNetwork(rmgpy.measure.network.Network):
    """
    A representation of a *partial* unimolecular reaction network. Each partial
    network has a single `source` isomer or reactant channel, and is responsible
    only for :math:`k(T,P)` values for net reactions with source as the
    reactant. Multiple partial networks can have the same source, but networks
    with the same source and any explored isomers must be combined.

    =================== ======================= ================================
    Attribute           Type                    Description
    =================== ======================= ================================
    `source`            ``list``                The isomer or reactant channel that acts as the source
    `explored`          ``list``                A list of the unimolecular isomers whose reactions have been fully explored
    =================== ======================= ================================

    """

    def __init__(self, index=-1, source=None):
        rmgpy.measure.network.Network.__init__(self, index=index)
        self.source = source
        self.explored = []

    def getLeakCoefficient(self, T, P):
        """
        Return the pressure-dependent rate coefficient :math:`k(T,P)` describing
        the total rate of "leak" from this network. This is defined as the sum
        of the :math:`k(T,P)` values for all net reactions to nonexplored
        unimolecular isomers.
        """
        k = 0.0
        if len(self.netReactions) == 0 and len(self.pathReactions) == 1:
            # The network is of the form A + B -> C* (with C* nonincluded)
            # For this special case we use the high-pressure limit k(T) to
            # ensure that we're estimating the total leak flux
            rxn = self.pathReactions[0]
            if rxn.kinetics is None:
                if rxn.reverse.kinetics is not None:
                    rxn = rxn.reverse
                else:
                    raise PressureDependenceError('Path reaction {0} with no high-pressure-limit kinetics encountered in PDepNetwork #{1:d} while evaluating leak flux.'.format(rxn, self.index))
            if rxn.products is self.source:
                k = rxn.getRateCoefficient(T,P) / rxn.getEquilibriumConstant(T)
            else:
                k = rxn.getRateCoefficient(T,P)
        else:
            # The network has at least one included isomer, so we can calculate
            # the leak flux normally
            for rxn in self.netReactions:
                if len(rxn.products) == 1 and rxn.products[0] not in self.explored:
                    k += rxn.getRateCoefficient(T,P)
        return k

    def getMaximumLeakSpecies(self, T, P):
        """
        Get the unexplored (unimolecular) isomer with the maximum leak flux.
        Note that the leak rate coefficients vary with temperature and
        pressure, so you must provide these in order to get a meaningful result.
        """
        # Choose species with maximum leak flux
        maxK = 0.0; maxSpecies = None
        if len(self.netReactions) == 0 and len(self.pathReactions) == 1:
            maxK = self.getLeakCoefficient(T,P)
            rxn = self.pathReactions[0]
            if rxn.products is self.source:
                assert len(rxn.reactants) == 1
                maxSpecies = rxn.reactants[0]
            else:
                assert len(rxn.products) == 1
                maxSpecies = rxn.products[0]
        else:
            for rxn in self.netReactions:
                if len(rxn.products) == 1 and rxn.products[0] not in self.explored:
                    k = rxn.getRateCoefficient(T,P)
                    if maxSpecies is None or k > maxK:
                        maxSpecies = rxn.products[0]
                        maxK = k

        # Make sure we've identified a species
        if maxSpecies is None:
            raise UnirxnNetworkException('No unimolecular isomers left to explore!')
        # Return the species
        return maxSpecies

    def getLeakBranchingRatios(self, T, P):
        """
        Return a dict with the unexplored isomers in the partial network as the
        keys and the fraction of the total leak coefficient as the values.
        """
        ratios = {}
        if len(self.netReactions) == 0 and len(self.pathReactions) == 1:
            rxn = self.pathReactions[0]
            if rxn.products is self.source:
                assert len(rxn.reactants) == 1
                ratios[rxn.reactants[0]] = 1.0
            else:
                assert len(rxn.products) == 1
                ratios[rxn.products[0]] = 1.0
        else:
            for rxn in self.netReactions:
                if len(rxn.products) == 1 and rxn.products[0] not in self.explored:
                    ratios[rxn.products[0]] = rxn.getRateCoefficient(T,P)

        kleak = sum(ratios.values())
        for spec in ratios:
            ratios[spec] /= kleak

        return ratios

    def exploreIsomer(self, isomer, reactionModel, database):
        """
        Explore a previously-unexplored unimolecular `isomer` in this partial
        network using the provided core-edge reaction model `reactionModel`,
        returning the new reactions and new species.
        """
        assert isomer not in self.explored
        assert isomer not in self.isomers
        assert isomer not in self.source

        logging.info('Exploring isomer {0} in pressure-dependent network #{1:d}'.format(isomer, self.index))
        self.explored.append(isomer)
        # Find reactions involving the found species as unimolecular
        # reactant or product (e.g. A <---> products)
        newReactionList = reactionModel.react(database, isomer)
        # Don't find reactions involving the new species as bimolecular
        # reactants or products with itself (e.g. A + A <---> products)
        # Don't find reactions involving the new species as bimolecular
        # reactants or products with other core species (e.g. A + B <---> products)

        return newReactionList

    def addPathReaction(self, newReaction, newSpecies):
        """
        Add a path reaction to the network. If the path reaction already exists,
        no action is taken.
        """
        # Add this reaction to that network if not already present
        found = False
        for rxn in self.pathReactions:
            if newReaction.reactants == rxn.reactants and newReaction.products == rxn.products:
                found = True
                break
            elif newReaction.products == rxn.reactants and newReaction.reactants == rxn.products:
                found = True
                break
        if not found:
            self.pathReactions.append(newReaction)
            self.invalidate()

    def merge(self, other):
        """
        Merge the partial network `other` into this network.
        """
        # Make sure the two partial networks have the same source configuration
        assert self.source == other.source

        # Merge isomers
        for isomer in other.isomers:
            if isomer not in self.isomers:
                self.isomers.append(isomer)
            else:
                # If the isomer is explored in either network, then it is explored in the combined network
                if isomer in other.explored and not isomer in self.explored:
                    self.explored.append(isomer)
        # Merge reactants
        for reactants in other.reactants:
            if reactants not in self.reactants:
                self.reactants.append(reactants)
        # Merge products
        for products in other.products:
            if isomer not in self.isomers:
                self.isomers.append(isomer)

        # Merge path reactions
        for reaction in other.pathReactions:
            found = False
            for rxn in self.pathReactions:
                if reaction.isEquivalent(rxn):
                    found = True
                    break
            if not found:
                self.pathReactions.append(reaction)

        # Also merge net reactions (so that when we update the network in the
        # future, we update the existing net reactions rather than making new ones)
        # Q: What to do when a net reaction exists in both networks being merged?
        for reaction in other.netReactions:
            found = False
            for rxn in self.netReactions:
                if reaction.isEquivalent(rxn):
                    found = True
                    break
            if not found:
                self.netReactions.append(reaction)

        # Mark this network as invalid
        self.valid = False

    def updateConfigurations(self):
        """
        Regenerate the :math:`k(T,P)` values for this partial network if the
        network is marked as invalid.
        """

        # Figure out which configurations are isomers, reactant channels, and product channels
        self.products = []
        if len(self.source) == 1:
            # The source is a unimolecular isomer
            self.isomers = self.explored[:]
            if self.source[0] not in self.isomers: self.isomers.insert(0, self.source[0])
            self.reactants = []
        else:
            # The source is a bimolecular reactant channel
            self.source.sort()
            self.isomers = self.explored[:]
            self.reactants = [self.source]
        for rxn in self.pathReactions:
            # Sort bimolecular configurations so that we always encounter them in the
            # same order
            # The actual order doesn't matter, as long as it is consistent
            rxn.reactants.sort()
            rxn.products.sort()
            # All reactants and products not already assigned as an isomer or
            # a reactant channel are product channels
            if rxn.reactants not in self.products:
                if len(rxn.reactants) == 1 and rxn.reactants[0] not in self.isomers:
                    self.products.append(rxn.reactants)
                elif len(rxn.reactants) > 1 and rxn.reactants not in self.reactants:
                    self.products.append(rxn.reactants)
            if rxn.products not in self.products:
                if len(rxn.products) == 1 and rxn.products[0] not in self.isomers:
                    self.products.append(rxn.products)
                elif len(rxn.products) > 1 and rxn.products not in self.reactants:
                    self.products.append(rxn.products)

    def update(self, reactionModel, database):
        """
        Regenerate the :math:`k(T,P)` values for this partial network if the
        network is marked as invalid.
        """
        from rmgpy.kinetics import Arrhenius, KineticsData
        from rmgpy.measure.collision import SingleExponentialDownModel
        from rmgpy.measure.reaction import fitInterpolationModel
        import rmgpy.measure.settings
        import rmgpy.measure.output

        # Get the parameters for the pressure dependence calculation
        method, Tmin, Tmax, Tlist, Pmin, Pmax, Plist, grainSize, numGrains, model = settings.pressureDependence

        # Figure out which configurations are isomers, reactant channels, and product channels
        self.updateConfigurations()

        # Make sure we have high-P kinetics for all path reactions
        for rxn in self.pathReactions:
            if rxn.kinetics is None and rxn.reverse.kinetics is None:
                raise PressureDependenceError('Path reaction {0} with no high-pressure-limit kinetics encountered in PDepNetwork #{1:d}.'.format(rxn, self.index))
            elif rxn.kinetics is not None and rxn.kinetics.isPressureDependent():
                raise PressureDependenceError('Pressure-dependent kinetics encountered for path reaction {0} in PDepNetwork #{1:d}.'.format(rxn, self.index))
        
        # Do nothing if the network is already valid
        if self.valid: return
        # Do nothing if there are no explored wells
        if len(self.explored) == 0 and len(self.source) > 1: return

        # Generate states data for unimolecular isomers and reactants if necessary
        for spec in self.isomers:
            if spec.states is None: spec.generateStatesData(database)
        for reactants in self.reactants:
            for spec in reactants:
                if spec.states is None: spec.generateStatesData(database)

        # Determine transition state energies on potential energy surface
        # In the absence of any better information, we simply set it to
        # be the reactant ground-state energy + the activation energy
        # Note that we need Arrhenius kinetics in order to do this
        for rxn in self.pathReactions:
            if rxn.kinetics is None:
                raise Exception('Path reaction "{0}" in PDepNetwork #{1:d} has no kinetics!'.format(rxn, self.index))
            elif isinstance(rxn.kinetics, KineticsData):
                if len(rxn.reactants) == 1:
                    kunits = 's^-1'
                elif len(rxn.reactants) == 2:
                    kunits = 'm^3/(mol*s)'
                elif len(rxn.reactants) == 3:
                    kunits = 'm^6/(mol^2*s)'
                else:
                    kunits = ''
                rxn.kinetics = Arrhenius().fitToData(Tlist=rxn.kinetics.Tdata.values, klist=rxn.kinetics.kdata.values, kunits=kunits)
            elif not isinstance(rxn.kinetics, Arrhenius):
                raise Exception('Path reaction "{0}" in PDepNetwork #{1:d} has invalid kinetics type "{2}".'.format(rxn, rxn.kinetics.__class__))
            rxn.transitionState = rmgpy.species.TransitionState(
                E0=((sum([spec.E0.value for spec in rxn.reactants]) + rxn.kinetics.Ea.value)/1000.,"kJ/mol"),
            )

        # Set collision model
        bathGas = [spec for spec in reactionModel.core.species if not spec.reactive]
        self.bathGas = {}
        for spec in bathGas:
            self.bathGas[spec] = 1.0 / len(bathGas)
        self.collisionModel = SingleExponentialDownModel(alpha0=4.86 * 4184)

        # Save input file
        rmgpy.measure.output.writeInput(os.path.join(settings.outputDirectory, 'pdep', 'network{0:d}_{1:d}.py'.format(self.index, len(self.isomers))),
            self, Tlist, Plist, (grainSize, numGrains), method, model)

        self.printSummary(level=logging.INFO)

        # Automatically choose a suitable set of energy grains if they were not
        # explicitly specified in the input file
        Elist = self.autoGenerateEnergyGrains(Tmax=Tmax, grainSize=grainSize, Ngrains=numGrains)
        
        # Calculate the rate coefficients
        K, p0 = self.calculateRateCoefficients(Tlist, Plist, Elist, method)

        # Generate PDepReaction objects
        configurations = []
        configurations.extend([[isom] for isom in self.isomers])
        configurations.extend([reactants for reactants in self.reactants])
        configurations.extend([products for products in self.products])
        j = configurations.index(self.source)

        for i in range(K.shape[2]):
            if i != j:
                # Find the path reaction
                netReaction = None
                for r in self.netReactions:
                    if r.hasTemplate(configurations[j], configurations[i]):
                        netReaction = r
                # If net reaction does not already exist, make a new one
                if netReaction is None:
                    netReaction = PDepReaction(
                        reactants=configurations[j],
                        products=configurations[i],
                        network=self,
                        kinetics=None
                    )
                    netReaction = reactionModel.makeNewPDepReaction(netReaction)
                    self.netReactions.append(netReaction)

                    # Place the net reaction in the core or edge if necessary
                    # Note that leak reactions are not placed in the edge
                    if all([s in reactionModel.core.species for s in netReaction.reactants]) and all([s in reactionModel.core.species for s in netReaction.products]):
                        reactionModel.addReactionToCore(netReaction)
                    else:
                        reactionModel.addReactionToEdge(netReaction)

                # Set/update the net reaction kinetics using interpolation model
                netReaction.kinetics = fitInterpolationModel(netReaction, Tlist, Plist, K[:,:,i,j], model, Tmin, Tmax, Pmin, Pmax)

        # We're done processing this network, so mark it as valid
        self.valid = True

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

class ReactionModel:
    """
    Represent a generic reaction model. A reaction model consists of `species`,
    a list of species, and `reactions`, a list of reactions.
    """

    def __init__(self, species=None, reactions=None):
        self.species = species or []
        self.reactions = reactions or []

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

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
    `absoluteTolerance`        The absolute tolerance used in the ODE/DAE solver
    `relativeTolerance`        The relative tolerance used in the ODE/DAE solver
    `fluxToleranceKeepInEdge`  The relative species flux below which species are discarded from the edge
    `fluxToleranceMoveToCore`  The relative species flux above which species are moved from the edge to the core
    `fluxToleranceInterrupt`   The relative species flux above which the simulation will halt
    `maximumEdgeSpecies`       The maximum number of edge species allowed at any time
    `termination`              A list of termination targets (i.e :class:`TerminationTime` and :class:`TerminationConversion` objects)
    `unirxnNetworks`           A list of unimolecular reaction networks (:class:`unirxn.network.Network` objects)
    `networkCount`             A counter for the number of unirxn networks created
    =========================  ==============================================================


    """

    def __init__(self, core=None, edge=None):
        if core is None:
            self.core = ReactionModel()
        else:
            self.core = core
        if edge is None:
            self.edge = ReactionModel()
        else:
            self.edge = edge
        # 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.fluxToleranceKeepInEdge = 0.0
        self.fluxToleranceMoveToCore = 1.0
        self.fluxToleranceInterrupt = 1.0
        self.absoluteTolerance = 1.0e-8
        self.relativeTolerance = 1.0e-4
        self.maximumEdgeSpecies = 1000000
        self.termination = []
        self.unirxnNetworks = []
        self.networkCount = 0
        self.speciesDict = {}
        self.reactionDict = {}
        self.speciesCache = [None for i in range(4)]
        self.speciesCounter = 0
        self.reactionCounter = 0
        self.newSpeciesList = []
        self.newReactionList = []

    def checkForExistingSpecies(self, molecule):
        """
        Check to see if an existing species contains the same
        :class:`structure.Structure` as `structure`. Returns ``True`` or
        ``False`` and the matched species (if found, or ``None`` if not).
        """

        # First check cache and return if species is found
        for i, spec in enumerate(self.speciesCache):
            if spec is not None:
                for mol in spec.molecule:
                    if molecule.isIsomorphic(mol):
                        self.speciesCache.pop(i)
                        self.speciesCache.insert(0, spec)
                        return True, spec

        # Return an existing species if a match is found
        formula = molecule.getFormula()
        try:
            for spec in self.speciesDict[formula]:
                for mol in spec.molecule:
                    if molecule.isIsomorphic(mol):
                        self.speciesCache.pop()
                        self.speciesCache.insert(0, spec)
                        return True, spec
        except KeyError: pass
        
        # At this point we can conclude that the structure does not exist
        return False, None

    def makeNewSpecies(self, object, label='', reactive=True, checkForExisting=True):
        """
        Formally create a new species from the specified `object`, which can be
        either a :class:`Molecule` object or an :class:`rmgpy.species.Species`
        object.
        """

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

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

        # Check that the structure is not forbidden

        # If we're here then we're ready to make the new species
        if label == '': label = molecule.getFormula()
        logging.debug('Creating new species {0}'.format(label))
        spec = Species(index=self.speciesCounter+1, label=label, molecule=[molecule], reactive=reactive)
        spec.coreSizeAtCreation = len(self.core.species)
        spec.generateResonanceIsomers()
        spec.molecularWeight = Quantity(spec.molecule[0].getMolecularWeight()*1000.,"g/mol")
        spec.generateLennardJonesParameters()
        formula = molecule.getFormula()
        if formula in self.speciesDict:
            self.speciesDict[formula].append(spec)
        else:
            self.speciesDict[formula] = [spec]

        # Store hydrogens implicitly to conserve memory and speed up isomorphism
        for mol in spec.molecule:
            #mol.updateConnectivityValues()
            mol.makeHydrogensImplicit()

        self.speciesCounter += 1

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

        return spec, True

    def checkForExistingReaction(self, rxn):
        """
        Check to see if an existing reaction has the same reactants, products, and
        family as `rxn`. Returns :data:`True` or :data:`False` and the matched
        reaction (if found).
        """

        # Make sure the reactant and product lists are sorted before performing the check
        rxn.reactants.sort()
        rxn.products.sort()

        # Get the short-list of reactions with the same family, reactant1 and reactant2
        r1 = rxn.reactants[0]
        if len(rxn.reactants)==1: r2 = None
        else: r2 = rxn.reactants[1]
        family = rxn.getSource()
        try:
            my_reactionList = self.reactionDict[family][r1][r2]
        except KeyError: # no such short-list: must be new, unless in seed.
            my_reactionList = []

        # Now use short-list to check for matches. All should be in same forward direction.
        for rxn0 in my_reactionList:
            if (rxn0.reactants == rxn.reactants and rxn0.products == rxn.products):
                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
        for family0 in self.reactionDict:
            if isinstance(family, KineticsLibrary) and family0 != family:

                # First check seed short-list in forward direction
                try:
                    my_reactionList = self.reactionDict[family0][r1][r2]
                except KeyError:
                    my_reactionList = []
                for rxn0 in my_reactionList:
                    if (rxn0.reactants == rxn.reactants and rxn0.products == rxn.products) or \
                        (rxn0.reactants == rxn.products and rxn0.products == rxn.reactants):
                        return True, rxn0
                # Now get the seed short-list of the reverse reaction
                r1 = rxn.products[0]
                if len(rxn.products)==1: r2 = None
                else: r2 = rxn.products[1]
                try:
                    my_reactionList = self.reactionDict[family0][r1][r2]
                except KeyError:
                    my_reactionList = []
                for rxn0 in my_reactionList:
                    if (rxn0.reactants == rxn.reactants and rxn0.products == rxn.products) or \
                        (rxn0.reactants == rxn.products and rxn0.products == rxn.reactants):
                        return True, rxn0

        return False, None

    def makeNewReaction(self, forward, checkExisting=True):
        """
        Make a new reaction given a :class:`Reaction` object `forward`. The kinetics
        of the reaction are estimated and 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.
        """

        # Determine the proper species objects for all reactants and products
        forward.reactants = [self.makeNewSpecies(reactant)[0] for reactant in forward.reactants]
        forward.products  = [self.makeNewSpecies(product)[0]  for product  in forward.products ]

        if checkExisting:
            found, rxn = self.checkForExistingReaction(forward)
            if found: return rxn, False

        # Note in the log
        if isinstance(rxn, TemplateReaction):
            logging.debug('Creating new {0} reaction {1}'.format(forward.family.label, forward))
        else:
            logging.debug('Creating new library reaction {0}'.format(forward))
        
        # Add to the global dict/list of existing reactions (a list broken down by family, r1, r2)
        # identify r1 and r2
        r1 = forward.reactants[0]
        r2 = None if len(forward.reactants) == 1 else forward.reactants[1]
        family = forward.getSource()
        # make dictionary entries if necessary
        if family not in self.reactionDict:
            self.reactionDict[family] = {}
        if not self.reactionDict[family].has_key(r1):
            self.reactionDict[family][r1] = dict()
        if not self.reactionDict[family][r1].has_key(r2):
            self.reactionDict[family][r1][r2] = list()
        # store this reaction at the top of the relevant short-list
        self.reactionDict[family][r1][r2].insert(0, forward)

        forward.index = self.reactionCounter + 1
        self.reactionCounter += 1

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

        # Return newly created reaction
        return forward, True

    def makeNewPDepReaction(self, forward):
        """
        Make a new pressure-dependent reaction based on a list of `reactants` and a
        list of `products`. The reaction belongs to the specified `network` and
        has pressure-dependent kinetics given by `kinetics`.

        No checking for existing reactions is made here. The returned PDepReaction
        object is not added to the global list of reactions, as that is intended
        to represent only the high-pressure-limit set. The reactionCounter is
        incremented, however, since the returned reaction can and will exist in
        the model edge and/or core.
        """

        # Don't create reverse reaction: all such reactions are treated as irreversible
        # The reverse direction will come from a different partial network
        # Note that this isn't guaranteed to satisfy thermodynamics (but will probably be close)
        forward.reverse = None
        forward.reversible = False

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

        return forward

    def react(self, database, speciesA, speciesB=None):
        reactionList = []
        if speciesB is None:
            for moleculeA in speciesA.molecule:
                reactionList.extend(database.kinetics.generateReactions([moleculeA]))
                moleculeA.clearLabeledAtoms()
        else:
            for moleculeA in speciesA.molecule:
                for moleculeB in speciesB.molecule:
                    reactionList.extend(database.kinetics.generateReactions([moleculeA, moleculeB]))
                    moleculeA.clearLabeledAtoms()
                    moleculeB.clearLabeledAtoms()
        return reactionList

    def enlarge(self, newObject, database):
        """
        Enlarge a reaction model by processing `newObject`. If `newObject` is a
        :class:`rmg.species.Species` object, then the species is moved from
        the edge to the core and reactions generated for that species, reacting
        with itself and with all other species in the model core. If `newObject`
        is a :class:`rmg.unirxn.network.Network` object, then reactions are
        generated for the species in the network with the largest leak flux.
        """

        self.newReactionList = []; self.newSpeciesList = []
        newReactions = []

        numOldCoreSpecies = len(self.core.species)
        numOldCoreReactions = len(self.core.reactions)
        numOldEdgeSpecies = len(self.edge.species)
        numOldEdgeReactions = len(self.edge.reactions)

        pdepNetwork = None

        if isinstance(newObject, Species):

            newSpecies = newObject

            if not newSpecies.reactive:
                logging.info('NOT generating reactions for unreactive species {0}'.format(newSpecies))
            else:
                logging.info('Adding species {0} to model core'.format(newSpecies))
                # Find reactions involving the new species as unimolecular reactant
                # or product (e.g. A <---> products)
                newReactions.extend(self.react(database, newSpecies))
                # Find reactions involving the new species as bimolecular reactants
                # or products with other core species (e.g. A + B <---> products)
                for coreSpecies in self.core.species:
                    if coreSpecies.reactive:
                        newReactions.extend(self.react(database, newSpecies, coreSpecies))
                # Find reactions involving the new species as bimolecular reactants
                # or products with itself (e.g. A + A <---> products)
                newReactions.extend(self.react(database, newSpecies, newSpecies))

            # Add new species
            self.addSpeciesToCore(newSpecies)

            # Process the new reactions
            self.processNewReactions(newReactions, newSpecies, pdepNetwork)

        elif isinstance(newObject, tuple) and isinstance(newObject[0], PDepNetwork) and settings.pressureDependence:

            pdepNetwork, newSpecies = newObject
            newReactions.extend(network.exploreIsomer(species, self, database))
            self.processNewReactions(newReactions, newSpecies, pdepNetwork)
<commit_msg>Fixed bug in exploring of partial pdep network with high leak flux.<commit_after>#!/usr/bin/python
# -*- coding: utf-8 -*-

################################################################################
#
#	RMG - Reaction Mechanism Generator
#
#	Copyright (c) 2002-2009 Prof. <NAME> (<EMAIL>) and the
#	RMG Team (<EMAIL>)
#
#	Permission is hereby granted, free of charge, to any person obtaining a
#	copy of this software and associated documentation files (the 'Software'),
#	to deal in the Software without restriction, including without limitation
#	the rights to use, copy, modify, merge, publish, distribute, sublicense,
#	and/or sell copies of the Software, and to permit persons to whom the
#	Software is furnished to do so, subject to the following conditions:
#
#	The above copyright notice and this permission notice shall be included in
#	all copies or substantial portions of the Software.
#
#	THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
#	FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
#	DEALINGS IN THE SOFTWARE.
#
################################################################################

"""
Contains classes for working with the reaction model generated by RMG.
"""

import logging
import math
import numpy
import os
import os.path

from rmgpy.quantity import Quantity, constants
import rmgpy.species
from rmgpy.thermo import Wilhoit, MultiNASA

from rmgpy.data.thermo import *
from rmgpy.data.kinetics import *
from rmgpy.data.states import *

import rmgpy.measure.network

import settings

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

class Species(rmgpy.species.Species):

    def __init__(self, index=-1, label='', thermo=None, states=None, molecule=None, E0=0.0, lennardJones=None, molecularWeight=0.0, reactive=True):
        rmgpy.species.Species.__init__(self, index, label, thermo, states, molecule, E0, lennardJones, molecularWeight, reactive)
        self.coreSizeAtCreation = 0

    def __reduce__(self):
        """
        A helper function used when pickling an object.
        """
        d = {
            'coreSizeAtCreation': self.coreSizeAtCreation,
        }
        return (Species, (self.index, self.label, self.thermo, self.states, self.molecule, self.E0, self.lennardJones, self.molecularWeight, self.reactive), d)

    def __setstate__(self, d):
        """
        A helper function used when unpickling an object.
        """
        self.coreSizeAtCreation = d['coreSizeAtCreation']

    def generateThermoData(self, database, thermoClass=MultiNASA):
        """
        Generate thermodynamic data for the species using the thermo database.

        Generates the thermo data for each structure (resonance isomer),
        picks that with lowest H298 value, and saves it to `self.thermoData`.
        """

        # Ensure molecules are using explicit hydrogens
        implicitH = [mol.implicitHydrogens for mol in self.molecule]
        for molecule in self.molecule:
            molecule.makeHydrogensExplicit()
        
        thermo = []
        for molecule in self.molecule:
            molecule.clearLabeledAtoms()
            molecule.updateAtomTypes()
            tdata = database.thermo.getThermoData(molecule)
            thermo.append(tdata)

        H298 = numpy.array([t.getEnthalpy(298.) for t in thermo])
        indices = H298.argsort()

        # If multiple resonance isomers are present, use the thermo data of
        # the most stable isomer (i.e. one with lowest enthalpy of formation)
        # as the thermo data of the species
        self.thermo = thermo[indices[0]]

        # Sort the structures in order of decreasing stability
        self.molecule = [self.molecule[ind] for ind in indices]
        implicitH = [implicitH[ind] for ind in indices]

        # Convert to desired thermo class
        thermo0 = self.thermo
        linear = self.molecule[0].isLinear()
        nRotors = self.molecule[0].countInternalRotors()
        nFreq = 3 * len(self.molecule[0].atoms) - (5 if linear else 6) - nRotors
        self.thermo = convertThermoModel(self.thermo, Wilhoit, linear=linear, nFreq=nFreq, nRotors=nRotors)
        self.E0 = Quantity(self.thermo.getEnthalpy(1.0)/1000.0,"kJ/mol")
        self.thermo = convertThermoModel(self.thermo, thermoClass, Tmin=100.0, Tmax=5000.0, Tint=1000.0)
        if self.thermo.__class__ != thermo0.__class__:
            # Compute RMS error of overall transformation
            Tlist = numpy.array([300.0, 400.0, 500.0, 600.0, 800.0, 1000.0, 1500.0], numpy.float64)
            err = math.sqrt(numpy.sum((self.thermo.getHeatCapacities(Tlist) - thermo0.getHeatCapacities(Tlist))**2))/constants.R/len(Tlist)
            logging.log(logging.WARNING if err > 0.1 else 0, 'Average RMS error in heat capacity fit to {0} = {1:g}*R'.format(self, err))

        # Restore implicit hydrogens if necessary
        for implicit, molecule in zip(implicitH, self.molecule):
            if implicit: molecule.makeHydrogensImplicit()

        return self.thermo

    def generateStatesData(self, database):
        """
        Generate molecular degree of freedom data for the species. You must
        have already provided a thermodynamics model using e.g.
        :meth:`generateThermoData()`.
        """
        if not self.thermo:
            raise Exception("Unable to determine states model for species {0}: No thermodynamics model found.".format(self))
        molecule = self.molecule[0]
        implicitH = molecule.implicitHydrogens
        molecule.makeHydrogensExplicit()
        self.states = database.states.getStatesData(molecule, self.thermo)
        if implicitH: molecule.makeHydrogensImplicit()

    def generateLennardJonesParameters(self):
        """
        Generate the Lennard-Jones parameters for the species. This "algorithm"
        is *very* much in need of improvement.
        """

        count = sum([1 for atom in self.molecule[0].vertices if atom.isNonHydrogen()])
        self.lennardJones = rmgpy.species.LennardJones()

        if count == 1:
            self.lennardJones.sigma = Quantity(3.758e-10,"m")
            self.lennardJones.epsilon = Quantity(148.6 * constants.kB,"J")
        elif count == 2:
            self.lennardJones.sigma = Quantity(4.443e-10,"m")
            self.lennardJones.epsilon = Quantity(110.7 * constants.kB,"J")
        elif count == 3:
            self.lennardJones.sigma = Quantity(5.118e-10,"m")
            self.lennardJones.epsilon = Quantity(237.1 * constants.kB,"J")
        elif count == 4:
            self.lennardJones.sigma = Quantity(4.687e-10,"m")
            self.lennardJones.epsilon = Quantity(531.4 * constants.kB,"J")
        elif count == 5:
            self.lennardJones.sigma = Quantity(5.784e-10,"m")
            self.lennardJones.epsilon = Quantity(341.1 * constants.kB,"J")
        else:
            self.lennardJones.sigma = Quantity(5.949e-10,"m")
            self.lennardJones.epsilon = Quantity(399.3 * constants.kB,"J")

class PDepReaction(rmgpy.reaction.Reaction):

    def __init__(self, index=-1, reactants=None, products=None, network=None, kinetics=None, reversible=True, transitionState=None, thirdBody=False):
        rmgpy.reaction.Reaction.__init__(self, index, reactants, products, kinetics, reversible, transitionState, thirdBody)
        self.network = network

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

class PressureDependenceError(Exception):
    """
    An exception class to use when an error involving pressure dependence is
    encountered. Pass a string describing the circumstances of the exceptional
    behavior.
    """
    pass

class PDepNetwork(rmgpy.measure.network.Network):
    """
    A representation of a *partial* unimolecular reaction network. Each partial
    network has a single `source` isomer or reactant channel, and is responsible
    only for :math:`k(T,P)` values for net reactions with source as the
    reactant. Multiple partial networks can have the same source, but networks
    with the same source and any explored isomers must be combined.

    =================== ======================= ================================
    Attribute           Type                    Description
    =================== ======================= ================================
    `source`            ``list``                The isomer or reactant channel that acts as the source
    `explored`          ``list``                A list of the unimolecular isomers whose reactions have been fully explored
    =================== ======================= ================================

    """

    def __init__(self, index=-1, source=None):
        rmgpy.measure.network.Network.__init__(self, index=index)
        self.source = source
        self.explored = []

    def getLeakCoefficient(self, T, P):
        """
        Return the pressure-dependent rate coefficient :math:`k(T,P)` describing
        the total rate of "leak" from this network. This is defined as the sum
        of the :math:`k(T,P)` values for all net reactions to nonexplored
        unimolecular isomers.
        """
        k = 0.0
        if len(self.netReactions) == 0 and len(self.pathReactions) == 1:
            # The network is of the form A + B -> C* (with C* nonincluded)
            # For this special case we use the high-pressure limit k(T) to
            # ensure that we're estimating the total leak flux
            rxn = self.pathReactions[0]
            if rxn.kinetics is None:
                if rxn.reverse.kinetics is not None:
                    rxn = rxn.reverse
                else:
                    raise PressureDependenceError('Path reaction {0} with no high-pressure-limit kinetics encountered in PDepNetwork #{1:d} while evaluating leak flux.'.format(rxn, self.index))
            if rxn.products is self.source:
                k = rxn.getRateCoefficient(T,P) / rxn.getEquilibriumConstant(T)
            else:
                k = rxn.getRateCoefficient(T,P)
        else:
            # The network has at least one included isomer, so we can calculate
            # the leak flux normally
            for rxn in self.netReactions:
                if len(rxn.products) == 1 and rxn.products[0] not in self.explored:
                    k += rxn.getRateCoefficient(T,P)
        return k

    def getMaximumLeakSpecies(self, T, P):
        """
        Get the unexplored (unimolecular) isomer with the maximum leak flux.
        Note that the leak rate coefficients vary with temperature and
        pressure, so you must provide these in order to get a meaningful result.
        """
        # Choose species with maximum leak flux
        maxK = 0.0; maxSpecies = None
        if len(self.netReactions) == 0 and len(self.pathReactions) == 1:
            maxK = self.getLeakCoefficient(T,P)
            rxn = self.pathReactions[0]
            if rxn.products is self.source:
                assert len(rxn.reactants) == 1
                maxSpecies = rxn.reactants[0]
            else:
                assert len(rxn.products) == 1
                maxSpecies = rxn.products[0]
        else:
            for rxn in self.netReactions:
                if len(rxn.products) == 1 and rxn.products[0] not in self.explored:
                    k = rxn.getRateCoefficient(T,P)
                    if maxSpecies is None or k > maxK:
                        maxSpecies = rxn.products[0]
                        maxK = k

        # Make sure we've identified a species
        if maxSpecies is None:
            raise UnirxnNetworkException('No unimolecular isomers left to explore!')
        # Return the species
        return maxSpecies

    def getLeakBranchingRatios(self, T, P):
        """
        Return a dict with the unexplored isomers in the partial network as the
        keys and the fraction of the total leak coefficient as the values.
        """
        ratios = {}
        if len(self.netReactions) == 0 and len(self.pathReactions) == 1:
            rxn = self.pathReactions[0]
            if rxn.products is self.source:
                assert len(rxn.reactants) == 1
                ratios[rxn.reactants[0]] = 1.0
            else:
                assert len(rxn.products) == 1
                ratios[rxn.products[0]] = 1.0
        else:
            for rxn in self.netReactions:
                if len(rxn.products) == 1 and rxn.products[0] not in self.explored:
                    ratios[rxn.products[0]] = rxn.getRateCoefficient(T,P)

        kleak = sum(ratios.values())
        for spec in ratios:
            ratios[spec] /= kleak

        return ratios

    def exploreIsomer(self, isomer, reactionModel, database):
        """
        Explore a previously-unexplored unimolecular `isomer` in this partial
        network using the provided core-edge reaction model `reactionModel`,
        returning the new reactions and new species.
        """
        assert isomer not in self.explored
        assert isomer not in self.isomers
        assert isomer not in self.source

        logging.info('Exploring isomer {0} in pressure-dependent network #{1:d}'.format(isomer, self.index))
        self.explored.append(isomer)
        # Find reactions involving the found species as unimolecular
        # reactant or product (e.g. A <---> products)
        newReactionList = reactionModel.react(database, isomer)
        # Don't find reactions involving the new species as bimolecular
        # reactants or products with itself (e.g. A + A <---> products)
        # Don't find reactions involving the new species as bimolecular
        # reactants or products with other core species (e.g. A + B <---> products)

        return newReactionList

    def addPathReaction(self, newReaction, newSpecies):
        """
        Add a path reaction to the network. If the path reaction already exists,
        no action is taken.
        """
        # Add this reaction to that network if not already present
        found = False
        for rxn in self.pathReactions:
            if newReaction.reactants == rxn.reactants and newReaction.products == rxn.products:
                found = True
                break
            elif newReaction.products == rxn.reactants and newReaction.reactants == rxn.products:
                found = True
                break
        if not found:
            self.pathReactions.append(newReaction)
            self.invalidate()

    def merge(self, other):
        """
        Merge the partial network `other` into this network.
        """
        # Make sure the two partial networks have the same source configuration
        assert self.source == other.source

        # Merge isomers
        for isomer in other.isomers:
            if isomer not in self.isomers:
                self.isomers.append(isomer)
            else:
                # If the isomer is explored in either network, then it is explored in the combined network
                if isomer in other.explored and not isomer in self.explored:
                    self.explored.append(isomer)
        # Merge reactants
        for reactants in other.reactants:
            if reactants not in self.reactants:
                self.reactants.append(reactants)
        # Merge products
        for products in other.products:
            if isomer not in self.isomers:
                self.isomers.append(isomer)

        # Merge path reactions
        for reaction in other.pathReactions:
            found = False
            for rxn in self.pathReactions:
                if reaction.isEquivalent(rxn):
                    found = True
                    break
            if not found:
                self.pathReactions.append(reaction)

        # Also merge net reactions (so that when we update the network in the
        # future, we update the existing net reactions rather than making new ones)
        # Q: What to do when a net reaction exists in both networks being merged?
        for reaction in other.netReactions:
            found = False
            for rxn in self.netReactions:
                if reaction.isEquivalent(rxn):
                    found = True
                    break
            if not found:
                self.netReactions.append(reaction)

        # Mark this network as invalid
        self.valid = False

    def updateConfigurations(self):
        """
        Regenerate the :math:`k(T,P)` values for this partial network if the
        network is marked as invalid.
        """

        # Figure out which configurations are isomers, reactant channels, and product channels
        self.products = []
        if len(self.source) == 1:
            # The source is a unimolecular isomer
            self.isomers = self.explored[:]
            if self.source[0] not in self.isomers: self.isomers.insert(0, self.source[0])
            self.reactants = []
        else:
            # The source is a bimolecular reactant channel
            self.source.sort()
            self.isomers = self.explored[:]
            self.reactants = [self.source]
        for rxn in self.pathReactions:
            # Sort bimolecular configurations so that we always encounter them in the
            # same order
            # The actual order doesn't matter, as long as it is consistent
            rxn.reactants.sort()
            rxn.products.sort()
            # All reactants and products not already assigned as an isomer or
            # a reactant channel are product channels
            if rxn.reactants not in self.products:
                if len(rxn.reactants) == 1 and rxn.reactants[0] not in self.isomers:
                    self.products.append(rxn.reactants)
                elif len(rxn.reactants) > 1 and rxn.reactants not in self.reactants:
                    self.products.append(rxn.reactants)
            if rxn.products not in self.products:
                if len(rxn.products) == 1 and rxn.products[0] not in self.isomers:
                    self.products.append(rxn.products)
                elif len(rxn.products) > 1 and rxn.products not in self.reactants:
                    self.products.append(rxn.products)

    def update(self, reactionModel, database):
        """
        Regenerate the :math:`k(T,P)` values for this partial network if the
        network is marked as invalid.
        """
        from rmgpy.kinetics import Arrhenius, KineticsData
        from rmgpy.measure.collision import SingleExponentialDownModel
        from rmgpy.measure.reaction import fitInterpolationModel
        import rmgpy.measure.settings
        import rmgpy.measure.output

        # Get the parameters for the pressure dependence calculation
        method, Tmin, Tmax, Tlist, Pmin, Pmax, Plist, grainSize, numGrains, model = settings.pressureDependence

        # Figure out which configurations are isomers, reactant channels, and product channels
        self.updateConfigurations()

        # Make sure we have high-P kinetics for all path reactions
        for rxn in self.pathReactions:
            if rxn.kinetics is None and rxn.reverse.kinetics is None:
                raise PressureDependenceError('Path reaction {0} with no high-pressure-limit kinetics encountered in PDepNetwork #{1:d}.'.format(rxn, self.index))
            elif rxn.kinetics is not None and rxn.kinetics.isPressureDependent():
                raise PressureDependenceError('Pressure-dependent kinetics encountered for path reaction {0} in PDepNetwork #{1:d}.'.format(rxn, self.index))
        
        # Do nothing if the network is already valid
        if self.valid: return
        # Do nothing if there are no explored wells
        if len(self.explored) == 0 and len(self.source) > 1: return

        # Generate states data for unimolecular isomers and reactants if necessary
        for spec in self.isomers:
            if spec.states is None: spec.generateStatesData(database)
        for reactants in self.reactants:
            for spec in reactants:
                if spec.states is None: spec.generateStatesData(database)

        # Determine transition state energies on potential energy surface
        # In the absence of any better information, we simply set it to
        # be the reactant ground-state energy + the activation energy
        # Note that we need Arrhenius kinetics in order to do this
        for rxn in self.pathReactions:
            if rxn.kinetics is None:
                raise Exception('Path reaction "{0}" in PDepNetwork #{1:d} has no kinetics!'.format(rxn, self.index))
            elif isinstance(rxn.kinetics, KineticsData):
                if len(rxn.reactants) == 1:
                    kunits = 's^-1'
                elif len(rxn.reactants) == 2:
                    kunits = 'm^3/(mol*s)'
                elif len(rxn.reactants) == 3:
                    kunits = 'm^6/(mol^2*s)'
                else:
                    kunits = ''
                rxn.kinetics = Arrhenius().fitToData(Tlist=rxn.kinetics.Tdata.values, klist=rxn.kinetics.kdata.values, kunits=kunits)
            elif not isinstance(rxn.kinetics, Arrhenius):
                raise Exception('Path reaction "{0}" in PDepNetwork #{1:d} has invalid kinetics type "{2}".'.format(rxn, rxn.kinetics.__class__))
            rxn.transitionState = rmgpy.species.TransitionState(
                E0=((sum([spec.E0.value for spec in rxn.reactants]) + rxn.kinetics.Ea.value)/1000.,"kJ/mol"),
            )

        # Set collision model
        bathGas = [spec for spec in reactionModel.core.species if not spec.reactive]
        self.bathGas = {}
        for spec in bathGas:
            self.bathGas[spec] = 1.0 / len(bathGas)
        self.collisionModel = SingleExponentialDownModel(alpha0=4.86 * 4184)

        # Save input file
        rmgpy.measure.output.writeInput(os.path.join(settings.outputDirectory, 'pdep', 'network{0:d}_{1:d}.py'.format(self.index, len(self.isomers))),
            self, Tlist, Plist, (grainSize, numGrains), method, model)

        self.printSummary(level=logging.INFO)

        # Automatically choose a suitable set of energy grains if they were not
        # explicitly specified in the input file
        Elist = self.autoGenerateEnergyGrains(Tmax=Tmax, grainSize=grainSize, Ngrains=numGrains)
        
        # Calculate the rate coefficients
        K, p0 = self.calculateRateCoefficients(Tlist, Plist, Elist, method)

        # Generate PDepReaction objects
        configurations = []
        configurations.extend([[isom] for isom in self.isomers])
        configurations.extend([reactants for reactants in self.reactants])
        configurations.extend([products for products in self.products])
        j = configurations.index(self.source)

        for i in range(K.shape[2]):
            if i != j:
                # Find the path reaction
                netReaction = None
                for r in self.netReactions:
                    if r.hasTemplate(configurations[j], configurations[i]):
                        netReaction = r
                # If net reaction does not already exist, make a new one
                if netReaction is None:
                    netReaction = PDepReaction(
                        reactants=configurations[j],
                        products=configurations[i],
                        network=self,
                        kinetics=None
                    )
                    netReaction = reactionModel.makeNewPDepReaction(netReaction)
                    self.netReactions.append(netReaction)

                    # Place the net reaction in the core or edge if necessary
                    # Note that leak reactions are not placed in the edge
                    if all([s in reactionModel.core.species for s in netReaction.reactants]) and all([s in reactionModel.core.species for s in netReaction.products]):
                        reactionModel.addReactionToCore(netReaction)
                    else:
                        reactionModel.addReactionToEdge(netReaction)

                # Set/update the net reaction kinetics using interpolation model
                netReaction.kinetics = fitInterpolationModel(netReaction, Tlist, Plist, K[:,:,i,j], model, Tmin, Tmax, Pmin, Pmax)

        # We're done processing this network, so mark it as valid
        self.valid = True

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

class ReactionModel:
    """
    Represent a generic reaction model. A reaction model consists of `species`,
    a list of species, and `reactions`, a list of reactions.
    """

    def __init__(self, species=None, reactions=None):
        self.species = species or []
        self.reactions = reactions or []

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

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
    `absoluteTolerance`        The absolute tolerance used in the ODE/DAE solver
    `relativeTolerance`        The relative tolerance used in the ODE/DAE solver
    `fluxToleranceKeepInEdge`  The relative species flux below which species are discarded from the edge
    `fluxToleranceMoveToCore`  The relative species flux above which species are moved from the edge to the core
    `fluxToleranceInterrupt`   The relative species flux above which the simulation will halt
    `maximumEdgeSpecies`       The maximum number of edge species allowed at any time
    `termination`              A list of termination targets (i.e :class:`TerminationTime` and :class:`TerminationConversion` objects)
    `unirxnNetworks`           A list of unimolecular reaction networks (:class:`unirxn.network.Network` objects)
    `networkCount`             A counter for the number of unirxn networks created
    =========================  ==============================================================


    """

    def __init__(self, core=None, edge=None):
        if core is None:
            self.core = ReactionModel()
        else:
            self.core = core
        if edge is None:
            self.edge = ReactionModel()
        else:
            self.edge = edge
        # 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.fluxToleranceKeepInEdge = 0.0
        self.fluxToleranceMoveToCore = 1.0
        self.fluxToleranceInterrupt = 1.0
        self.absoluteTolerance = 1.0e-8
        self.relativeTolerance = 1.0e-4
        self.maximumEdgeSpecies = 1000000
        self.termination = []
        self.unirxnNetworks = []
        self.networkCount = 0
        self.speciesDict = {}
        self.reactionDict = {}
        self.speciesCache = [None for i in range(4)]
        self.speciesCounter = 0
        self.reactionCounter = 0
        self.newSpeciesList = []
        self.newReactionList = []

    def checkForExistingSpecies(self, molecule):
        """
        Check to see if an existing species contains the same
        :class:`structure.Structure` as `structure`. Returns ``True`` or
        ``False`` and the matched species (if found, or ``None`` if not).
        """

        # First check cache and return if species is found
        for i, spec in enumerate(self.speciesCache):
            if spec is not None:
                for mol in spec.molecule:
                    if molecule.isIsomorphic(mol):
                        self.speciesCache.pop(i)
                        self.speciesCache.insert(0, spec)
                        return True, spec

        # Return an existing species if a match is found
        formula = molecule.getFormula()
        try:
            for spec in self.speciesDict[formula]:
                for mol in spec.molecule:
                    if molecule.isIsomorphic(mol):
                        self.speciesCache.pop()
                        self.speciesCache.insert(0, spec)
                        return True, spec
        except KeyError: pass
        
        # At this point we can conclude that the structure does not exist
        return False, None

    def makeNewSpecies(self, object, label='', reactive=True, checkForExisting=True):
        """
        Formally create a new species from the specified `object`, which can be
        either a :class:`Molecule` object or an :class:`rmgpy.species.Species`
        object.
        """

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

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

        # Check that the structure is not forbidden

        # If we're here then we're ready to make the new species
        if label == '': label = molecule.getFormula()
        logging.debug('Creating new species {0}'.format(label))
        spec = Species(index=self.speciesCounter+1, label=label, molecule=[molecule], reactive=reactive)
        spec.coreSizeAtCreation = len(self.core.species)
        spec.generateResonanceIsomers()
        spec.molecularWeight = Quantity(spec.molecule[0].getMolecularWeight()*1000.,"g/mol")
        spec.generateLennardJonesParameters()
        formula = molecule.getFormula()
        if formula in self.speciesDict:
            self.speciesDict[formula].append(spec)
        else:
            self.speciesDict[formula] = [spec]

        # Store hydrogens implicitly to conserve memory and speed up isomorphism
        for mol in spec.molecule:
            #mol.updateConnectivityValues()
            mol.makeHydrogensImplicit()

        self.speciesCounter += 1

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

        return spec, True

    def checkForExistingReaction(self, rxn):
        """
        Check to see if an existing reaction has the same reactants, products, and
        family as `rxn`. Returns :data:`True` or :data:`False` and the matched
        reaction (if found).
        """

        # Make sure the reactant and product lists are sorted before performing the check
        rxn.reactants.sort()
        rxn.products.sort()

        # Get the short-list of reactions with the same family, reactant1 and reactant2
        r1 = rxn.reactants[0]
        if len(rxn.reactants)==1: r2 = None
        else: r2 = rxn.reactants[1]
        family = rxn.getSource()
        try:
            my_reactionList = self.reactionDict[family][r1][r2]
        except KeyError: # no such short-list: must be new, unless in seed.
            my_reactionList = []

        # Now use short-list to check for matches. All should be in same forward direction.
        for rxn0 in my_reactionList:
            if (rxn0.reactants == rxn.reactants and rxn0.products == rxn.products):
                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
        for family0 in self.reactionDict:
            if isinstance(family, KineticsLibrary) and family0 != family:

                # First check seed short-list in forward direction
                try:
                    my_reactionList = self.reactionDict[family0][r1][r2]
                except KeyError:
                    my_reactionList = []
                for rxn0 in my_reactionList:
                    if (rxn0.reactants == rxn.reactants and rxn0.products == rxn.products) or \
                        (rxn0.reactants == rxn.products and rxn0.products == rxn.reactants):
                        return True, rxn0
                # Now get the seed short-list of the reverse reaction
                r1 = rxn.products[0]
                if len(rxn.products)==1: r2 = None
                else: r2 = rxn.products[1]
                try:
                    my_reactionList = self.reactionDict[family0][r1][r2]
                except KeyError:
                    my_reactionList = []
                for rxn0 in my_reactionList:
                    if (rxn0.reactants == rxn.reactants and rxn0.products == rxn.products) or \
                        (rxn0.reactants == rxn.products and rxn0.products == rxn.reactants):
                        return True, rxn0

        return False, None

    def makeNewReaction(self, forward, checkExisting=True):
        """
        Make a new reaction given a :class:`Reaction` object `forward`. The kinetics
        of the reaction are estimated and 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.
        """

        # Determine the proper species objects for all reactants and products
        forward.reactants = [self.makeNewSpecies(reactant)[0] for reactant in forward.reactants]
        forward.products  = [self.makeNewSpecies(product)[0]  for product  in forward.products ]

        if checkExisting:
            found, rxn = self.checkForExistingReaction(forward)
            if found: return rxn, False

        # Note in the log
        if isinstance(rxn, TemplateReaction):
            logging.debug('Creating new {0} reaction {1}'.format(forward.family.label, forward))
        else:
            logging.debug('Creating new library reaction {0}'.format(forward))
        
        # Add to the global dict/list of existing reactions (a list broken down by family, r1, r2)
        # identify r1 and r2
        r1 = forward.reactants[0]
        r2 = None if len(forward.reactants) == 1 else forward.reactants[1]
        family = forward.getSource()
        # make dictionary entries if necessary
        if family not in self.reactionDict:
            self.reactionDict[family] = {}
        if not self.reactionDict[family].has_key(r1):
            self.reactionDict[family][r1] = dict()
        if not self.reactionDict[family][r1].has_key(r2):
            self.reactionDict[family][r1][r2] = list()
        # store this reaction at the top of the relevant short-list
        self.reactionDict[family][r1][r2].insert(0, forward)

        forward.index = self.reactionCounter + 1
        self.reactionCounter += 1

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

        # Return newly created reaction
        return forward, True

    def makeNewPDepReaction(self, forward):
        """
        Make a new pressure-dependent reaction based on a list of `reactants` and a
        list of `products`. The reaction belongs to the specified `network` and
        has pressure-dependent kinetics given by `kinetics`.

        No checking for existing reactions is made here. The returned PDepReaction
        object is not added to the global list of reactions, as that is intended
        to represent only the high-pressure-limit set. The reactionCounter is
        incremented, however, since the returned reaction can and will exist in
        the model edge and/or core.
        """

        # Don't create reverse reaction: all such reactions are treated as irreversible
        # The reverse direction will come from a different partial network
        # Note that this isn't guaranteed to satisfy thermodynamics (but will probably be close)
        forward.reverse = None
        forward.reversible = False

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

        return forward

    def react(self, database, speciesA, speciesB=None):
        reactionList = []
        if speciesB is None:
            for moleculeA in speciesA.molecule:
                reactionList.extend(database.kinetics.generateReactions([moleculeA]))
                moleculeA.clearLabeledAtoms()
        else:
            for moleculeA in speciesA.molecule:
                for moleculeB in speciesB.molecule:
                    reactionList.extend(database.kinetics.generateReactions([moleculeA, moleculeB]))
                    moleculeA.clearLabeledAtoms()
                    moleculeB.clearLabeledAtoms()
        return reactionList

    def enlarge(self, newObject, database):
        """
        Enlarge a reaction model by processing `newObject`. If `newObject` is a
        :class:`rmg.species.Species` object, then the species is moved from
        the edge to the core and reactions generated for that species, reacting
        with itself and with all other species in the model core. If `newObject`
        is a :class:`rmg.unirxn.network.Network` object, then reactions are
        generated for the species in the network with the largest leak flux.
        """

        self.newReactionList = []; self.newSpeciesList = []
        newReactions = []

        numOldCoreSpecies = len(self.core.species)
        numOldCoreReactions = len(self.core.reactions)
        numOldEdgeSpecies = len(self.edge.species)
        numOldEdgeReactions = len(self.edge.reactions)

        pdepNetwork = None

        if isinstance(newObject, Species):

            newSpecies = newObject

            if not newSpecies.reactive:
                logging.info('NOT generating reactions for unreactive species {0}'.format(newSpecies))
            else:
                logging.info('Adding species {0} to model core'.format(newSpecies))
                # Find reactions involving the new species as unimolecular reactant
                # or product (e.g. A <---> products)
                newReactions.extend(self.react(database, newSpecies))
                # Find reactions involving the new species as bimolecular reactants
                # or products with other core species (e.g. A + B <---> products)
                for coreSpecies in self.core.species:
                    if coreSpecies.reactive:
                        newReactions.extend(self.react(database, newSpecies, coreSpecies))
                # Find reactions involving the new species as bimolecular reactants
                # or products with itself (e.g. A + A <---> products)
                newReactions.extend(self.react(database, newSpecies, newSpecies))

            # Add new species
            self.addSpeciesToCore(newSpecies)

            # Process the new reactions
            self.processNewReactions(newReactions, newSpecies, pdepNetwork)

        elif isinstance(newObject, tuple) and isinstance(newObject[0], PDepNetwork) and settings.pressureDependence:

            pdepNetwork, newSpecies = newObject
            newReactions.extend(pdepNetwork.exploreIsomer(newSpecies, self, database))
            self.processNewReactions(newReactions, newSpecies, pdepNetwork)

{'fraction_non_alphanumeric': 0.10116572109595164, 'fraction_numerical': 0.009012923489466302, 'mean_word_length': 2.245577984608566, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4509400', 'n_tokens_mistral': 22290, 'n_tokens_neox': 20586, 'n_words': 6913}