From accc60902b9e1fd4dac7cdb8c6e7cc89a3c9999e Mon Sep 17 00:00:00 2001 From: Richard West Date: Sat, 18 Oct 2025 23:39:54 -0400 Subject: [PATCH 1/2] Cythonize Molecule.get_element_count() Probably won't speed it up much, but we might as well. --- rmgpy/molecule/molecule.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) 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 From 242798e3f66e1b639f120b3a5192cecd2ce0729e Mon Sep 17 00:00:00 2001 From: Richard West Date: Sat, 18 Oct 2025 23:46:55 -0400 Subject: [PATCH 2/2] Added completedNetworks option to prevent exploring PDep networks that are complete. Sometimes you have a pressure-dependent network that has been studied in great detail and put in a seed mechanism or reaction library, (such as CH2O2) and you don't want RMG to add reactions to the network and your model, because doing so will make it worse. Any reactions net reactions that are in your seed mechanism won't be changed, but currently if there's a new reaction RMG would add it. This option lets you specify in the input file that you don't want RMG to expand certain networks further. You specify them with a line in your pressureDependence block of your input file like completedNetworks = ['CH2O2', 'C2H6'], --- examples/rmg/ethane-oxidation/input.py | 1 + rmgpy/rmg/input.py | 10 +++++++ rmgpy/rmg/model.py | 36 ++++++++++++++++++++++++++ 3 files changed, 47 insertions(+) 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/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: