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'],

Author: rwest

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/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: