Invalid k(E) values computed for path reaction. #77
Comments
|
I am still getting this after a trial merge of the 'molecule' branch. Merging PDepNetwork #41809 and PDepNetwork #41810
Updating 2 modified unimolecular reaction networks...
========================================================================
Network Information (Network #41863)
-------------------
Isomers:
CC=C=[CH](13976) 315.344 kJ/mol
C(C)C#[C](5134) 486.735 kJ/mol
Reactant channels:
CH3(10) + C3H2(69) 628.787 kJ/mol
H(17) + C(=C)C#C(197) 498.607 kJ/mol
CH3(10) + [CH2]C#[C](2305) 681.976 kJ/mol
C2(40) + C2H5(27) 943.022 kJ/mol
Product channels:
H(17) + C(=C=[CH])[CH2](15712) 664.527 kJ/mol
H(17) + C[C]=C=[CH](13984) 750.731 kJ/mol
H(17) + C[CH]C#[C](13980) 720.442 kJ/mol
[C@@H]1(C)C=[C]1(31851) 478.812 kJ/mol
C([CH2])C#C(531) 354.849 kJ/mol
C([CH2])C#[C](511) + H(17) 903.775 kJ/mol
Path reactions:
CH3(10) + C3H2(69) <=> CC=C=[CH](13976) 628.787 kJ/mol
H(17) + C(=C=[CH])[CH2](15712) <=> CC=C=[CH](13976) 664.527 kJ/mol
H(17) + C[C]=C=[CH](13984) <=> CC=C=[CH](13976) 750.731 kJ/mol
H(17) + C(=C)C#C(197) <=> CC=C=[CH](13976) 502.791 kJ/mol
H(17) + C[CH]C#[C](13980) <=> CC=C=[CH](13976) 720.442 kJ/mol
CC=C=[CH](13976) <=> [C@@H]1(C)C=[C]1(31851) 541.28 kJ/mol
C([CH2])C#C(531) <=> CC=C=[CH](13976) 471.583 kJ/mol
C(C)C#[C](5134) <=> CC=C=[CH](13976) 631.545 kJ/mol
CH3(10) + [CH2]C#[C](2305) <=> C(C)C#[C](5134) 681.976 kJ/mol
C2(40) + C2H5(27) <=> C(C)C#[C](5134) 943.022 kJ/mol
H(17) + C[CH]C#[C](13980) <=> C(C)C#[C](5134) 720.442 kJ/mol
C([CH2])C#[C](511) + H(17) <=> C(C)C#[C](5134) 903.775 kJ/mol
C(C)C#[C](5134) <=> C([CH2])C#C(531) 631.545 kJ/mol
========================================================================
Calculating densities of states for Network #41863...
Using 463 grains from 315.34 to 1998.47 kJ/mol in steps of 3.64 kJ/mol to compute densities of states
Calculating phenomenological rate coefficients for Network #41863...
Using 200 grains from 315.34 to 1040.33 kJ/mol in steps of 3.64 kJ/mol to compute the k(T,P) values at 292.578 K
Using 200 grains from 315.34 to 1047.55 kJ/mol in steps of 3.68 kJ/mol to compute the k(T,P) values at 314.289 K
Using 200 grains from 315.34 to 1064.16 kJ/mol in steps of 3.76 kJ/mol to compute the k(T,P) values at 364.231 K
Using 200 grains from 315.34 to 1095.90 kJ/mol in steps of 3.92 kJ/mol to compute the k(T,P) values at 459.667 K
Using 203 grains from 315.34 to 1160.51 kJ/mol in steps of 4.18 kJ/mol to compute the k(T,P) values at 641.642 K
Using 232 grains from 315.34 to 1281.85 kJ/mol in steps of 4.18 kJ/mol to compute the k(T,P) values at 1011.65 K
Using 295 grains from 315.34 to 1545.44 kJ/mol in steps of 4.18 kJ/mol to compute the k(T,P) values at 1810.91 K
Using 403 grains from 315.34 to 1997.31 kJ/mol in steps of 4.18 kJ/mol to compute the k(T,P) values at 3163.57 K
========================================================================
Network Information (Network #41809)
-------------------
Isomers:
H2CCCH(68) 347.693 kJ/mol
CC#[C](459) 509.242 kJ/mol
C1C=[C]1(4505) 513.113 kJ/mol
Reactant channels:
C3H2(69) + H(17) 703.529 kJ/mol
H(17) + [CH2]C#[C](2305) 756.717 kJ/mol
CH3(10) + C2(40) 968.493 kJ/mol
Product channels:
H(17) + C1=C=C1(19322) 706.706 kJ/mol
H(17) + C1=[C][CH]1(19323) 826.518 kJ/mol
H(17) + C1[C]=[C]1(19324) 988.462 kJ/mol
[CH]1C=C1(19325) 377.133 kJ/mol
Path reactions:
H2CCCH(68) <=> C3H2(69) + H(17) 675.915 kJ/mol
H(17) + [CH2]C#[C](2305) <=> H2CCCH(68) 756.717 kJ/mol
H2CCCH(68) <=> C1C=[C]1(4505) 573.629 kJ/mol
CC#[C](459) <=> H2CCCH(68) 654.052 kJ/mol
CH3(10) + C2(40) <=> CC#[C](459) 968.493 kJ/mol
H(17) + [CH2]C#[C](2305) <=> CC#[C](459) 756.717 kJ/mol
H(17) + C1=C=C1(19322) <=> C1C=[C]1(4505) 718.839 kJ/mol
H(17) + C1=[C][CH]1(19323) <=> C1C=[C]1(4505) 826.518 kJ/mol
H(17) + C1[C]=[C]1(19324) <=> C1C=[C]1(4505) 988.462 kJ/mol
C1C=[C]1(4505) <=> [CH]1C=C1(19325) 657.923 kJ/mol
========================================================================
Calculating densities of states for Network #41809...
Using 458 grains from 347.69 to 2042.67 kJ/mol in steps of 3.71 kJ/mol to compute densities of states
Calculating phenomenological rate coefficients for Network #41809...
Using 200 grains from 347.69 to 1085.77 kJ/mol in steps of 3.71 kJ/mol to compute the k(T,P) values at 292.578 K
Error: For path reaction H2CCCH(68) <=> C3H2(69) + H(17):
Error: Expected kf(292.578 K) = 1.31499e-46
Error: Actual kf(292.578 K) = 7.44885e-49
Error: Expected Keq(292.578 K) = 2.85228e-57
Error: Actual Keq(292.578 K) = 2.85228e-57
Traceback (most recent call last):
File "rmg.py", line 134, in <module>
cProfile.runctx(command, global_vars, local_vars, stats_file)
File "/usr/local/Cellar/python/2.7.3/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.3/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 "/Users/rwest/Code/rmgpy/RMG-Py/rmgpy/rmg/main.py", line 391, in execute
self.reactionModel.enlarge(objectsToEnlarge)
File "/Users/rwest/Code/rmgpy/RMG-Py/rmgpy/rmg/model.py", line 615, in enlarge
self.updateUnimolecularReactionNetworks(database)
File "/Users/rwest/Code/rmgpy/RMG-Py/rmgpy/rmg/model.py", line 1291, in updateUnimolecularReactionNetworks
network.update(self, database, self.pressureDependence)
File "/Users/rwest/Code/rmgpy/RMG-Py/rmgpy/rmg/pdep.py", line 443, in update
K = self.calculateRateCoefficients(Tlist, Plist, method, grainSize=grainSize, grainCount=grainCount)
File "/Users/rwest/Code/rmgpy/RMG-Py/rmgpy/measure/network.py", line 852, in calculateRateCoefficients
self.setConditions(T, P)
File "/Users/rwest/Code/rmgpy/RMG-Py/rmgpy/measure/network.py", line 761, in setConditions
self.calculateMicrocanonicalRates()
File "/Users/rwest/Code/rmgpy/RMG-Py/rmgpy/measure/network.py", line 613, in calculateMicrocanonicalRates
raise NetworkError('Invalid k(E) values computed for path reaction "{0}" do not satisfy'.format(rxn))
rmgpy.measure.network.NetworkError: Invalid k(E) values computed for path reaction "H2CCCH(68) <=> C3H2(69) + H(17)" do not satisfy(But with the faster pdep it now only takes 6 hours to reach that point) |
So the kf values are an order of magnitude apart, but those are very small numbers (whatever the units may be). What is the machine precision? Are they really different? |
|
Yes, I still see this more often than I would like. Some of the error is probably just due to discretization, but I don't think that explains a full order of magnitude. This is especially true when using ILT to go from k(T) -> k(E), which by definition should give you the same k(T) at the end. (This assumes that Ea and the barrier height are the same, which is what we're usually assuming in RMG jobs.) Feel free to send me a few networks that are known to have the issue. I've never really thought about machine precision being an issue. Now that you mention it, that could be problematic, since the formula is k(E) = A * rho(E - Ea) / rho(E), where rho(E) is the reactant density of states. For large Ea - as in this case - the ratio rho(E - Ea) / rho(E) is dividing a small number by a big number, which could be troublesome. I'll think about this some more; it may be that we don't care about this issue when the kf is negligible. Looking at the first example in the thread, I see that the barrier height is negative. Not sure if RMG-Py handles this correctly, especially because I'm not really sure how it should be handled. |
|
The first example in the thread is also very small kf: 1.31499e-46. species(
label='Hyad',
reactive=True,
structure=SMILES("C(C=O)O"),
)but I'll send you an input file. |
|
Just realized that although in all these cases the kf is very small (which still hints at numerical error accumulation?) a very small Keq means that the reverse rate could in fact be significant and important to get right? Implementing an "ignore this problem if both kf and kr are less then 1e30" allows the Hyad example to get past the first problem, and run straight into another one: But what is "small" for a rate coefficient anyway? Are these even always in the same units (eg. all bimolecular?) |
|
The 'sean-v8@26fbe21' still not working, I got the same error: I changed cutoff to e-15, but again job crashed for the same reason: |
|
Yes, even allowing relatively large discrepancies only seems to delay the inevitable |
…aceTransformMethod. Hopefully this is the nail in the coffin of issue #77 (The numerics seem to be worst for low n)
|
In commit 22f4dd3 I made it so that if it gets an invalid k(E) error, it tries again with twice as many grains. This seemed to help a lot, but running the methylformate example, it runs for three hours, then this happens: Those attempts take about another 3 hours, then it crashes. We need a strategy for when adding more grains doesn't help. |
|
Clearly increasing the number of grains is not helping in this system, as the actual kf barely changes at all. This suggests a real (not a numerical) error. If you can send me the network I'll see if I can identify what's going on. |
|
The problem with this particular network is that the reaction H2CCCH(65)<=>C3H2(60)+H(17), coming from Glarborg/C2, has an activation energy that is ~20 kcal/mol below the barrier height (as estimated from the enthalpy of reaction). Using the ILT method to get k(E) requires that the activation energy and the barrier height be the same. Since the k(T) estimate came from a library instead of an estimate, its Ea was not adjusted to be at least dHrxn. I'm not quite sure what the appropriate course of action is. |
|
One thing we can do is ensure that increasing the number of grains is getting us closer to an acceptable k. What to do with such an exception, or how to avoid it in the first place, however, is still unclear. |
|
At the meeting Bill suggested that we go ahead and raise the barrier height to match the enthalpy of reaction for kinetics from reaction libraries, just like we currently do for rate rules. This would at least allow things to continue. Otherwise you would just have to remove/comment out the offending reaction from the library. Even though the plot above is shifted due to the systematic error, the rate of convergence looks rather typical of what I've seen. Clearly adding additional grains is having an effect, even if the convergence is to the wrong value in this instance. On the flip side, there is clearly convergence, so we can probably stop trying smaller discretizations before 10^6 grains. Incidentally, when you get a network that doesn't work you should be able to run CanTherm on it and reproduce the same error, but also get a drawing of the PES to study. This would perhaps more readily help identify issues such as this reaction. That way you might not have to wait for me to get around to it. |
|
We just ran into this again (exactly the same reaction) |
|
Do we want to raise the barrier just for the sake of PDep, or raise it for all reactions always (eg. as soon as we look it up)? |
|
Barrier height problem hopefully addressed in f477d79 |
|
Looks like f477d79 didn't solve it. I've pared down an input file that hangs on its first network, to help debugging: # Data sources
database(
thermoLibraries = ['primaryThermoLibrary','DFT_QCI_thermo','GRI-Mech3.0'],
reactionLibraries = [('Methylformate',False),('Glarborg/highP',False)],
seedMechanisms = [],
kineticsDepositories = ['training'],
kineticsFamilies = ['!Intra_Disproportionation'],
kineticsEstimator = 'rate rules',
)
# List of species
species(
label='H2CCCH',
reactive=True,
structure=SMILES("C#C[CH2]"),
)
species(
label='[CH]=C=[CH]',
reactive=True,
structure=SMILES("[CH]=C=[CH]"),
)
species(
label='[H]',
reactive=True,
structure=SMILES("[H]"),
)
# Bath gas
species(
label='Ar',
reactive=False,
structure=InChI("InChI=1S/Ar"),
)
# Reaction systems
simpleReactor(
temperature=(650,'K'),
pressure=(1.0,'bar'),
initialMoleFractions={
"H2CCCH": 0.01,
"Ar": 0.08,
},
terminationTime=(0.5,'s'),
)
simpleReactor(
temperature=(1350,'K'),
pressure=(3.0,'bar'),
initialMoleFractions={
"H2CCCH": 0.01,
"Ar": 0.97,
},
terminationTime=(0.5,'s'),
)
simpleReactor(
temperature=(1950,'K'),
pressure=(10.0,'bar'),
initialMoleFractions={
"H2CCCH": 0.01,
"Ar": 0.97,
},
terminationTime=(0.5,'s'),
)
simulator(
atol=1e-22,
rtol=1e-8,
)
model(
toleranceKeepInEdge=0.0,
toleranceMoveToCore=0.0005,
toleranceInterruptSimulation=1.0,
maximumEdgeSpecies=100000
)
pressureDependence(
method='modified strong collision', # 'reservoir state'
maximumGrainSize=(1.0,'kcal/mol'),
minimumNumberOfGrains=200,
temperatures=(290,3500,'K',8),
pressures=(0.02,100,'bar',5),
interpolation=('Chebyshev', 6, 4),
)
options(
units='si',
saveRestartPeriod=None,
drawMolecules=True,
generatePlots=False,
saveConcentrationProfiles=False,
)Gives this network: |
|
In this case it looks like the reaction at issue is given in the |
|
I thought the Ea >= dHrxn was always checked for endothermic Arrhenius reactions. |
|
What should we do with |
|
If that reaction is not going to be treated as pressure-dependent, then we can leave it as-is. For pressure dependence I think we'll have to fit a single Arrhenius so that we can have a single barrier height. So I guess do that and then adjust the resulting Ea if necessary? |
|
If it's in a Reaction Library, it'll be forced through PDep. Do we turn everything into Arrhenius at that point? (I guess we could pre-emptively do it when loading the reaction library). I see around line 473 of rmgpy/rmg/pdep.py elif isinstance(rxn.kinetics, KineticsData):
...
elif not isinstance(rxn.kinetics, Arrhenius):
raise Exception('Path reaction "{0}" in PDepNetwork #{1:d} has invalid kinetics type "{2!s}".'.format(rxn, self.index, rxn.kinetics.__class__))So it looks like everything is expected to be KineticsData or Arrhenius by that point, but I'm not sure where it'd be converted if it started as MultiArrhenius. Do seed mechanisms avoid the PDep code entirely? (i.e. should this checking not be done for them?) |
|
Commit 8aa385d checks for convergence inside the "add more grains" loop. |
|
I've found what happens to MultiArrhenius reactions in a reaction library... They die: File "/nics/d/home/rwest/Code/RMG-Py/rmgpy/rmg/main.py", line 428, in execute
self.reactionModel.enlarge(objectsToEnlarge)
File "/nics/d/home/rwest/Code/RMG-Py/rmgpy/rmg/model.py", line 655, in enlarge
self.updateUnimolecularReactionNetworks(database)
File "/nics/d/home/rwest/Code/RMG-Py/rmgpy/rmg/model.py", line 1389, in updateUnimolecularReactionNetworks
network.update(self, database, self.pressureDependence)
File "/nics/d/home/rwest/Code/RMG-Py/rmgpy/rmg/pdep.py", line 490, in update
raise Exception('Path reaction "{0}" in PDepNetwork #{1:d} has invalid kinetics type "{2!s}".'.format(rxn, self.index, rxn.kinetics.__class__))
Exception: Path reaction "OH(11) + C2H2(41) <=> CHCHOH(43)" in PDepNetwork #384 has invalid kinetics type "<type 'rmgpy.kinetics.arrhenius.MultiArrhenius'>". |
|
I've made a .toArrhenius() method for MultiArrhenius kinetics. When do you think this should be called, @jwallen ? The place I had put the call to fixBarrierHeight() when adding reaction libraries to the edge does not seem to catch any MultiArrhenius kinetics. |
|
Think I have fixed this. Decided to implement it within the PDepNetwork code. |
…code. Previously, kinetics coming from a reaction library (such as Glarborg/HighP) were exempt from the fixBarrierHeight check. Now, if we are going to run them through a PDep calculation (PDep is on and the reaction is unimolecular) then we fix the barrier height. Should address #77, or at least the comment #77 (comment)
Still addressing issue #77. Commit f477d79 added a lot of what's needed, but the fixBarrierHeight was not being called on Seed Mechanisms or Reaction Libraries, because these avoid the model.enlarge() method. This commit adds that, calculating the thermo when necessary, and also limits it to only working on Arrhenius reactions, which have an Ea to change. Remaining potential problem: Some high-P libraries have things like MultiArrhenius reactions.
…calculation. This should stop the endless loop that fills memory with 100 million grains and fails, when a network is fundamentally flawed. (See issue #77)
pierrelb
added
the
Status: Stale
|
This issue is being automatically marked as stale because it has not received any interaction in the last 90 days. Please leave a comment if this is still a relevant issue, otherwise it will automatically be closed in 30 days. |
Not sure what's going on here, but the methylformate example job dies with this:
Measure input file for network 41809 is at https://gist.github.com/2910382
The text was updated successfully, but these errors were encountered: