diff --git a/documentation/source/users/rmg/releaseNotes.rst b/documentation/source/users/rmg/releaseNotes.rst index cc55fca45e..7134c58430 100644 --- a/documentation/source/users/rmg/releaseNotes.rst +++ b/documentation/source/users/rmg/releaseNotes.rst @@ -4,6 +4,22 @@ Release Notes ************* +RMG-Py Version 2.2.1 +==================== +Date July 23, 2018 + +This release is minor patch which fixes a number of issues discovered after 2.2.0. + +- Collision limit checking: + - RMG will now output a list of collision limit violations for the generated model + +- Fixes: + - Ambiguous chemical formulas in SMILES lookup leading to incorrect SMILES generation + - Fixed issue with reading geometries from QChem output files + - React flags for reaction filter were not properly updated on each iteration + - Fixed issue with inconsistent symmetry number calculation + + RMG-Py Version 2.2.0 ==================== Date: July 5, 2018 diff --git a/rmgpy/cantherm/qchem.py b/rmgpy/cantherm/qchem.py index f5227a4ce0..a975b41bd3 100644 --- a/rmgpy/cantherm/qchem.py +++ b/rmgpy/cantherm/qchem.py @@ -33,6 +33,7 @@ import logging import os.path import rmgpy.constants as constants +from rmgpy.exceptions import InputError from rmgpy.cantherm.common import checkConformerEnergy from rmgpy.statmech import IdealGasTranslation, NonlinearRotor, LinearRotor, HarmonicOscillator, Conformer ################################################################################ @@ -115,35 +116,41 @@ def loadGeometry(self): """ atom = []; coord = []; number = []; - f = open(self.path, 'r') - line = f.readline() - while line != '': - if 'Final energy is' in line: - print 'found a sucessfully completed Qchem Geometry Optimization Job' - line = f.readline() - atom = []; coord = [] + + with open(self.path) as f: + log = f.read().splitlines() + + #First check that the Qchem job file (not necessarily a geometry optimization) + #has successfully completed, if not an error is thrown + completed_job = False + for line in reversed(log): + if 'Total job time:' in line: + logging.debug('Found a sucessfully completed Qchem Job') + completed_job = True break - line = f.readline() - found = 0 - while line != '': + + if not completed_job: + raise InputError('Could not find a successfully completed Qchem job in Qchem output file {0}'.format(self.path)) + + #Now look for the geometry. + #Will return the final geometry in the file under Standard Nuclear Orientation. + geometry_flag = False + for i in reversed(xrange(len(log))): + line = log[i] if 'Standard Nuclear Orientation' in line: - found += 1 - for i in range(3): line = f.readline() # skip lines - while '----------------------------------------------------' not in line: - data = line.split() - atom.append((data[1])) - coord.append([float(data[2]), float(data[3]), float(data[4])]) - line = f.readline() - # Read the next line in the file - line = f.readline() - # Read the next line in the file - line = f.readline() - if found ==1: break - line = f.readline() - #print coord - f.close() + atom, coord, number = [], [], [] + for line in log[(i+3):]: + if '------------' not in line: + data = line.split() + atom.append(data[1]) + coord.append([float(c) for c in data [2:]]) + geometry_flag = True + else: + break + if geometry_flag: + break + coord = numpy.array(coord, numpy.float64) - mass = numpy.array(coord, numpy.float64) # Assign appropriate mass to each atom in molecule # These values were taken from "Atomic Weights and Isotopic Compositions" v3.0 (July 2010) from NIST @@ -173,7 +180,9 @@ def loadGeometry(self): number.append('17') else: raise NotImplementedError('Atomic atom {0:d} not yet supported in loadGeometry().'.format(atom[i])) - number = numpy.array(number, numpy.int) + number = numpy.array(number, numpy.int) + if len(number) == 0 or len(coord) == 0 or len(mass) == 0: + raise InputError('Unable to read the numbers and types of atoms from Qchem output file {0}'.format(self.path)) return coord, number, mass def loadConformer(self, symmetry=None, spinMultiplicity=0, opticalIsomers=1, symfromlog=None, label=''): @@ -301,7 +310,7 @@ def loadEnergy(self,frequencyScaleFactor=1.): if E0 is not None: return E0 else: - raise Exception('Unable to find energy in Qchem output file.') + raise InputError('Unable to find energy in Qchem output file.') def loadZeroPointEnergy(self,frequencyScaleFactor=1.): """ @@ -331,7 +340,7 @@ def loadZeroPointEnergy(self,frequencyScaleFactor=1.): if ZPE is not None: return ZPE else: - raise Exception('Unable to find zero-point energy in Qchem output file.') + raise InputError('Unable to find zero-point energy in Qchem output file.') def loadScanEnergies(self): """ @@ -394,4 +403,4 @@ def loadNegativeFrequency(self): if frequency < 0: return frequency else: - raise Exception('Unable to find imaginary frequency in QChem output file {0}'.format(self.path)) + raise InputError('Unable to find imaginary frequency in QChem output file {0}'.format(self.path)) diff --git a/rmgpy/data/thermoTest.py b/rmgpy/data/thermoTest.py index f89ee910dd..194473d90b 100644 --- a/rmgpy/data/thermoTest.py +++ b/rmgpy/data/thermoTest.py @@ -424,7 +424,6 @@ def testNewThermoGeneration(self): self.assertAlmostEqual(Cp, thermoData.getHeatCapacity(T) / 4.184, places=1, msg="Cp{3} error for {0}. Expected {1} but calculated {2}.".format(smiles, Cp, thermoData.getHeatCapacity(T) / 4.184, T)) - @work_in_progress def testSymmetryNumberGeneration(self): """ Test we generate symmetry numbers correctly. diff --git a/rmgpy/molecule/vf2.pyx b/rmgpy/molecule/vf2.pyx index 44456637e2..1158e8f690 100644 --- a/rmgpy/molecule/vf2.pyx +++ b/rmgpy/molecule/vf2.pyx @@ -166,6 +166,14 @@ cdef class VF2: graph1.restore_vertex_order() graph2.restore_vertex_order() + # We're done, so clear the mappings to prevent downstream effects + for vertex1 in graph1.vertices: + vertex1.mapping = None + vertex1.terminal = False + for vertex2 in graph2.vertices: + vertex2.mapping = None + vertex2.terminal = False + cdef bint match(self, int callDepth) except -2: """ Recursively search for pairs of vertices to match, until all vertices diff --git a/rmgpy/molecule/vf2Test.py b/rmgpy/molecule/vf2Test.py index 04846cf178..7e1aff6f2a 100644 --- a/rmgpy/molecule/vf2Test.py +++ b/rmgpy/molecule/vf2Test.py @@ -81,6 +81,15 @@ def test_feasible(self): self.assertTrue(self.vf2.feasible(atom1, atom2)) else: # different connectivity values should return false self.assertFalse(self.vf2.feasible(atom1, atom2)) + + def test_clear_mapping(self): + """Test that vertex mapping is cleared after isomorphism.""" + self.vf2.isIsomorphic(self.mol, self.mol2, None) + + for atom in self.mol.atoms: + self.assertIsNone(atom.mapping) + self.assertFalse(atom.terminal) + ################################################################################ if __name__ == '__main__': diff --git a/rmgpy/reaction.pxd b/rmgpy/reaction.pxd index 8434d03b79..87bbec244b 100644 --- a/rmgpy/reaction.pxd +++ b/rmgpy/reaction.pxd @@ -113,4 +113,12 @@ cdef class Reaction: cpdef ensure_species(self, bint reactant_resonance=?, bint product_resonance=?) + cpdef list check_collision_limit_violation(self, float t_min, float t_max, float p_min, float p_max) + + cpdef calculate_coll_limit(self, float temp, bint reverse=?) + + cpdef get_reduced_mass(self, bint reverse=?) + + cpdef get_mean_sigma_and_epsilon(self, bint reverse=?) + cpdef bint isomorphic_species_lists(list list1, list list2, bint check_identical=?, bint only_check_label=?) diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index b3eed41bf9..8222263809 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -1146,6 +1146,106 @@ def ensure_species(self, reactant_resonance=False, product_resonance=True): except AttributeError: pass + def check_collision_limit_violation(self, t_min, t_max, p_min, p_max): + """ + Warn if a core reaction violates the collision limit rate in either the forward or reverse direction + at the relevant extreme T/P conditions. Assuming a monotonic behaviour of the kinetics. + Returns a list with the reaction object and the direction in which the violation was detected. + """ + conditions = [[t_min, p_min]] + if t_min != t_max: + conditions.append([t_max, p_min]) + if self.kinetics.isPressureDependent() and p_max != p_min: + conditions.append([t_min, p_max]) + if t_min != t_max: + conditions.append([t_max, p_max]) + logging.debug("Checking whether reaction {0} violates the collision rate limit...".format(self)) + violator_list = [] + kf_list = [] + kr_list = [] + collision_limit_f = [] + collision_limit_r = [] + for condition in conditions: + if len(self.reactants) >= 2: + try: + collision_limit_f.append(self.calculate_coll_limit(temp=condition[0], reverse=False)) + except ValueError: + continue + else: + kf_list.append(self.getRateCoefficient(condition[0], condition[1])) + if len(self.products) >= 2: + try: + collision_limit_r.append(self.calculate_coll_limit(temp=condition[0], reverse=True)) + except ValueError: + continue + else: + kr_list.append(self.generateReverseRateCoefficient().getRateCoefficient(condition[0], condition[1])) + if len(self.reactants) >= 2: + for i, k in enumerate(kf_list): + if k > collision_limit_f[i]: + ratio = k / collision_limit_f[i] + condition = '{0} K, {1:.1f} bar'.format(conditions[i][0], conditions[i][1]/1e5) + violator_list.append([self, 'forward', ratio, condition]) + if len(self.products) >= 2: + for i, k in enumerate(kr_list): + if k > collision_limit_r[i]: + ratio = k / collision_limit_r[i] + condition = '{0} K, {1:.1f} bar'.format(conditions[i][0], conditions[i][1]/1e5) + violator_list.append([self, 'reverse', ratio, condition]) + return violator_list + + def calculate_coll_limit(self, temp, reverse=False): + """ + Calculate the collision limit rate for the given temperature + implemented as recommended in Wang et al. doi 10.1016/j.combustflame.2017.08.005 (Eq. 1) + """ + reduced_mass = self.get_reduced_mass(reverse) + sigma, epsilon = self.get_mean_sigma_and_epsilon(reverse) + Tr = temp * constants.kB * constants.Na / epsilon + reduced_coll_integral = 1.16145 * Tr ** (-0.14874) + 0.52487 * math.exp(-0.7732 * Tr) + 2.16178 * math.exp( + -2.437887 * Tr) + k_coll = (math.sqrt(8 * math.pi * constants.kB * temp * constants.Na / reduced_mass) * sigma ** 2 + * reduced_coll_integral * constants.Na) + return k_coll + + def get_reduced_mass(self, reverse=False): + """ + Returns the reduced mass of the reactants if reverse is ``False`` + Returns the reduced mass of the products if reverse is ``True`` + """ + if reverse: + mass_list = [spc.molecule[0].getMolecularWeight() for spc in self.products] + else: + mass_list = [spc.molecule[0].getMolecularWeight() for spc in self.reactants] + reduced_mass = reduce((lambda x, y: x * y), mass_list) / sum(mass_list) + return reduced_mass + + def get_mean_sigma_and_epsilon(self, reverse=False): + """ + Calculates the collision diameter (sigma) using an arithmetic mean + Calculates the well depth (epsilon) using a geometric mean + If reverse is ``False`` the above is calculated for the reactants, otherwise for the products + """ + sigmas = [] + epsilons = [] + if reverse: + for spc in self.products: + trans = spc.getTransportData() + sigmas.append(trans.sigma.value_si) + epsilons.append(trans.epsilon.value_si) + num_of_spcs = len(self.products) + else: + for spc in self.reactants: + trans = spc.getTransportData() + sigmas.append(trans.sigma.value_si) + epsilons.append(trans.epsilon.value_si) + num_of_spcs = len(self.reactants) + if any([x == 0 for x in sigmas + epsilons]): + raise ValueError + mean_sigmas = sum(sigmas) / num_of_spcs + mean_epsilons = reduce((lambda x, y: x * y), epsilons) ** (1 / len(epsilons)) + return mean_sigmas, mean_epsilons + def isomorphic_species_lists(list1, list2, check_identical=False, only_check_label=False): """ This method compares whether lists of species or molecules are isomorphic diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index 01f1e39787..46d5647f74 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -33,19 +33,19 @@ Generator (RMG). """ -import os.path import sys -import logging import warnings import time +import logging +import os import shutil -import numpy -np = numpy + +import numpy as np import gc import copy - from copy import deepcopy from scipy.optimize import brute +from cantera import ck2cti from rmgpy.rmg.settings import ModelSettings from rmgpy.constraints import failsSpeciesConstraints @@ -53,9 +53,10 @@ from rmgpy.solver.base import TerminationTime, TerminationConversion from rmgpy.solver.simple import SimpleReactor from rmgpy.data.rmg import RMGDatabase -from rmgpy.exceptions import ForbiddenStructureException, DatabaseError +from rmgpy.exceptions import ForbiddenStructureException, DatabaseError, CoreError from rmgpy.data.kinetics.library import KineticsLibrary, LibraryReaction from rmgpy.data.kinetics.family import KineticsFamily, TemplateReaction +from rmgpy.rmg.pdep import PDepReaction from rmgpy.data.thermo import ThermoLibrary from rmgpy.data.base import Entry @@ -76,8 +77,6 @@ from rmgpy.stats import ExecutionStatsWriter from rmgpy.thermo.thermoengine import submit from rmgpy.tools.simulate import plot_sensitivity -from cantera import ck2cti -from rmgpy.exceptions import CoreError ################################################################################ solvent = None @@ -147,6 +146,10 @@ def __init__(self, inputFile=None, outputDirectory=None): self.clear() self.modelSettingsList = [] self.simulatorSettingsList = [] + self.Tmin = 0.0 + self.Tmax = 0.0 + self.Pmin = 0.0 + self.Pmax = 0.0 def clear(self): """ @@ -599,15 +602,16 @@ def execute(self, **kwargs): self.register_listeners() self.done = False - - Tlist = [] - for x in self.reactionSystems: - if x.Trange: - Tlist.append(x.Trange[1].value_si) - elif x.T: - Tlist.append(x.T.value_si) - - self.Tmax = max(Tlist) + + # determine min and max values for T and P (don't determine P values for liquid reactors) + self.Tmin = min([r_sys.Trange[0] if r_sys.Trange else r_sys.T for r_sys in self.reactionSystems]).value_si + self.Tmax = max([r_sys.Trange[1] if r_sys.Trange else r_sys.T for r_sys in self.reactionSystems]).value_si + try: + self.Pmin = min([x.Prange[0] if x.Prange else x.P for x in self.reactionSystems]).value_si + self.Pmax = max([x.Prange[1] if x.Prange else x.P for x in self.reactionSystems]).value_si + except AttributeError: + # For LiquidReactor, Pmin and Pmax remain with the default value of `None` + pass self.rmg_memories = [] @@ -641,7 +645,7 @@ def execute(self, **kwargs): self.rmg_memories[index].generate_cond() log_conditions(self.rmg_memories,index) - if not numpy.isinf(self.modelSettingsList[0].toleranceThermoKeepSpeciesInEdge): + if not np.isinf(self.modelSettingsList[0].toleranceThermoKeepSpeciesInEdge): self.reactionModel.setThermodynamicFilteringParameters(self.Tmax,toleranceThermoKeepSpeciesInEdge=self.modelSettingsList[0].toleranceThermoKeepSpeciesInEdge, minCoreSizeForPrune=self.modelSettingsList[0].minCoreSizeForPrune, maximumEdgeSpecies =self.modelSettingsList[0].maximumEdgeSpecies, @@ -652,7 +656,7 @@ def execute(self, **kwargs): bimolecularReact=self.bimolecularReact, trimolecularReact=self.trimolecularReact) - if not numpy.isinf(self.modelSettingsList[0].toleranceThermoKeepSpeciesInEdge): + if not np.isinf(self.modelSettingsList[0].toleranceThermoKeepSpeciesInEdge): self.reactionModel.thermoFilterDown(maximumEdgeSpecies=self.modelSettingsList[0].maximumEdgeSpecies) logging.info('Completed initial enlarge edge step...') @@ -762,54 +766,52 @@ def execute(self, **kwargs): for objectToEnlarge in objectsToEnlarge: self.reactionModel.enlarge(objectToEnlarge) - if len(self.reactionModel.core.species) > numCoreSpecies or self.reactionModel.iterationNum == 1: + if modelSettings.filterReactions: + # Run a raw simulation to get updated reaction system threshold values + # Run with the same conditions as with pruning off tempModelSettings = deepcopy(modelSettings) tempModelSettings.fluxToleranceKeepInEdge = 0 - # If there were core species added, then react the edge - # If there were no new core species, it means the pdep network needs be updated through another enlarge core step - if modelSettings.filterReactions: - # Run a raw simulation to get updated reaction system threshold values - # Run with the same conditions as with pruning off - if not resurrected: - try: - reactionSystem.simulate( - coreSpecies = self.reactionModel.core.species, - coreReactions = self.reactionModel.core.reactions, - edgeSpecies = [], - edgeReactions = [], - surfaceSpecies = self.reactionModel.surface.species, - surfaceReactions = self.reactionModel.surface.reactions, - pdepNetworks = self.reactionModel.networkList, - modelSettings = tempModelSettings, - simulatorSettings = simulatorSettings, - conditions = self.rmg_memories[index].get_cond() - ) - except: - self.updateReactionThresholdAndReactFlags( - rxnSysUnimolecularThreshold = reactionSystem.unimolecularThreshold, - rxnSysBimolecularThreshold = reactionSystem.bimolecularThreshold, skipUpdate=True) - logging.warn('Reaction thresholds/flags for Reaction System {0} was not updated due to simulation failure'.format(index+1)) - else: - self.updateReactionThresholdAndReactFlags( - rxnSysUnimolecularThreshold = reactionSystem.unimolecularThreshold, - rxnSysBimolecularThreshold = reactionSystem.bimolecularThreshold, - rxnSysTrimolecularThreshold = reactionSystem.trimolecularThreshold - ) - else: + if not resurrected: + try: + reactionSystem.simulate( + coreSpecies = self.reactionModel.core.species, + coreReactions = self.reactionModel.core.reactions, + edgeSpecies = [], + edgeReactions = [], + surfaceSpecies = self.reactionModel.surface.species, + surfaceReactions = self.reactionModel.surface.reactions, + pdepNetworks = self.reactionModel.networkList, + modelSettings = tempModelSettings, + simulatorSettings = simulatorSettings, + conditions = self.rmg_memories[index].get_cond() + ) + except: self.updateReactionThresholdAndReactFlags( rxnSysUnimolecularThreshold = reactionSystem.unimolecularThreshold, rxnSysBimolecularThreshold = reactionSystem.bimolecularThreshold, rxnSysTrimolecularThreshold = reactionSystem.trimolecularThreshold, - skipUpdate = True + skipUpdate=True) + logging.warn('Reaction thresholds/flags for Reaction System {0} was not updated due to simulation failure'.format(index+1)) + else: + self.updateReactionThresholdAndReactFlags( + rxnSysUnimolecularThreshold = reactionSystem.unimolecularThreshold, + rxnSysBimolecularThreshold = reactionSystem.bimolecularThreshold, + rxnSysTrimolecularThreshold = reactionSystem.trimolecularThreshold ) - logging.warn('Reaction thresholds/flags for Reaction System {0} was not updated due to resurrection'.format(index+1)) - - - logging.info('') else: - self.updateReactionThresholdAndReactFlags() - - if not numpy.isinf(modelSettings.toleranceThermoKeepSpeciesInEdge): + self.updateReactionThresholdAndReactFlags( + rxnSysUnimolecularThreshold = reactionSystem.unimolecularThreshold, + rxnSysBimolecularThreshold = reactionSystem.bimolecularThreshold, + rxnSysTrimolecularThreshold = reactionSystem.trimolecularThreshold, + skipUpdate = True + ) + logging.warn('Reaction thresholds/flags for Reaction System {0} was not updated due to resurrection'.format(index+1)) + + logging.info('') + else: + self.updateReactionThresholdAndReactFlags() + + if not np.isinf(modelSettings.toleranceThermoKeepSpeciesInEdge): self.reactionModel.setThermodynamicFilteringParameters(self.Tmax, toleranceThermoKeepSpeciesInEdge=modelSettings.toleranceThermoKeepSpeciesInEdge, minCoreSizeForPrune=modelSettings.minCoreSizeForPrune, maximumEdgeSpecies=modelSettings.maximumEdgeSpecies, @@ -825,7 +827,7 @@ def execute(self, **kwargs): if oldEdgeSize != len(self.reactionModel.edge.reactions) or oldCoreSize != len(self.reactionModel.core.reactions): reactorDone = False - if not numpy.isinf(self.modelSettingsList[0].toleranceThermoKeepSpeciesInEdge): + if not np.isinf(self.modelSettingsList[0].toleranceThermoKeepSpeciesInEdge): self.reactionModel.thermoFilterDown(maximumEdgeSpecies=modelSettings.maximumEdgeSpecies) maxNumSpcsHit = len(self.reactionModel.core.species) >= modelSettings.maxNumSpecies @@ -952,6 +954,53 @@ def check_model(self): '\n{2}\n{3}'.format(spc.label, spc2.label, spc.toAdjacencyList(), spc2.toAdjacencyList()) ) + # Check all core reactions (in both directions) for collision limit violation + violators = [] + num_rxn_violators = 0 + for rxn in self.reactionModel.core.reactions: + violator_list = rxn.check_collision_limit_violation(t_min=self.Tmin, t_max=self.Tmax, + p_min=self.Pmin, p_max=self.Pmax) + if violator_list: + violators.extend(violator_list) + num_rxn_violators += 1 + # Whether or not violators were found, rename 'collision_rate_violators.log' if it exists + if os.path.isfile('collision_rate_violators.log'): + # If there are no violators, yet the violators log exists (probably from a previous run + # in the same folder), rename it. + if os.path.isfile('collision_rate_violators_OLD.log'): + os.remove('collision_rate_violators_OLD.log') + os.rename('collision_rate_violators.log', 'collision_rate_violators_OLD.log') + if violators: + logging.info("\n") + logging.warning("{0} CORE reactions violate the collision rate limit!" + "\nSee the 'collision_rate_violators.log' for details.\n\n".format(num_rxn_violators)) + with open('collision_rate_violators.log', 'w') as violators_f: + violators_f.write('*** Collision rate limit violators report ***\n' + '"Violation factor" is the ratio of the rate coefficient to the collision limit' + ' rate at the relevant conditions\n\n') + for violator in violators: + rxn_string = str(violator[0]) + kinetics = violator[0].kinetics + comment='' + if isinstance(violator[0], TemplateReaction): + comment = violator[0].kinetics.comment + violator[0].kinetics.comment = '' # the comment is printed better when outside of the object + if isinstance(violator[0], LibraryReaction): + comment = 'Kinetic library: {0}'.format(violator[0].library) + if isinstance(violator[0], PDepReaction): + comment = 'Network #{0}'.format(violator[0].network) + direction = violator[1] + ratio = violator[2] + condition = violator[3] + violators_f.write('{0}\n{1}\n{2}\nDirection: {3}\nViolation factor: {4:.2f}\n' + 'Violation condition: {5}\n\n'.format( + rxn_string, kinetics, comment, direction, ratio, condition)) + if isinstance(violator[0], TemplateReaction): + # although this is the end of the run, restore the original comment + violator[0].kinetics.comment = comment + else: + logging.info("No collision rate violators found.") + def makeSeedMech(self,firstTime=False): """ causes RMG to make a seed mechanism out of the current chem_annotated.inp and species_dictionary.txt @@ -1169,19 +1218,19 @@ def generateCanteraFiles(self, chemkinFile, **kwargs): def initializeReactionThresholdAndReactFlags(self): numCoreSpecies = len(self.reactionModel.core.species) if self.filterReactions: - self.unimolecularReact = numpy.zeros((numCoreSpecies),bool) - self.bimolecularReact = numpy.zeros((numCoreSpecies, numCoreSpecies),bool) - self.unimolecularThreshold = numpy.zeros((numCoreSpecies),bool) - self.bimolecularThreshold = numpy.zeros((numCoreSpecies, numCoreSpecies),bool) + self.unimolecularReact = np.zeros((numCoreSpecies),bool) + self.bimolecularReact = np.zeros((numCoreSpecies, numCoreSpecies),bool) + self.unimolecularThreshold = np.zeros((numCoreSpecies),bool) + self.bimolecularThreshold = np.zeros((numCoreSpecies, numCoreSpecies),bool) if self.trimolecular: - self.trimolecularReact = numpy.zeros((numCoreSpecies, numCoreSpecies, numCoreSpecies), bool) - self.trimolecularThreshold = numpy.zeros((numCoreSpecies, numCoreSpecies, numCoreSpecies), bool) + self.trimolecularReact = np.zeros((numCoreSpecies, numCoreSpecies, numCoreSpecies), bool) + self.trimolecularThreshold = np.zeros((numCoreSpecies, numCoreSpecies, numCoreSpecies), bool) else: # By default, react everything - self.unimolecularReact = numpy.ones((numCoreSpecies),bool) - self.bimolecularReact = numpy.ones((numCoreSpecies, numCoreSpecies),bool) + self.unimolecularReact = np.ones((numCoreSpecies),bool) + self.bimolecularReact = np.ones((numCoreSpecies, numCoreSpecies),bool) if self.trimolecular: - self.trimolecularReact = numpy.ones((numCoreSpecies, numCoreSpecies, numCoreSpecies),bool) + self.trimolecularReact = np.ones((numCoreSpecies, numCoreSpecies, numCoreSpecies),bool) # No need to initialize reaction threshold arrays in this case def updateReactionThresholdAndReactFlags(self, @@ -1194,18 +1243,19 @@ def updateReactionThresholdAndReactFlags(self, """ numCoreSpecies = len(self.reactionModel.core.species) prevNumCoreSpecies = len(self.unimolecularReact) - stale = True if numCoreSpecies > prevNumCoreSpecies else False - - - if self.filterReactions: - if stale: - # Reset and expand the react arrays if there were new core species added - self.unimolecularReact = numpy.zeros((numCoreSpecies), bool) - self.bimolecularReact = numpy.zeros((numCoreSpecies, numCoreSpecies), bool) + new_core_species = numCoreSpecies > prevNumCoreSpecies + + # Always reset the react arrays from prior iterations + self.unimolecularReact = np.zeros((numCoreSpecies), bool) + self.bimolecularReact = np.zeros((numCoreSpecies, numCoreSpecies), bool) + if self.trimolecular: + self.trimolecularReact = np.zeros((numCoreSpecies, numCoreSpecies, numCoreSpecies), bool) + if self.filterReactions: + if new_core_species: # Expand the threshold arrays if there were new core species added - unimolecularThreshold = numpy.zeros((numCoreSpecies), bool) - bimolecularThreshold = numpy.zeros((numCoreSpecies, numCoreSpecies), bool) + unimolecularThreshold = np.zeros((numCoreSpecies), bool) + bimolecularThreshold = np.zeros((numCoreSpecies, numCoreSpecies), bool) # Broadcast original thresholds unimolecularThreshold[:prevNumCoreSpecies] = self.unimolecularThreshold @@ -1214,8 +1264,7 @@ def updateReactionThresholdAndReactFlags(self, self.bimolecularThreshold = bimolecularThreshold if self.trimolecular: - self.trimolecularReact = numpy.zeros((numCoreSpecies, numCoreSpecies, numCoreSpecies), bool) - trimolecularThreshold = numpy.zeros((numCoreSpecies, numCoreSpecies, numCoreSpecies), bool) + trimolecularThreshold = np.zeros((numCoreSpecies, numCoreSpecies, numCoreSpecies), bool) trimolecularThreshold[:prevNumCoreSpecies, :prevNumCoreSpecies, :prevNumCoreSpecies] = self.trimolecularThreshold @@ -1248,13 +1297,7 @@ def updateReactionThresholdAndReactFlags(self, self.trimolecularThreshold[i,j,k] = True else: # We are not filtering reactions - if stale: - # Reset and expand the react arrays if there were new core species added - self.unimolecularReact = numpy.zeros((numCoreSpecies), bool) - self.bimolecularReact = numpy.zeros((numCoreSpecies, numCoreSpecies), bool) - if self.trimolecular: - self.trimolecularReact = numpy.zeros((numCoreSpecies, numCoreSpecies, numCoreSpecies), bool) - + if new_core_species: # React all the new core species unimolecularly for i in xrange(prevNumCoreSpecies, numCoreSpecies): self.unimolecularReact[i] = True @@ -1585,7 +1628,7 @@ def loadRMGJavaInput(self, path): assert len(Tlist) > 0 assert len(Plist) > 0 - concentrationList = numpy.array(concentrationList) + concentrationList = np.array(concentrationList) assert concentrationList.shape[1] > 0 # An arbitrary number of concentrations is acceptable, and should be run for each reactor system # Make a reaction system for each (T,P) combination @@ -1593,7 +1636,7 @@ def loadRMGJavaInput(self, path): for P in Plist: for i in range(concentrationList.shape[1]): concentrations = concentrationList[:,i] - totalConc = numpy.sum(concentrations) + totalConc = np.sum(concentrations) initialMoleFractions = dict([(self.initialSpecies[i], concentrations[i] / totalConc) for i in range(len(self.initialSpecies))]) reactionSystem = SimpleReactor(T, P, initialMoleFractions=initialMoleFractions, termination=termination) self.reactionSystems.append(reactionSystem) diff --git a/rmgpy/version.py b/rmgpy/version.py index 55fe2afe7c..ad3d8c090e 100644 --- a/rmgpy/version.py +++ b/rmgpy/version.py @@ -34,4 +34,4 @@ This value can be accessed via `rmgpy.__version__`. """ -__version__ = '2.2.0' +__version__ = '2.2.1'