Added completedNetworks option to prevent exploring PDep networks that are complete.

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(

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
 

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

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: