Skip to content

Commit 4cf074d

Browse files
KEHANGconnie
KEHANG
authored and
connie
committed
created an ipynb showing the effectiveness of kinetics improvement
1 parent 70daf53 commit 4cf074d

File tree

1 file changed

+282
-0
lines changed

1 file changed

+282
-0
lines changed

‎scripts/showNewTrainingRxnsImprovements.ipynb‎

Lines changed: 282 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,282 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": null,
6+
"metadata": {
7+
"collapsed": true
8+
},
9+
"outputs": [],
10+
"source": [
11+
"from rmgpy.data.rmg import RMGDatabase\n",
12+
"from rmgpy.chemkin import saveChemkinFile, saveSpeciesDictionary\n",
13+
"from rmgpy.rmg.model import Species, getFamilyLibraryObject, CoreEdgeReactionModel\n",
14+
"from rmgpy import settings\n",
15+
"from convertKineticsLibraryToTrainingReactions import addAtomLabelsForReaction"
16+
]
17+
},
18+
{
19+
"cell_type": "code",
20+
"execution_count": null,
21+
"metadata": {
22+
"collapsed": false
23+
},
24+
"outputs": [],
25+
"source": [
26+
"database = RMGDatabase()\n",
27+
"libraries = ['C3']\n",
28+
"database.load(settings['database.directory'], kineticsFamilies='all', reactionLibraries = libraries, kineticsDepositories='all')"
29+
]
30+
},
31+
{
32+
"cell_type": "markdown",
33+
"metadata": {},
34+
"source": [
35+
"## step1: find fam_rxn for each lib_rxn"
36+
]
37+
},
38+
{
39+
"cell_type": "code",
40+
"execution_count": null,
41+
"metadata": {
42+
"collapsed": false
43+
},
44+
"outputs": [],
45+
"source": [
46+
"reactionDict = {}\n",
47+
"for libraryName in libraries:\n",
48+
" kineticLibrary = database.kinetics.libraries[libraryName]\n",
49+
" for index, entry in kineticLibrary.entries.iteritems():\n",
50+
" lib_rxn = entry.item\n",
51+
" lib_rxn.kinetics = entry.data \n",
52+
" lib_rxn.index = entry.index\n",
53+
" # Let's make RMG generate this reaction from the families!\n",
54+
" fam_rxn_list = []\n",
55+
" rxt_mol_mutation_num = 1\n",
56+
" pdt_mol_mutation_num = 1\n",
57+
" for reactant in lib_rxn.reactants:\n",
58+
" rxt_mol_mutation_num *= len(reactant.molecule)\n",
59+
"\n",
60+
" for product in lib_rxn.products:\n",
61+
" pdt_mol_mutation_num *= len(product.molecule)\n",
62+
"\n",
63+
" for mutation_i in range(rxt_mol_mutation_num):\n",
64+
" rxts_mol = [spc.molecule[mutation_i%(len(spc.molecule))] for spc in lib_rxn.reactants]\n",
65+
" pdts_mol = [spc.molecule[0] for spc in lib_rxn.products]\n",
66+
" fam_rxn_list.extend(database.kinetics.generateReactionsFromFamilies(\n",
67+
" reactants=rxts_mol, products=pdts_mol))\n",
68+
"\n",
69+
" if len(fam_rxn_list) == 1:\n",
70+
" fam_rxn = fam_rxn_list[0] \n",
71+
" lib_reactants = [r for r in lib_rxn.reactants] \n",
72+
" fam_reactants = [r for r in fam_rxn.reactants]\n",
73+
" for lib_reactant in lib_reactants:\n",
74+
" for fam_reactant in fam_reactants:\n",
75+
" if lib_reactant.isIsomorphic(fam_reactant):\n",
76+
" fam_reactants.remove(fam_reactant)\n",
77+
" break\n",
78+
"\n",
79+
" lib_products = [r for r in lib_rxn.products] \n",
80+
" fam_products = [r for r in fam_rxn.products]\n",
81+
" for lib_product in lib_products:\n",
82+
" for fam_product in fam_products:\n",
83+
" if lib_product.isIsomorphic(fam_product):\n",
84+
" fam_products.remove(fam_product)\n",
85+
" break\n",
86+
"\n",
87+
" forward = not (len(fam_reactants) != 0 or len(fam_products) != 0)\n",
88+
" # find the labeled atoms using family and reactants & products from fam_rxn \n",
89+
" addAtomLabelsForReaction(fam_rxn, database)\n",
90+
" fam_rxn.index = lib_rxn.index\n",
91+
" reactionDict[fam_rxn.family] = [fam_rxn]"
92+
]
93+
},
94+
{
95+
"cell_type": "code",
96+
"execution_count": null,
97+
"metadata": {
98+
"collapsed": false
99+
},
100+
"outputs": [],
101+
"source": [
102+
"from IPython.display import display\n",
103+
"for fam_rxn in reactionDict['Intra_R_Add_Endocyclic']:\n",
104+
" print fam_rxn.index\n",
105+
" display(fam_rxn)\n",
106+
" for spec in fam_rxn.reactants + fam_rxn.products:\n",
107+
" print spec.molecule[0].toSMILES()"
108+
]
109+
},
110+
{
111+
"cell_type": "code",
112+
"execution_count": null,
113+
"metadata": {
114+
"collapsed": false
115+
},
116+
"outputs": [],
117+
"source": [
118+
"for index, entry in kineticLibrary.entries.iteritems():\n",
119+
" if entry.index == fam_rxn.index:\n",
120+
" lib_rxn = entry.item\n",
121+
" lib_rxn.kinetics = entry.data \n",
122+
" lib_rxn.index = entry.index\n",
123+
" break\n",
124+
"lib_rxn"
125+
]
126+
},
127+
{
128+
"cell_type": "code",
129+
"execution_count": null,
130+
"metadata": {
131+
"collapsed": false
132+
},
133+
"outputs": [],
134+
"source": [
135+
"id(fam_rxn.reactants[0].molecule[0]), id(lib_rxn.reactants[0].molecule[0])"
136+
]
137+
},
138+
{
139+
"cell_type": "markdown",
140+
"metadata": {},
141+
"source": [
142+
"## step2: get fam_rxn's kinetics"
143+
]
144+
},
145+
{
146+
"cell_type": "markdown",
147+
"metadata": {},
148+
"source": [
149+
"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)."
150+
]
151+
},
152+
{
153+
"cell_type": "markdown",
154+
"metadata": {},
155+
"source": [
156+
"## step3: after training get fam_rxn's kinetics"
157+
]
158+
},
159+
{
160+
"cell_type": "code",
161+
"execution_count": null,
162+
"metadata": {
163+
"collapsed": true
164+
},
165+
"outputs": [],
166+
"source": [
167+
"cem = CoreEdgeReactionModel()\n",
168+
"cem.kineticsEstimator = 'rate rules'\n",
169+
"cem.verboseComments = True"
170+
]
171+
},
172+
{
173+
"cell_type": "code",
174+
"execution_count": null,
175+
"metadata": {
176+
"collapsed": false
177+
},
178+
"outputs": [],
179+
"source": [
180+
"from rmgpy.kinetics.kineticsdata import KineticsData\n",
181+
"from rmgpy.data.kinetics.family import TemplateReaction\n",
182+
"from rmgpy.data.kinetics.depository import DepositoryReaction\n",
183+
"\n",
184+
"for idx, spec in enumerate(fam_rxn.reactants):\n",
185+
" spec = Species(label=spec.label, molecule=spec.molecule)\n",
186+
" spec.generateThermoData(database)\n",
187+
" fam_rxn.reactants[idx] = spec\n",
188+
"for idx, spec in enumerate(fam_rxn.products):\n",
189+
" spec = Species(label=spec.label, molecule=spec.molecule)\n",
190+
" spec.generateThermoData(database)\n",
191+
" fam_rxn.products[idx] = spec\n",
192+
"\n",
193+
"family = getFamilyLibraryObject(fam_rxn.family)\n",
194+
"\n",
195+
"# If the reaction already has kinetics (e.g. from a library),\n",
196+
"# assume the kinetics are satisfactory\n",
197+
"if fam_rxn.kinetics is None:\n",
198+
" # Set the reaction kinetics\n",
199+
" kinetics, source, entry, isForward = cem.generateKinetics(fam_rxn)\n",
200+
" fam_rxn.kinetics = kinetics\n",
201+
" # Flip the reaction direction if the kinetics are defined in the reverse direction\n",
202+
" if not isForward:\n",
203+
" fam_rxn.reactants, fam_rxn.products = fam_rxn.products, fam_rxn.reactants\n",
204+
" fam_rxn.pairs = [(p,r) for r,p in fam_rxn.pairs]\n",
205+
" if family.ownReverse and hasattr(fam_rxn,'reverse'):\n",
206+
" if not isForward:\n",
207+
" fam_rxn.template = fam_rxn.reverse.template\n",
208+
" # We're done with the \"reverse\" attribute, so delete it to save a bit of memory\n",
209+
" delattr(fam_rxn,'reverse')\n",
210+
"\n",
211+
"# convert KineticsData to Arrhenius forms\n",
212+
"if isinstance(fam_rxn.kinetics, KineticsData):\n",
213+
" fam_rxn.kinetics = fam_rxn.kinetics.toArrhenius()\n",
214+
"# correct barrier heights of estimated kinetics\n",
215+
"if isinstance(fam_rxn,TemplateReaction) or isinstance(fam_rxn,DepositoryReaction): # i.e. not LibraryReaction\n",
216+
" fam_rxn.fixBarrierHeight() # also converts ArrheniusEP to Arrhenius.\n",
217+
"\n",
218+
"if cem.pressureDependence and fam_rxn.isUnimolecular():\n",
219+
" # If this is going to be run through pressure dependence code,\n",
220+
" # we need to make sure the barrier is positive.\n",
221+
" fam_rxn.fixBarrierHeight(forcePositive=True)"
222+
]
223+
},
224+
{
225+
"cell_type": "code",
226+
"execution_count": null,
227+
"metadata": {
228+
"collapsed": false
229+
},
230+
"outputs": [],
231+
"source": [
232+
"fam_rxn.kinetics"
233+
]
234+
},
235+
{
236+
"cell_type": "markdown",
237+
"metadata": {},
238+
"source": [
239+
"## step4: compare with lib_rxn's kinetics"
240+
]
241+
},
242+
{
243+
"cell_type": "code",
244+
"execution_count": null,
245+
"metadata": {
246+
"collapsed": false
247+
},
248+
"outputs": [],
249+
"source": [
250+
"lib_rxn.kinetics"
251+
]
252+
},
253+
{
254+
"cell_type": "markdown",
255+
"metadata": {},
256+
"source": [
257+
"## Conclusion: it improves the kinetics by factor of 10,000 at 673 for this reaction"
258+
]
259+
}
260+
],
261+
"metadata": {
262+
"kernelspec": {
263+
"display_name": "Python 2",
264+
"language": "python",
265+
"name": "python2"
266+
},
267+
"language_info": {
268+
"codemirror_mode": {
269+
"name": "ipython",
270+
"version": 2
271+
},
272+
"file_extension": ".py",
273+
"mimetype": "text/x-python",
274+
"name": "python",
275+
"nbconvert_exporter": "python",
276+
"pygments_lexer": "ipython2",
277+
"version": "2.7.11"
278+
}
279+
},
280+
"nbformat": 4,
281+
"nbformat_minor": 0
282+
}

0 commit comments

Comments
 (0)