Add a new script to plot rxn path diagram with nodes labeled with species images

model_analyzer/plotRxnPathDiagram/Readme.txt

@@ -0,0 +1,23 @@
+The script 'plotRxnPath.py' in this folder is used to plot the reaction path diagram in a constant volume, adiabatic kinetics simulation.
+
+The simulation is performed using Cantera.
+
+'chem_annotated.cti' is the mech file.
+
+The dir 'structure' contains the species images which could be obtained from an RMG job.
+
+The script 'resizePNGSize.py' could be used to resize the structure images if you find they are not in proper size in the rxn path diagram.
+
+
+Previously the pathway analysis tool in cantera could generate a rxn pathway diagram with the species labeled with species names. 
+It's a boring work to indentify the major species in this diagram.
+Here I am trying to replace the text labels with the species images.
+Hopefully it will make the model analysis process less painful.
+
+It could be modified to simulate a constant pressure, constant temperature reactor by setting the wall expansion_rate_coeff and heat transfer coeff large enough. 
+If you did that, please push it to the RMG/WIP_script repo.
+
+Peng
+Nov-14-2016
+
+

model_analyzer/plotRxnPathDiagram/plotRxnPath.py

@@ -0,0 +1,116 @@
+"""
+Constant volume, adiabatic kinetics simulation. Plot rxn path diagram.
+
+Previously the pathway analysis tool in cantera could generate a rxn pathway diagram with the species labeled with
+species names. It's a boring work to indentify the major species in this diagram.
+Here I am trying to replace the text label with the species IMG.
+Hopefully it will make the model analysis process less painful.
+
+It could be modified to simulate a constant pressure, constant temperature reactor by setting the wall expansion_rate_coeff
+and heat transfer coeff large enough. If you did that, please push it to the RMG/WIP_script repo.
+
+Peng Zhang
+Nov-14-2016
+"""
+
+import sys
+import os
+import csv
+import numpy as np
+
+import cantera as ct
+
+###########################
+# Input section
+gas = ct.Solution('chem_annotated.cti') # Input mech file.
+air = ct.Solution('air.xml')
+
+P = 20      # Initial pressure. Unit: bar 
+T = 700     # Initial temperature. Unit: K
+dt = 1e-5   # Time step.
+n_steps = int(1.0/dt) # Total time steps. Set the end time long enough. Here we set it to be 1 second. 
+
+gas.TPX = T, P*1e5, 'butane(5):10, O2(2):68.5, N2:274'  # Mixture composition. C8H10O(1):1, butane(5):9, 
+###########################
+# Cantera reaction network setting
+r = ct.IdealGasReactor(gas)
+env = ct.Reservoir(air)
+w = ct.Wall(r, env)
+w.expansion_rate_coeff = 0.0
+w.area = 1.0
+sim = ct.ReactorNet([r])
+###########################
+time = 0.0 # Initialize time
+for n in range(n_steps):
+    time += dt
+    sim.advance(time)
+    # print r.T    # Uncomment this line to moniter the calculation progress.
+
+    if r.T > T + 50:
+        print 'Initial temperature: {0} K'.format(T)
+        print 'End temperature:     {0:.0f} K'.format(r.T)
+        print 'End time:            {0:.2f} ms'.format(time*1000)
+        ###########################################
+        # Reaction path diagram setting
+        element = 'C'
+        diagram = ct.ReactionPathDiagram(gas, element)
+        diagram.threshold = 0.05    # Default value is 0.005
+        diagram.scale = bool(0)
+        diagram.label_threshold = 1/100000000 # For now, I'm not sure what this is used for.
+        diagram.show_details = bool(0) # Make it 1 to turn on details.
+        ###########################################
+
+        dot_file = 'rxnpath.dot'
+        diagram.write_dot(dot_file)
+
+        #---------------------------------------------------------------------------------------
+        # Plot the rxn path diagram and label nodes with spc structures.
+        dotFileOriginal = dot_file
+        dotFileNew = 'rxnpath_Structure.dot'    # This is an intermediate file, which will be deleted in the end.
+        imgFileNew = 'rxnpath_Structure.png'    # This is the rxn path pic. 
+        structureDir = 'structure'              # This dir stores the structure of the species. 
+
+        structureDict = {}
+        for filename in os.listdir(structureDir):
+            if filename.endswith(').png'): 
+                fnLeftParentheses = filename.rfind('(')
+                fnRightParenthese = filename.rfind(')')
+                structureDict[filename[fnLeftParentheses:fnRightParenthese+1]] = filename
+
+        with open(dotFileOriginal, 'r') as f:
+            stream = f.readlines()
+        for i in range(len(stream)):
+            indexLabel = stream[i].find('label')
+            if indexLabel > 0:
+                indexLeftQuote = stream[i].find('="', indexLabel)
+                indexRightQuote = stream[i].find('"];', indexLabel)
+                stripLabelQuote = stream[i][indexLeftQuote+2:indexRightQuote]
+                if 'fwd:' not in stripLabelQuote and '(' in stripLabelQuote:
+                    indexLeftParentheses = stripLabelQuote.find('(')
+                    indexRightParentheses = stripLabelQuote.find(')')
+                    spcIndex = stripLabelQuote[indexLeftParentheses: indexRightParentheses+1]
+                    if spcIndex in structureDict:
+                        newLabel = """penwidth=0, shape=box, margin=0, label=<<TABLE border="0">
+        <TR><TD><IMG scale="True" SRC="{0}"/></TD></TR>
+        <TR><TD>{1}</TD></TR>
+        </TABLE>>];
+""".format(os.path.join(structureDir, structureDict[spcIndex]), stripLabelQuote)
+                        stream[i] = stream[i][:indexLabel]+newLabel
+                    else:
+                        print 'Warning: {0} is not in the structure dictionary! Use the original label.'.format(stripLabelQuote)
+                        stream[i] = stream[i][:indexLabel] + 'penwidth=0, shape=box, margin=0, ' + stream[i][indexLabel:]
+
+        with open(dotFileNew, 'w') as f:
+            f.writelines(stream)
+
+        os.system('dot {0} -Tpng -o{1} -Gdpi=300'.format(dotFileNew, imgFileNew))
+        print("Wrote output file to '{0}'.".format(imgFileNew))
+
+        # Remove the intermediate dot file.
+        os.remove(dotFileOriginal)
+        os.remove(dotFileNew)
+
+        break
+# If the for loop didn't break, that means the temperature rise is less than 50 K. 
+else:
+    print 'Temperature rise is less than 50 K. Please reset the initial condition.'

model_analyzer/plotRxnPathDiagram/resizePNGSize.py

@@ -0,0 +1,20 @@
+"""
+Resize the structure png file to make it visible in the reaction pathway diagram.
+Please adjust the factor `f` accordingly.
+Peng Zhang
+Nov-14-2016
+"""
+
+import os
+from PIL import Image
+
+structureDir = 'structure'
+for filename in os.listdir(structureDir):
+    if filename.endswith(').png'):
+        filename = os.path.join(structureDir, filename)
+        im = Image.open(filename)
+        # f = 200.0/max(im.width, im.height)
+        f = 3
+        im = im.resize((int(im.width*f), int(im.height*f)))
+        im.save(filename)
+