diff --git a/examples/rmg/ethane-oxidation/input.py b/examples/rmg/ethane-oxidation/input.py index f5c378e069..c90797ec98 100644 --- a/examples/rmg/ethane-oxidation/input.py +++ b/examples/rmg/ethane-oxidation/input.py @@ -61,6 +61,7 @@ temperatures=(300,3000,'K',8), pressures=(0.001,100,'bar',5), interpolation=('Chebyshev', 6, 4), + completedNetworks=['C2H6'], ) options( diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 1a32f07add..26d08b5e84 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -1590,15 +1590,14 @@ def get_element_count(self): """ Returns the element count for the molecule as a dictionary. """ + cython.declare(atom=Atom, element_count=dict, symbol=str, key=str) element_count = {} for atom in self.atoms: symbol = atom.element.symbol - isotope = atom.element.isotope - key = symbol - if key in element_count: - element_count[key] += 1 + if symbol in element_count: + element_count[symbol] += 1 else: - element_count[key] = 1 + element_count[symbol] = 1 return element_count diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index 8c34e815d4..15a158f5f0 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -1326,6 +1326,7 @@ def pressure_dependence( minimumNumberOfGrains=0, interpolation=None, maximumAtoms=None, + completedNetworks=None, ): from arkane.pdep import PressureDependenceJob @@ -1367,6 +1368,10 @@ def pressure_dependence( rmg.pressure_dependence.active_k_rotor = True rmg.pressure_dependence.rmgmode = True + if completedNetworks: + for formula in completedNetworks: + rmg.reaction_model.add_completed_pdep_network(formula) + def options(name='Seed', generateSeedEachIteration=True, saveSeedToDatabase=False, units='si', saveRestartPeriod=None, generateOutputHTML=False, generatePlots=False, generatePESDiagrams=False, saveSimulationProfiles=False, verboseComments=False, @@ -1804,6 +1809,11 @@ def save_input_file(path, rmg): )) f.write(' interpolation = {0},\n'.format(rmg.pressure_dependence.interpolation_model)) f.write(' maximumAtoms = {0}, \n'.format(rmg.pressure_dependence.maximum_atoms)) + if rmg.reaction_model.completed_pdep_networks: + def formula(elements): + return ''.join(f'{el}{count}' if count > 1 else f'{el}' for el, count in elements) + f.write(' completedNetworks = {0},\n'.format( + [formula(net) for net in rmg.reaction_model.completed_pdep_networks])) f.write(')\n\n') # Quantum Mechanics diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index d677755b07..e1afe03624 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -35,6 +35,7 @@ import itertools import logging import os +import re import numpy as np @@ -248,6 +249,27 @@ def __init__(self, core=None, edge=None, surface=None): """ ) ] + self.completed_pdep_networks = set() + + def add_completed_pdep_network(self, formula): + """ + Add a completed pressure-dependent network formula to the set. + """ + # turn C2H4 into {'C':2,'H':4} + if not isinstance(formula, str): + raise TypeError("Expected string for formula, got {0}".format(formula.__class__)) + pattern = r'([A-Z][a-z]?)(\d*)' + element_count = {} + + for match in re.finditer(pattern, formula): + element = match.group(1) + count = int(match.group(2)) if match.group(2) else 1 + element_count[element] = count + # must be hashable and match what is done in add_reaction_to_unimolecular_networks + key = tuple(sorted(element_count.items())) + self.completed_pdep_networks.add(key) + logging.info(f"Added {formula} to list of completed PDep networks that will not be further explored.") + def check_for_existing_species(self, molecule): """ @@ -1893,6 +1915,20 @@ def add_reaction_to_unimolecular_networks(self, newReaction, new_species, networ source = tuple(reactants) + if len(reactants) == 1: + elements = reactants[0].molecule[0].get_element_count() + elif len(products) == 1: + elements = products[0].molecule[0].get_element_count() + else: + raise ValueError("Unimolecular reaction networks can only be formed for unimolecular reactions or isomerizations.") + # make a hashable key from the elements dict + elements_key = tuple(sorted(elements.items())) + if elements_key in self.completed_pdep_networks: + formula = ''.join(f'{el}{count}' if count>1 else el for el, count in elements_key) + logging.info(f"Not adding reaction {newReaction} to unimolecular networks because the network for {formula} is marked as completed.") + return + + # Only search for a network if we don't specify it as a parameter if network is None: if len(reactants) == 1: