diff --git a/scripts/convertKineticsLibraryToTrainingReactions.ipynb b/scripts/convertKineticsLibraryToTrainingReactions.ipynb index 6b88443bb5..5e1dfffbf7 100644 --- a/scripts/convertKineticsLibraryToTrainingReactions.ipynb +++ b/scripts/convertKineticsLibraryToTrainingReactions.ipynb @@ -19,7 +19,8 @@ }, "outputs": [], "source": [ - "libraryName = 'vinylCPD_H'" + "libraryName = 'vinylCPD_H'\n", + "compareKinetics = True" ] }, { @@ -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" ] }, { @@ -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)" ] }, { @@ -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", @@ -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 = ['']\n", + " html += [''.format(full)]\n", + " html += ['']\n", + " html += [''.format(full, b64encode(fam_rxn._repr_png_()))]\n", + " html += ['']\n", + " html += [''.format(half)]\n", + " html += [''.format(half, ' + '.join([reactant.molecule[0].toSMILES() for reactant in lib_rxn.reactants]))]\n", + " html += ['']\n", + " html += [''.format(half)]\n", + " html += [''.format(half, ' + '.join([product.molecule[0].toSMILES() for product in lib_rxn.products]))]\n", + " html += ['']\n", + " html += [''.format(half)]\n", + " html += [''.format(half, fam_rxn.family)]\n", + " html += ['']\n", + " html += [''.format(half)]\n", + " html += [''.format(half, [entry.label for entry in template])]\n", + " html += ['']\n", " for entry in template:\n", - " group = entry.item\n", - " print entry.label\n", - " display(group)\n", - " print '=================================================='\n", - "\n", + " html += [''.format(full/templateSize, entry.label)]\n", + " html += ['']\n", + " for entry in template:\n", + " html += [''.format(full/templateSize, b64encode(entry.item._repr_png_()))]\n", + " if compareKinetics:\n", + " if not forward:\n", + " html += ['']\n", + " html += ['']\n", + " html += [''.format(half, newKinetics, oldKinetics)]\n", + " html += [''.format(half, figdata)]\n", + " html += ['
One RMG match for this reaction
Reactant SMILES{1}
Product SMILES{1}
Matched Family{1}
Matched Template{1}
{1}
Note: Training reaction written in opposite direction from reaction family.'.format(full)]\n", + " html += ['
New Kinetics:
{1}

Current Kinetics
{2}
']\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 = ['']\n", + " html += ['']\n", + " html += ['']\n", + " html += [''.format(str(lib_rxn))]\n", + " html += ['']\n", + " html += [''.format(b64encode(lib_rxn._repr_png_()))]\n", + " html += ['']\n", + " html += ['']\n", + " html += [''.format(' + '.join([reactant.molecule[0].toSMILES() for reactant in lib_rxn.reactants]))]\n", + " html += ['']\n", + " html += ['']\n", + " html += [''.format(' + '.join([product.molecule[0].toSMILES() for product in lib_rxn.products]))]\n", + " html += ['
Sad :( There are no matches. This is a magic reaction or has chemistry that should be made into a new reaction family
{0}
Reactant SMILES{0}
Product SMILES{0}
']\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 = ['']\n", + " html += [''.format(full)]\n", + " html += ['']\n", + " html += [''.format(full, str(lib_rxn))]\n", + " html += ['']\n", + " html += [''.format(half)]\n", + " html += [''.format(half, ' + '.join([reactant.molecule[0].toSMILES() for reactant in lib_rxn.reactants]))]\n", + " html += ['']\n", + " html += [''.format(half)]\n", + " html += [''.format(half, ' + '.join([product.molecule[0].toSMILES() for product in lib_rxn.products]))]\n", + " html += ['']\n", + " \n", + " html += [''.format(full, i + 1)]\n", + " html += ['']\n", + " html += [''.format(full, b64encode(rxn._repr_png_()))]\n", + " html += ['']\n", + " html += [''.format(half)]\n", + " html += [''.format(half, rxn.family)]\n", + " html += ['']\n", + " html += [''.format(half)]\n", + " html += [''.format(half, [entry.label for entry in template])]\n", + " html += ['']\n", + " for entry in template:\n", + " html += [''.format(full/templateSize, entry.label)]\n", + " html += ['']\n", + " for entry in template:\n", + " html += [''.format(full/templateSize, b64encode(entry.item._repr_png_()))]\n", + " html += ['']\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 += ['']\n", + " html += [''.format(half, figdata)]\n", + " \n", + " html += ['
There are multiple RMG matches for this reaction. You have to manually create this training reaction.
{1}
Reactant SMILES{1}
Product SMILES{1}
Match #{1} - For the following resonance form of the reaction:
Matched Family{1}
Matched Template{1}
{1}
Note: Training reaction written in opposite direction from reaction family.'.format(full)]\n", + " html += ['
'.format(half)]\n", + " html += ['New Kinetics:
{0}'.format(newKinetics)]\n", + " for i, kinetics in enumerate(oldKinetics):\n", + " html += ['

Match #{0} Kinetics:
{1}'.format(i + 1, kinetics)]\n", + " html += ['
']\n", + " \n", + " display(HTML(''.join(html)))\n" ] }, { @@ -382,8 +528,9 @@ } ], "metadata": { + "anaconda-cloud": {}, "kernelspec": { - "display_name": "Python 2", + "display_name": "Python [default]", "language": "python", "name": "python2" }, @@ -397,7 +544,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", - "version": "2.7.11" + "version": "2.7.12" } }, "nbformat": 4,