Skip to content

Commit

Permalink
Merge pull request #157 from ReactionMechanismGenerator/detailed_outp…
Browse files Browse the repository at this point in the history 
…ut_in_convertKinetics_ipynb

Modified convertKineticsLibrarytoTrainingReactions.ipynb to output more details
  • Loading branch information
@nyee
nyee authored Jan 23, 2017
2 parents da2c46c + 7ee1d67 commit 63c3f0e
Showing 1 changed file with 189 additions and 42 deletions.
231 changes: 189 additions & 42 deletions scripts/convertKineticsLibraryToTrainingReactions.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,8 @@
},
"outputs": [],
"source": [
"libraryName = 'vinylCPD_H'"
"libraryName = 'vinylCPD_H'\n",
"compareKinetics = True"
]
},
{
Expand All @@ -34,7 +35,12 @@
"from rmgpy.chemkin import saveChemkinFile, saveSpeciesDictionary\n",
"from rmgpy.rmg.model import Species\n",
"from rmgpy import settings\n",
"from IPython.display import display"
"from IPython.display import display, HTML, Image\n",
"from base64 import b64encode\n",
"if compareKinetics:\n",
" import numpy as np\n",
" import matplotlib.pyplot as plt\n",
" from io import BytesIO"
]
},
{
Expand All @@ -54,7 +60,18 @@
"outputs": [],
"source": [
"database = RMGDatabase()\n",
"database.load(settings['database.directory'], kineticsFamilies='all', reactionLibraries = [libraryName], kineticsDepositories='all')"
"database.load(\n",
" path = settings['database.directory'], \n",
" thermoLibraries = ['primaryThermoLibrary'], # can add others if necessary\n",
" kineticsFamilies = 'all', \n",
" reactionLibraries = [libraryName], \n",
" kineticsDepositories = 'all'\n",
")\n",
"# if we want accurate kinetics comparison, add existing training reactions and fill tree by averaging\n",
"if compareKinetics:\n",
" for family in database.kinetics.families.values():\n",
" family.addKineticsRulesFromTrainingSet(thermoDatabase=database.thermo)\n",
" family.fillKineticsRulesByAveragingUp(verbose=True)"
]
},
{
Expand All @@ -69,12 +86,17 @@
"execution_count": null,
"metadata": {
"collapsed": false,
"scrolled": false
"scrolled": true
},
"outputs": [],
"source": [
"reactionDict = {}\n",
"kineticLibrary = database.kinetics.libraries[libraryName]\n",
"\n",
"# html table settings\n",
"full = 12\n",
"half = full / 2\n",
"\n",
"for index, entry in kineticLibrary.entries.iteritems():\n",
" lib_rxn = entry.item\n",
" lib_rxn.kinetics = entry.data \n",
Expand Down Expand Up @@ -136,46 +158,170 @@
" else:\n",
" reactionDict[fam_rxn.family] = [lib_rxn]\n",
" \n",
" \n",
" template = database.kinetics.families[fam_rxn.family].getReactionTemplate(fam_rxn)\n",
" display(fam_rxn)\n",
" print 'Matched Template: {}'.format([entry.label for entry in template])\n",
" \n",
" templateSize = len(template)\n",
" # html table currently uses a 12 column setup, so templates with 5 groups will break the table\n",
" assert templateSize < 5\n",
" \n",
" if compareKinetics:\n",
" # check what the current kinetics for this template are\n",
" newKinetics = lib_rxn.kinetics\n",
" oldKinetics = database.kinetics.families[fam_rxn.family].getKineticsForTemplate(template, degeneracy=fam_rxn.degeneracy)[0]\n",
" # evaluate kinetics\n",
" tlistinv = np.linspace(1000/1500, 1000/300, num=10)\n",
" tlist = 1000 * np.reciprocal(tlistinv)\n",
" newklist = np.log10(np.array([newKinetics.getRateCoefficient(t) for t in tlist]))\n",
" oldklist = np.log10(np.array([oldKinetics.getRateCoefficient(t) for t in tlist]))\n",
" # create plot\n",
" plt.cla()\n",
" plt.plot(tlistinv, newklist, label='New')\n",
" plt.plot(tlistinv, oldklist, label='Current')\n",
" plt.legend()\n",
" plt.xlabel('1000/T')\n",
" plt.ylabel('log(k)')\n",
" fig = BytesIO()\n",
" plt.savefig(fig, format='png')\n",
" fig.seek(0)\n",
" figdata = b64encode(fig.getvalue())\n",
" fig.close()\n",
" \n",
" # format output using html\n",
" html = ['<table style=\"width:100%;table-layout:fixed;\"><tr>']\n",
" html += ['<th colspan=\"{0}\" style=\"color:green\">One RMG match for this reaction</th>'.format(full)]\n",
" html += ['</tr><tr>']\n",
" html += ['<td colspan=\"{0}\"><img src=\"data:image/png;base64,{1}\"></td>'.format(full, b64encode(fam_rxn._repr_png_()))]\n",
" html += ['</tr><tr>']\n",
" html += ['<th colspan=\"{0}\">Reactant SMILES</th>'.format(half)]\n",
" html += ['<td colspan=\"{0}\">{1}</td>'.format(half, ' + '.join([reactant.molecule[0].toSMILES() for reactant in lib_rxn.reactants]))]\n",
" html += ['</tr><tr>']\n",
" html += ['<th colspan=\"{0}\">Product SMILES</th>'.format(half)]\n",
" html += ['<td colspan=\"{0}\">{1}</td>'.format(half, ' + '.join([product.molecule[0].toSMILES() for product in lib_rxn.products]))]\n",
" html += ['</tr><tr>']\n",
" html += ['<th colspan=\"{0}\">Matched Family</th>'.format(half)]\n",
" html += ['<td colspan=\"{0}\">{1}</td>'.format(half, fam_rxn.family)]\n",
" html += ['</tr><tr>']\n",
" html += ['<th colspan=\"{0}\">Matched Template</th>'.format(half)]\n",
" html += ['<td colspan=\"{0}\">{1}</td>'.format(half, [entry.label for entry in template])]\n",
" html += ['</tr><tr>']\n",
" for entry in template:\n",
" group = entry.item\n",
" print entry.label\n",
" display(group)\n",
" print '=================================================='\n",
"\n",
" html += ['<td colspan=\"{0}\">{1}</td>'.format(full/templateSize, entry.label)]\n",
" html += ['</tr><tr>']\n",
" for entry in template:\n",
" html += ['<td colspan=\"{0}\"><img src=\"data:image/png;base64,{1}\"></td>'.format(full/templateSize, b64encode(entry.item._repr_png_()))]\n",
" if compareKinetics:\n",
" if not forward:\n",
" html += ['</tr><tr>']\n",
" html += ['<th colspan=\"{0}\" style=\"color:red;text-align:center\">Note: Training reaction written in opposite direction from reaction family.'.format(full)]\n",
" html += ['</tr><tr>']\n",
" html += ['<td colspan=\"{0}\"><strong>New Kinetics:</strong><br>{1}<br><br><strong>Current Kinetics</strong><br>{2}</td>'.format(half, newKinetics, oldKinetics)]\n",
" html += ['<td colspan=\"{0}\"><img src=\"data:image/png;base64,{1}\"></td>'.format(half, figdata)]\n",
" html += ['</tr></table>']\n",
" \n",
" display(HTML(''.join(html)))\n",
" elif len(fam_rxn_list) == 0:\n",
" print \"Sad :( There are no matches. This is a magic reaction or has chemistry that should be made into a new reaction family\"\n",
" print ''\n",
" print lib_rxn\n",
" print ''\n",
" print 'Reactant SMILES:'\n",
" for reactant in lib_rxn.reactants:\n",
" print reactant.molecule[0].toSMILES()\n",
" print 'Product SMILES:'\n",
" for product in lib_rxn.products:\n",
" print product.molecule[0].toSMILES()\n",
" print '=================================================='\n",
" html = ['<table style=\"width:100%;table-layout:fixed;\"><tr>']\n",
" html += ['<th colspan=\"2\" style=\"color:red\">Sad :( There are no matches. This is a magic reaction or has chemistry that should be made into a new reaction family</th>']\n",
" html += ['</tr><tr>']\n",
" html += ['<td colspan=\"2\">{0}</td>'.format(str(lib_rxn))]\n",
" html += ['</tr><tr>']\n",
" html += ['<td colspan=\"2\"><img src=\"data:image/png;base64,{0}\"></td>'.format(b64encode(lib_rxn._repr_png_()))]\n",
" html += ['</tr><tr>']\n",
" html += ['<th>Reactant SMILES</th>']\n",
" html += ['<td>{0}</td>'.format(' + '.join([reactant.molecule[0].toSMILES() for reactant in lib_rxn.reactants]))]\n",
" html += ['</tr><tr>']\n",
" html += ['<th>Product SMILES</th>']\n",
" html += ['<td>{0}</td>'.format(' + '.join([product.molecule[0].toSMILES() for product in lib_rxn.products]))]\n",
" html += ['</tr></table>']\n",
" \n",
" display(HTML(''.join(html)))\n",
" else:\n",
" print \"There are multiple RMG matches for this reaction. You have to manually create this training reaction\"\n",
" print ''\n",
" print lib_rxn\n",
" print ''\n",
" print 'Reactant SMILES:'\n",
" for reactant in lib_rxn.reactants:\n",
" print reactant.molecule[0].toSMILES()\n",
" print 'Product SMILES'\n",
" for product in lib_rxn.products:\n",
" print product.molecule[0].toSMILES()\n",
" print ''\n",
" print 'The following families were matched:'\n",
" for rxn in fam_rxn_list:\n",
" print rxn.family\n",
" print '=================================================='\n",
"\n",
"\n"
" if compareKinetics: oldKinetics = []\n",
" lib_reactant = lib_rxn.reactants[0]\n",
" forward = True\n",
" for i, rxn in enumerate(fam_rxn_list):\n",
" # determine direction of family reaction\n",
" fam_reactants = list(rxn.reactants)\n",
" for fam_reactant in fam_reactants:\n",
" if lib_reactant.isIsomorphic(fam_reactant):\n",
" break\n",
" else:\n",
" forward = False\n",
" \n",
" # find the labeled atoms using family and reactants & products from fam_rxn\n",
" family_database = database.kinetics.families[rxn.family]\n",
" family_database.addAtomLabelsForReaction(rxn)\n",
" template = database.kinetics.families[rxn.family].getReactionTemplate(rxn)\n",
" templateSize = len(template)\n",
" \n",
" if compareKinetics:\n",
" oldKinetics.append(database.kinetics.families[rxn.family].getKineticsForTemplate(template, degeneracy=rxn.degeneracy)[0])\n",
" \n",
" if i == 0: \n",
" html = ['<table style=\"width:100%;table-layout:fixed;\"><tr>']\n",
" html += ['<th colspan=\"{0}\" style=\"color:blue\">There are multiple RMG matches for this reaction. You have to manually create this training reaction.</th>'.format(full)]\n",
" html += ['</tr><tr>']\n",
" html += ['<td colspan=\"{0}\">{1}</td>'.format(full, str(lib_rxn))]\n",
" html += ['</tr><tr>']\n",
" html += ['<th colspan=\"{0}\">Reactant SMILES</th>'.format(half)]\n",
" html += ['<td colspan=\"{0}\">{1}</td>'.format(half, ' + '.join([reactant.molecule[0].toSMILES() for reactant in lib_rxn.reactants]))]\n",
" html += ['</tr><tr>']\n",
" html += ['<th colspan=\"{0}\">Product SMILES</th>'.format(half)]\n",
" html += ['<td colspan=\"{0}\">{1}</td>'.format(half, ' + '.join([product.molecule[0].toSMILES() for product in lib_rxn.products]))]\n",
" html += ['</tr><tr>']\n",
" \n",
" html += ['<th colspan=\"{0}\">Match #{1} - For the following resonance form of the reaction:</th>'.format(full, i + 1)]\n",
" html += ['</tr><tr>']\n",
" html += ['<td colspan=\"{0}\"><img src=\"data:image/png;base64,{1}\"></td>'.format(full, b64encode(rxn._repr_png_()))]\n",
" html += ['</tr><tr>']\n",
" html += ['<th colspan=\"{0}\">Matched Family</th>'.format(half)]\n",
" html += ['<td colspan=\"{0}\">{1}</td>'.format(half, rxn.family)]\n",
" html += ['</tr><tr>']\n",
" html += ['<th colspan=\"{0}\">Matched Template</th>'.format(half)]\n",
" html += ['<td colspan=\"{0}\">{1}</td>'.format(half, [entry.label for entry in template])]\n",
" html += ['</tr><tr>']\n",
" for entry in template:\n",
" html += ['<td colspan=\"{0}\">{1}</td>'.format(full/templateSize, entry.label)]\n",
" html += ['</tr><tr>']\n",
" for entry in template:\n",
" html += ['<td colspan=\"{0}\"><img src=\"data:image/png;base64,{1}\"></td>'.format(full/templateSize, b64encode(entry.item._repr_png_()))]\n",
" html += ['</tr><tr>']\n",
" \n",
" if compareKinetics:\n",
" newKinetics = lib_rxn.kinetics\n",
" # evaluate kinetics\n",
" tlistinv = np.linspace(1000/1500, 1000/300, num=10)\n",
" tlist = 1000 * np.reciprocal(tlistinv)\n",
" newklist = np.log10(np.array([newKinetics.getRateCoefficient(t) for t in tlist]))\n",
" oldklist = []\n",
" for kinetics in oldKinetics:\n",
" oldklist.append(np.log10(np.array([kinetics.getRateCoefficient(t) for t in tlist])))\n",
" # create plot\n",
" plt.cla()\n",
" plt.plot(tlistinv, newklist, label='New')\n",
" for i, k in enumerate(oldklist):\n",
" plt.plot(tlistinv, k, label='Match #{0}'.format(i + 1))\n",
" plt.legend()\n",
" plt.xlabel('1000/T')\n",
" plt.ylabel('log(k)')\n",
" fig = BytesIO()\n",
" plt.savefig(fig, format='png')\n",
" fig.seek(0)\n",
" figdata = b64encode(fig.getvalue())\n",
" fig.close()\n",
" \n",
" if not forward:\n",
" html += ['</tr><tr>']\n",
" html += ['<th colspan=\"{0}\" style=\"color:red;text-align:center\">Note: Training reaction written in opposite direction from reaction family.'.format(full)]\n",
" html += ['</tr><tr><td colspan=\"{0}\">'.format(half)]\n",
" html += ['<strong>New Kinetics:</strong><br>{0}'.format(newKinetics)]\n",
" for i, kinetics in enumerate(oldKinetics):\n",
" html += ['<br><br><strong>Match #{0} Kinetics:</strong><br>{1}'.format(i + 1, kinetics)]\n",
" html += ['</td><td colspan=\"{0}\"><img src=\"data:image/png;base64,{1}\"></td>'.format(half, figdata)]\n",
" \n",
" html += ['</tr></table>']\n",
" \n",
" display(HTML(''.join(html)))\n"
]
},
{
Expand Down Expand Up @@ -382,8 +528,9 @@
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 2",
"display_name": "Python [default]",
"language": "python",
"name": "python2"
},
Expand All @@ -397,7 +544,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.11"
"version": "2.7.12"
}
},
"nbformat": 4,
Expand Down

0 comments on commit 63c3f0e

Please sign in to comment.