From 3b83153f841e7c0f6fbdac40835d658f09734569 Mon Sep 17 00:00:00 2001 From: alongd Date: Wed, 9 May 2018 13:41:33 -0400 Subject: [PATCH 01/12] Added methods for collision rate limit checks to Reaction --- rmgpy/reaction.pxd | 8 ++++ rmgpy/reaction.py | 100 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 108 insertions(+) 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 From 16c1b9a5f6e509dc42b9bfc096c05e95b8979d21 Mon Sep 17 00:00:00 2001 From: alongd Date: Wed, 9 May 2018 13:42:12 -0400 Subject: [PATCH 02/12] Check collision limit violators among core reactions --- rmgpy/rmg/main.py | 74 ++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 63 insertions(+), 11 deletions(-) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index 01f1e39787..a6d479d8a7 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -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 @@ -77,7 +78,6 @@ 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 +147,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 +603,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 = [] @@ -952,6 +957,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 From a62c3cec4ca6af10ebe1320a702d751215fe7f70 Mon Sep 17 00:00:00 2001 From: alongd Date: Thu, 28 Jun 2018 10:09:12 -0400 Subject: [PATCH 03/12] Improved module imports in rmg/main.py --- rmgpy/rmg/main.py | 59 +++++++++++++++++++++++------------------------ 1 file changed, 29 insertions(+), 30 deletions(-) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index a6d479d8a7..1feb613159 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 @@ -77,7 +77,6 @@ from rmgpy.stats import ExecutionStatsWriter from rmgpy.thermo.thermoengine import submit from rmgpy.tools.simulate import plot_sensitivity -from cantera import ck2cti ################################################################################ solvent = None @@ -646,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, @@ -657,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...') @@ -814,7 +813,7 @@ def execute(self, **kwargs): else: self.updateReactionThresholdAndReactFlags() - if not numpy.isinf(modelSettings.toleranceThermoKeepSpeciesInEdge): + if not np.isinf(modelSettings.toleranceThermoKeepSpeciesInEdge): self.reactionModel.setThermodynamicFilteringParameters(self.Tmax, toleranceThermoKeepSpeciesInEdge=modelSettings.toleranceThermoKeepSpeciesInEdge, minCoreSizeForPrune=modelSettings.minCoreSizeForPrune, maximumEdgeSpecies=modelSettings.maximumEdgeSpecies, @@ -830,7 +829,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 @@ -1221,19 +1220,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, @@ -1252,12 +1251,12 @@ def updateReactionThresholdAndReactFlags(self, 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) + self.unimolecularReact = np.zeros((numCoreSpecies), bool) + self.bimolecularReact = np.zeros((numCoreSpecies, numCoreSpecies), bool) # 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 @@ -1266,8 +1265,8 @@ def updateReactionThresholdAndReactFlags(self, self.bimolecularThreshold = bimolecularThreshold if self.trimolecular: - self.trimolecularReact = numpy.zeros((numCoreSpecies, numCoreSpecies, numCoreSpecies), bool) - trimolecularThreshold = numpy.zeros((numCoreSpecies, numCoreSpecies, numCoreSpecies), bool) + self.trimolecularReact = np.zeros((numCoreSpecies, numCoreSpecies, numCoreSpecies), bool) + trimolecularThreshold = np.zeros((numCoreSpecies, numCoreSpecies, numCoreSpecies), bool) trimolecularThreshold[:prevNumCoreSpecies, :prevNumCoreSpecies, :prevNumCoreSpecies] = self.trimolecularThreshold @@ -1302,10 +1301,10 @@ def updateReactionThresholdAndReactFlags(self, # 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) + self.unimolecularReact = np.zeros((numCoreSpecies), bool) + self.bimolecularReact = np.zeros((numCoreSpecies, numCoreSpecies), bool) if self.trimolecular: - self.trimolecularReact = numpy.zeros((numCoreSpecies, numCoreSpecies, numCoreSpecies), bool) + self.trimolecularReact = np.zeros((numCoreSpecies, numCoreSpecies, numCoreSpecies), bool) # React all the new core species unimolecularly for i in xrange(prevNumCoreSpecies, numCoreSpecies): @@ -1637,7 +1636,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 @@ -1645,7 +1644,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) From 061734e7736f59a78de2d34b79378141f5f3f7cf Mon Sep 17 00:00:00 2001 From: rgillis8 Date: Fri, 13 Jul 2018 16:37:43 -0400 Subject: [PATCH 04/12] Corrects a geometry reading error for Qchem files --- rmgpy/cantherm/qchem.py | 69 +++++++++++++++++++++++------------------ 1 file changed, 39 insertions(+), 30 deletions(-) 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)) From 90a842991b1f98ec710501150b089d5255dc2174 Mon Sep 17 00:00:00 2001 From: Max Liu Date: Mon, 16 Jul 2018 12:08:26 -0400 Subject: [PATCH 05/12] Change conda build test to use superminimal example --- meta.yaml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/meta.yaml b/meta.yaml index dc72f6120d..142f12d0ce 100644 --- a/meta.yaml +++ b/meta.yaml @@ -70,12 +70,12 @@ requirements: test: source_files: - - 'examples/rmg/minimal' + - 'examples/rmg/superminimal' imports: - rmgpy commands: - - rmg.py examples/rmg/minimal/input.py # [unix] - - python %SCRIPTS%\rmg.py examples\rmg\minimal\input.py # [win] + - rmg.py examples/rmg/superminimal/input.py # [unix] + - python %SCRIPTS%\rmg.py examples\rmg\superminimal\input.py # [win] about: home: http://github.com/ReactionMechanismGenerator/RMG-Py From 5db3e7c5c1fff235cba938e601156fc57a1d1ce8 Mon Sep 17 00:00:00 2001 From: Colin Grambow Date: Mon, 16 Jul 2018 14:32:55 -0400 Subject: [PATCH 06/12] Add missing trimolecularThreshold --- rmgpy/rmg/main.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index 1feb613159..a5c324c8d1 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -791,7 +791,9 @@ def execute(self, **kwargs): except: self.updateReactionThresholdAndReactFlags( rxnSysUnimolecularThreshold = reactionSystem.unimolecularThreshold, - rxnSysBimolecularThreshold = reactionSystem.bimolecularThreshold, skipUpdate=True) + rxnSysBimolecularThreshold = reactionSystem.bimolecularThreshold, + rxnSysTrimolecularThreshold = reactionSystem.trimolecularThreshold, + skipUpdate=True) logging.warn('Reaction thresholds/flags for Reaction System {0} was not updated due to simulation failure'.format(index+1)) else: self.updateReactionThresholdAndReactFlags( From 94a5e9f192cce711b1be724d244353dd0fdcfb47 Mon Sep 17 00:00:00 2001 From: Max Liu Date: Thu, 19 Jul 2018 15:56:08 -0400 Subject: [PATCH 07/12] Adjust updating and resetting of react flags in execute Modify updateReactionThresholdAndReactFlags to always reset the react flags before updating. Also, now the reaction threshold and react flags are updated every iteration, regardless of whether or not core species are added. --- rmgpy/rmg/main.py | 104 +++++++++++++++++++++------------------------- 1 file changed, 47 insertions(+), 57 deletions(-) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index a5c324c8d1..46d5647f74 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -766,55 +766,51 @@ 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, - rxnSysTrimolecularThreshold = reactionSystem.trimolecularThreshold, - 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() - + 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, @@ -1247,15 +1243,16 @@ 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 = np.zeros((numCoreSpecies), bool) - self.bimolecularReact = np.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 = np.zeros((numCoreSpecies), bool) bimolecularThreshold = np.zeros((numCoreSpecies, numCoreSpecies), bool) @@ -1267,7 +1264,6 @@ def updateReactionThresholdAndReactFlags(self, self.bimolecularThreshold = bimolecularThreshold if self.trimolecular: - self.trimolecularReact = np.zeros((numCoreSpecies, numCoreSpecies, numCoreSpecies), bool) trimolecularThreshold = np.zeros((numCoreSpecies, numCoreSpecies, numCoreSpecies), bool) trimolecularThreshold[:prevNumCoreSpecies, :prevNumCoreSpecies, @@ -1301,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 = np.zeros((numCoreSpecies), bool) - self.bimolecularReact = np.zeros((numCoreSpecies, numCoreSpecies), bool) - if self.trimolecular: - self.trimolecularReact = np.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 From 5e144c02e04cf22603b0a197e070102b267acaa2 Mon Sep 17 00:00:00 2001 From: Max Liu Date: Fri, 20 Jul 2018 15:30:59 -0400 Subject: [PATCH 08/12] Clear vertex mapping attribute after running isomorphism Previously, the mapping attribute was causing the symmetry algorithm to return different results. Since nothing uses the mapping attribute aside from isomorphism, we can clear it once we're done with it. Fixes #1344 --- rmgpy/molecule/vf2.pyx | 8 ++++++++ 1 file changed, 8 insertions(+) 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 From 3e3e440315928d049f5e7fa14709d4998b028bc8 Mon Sep 17 00:00:00 2001 From: Max Liu Date: Fri, 20 Jul 2018 15:43:53 -0400 Subject: [PATCH 09/12] Add test for whether vertex mappings are cleared --- rmgpy/molecule/vf2Test.py | 9 +++++++++ 1 file changed, 9 insertions(+) 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__': From f64e5c85d4fa56d1d68bb46bf49beb1b03f327c0 Mon Sep 17 00:00:00 2001 From: Max Liu Date: Mon, 23 Jul 2018 12:55:28 -0400 Subject: [PATCH 10/12] Remove work_in_progress decorator from symmetry test RMG now gives the correct symmetry number after thermo generation --- rmgpy/data/thermoTest.py | 1 - 1 file changed, 1 deletion(-) 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. From 9b7cbfc3fe9dcf1f3b2c9a278da0ed0f2006f595 Mon Sep 17 00:00:00 2001 From: Max Liu Date: Mon, 23 Jul 2018 15:03:31 -0400 Subject: [PATCH 11/12] Update version number to 2.2.1 --- rmgpy/version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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' From c80b0288ec3aa5ca0c10a193154fa7b982ee0ca8 Mon Sep 17 00:00:00 2001 From: Max Liu Date: Mon, 23 Jul 2018 15:03:40 -0400 Subject: [PATCH 12/12] Add release notes for v2.2.1 --- documentation/source/users/rmg/releaseNotes.rst | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) 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