diff --git a/arkane/encorr/corr.py b/arkane/encorr/corr.py index 49795c5da58..e3770475dc4 100644 --- a/arkane/encorr/corr.py +++ b/arkane/encorr/corr.py @@ -120,7 +120,7 @@ def get_atom_correction(level_of_theory: Union[LevelOfTheory, CompositeLevelOfTh try: atom_energies = data.atom_energies[energy_level.simple()] except KeyError: - raise AtomEnergyCorrectionError(f'Missing atom energies for {energy_level}') + return corr else: logging.warning(f'No exact atom energy match found for {energy_level}.' f' Using {energy_level.simple()} instead.') diff --git a/arkane/ess/gaussian.py b/arkane/ess/gaussian.py index b421ec02831..121765b6174 100644 --- a/arkane/ess/gaussian.py +++ b/arkane/ess/gaussian.py @@ -257,7 +257,7 @@ def load_conformer(self, symmetry=None, spin_multiplicity=0, optical_isomers=Non return Conformer(E0=(e0 * 0.001, "kJ/mol"), modes=modes, spin_multiplicity=spin_multiplicity, optical_isomers=optical_isomers), unscaled_frequencies - def load_energy(self, zpe_scale_factor=1.): + def load_energy(self, zpe_scale_factor=1. , remove_empirical = False): """ Load the energy in J/mol from a Gaussian log file. The file is checked for a complete basis set extrapolation; if found, that value is @@ -299,6 +299,13 @@ def load_energy(self, zpe_scale_factor=1.): elif 'E(ZPE)' in line: scaled_zpe = float(line.split()[1]) * constants.E_h * constants.Na + elif ' DE(CBS)=' in line: + cbs = float(line.split()[1]) * constants.E_h * constants.Na + elif ' DE(Int)=' in line: + cbs_int = float(line.split()[-1]) * constants.E_h * constants.Na + elif ' DE(Empirical)=' in line: + empirical = float(line.split()[1]) * constants.E_h * constants.Na + elif '\\ZeroPoint=' in line: line = line.strip() + f.readline().strip() start = line.find('\\ZeroPoint=') + 11 @@ -307,16 +314,19 @@ def load_energy(self, zpe_scale_factor=1.): # Read the next line in the file line = f.readline() - if e0_composite is not None: - logging.debug("Using the composite energy from the gaussian output file") - if scaled_zpe is None: - raise LogError('Unable to find zero-point energy in Gaussian log file.') - return e0_composite - scaled_zpe - elif e_elect is not None: - logging.debug("Using the {0} energy from the gaussian output file".format(elect_energy_source)) - return e_elect - else: - raise LogError('Unable to find energy in Gaussian log file.') + if e0_composite is not None and remove_empirical: + logging.debug("removing cbs, interference, and empirical correction for the composite energy") + return e0_composite - scaled_zpe - (cbs + cbs_int + empirical) + elif e0_composite is not None: + logging.debug("Using the composite energy from the gaussian output file") + if scaled_zpe is None: + raise LogError('Unable to find zero-point energy in Gaussian log file.') + return e0_composite - scaled_zpe + elif e_elect is not None: + logging.debug("Using the {0} energy from the gaussian output file".format(elect_energy_source)) + return e_elect + else: + raise LogError('Unable to find energy in Gaussian log file.') def load_zero_point_energy(self): """ diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index 1654331fa25..2ab89b5bac8 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -690,18 +690,20 @@ def read_reaction_comments(reaction, comments, read=True): if reactant.label == reac_str: break else: - raise ChemkinError('Unexpected species identifier {0} encountered in flux pairs ' - 'for reaction {1}.'.format(reac_str, reaction)) + logging.error('Unexpected species identifier {0} encountered in flux pairs ' + 'for reaction {1}.'.format(reac_str, reaction)) + continue if prod_str[-1] == ';': prod_str = prod_str[:-1] for product in reaction.products: if product.label == prod_str: break else: - raise ChemkinError('Unexpected species identifier {0} encountered in flux pairs ' - 'for reaction {1}.'.format(prod_str, reaction)) + logging.error('Unexpected species identifier {0} encountered in flux pairs ' + 'for reaction {1}.'.format(prod_str, reaction)) + continue reaction.pairs.append((reactant, product)) - assert len(reaction.pairs) == max(len(reaction.reactants), len(reaction.products)) + # assert len(reaction.pairs) == max(len(reaction.reactants), len(reaction.products)) elif isinstance(reaction, TemplateReaction) and 'rate rule ' in line: bracketed_rule = tokens[-1] diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 8ea1b553d33..0b905aeeccc 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -430,6 +430,8 @@ def __init__(self, forward_recipe=None, reverse_template=None, reverse_recipe=None, + only_forward=False, + only_reverse=False, forbidden=None, boundary_atoms=None, tree_distances=None, @@ -442,6 +444,8 @@ def __init__(self, self.forward_recipe = forward_recipe self.reverse_template = reverse_template self.reverse_recipe = reverse_recipe + self.only_forward = only_forward + self.only_reverse = only_reverse if self.reversible else False self.forbidden = forbidden self.own_reverse = forward_template is not None and reverse_template is None self.boundary_atoms = boundary_atoms @@ -658,6 +662,8 @@ def load(self, path, local_context=None, global_context=None, depository_labels= local_context['False'] = False local_context['reverse'] = None local_context['reversible'] = None + local_context['only_forward'] = None + local_context['only_reverse'] = None local_context['boundaryAtoms'] = None local_context['treeDistances'] = None local_context['reverseMap'] = None @@ -692,6 +698,10 @@ def load(self, path, local_context=None, global_context=None, depository_labels= self.reverse = local_context.get('reverse', None) self.reversible = True if local_context.get('reversible', None) is None else local_context.get('reversible', None) self.forward_template.products = self.generate_product_template(self.forward_template.reactants) + self.only_forward = False if local_context.get('only_forward', None) is None\ + else local_context.get('only_forward', None) + self.only_reverse = False if local_context.get('only_reverse', None) is None or not self.reversible\ + else local_context.get('only_reverse', None) if self.reversible: self.reverse_template = Reaction(reactants=self.forward_template.products, products=self.forward_template.reactants) @@ -1233,40 +1243,43 @@ def add_rules_from_training(self, thermo_database=None, train_indices=None): template = self.get_reaction_template(item) item.template = self.get_reaction_template_labels(item) - new_degeneracy = self.calculate_degeneracy(item) + try: + new_degeneracy = self.calculate_degeneracy(item) + + if isinstance(entry.data, SurfaceArrhenius): + data = SurfaceArrheniusBEP( + # analogous to Arrhenius.to_arrhenius_ep + A=deepcopy(data.A), + n=deepcopy(data.n), + alpha=0, + E0=deepcopy(data.Ea), + Tmin=deepcopy(data.Tmin), + Tmax=deepcopy(data.Tmax) + ) + else: + data = data.to_arrhenius_ep() - if isinstance(entry.data, SurfaceArrhenius): - data = SurfaceArrheniusBEP( - # analogous to Arrhenius.to_arrhenius_ep - A=deepcopy(data.A), - n=deepcopy(data.n), - alpha=0, - E0=deepcopy(data.Ea), - Tmin=deepcopy(data.Tmin), - Tmax=deepcopy(data.Tmax) + new_entry = Entry( + index=index, + label=';'.join([g.label for g in template]), + item=Reaction(reactants=[g.item for g in template], + products=[]), + data=data, + rank=entry.rank, + reference=entry.reference, + short_desc="Rate rule generated from training reaction {0}. ".format(entry.index) + entry.short_desc, + long_desc="Rate rule generated from training reaction {0}. ".format(entry.index) + entry.long_desc, ) - else: - data = data.to_arrhenius_ep() - - new_entry = Entry( - index=index, - label=';'.join([g.label for g in template]), - item=Reaction(reactants=[g.item for g in template], - products=[]), - data=data, - rank=entry.rank, - reference=entry.reference, - short_desc="Rate rule generated from training reaction {0}. ".format(entry.index) + entry.short_desc, - long_desc="Rate rule generated from training reaction {0}. ".format(entry.index) + entry.long_desc, - ) - new_entry.data.comment = "From training reaction {1} used for {0}".format(';'.join([g.label for g in template]), entry.index) + new_entry.data.comment = "From training reaction {1} used for {0}".format(';'.join([g.label for g in template]), entry.index) - new_entry.data.A.value_si /= new_degeneracy - try: - self.rules.entries[new_entry.label].append(new_entry) - except KeyError: - self.rules.entries[new_entry.label] = [new_entry] - index += 1 + new_entry.data.A.value_si /= new_degeneracy + try: + self.rules.entries[new_entry.label].append(new_entry) + except KeyError: + self.rules.entries[new_entry.label] = [new_entry] + index += 1 + except: + pass def get_root_template(self): """ @@ -1620,17 +1633,23 @@ def is_molecule_forbidden(self, molecule): return False - def _create_reaction(self, reactants, products, is_forward): + def _create_reaction(self, reactants, products, is_forward, consider_only_direction=False): """ Create and return a new :class:`Reaction` object containing the provided `reactants` and `products` as lists of :class:`Molecule` objects. + + If `consider_only_direction` is True, a reaction will only be identified in the direction + the family allows using the only_forward and only_reverse flags. """ # Make sure the products are in fact different than the reactants if same_species_lists(reactants, products): return None + if consider_only_direction and ((self.only_reverse and is_forward) or (self.only_forward and not is_forward)): + return None + # Create and return template reaction object reaction = TemplateReaction( reactants=reactants if is_forward else products, @@ -1707,16 +1726,23 @@ def generate_reactions(self, reactants, products=None, prod_resonance=True): # Forward direction (the direction in which kinetics is defined) reaction_list.extend( - self._generate_reactions(reactants, products=products, forward=True, prod_resonance=prod_resonance)) + self._generate_reactions(reactants, products=products, forward=True, prod_resonance=prod_resonance, + consider_only_direction=True)) if not self.own_reverse and self.reversible: # Reverse direction (the direction in which kinetics is not defined) reaction_list.extend( - self._generate_reactions(reactants, products=products, forward=False, prod_resonance=prod_resonance)) + self._generate_reactions(reactants, products=products, forward=False, prod_resonance=prod_resonance, + consider_only_direction=True)) + + # filter out reactions that have number of reactants or products that do not match the family + if self.reactant_num and self.product_num: + reaction_list = [rxn for rxn in reaction_list if + len(rxn.reactants) == self.reactant_num and len(rxn.products) == self.product_num] return reaction_list - def add_reverse_attribute(self, rxn, react_non_reactive=True): + def add_reverse_attribute(self, rxn, react_non_reactive=True, consider_only_direction=False): """ For rxn (with species' objects) from families with ownReverse, this method adds a `reverse` attribute that contains the reverse reaction information (like degeneracy) @@ -1743,7 +1769,8 @@ def add_reverse_attribute(self, rxn, react_non_reactive=True): reaction_list = self._generate_reactions([spc.molecule for spc in rxn.products], products=rxn.reactants, forward=True, - react_non_reactive=react_non_reactive) + react_non_reactive=react_non_reactive, + consider_only_direction=False) reactions = find_degenerate_reactions(reaction_list, same_reactants, kinetics_family=self) if len(reactions) == 0: logging.error("Expecting one matching reverse reaction, not zero in reaction family {0} for " @@ -1766,7 +1793,8 @@ def add_reverse_attribute(self, rxn, react_non_reactive=True): try: reaction_list = self._generate_reactions([spc.molecule for spc in rxn.products], products=rxn.reactants, forward=True, - react_non_reactive=react_non_reactive) + react_non_reactive=react_non_reactive, + consider_only_direction=False) reactions = find_degenerate_reactions(reaction_list, same_reactants, kinetics_family=self) finally: self.forbidden = temp_object @@ -1859,7 +1887,7 @@ def calculate_degeneracy(self, reaction): reactions = [] for combo in molecule_combos: reactions.extend(self._generate_reactions(combo, products=reaction.products, forward=True, - react_non_reactive=True)) + react_non_reactive=True, consider_only_direction=False)) # remove degenerate reactions reactions = find_degenerate_reactions(reactions, same_reactants, template=reaction.template, @@ -1868,16 +1896,16 @@ def calculate_degeneracy(self, reaction): # log issues if len(reactions) != 1: for reactant in reaction.reactants: - logging.error("Reactant: {0!r}".format(reactant)) + logging.debug("Reactant: {0!r}".format(reactant)) for product in reaction.products: - logging.error("Product: {0!r}".format(product)) + logging.debug("Product: {0!r}".format(product)) raise KineticsError(('Unable to calculate degeneracy for reaction {0} ' 'in reaction family {1}. Expected 1 reaction ' 'but generated {2}').format(reaction, self.label, len(reactions))) return reactions[0].degeneracy def _generate_reactions(self, reactants, products=None, forward=True, prod_resonance=True, - react_non_reactive=False): + react_non_reactive=False, consider_only_direction=False): """ Generate a list of all the possible reactions of this family between the list of `reactants`. The number of reactants provided must match @@ -1954,7 +1982,8 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson pass else: if product_structures is not None: - rxn = self._create_reaction(reactant_structures, product_structures, forward) + rxn = self._create_reaction(reactant_structures, product_structures, forward, + consider_only_direction=consider_only_direction) if rxn: rxn_list.append(rxn) # Bimolecular reactants: A + B --> products @@ -1995,7 +2024,8 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson pass else: if product_structures is not None: - rxn = self._create_reaction(reactant_structures, product_structures, forward) + rxn = self._create_reaction(reactant_structures, product_structures, forward, + consider_only_direction=consider_only_direction) if rxn: rxn_list.append(rxn) # Only check for swapped reactants if they are different @@ -2017,7 +2047,8 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson else: if product_structures is not None: rxn = self._create_reaction(reactant_structures, product_structures, - forward) + forward, + consider_only_direction=consider_only_direction) if rxn: rxn_list.append(rxn) @@ -2073,7 +2104,8 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson pass else: if product_structures is not None: - rxn = self._create_reaction(reactant_structures, product_structures, forward) + rxn = self._create_reaction(reactant_structures, product_structures, forward, + consider_only_direction=consider_only_direction) if rxn: rxn_list.append(rxn) else: @@ -2140,7 +2172,8 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson pass else: if product_structures is not None: - rxn = self._create_reaction(reactant_structures, product_structures, forward) + rxn = self._create_reaction(reactant_structures, product_structures, forward, + consider_only_direction=consider_only_direction) if rxn: rxn_list.append(rxn) @@ -2187,7 +2220,8 @@ def generate_products_and_reactions(order): if _productStructures is not None: _rxn = self._create_reaction(_reactantStructures, _productStructures, - forward) + forward, + consider_only_direction) if _rxn: rxn_list.append(_rxn) diff --git a/rmgpy/data/kinetics/library.py b/rmgpy/data/kinetics/library.py index 74f7f4f4a30..99b6f52bfc1 100644 --- a/rmgpy/data/kinetics/library.py +++ b/rmgpy/data/kinetics/library.py @@ -95,6 +95,20 @@ def __init__(self, self.library = library self.family = library self.entry = entry + if self.kinetics.is_pressure_dependent(): + success = self.generate_high_p_limit_kinetics() + if success: + self.kinetics = self.network_kinetics + logging.info(f'Replaced PDep kinetics for library reaction {self} with the high pressure limit rate.') + else: + logging.error(f'library reaction {self} is pressure dependent, and its high-pressure limit rate ' + f'could not be extracted. Ignoring this library reaction (the rate could still be ' + f'estimated from families).') + self.kinetics = None + self.reactants = list() + self.products = list() + self.entry = None + def __reduce__(self): """ @@ -170,8 +184,8 @@ def generate_high_p_limit_kinetics(self): if isinstance(self.kinetics, (Lindemann, Troe)): self.network_kinetics = self.kinetics.arrheniusHigh self.network_kinetics.comment = self.kinetics.comment - self.network_kinetics.comment = "Kinetics taken from the arrheniusHigh attribute of a" \ - " Troe/Lindemann exprssion. Originally from reaction library {0}".format(self.library) + # self.network_kinetics.comment = "Kinetics taken from the arrheniusHigh attribute of a" \ + # " Troe/Lindemann exprssion. Originally from reaction library {0}".format(self.library) return True if isinstance(self.kinetics, PDepArrhenius): if self.kinetics.pressures.value_si[-1] >= 9000000: # Pa units diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx index 802cc12db9a..883556ca54f 100644 --- a/rmgpy/kinetics/arrhenius.pyx +++ b/rmgpy/kinetics/arrhenius.pyx @@ -398,14 +398,15 @@ cdef class ArrheniusEP(KineticsModel): Return the activation energy in J/mol corresponding to the given enthalpy of reaction `dHrxn` in J/mol. """ - cdef double Ea - Ea = self._alpha.value_si * dHrxn + self._E0.value_si - if self._E0.value_si > 0: - if dHrxn < 0.0 and Ea < 0.0: - Ea = 0.0 - elif dHrxn > 0.0 and Ea < dHrxn: - Ea = dHrxn - return Ea + return self._E0.value_si + # cdef double Ea + # Ea = self._alpha.value_si * dHrxn + self._E0.value_si + # if self._E0.value_si > 0: + # if dHrxn < 0.0 and Ea < 0.0: + # Ea = 0.0 + # elif dHrxn > 0.0 and Ea < dHrxn: + # Ea = dHrxn + # return Ea cpdef Arrhenius to_arrhenius(self, double dHrxn): """ diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index ed4b9690a33..5fd147b46ec 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -1350,7 +1350,7 @@ def update_atomtypes(self, log_species=True, raise_exception=True): atom.atomtype = get_atomtype(atom, atom.edges) except AtomTypeError: if log_species: - logging.error("Could not update atomtypes for this molecule:\n{0}".format(self.to_adjacency_list())) + logging.debug("Could not update atomtypes for this molecule:\n{0}".format(self.to_adjacency_list())) if raise_exception: raise atom.atomtype = ATOMTYPES['R'] diff --git a/rmgpy/molecule/resonance.py b/rmgpy/molecule/resonance.py index d9233084441..bc9b3914d53 100644 --- a/rmgpy/molecule/resonance.py +++ b/rmgpy/molecule/resonance.py @@ -221,8 +221,8 @@ def generate_resonance_structures(mol, clar_structures=True, keep_isomorphic=Fal if features['is_radical'] and not features['is_aryl_radical']: _generate_resonance_structures(mol_list, [generate_kekule_structure], keep_isomorphic=keep_isomorphic, filter_structures=filter_structures) - _generate_resonance_structures(mol_list, [generate_allyl_delocalization_resonance_structures], - keep_isomorphic=keep_isomorphic, filter_structures=filter_structures) + # _generate_resonance_structures(mol_list, [generate_allyl_delocalization_resonance_structures], + # keep_isomorphic=keep_isomorphic, filter_structures=filter_structures) if features['isPolycyclicAromatic'] and clar_structures: _generate_resonance_structures(mol_list, [generate_clar_structures], keep_isomorphic=keep_isomorphic, filter_structures=filter_structures) diff --git a/rmgpy/rmg/decay.py b/rmgpy/rmg/decay.py new file mode 100644 index 00000000000..646e3adf0a4 --- /dev/null +++ b/rmgpy/rmg/decay.py @@ -0,0 +1,111 @@ +#!/usr/bin/env python3 + +############################################################################### +# # +# RMG - Reaction Mechanism Generator # +# # +# Copyright (c) 2002-2019 Prof. William H. Green (whgreen@mit.edu), # +# Prof. Richard H. West (r.west@neu.edu) and the RMG Team (rmg_dev@mit.edu) # +# # +# Permission is hereby granted, free of charge, to any person obtaining a # +# copy of this software and associated documentation files (the 'Software'), # +# to deal in the Software without restriction, including without limitation # +# the rights to use, copy, modify, merge, publish, distribute, sublicense, # +# and/or sell copies of the Software, and to permit persons to whom the # +# Software is furnished to do so, subject to the following conditions: # +# # +# The above copyright notice and this permission notice shall be included in # +# all copies or substantial portions of the Software. # +# # +# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING # +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER # +# DEALINGS IN THE SOFTWARE. # +# # +############################################################################### +""" +Contains classes for automatically breaking down unstable or non-existent species +""" +import logging +from rmgpy.molecule.group import Group +from rmgpy.species import Species +from rmgpy.molecule.molecule import Molecule +from rmgpy.data.kinetics.family import ReactionRecipe + +decay_group_recipe_pairs = [(Group().from_adjacency_list(""" + 1 *3 O u0 p2 c0 {2,S} {4,S} + 2 *2 O u0 p2 c0 {1,S} {3,S} + 3 *1 R!H u1 px c0 {2,S} + 4 H u0 p0 c0 {1,S} + """), + ReactionRecipe(actions=[ + ['BREAK_BOND', '*3', 1, '*2'], + ['CHANGE_BOND', '*2', 1, '*1'], + ['LOSE_RADICAL', '*1','1'], + ['GAIN_RADICAL','*3','1'] + ])), + ] + +def decay_species(spc,check_deterministic=True): + """ + recursively decays a species object as long as there are valid decay recipes + raises a warning if the decays could've been done in a different order + """ + decay = None + for mol in spc.molecule: + d = decay_molecule(mol,check_deterministic=check_deterministic) + if d != mol and decay is None: + decay = d + if not check_deterministic: + break + elif d != mol: + logging.warning("found more than one possible decay for {}, may not be able to deterministically determine decay".format(spc.molecule[0].to_smiles())) + break + + if type(decay) == list: + dlist = [] + for dmol in decay: + dspc = Species(molecule=[dmol]) + dspc.generate_resonance_structures() + decay_spc = decay_species(dspc) + dlist.extend(decay_spc) + return dlist + else: + return [spc] + +def decay_molecule(mol,check_deterministic=True): + """ + breaks down a molecule object down one step using the first valid recipe and atom mappings + raises a warning if there were multiple valid recipes/atom mappings + """ + ind = None + for i,(grp,recipe) in enumerate(decay_group_recipe_pairs): + if mol.is_subgraph_isomorphic(grp): #check is True at most twice + if ind is None: + ind = i + if not check_deterministic: + break + else: + logging.warning("found more than one possible decay for {}, may not be able to deterministically determine decay".format(mol.to_smiles())) + break + + if ind is None: + return mol + else: + grp,recipe = decay_group_recipe_pairs[ind] + mol2 = mol.copy(deep=True) + subs = mol2.find_subgraph_isomorphisms(grp) + + if len(subs)>1: + logging.warning("found more than one possible decay for {}, may not be able to deterministically determine decay".format(mol.to_smiles())) + sub = subs[0] + for key,item in sub.items(): + if item.label: + key.label = item.label + recipe.apply_forward(mol2) + mol2.clear_labeled_atoms() + return mol2.split() + \ No newline at end of file diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index ff0c3de9260..305af105723 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -54,6 +54,7 @@ from rmgpy.rmg.react import react_all from rmgpy.species import Species from rmgpy.thermo.thermoengine import submit +from rmgpy.rmg.decay import decay_species ################################################################################ @@ -257,7 +258,7 @@ def check_for_existing_species(self, molecule): # At this point we can conclude that the species is new return None - def make_new_species(self, object, label='', reactive=True, check_existing=True, generate_thermo=True): + def make_new_species(self, object, label='', reactive=True, check_existing=True, generate_thermo=True, check_decay=False): """ Formally create a new species from the specified `object`, which can be either a :class:`Molecule` object or an :class:`rmgpy.species.Species` @@ -280,21 +281,29 @@ def make_new_species(self, object, label='', reactive=True, check_existing=True, spec = self.check_for_existing_species(molecule) if spec is not None: return spec, False - + # If we're here then we're ready to make the new species + try: + spec = Species(label=label,molecule=[molecule],reactive=reactive,thermo=object.thermo, transport_data=object.transport_data) + except AttributeError: + spec = Species(label=label, molecule=[molecule], reactive=reactive) + + spec.generate_resonance_structures() + + if check_decay: + spcs = decay_species(spec) + if len(spcs) == 1: + spec = spcs[0] + else: + return [self.make_new_species(spc) for spc in spcs] + if reactive: self.species_counter += 1 # count only reactive species - species_index = self.species_counter + spec.index = self.species_counter else: - species_index = -1 - try: - spec = Species(index=species_index, label=label, molecule=[molecule], reactive=reactive, - thermo=object.thermo, transport_data=object.transport_data) - except AttributeError: - spec = Species(index=species_index, label=label, molecule=[molecule], reactive=reactive) + spec.index = -1 spec.creation_iteration = self.iteration_num - spec.generate_resonance_structures() spec.molecular_weight = Quantity(spec.molecule[0].get_molecular_weight() * 1000., "amu") if generate_thermo: @@ -418,8 +427,19 @@ def make_new_reaction(self, forward, check_existing=True, generate_thermo=True): """ # Determine the proper species objects for all reactants and products - reactants = [self.make_new_species(reactant, generate_thermo=generate_thermo)[0] for reactant in forward.reactants] - products = [self.make_new_species(product, generate_thermo=generate_thermo)[0] for product in forward.products] + if forward.family and forward.is_forward: + reactants = [self.make_new_species(reactant, generate_thermo=generate_thermo)[0] for reactant in forward.reactants] + products = [] + for product in forward.products: + spcs = self.make_new_species(product, generate_thermo=generate_thermo,check_decay=True) + if type(spcs) == tuple: + products.append(spcs[0]) + elif type(spcs) == list: + products.extend([spc[0] for spc in spcs]) + else: + reactants = [self.make_new_species(reactant, generate_thermo=generate_thermo)[0] for reactant in forward.reactants] + products = [self.make_new_species(product, generate_thermo=generate_thermo)[0] for product in forward.products] + if forward.specific_collider is not None: forward.specific_collider = self.make_new_species(forward.specific_collider)[0] @@ -872,7 +892,7 @@ def generate_kinetics(self, reaction): G298 = reaction.get_free_energy_of_reaction(298) gibbs_is_positive = G298 > -1e-8 - if family.own_reverse and hasattr(reaction, 'reverse'): + if family.own_reverse and len(reaction.products)==len(reaction.reactants) and hasattr(reaction, 'reverse'): if reaction.reverse: # The kinetics family is its own reverse, so we could estimate kinetics in either direction @@ -1357,6 +1377,10 @@ def add_reaction_to_core(self, rxn): ensure it is supposed to be a core reaction (i.e. all of its reactants AND all of its products are in the list of core species). """ + if any([spc.label == 'API' for spc in rxn.reactants + rxn.products]) \ + and not isinstance(rxn, LibraryReaction): + logging.info(f'\n\n!!!!!\nNot considering reaction\n{rxn}\n\n') + return None if rxn not in self.core.reactions: self.core.reactions.append(rxn) if rxn in self.edge.reactions: diff --git a/rmgpy/thermo/thermoengine.py b/rmgpy/thermo/thermoengine.py index c0ce0d76d23..f317318206e 100644 --- a/rmgpy/thermo/thermoengine.py +++ b/rmgpy/thermo/thermoengine.py @@ -60,15 +60,23 @@ def process_thermo_data(spc, thermo0, thermo_class=NASA, solvent_name=''): logging.debug('Solvent database or solvent_name not found. Solvent effect was not utilized') solvent_data = None else: - solvent_data = solvation_database.get_solvent_data(solvent_name) + try: + solvent_data = solvation_database.get_solvent_data(solvent_name) + except Exception as e: + logging.debug(f'Could not determining solvent thermo corrections for {spc.label}. Got:\n{e}') + solvent_data = None if solvent_data and not "Liquid thermo library" in thermo0.comment: solvation_database = get_db('solvation') - solute_data = solvation_database.get_solute_data(spc) - solvation_correction = solvation_database.get_solvation_correction(solute_data, solvent_data) - # correction is added to the entropy and enthalpy - wilhoit.S0.value_si = (wilhoit.S0.value_si + solvation_correction.entropy) - wilhoit.H0.value_si = (wilhoit.H0.value_si + solvation_correction.enthalpy) - wilhoit.comment += ' + Solvation correction with {} as solvent and solute estimated using {}'.format(solvent_name, solute_data.comment) + try: + solute_data = solvation_database.get_solute_data(spc) + except Exception as e: + logging.debug(f'Could not determining solvent thermo corrections for {spc.label}. Got:\n{e}') + else: + solvation_correction = solvation_database.get_solvation_correction(solute_data, solvent_data) + # correction is added to the entropy and enthalpy + wilhoit.S0.value_si = (wilhoit.S0.value_si + solvation_correction.entropy) + wilhoit.H0.value_si = (wilhoit.H0.value_si + solvation_correction.enthalpy) + wilhoit.comment += ' + Solvation correction with {} as solvent and solute estimated using {}'.format(solvent_name, solute_data.comment) # Compute E0 by extrapolation to 0 K if spc.conformer is None: