More helpful debugging for PDep networks, and option to always generatePESDiagrams.

arkane/pdep.py

@@ -618,14 +618,18 @@ def plot(self, output_directory):
                 plt.savefig(os.path.join(output_directory, 'plots', 'kinetics_{0:d}.pdf'.format(count)))
                 plt.close()
 
-    def draw(self, output_directory, file_format='pdf'):
+    def draw(self, output_directory, file_format='pdf', filename_stem='network'):
         """
         Generate a PDF drawing of the pressure-dependent reaction network.
         This requires that Cairo and its Python wrapper be available; if not,
         the drawing is not generated.
 
         You may also generate different formats of drawings, by changing format to
         one of the following: `pdf`, `svg`, `png`.
+
+        The default filename stem is 'network', which will result in a file
+        named 'network.pdf' in the specified output directory. You can change
+        this by passing a different `filename_stem` argument.
         """
 
         # Skip this step if cairo is not installed
@@ -639,7 +643,7 @@ def draw(self, output_directory, file_format='pdf'):
 
         from rmgpy.pdep.draw import NetworkDrawer
 
-        path = os.path.join(output_directory, 'network.' + file_format)
+        path = os.path.join(output_directory, f'{filename_stem}.{file_format}')
 
         NetworkDrawer().draw(self.network, file_format=file_format, path=path)
 

documentation/source/users/rmg/input.rst

@@ -965,6 +965,7 @@ Miscellaneous options::
         units='si',
         generateOutputHTML=True,
         generatePlots=False,
+		generatePESDiagrams=False,
         saveSimulationProfiles=True,
         verboseComments=False,
         saveEdgeSpecies=True,
@@ -986,6 +987,8 @@ HTML file for your model containing all the species and reactions.  Turning this
 
 Setting ``generatePlots`` to ``True`` will generate a number of plots describing the statistics of the RMG job, including the reaction model core and edge size and memory use versus  execution time. These will be placed in the output directory in the plot/ folder.
 
+Setting ``generatePESDiagrams`` to ``True`` will generate potential energy surface diagrams for each pressure dependent network in the model.  These diagrams will be saved in the ``pdep/`` folder in the output directory. Only applicable if pressure dependence is enabled.
+
 Setting ``saveSimulationProfiles`` to ``True`` will make RMG save csv files of the simulation in .csv files in the ``solver/`` folder.  The filename will be ``simulation_1_26.csv`` where the first number corresponds to the reaciton system, and the second number corresponds to the total number of species at the point of the simulation.  Therefore, the highest second number will indicate the latest simulation that RMG has complete while enlarging the core model.  The information inside the csv file will provide the time, reactor volume in m^3, as well as mole fractions of the individual species.
 
 Setting ``verboseComments`` to ``True`` will make RMG generate chemkin files with complete verbose commentary for the kinetic and thermo parameters.  This will be helpful in debugging what values are being averaged for the kinetics.  Note that this may produce very large files.

examples/rmg/1,3-hexadiene/input.py

@@ -89,4 +89,5 @@
     units='si',
     generateOutputHTML=False,
     generatePlots=False,
+    generatePESDiagrams=True,
 )

examples/rmg/commented/input.py

@@ -217,6 +217,8 @@
     generateOutputHTML=True,
     # generates plots of the RMG's performance statistics. Not helpful if you just want a model.
     generatePlots=False,
+    # generates potential energy surface diagrams for pressure dependent networks in the model.
+    generatePESDiagrams=False,
     # saves mole fraction of species in 'solver/' to help you create plots
     saveSimulationProfiles=False,
     # gets RMG to output comments on where kinetics were obtained in the chemkin file.

examples/rmg/ethane-oxidation/input.py

@@ -67,4 +67,5 @@
     units='si',
     generateOutputHTML=False,
     generatePlots=False,
+    generatePESDiagrams=True,
 )

rmgpy/data/kinetics/family.py

@@ -4487,6 +4487,7 @@ def extract_source_from_comments(self, reaction):
         template_matches = re.search(template_pattern, full_comment_string)
         if autogen_node_matches is not None:  # autogenerated trees
             template_str = autogen_node_matches.group(1).split('Multiplied by reaction path degeneracy')[0].strip()
+            template_str = template_str.split('in family')[0].strip()
             tokens = template_str.split()
             if len(tokens) == 2:  # The node was probably split because wordwrap was turned off
                 assert len(template_str) > 115, 'The node name is too short to have been broken up by the chemkin writer'

rmgpy/data/kinetics/rules.py

@@ -367,20 +367,18 @@ def estimate_kinetics(self, template, degeneracy=1):
 
             kinetics = deepcopy(entry.data)
             if entry0 == entry:
-                kinetics.comment = "Estimated from node {}".format(entry.label)
+                kinetics.comment = f"Estimated from node {entry.label} in family {self.label.replace('/rules','')}."
                 kinetics.A.value_si *= degeneracy
                 if degeneracy > 1:
-                    kinetics.comment += "\n"
-                    kinetics.comment += "Multiplied by reaction path degeneracy {0}".format(degeneracy)
-                return kinetics,entry
+                    kinetics.comment += f"\nMultiplied by reaction path degeneracy {degeneracy}"
+                return kinetics, entry
             else:
-                kinetics.comment = "Matched node {}\n".format(entry0.label)
-                kinetics.comment += "Estimated from node {}".format(entry.label)
+                kinetics.comment = f"Matched node {entry0.label}\n"
+                kinetics.comment += f"Estimated from node {entry.label} in family {self.label.replace('/rules','')}."
                 kinetics.A.value_si *= degeneracy
                 if degeneracy > 1:
-                    kinetics.comment += "\n"
-                    kinetics.comment += "Multiplied by reaction path degeneracy {0}".format(degeneracy)
-                return kinetics,None
+                    kinetics.comment += f"\nMultiplied by reaction path degeneracy {degeneracy}"
+                return kinetics, None
 
         original_leaves = get_template_label(template)
         template_list = [template]

rmgpy/pdep/network.py

@@ -408,6 +408,7 @@ def set_conditions(self, T, P, ymB=None):
                             logging.error("Increasing number of grains did not decrease error enough "
                                           "(Current badness: {0:.1f}, previous {1:.1f}). Something must be wrong with "
                                           "network {2}".format(badness, previous_error.badness(), self.label))
+                            self.log_summary()
                             raise error
                     previous_error = error
                     success = False

rmgpy/rmg/input.py

@@ -1369,7 +1369,7 @@ def pressure_dependence(
 
 
 def options(name='Seed', generateSeedEachIteration=True, saveSeedToDatabase=False, units='si', saveRestartPeriod=None,
-            generateOutputHTML=False, generatePlots=False, saveSimulationProfiles=False, verboseComments=False,
+            generateOutputHTML=False, generatePlots=False, generatePESDiagrams=False, saveSimulationProfiles=False, verboseComments=False,
             saveEdgeSpecies=False, keepIrreversible=False, trimolecularProductReversible=True, wallTime='00:00:00:00',
             saveSeedModulus=-1):
     if saveRestartPeriod:
@@ -1386,6 +1386,9 @@ def options(name='Seed', generateSeedEachIteration=True, saveSeedToDatabase=Fals
         logging.warning('Generate Output HTML option was turned on. Note that this will slow down model generation.')
     rmg.generate_output_html = generateOutputHTML
     rmg.generate_plots = generatePlots
+    rmg.generate_PES_diagrams = generatePESDiagrams
+    if generatePESDiagrams:
+        logging.info('Potential Energy Surface diagrams will be generated in the "pdep" folder.')
     rmg.save_simulation_profiles = saveSimulationProfiles
     rmg.verbose_comments = verboseComments
     if saveEdgeSpecies:
@@ -1835,6 +1838,7 @@ def save_input_file(path, rmg):
     f.write('    units = "{0}",\n'.format(rmg.units))
     f.write('    generateOutputHTML = {0},\n'.format(rmg.generate_output_html))
     f.write('    generatePlots = {0},\n'.format(rmg.generate_plots))
+    f.write('    generatePESDiagrams = {0},\n'.format(rmg.generate_PES_diagrams))
     f.write('    saveSimulationProfiles = {0},\n'.format(rmg.save_simulation_profiles))
     f.write('    saveEdgeSpecies = {0},\n'.format(rmg.save_edge_species))
     f.write('    keepIrreversible = {0},\n'.format(rmg.keep_irreversible))

rmgpy/rmg/main.py

@@ -143,6 +143,7 @@ class RMG(util.Subject):
     `units`                                                    The unit system to use to save output files (currently must be 'si')
     `generate_output_html`                                     ``True`` to draw pictures of the species and reactions, saving a visualized model in an output HTML file.  ``False`` otherwise
     `generate_plots`                                           ``True`` to generate plots of the job execution statistics after each iteration, ``False`` otherwise
+    `generate_PES_diagrams`                                    ``True`` to generate potential energy surface diagrams for pressure dependent networks in the model, ``False`` otherwise
     `verbose_comments`                                         ``True`` to keep the verbose comments for database estimates, ``False`` otherwise
     `save_edge_species`                                        ``True`` to save chemkin and HTML files of the edge species, ``False`` otherwise
     `keep_irreversible`                                        ``True`` to keep ireversibility of library reactions as is ('<=>' or '=>'). ``False`` (default) to force all library reactions to be reversible ('<=>')
@@ -222,6 +223,7 @@ def clear(self):
         self.units = "si"
         self.generate_output_html = None
         self.generate_plots = None
+        self.generate_PES_diagrams = None
         self.save_simulation_profiles = None
         self.verbose_comments = None
         self.save_edge_species = None
@@ -271,6 +273,7 @@ def load_input(self, path=None):
         if self.pressure_dependence:
             self.pressure_dependence.output_file = self.output_directory
             self.reaction_model.pressure_dependence = self.pressure_dependence
+            self.pressure_dependence.generate_PES_diagrams = self.generate_PES_diagrams
         if self.solvent:
             self.reaction_model.solvent_name = self.solvent
 

rmgpy/rmg/pdep.py

@@ -44,7 +44,7 @@
 import rmgpy.reaction
 from rmgpy.constants import R
 from rmgpy.data.kinetics.library import LibraryReaction
-from rmgpy.exceptions import PressureDependenceError, NetworkError
+from rmgpy.exceptions import PressureDependenceError, NetworkError, InvalidMicrocanonicalRateError
 from rmgpy.pdep import Configuration
 from rmgpy.rmg.react import react_species
 from rmgpy.statmech import Conformer
@@ -872,11 +872,20 @@ def update(self, reaction_model, pdep_settings, requires_rms=False):
         if output_directory:
             job.save_input_file(
                 os.path.join(output_directory, 'pdep', 'network{0:d}_{1:d}.py'.format(self.index, len(self.isomers))))
+            if getattr(pdep_settings, 'generate_PES_diagrams', False):
+                job.draw(os.path.join(output_directory, 'pdep'), filename_stem=f'network{self.index:d}_{len(self.isomers):d}', file_format='pdf')
 
         # Calculate the rate coefficients
         self.initialize(Tmin, Tmax, Pmin, Pmax, maximum_grain_size, minimum_grain_count, active_j_rotor, active_k_rotor,
                         rmgmode)
-        K = self.calculate_rate_coefficients(Tlist, Plist, method)
+        try:
+            K = self.calculate_rate_coefficients(Tlist, Plist, method)
+        except InvalidMicrocanonicalRateError:
+            if output_directory:
+                filename_stem = f'network{self.index:d}_{len(self.isomers):d}'
+                job.draw(output_directory, filename_stem=filename_stem, file_format='pdf')
+                logging.info(f"Network {self.index} has been drawn and saved as {filename_stem}.pdf in {output_directory} to aid debugging.")
+            raise
 
         # Generate PDepReaction objects
         configurations = []

test/rmgpy/test_data/parsing_data/chem_annotated.inp

@@ -259,7 +259,7 @@ HCCO(23)(+M)=O(2)+C2H(21)(+M)                       1.000e+00 0.000     0.000
 ! Reaction index: Chemkin #25; RMG #97
 ! Template reaction: Disproportionation
 ! Flux pairs: CH3CHCH3, C3H6(28); C2H5(27), C2H6; 
-! Estimated from node Root_Ext-2R!H-R_2R!H->C_4R->C
+! Estimated from node Root_Ext-2R!H-R_2R!H->C_4R->C in family Disproportionation.
 ! Multiplied by reaction path degeneracy 6.0
 C2H5(27)+CH3CHCH3<=>C2H6+C3H6(28)                   3.000e+11 0.000     0.000