From 42dae5ee5cebf738fc7c9368a207b83e8e477dff Mon Sep 17 00:00:00 2001 From: zjburas Date: Mon, 18 Sep 2017 22:42:02 -0400 Subject: [PATCH 1/3] Created a new reactor type MBSampledReactor For modeling molecular beam sampling of species from an isothermal, isobaric, homogenous batch reactor. The sampling process is modeled as two first-order unimolecular reactions, following the approach developed by Baeza-Romero et al., in their 2011 IJCK paper. Currently, the new reactor is only intended for use with the simulate.py script, after a conventional RMG mechanism has already been made. --- rmgpy/chemkin.pyx | 12 +- rmgpy/rmg/input.py | 67 ++++- rmgpy/solver/mbSampled.pyx | 494 +++++++++++++++++++++++++++++++++++++ rmgpy/tools/loader.py | 23 +- setup.py | 1 + 5 files changed, 593 insertions(+), 4 deletions(-) create mode 100644 rmgpy/solver/mbSampled.pyx diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index 5848fa5a8a..0fe1a80b96 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -1411,7 +1411,11 @@ def getSpeciesIdentifier(species): # Try the chemical formula name = '{0}({1:d})'.format(species.molecule[0].getFormula(), species.index) if len(name) <= 10: - return name + if 'obs' in label: + # For MBSampledReactor, keep observed species tag + return name + '_obs' + else: + return name # As a last resort, just use the index if species.index >= 0: @@ -1421,7 +1425,11 @@ def getSpeciesIdentifier(species): else: name = 'S({0:d})'.format(species.index) if len(name) <= 10: - return name + if 'obs' in label: + # For MBSampledReactor, keep observed species tag + return name + '_obs' + else: + return name # If we're here then we just can't come up with a valid Chemkin name # for this species, so raise an exception diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index 009ea77789..fa4cb0ae16 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -38,10 +38,11 @@ from rmgpy import settings from rmgpy.molecule import Molecule -from rmgpy.quantity import Quantity, Energy, SurfaceConcentration +from rmgpy.quantity import Quantity, Energy, RateCoefficient, SurfaceConcentration from rmgpy.solver.base import TerminationTime, TerminationConversion, TerminationRateRatio from rmgpy.solver.simple import SimpleReactor from rmgpy.solver.liquid import LiquidReactor +from rmgpy.solver.mbSampled import MBSampledReactor from rmgpy.solver.surface import SurfaceReactor from rmgpy.rmg.settings import ModelSettings, SimulatorSettings from model import CoreEdgeReactionModel @@ -447,6 +448,69 @@ def surfaceReactor(temperature, rmg.reactionSystems.append(system) system.log_initial_conditions(number=len(rmg.reactionSystems)) + +# Reaction systems +def mbsampledReactor(temperature, + pressure, + initialMoleFractions, + mbsamplingRate, + terminationConversion=None, + terminationTime=None, + sensitivity=None, + sensitivityThreshold=1e-3, + constantSpecies=None, + ): + logging.debug('Found MBSampledReactor reaction system') + + for value in initialMoleFractions.values(): + if value < 0: + raise InputError('Initial mole fractions cannot be negative.') + + for spec in initialMoleFractions: + initialMoleFractions[spec] = float(initialMoleFractions[spec]) + + totalInitialMoles = sum(initialMoleFractions.values()) + if totalInitialMoles != 1: + logging.warning('Initial mole fractions do not sum to one; normalizing.') + logging.info('') + logging.info('Original composition:') + for spec, molfrac in initialMoleFractions.iteritems(): + logging.info("{0} = {1}".format(spec, molfrac)) + for spec in initialMoleFractions: + initialMoleFractions[spec] /= totalInitialMoles + logging.info('') + logging.info('Normalized mole fractions:') + for spec, molfrac in initialMoleFractions.iteritems(): + logging.info("{0} = {1}".format(spec, molfrac)) + + T = Quantity(temperature) + P = Quantity(pressure) + + k_sampling = RateCoefficient(mbsamplingRate, 's^-1') + + constantSpeciesList = [] + + for spec in constantSpecies: + constantSpeciesList.append(speciesDict[spec]) + + termination = [] + if terminationConversion is not None: + for spec, conv in terminationConversion.iteritems(): + termination.append(TerminationConversion(speciesDict[spec], conv)) + if terminationTime is not None: + termination.append(TerminationTime(Quantity(terminationTime))) + if len(termination) == 0: + raise InputError( + 'No termination conditions specified for reaction system #{0}.'.format(len(rmg.reactionSystems) + 2)) + + sensitiveSpecies = [] + if sensitivity: + if isinstance(sensitivity, str): sensitivity = [sensitivity] + for spec in sensitivity: + sensitiveSpecies.append(speciesDict[spec]) + system = MBSampledReactor(T, P, initialMoleFractions, k_sampling, constantSpeciesList, termination, sensitiveSpecies, sensitivityThreshold) + rmg.reactionSystems.append(system) + def simulator(atol, rtol, sens_atol=1e-6, sens_rtol=1e-4): rmg.simulatorSettingsList.append(SimulatorSettings(atol, rtol, sens_atol, sens_rtol)) @@ -756,6 +820,7 @@ def readInputFile(path, rmg0): 'simpleReactor': simpleReactor, 'liquidReactor': liquidReactor, 'surfaceReactor': surfaceReactor, + 'mbsampledReactor': mbsampledReactor, 'simulator': simulator, 'solvation': solvation, 'model': model, diff --git a/rmgpy/solver/mbSampled.pyx b/rmgpy/solver/mbSampled.pyx new file mode 100644 index 0000000000..c608fbdd69 --- /dev/null +++ b/rmgpy/solver/mbSampled.pyx @@ -0,0 +1,494 @@ +################################################################################ +# +# RMG - Reaction Mechanism Generator +# +# Copyright (c) 2002-2017 Prof. William H. Green (whgreen@mit.edu), +# Prof. Richard H. West (r.west@neu.edu) and the RMG Team (rmg_dev@mit.edu) +# +# Permission is hereby granted, free of charge, to any person obtaining a +# copy of this software and associated documentation files (the 'Software'), +# to deal in the Software without restriction, including without limitation +# the rights to use, copy, modify, merge, publish, distribute, sublicense, +# and/or sell copies of the Software, and to permit persons to whom the +# Software is furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +# DEALINGS IN THE SOFTWARE. +# +################################################################################ + +""" +Contains the :class:`MBSampledReactor` class, providing a reaction system +consisting of a homogeneous, isothermal, isobaric batch reactor that is being +sampled by a molecular beam (MB). +""" + +import numpy, logging +cimport numpy + +import itertools + +from base cimport ReactionSystem +cimport cython + +import rmgpy.constants as constants +cimport rmgpy.constants as constants +from rmgpy.quantity import Quantity, RateCoefficient +from rmgpy.quantity cimport ScalarQuantity, ArrayQuantity + +cdef class MBSampledReactor(ReactionSystem): + """ + A reaction system consisting of a homogeneous, isothermal, isobaric batch + reactor that is being sample by a molecular beam. + The sampling process is modeled as a unimolecular reaction. + These assumptions allow for a number of optimizations that enable + this solver to complete very rapidly, even for large kinetic models. + + This is currently only intended for use with the ``simulate.py``script, + and cannot be used for a standard RMG job. + """ + + cdef public ScalarQuantity T + cdef public ScalarQuantity P + cdef public ScalarQuantity k_sampling + cdef public double V + cdef public bint constantVolume + cdef public dict initialMoleFractions + cdef public list constantSpeciesList + + # collider variables + + """ + pdepColliderKinetics: + an array that contains a reference to the kinetics object of the reaction + that has pressure dependent kinetics. + """ + cdef public list pdepColliderKinetics + + """ + colliderEfficiencies: + an array consisting of array elements, each element corresponding to a reaction. + Each element is an array with each position in the array corresponding to the collider efficiency + of the core species. The collider efficiency is set to 1 if the species was not found in the list + of colliders. + """ + cdef public numpy.ndarray colliderEfficiencies + + """ + pdepColliderReactionIndices: + array that contains the indices of those reactions that + have pressure dependent kinetics. E.g. [4, 10, 2, 123] + """ + cdef public numpy.ndarray pdepColliderReactionIndices + + """ + pdepSpecificColliderKinetics: + an array that contains a reference to the kinetics object of the reaction + that has pressure dependent kinetics with a specific species as a third body collider. + """ + cdef public list pdepSpecificColliderKinetics + + """ + specificColliderSpecies: + a list that contains object references to species which are specific third body colliders + in the respective reactions in pdepSpecificColliderReactionIndices. + """ + cdef public list specificColliderSpecies + + """ + pdepSpecificColliderReactionIndices: + an array that contains the indices of reactions that have + a specifcCollider attribyte. E.g. [16, 155, 90] + """ + cdef public numpy.ndarray pdepSpecificColliderReactionIndices + + def __init__(self, T, P, initialMoleFractions, k_sampling, constantSpeciesList, termination, sensitiveSpecies=None, sensitivityThreshold=1e-3): + ReactionSystem.__init__(self, termination, sensitiveSpecies, sensitivityThreshold) + self.T = Quantity(T) + self.P = Quantity(P) + self.initialMoleFractions = initialMoleFractions + self.k_sampling = RateCoefficient(k_sampling) + self.constantSpeciesList = constantSpeciesList + + self.V = 0 # will be set in initializeModel + self.constantVolume = False + + self.pdepColliderReactionIndices = None + self.pdepColliderKinetics = None + self.colliderEfficiencies = None + self.pdepSpecificColliderReactionIndices = None + self.pdepSpecificColliderKinetics = None + self.specificColliderSpecies = None + + def __reduce__(self): + """ + A helper function used when pickling an object. + """ + return (self.__class__, + (self.T, self.P, self.initialMoleFractions, self.termination)) + + def convertInitialKeysToSpeciesObjects(self, speciesDict): + """ + Convert the initialMoleFractions dictionary from species names into species objects, + using the given dictionary of species. + """ + initialMoleFractions = {} + for label, moleFrac in self.initialMoleFractions.iteritems(): + initialMoleFractions[speciesDict[label]] = moleFrac + self.initialMoleFractions = initialMoleFractions + + cpdef initializeModel(self, list coreSpecies, list coreReactions, list edgeSpecies, list edgeReactions, + list surfaceSpecies=None, list surfaceReactions=None, list pdepNetworks=None, + atol=1e-16, rtol=1e-8, sensitivity=False, sens_atol=1e-6, sens_rtol=1e-4, + filterReactions=False, dict conditions=None): + """ + Initialize a simulation of the reaction system using the provided + kinetic model. You will probably want to create your own version of this + method in the derived class; don't forget to also call the base class + version, too. + """ + if surfaceSpecies is None: + surfaceSpecies = [] + if surfaceReactions is None: + surfaceReactions = [] + + # First call the base class version of the method + # This initializes the attributes declared in the base class + ReactionSystem.initializeModel(self, coreSpecies=coreSpecies, coreReactions=coreReactions, edgeSpecies=edgeSpecies, + edgeReactions=edgeReactions, surfaceSpecies=surfaceSpecies, surfaceReactions=surfaceReactions, + pdepNetworks=pdepNetworks, atol=atol, rtol=rtol, sensitivity=sensitivity, sens_atol=sens_atol, + sens_rtol=sens_rtol) + + # Set initial conditions + self.set_initial_conditions() + + # Compute reaction thresholds if reaction filtering is turned on + if filterReactions: + ReactionSystem.set_initial_reaction_thresholds(self) + + self.set_colliders(coreReactions, edgeReactions, coreSpecies) + + ReactionSystem.compute_network_variables(self, pdepNetworks) + + # Generate forward and reverse rate coefficients k(T,P) + self.generate_rate_coefficients(coreReactions, edgeReactions) + + ReactionSystem.set_initial_derivative(self) + # Initialize the model + ReactionSystem.initialize_solver(self) + + def calculate_effective_pressure(self, rxn): + """ + Computes the effective pressure for a reaction as: + Peff = P * sum(yi * effi / sum(y)) + with: + - P the pressure of the reactor, + - y the array of initial moles of the core species + or as: + Peff = P * y_specificCollider / sum(y) + if a secificCollider is mentioned + """ + + y0_coreSpecies = self.y0[:self.numCoreSpecies] + sum_core_species = numpy.sum(y0_coreSpecies) + + j = self.reactionIndex[rxn] + for i in xrange(self.pdepColliderReactionIndices.shape[0]): + if j == self.pdepColliderReactionIndices[i]: + # Calculate effective pressure + if rxn.specificCollider is None: + Peff = self.P.value_si * numpy.sum(self.colliderEfficiencies[i]*y0_coreSpecies / sum_core_species) + else: + logging.debug("Calculating Peff using {0} as a specificCollider".format(rxn.specificCollider)) + Peff = self.P.value_si * self.y0[self.speciesIndex[rxn.specificCollider]] / sum_core_species + return Peff + return self.P.value_si + + def generate_rate_coefficients(self, coreReactions, edgeReactions): + """ + Populates the forward rate coefficients (kf), reverse rate coefficients (kb) + and equilibrium constants (Keq) arrays with the values computed at the temperature + and (effective) pressure of the reaction system. + """ + + for rxn in itertools.chain(coreReactions, edgeReactions): + j = self.reactionIndex[rxn] + Peff = self.calculate_effective_pressure(rxn) + self.kf[j] = rxn.getRateCoefficient(self.T.value_si, Peff) + + if rxn.reversible: + self.Keq[j] = rxn.getEquilibriumConstant(self.T.value_si) + self.kb[j] = self.kf[j] / self.Keq[j] + + def set_colliders(self, coreReactions, edgeReactions, coreSpecies): + """ + Store collider efficiencies and reaction indices for pdep reactions that have collider efficiencies, + and store specific collider indices + """ + pdepColliderReactionIndices = [] + self.pdepColliderKinetics = [] + colliderEfficiencies = [] + pdepSpecificColliderReactionIndices = [] + self.pdepSpecificColliderKinetics = [] + self.specificColliderSpecies = [] + + for rxn in itertools.chain(coreReactions, edgeReactions): + if rxn.kinetics.isPressureDependent(): + if rxn.kinetics.efficiencies: + j = self.reactionIndex[rxn] + pdepColliderReactionIndices.append(j) + self.pdepColliderKinetics.append(rxn.kinetics) + colliderEfficiencies.append(rxn.kinetics.getEffectiveColliderEfficiencies(coreSpecies)) + if rxn.specificCollider: + pdepSpecificColliderReactionIndices.append(self.reactionIndex[rxn]) + self.pdepSpecificColliderKinetics.append(rxn.kinetics) + self.specificColliderSpecies.append(rxn.specificCollider) + + self.pdepColliderReactionIndices = numpy.array(pdepColliderReactionIndices, numpy.int) + self.colliderEfficiencies = numpy.array(colliderEfficiencies, numpy.float64) + self.pdepSpecificColliderReactionIndices = numpy.array(pdepSpecificColliderReactionIndices, numpy.int) + + def set_initial_conditions(self): + """ + Sets the initial conditions of the rate equations that represent the + current reactor model. + + The volume is set to the value derived from the ideal gas law, using the + user-defined pressure, temperature, and the number of moles of initial species. + + The species moles array (y0) is set to the values stored in the + initial mole fractions dictionary. + + The initial species concentration is computed and stored in the + coreSpeciesConcentrations array. + + """ + + ReactionSystem.set_initial_conditions(self) + + for spec, moleFrac in self.initialMoleFractions.iteritems(): + i = self.get_species_index(spec) + self.y0[i] = moleFrac + + # Use ideal gas law to compute volume + self.V = constants.R * self.T.value_si * numpy.sum(self.y0[:self.numCoreSpecies]) / self.P.value_si# volume in m^3 + for j in xrange(self.numCoreSpecies): + self.coreSpeciesConcentrations[j] = self.y0[j] / self.V + + @cython.boundscheck(False) + def residual(self, double t, numpy.ndarray[numpy.float64_t, ndim=1] y, numpy.ndarray[numpy.float64_t, ndim=1] dydt, + numpy.ndarray[numpy.float64_t, ndim=1] senpar = numpy.zeros(1, numpy.float64)): + """ + Return the residual function for the governing DAE system for the + simple reaction system. + """ + cdef numpy.ndarray[numpy.int_t, ndim=2] ir, ip, inet + cdef numpy.ndarray[numpy.float64_t, ndim=1] res, kf, kr, knet, delta, equilibriumConstants + cdef int numCoreSpecies, numCoreReactions, numEdgeSpecies, numEdgeReactions, numPdepNetworks + cdef int i, j, z, first, second, third, realspeciesIndex + cdef double k, V, reactionRate, revReactionRate, T, P, Peff, coreSpeciesRate + cdef numpy.ndarray[numpy.float64_t, ndim=1] coreSpeciesConcentrations, coreSpeciesRates, coreReactionRates, edgeSpeciesRates, edgeReactionRates, networkLeakRates, coreSpeciesConsumptionRates, coreSpeciesProductionRates + cdef numpy.ndarray[numpy.float64_t, ndim=1] C, y_coreSpecies + cdef numpy.ndarray[numpy.float64_t, ndim=2] jacobian, dgdk, colliderEfficiencies + cdef numpy.ndarray[numpy.int_t, ndim=1] pdepColliderReactionIndices, pdepSpecificColliderReactionIndices + cdef list pdepColliderKinetics, pdepSpecificColliderKinetics + + ir = self.reactantIndices + ip = self.productIndices + + numCoreSpecies = len(self.coreSpeciesRates) + numCoreReactions = len(self.coreReactionRates) + numEdgeSpecies = len(self.edgeSpeciesRates) + numEdgeReactions = len(self.edgeReactionRates) + numPdepNetworks = len(self.networkLeakRates) + kf = self.kf + kr = self.kb + + y_coreSpecies = y[:numCoreSpecies] + + # Recalculate any forward and reverse rate coefficients that involve pdep collision efficiencies + if self.pdepColliderReactionIndices.shape[0] != 0: + T = self.T.value_si + P = self.P.value_si + equilibriumConstants = self.Keq + pdepColliderReactionIndices = self.pdepColliderReactionIndices + pdepColliderKinetics = self.pdepColliderKinetics + colliderEfficiencies = self.colliderEfficiencies + for i in xrange(pdepColliderReactionIndices.shape[0]): + # Calculate effective pressure + Peff = P*numpy.sum(colliderEfficiencies[i]*y_coreSpecies / numpy.sum(y_coreSpecies)) + j = pdepColliderReactionIndices[i] + kf[j] = pdepColliderKinetics[i].getRateCoefficient(T, Peff) + kr[j] = kf[j] / equilibriumConstants[j] + if self.pdepSpecificColliderReactionIndices.shape[0] != 0: + T = self.T.value_si + P = self.P.value_si + equilibriumConstants = self.Keq + pdepSpecificColliderReactionIndices = self.pdepSpecificColliderReactionIndices + pdepSpecificColliderKinetics = self.pdepSpecificColliderKinetics + specificColliderSpecies = self.specificColliderSpecies + for i in xrange(pdepSpecificColliderReactionIndices.shape[0]): + # Calculate effective pressure + Peff = P * y[self.speciesIndex[specificColliderSpecies[i]]] / numpy.sum(y_coreSpecies) + j = pdepSpecificColliderReactionIndices[i] + kf[j] = pdepSpecificColliderKinetics[i].getRateCoefficient(T, Peff) + kr[j] = kf[j] / equilibriumConstants[j] + + inet = self.networkIndices + knet = self.networkLeakCoefficients + + + res = numpy.zeros(numCoreSpecies, numpy.float64) + + coreSpeciesConcentrations = numpy.zeros_like(self.coreSpeciesConcentrations) + coreSpeciesRates = numpy.zeros_like(self.coreSpeciesRates) + coreReactionRates = numpy.zeros_like(self.coreReactionRates) + coreSpeciesConsumptionRates = numpy.zeros_like(self.coreSpeciesConsumptionRates) + coreSpeciesProductionRates = numpy.zeros_like(self.coreSpeciesProductionRates) + edgeSpeciesRates = numpy.zeros_like(self.edgeSpeciesRates) + edgeReactionRates = numpy.zeros_like(self.edgeReactionRates) + networkLeakRates = numpy.zeros_like(self.networkLeakRates) + + C = numpy.zeros_like(self.coreSpeciesConcentrations) + + # Use ideal gas law to compute volume + V = constants.R * self.T.value_si * numpy.sum(y_coreSpecies) / self.P.value_si + self.V = V + + for j in xrange(numCoreSpecies): + C[j] = y[j] / V + coreSpeciesConcentrations[j] = C[j] + + for j in xrange(ir.shape[0]): + k = kf[j] + if ir[j,0] >= numCoreSpecies or ir[j,1] >= numCoreSpecies or ir[j,2] >= numCoreSpecies: + fReactionRate = 0.0 + elif ir[j,1] == -1: # only one reactant + fReactionRate = k * C[ir[j,0]] + elif ir[j,2] == -1: # only two reactants + fReactionRate = k * C[ir[j,0]] * C[ir[j,1]] + else: # three reactants!! (really?) + fReactionRate = k * C[ir[j,0]] * C[ir[j,1]] * C[ir[j,2]] + k = kr[j] + if ip[j,0] >= numCoreSpecies or ip[j,1] >= numCoreSpecies or ip[j,2] >= numCoreSpecies: + revReactionRate = 0.0 + elif ip[j,1] == -1: # only one reactant + revReactionRate = k * C[ip[j,0]] + elif ip[j,2] == -1: # only two reactants + revReactionRate = k * C[ip[j,0]] * C[ip[j,1]] + else: # three reactants!! (really?) + revReactionRate = k * C[ip[j,0]] * C[ip[j,1]] * C[ip[j,2]] + + reactionRate = fReactionRate-revReactionRate + + # Set the reaction and species rates + if j < numCoreReactions: + # The reaction is a core reaction + coreReactionRates[j] = reactionRate + + # Add/substract the total reaction rate from each species rate + # Since it's a core reaction we know that all of its reactants + # and products are core species + first = ir[j,0] + coreSpeciesRates[first] -= reactionRate + coreSpeciesConsumptionRates[first] += fReactionRate + coreSpeciesProductionRates[first] += revReactionRate + second = ir[j,1] + if second != -1: + coreSpeciesRates[second] -= reactionRate + coreSpeciesConsumptionRates[second] += fReactionRate + coreSpeciesProductionRates[second] += revReactionRate + third = ir[j,2] + if third != -1: + coreSpeciesRates[third] -= reactionRate + coreSpeciesConsumptionRates[third] += fReactionRate + coreSpeciesProductionRates[third] += revReactionRate + first = ip[j,0] + coreSpeciesRates[first] += reactionRate + coreSpeciesProductionRates[first] += fReactionRate + coreSpeciesConsumptionRates[first] += revReactionRate + second = ip[j,1] + if second != -1: + coreSpeciesRates[second] += reactionRate + coreSpeciesProductionRates[second] += fReactionRate + coreSpeciesConsumptionRates[second] += revReactionRate + third = ip[j,2] + if third != -1: + coreSpeciesRates[third] += reactionRate + coreSpeciesProductionRates[third] += fReactionRate + coreSpeciesConsumptionRates[third] += revReactionRate + + else: + # The reaction is an edge reaction + edgeReactionRates[j-numCoreReactions] = reactionRate + + # Add/substract the total reaction rate from each species rate + # Since it's an edge reaction its reactants and products could + # be either core or edge species + # We're only interested in the edge species + first = ir[j,0] + if first >= numCoreSpecies: edgeSpeciesRates[first-numCoreSpecies] -= reactionRate + second = ir[j,1] + if second != -1: + if second >= numCoreSpecies: edgeSpeciesRates[second-numCoreSpecies] -= reactionRate + third = ir[j,2] + if third != -1: + if third >= numCoreSpecies: edgeSpeciesRates[third-numCoreSpecies] -= reactionRate + first = ip[j,0] + if first >= numCoreSpecies: edgeSpeciesRates[first-numCoreSpecies] += reactionRate + second = ip[j,1] + if second != -1: + if second >= numCoreSpecies: edgeSpeciesRates[second-numCoreSpecies] += reactionRate + third = ip[j,2] + if third != -1: + if third >= numCoreSpecies: edgeSpeciesRates[third-numCoreSpecies] += reactionRate + + for j in xrange(inet.shape[0]): + k = knet[j] + if inet[j,1] == -1: # only one reactant + reactionRate = k * C[inet[j,0]] + elif inet[j,2] == -1: # only two reactants + reactionRate = k * C[inet[j,0]] * C[inet[j,1]] + else: # three reactants!! (really?) + reactionRate = k * C[inet[j,0]] * C[inet[j,1]] * C[inet[j,2]] + networkLeakRates[j] = reactionRate + + #Add unimolecular reactions describing the transport of sampled molecules to the measurement zone (i.e., + #ionization region of a Mass Spectrometer, intersection with probing laser, etc.) + for j, coreSpeciesRate in enumerate(coreSpeciesRates): + species = self.speciesIndex.keys()[self.speciesIndex.values().index(j)] + if '_obs' in species.label: + for species2 in self.speciesIndex.keys(): + if species2.label == species.label.replace('_obs',''): + realspeciesIndex = self.speciesIndex[species2] + if realspeciesIndex > len(C) - 1: + continue + else: + coreSpeciesRates[j] = self.k_sampling.value_si * C[realspeciesIndex] - self.k_sampling.value_si * C[j] + coreSpeciesConsumptionRates[j] = self.k_sampling.value_si * C[j] + coreSpeciesProductionRates[j] = self.k_sampling.value_si * C[realspeciesIndex] + break + + self.coreSpeciesConcentrations = coreSpeciesConcentrations + self.coreSpeciesRates = coreSpeciesRates + self.coreSpeciesProductionRates = coreSpeciesProductionRates + self.coreSpeciesConsumptionRates = coreSpeciesConsumptionRates + self.coreReactionRates = coreReactionRates + self.edgeSpeciesRates = edgeSpeciesRates + self.edgeReactionRates = edgeReactionRates + self.networkLeakRates = networkLeakRates + + res = coreSpeciesRates * V + + delta = res - dydt + + # Return DELTA, IRES. IRES is set to 1 in order to tell DASPK to evaluate the sensitivity residuals + return delta, 1 diff --git a/rmgpy/tools/loader.py b/rmgpy/tools/loader.py index b600381a7e..d666d19b42 100644 --- a/rmgpy/tools/loader.py +++ b/rmgpy/tools/loader.py @@ -37,6 +37,7 @@ import warnings from rmgpy.chemkin import loadChemkinFile from rmgpy.solver.liquid import LiquidReactor +from rmgpy.solver.mbSampled import MBSampledReactor from rmgpy.solver.surface import SurfaceReactor from rmgpy.solver.base import TerminationConversion @@ -77,6 +78,27 @@ def loadRMGPyJob(inputFile, chemkinFile=None, speciesDict=None, generateImages=T speciesList, reactionList = loadChemkinFile(chemkinFile, speciesDict, useChemkinNames=useChemkinNames, checkDuplicates=checkDuplicates) + # Created "observed" versions of all reactive species that are not explicitly + # identified as "constant" species + for reactionSystem in rmg.reactionSystems: + if isinstance(reactionSystem, MBSampledReactor): + observedspeciesList = [] + for species in speciesList: + if '_obs' not in species.label and species.reactive: + for constantSpecies in reactionSystem.constantSpeciesList: + if species.isIsomorphic(constantSpecies): + break + else: + for species2 in speciesList: + if species2.label == species.label + '_obs': + break + else: + observedspecies = species.copy(deep=True) + observedspecies.label = species.label + '_obs' + observedspeciesList.append(observedspecies) + + speciesList.extend(observedspeciesList) + # Map species in input file to corresponding species in Chemkin file speciesDict = {} for spec0 in rmg.initialSpecies: @@ -210,4 +232,3 @@ def loadRMGJavaJob(inputFile, chemkinFile=None, speciesDict=None, generateImages species.molecule[0].draw(str(path)) return rmg - diff --git a/setup.py b/setup.py index cdd4405435..10b7ef2af5 100644 --- a/setup.py +++ b/setup.py @@ -122,6 +122,7 @@ def getSolverExtensionModules(): Extension('rmgpy.solver.base', ['rmgpy/solver/base.pyx'], include_dirs=['.']), Extension('rmgpy.solver.simple', ['rmgpy/solver/simple.pyx'], include_dirs=['.']), Extension('rmgpy.solver.liquid', ['rmgpy/solver/liquid.pyx'], include_dirs=['.']), + Extension('rmgpy.solver.mbSampled', ['rmgpy/solver/mbSampled.pyx'], include_dirs=['.']), Extension('rmgpy.solver.surface', ['rmgpy/solver/surface.pyx'], include_dirs=['.']), ] From d5fc97d2f279de5994c8c83f9b28f1df7006533d Mon Sep 17 00:00:00 2001 From: Max Liu Date: Sun, 4 Aug 2019 22:21:10 -0400 Subject: [PATCH 2/3] Add unit test for simulation of MBSampledReactor --- rmgpy/tools/data/sim/mbSampled/chem.inp | 440 ++++++++++++++++++ rmgpy/tools/data/sim/mbSampled/input.py | 87 ++++ .../data/sim/mbSampled/species_dictionary.txt | 249 ++++++++++ rmgpy/tools/simulateTest.py | 23 +- 4 files changed, 797 insertions(+), 2 deletions(-) create mode 100644 rmgpy/tools/data/sim/mbSampled/chem.inp create mode 100644 rmgpy/tools/data/sim/mbSampled/input.py create mode 100644 rmgpy/tools/data/sim/mbSampled/species_dictionary.txt diff --git a/rmgpy/tools/data/sim/mbSampled/chem.inp b/rmgpy/tools/data/sim/mbSampled/chem.inp new file mode 100644 index 0000000000..13a454f23d --- /dev/null +++ b/rmgpy/tools/data/sim/mbSampled/chem.inp @@ -0,0 +1,440 @@ +ELEMENTS + H + D /2.014/ + T /3.016/ + C + CI /13.003/ + O + OI /18.000/ + N + Ne + Ar + He + Si + S + Cl +END + +SPECIES + C6H5 + C2H4 + C6H6 + C2H3 + C6H5C2H3 + H + i1 + i2 + i3 + i4 + i5 + i6 + i7 + i8 + i9 + He +END + + + +THERM ALL + 300.000 1000.000 5000.000 + +! Narayanaswamy +C6H5 H 5C 6 G 100.000 5000.000 845.04 1 + 7.43177438E+00 2.52167887E-02-1.11182014E-05 2.22179629E-09-1.64849859E-13 2 + 3.69251718E+04-1.62469857E+01-2.38126808E-01 2.56754111E-02 5.16978224E-05 3 +-9.75334271E-08 4.41979860E-11 3.95013542E+04 2.70374991E+01 4 + +! Narayanaswamy +C2H4 H 4C 2 G 100.000 5000.000 889.17 1 + 2.57579789E+00 1.32188903E-02-5.48150328E-06 1.04904217E-09-7.54583209E-14 2 + 4.80469837E+03 7.51939203E+00 3.59281423E+00-3.75064618E-03 4.40545192E-05 3 +-5.17680301E-08 1.91823593E-11 5.11380306E+03 5.48752425E+00 4 + +! Narayanaswamy +C6H6 H 6C 6 G 100.000 5000.000 846.39 1 + 6.88981564E+00 2.87389241E-02-1.25905881E-05 2.50537788E-09-1.85350036E-13 2 + 6.00007188E+03-1.66565842E+01-3.83583688E-01 2.39231158E-02 6.53972086E-05 3 +-1.13627777E-07 5.02755970E-11 8.63502031E+03 2.55127510E+01 4 + +! Narayanaswamy +C2H3 H 3C 2 G 100.000 5000.000 901.39 1 + 3.44360414E+00 9.27038173E-03-3.80263060E-06 7.22486378E-10-5.17120430E-14 2 + 3.44940404E+04 5.54954191E+00 2.98095681E+00 4.01214338E-03 1.71141801E-05 3 +-2.37458916E-08 9.23024608E-12 3.48744695E+04 9.38110056E+00 4 + +! Narayanaswamy +C6H5C2H3 H 8C 8 G 100.000 5000.000 1047.78 1 + 1.40832380E+01 3.09758115E-02-1.27931070E-05 2.42245720E-09-1.73876461E-13 2 + 1.05677156E+04-5.11291054E+01 1.41597467E+00 4.24124721E-02 2.36912100E-05 3 +-5.44222614E-08 2.14135715E-11 1.52489424E+04 2.02458503E+01 4 + +! primaryThermoLibrary +H H 1 G 100.000 5000.000 3142.77 1 + 2.50000000E+00 3.32026413E-12-1.37161245E-15 2.48094222E-19-1.65799619E-23 2 + 2.54742178E+04-4.44972878E-01 2.50000000E+00 9.62838494E-15-1.50167895E-17 3 + 7.73923008E-21-1.23210440E-24 2.54742178E+04-4.44972896E-01 4 + +! Narayanaswamy +i1 H 9C 8 G 100.000 5000.000 1027.73 1 + 1.29602849E+01 3.44256941E-02-1.41379607E-05 2.66396201E-09-1.88898251E-13 2 + 2.20240268E+04-4.26574166E+01 1.06920806E+00 5.13336312E-02 4.05528157E-06 3 +-3.69471051E-08 1.62114832E-11 2.60194156E+04 2.25822720E+01 4 + +! CCSD(T)-F12/cc-pVTZ-F12 calculation, with 1-D hindered rotor +i2 H 9C 8 G 10.000 3000.000 494.35 1 +-2.20727418E+00 7.18030212E-02-4.58339987E-05 1.40113251E-08-1.64310381E-12 2 + 2.13352323E+04 3.68838437E+01 3.90724086E+00 5.76073074E-03 2.04827457E-04 3 +-3.91817793E-07 2.37875592E-10 2.09331219E+04 1.37400593E+01 4 + +! Thermo group additivity estimation: group(Cb-H) + group(Cb-H) + group(Cb-H) + group(Cb-H) + group(Cb-H) + group(Cb-Cs) + group(Cs-CsHHH) + group(Cs- +! CbCsHH) + ring(Benzene) + radical(CbJ) +i3 H 9C 8 G 100.000 5000.000 982.26 1 + 1.46961187E+01 2.95647284E-02-1.09249131E-05 2.04710319E-09-1.48900204E-13 2 + 2.70636886E+04-5.07831711E+01 1.28416603E+00 4.25780751E-02 3.27346770E-05 3 +-7.07047023E-08 2.93421648E-11 3.17055107E+04 2.38980107E+01 4 + +! Thermo group additivity estimation: group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cs-(Cds-Cds)(Cds-Cds)HH) + +! group(Cs-CsCsHH) + group(Cs-CsCsHH) + group(Cs-(Cds-Cds)(Cds-Cds)CsCs) + polycyclic(s1_3_6_diene_1_4) + radical(C=CCJC=C-cyclohexadiene) +i4 H 9C 8 G 100.000 5000.000 1421.74 1 + 8.28198008E+00 1.84893779E-02 5.58691374E-06-2.15788590E-09 1.85746382E-13 2 + 3.40548623E+04-1.17200135E+01-4.28932547E-01 8.12428030E-02-1.00971926E-04 3 + 6.67296285E-08-1.52546181E-11 3.26664299E+04 1.97764999E+01 4 + +! Thermo group additivity estimation: group(Cds-Cds(Cds-Cds)H) + group(Cds-Cds(Cds-Cds)H) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cs-CsCsHH) + +! group(Cs-CsCsHH) + group(Cs-(Cds-Cds)CsCsH) + group(Cs-(Cds-Cds)CsCsH) + polycyclic(s2_4_6_diene_1_3) + radical(Allyl_T) +i5 H 9C 8 G 100.000 5000.000 982.89 1 + 1.23035693E+01 3.46969352E-02-1.32391077E-05 2.54167765E-09-1.87713284E-13 2 + 3.50570935E+04-4.58231453E+01 2.02886321E+00 2.15595261E-02 9.06732887E-05 3 +-1.24821832E-07 4.66758048E-11 3.97312139E+04 1.70718461E+01 4 + +! Thermo group additivity estimation: group(Cds-CdsHH) + group(Cds-Cds(Cds-Cds)H) + group(Cds-Cds(Cds-Cds)H) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + +! group(Cds-CdsCsH) + group(Cs-(Cds-Cds)CsHH) + group(Cs-(Cds-Cds)(Cds-Cds)CsH) + ring(1,3-Cyclohexadiene) + radical(C=CCJC=C-cyclohexadiene) +i6 H 9C 8 G 100.000 5000.000 1368.12 1 + 1.00899997E+01 1.51499493E-02 7.45552368E-06-2.56904073E-09 2.17408157E-13 2 + 3.10380015E+04-1.86048699E+01-6.57652868E-01 8.92974260E-02-1.20682400E-04 3 + 8.26969458E-08-1.95345474E-11 2.99803671E+04 2.20020859E+01 4 + +! Thermo group additivity estimation: group(Cds-Cds(Cds-Cds)H) + group(Cds-Cds(Cds-Cds)H) + group(Cds-Cds(Cds-Cds)H) + group(Cds-CdsCsH) + group(Cds- +! CdsCsH) + group(Cds-Cds(Cds-Cds)Cs) + group(Cs-(Cds-Cds)HHH) + group(Cs-(Cds-Cds)(Cds-Cds)HH) + ring(13cyclohexadiene5methylene) + radical(C=CC=CCJ) +i7 H 9C 8 G 100.000 5000.000 965.91 1 + 1.49840829E+01 3.00914122E-02-1.04911179E-05 1.95037002E-09-1.43618297E-13 2 + 2.03234173E+04-5.36158193E+01 1.51585120E+00 3.36849031E-02 6.49627950E-05 3 +-1.06057678E-07 4.22875321E-11 2.53594006E+04 2.34938970E+01 4 + +! Thermo group additivity estimation: group(Cds-CdsHH) + group(Cds-Cds(Cds-Cds)H) + group(Cds-Cds(Cds-Cds)H) + group(Cds-Cds(Cds-Cds)H) + group(Cds- +! CdsCsH) + group(Cds-Cds(Cds-Cds)Cs) + group(Cs-(Cds-Cds)CsHH) + group(Cs-(Cds-Cds)CsHH) + ring(1,3-Cyclohexadiene) + radical(C=CCJC=C-cyclohexadiene) +i8 H 9C 8 G 100.000 5000.000 1476.19 1 + 1.17991106E+01 1.09901877E-02 1.02947402E-05-3.09339181E-09 2.49280097E-13 2 + 2.85785491E+04-2.65800665E+01-1.38679798E+00 9.29739482E-02-1.20011913E-04 3 + 7.69809924E-08-1.69064929E-11 2.74318527E+04 2.50998255E+01 4 + +! Thermo group additivity estimation: group(Cds-Cds(Cds-Cds)H) + group(Cds-Cds(Cds-Cds)H) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + +! group(Cds-Cds(Cds-Cds)(Cds-Cds)) + group(Cs-(Cds-Cds)HHH) + group(Cs-(Cds-Cds)(Cds-Cds)HH) + ring(14cyclohexadiene3methylene) + radical(C=CC=CCJ) +i9 H 9C 8 G 100.000 5000.000 966.00 1 + 1.48486579E+01 3.02522264E-02-1.05601334E-05 1.96473718E-09-1.44728549E-13 2 + 2.38763847E+04-5.35888431E+01 1.55095799E+00 3.27961899E-02 6.70410574E-05 3 +-1.07871615E-07 4.28464394E-11 2.88959130E+04 2.27881998E+01 4 + +! primaryThermoLibrary +He He 1 G 100.000 5000.000 3142.77 1 + 2.50000000E+00 3.32026413E-12-1.37161245E-15 2.48094222E-19-1.65799619E-23 2 +-7.45374998E+02 9.28724018E-01 2.50000000E+00 9.62838494E-15-1.50167895E-17 3 + 7.73923008E-21-1.23210440E-24-7.45375000E+02 9.28724000E-01 4 + +END + + + +REACTIONS KCAL/MOLE MOLES + +C6H5+C2H4=C2H3+C6H6 8.156e-03 4.488 4.420 + +i2(+M)=i1(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ -6.278e+00 1.679e-01 -7.095e-02 1.265e-02 / + CHEB/ 1.205e+01 3.057e-01 -1.262e-01 2.032e-02 / + CHEB/ 5.054e-01 2.315e-01 -8.872e-02 9.278e-03 / + CHEB/ -2.317e-01 1.477e-01 -4.884e-02 -5.245e-04 / + CHEB/ -3.734e-02 8.072e-02 -2.025e-02 -5.082e-03 / + CHEB/ -2.160e-02 3.852e-02 -5.418e-03 -5.102e-03 / + +i3(+M)=i1(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 4.520e+00 3.354e-01 -1.347e-01 1.901e-02 / + CHEB/ 4.051e+00 5.840e-01 -2.243e-01 2.396e-02 / + CHEB/ 1.308e-01 3.880e-01 -1.277e-01 -2.170e-03 / + CHEB/ -1.818e-01 2.003e-01 -4.544e-02 -1.630e-02 / + CHEB/ -5.626e-02 8.272e-02 -5.123e-03 -1.439e-02 / + CHEB/ -2.322e-02 2.903e-02 4.798e-03 -6.779e-03 / + +i4(+M)=i1(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 9.423e+00 2.164e+00 -1.866e-01 -4.970e-02 / + CHEB/ 1.107e-01 5.200e-01 1.214e-01 4.322e-02 / + CHEB/ -1.811e-01 1.259e-02 -2.018e-02 2.488e-02 / + CHEB/ -4.833e-02 3.289e-02 -3.999e-02 1.099e-03 / + CHEB/ -2.664e-02 4.138e-02 -2.038e-02 -1.498e-03 / + CHEB/ -1.821e-02 2.840e-02 -8.859e-03 -7.638e-05 / + +i5(+M)=i1(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 6.961e+00 1.064e+00 -2.779e-01 1.167e-02 / + CHEB/ 2.681e+00 1.066e+00 -1.093e-01 -6.232e-02 / + CHEB/ -3.724e-01 3.862e-01 2.684e-02 -2.151e-02 / + CHEB/ -1.668e-01 1.283e-01 1.985e-02 6.260e-03 / + CHEB/ -5.796e-02 3.297e-02 3.832e-03 1.024e-02 / + CHEB/ -2.012e-02 1.295e-02 -8.162e-03 6.745e-03 / + +C2H4+C6H5(+M)=i1(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 1.014e+01 1.760e-01 -7.643e-02 1.540e-02 / + CHEB/ 1.096e+00 3.047e-01 -1.265e-01 2.077e-02 / + CHEB/ 2.714e-02 2.105e-01 -7.714e-02 4.991e-03 / + CHEB/ -4.521e-02 1.237e-01 -3.633e-02 -4.467e-03 / + CHEB/ -3.700e-02 6.492e-02 -1.301e-02 -6.599e-03 / + CHEB/ -2.166e-02 3.120e-02 -2.768e-03 -5.127e-03 / + +i3(+M)=i2(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ -7.113e-01 -1.657e+00 -1.837e-01 7.714e-02 / + CHEB/ 7.840e+00 1.223e+00 -3.203e-01 3.893e-02 / + CHEB/ 1.768e-01 7.042e-01 -7.239e-02 -5.347e-03 / + CHEB/ -9.309e-02 2.246e-01 -3.458e-02 -4.545e-02 / + CHEB/ -1.201e-02 7.433e-02 -1.653e-03 -1.892e-02 / + CHEB/ -1.951e-02 1.071e-02 -8.229e-03 -8.536e-03 / + +i4(+M)=i2(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 1.928e+00 1.243e-01 -1.627e-01 -2.323e-02 / + CHEB/ 5.680e+00 5.201e-01 -8.736e-02 5.769e-02 / + CHEB/ 1.531e-01 2.714e-01 -1.206e-01 2.968e-03 / + CHEB/ -1.328e-01 1.869e-01 -4.896e-02 -4.618e-03 / + CHEB/ 1.054e-02 8.764e-02 -1.855e-02 -8.371e-03 / + CHEB/ -5.357e-03 3.991e-02 -2.340e-03 -6.708e-03 / + +i5(+M)=i2(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 2.203e+00 -7.543e-01 -4.172e-01 -7.724e-03 / + CHEB/ 5.529e+00 1.291e+00 -9.530e-02 -4.903e-02 / + CHEB/ 3.716e-02 4.048e-01 -1.367e-02 1.894e-02 / + CHEB/ -1.108e-01 2.111e-01 -3.158e-02 1.021e-02 / + CHEB/ -2.260e-02 8.598e-02 -1.233e-02 -2.560e-03 / + CHEB/ -7.484e-03 3.333e-02 -4.060e-03 -5.749e-03 / + +C2H4+C6H5(+M)=i2(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 6.745e+00 -2.028e+00 -1.830e-01 2.657e-02 / + CHEB/ 3.061e+00 7.777e-01 -2.811e-01 2.628e-02 / + CHEB/ 2.582e-01 4.596e-01 -1.314e-01 -1.076e-02 / + CHEB/ 1.210e-02 2.086e-01 -3.224e-02 -2.190e-02 / + CHEB/ 1.502e-02 8.083e-02 2.033e-03 -1.434e-02 / + CHEB/ 2.508e-04 3.186e-02 5.366e-03 -5.546e-03 / + +i4(+M)=i3(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 2.584e+00 3.144e-01 -2.388e-01 -2.004e-02 / + CHEB/ 4.503e+00 8.348e-01 -1.889e-01 5.922e-02 / + CHEB/ -9.833e-02 4.377e-01 -1.453e-01 -1.274e-02 / + CHEB/ -1.301e-01 2.198e-01 -3.823e-02 -2.406e-02 / + CHEB/ -6.149e-04 8.229e-02 8.483e-04 -1.609e-02 / + CHEB/ -1.299e-03 2.863e-02 6.730e-03 -5.707e-03 / + +i5(+M)=i3(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 2.599e+00 -6.402e-01 -4.842e-01 8.263e-03 / + CHEB/ 4.544e+00 1.694e+00 -2.304e-01 -5.849e-02 / + CHEB/ -1.752e-01 5.888e-01 -2.815e-02 -1.228e-02 / + CHEB/ -1.382e-01 2.280e-01 -2.206e-03 -7.969e-03 / + CHEB/ -2.132e-02 7.165e-02 7.363e-03 -6.171e-03 / + CHEB/ -3.574e-03 2.313e-02 1.852e-03 -2.575e-03 / + +C2H4+C6H5(+M)=i3(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 6.795e+00 -1.703e+00 -3.088e-01 4.683e-02 / + CHEB/ 2.282e+00 1.107e+00 -3.466e-01 -8.891e-03 / + CHEB/ 1.044e-01 5.243e-01 -9.293e-02 -5.128e-02 / + CHEB/ -4.476e-04 1.977e-01 1.262e-03 -3.236e-02 / + CHEB/ 1.711e-02 7.082e-02 1.243e-02 -1.057e-02 / + CHEB/ 1.744e-03 2.754e-02 6.984e-03 -2.387e-03 / + +i5(+M)=i4(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 3.415e+00 1.128e+00 -3.593e-01 -2.407e-02 / + CHEB/ 2.944e+00 1.221e+00 -5.614e-02 -3.723e-02 / + CHEB/ -3.917e-01 3.817e-01 3.144e-02 -1.093e-02 / + CHEB/ -1.001e-01 1.233e-01 9.701e-03 1.142e-02 / + CHEB/ 2.974e-03 4.147e-02 -5.392e-03 1.043e-02 / + CHEB/ 1.358e-03 2.264e-02 -1.217e-02 3.880e-03 / + +C2H4+C6H5(+M)=i4(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 6.968e+00 1.605e-01 -1.340e-01 -9.209e-03 / + CHEB/ 1.233e+00 4.215e-01 -1.141e-01 3.872e-02 / + CHEB/ -4.120e-02 2.508e-01 -8.162e-02 1.816e-03 / + CHEB/ 1.010e-02 1.358e-01 -3.447e-02 -9.680e-03 / + CHEB/ 2.660e-02 6.662e-02 -8.606e-03 -9.539e-03 / + CHEB/ 4.799e-03 3.002e-02 4.891e-04 -6.158e-03 / + +C2H4+C6H5(+M)=i5(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 4.784e+00 -4.969e-01 -4.367e-01 -4.174e-02 / + CHEB/ 2.415e+00 9.140e-01 -3.456e-03 -1.290e-02 / + CHEB/ 1.125e-01 3.727e-01 -2.008e-02 1.470e-02 / + CHEB/ -3.277e-02 1.691e-01 -7.886e-03 -5.592e-04 / + CHEB/ 1.486e-02 5.485e-02 -2.039e-03 -3.572e-03 / + CHEB/ 3.180e-03 2.435e-02 -3.211e-04 -5.944e-03 / + +i1(+M)=C6H5C2H3+H(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ -2.715e+00 1.275e-01 -5.617e-02 1.160e-02 / + CHEB/ 1.035e+01 2.297e-01 -9.846e-02 1.809e-02 / + CHEB/ 3.876e-02 1.735e-01 -6.875e-02 8.382e-03 / + CHEB/ -5.629e-02 1.146e-01 -3.987e-02 7.019e-04 / + CHEB/ -4.325e-02 6.757e-02 -1.904e-02 -3.179e-03 / + CHEB/ -2.559e-02 3.507e-02 -6.472e-03 -4.199e-03 / + +i2(+M)=C6H5C2H3+H(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ -6.660e+00 8.943e-02 -4.331e-02 1.078e-02 / + CHEB/ 1.342e+01 1.647e-01 -7.871e-02 1.872e-02 / + CHEB/ 4.959e-02 1.283e-01 -5.870e-02 1.179e-02 / + CHEB/ -3.627e-02 8.401e-02 -3.518e-02 4.341e-03 / + CHEB/ -2.835e-02 4.586e-02 -1.619e-02 -6.028e-04 / + CHEB/ -1.465e-02 2.077e-02 -5.121e-03 -2.266e-03 / + +i3(+M)=C6H5C2H3+H(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ -1.498e+00 -1.799e+00 -2.673e-01 3.650e-02 / + CHEB/ 8.721e+00 1.034e+00 -3.355e-01 1.490e-03 / + CHEB/ 1.055e-01 5.100e-01 -1.009e-01 -4.421e-02 / + CHEB/ -2.628e-02 1.994e-01 -5.147e-03 -3.080e-02 / + CHEB/ 7.188e-03 7.577e-02 8.064e-03 -9.697e-03 / + CHEB/ -2.399e-03 3.109e-02 4.757e-03 -1.898e-03 / + +i4(+M)=C6H5C2H3+H(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 7.516e-01 9.729e-02 -1.188e-01 -1.553e-02 / + CHEB/ 7.230e+00 3.499e-01 -7.953e-02 4.181e-02 / + CHEB/ -7.376e-02 2.134e-01 -7.795e-02 6.983e-03 / + CHEB/ -2.357e-02 1.311e-01 -4.231e-02 -3.547e-03 / + CHEB/ 1.222e-02 7.315e-02 -1.663e-02 -6.210e-03 / + CHEB/ -1.523e-03 3.524e-02 -3.321e-03 -5.835e-03 / + +i5(+M)=C6H5C2H3+H(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 1.507e+00 -6.809e-01 -4.142e-01 -2.365e-02 / + CHEB/ 6.458e+00 9.730e-01 7.907e-03 -3.891e-02 / + CHEB/ 7.707e-02 3.211e-01 -5.213e-03 2.434e-02 / + CHEB/ -7.142e-02 1.649e-01 -1.925e-02 8.278e-03 / + CHEB/ -3.003e-03 6.393e-02 -1.043e-02 -4.446e-04 / + CHEB/ -3.335e-03 2.965e-02 -4.796e-03 -5.820e-03 / + +i6(+M)=C6H5C2H3+H(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 1.314e+00 3.966e-01 -1.179e-01 3.458e-03 / + CHEB/ 7.106e+00 6.549e-01 -1.777e-01 -4.885e-03 / + CHEB/ -1.653e-01 3.663e-01 -6.690e-02 -2.065e-02 / + CHEB/ -1.175e-01 1.341e-01 3.223e-03 -1.905e-02 / + CHEB/ -3.739e-02 3.057e-02 1.641e-02 -7.436e-03 / + CHEB/ -1.319e-02 6.151e-03 8.011e-03 -3.797e-04 / + +i7(+M)=C6H5C2H3+H(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ -2.939e-01 3.408e-01 -1.078e-01 5.609e-03 / + CHEB/ 8.839e+00 5.791e-01 -1.717e-01 1.860e-03 / + CHEB/ -1.733e-01 3.499e-01 -7.953e-02 -1.365e-02 / + CHEB/ -1.177e-01 1.398e-01 -7.868e-03 -1.757e-02 / + CHEB/ -3.637e-02 2.784e-02 1.599e-02 -9.764e-03 / + CHEB/ -7.315e-03 -1.716e-03 1.039e-02 -1.414e-03 / + +i8(+M)=C6H5C2H3+H(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 1.159e+00 3.857e-01 -1.209e-01 5.791e-03 / + CHEB/ 7.668e+00 6.367e-01 -1.833e-01 -1.291e-03 / + CHEB/ -1.716e-01 3.539e-01 -7.046e-02 -1.946e-02 / + CHEB/ -1.064e-01 1.242e-01 3.019e-03 -1.973e-02 / + CHEB/ -3.078e-02 2.281e-02 1.761e-02 -8.175e-03 / + CHEB/ -9.831e-03 3.214e-03 7.929e-03 -3.410e-04 / + +i9(+M)=C6H5C2H3+H(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ -6.135e-01 2.976e-01 -1.024e-01 8.157e-03 / + CHEB/ 9.009e+00 5.109e-01 -1.669e-01 7.470e-03 / + CHEB/ -1.284e-01 3.191e-01 -8.493e-02 -8.496e-03 / + CHEB/ -1.094e-01 1.373e-01 -1.685e-02 -1.525e-02 / + CHEB/ -3.694e-02 3.444e-02 1.041e-02 -1.041e-02 / + CHEB/ -9.261e-03 2.730e-03 9.133e-03 -2.879e-03 / + +C2H4+C6H5(+M)=C6H5C2H3+H(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 6.435e+00 -2.199e+00 -1.217e-01 1.970e-02 / + CHEB/ 3.742e+00 4.830e-01 -1.788e-01 1.675e-02 / + CHEB/ 4.195e-01 2.796e-01 -7.856e-02 -1.040e-02 / + CHEB/ 6.611e-02 1.380e-01 -2.118e-02 -1.683e-02 / + CHEB/ 3.216e-02 6.435e-02 -3.435e-04 -1.241e-02 / + CHEB/ 6.731e-03 2.484e-02 6.690e-03 -7.443e-03 / + +i6(+M)=C6H6+C2H3(+M) 1.000e+00 0.000 0.000 + TCHEB/ 298.000 1500.000 / + PCHEB/ 0.001 131.579 / + CHEB/ 6 4/ + CHEB/ 1.115e+00 4.362e-01 -1.335e-01 6.663e-03 / + CHEB/ 7.898e+00 6.921e-01 -1.869e-01 -6.175e-03 / + CHEB/ -2.637e-01 3.503e-01 -5.343e-02 -2.639e-02 / + CHEB/ -1.110e-01 1.115e-01 1.358e-02 -2.023e-02 / + CHEB/ -3.564e-02 2.386e-02 1.686e-02 -5.660e-03 / + CHEB/ -1.408e-02 7.779e-03 6.060e-03 3.437e-04 / + +END + diff --git a/rmgpy/tools/data/sim/mbSampled/input.py b/rmgpy/tools/data/sim/mbSampled/input.py new file mode 100644 index 0000000000..a7355b2370 --- /dev/null +++ b/rmgpy/tools/data/sim/mbSampled/input.py @@ -0,0 +1,87 @@ +# Data sources +database( + thermoLibraries=['BurkeH2O2', 'primaryThermoLibrary', 'DFT_QCI_thermo', 'CBS_QB3_1dHR', 'FFCM1(-)'], + reactionLibraries=[('BurkeH2O2inArHe', False), ('FFCM1(-)', False)], + seedMechanisms=[], + kineticsDepositories=['training'], + kineticsFamilies='default', + kineticsEstimator='rate rules', +) + +generatedSpeciesConstraints( + allowed=['seed mechanisms', 'input species', 'reaction libraries'], + maximumOxygenAtoms=4, + maximumCarbonAtoms=18, + maximumRadicalElectrons=1, + allowSingletO2=False, +) + +# List of species +species( + label='C6H5', + reactive=True, + structure=adjacencyList( + """ + multiplicity 2 + 1 C u0 p0 c0 {2,B} {3,B} {8,S} + 2 C u0 p0 c0 {1,B} {4,B} {7,S} + 3 C u0 p0 c0 {1,B} {5,B} {9,S} + 4 C u0 p0 c0 {2,B} {6,B} {10,S} + 5 C u0 p0 c0 {3,B} {6,B} {11,S} + 6 C u1 p0 c0 {4,B} {5,B} + 7 H u0 p0 c0 {2,S} + 8 H u0 p0 c0 {1,S} + 9 H u0 p0 c0 {3,S} + 10 H u0 p0 c0 {4,S} + 11 H u0 p0 c0 {5,S} + """), +) + +species( + label='C2H4', + reactive=True, + structure=SMILES('C=C'), +) + +species( + label='He', + reactive=False, + structure=SMILES("[He]"), +) + +mbsampledReactor( + temperature=(796, 'K'), + pressure=(10, 'torr'), + initialMoleFractions={ + "C6H5": 0.7e+12, + "C2H4": 2.68e+16, + "He": 9.38e+16, + }, + mbsamplingRate=3500, + terminationTime=(0.01, 's'), + constantSpecies=['C2H4', 'He'] +) + +simulator( + atol=1e-10, + rtol=1e-6, + sens_atol=1e-6, + sens_rtol=1e-4, +) + +model( + toleranceKeepInEdge=0.01, + toleranceMoveToCore=0.1, + toleranceInterruptSimulation=1E8, + maximumEdgeSpecies=100000, + filterReactions=True, +) + +options( + units='si', + saveRestartPeriod=None, + saveSimulationProfiles=True, + generateOutputHTML=False, + generatePlots=False, + saveEdgeSpecies=False, +) diff --git a/rmgpy/tools/data/sim/mbSampled/species_dictionary.txt b/rmgpy/tools/data/sim/mbSampled/species_dictionary.txt new file mode 100644 index 0000000000..1a2efe9530 --- /dev/null +++ b/rmgpy/tools/data/sim/mbSampled/species_dictionary.txt @@ -0,0 +1,249 @@ +C2H4 +1 C u0 p0 c0 {2,D} {3,S} {4,S} +2 C u0 p0 c0 {1,D} {5,S} {6,S} +3 H u0 p0 c0 {1,S} +4 H u0 p0 c0 {1,S} +5 H u0 p0 c0 {2,S} +6 H u0 p0 c0 {2,S} + +C6H5 +multiplicity 2 +1 C u0 p0 c0 {2,D} {3,S} {8,S} +2 C u0 p0 c0 {1,D} {4,S} {7,S} +3 C u0 p0 c0 {1,S} {5,D} {9,S} +4 C u0 p0 c0 {2,S} {6,D} {10,S} +5 C u0 p0 c0 {3,D} {6,S} {11,S} +6 C u1 p0 c0 {4,D} {5,S} +7 H u0 p0 c0 {2,S} +8 H u0 p0 c0 {1,S} +9 H u0 p0 c0 {3,S} +10 H u0 p0 c0 {4,S} +11 H u0 p0 c0 {5,S} + +C2H3 +multiplicity 2 +1 C u0 p0 c0 {2,D} {3,S} {4,S} +2 C u1 p0 c0 {1,D} {5,S} +3 H u0 p0 c0 {1,S} +4 H u0 p0 c0 {1,S} +5 H u0 p0 c0 {2,S} + +C6H6 +1 C u0 p0 c0 {2,D} {6,S} {7,S} +2 C u0 p0 c0 {1,D} {3,S} {8,S} +3 C u0 p0 c0 {2,S} {4,D} {9,S} +4 C u0 p0 c0 {3,D} {5,S} {10,S} +5 C u0 p0 c0 {4,S} {6,D} {11,S} +6 C u0 p0 c0 {1,S} {5,D} {12,S} +7 H u0 p0 c0 {1,S} +8 H u0 p0 c0 {2,S} +9 H u0 p0 c0 {3,S} +10 H u0 p0 c0 {4,S} +11 H u0 p0 c0 {5,S} +12 H u0 p0 c0 {6,S} + +C6H5C2H3 +1 C u0 p0 c0 {2,S} {3,D} {4,S} +2 C u0 p0 c0 {1,S} {5,D} {10,S} +3 C u0 p0 c0 {1,D} {7,S} {14,S} +4 C u0 p0 c0 {1,S} {8,D} {9,S} +5 C u0 p0 c0 {2,D} {6,S} {11,S} +6 C u0 p0 c0 {5,S} {7,D} {12,S} +7 C u0 p0 c0 {3,S} {6,D} {13,S} +8 C u0 p0 c0 {4,D} {15,S} {16,S} +9 H u0 p0 c0 {4,S} +10 H u0 p0 c0 {2,S} +11 H u0 p0 c0 {5,S} +12 H u0 p0 c0 {6,S} +13 H u0 p0 c0 {7,S} +14 H u0 p0 c0 {3,S} +15 H u0 p0 c0 {8,S} +16 H u0 p0 c0 {8,S} + +H +multiplicity 2 +1 H u1 p0 c0 + +i1 +multiplicity 2 +1 C u0 p0 c0 {2,S} {8,S} {9,S} {10,S} +2 C u0 p0 c0 {1,S} {3,S} {4,D} +3 C u0 p0 c0 {2,S} {5,D} {11,S} +4 C u0 p0 c0 {2,D} {7,S} {15,S} +5 C u0 p0 c0 {3,D} {6,S} {12,S} +6 C u0 p0 c0 {5,S} {7,D} {13,S} +7 C u0 p0 c0 {4,S} {6,D} {14,S} +8 C u1 p0 c0 {1,S} {16,S} {17,S} +9 H u0 p0 c0 {1,S} +10 H u0 p0 c0 {1,S} +11 H u0 p0 c0 {3,S} +12 H u0 p0 c0 {5,S} +13 H u0 p0 c0 {6,S} +14 H u0 p0 c0 {7,S} +15 H u0 p0 c0 {4,S} +16 H u0 p0 c0 {8,S} +17 H u0 p0 c0 {8,S} + +i2 +multiplicity 2 +1 C u0 p0 c0 {3,S} {9,S} {10,S} {11,S} +2 C u0 p0 c0 {3,S} {4,S} {5,D} +3 C u1 p0 c0 {1,S} {2,S} {12,S} +4 C u0 p0 c0 {2,S} {6,D} {13,S} +5 C u0 p0 c0 {2,D} {8,S} {17,S} +6 C u0 p0 c0 {4,D} {7,S} {14,S} +7 C u0 p0 c0 {6,S} {8,D} {15,S} +8 C u0 p0 c0 {5,S} {7,D} {16,S} +9 H u0 p0 c0 {1,S} +10 H u0 p0 c0 {1,S} +11 H u0 p0 c0 {1,S} +12 H u0 p0 c0 {3,S} +13 H u0 p0 c0 {4,S} +14 H u0 p0 c0 {6,S} +15 H u0 p0 c0 {7,S} +16 H u0 p0 c0 {8,S} +17 H u0 p0 c0 {5,S} + +i3 +multiplicity 2 +1 C u0 p0 c0 {2,S} {3,S} {9,S} {10,S} +2 C u0 p0 c0 {1,S} {11,S} {12,S} {13,S} +3 C u0 p0 c0 {1,S} {4,D} {8,S} +4 C u0 p0 c0 {3,D} {5,S} {16,S} +5 C u0 p0 c0 {4,S} {6,D} {15,S} +6 C u0 p0 c0 {5,D} {7,S} {14,S} +7 C u0 p0 c0 {6,S} {8,D} {17,S} +8 C u1 p0 c0 {3,S} {7,D} +9 H u0 p0 c0 {1,S} +10 H u0 p0 c0 {1,S} +11 H u0 p0 c0 {2,S} +12 H u0 p0 c0 {2,S} +13 H u0 p0 c0 {2,S} +14 H u0 p0 c0 {6,S} +15 H u0 p0 c0 {5,S} +16 H u0 p0 c0 {4,S} +17 H u0 p0 c0 {7,S} + +i4 +multiplicity 2 +1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S} +2 C u0 p0 c0 {1,S} {3,S} {9,S} {10,S} +3 C u0 p0 c0 {1,S} {2,S} {11,S} {12,S} +4 C u0 p0 c0 {1,S} {6,D} {13,S} +5 C u0 p0 c0 {1,S} {7,D} {14,S} +6 C u0 p0 c0 {4,D} {8,S} {16,S} +7 C u0 p0 c0 {5,D} {8,S} {17,S} +8 C u1 p0 c0 {6,S} {7,S} {15,S} +9 H u0 p0 c0 {2,S} +10 H u0 p0 c0 {2,S} +11 H u0 p0 c0 {3,S} +12 H u0 p0 c0 {3,S} +13 H u0 p0 c0 {4,S} +14 H u0 p0 c0 {5,S} +15 H u0 p0 c0 {8,S} +16 H u0 p0 c0 {6,S} +17 H u0 p0 c0 {7,S} + +i5 +multiplicity 2 +1 C u0 p0 c0 {2,S} {4,S} {5,S} {9,S} +2 C u0 p0 c0 {1,S} {3,S} {10,S} {11,S} +3 C u0 p0 c0 {2,S} {4,S} {12,S} {13,S} +4 C u0 p0 c0 {1,S} {3,S} {6,D} +5 C u0 p0 c0 {1,S} {7,D} {14,S} +6 C u0 p0 c0 {4,D} {8,S} {17,S} +7 C u0 p0 c0 {5,D} {8,S} {16,S} +8 C u1 p0 c0 {6,S} {7,S} {15,S} +9 H u0 p0 c0 {1,S} +10 H u0 p0 c0 {2,S} +11 H u0 p0 c0 {2,S} +12 H u0 p0 c0 {3,S} +13 H u0 p0 c0 {3,S} +14 H u0 p0 c0 {5,S} +15 H u0 p0 c0 {8,S} +16 H u0 p0 c0 {7,S} +17 H u0 p0 c0 {6,S} + +i6 +multiplicity 2 +1 C u0 p0 c0 {2,S} {3,S} {4,S} {9,S} +2 C u0 p0 c0 {1,S} {5,D} {11,S} +3 C u0 p0 c0 {1,S} {6,D} {12,S} +4 C u0 p0 c0 {1,S} {8,D} {10,S} +5 C u0 p0 c0 {2,D} {7,S} {14,S} +6 C u0 p0 c0 {3,D} {7,S} {15,S} +7 C u1 p0 c0 {5,S} {6,S} {13,S} +8 C u0 p0 c0 {4,D} {16,S} {17,S} +9 H u0 p0 c0 {1,S} +10 H u0 p0 c0 {4,S} +11 H u0 p0 c0 {2,S} +12 H u0 p0 c0 {3,S} +13 H u0 p0 c0 {7,S} +14 H u0 p0 c0 {5,S} +15 H u0 p0 c0 {6,S} +16 H u0 p0 c0 {8,S} +17 H u0 p0 c0 {8,S} + +i7 +multiplicity 2 +1 C u0 p0 c0 {2,S} {3,S} {9,S} {10,S} +2 C u0 p0 c0 {1,S} {4,D} {5,S} +3 C u0 p0 c0 {1,S} {6,D} {11,S} +4 C u0 p0 c0 {2,D} {7,S} {14,S} +5 C u0 p0 c0 {2,S} {8,D} {13,S} +6 C u0 p0 c0 {3,D} {7,S} {15,S} +7 C u1 p0 c0 {4,S} {6,S} {12,S} +8 C u0 p0 c0 {5,D} {16,S} {17,S} +9 H u0 p0 c0 {1,S} +10 H u0 p0 c0 {1,S} +11 H u0 p0 c0 {3,S} +12 H u0 p0 c0 {7,S} +13 H u0 p0 c0 {5,S} +14 H u0 p0 c0 {4,S} +15 H u0 p0 c0 {6,S} +16 H u0 p0 c0 {8,S} +17 H u0 p0 c0 {8,S} + +i8 +multiplicity 2 +1 C u0 p0 c0 {3,S} {4,S} {9,S} {10,S} +2 C u0 p0 c0 {3,D} {5,S} {7,S} +3 C u0 p0 c0 {1,S} {2,D} {12,S} +4 C u0 p0 c0 {1,S} {6,D} {11,S} +5 C u1 p0 c0 {2,S} {6,S} {13,S} +6 C u0 p0 c0 {4,D} {5,S} {15,S} +7 C u0 p0 c0 {2,S} {8,D} {14,S} +8 C u0 p0 c0 {7,D} {16,S} {17,S} +9 H u0 p0 c0 {1,S} +10 H u0 p0 c0 {1,S} +11 H u0 p0 c0 {4,S} +12 H u0 p0 c0 {3,S} +13 H u0 p0 c0 {5,S} +14 H u0 p0 c0 {7,S} +15 H u0 p0 c0 {6,S} +16 H u0 p0 c0 {8,S} +17 H u0 p0 c0 {8,S} + +i9 +multiplicity 2 +1 C u0 p0 c0 {3,S} {4,S} {9,S} {10,S} +2 C u1 p0 c0 {5,S} {6,S} {7,S} +3 C u0 p0 c0 {1,S} {5,D} {11,S} +4 C u0 p0 c0 {1,S} {6,D} {12,S} +5 C u0 p0 c0 {2,S} {3,D} {14,S} +6 C u0 p0 c0 {2,S} {4,D} {15,S} +7 C u0 p0 c0 {2,S} {8,D} {13,S} +8 C u0 p0 c0 {7,D} {16,S} {17,S} +9 H u0 p0 c0 {1,S} +10 H u0 p0 c0 {1,S} +11 H u0 p0 c0 {3,S} +12 H u0 p0 c0 {4,S} +13 H u0 p0 c0 {7,S} +14 H u0 p0 c0 {5,S} +15 H u0 p0 c0 {6,S} +16 H u0 p0 c0 {8,S} +17 H u0 p0 c0 {8,S} + +He +1 He u0 p1 c0 + diff --git a/rmgpy/tools/simulateTest.py b/rmgpy/tools/simulateTest.py index b9190da7b2..7ad3975ece 100644 --- a/rmgpy/tools/simulateTest.py +++ b/rmgpy/tools/simulateTest.py @@ -39,7 +39,8 @@ class SimulateTest(unittest.TestCase): def test_minimal(self): - folder = os.path.join(os.path.dirname(rmgpy.__file__), 'tools/data/sim/simple') + """Test that we can simlulate a SimpleReactor with sensitivity""" + folder = os.path.join(os.path.dirname(rmgpy.__file__), 'tools', 'data', 'sim', 'simple') inputFile = os.path.join(folder, 'input.py') chemkinFile = os.path.join(folder, 'chem.inp') @@ -57,7 +58,8 @@ def test_minimal(self): os.remove(os.path.join(folder, 'simulate.log')) def test_liquid(self): - folder = os.path.join(os.path.dirname(rmgpy.__file__), 'tools/data/sim/liquid') + """Test that we can simulate a LiquidReactor with sensitivity""" + folder = os.path.join(os.path.dirname(rmgpy.__file__), 'tools', 'data', 'sim', 'liquid') inputFile = os.path.join(folder, 'input.py') chemkinFile = os.path.join(folder, 'chem.inp') @@ -74,6 +76,23 @@ def test_liquid(self): shutil.rmtree(os.path.join(folder, 'solver')) os.remove(os.path.join(folder, 'simulate.log')) + def test_mbSampled(self): + """Test that we can simulate an MBSampledReactor""" + folder = os.path.join(os.path.dirname(rmgpy.__file__), 'tools', 'data', 'sim', 'mbSampled') + + inputFile = os.path.join(folder, 'input.py') + chemkinFile = os.path.join(folder, 'chem.inp') + dictFile = os.path.join(folder, 'species_dictionary.txt') + + run_simulation(inputFile, chemkinFile, dictFile) + + simfile = os.path.join(folder, 'solver', 'simulation_1_30.csv') + + self.assertTrue(os.path.isfile(simfile)) + + shutil.rmtree(os.path.join(folder, 'solver')) + os.remove(os.path.join(folder, 'simulate.log')) + def tearDown(self): import rmgpy.data.rmg rmgpy.data.rmg.database = None From 206f149dd939cee5ae325a148702f8240bc2508b Mon Sep 17 00:00:00 2001 From: Max Liu Date: Sun, 4 Aug 2019 22:25:10 -0400 Subject: [PATCH 3/3] Disable logging in rmgpy.tools.simulateTest --- rmgpy/tools/simulateTest.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/rmgpy/tools/simulateTest.py b/rmgpy/tools/simulateTest.py index 7ad3975ece..d855ffb787 100644 --- a/rmgpy/tools/simulateTest.py +++ b/rmgpy/tools/simulateTest.py @@ -28,6 +28,7 @@ # # ############################################################################### +import logging import unittest import os import os.path @@ -36,8 +37,14 @@ from rmgpy.tools.simulate import run_simulation import rmgpy + class SimulateTest(unittest.TestCase): + def setUp(self): + """This method is run once before each unit test""" + # Disable logging + logging.disable(logging.WARNING) + def test_minimal(self): """Test that we can simlulate a SimpleReactor with sensitivity""" folder = os.path.join(os.path.dirname(rmgpy.__file__), 'tools', 'data', 'sim', 'simple') @@ -96,3 +103,6 @@ def test_mbSampled(self): def tearDown(self): import rmgpy.data.rmg rmgpy.data.rmg.database = None + + # Reset logging + logging.disable(logging.NOTSET)