Skip to content

Commit

Permalink
Merge pull request #112 from ReactionMechanismGenerator/fixTrainingSc…
Browse files Browse the repository at this point in the history 
…ript

Fix Kinetics Library to Training Reaction scripts
  • Loading branch information
@connie
connie authored Jul 6, 2016
2 parents f923ba3 + 118f22a commit 8b616a8
Show file tree
Hide file tree
Showing 2 changed files with 97 additions and 131 deletions.
12 changes: 9 additions & 3 deletions scripts/convertKineticsLibraryToTrainingReactions.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@
"from rmgpy.chemkin import saveChemkinFile, saveSpeciesDictionary\n",
"from rmgpy.rmg.model import Species\n",
"from rmgpy import settings\n",
"from convertKineticsLibraryToTrainingReactions import addAtomLabelsForReaction"
"from IPython.display import display"
]
},
{
Expand Down Expand Up @@ -121,8 +121,9 @@
" break\n",
"\n",
" forward = not (len(fam_reactants) != 0 or len(fam_products) != 0)\n",
" # find the labeled atoms using family and reactants & products from fam_rxn \n",
" addAtomLabelsForReaction(fam_rxn, database)\n",
" # find the labeled atoms using family and reactants & products from fam_rxn\n",
" family_database = database.kinetics.families[fam_rxn.family]\n",
" family_database.addAtomLabelsForReaction(fam_rxn)\n",
" # species replacement so that labeledAtoms is retored\n",
" if forward:\n",
" lib_rxn.reactants = fam_rxn.reactants\n",
Expand All @@ -134,6 +135,11 @@
" reactionDict[fam_rxn.family].append(lib_rxn)\n",
" 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",
" 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",
Expand Down
216 changes: 88 additions & 128 deletions scripts/showNewTrainingRxnsImprovements.ipynb
View file
Original file line number Diff line number Diff line change
@@ -1,18 +1,31 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Show Training Reaction Improvements\n",
"Run this script once before adding training reactions, and once after"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
"collapsed": false
},
"outputs": [],
"source": [
"from rmgpy.data.rmg import RMGDatabase\n",
"from rmgpy.chemkin import saveChemkinFile, saveSpeciesDictionary\n",
"from rmgpy.rmg.model import Species, getFamilyLibraryObject, CoreEdgeReactionModel\n",
"from rmgpy import settings\n",
"from convertKineticsLibraryToTrainingReactions import addAtomLabelsForReaction"
"from IPython.display import display\n",
"from rmgpy.cantherm.output import prettify\n",
"from rmgpy.kinetics.kineticsdata import KineticsData\n",
"from rmgpy.data.kinetics.family import TemplateReaction\n",
"from rmgpy.data.kinetics.depository import DepositoryReaction\n",
"from copy import copy"
]
},
{
Expand All @@ -24,7 +37,7 @@
"outputs": [],
"source": [
"database = RMGDatabase()\n",
"libraries = ['C3']\n",
"libraries = ['vinylCPD_H']\n",
"database.load(settings['database.directory'], kineticsFamilies='all', reactionLibraries = libraries, kineticsDepositories='all')"
]
},
Expand All @@ -48,8 +61,7 @@
" kineticLibrary = database.kinetics.libraries[libraryName]\n",
" for index, entry in kineticLibrary.entries.iteritems():\n",
" lib_rxn = entry.item\n",
" lib_rxn.kinetics = entry.data \n",
" lib_rxn.index = entry.index\n",
" display(lib_rxn)\n",
" # Let's make RMG generate this reaction from the families!\n",
" fam_rxn_list = []\n",
" rxt_mol_mutation_num = 1\n",
Expand Down Expand Up @@ -86,53 +98,14 @@
"\n",
" forward = not (len(fam_reactants) != 0 or len(fam_products) != 0)\n",
" # find the labeled atoms using family and reactants & products from fam_rxn \n",
" addAtomLabelsForReaction(fam_rxn, database)\n",
" family_database = database.kinetics.families[fam_rxn.family]\n",
" family_database.addAtomLabelsForReaction(fam_rxn)\n",
" fam_rxn.index = lib_rxn.index\n",
" reactionDict[fam_rxn.family] = [fam_rxn]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from IPython.display import display\n",
"for fam_rxn in reactionDict['Intra_R_Add_Endocyclic']:\n",
" print fam_rxn.index\n",
" display(fam_rxn)\n",
" for spec in fam_rxn.reactants + fam_rxn.products:\n",
" print spec.molecule[0].toSMILES()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"for index, entry in kineticLibrary.entries.iteritems():\n",
" if entry.index == fam_rxn.index:\n",
" lib_rxn = entry.item\n",
" lib_rxn.kinetics = entry.data \n",
" lib_rxn.index = entry.index\n",
" break\n",
"lib_rxn"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"id(fam_rxn.reactants[0].molecule[0]), id(lib_rxn.reactants[0].molecule[0])"
" reactionDict[lib_rxn]=fam_rxn\n",
" for spec in fam_rxn.reactants + fam_rxn.products:\n",
" print spec.molecule[0].toSMILES()\n",
" else:\n",
" print 'There was an issue finding a single reaction family for this reaction. Please investigate.'"
]
},
{
Expand All @@ -146,14 +119,19 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Before training RMG estimates fam_rxn's kinetics as $ A = 10^9, n = 0.19, E_a = 83.68 kJ/mol $ at [here](http://rmg.mit.edu/database/kinetics/families/Intra_R_Add_Endocyclic/rate_rules/reactant1=multiplicity%25202%250A1%2520%2520C%2520u0%2520p0%2520c0%2520%257B2%252CS%257D%2520%257B3%252CS%257D%2520%257B8%252CS%257D%2520%257B9%252CS%257D%250A2%2520%2520C%2520u0%2520p0%2520c0%2520%257B1%252CS%257D%2520%257B4%252CS%257D%2520%257B5%252CS%257D%2520%257B6%252CS%257D%250A3%2520%2520C%2520u1%2520p0%2520c0%2520%257B1%252CS%257D%2520%257B7%252CD%257D%250A4%2520%2520C%2520u0%2520p0%2520c0%2520%257B2%252CS%257D%2520%257B10%252CT%257D%250A5%2520%2520H%2520u0%2520p0%2520c0%2520%257B2%252CS%257D%250A6%2520%2520H%2520u0%2520p0%2520c0%2520%257B2%252CS%257D%250A7%2520%2520C%2520u0%2520p0%2520c0%2520%257B3%252CD%257D%2520%257B11%252CS%257D%2520%257B12%252CS%257D%250A8%2520%2520H%2520u0%2520p0%2520c0%2520%257B1%252CS%257D%250A9%2520%2520H%2520u0%2520p0%2520c0%2520%257B1%252CS%257D%250A10%2520C%2520u0%2520p0%2520c0%2520%257B4%252CT%257D%2520%257B13%252CS%257D%250A11%2520H%2520u0%2520p0%2520c0%2520%257B7%252CS%257D%250A12%2520H%2520u0%2520p0%2520c0%2520%257B7%252CS%257D%250A13%2520H%2520u0%2520p0%2520c0%2520%257B10%252CS%257D%250A__product1=multiplicity%25202%250A1%2520%2520C%2520u0%2520p0%2520c0%2520%257B2%252CS%257D%2520%257B3%252CS%257D%2520%257B7%252CS%257D%2520%257B8%252CS%257D%250A2%2520%2520C%2520u0%2520p0%2520c0%2520%257B1%252CS%257D%2520%257B4%252CS%257D%2520%257B9%252CS%257D%2520%257B10%252CS%257D%250A3%2520%2520C%2520u0%2520p0%2520c0%2520%257B1%252CS%257D%2520%257B5%252CS%257D%2520%257B6%252CD%257D%250A4%2520%2520C%2520u1%2520p0%2520c0%2520%257B2%252CS%257D%2520%257B5%252CD%257D%250A5%2520%2520C%2520u0%2520p0%2520c0%2520%257B3%252CS%257D%2520%257B4%252CD%257D%2520%257B11%252CS%257D%250A6%2520%2520C%2520u0%2520p0%2520c0%2520%257B3%252CD%257D%2520%257B12%252CS%257D%2520%257B13%252CS%257D%250A7%2520%2520H%2520u0%2520p0%2520c0%2520%257B1%252CS%257D%250A8%2520%2520H%2520u0%2520p0%2520c0%2520%257B1%252CS%257D%250A9%2520%2520H%2520u0%2520p0%2520c0%2520%257B2%252CS%257D%250A10%2520H%2520u0%2520p0%2520c0%2520%257B2%252CS%257D%250A11%2520H%2520u0%2520p0%2520c0%2520%257B5%252CS%257D%250A12%2520H%2520u0%2520p0%2520c0%2520%257B6%252CS%257D%250A13%2520H%2520u0%2520p0%2520c0%2520%257B6%252CS%257D%250A)."
"Check what the RMG estimated kinetics are"
]
},
{
"cell_type": "markdown",
"metadata": {},
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"## step3: after training get fam_rxn's kinetics"
"# Pick a temperature to evaluate the kinetics\n",
"T = 623.0 "
]
},
{
Expand All @@ -177,84 +155,66 @@
},
"outputs": [],
"source": [
"from rmgpy.kinetics.kineticsdata import KineticsData\n",
"from rmgpy.data.kinetics.family import TemplateReaction\n",
"from rmgpy.data.kinetics.depository import DepositoryReaction\n",
"\n",
"for idx, spec in enumerate(fam_rxn.reactants):\n",
" spec = Species(label=spec.label, molecule=spec.molecule)\n",
" spec.generateThermoData(database)\n",
" fam_rxn.reactants[idx] = spec\n",
"for idx, spec in enumerate(fam_rxn.products):\n",
" spec = Species(label=spec.label, molecule=spec.molecule)\n",
" spec.generateThermoData(database)\n",
" fam_rxn.products[idx] = spec\n",
"for libraryName in libraries:\n",
" kineticLibrary = database.kinetics.libraries[libraryName]\n",
" for index, entry in kineticLibrary.entries.iteritems():\n",
" lib_rxn = entry.item\n",
" try:\n",
" fam_rxn = reactionDict[lib_rxn]\n",
" except:\n",
" print 'There was an issue finding a single reaction family for this reaction in step 1. Skipping.'\n",
" continue\n",
" for idx, spec in enumerate(fam_rxn.reactants):\n",
" spec = Species(label=spec.label, molecule=spec.molecule)\n",
" spec.generateThermoData(database)\n",
" fam_rxn.reactants[idx] = spec\n",
" for idx, spec in enumerate(fam_rxn.products):\n",
" spec = Species(label=spec.label, molecule=spec.molecule)\n",
" spec.generateThermoData(database)\n",
" fam_rxn.products[idx] = spec\n",
"\n",
"family = getFamilyLibraryObject(fam_rxn.family)\n",
" family = getFamilyLibraryObject(fam_rxn.family)\n",
"\n",
"# If the reaction already has kinetics (e.g. from a library),\n",
"# assume the kinetics are satisfactory\n",
"if fam_rxn.kinetics is None:\n",
" # Set the reaction kinetics\n",
" kinetics, source, entry, isForward = cem.generateKinetics(fam_rxn)\n",
" fam_rxn.kinetics = kinetics\n",
" # Flip the reaction direction if the kinetics are defined in the reverse direction\n",
" if not isForward:\n",
" fam_rxn.reactants, fam_rxn.products = fam_rxn.products, fam_rxn.reactants\n",
" fam_rxn.pairs = [(p,r) for r,p in fam_rxn.pairs]\n",
" if family.ownReverse and hasattr(fam_rxn,'reverse'):\n",
" # Set the reaction kinetics\n",
" kinetics, source, entryFamily, isForward = cem.generateKinetics(fam_rxn)\n",
" fam_rxn.kinetics = kinetics\n",
" # Flip the reaction direction if the kinetics are defined in the reverse direction\n",
" if not isForward:\n",
" fam_rxn.template = fam_rxn.reverse.template\n",
" # We're done with the \"reverse\" attribute, so delete it to save a bit of memory\n",
" delattr(fam_rxn,'reverse')\n",
" fam_rxn.reactants, fam_rxn.products = fam_rxn.products, fam_rxn.reactants\n",
" fam_rxn.pairs = [(p,r) for r,p in fam_rxn.pairs]\n",
" if family.ownReverse and hasattr(fam_rxn,'reverse'):\n",
" if not isForward:\n",
" fam_rxn.template = fam_rxn.reverse.template\n",
" # We're done with the \"reverse\" attribute, so delete it to save a bit of memory\n",
" delattr(fam_rxn,'reverse')\n",
"\n",
"# convert KineticsData to Arrhenius forms\n",
"if isinstance(fam_rxn.kinetics, KineticsData):\n",
" fam_rxn.kinetics = fam_rxn.kinetics.toArrhenius()\n",
"# correct barrier heights of estimated kinetics\n",
"if isinstance(fam_rxn,TemplateReaction) or isinstance(fam_rxn,DepositoryReaction): # i.e. not LibraryReaction\n",
" fam_rxn.fixBarrierHeight() # also converts ArrheniusEP to Arrhenius.\n",
" # convert KineticsData to Arrhenius forms\n",
" if isinstance(fam_rxn.kinetics, KineticsData):\n",
" fam_rxn.kinetics = fam_rxn.kinetics.toArrhenius()\n",
" # correct barrier heights of estimated kinetics\n",
" if isinstance(fam_rxn,TemplateReaction) or isinstance(fam_rxn,DepositoryReaction): # i.e. not LibraryReaction\n",
" fam_rxn.fixBarrierHeight() # also converts ArrheniusEP to Arrhenius.\n",
"\n",
"if cem.pressureDependence and fam_rxn.isUnimolecular():\n",
" # If this is going to be run through pressure dependence code,\n",
" # we need to make sure the barrier is positive.\n",
" fam_rxn.fixBarrierHeight(forcePositive=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fam_rxn.kinetics"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## step4: compare with lib_rxn's kinetics"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"lib_rxn.kinetics"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Conclusion: it improves the kinetics by factor of 10,000 at 673 for this reaction"
" if cem.pressureDependence and fam_rxn.isUnimolecular():\n",
" # If this is going to be run through pressure dependence code,\n",
" # we need to make sure the barrier is positive.\n",
" fam_rxn.fixBarrierHeight(forcePositive=True)\n",
" print '==============='\n",
" print index\n",
" display(lib_rxn)\n",
" print ''\n",
" print 'Library Kinetics'\n",
" print prettify(repr(entry.data))\n",
" \n",
" print ''\n",
" print 'Reaction Family Kinetics After Training'\n",
" print 'Family: {}'.format(fam_rxn.family)\n",
" print prettify(repr(fam_rxn.kinetics))\n",
" \n",
" print ''\n",
" print 'Kinetics evaluated at {} K'.format(T)\n",
" print 'Library Kinetics: {:.2E}'.format(entry.data.getRateCoefficient(T))\n",
" print 'Reaction Family Kinetics After Training: {:.2E}'.format(fam_rxn.kinetics.getRateCoefficient(T))"
]
}
],
Expand Down

0 comments on commit 8b616a8

Please sign in to comment.