reactionmechanismgenerator / rmg-py Goto Github PK
View Code? Open in Web Editor NEWPython version of the amazing Reaction Mechanism Generator (RMG).
Home Page: http://reactionmechanismgenerator.github.io/RMG-Py/
License: Other
rmg-py's People
Contributors
Stargazers
Watchers
Forkers
jwallen gmagoon athlonshi shamelmerchant pierrelb faribas nickvandewiele sean-v8 calebclass comocheng nyee bslakman jbarlow3 bbuesser alimalek raymondlanger brianwolfe agvandeputte stavlebet enochd chatelak wking rwest asilvi jabrownmit kehang alaraen cainja olusade zhugeyicixin mliu49 goldmanm connorcoley jimchu10 thackray kyleniemeyer dgreen18 alui1 alexpaquin cfgoldsmith pwstmrlnd alongd yunsiechung rgillis8 iyhsieh aelong cgrambow zjburas lyle-zhang soumsrani ionela89 pengzhang13 phalgunlolur alborzi wordplayer yepunaci ouchemical hamza-khan adeeljamal nateharms ngvjai wangyuran ray8201131 gitter-badger dbaugher jb-delafosse tafmuta fredwbacon mjohnson541 keceli yplitw wayneyann canh-thu c-underkoffler zhuzhu0713 solertis tybiot pattanaikl ajocher lily90502 nseyke mcs16 sarahkha chemphy zach-youssef kblondal bellamkondaprakash gyg16 nsm120 amarkpayne aspirincode ajuluc andreachlebikova drdew2 speth sleepyguyallday oscarwumit z5476t4508 minghao2016 dranasinghermg-py's Issues
Derivative continuity constraints between low- and high-temperature NASA polynomials
Currently, NASA polynomials created from Wilhoit polynomials are constrained to have continuous first and second derivatives at the intermediate temperature. Perhaps we should provide option(s) for relaxation of these constraints to allow discontinuous second-derivatives, and perhaps even discontinuous first derivatives, allowing a better fit.
Parallel Processing in RMG-Py
All the computers we will be running RMG on from now on have several cores, and many of them are even clusters. It would be nice to run code on all of them at once to reduce the wall-time of our calculations.
This issue is to keep a record of work towards parallelizing [parts of] RMG-Py.
output.html links to image names with # in, which don't load
All the species with triple bonds, the images don't load in output.html.
# should be replaced with %23 in the <img src
Singular matrix in rmgpy.measure.msc.applyModifiedStrongCollisionMethod
My methylformate job now dies with this:
Traceback (most recent call last):
File "/home/rwest/code/RMG-Py/rmg.py", line 134, in <module>
cProfile.runctx(command, global_vars, local_vars, stats_file)
File "/usr/lib/python2.6/cProfile.py", line 49, in runctx
prof = prof.runctx(statement, globals, locals)
File "/usr/lib/python2.6/cProfile.py", line 140, in runctx
exec cmd in globals, locals
File "<string>", line 1, in <module>
File "/home/rwest/code/RMG-Py/rmgpy/rmg/main.py", line 333, in execute
self.reactionModel.enlarge(object)
File "/home/rwest/code/RMG-Py/rmgpy/rmg/model.py", line 577, in enlarge
self.updateUnimolecularReactionNetworks(database)
File "/home/rwest/code/RMG-Py/rmgpy/rmg/model.py", line 1038, in updateUnimolecularReactionNetworks
network.update(self, database, self.pressureDependence)
File "/home/rwest/code/RMG-Py/rmgpy/rmg/pdep.py", line 400, in update
K, p0 = self.calculateRateCoefficients(Tlist, Plist, Elist, method)
File "/home/rwest/code/RMG-Py/rmgpy/measure/network.py", line 591, in calculateRateCoefficients
K[t,p,:,:], p0[t,p,:,:,:] = msc.applyModifiedStrongCollisionMethod(T, P, Elist, densStates, collFreq, Kij, Fim, Gnj, Ereac, Nisom, Nreac, Nprod)
File "msc.pyx", line 131, in rmgpy.measure.msc.applyModifiedStrongCollisionMethod (build/pyrex/rmgpy/measure/msc.c:1914)
File "/usr/lib/python2.6/dist-packages/numpy/linalg/linalg.py", line 254, in solve
raise LinAlgError, 'Singular matrix'
numpy.linalg.linalg.LinAlgError: Singular matrixThe MEASURE file in question is at http://rmg.mit.edu/measure/networks/0a245b5d26cbcc95f466c7d3d378f8ec
LinAlgError, 'Singular matrix' in rmgpy.measure.msc.applyModifiedStrongCollisionMethod
Traceback (most recent call last):
File "/home/rwest/code/RMG-Py/rmg.py", line 134, in <module>
cProfile.runctx(command, global_vars, local_vars, stats_file)
File "/usr/lib/python2.6/cProfile.py", line 49, in runctx
prof = prof.runctx(statement, globals, locals)
File "/usr/lib/python2.6/cProfile.py", line 140, in runctx
exec cmd in globals, locals
File "<string>", line 1, in <module>
File "/home/rwest/code/RMG-Py/rmgpy/rmg/main.py", line 337, in execute
self.reactionModel.enlarge(object)
File "/home/rwest/code/RMG-Py/rmgpy/rmg/model.py", line 600, in enlarge
self.updateUnimolecularReactionNetworks(database)
File "/home/rwest/code/RMG-Py/rmgpy/rmg/model.py", line 1106, in updateUnimolecularReactionNetworks
network.update(self, database, self.pressureDependence)
File "/home/rwest/code/RMG-Py/rmgpy/rmg/pdep.py", line 413, in update
K, p0 = self.calculateRateCoefficients(Tlist, Plist, Elist, method)
File "/home/rwest/code/RMG-Py/rmgpy/measure/network.py", line 591, in calculateRateCoefficients
K[t,p,:,:], p0[t,p,:,:,:] = msc.applyModifiedStrongCollisionMethod(T, P, Elist, densStates, collFreq, Kij, Fim, Gnj, Ereac, Nisom, Nreac, Nprod)
File "msc.pyx", line 131, in rmgpy.measure.msc.applyModifiedStrongCollisionMethod (build/pyrex/rmgpy/measure/msc.c:1914)
File "/usr/lib/python2.6/dist-packages/numpy/linalg/linalg.py", line 254, in solve
raise LinAlgError, 'Singular matrix'The network in question is at
http://rmg.mit.edu/measure/networks/7977a99b58d3185814d1a62926b68241
Here is the summary:
========================================================================
Network Information (Network #5485)
-------------------
Isomers:
CC=C=C=O(5168) 123.55 kJ/mol
C(=C[C]=O)[CH2](2225) 145.259 kJ/mol
C(=C=C[O])[CH2](5407) 215.06 kJ/mol
C[C]=C[C]=O(2112) 231.451 kJ/mol
C([CH2])C#C[O](503) 337.193 kJ/mol
C[C]=C=C[O](5406) 300.862 kJ/mol
C(=C=C)C=O(2357) 41.2816 kJ/mol
Reactant channels:
CH3(10) + C(#C)[C]=O(313) 375.915 kJ/mol
Product channels:
H(17) + C(=C=C=O)[CH2](5404) 466.065 kJ/mol
H(17) + C[C]=C=C=O(5405) 407.34 kJ/mol
H(17) + C(=C=C)[C]=O(2660) 406.73 kJ/mol
H(17) + C(=[C][CH2])[C]=O(2657) 593.893 kJ/mol
H(17) + C(=[C][C]=O)[CH2](5008) 554.88 kJ/mol
CC=[C][C]=O(5166) 192.299 kJ/mol
C(=[C][CH2])C=O(2501) 254.042 kJ/mol
[CH]([CH2])[C]=C[O](5408) 539.5 kJ/mol
C(=C)C=C=O(5409) 17.5958 kJ/mol
H(17) + [CH](C=[CH])[C]=O(3647) 575.177 kJ/mol
C([C]=C)[C]=O(2358) 300.62 kJ/mol
C(C=[CH])[C]=O(2224) 309.905 kJ/mol
C2H3(33) + HCCO(38) 437.405 kJ/mol
C2H3(33) + [C](=[CH])[O](5410) 718.543 kJ/mol
C([CH2])[C]=[C][O](4403) 581.781 kJ/mol
C(=C[O])C=[CH](5411) 235.33 kJ/mol
[CH](C=[C][O])[CH2](5412) 494.254 kJ/mol
H(17) + [CH2]C#CC=O(2895) 408.851 kJ/mol
C(=O)C#C(177) + CH2(24) 509.337 kJ/mol
H(17) + C(=O)[C]=[C][CH2](5413) 663.54 kJ/mol
CC#CC=O(5414) 41.9879 kJ/mol
C(=C)C#C(183) + O(18) 529.112 kJ/mol
H(17) + C(=[CH])[C]=C[O](5415) 644.738 kJ/mol
C(=C)C#CO(5416) 129.269 kJ/mol
C(=C=C=C)O(5417) 144.634 kJ/mol
C(=C)[C]=[C]O(5418) 359.554 kJ/mol
O(18) + C(=C)[C]=[CH](3084) 694.995 kJ/mol
C(=[C][CH2])[CH][O](5419) 594.688 kJ/mol
H(17) + C[C]=[C][C]=O(5420) 640.665 kJ/mol
CO(7) + C[C]=[CH](2053) 371.665 kJ/mol
CH2(24) + C=C=C=O(4985) 507.188 kJ/mol
[CH2][CH2](80) + C2O(41) 693.815 kJ/mol
O(18) + C([CH2])C#[C](502) 934.281 kJ/mol
C[C]=C=[C]O(5421) 444.533 kJ/mol
O(18) + C[C]=C=[CH](3082) 781.237 kJ/mol
CH3(10) + [CH](C#[C])[O](5422) 598.523 kJ/mol
C(C#C[CH2])[O](2991) 357.799 kJ/mol
H(17) + C(=C=[CH])C=O(2658) 406.437 kJ/mol
C(C=C=[CH])[O](5423) 361.658 kJ/mol
Path reactions:
CH3(10) + C(#C)[C]=O(313) <=> CC=C=C=O(5168) 375.915 kJ/mol
H(17) + C(=C=C=O)[CH2](5404) <=> CC=C=C=O(5168) 465.727 kJ/mol
H(17) + C[C]=C=C=O(5405) <=> CC=C=C=O(5168) 407.456 kJ/mol
C([CH2])C#C[O](503) <=> CC=C=C=O(5168) 376.941 kJ/mol
C[C]=C[C]=O(2112) <=> CC=C=C=O(5168) 271.199 kJ/mol
C(=C[C]=O)[CH2](2225) <=> CC=C=C=O(5168) 213.458 kJ/mol
C[C]=C=C[O](5406) <=> CC=C=C=O(5168) 369.061 kJ/mol
C(=C=C[O])[CH2](5407) <=> CC=C=C=O(5168) 247.486 kJ/mol
H(17) + C(=C=C)[C]=O(2660) <=> C(=C[C]=O)[CH2](2225) 419.282 kJ/mol
H(17) + C(=C=C=O)[CH2](5404) <=> C(=C[C]=O)[CH2](2225) 479.881 kJ/mol
H(17) + C(=[C][CH2])[C]=O(2657) <=> C(=C[C]=O)[CH2](2225) 593.893 kJ/mol
H(17) + C(=[C][C]=O)[CH2](5008) <=> C(=C[C]=O)[CH2](2225) 554.88 kJ/mol
C(=C[C]=O)[CH2](2225) <=> C(=C=C)C=O(2357) 213.458 kJ/mol
C[C]=C[C]=O(2112) <=> C(=C[C]=O)[CH2](2225) 402.312 kJ/mol
C(=C=C[O])[CH2](5407) <=> C(=C[C]=O)[CH2](2225) 382.282 kJ/mol
CC=[C][C]=O(5166) <=> C(=C[C]=O)[CH2](2225) 338.214 kJ/mol
C(=[C][CH2])C=O(2501) <=> C(=C[C]=O)[CH2](2225) 403.91 kJ/mol
[CH]([CH2])[C]=C[O](5408) <=> C(=C[C]=O)[CH2](2225) 579.248 kJ/mol
C([CH2])C#C[O](503) <=> C(=C[C]=O)[CH2](2225) 495.803 kJ/mol
C(=C[C]=O)[CH2](2225) <=> C(=C)C=C=O(5409) 145.259 kJ/mol
H(17) + [CH](C=[CH])[C]=O(3647) <=> C(=C[C]=O)[CH2](2225) 575.177 kJ/mol
C([C]=C)[C]=O(2358) <=> C(=C[C]=O)[CH2](2225) 467.842 kJ/mol
C(C=[CH])[C]=O(2224) <=> C(=C[C]=O)[CH2](2225) 434.328 kJ/mol
C2H3(33) + HCCO(38) <=> C(=C[C]=O)[CH2](2225) 449.68 kJ/mol
C2H3(33) + [C](=[CH])[O](5410) <=> C(=C[C]=O)[CH2](2225) 719.054 kJ/mol
C([CH2])[C]=[C][O](4403) <=> C(=C[C]=O)[CH2](2225) 621.529 kJ/mol
C(=C[O])C=[CH](5411) <=> C(=C[C]=O)[CH2](2225) 384.034 kJ/mol
[CH](C=[C][O])[CH2](5412) <=> C(=C[C]=O)[CH2](2225) 494.254 kJ/mol
H(17) + [CH2]C#CC=O(2895) <=> C(=C=C[O])[CH2](5407) 418.278 kJ/mol
H(17) + C(=C=C=O)[CH2](5404) <=> C(=C=C[O])[CH2](5407) 479.881 kJ/mol
C(=O)C#C(177) + CH2(24) <=> C(=C=C[O])[CH2](5407) 528.499 kJ/mol
H(17) + C(=O)[C]=[C][CH2](5413) <=> C(=C=C[O])[CH2](5407) 663.54 kJ/mol
H(17) + C(=[C][C]=O)[CH2](5008) <=> C(=C=C[O])[CH2](5407) 554.88 kJ/mol
C(=C=C[O])[CH2](5407) <=> CC#CC=O(5414) 254.808 kJ/mol
C(=C=C[O])[CH2](5407) <=> C(=C=C)C=O(2357) 254.808 kJ/mol
C[C]=C=C[O](5406) <=> C(=C=C[O])[CH2](5407) 463.111 kJ/mol
C(=[C][CH2])C=O(2501) <=> C(=C=C[O])[CH2](5407) 504.303 kJ/mol
C(=C)C#C(183) + O(18) <=> C(=C=C[O])[CH2](5407) 548.274 kJ/mol
H(17) + C(=[CH])[C]=C[O](5415) <=> C(=C=C[O])[CH2](5407) 644.738 kJ/mol
C(=C=C[O])[CH2](5407) <=> C(=C)C=C=O(5409) 254.808 kJ/mol
C(=C=C[O])[CH2](5407) <=> C(=C)C#CO(5416) 254.808 kJ/mol
C(=C=C[O])[CH2](5407) <=> C(=C=C=C)O(5417) 247.486 kJ/mol
C(=C)[C]=[C]O(5418) <=> C(=C=C[O])[CH2](5407) 518.164 kJ/mol
C(=C[O])C=[CH](5411) <=> C(=C=C[O])[CH2](5407) 394.829 kJ/mol
O(18) + C(=C)[C]=[CH](3084) <=> C(=C=C[O])[CH2](5407) 694.995 kJ/mol
C(=[C][CH2])[CH][O](5419) <=> C(=C=C[O])[CH2](5407) 634.436 kJ/mol
[CH](C=[C][O])[CH2](5412) <=> C(=C=C[O])[CH2](5407) 534.002 kJ/mol
H(17) + C(=C=C)[C]=O(2660) <=> C[C]=C[C]=O(2112) 411.751 kJ/mol
H(17) + C[C]=C=C=O(5405) <=> C[C]=C[C]=O(2112) 417.704 kJ/mol
H(17) + C(=[C][CH2])[C]=O(2657) <=> C[C]=C[C]=O(2112) 593.893 kJ/mol
H(17) + C[C]=[C][C]=O(5420) <=> C[C]=C[C]=O(2112) 640.665 kJ/mol
C[C]=C[C]=O(2112) <=> CC#CC=O(5414) 271.199 kJ/mol
C[C]=C[C]=O(2112) <=> C(=C=C)C=O(2357) 263.877 kJ/mol
C[C]=C[C]=O(2112) <=> CC=[C][C]=O(5166) 385.786 kJ/mol
C[C]=C=C[O](5406) <=> C[C]=C[C]=O(2112) 502.064 kJ/mol
CO(7) + C[C]=[CH](2053) <=> C[C]=C[C]=O(2112) 392.585 kJ/mol
H(17) + C(=C=C=O)[CH2](5404) <=> C([CH2])C#C[O](503) 480.925 kJ/mol
CH2(24) + C=C=C=O(4985) <=> C([CH2])C#C[O](503) 530.497 kJ/mol
H(17) + C(=[C][C]=O)[CH2](5008) <=> C([CH2])C#C[O](503) 554.88 kJ/mol
C([CH2])C#C[O](503) <=> C(=C)C=C=O(5409) 376.941 kJ/mol
C([CH2])C#C[O](503) <=> CC=[C][C]=O(5166) 507.066 kJ/mol
[CH2][CH2](80) + C2O(41) <=> C([CH2])C#C[O](503) 693.745 kJ/mol
C([CH2])C#C[O](503) <=> C(=C)C#CO(5416) 369.619 kJ/mol
[CH](C=[C][O])[CH2](5412) <=> C([CH2])C#C[O](503) 534.002 kJ/mol
C(=C)[C]=[C]O(5418) <=> C([CH2])C#C[O](503) 523.187 kJ/mol
O(18) + C([CH2])C#[C](502) <=> C([CH2])C#C[O](503) 934.281 kJ/mol
H(17) + [CH2]C#CC=O(2895) <=> C[C]=C=C[O](5406) 419.329 kJ/mol
H(17) + C[C]=C=C=O(5405) <=> C[C]=C=C[O](5406) 421.156 kJ/mol
H(17) + C(=O)[C]=[C][CH2](5413) <=> C[C]=C=C[O](5406) 663.54 kJ/mol
H(17) + C[C]=[C][C]=O(5420) <=> C[C]=C=C[O](5406) 640.665 kJ/mol
C[C]=C=C[O](5406) <=> C(=C=C=C)O(5417) 316.971 kJ/mol
C[C]=C=[C]O(5421) <=> C[C]=C=C[O](5406) 603.143 kJ/mol
O(18) + C[C]=C=[CH](3082) <=> C[C]=C=C[O](5406) 781.237 kJ/mol
CH3(10) + [CH](C#[C])[O](5422) <=> C[C]=C=C[O](5406) 598.888 kJ/mol
[CH]([CH2])[C]=C[O](5408) <=> C[C]=C=C[O](5406) 579.248 kJ/mol
C(C#C[CH2])[O](2991) <=> C[C]=C=C[O](5406) 454.728 kJ/mol
H(17) + [CH2]C#CC=O(2895) <=> C(=C=C)C=O(2357) 408.967 kJ/mol
H(17) + C(=C=C)[C]=O(2660) <=> C(=C=C)C=O(2357) 406.409 kJ/mol
H(17) + C(=C=[CH])C=O(2658) <=> C(=C=C)C=O(2357) 406.437 kJ/mol
C(C#C[CH2])[O](2991) <=> C(=C=C)C=O(2357) 397.547 kJ/mol
C([C]=C)[C]=O(2358) <=> C(=C=C)C=O(2357) 340.368 kJ/mol
C(=C[O])C=[CH](5411) <=> C(=C=C)C=O(2357) 275.078 kJ/mol
C(C=C=[CH])[O](5423) <=> C(=C=C)C=O(2357) 394.084 kJ/mol
========================================================================
Using 455 grains from 41.28 to 1940.82 kJ/mol in steps of 4.18 kJ/mol
Calculating densities of states for Network #5485...
Calculating phenomenological rate coefficients for Network #5485...
Reaction flux diagrams
Reaction path / flux diagrams would be nice.
I think you can ask something like this of Cantera:
# Write reaction path diagram
d = rxnpath.PathDiagram(threshold = 0.10, detailed = 'true',flow_type='both')
rxnpath.write(gas,'C','filename.dot',d)
Misinterprets SMILES string '[CH]' as having 0 radical Electrons.
When making a molecule from the SMILES [CH] the carbon ends up with 0 radicalElectrons instead of 3.
(There does not seem to be a problem with [CH2] or [CH3])
I have added a (failing) unit test in ca3b3bc
I notice that ChemRes make a complete hash of these radicals: http://cactus.nci.nih.gov/chemical/structure/%5BCH%5D/formula
But the syntax at http://www.opensmiles.org/spec/open-smiles-3-input.html is clear enough.
Molecule(SMILES='[H]').toAdjacencyList() gives 1 H 0
the output of Molecule(SMILES='[H]').toAdjacencyList() should be 1 H 1 (note the radical) not 1 H 0
put species.speciesCounter somewhere where it'll get pickled
currently species.speciesCounter is reset to be the length of the speciesList after a restart,
which is not safe.
Can't cope with missing bath gas when running in PDep mode.
I am running with:
# Data sources
database(
'../RMG-database/input',
thermoLibraries = ['primaryThermoLibrary','DFT_QCI_thermo','GRI-Mech3.0'],
reactionLibraries = ['Methylformate','Glarborg/highP'],
seedMechanisms = ['GRI-Mech3.0'],
)
pressureDependence(
method='modified strong collision',
minimumGrainSize=(2.0,'kcal/mol'),
minimumNumberOfGrains=200,
temperatures=(300,'K',2000,'K',8),
pressures=(0.01,'bar',100,'bar',5),
interpolation=('Chebyshev', 6, 4),
)I have a lot of lines saying:
Calculating densities of states for Network #1...
Calculating phenomenological rate coefficients for Network #1...
Exception ZeroDivisionError: 'float division' in 'rmgpy.measure.collision.calculateCollisionFrequency' ignored
Exception ZeroDivisionError: 'float division' in 'rmgpy.measure.collision.calculateCollisionFrequency' ignored
Exception ZeroDivisionError: 'float division' in 'rmgpy.measure.collision.calculateCollisionFrequency' ignored
Exception ZeroDivisionError: 'float division' in 'rmgpy.measure.collision.calculateCollisionFrequency' ignored
...etc...
then it dies with:
Traceback (most recent call last):
File "rmg.py", line 711, in <module>
execute(args)
File "rmg.py", line 275, in execute
reactionModel.enlarge(spec, database)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 938, in enlarge
self.updateUnimolecularReactionNetworks(database)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 1399, in updateUnimolecularReactionNetworks
network.update(self, database)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 532, in update
netReaction.kinetics = fitInterpolationModel(netReaction, Tlist, Plist, K[:,:,i,j], model, Tmin, Tmax, Pmin, Pmax)
File "reaction.pyx", line 332, in rmgpy.measure.reaction.fitInterpolationModel (build/pyrex/rmgpy/measure/reaction.c:5141)
File "kinetics.py", line 732, in rmgpy.kinetics.Chebyshev.fitToData (build/pyrex/rmgpy/kinetics.c:12744)
File "kinetics.py", line 767, in rmgpy.kinetics.Chebyshev.fitToData (build/pyrex/rmgpy/kinetics.c:12358)
ValueError: math domain errorA little up the stack trace is :
netReaction.kinetics = fitInterpolationModel(netReaction, Tlist, Plist, K[:,:,i,j], model, Tmin, Tmax, Pmin, Pmax)at which point the matrix K is entirely zero. If it helps, this is the netReaction:
Reaction(index=334, reactants=[Species(index=1, label="methylformate", thermo=MultiNASA(Tmin=100, Tmax=5000, polynomials=[NASA(Tmin=(100,"K"), Tmax=(1133.91,"K"), coeffs=[3.26245,0.0112056,1.91547e-05,-2.51712e-08,8.14166e-12,-44760.2,12.148], comment="""Low temperature range polynomial"""),NASA(Tmin=(1133.91,"K"), Tmax=(5000,"K"), coeffs=[5.13303,0.0182987,-8.34096e-06,1.64356e-09,-1.18236e-13,-46064.6,-0.992883], comment="""High temperature range polynomial""")], comment="""NASA function fitted to Wilhoit function. Weighted RMS error = 0.042_R;(Unweighted) RMS error = 0.046_R;"""), states=StatesModel(modes=[HarmonicOscillator(frequencies=([2782.5,750,1395,475,1775,1000,2750,2800,2850,1350,1500,750,1050,1375,1000,425],"cm^-1")), HinderedRotor(inertia=(0.0877936,"amu_angstrom^2"), barrier=(2.01855,"kJ/mol"), symmetry=1), HinderedRotor(inertia=(0.0654664,"amu_angstrom^2"), barrier=(47.6296,"kJ/mol"), symmetry=1)], spinMultiplicity=1), molecule=[Molecule(SMILES="COC=O")], E0=(-372.379,"kJ/mol"), lennardJones=LennardJones(sigma=(4.687e-10,"m"), epsilon=(7.33678e-21,"J")), molecularWeight=(60.052,"g/mol"))], products=[Species(index=16, label="CH3OH", thermo=MultiNASA(Tmin=100, Tmax=5000, polynomials=[NASA(Tmin=(100,"K"), Tmax=(1143.22,"K"), coeffs=[3.89247,-0.000484239,2.42754e-05,-2.44333e-08,7.4633e-12,-25699.7,5.86289], comment="""Low temperature range polynomial"""),NASA(Tmin=(1143.22,"K"), Tmax=(5000,"K"), coeffs=[5.08762,0.00759196,-2.40448e-06,5.03905e-10,-4.09468e-14,-26774,-3.56642], comment="""High temperature range polynomial""")], comment="""NASA function fitted to Wilhoit function. Weighted RMS error = 0.028_R;(Unweighted) RMS error = 0.030_R;"""), molecule=[Molecule(SMILES="CO")], E0=(-213.698,"kJ/mol"), lennardJones=LennardJones(sigma=(4.443e-10,"m"), epsilon=(1.52838e-21,"J")), molecularWeight=(32.0419,"g/mol")), Species(index=25, label="CO", thermo=MultiNASA(Tmin=100, Tmax=5000, polynomials=[NASA(Tmin=(100,"K"), Tmax=(1552.61,"K"), coeffs=[3.65108,-0.00132484,3.28978e-06,-2.09102e-09,4.32436e-13,-14359.9,3.18576], comment="""Low temperature range polynomial"""),NASA(Tmin=(1552.61,"K"), Tmax=(5000,"K"), coeffs=[3.09568,0.00139239,-5.78124e-07,1.03416e-10,-6.83094e-15,-14342.5,5.61029], comment="""High temperature range polynomial""")], comment="""NASA function fitted to Wilhoit function. Weighted RMS error = 0.006_R;(Unweighted) RMS error = 0.003_R;"""), molecule=[Molecule(SMILES="[C]=O")], E0=(-119.306,"kJ/mol"), lennardJones=LennardJones(sigma=(4.443e-10,"m"), epsilon=(1.52838e-21,"J")), molecularWeight=(28.0101,"g/mol"))], reversible=False)
(edit: title changed after comment below added)
Modified Strong Collision PDep gives negative steady-state concentration (and dies)
When I run my methyl formate test case with modified strong collision pressure dependence, it dies in the following way:
========================================================================
Network Information (Network #126)
-------------------
Isomers:
C3H4O3(167) -133.858 kJ/mol
Reactant channels:
Mfmt(1) + CO(7) -491.686 kJ/mol
Product channels:
CH3(10) + C2HO3(180) -143.6 kJ/mol
CO2(8) + C2H4O(542) -138.391 kJ/mol
CH3(10) + C2HO3(520) 204.839 kJ/mol
H(30) + C3H3O3(442) 264.895 kJ/mol
C3H4O3(543) -151.111 kJ/mol
C3H4O3(544) -384.873 kJ/mol
CO(7) + C2H4O2(523) -74.267 kJ/mol
Path reactions:
Mfmt(1) + CO(7) <=> C3H4O3(167) -462.821 kJ/mol
CH3(10) + C2HO3(180) <=> C3H4O3(167) -107.217 kJ/mol
CO2(8) + C2H4O(542) <=> C3H4O3(167) -109.526 kJ/mol
CH3(10) + C2HO3(520) <=> C3H4O3(167) 205.413 kJ/mol
H(30) + C3H3O3(442) <=> C3H4O3(167) 264.557 kJ/mol
C3H4O3(167) <=> C3H4O3(543) 12.5529 kJ/mol
C3H4O3(167) <=> C3H4O3(544) -120.812 kJ/mol
CO(7) + C2H4O2(523) <=> C3H4O3(167) -61.715 kJ/mol
========================================================================
Using 226 grains from -491.69 to 1391.11 kJ/mol in steps of 8.37 kJ/mol
Calculating densities of states for Network #126...
Calculating phenomenological rate coefficients for Network #126...
Warning: invalid value encountered in divide
Traceback (most recent call last):
File "rmg.py", line 732, in <module>
cProfile.runctx(command, global_vars, local_vars, stats_file)
File "/usr/local/Cellar/python/2.7.1/Frameworks/Python.framework/Versions/2.7/lib/python2.7/cProfile.py", line 49, in runctx
prof = prof.runctx(statement, globals, locals)
File "/usr/local/Cellar/python/2.7.1/Frameworks/Python.framework/Versions/2.7/lib/python2.7/cProfile.py", line 140, in runctx
exec cmd in globals, locals
File "<string>", line 1, in <module>
File "rmg.py", line 359, in execute
reactionModel.enlarge(object, database)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 942, in enlarge
self.updateUnimolecularReactionNetworks(database)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 1403, in updateUnimolecularReactionNetworks
network.update(self, database)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 501, in update
K, p0 = self.calculateRateCoefficients(Tlist, Plist, Elist, method)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/measure/network.py", line 529, in calculateRateCoefficients
K[t,p,:,:], p0[t,p,:,:,:] = msc.applyModifiedStrongCollisionMethod(T, P, Elist, densStates, collFreq, Kij, Fim, Gnj, Ereac, Nisom, Nreac, Nprod)
File "msc.pyx", line 139, in rmgpy.measure.msc.applyModifiedStrongCollisionMethod (build/pyrex/rmgpy/measure/msc.c:2023)
rmgpy.measure.msc.ModifiedStrongCollisionError: A negative steady-state concentration was encountered.Poor documentation for installation and setup
Especially given how important cythonising is for speeding up the thermo,
some documentation on how to do this would be great (especially for windows users)
Extrapolating Cp to high temperatures inaccurate with linear interpolation
Would be better to use Wilhoit form.
ctml_writer.py interferes with cantera installation
I can't remember the details, but there's a reason we moved aside ctml_writer.py on greg's computer. I think it was being loaded by mistake when he tried to use cantera just as cantera, and didn't want our modified version. We should make sure they don't collide; perhaps rename ours?
Perhaps Greg can play around with cantera and figure out what was wrong.
Subgraph and equivalent() is wrong and confusing
Currently the equivalent() methods match for
[Cd,Cs].equivalent([Cd,Cb,H]) = True
[Cd,Cb,H].equivalent([Cd,Cs]) = True
[R].equivalent([Cs]) = True
[Cs].equivalent([R]) = True
etc.
and the subgraph matching is confusingly named too.
I propose methods called self.isSpecificCaseOf(other) that is true if they are identical or self is fully contained within other (other is equal to or more general than self).
We may also need exact matching, such as self.isIdenticalTo(other) that is true only if they match exactly.
Installing RMGpy (and dependencies) is a pain!
Cantera 1.8 is especially troublesome.
Should we just put pressure on the cantera devs for a binary?
Perhaps we can build our own and host it somewhere.
(certainly will need to if we modify cantera)
DDASSL fails with IDID = -8 (matrix is singular?)
After running for half an hour (on pharos) my methylformate job in RMG-Py dies as follows:
After model enlargement:
The model core has 67 species and 1055 reactions
The model edge has 2095 species and 4117 reactions
Saving latest model core to HTML file...
Saving latest model core to Chemkin file...
Chemkin file contains 1110 reactions.
Updating RMG execution statistics...
Execution time (HH:MM:SS): 00:32:23
Memory used: 250.63 MB
Conducting simulation of reaction system 1...
dbolsm(). more than (i1)=itmax iterations solving bounded least squares problem.
error number = 22, message level = 1
i1 = 50
dbolsm(). more than (i1)=itmax iterations solving bounded least squares problem.
error number = 22, message level = 1
i1 = 50
dbolsm(). more than (i1)=itmax iterations solving bounded least squares problem.
error number = 22, message level = 1
i1 = 50
***MESSAGE FROM ROUTINE DDASSL IN LIBRARY SLATEC.
***POTENTIALLY RECOVERABLE ERROR, PROGRAM CONTINUES.
* AT T = 4.701195D-01 AND STEPSIZE H = 1.185978D-04 THE ITERATION
* MATRIX IS SINGULAR
* ERROR NUMBER = -8
*
***END OF MESSAGE
***MESSAGE FROM ROUTINE DDASSL IN LIBRARY SLATEC.
***FATAL ERROR, PROGRAM ABORTED.
* THE LAST STEP TERMINATED WITH A NEGATIVE VALUE OF IDID = -8 AND NO
* APPROPRIATE ACTION WAS TAKEN. RUN TERMINATED
* ERROR NUMBER = -998
*
***END OF MESSAGE
***JOB ABORT DUE TO FATAL ERROR.
The first message suggests it's potentially recoverable, but we do nothing to recover from it.
Perhaps we have just made some super-stiff ODEs by messing up rate estimates. Either way, it's a problem..
Reduced numerical accuracy in cythonized thermo.py integrals
Numerical error seems to accumulate more readily within cythonized integrals and/or using computations based on results of cythonized integrals. As a result, when run in cython mode, the integral of squared errors can sometimes be substantially negative, when, in fact, the result must be non-negative. Example cases that cause such problems appear in thermotest.py unit tests. These cases do not cause such problems when run in pure-python mode. (There is also some evidence that the cython results are platform-dependent.)
KeyError: 'Node not found in database.' when importing GRI-Mech seed
A job with this:
# Data sources
database(
'../RMG-database/input',
thermoLibraries = ['primaryThermoLibrary'],
reactionLibraries = ['Methylformate'],
seedMechanisms = ['GRI-Mech3.0'],
)fails with this:
Adding seed mechanism GRI-Mech3.0 to model core...
Traceback (most recent call last):
File "rmg.py", line 673, in <module>
execute(args)
File "rmg.py", line 262, in execute
reactionModel.addSeedMechanismToCore(database, seedMechanism, react=False)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 1317, in addSeedMechanismToCore
if spec.reactive: spec.generateThermoData(database)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 93, in generateThermoData
tdata = database.thermo.getThermoData(molecule)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/data/thermo.py", line 518, in getThermoData
thermoData = self.getThermoDataFromGroups(molecule)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/data/thermo.py", line 629, in getThermoDataFromGroups
thermoData = self.__addThermoData(thermoData, self.__getGroupThermoData(self.groups['radical'], saturatedStruct, {'*':atom}))
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/data/thermo.py", line 720, in __getGroupThermoData
raise KeyError('Node not found in database.')
KeyError: 'Node not found in database.'toChemkin() in thermo is formatted incorrectly when run on Windows
On my system, the toChemkin() command will display exponents with three digits (e+000) rather than the desired two digits (e+00). On other (Linux) systems, it seems to work as desired.
PDep reactions: both directions are not getting into the core
PDep reactions are now irreversible, but when things get into the core you would expect therefore to have both directions listed. This seems not to be the case.
Images can't distinguish between spin state of biradicals (2S or 2T)
The images don't distinguish between 2T and 2S electron states.
This picture:
http://rmg.mit.edu/molecule/1%20C%202T;
looks the same as this picture:
http://rmg.mit.edu/molecule/1%20C%202S;
Overall symmetry number errors
Symmetry number unit tests (9b73cdd) for ethylene and formaldehyde fail (results are too low by a factor of two).
cython is not creating .html
the -a flag for the cython command line option should still work, but it'd be nice to have it in our setup.py
I think something like this is what we need, but have not been able to get it working myself yet:
build_ext.pyrex_default_options['annotate']=True
cantera 1.7 installation warning
We should put in the installation structions a warning about the cantera 1.7 screwup
A quote from cantera mailing list in December 2007:
This line out of Makfile.in
rm -fR @ct_libdir@/*
cause some MINOR problems, if you type
./configure --prefix=/usr
because in configure the line
ct_libdir=${prefix}/lib
isn't that optimal.
If you didn't get it till now: I lost all user libraries!
e85 job causes AttributeError: 'rmg.chem.Atom' object has no attribute 'breakBond'
This is with 16f8145 ::
RMG execution initiated at Thu Nov 12 17:58:48 2009
################################################################
# #
# RMG - Reaction Mechanism Generator #
# Python Version 0.0.1 #
# 14 May 2009 #
# #
# http://rmg.sourceforge.net/ #
# #
################################################################
Found 1 database
Traceback (most recent call last):
File "rmg.py", line 149, in <module>
rmg.main.execute( args[0], options )
File "/Users/rwest/XCodeProjects/RMGpy/source/rmg/main.py", line 96, in execute
reactionModel, coreSpecies, reactionSystems = io.readInputFile(inputFile)
File "/Users/rwest/XCodeProjects/RMGpy/source/rmg/io.py", line 295, in readInputFile
reaction.kineticsDatabase.load(database[2] + os.sep)
File "/Users/rwest/XCodeProjects/RMGpy/source/rmg/reaction.py", line 1650, in load
family.load(path)
File "/Users/rwest/XCodeProjects/RMGpy/source/rmg/reaction.py", line 636, in load
self.loadTemplate(tempstr)
File "/Users/rwest/XCodeProjects/RMGpy/source/rmg/reaction.py", line 1156, in loadTemplate
self.generateProductTemplate()
File "/Users/rwest/XCodeProjects/RMGpy/source/rmg/reaction.py", line 1191, in generateProductTemplate
productStructure = self.applyRecipe(reactantStructure, unique=False)
File "/Users/rwest/XCodeProjects/RMGpy/source/rmg/reaction.py", line 1277, in applyRecipe
if not self.recipe.applyForward(reactantStructure, unique):
File "/Users/rwest/XCodeProjects/RMGpy/source/rmg/reaction.py", line 529, in applyForward
return self.__apply(struct, True, unique)
File "/Users/rwest/XCodeProjects/RMGpy/source/rmg/reaction.py", line 490, in __apply
atom1.breakBond(struct.getBonds(atom1), unique)
AttributeError: 'rmg.chem.Atom' object has no attribute 'breakBond'
Many undeclared duplicate reactions in chemkin file
After 2.5 hours, my methylformate job reported:
The model core has 101 species and 4708 reactions
The model edge has 11324 species and 17960 reactions
When I try to run it through chemkin it finds 1230 undeclared duplicates. The ones I've checked all seem to be in the reverse of the direction which they duplicate.
Interestingly there are only 210 chebyshev reactions (a small fraction of 4708), and these all printed as irreversible (with => arrows).
Simulation (ODE solving) is too slow
something broken
this is a bad message - it's not specific
NetworkError: Unexpected type of path reaction in measure/network.py(382)
When running RMG with PDep for my methyl formate test case, I end up with a NetworkError:
NetworkError: Unexpected type of path reaction "C2H3O(74) <=> CO(7) + CH3(10)"
in rmgpy/measure/network.py(382)calculateMicrocanonicalRates()The network in question is
========================================================================
Network Information
-------------------
Isomers:
CH2CHO(34) 8.30956 kJ/mol
Reactant channels:
CO(7) + CH3(10) 17.7746 kJ/mol
Product channels:
cC2H3O(239) 158.188 kJ/mol
H(30) + CH2CO(29) 153.561 kJ/mol
H(30) + C2H2O(83) 394.282 kJ/mol
C2H3O(74) -19.9272 kJ/mol
H(30) + C2H2O(67) 362.957 kJ/mol
C2H3O(75) 107.003 kJ/mol
C2H3O(241) 122.318 kJ/mol
O(17) + C2H3(25) 534.417 kJ/mol
Path reactions:
CH2CHO(34) <=> CO(7) + CH3(10) 177.034 kJ/mol
cC2H3O(239) <=> CH2CHO(34) 220.923 kJ/mol
CH2CHO(34) <=> H(30) + CH2CO(29) 199.125 kJ/mol
H(30) + C2H2O(83) <=> CH2CHO(34) 394.026 kJ/mol
CH2CHO(34) <=> C2H3O(74) 178.183 kJ/mol
H(30) + C2H2O(67) <=> CH2CHO(34) 363.073 kJ/mol
CH2CHO(34) <=> C2H3O(75) 217.696 kJ/mol
CH2CHO(34) <=> C2H3O(241) 209.758 kJ/mol
O(17) + C2H3(25) <=> CH2CHO(34) 534.417 kJ/mol
C2H3O(74) <=> CO(7) + CH3(10) 50.7615 kJ/mol
H(30) + CH2CO(29) <=> C2H3O(74) 164.552 kJ/mol
H(30) + C2H2O(83) <=> C2H3O(74) 394.026 kJ/mol
C2H3O(74) <=> CH2CHO(34) 150.934 kJ/mol
========================================================================
When I try to run the generated network__.py file in stand-alone measure.py it complains about two different reactions:
ERROR: The input file "examples/rmg/methylformate/pdep/network216_1.py" was invalid:
For path reaction "H(30) + CH2CO(29) <=> C2H3O(74)":
* Both the reactants and the products are product channels.
For path reaction "H(30) + C2H2O(83) <=> C2H3O(74)":
* Both the reactants and the products are product channels.
Unrealistic spinMultiplicity in output.xml
Output.xml seems to contain unrealistic values for spinMultiplicity. Examples include:
<?xml version="1.0" ?>
<rmgoutput>
<speciesList>
<species id="1" label="HXD13" reactive="yes">
<cml>
<molecule spinMultiplicity="105174752">and:
<species id="2" label="CH4" reactive="yes">
<cml>
<molecule spinMultiplicity="-1644429672">f2py cannot build Fortran with g95 on MacOS X
Running the makefile on my Mac OS 10.5 gives me (among other things):
compiling Fortran sources
Fortran f77 compiler: /usr/bin/g95 -ffixed-form -fno-second-underscore -O
Fortran f90 compiler: /usr/bin/g95 -fno-second-underscore -O
Fortran fix compiler: /usr/bin/g95 -ffixed-form -fno-second-underscore -O
compile options: '-I/var/folders/ie/ie3Fq1qC2RaKbE+BYnW87U+++TM/-Tmp-/tmp024OtW/src.macosx-10.3-fat-2.5 -I/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/numpy-1.3.0-py2.5-macosx-10.3-fat.egg/numpy/core/include -I/Library/Frameworks/Python.framework/Versions/2.5/include/python2.5 -c'
g95:f90: spectralfit.f90
/usr/bin/g95 -shared -shared /var/folders/ie/ie3Fq1qC2RaKbE+BYnW87U+++TM/-Tmp-/tmp024OtW/var/folders/ie/ie3Fq1qC2RaKbE+BYnW87U+++TM/-Tmp-/tmp024OtW/src.macosx-10.3-fat-2.5/spectralfitmodule.o /var/folders/ie/ie3Fq1qC2RaKbE+BYnW87U+++TM/-Tmp-/tmp024OtW/var/folders/ie/ie3Fq1qC2RaKbE+BYnW87U+++TM/-Tmp-/tmp024OtW/src.macosx-10.3-fat-2.5/fortranobject.o /var/folders/ie/ie3Fq1qC2RaKbE+BYnW87U+++TM/-Tmp-/tmp024OtW/spectralfit.o calc_freq_code.o dqed.o -lblas -llapack -o ./spectralfit.so
g95: unrecognized option '-shared'
g95: unrecognized option '-shared'
Undefined symbols:
"_PyModule_GetDict", referenced from:
_initspectralfit in spectralfitmodule.o
etc. etc.
I think it may be related to
http://mail.scipy.org/pipermail/numpy-discussion/2007-June/028174.html
http://mail.scipy.org/pipermail/scipy-user/2008-January/015300.html
I'll try switching to Gfortran
'rmgpy.species.Species' has no 'clearLabeledAtoms' in makeNewSpecies
When I run with:
# Data sources
database(
'../RMG-database/input',
thermoLibraries = ['primaryThermoLibrary', 'GRI-Mech3.0'],
reactionLibraries = ['GRI-Mech3.0'],
seedMechanisms = [],
)it fails with:
Adding seed mechanism GRI-Mech3.0 to model core...
Traceback (most recent call last):
File "rmg.py", line 673, in <module>
execute(args)
File "rmg.py", line 262, in execute
reactionModel.addSeedMechanismToCore(database, seedMechanism, react=False)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 1305, in addSeedMechanismToCore
r, isNew = self.makeNewReaction(rxn)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 748, in makeNewReaction
forward.reactants = [self.makeNewSpecies(molecule)[0] for molecule in forward.reactants]
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 641, in makeNewSpecies
molecule.clearLabeledAtoms()
AttributeError: 'rmgpy.species.Species' object has no attribute 'clearLabeledAtoms'Run out of memory allocating stoichiometry matrix
This seems to be the limiting factor on how large a model you can build
Pressure-dependent kinetics in Cantera
Cantera has a number of falloff expressions for representing pressure-dependent kinetics. However, all of these are not much use in unimolecular reaction networks containing multiple wells. A workaround has been written for isothermal, isobaric systems. This will have to suffice until appropriate kinetic models (e.g. Chebyshev polynomials) are available within Cantera.
Negative transition state energies in PDep networks.
When I run RMG-Py on my methyl formate example with MSC pressure-dependence then inspect the network__.py files in the pdep folder, I find quite a few quite negative transition state energies.
Is this a problem? Seems a little worrisome to me, but I don't know the details of how they're estimated or used.
Here is an example:
isomer("C2H5O3")
reactants("methylformate", "OH")
################################################################################
reaction(
reactants=["methylformate(1)", "OH(31)"],
products=["C2H5O3(94)"],
kinetics=Arrhenius(A=(1.45049e+06,"m^3/(mol*s)"), n=0.49959, Ea=(3214.93,"J/mol"), T0=(300,"K")),
transitionState=TransitionState(
E0=(-338.664,"kJ/mol"),
),
)
reaction(
reactants=["CH3(11)", "CH2O3(189)"],
products=["C2H5O3(94)"],
kinetics=Arrhenius(A=(1.02726e+06,"m^3/(mol*s)"), n=0.119613, Ea=(36383.1,"J/mol"), T0=(300,"K")),
transitionState=TransitionState(
E0=(-129.905,"kJ/mol"),
),
)Wilhoit thermo enforces constant Cp for monoatomic species
Things like O atom will (probably?) have incorrect thermo once they have been converted to Wilhoit form (and then onwards to NASA). Wilhoit uses the Cp0 and Cp_inf limits, which will both be 2.5*R for a monoatomic species.
Somebody should check our element atom enthalpy curves. (C, H, O, etc.)
Primary Reaction (as opposed to kinetics) Library
The thing currently called a "Primary Reaction Library" is in fact a primary kinetics library (PKL): once RMG creates a reaction using its templates, we check the PKL to get the rate (if it exists).
What would be useful is actually a Primary Reaction Library (PRL): In the 'generate reactions' step, before running through reaction families and applying templates to the given molecule(s) we would first check them against reactants/products of reactions in the PRL. This way the PRL could contain reactions that are not created by a template. Currently the only way to include such reactions is to add them to the seed mechanism, which forces them to be in the core even if they're irrelevant.
structuretest fails (rmg-py)
I first noticed that Species(SMILES="CCC") was broken, then tracked it down to Structure, then decided to just try the structuretest.py, which fails:
Exception exceptions.TypeError: "'NoneType' object is not iterable" in 'rmg.chem.Atom.isSpecificCaseOf' ignored
Traceback (most recent call last):
File "structuretest.py", line 415, in
StructureCheck('testSubgraphIsomorphismManyLabels').debug()
File "/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/unittest.py", line 286, in debug
getattr(self, self._testMethodName)()
File "structuretest.py", line 124, in testSubgraphIsomorphismManyLabels
""")
File "/Users/rwest/XCodeProjects/RMGpy/source/rmg/structure.py", line 396, in fromAdjacencyList
raise InvalidAdjacencyListException(label + ". Invalid because of "+e.message)
rmg.structure.InvalidAdjacencyListException: Invalid adjacency list: . Invalid because of NULL result without error in PyObject_Call
Species from input file are renamed if already exist in seed mechanism
This could be confusing if I put in my input file a species "METHANE" but my output mechanism doesn't have such a species (I'd have to deduce that 'CH4' was what I wanted, for example).
See 89bbbee
Options:
- add the new name as a comment to the species
- log a warning
- rename the species with the name from the input file
math domain error in rmgpy.kinetics.Chebyshev.fitToData
This may be a problem in the rates that it's trying to fit a Chebyshev to (check the energies in the network - they're all 0kJ/mol!) but the exception occurs in the Chebyshev fitting.
Adding species C(#[C])C#[C](2810) to model core
Generating thermodynamics for new species...
Warning: Average RMS error in heat capacity fit to C(=CC#[C])C#CC#[C](3070) = 0.103363*R
Warning: Average RMS error in heat capacity fit to C(=C=[CH])([CH2])C#CC#[C](3056) = 0.138296*R
Warning: Average RMS error in heat capacity fit to C(=CC#CC#[C])[C]=[CH](3062) = 0.136091*R
Warning: Average RMS error in heat capacity fit to C(=[CH])(C#CC#[C])C#[C](3072) = 0.107614*R
Warning: Average RMS error in heat capacity fit to C(=C)(C#CC#[C])C#[C](3066) = 0.106151*R
For reaction CO2(8) + C(#[C])C#[C](2810) <=> C(#CC#[C])O[C]=O(3016), Ea raised from 2.6 to 56.1 kJ/mol to match endothermicity of reaction.
Updating 1 modified unimolecular reaction networks...
Warning: States node <Entry index=-1 label="CtCr"> and all its parents have data=None
Warning: States node <Entry index=-1 label="CtCr"> and all its parents have data=None
========================================================================
Network Information (Network #3233)
-------------------
Isomers:
C(#[C])C#[C](2810) 0 kJ/mol
Reactant channels:
Product channels:
[C]#[C](102) + [C]#[C](102) 0 kJ/mol
Path reactions:
[C]#[C](102) + [C]#[C](102) <=> C(#[C])C#[C](2810) 0 kJ/mol
========================================================================
Using 200 grains from 0.00 to 563.39 kJ/mol in steps of 2.83 kJ/mol
Calculating densities of states for Network #3233...
Calculating phenomenological rate coefficients for Network #3233...
Traceback (most recent call last):
File "/home/rwest/code/RMG-Py/rmg.py", line 732, in <module>
cProfile.runctx(command, global_vars, local_vars, stats_file)
File "/usr/lib/python2.6/cProfile.py", line 49, in runctx
prof = prof.runctx(statement, globals, locals)
File "/usr/lib/python2.6/cProfile.py", line 140, in runctx
exec cmd in globals, locals
File "<string>", line 1, in <module>
File "/home/rwest/code/RMG-Py/rmg.py", line 359, in execute
reactionModel.enlarge(object)
File "/home/rwest/code/RMG-Py/rmgpy/rmg/model.py", line 1009, in enlarge
self.updateUnimolecularReactionNetworks(database)
File "/home/rwest/code/RMG-Py/rmgpy/rmg/model.py", line 1467, in updateUnimolecularReactionNetworks
network.update(self, database)
File "/home/rwest/code/RMG-Py/rmgpy/rmg/model.py", line 560, in update
netReaction.kinetics = fitInterpolationModel(netReaction, Tlist, Plist, K[:,:,i,j], model, Tmin, Tmax, Pmin, Pmax)
File "reaction.pyx", line 332, in rmgpy.measure.reaction.fitInterpolationModel (build/pyrex/rmgpy/measure/reaction.c:5135)
File "kinetics.py", line 774, in rmgpy.kinetics.Chebyshev.fitToData (build/pyrex/rmgpy/kinetics.c:13455)
File "kinetics.py", line 809, in rmgpy.kinetics.Chebyshev.fitToData (build/pyrex/rmgpy/kinetics.c:13069)
ValueError: math domain errorThis occurs when I change the examples/rmg/methylformate input file to have:
temperature=(1950,'K'),
pressure=(10.0,'bar'),
bond symmetry number unreliable
Commit 0d83783 hid a bug which should be looked at further.
This is the commit message:
The function Structure.split() returns a list of structures that result from
the splitting. In calculateBondSymmetryNumber(), this function was assumed to
return exactly two structures, which it should if things are going normally.
However, this behavior is not guaranteed, so a check was added to ensure that
exactly two structures are returned from split(). If not, the function returns
1, which is as if it was not called for that bond.
Arithmetic underflow in spectral fortran code
_modes.f90:58.44:
qt = ((2 * 3.141592654 * mass) / (6.626e-34 * 6.626e-34))**(dim/2.0) * V
1
Warning: Arithmetic underflow at (1)
this is with gnu fortran gfortran, and may be because the -r8 flag is not recognized?
Chemkin-formatted mechanism
It'd be nice to be able to run the final results in chemkin.
README.rst is inaccurate
This is on the front page of the github repository, and should at least include useful links to where to find setup information, if not the information itself.
Reservoir State PDep gives Negative rate coefficient (and dies)
When I run my methyl formate test case with reservoir state pressure dependence, it dies on the following network:
========================================================================
Network Information (Network #73)
-------------------
Isomers:
HCCOH(31) 84.0013 kJ/mol
Reactant channels:
C2H(3) + OH(32) 593.111 kJ/mol
Product channels:
C2HO(78) + H(30) 620.468 kJ/mol
HCCO(28) + H(30) 356.102 kJ/mol
C2H2O(67) 152.171 kJ/mol
Path reactions:
C2H(3) + OH(32) <=> HCCOH(31) 592.991 kJ/mol
C2HO(78) + H(30) <=> HCCOH(31) 620.003 kJ/mol
HCCO(28) + H(30) <=> HCCOH(31) 355.709 kJ/mol
C2H2O(67) <=> HCCOH(31) 191.919 kJ/mol
========================================================================
Using 200 grains from 84.00 to 1237.87 kJ/mol in steps of 5.80 kJ/mol
Calculating densities of states for Network #73...
Calculating phenomenological rate coefficients for Network #73...
Error: Negative rate coefficient generated; rejecting result.
[[ -1.63573146e-004 -4.14035279e-011]
[ 3.45545096e-096 -3.48219201e+005]
[ 2.13680492e-104 5.07162335e-015]
[ 7.30616588e-041 4.59983034e-002]
[ 1.63573146e-004 3.48219155e+005]]
Error: Negative rate coefficient generated; rejecting result.
[[ -5.07420490e-05 -2.42014475e-11]
[ 1.18795415e-91 -3.53078192e+05]
[ 1.27366817e-99 9.30989438e-15]
[ 2.23810724e-39 4.64161293e-02]
[ 5.07420490e-05 3.53078146e+05]]
Error: Negative rate coefficient generated; rejecting result.
[[ -2.89164833e-03 -1.86363775e-11]
[ 6.86206389e-90 -3.53078192e+05]
[ 7.35717989e-98 9.30989438e-15]
[ 1.29281280e-37 4.64161293e-02]
[ 2.89164833e-03 3.53078146e+05]]
Error: Negative rate coefficient generated; rejecting result.
[[ -1.72916260e-04 -7.55625557e-12]
[ 8.89709669e-82 -3.65793841e+05]
[ 2.57594827e-89 2.79487200e-14]
[ 3.78890937e-35 4.75234804e-02]
[ 1.72916260e-04 3.65793794e+05]]
Error: Negative rate coefficient generated; rejecting result.
[[ -7.29237900e-01 -4.08448180e-11]
[ 6.32541264e-66 -3.94011049e+05]
[ 6.34468135e-73 1.10514814e-13]
[ 1.07233778e-26 4.99441807e-02]
[ 7.29237900e-01 3.94010999e+05]]
Error: Negative rate coefficient generated; rejecting result.
[[ -8.91118317e+00 -4.09657459e-11]
[ 7.73071711e-65 -3.94011049e+05]
[ 7.75426671e-72 1.10514814e-13]
[ 1.31057695e-25 4.99441807e-02]
[ 8.91118317e+00 3.94010999e+05]]
Error: Negative rate coefficient generated; rejecting result.
[[ -8.98194965e+02 -2.32132280e-11]
[ 1.42645769e-49 -4.54384015e+05]
[ 5.04862197e-56 4.29836825e-13]
[ 2.18188589e-18 5.46193148e-02]
[ 8.98194965e+02 4.54383961e+05]]
Traceback (most recent call last):
File "rmg.py", line 732, in <module>
cProfile.runctx(command, global_vars, local_vars, stats_file)
File "/usr/local/Cellar/python/2.7.1/Frameworks/Python.framework/Versions/2.7/lib/python2.7/cProfile.py", line 49, in runctx
prof = prof.runctx(statement, globals, locals)
File "/usr/local/Cellar/python/2.7.1/Frameworks/Python.framework/Versions/2.7/lib/python2.7/cProfile.py", line 140, in runctx
exec cmd in globals, locals
File "<string>", line 1, in <module>
File "rmg.py", line 288, in execute
reactionModel.enlarge(spec, database)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 942, in enlarge
self.updateUnimolecularReactionNetworks(database)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 1403, in updateUnimolecularReactionNetworks
network.update(self, database)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 536, in update
netReaction.kinetics = fitInterpolationModel(netReaction, Tlist, Plist, K[:,:,i,j], model, Tmin, Tmax, Pmin, Pmax)
File "reaction.pyx", line 332, in rmgpy.measure.reaction.fitInterpolationModel (build/pyrex/rmgpy/measure/reaction.c:5141)
File "kinetics.py", line 760, in rmgpy.kinetics.Chebyshev.fitToData (build/pyrex/rmgpy/kinetics.c:13234)
File "kinetics.py", line 795, in rmgpy.kinetics.Chebyshev.fitToData (build/pyrex/rmgpy/kinetics.c:12848)
ValueError: math domain errorgetAllCycles is called far too often
To find a species symmetry number we have to find out if each bond is in a cycle.
This is currently done, for each bond, by finding all the cycles from the atom at one end of the bond.
Perhaps this information can be cached.
Restart doesn't work
Don't think the pickling is working fully
ZeroDivisionError in Arrhenius.changeT0 called from Reaction.generateReverseRateCoefficient
When I add pressure dependence to my (otherwise functioning) methylformate test case, it dies with the following:
========================================================================
Network Information
-------------------
Isomers:
methylformate(1) -368.115 kJ/mol
Reactant channels:
Product channels:
CO(3) + CH4O(4) -330.855 kJ/mol
CHO2(5) + CH3(6) -34.8517 kJ/mol
C2H3O2(7) + H(8) 46.279 kJ/mol
CHO(9) + CH3O(10) 41.2276 kJ/mol
H(8) + C2H3O2(11) 42.8663 kJ/mol
CO2(12) + CH4(13) -487.574 kJ/mol
C2H4O2(14) 27.7775 kJ/mol
Path reactions:
CO(3) + CH4O(4) <=> methylformate(1) -107.597 kJ/mol
CHO2(5) + CH3(6) <=> methylformate(1) -34.4148 kJ/mol
C2H3O2(7) + H(8) <=> methylformate(1) 45.9413 kJ/mol
CHO(9) + CH3O(10) <=> methylformate(1) 39.762 kJ/mol
H(8) + C2H3O2(11) <=> methylformate(1) 42.5453 kJ/mol
CO2(12) + CH4(13) <=> methylformate(1) -156.38 kJ/mol
C2H4O2(14) <=> methylformate(1) 95.9767 kJ/mol
========================================================================
Using 200 grains from -368.12 to 283.99 kJ/mol in steps of 3.28 kJ/mol
Calculating densities of states for Network #1...
Calculating phenomenological rate coefficients for Network #1...
Traceback (most recent call last):
File "rmg.py", line 673, in <module>
execute(args)
File "rmg.py", line 275, in execute
reactionModel.enlarge(spec, database)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 929, in enlarge
self.updateUnimolecularReactionNetworks(database)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 1390, in updateUnimolecularReactionNetworks
network.update(self, database)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 497, in update
K, p0 = self.calculateRateCoefficients(Tlist, Plist, Elist, method)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/measure/network.py", line 465, in calculateRateCoefficients
Kij, Gnj, Fim = self.calculateMicrocanonicalRates(Elist, densStates0, T)
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/measure/network.py", line 361, in calculateMicrocanonicalRates
dummy, Gnj[reac,prod,:] = calculateMicrocanonicalRateCoefficient(rxn, Elist, None, densStates[prod,:], T)
File "reaction.pyx", line 143, in rmgpy.measure.reaction.calculateMicrocanonicalRateCoefficient (build/pyrex/rmgpy/measure/reaction.c:2148)
File "reaction.py", line 322, in rmgpy.reaction.Reaction.generateReverseRateCoefficient (build/pyrex/rmgpy/reaction.c:6740)
File "reaction.py", line 349, in rmgpy.reaction.Reaction.generateReverseRateCoefficient (build/pyrex/rmgpy/reaction.c:6684)
File "kinetics.py", line 320, in rmgpy.kinetics.Arrhenius.changeT0 (build/pyrex/rmgpy/kinetics.c:5276)
ZeroDivisionError: float divisionThe latter part of this stack trace is recreated by one of the unit tests:
======================================================================
ERROR: testGenerateReverseRateCoefficient (reactionTest.TestReaction)
----------------------------------------------------------------------
Traceback (most recent call last):
File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/unittest/reactionTest.py", line 236, in testGenerateReverseRateCoefficient
reverseKinetics = self.reaction2.generateReverseRateCoefficient(Tlist)
File "reaction.py", line 322, in rmgpy.reaction.Reaction.generateReverseRateCoefficient (build/pyrex/rmgpy/reaction.c:6740)
File "reaction.py", line 349, in rmgpy.reaction.Reaction.generateReverseRateCoefficient (build/pyrex/rmgpy/reaction.c:6684)
File "kinetics.py", line 320, in rmgpy.kinetics.Arrhenius.changeT0 (build/pyrex/rmgpy/kinetics.c:5276)
ZeroDivisionError: float divisionReaction path degeneracy not properly treated
The kinetics in the database are given on a per-reaction site basis because the reaction path degeneracy varies with the species involved. To correct for this, we need to multiply the rate from the database by the reaction path degeneracy. But how do we determine this degeneracy?
It may be possible to do this at the reaction generation step, as RMG will automatically generate all of the degenerate reactions. However, we must be careful to not generate the reactions again in a later step. For example, the reaction A -> B can be generated either by reaction of A or reaction of B.
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.