Modified Strong Collision PDep gives negative steady-state concentration (and dies) #44

rwest · 2011-06-18T00:12:29Z

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.

The text was updated successfully, but these errors were encountered:

rwest · 2011-07-19T15:53:30Z

Roll forward a month (and some changes) and this is still my problem:

========================================================================
Network Information (Network #145)
-------------------
Isomers:
    CO[CH]O[C]=O(170)                                    -133.858 kJ/mol
Reactant channels:
    Mfmt(1) + CO(7)                                      -491.686 kJ/mol
Product channels:
    CH3(10) + C(=O)O[C]=O(184)                             -143.6 kJ/mol
    CO2(8) + CO[CH](419)                                 -138.391 kJ/mol
    CH3(10) + [CH](O[C]=O)[O](420)                        204.839 kJ/mol
    H(30) + [CH](O[CH2])O[C]=O(421)                       264.895 kJ/mol
    C(O[CH2])O[C]=O(422)                                 -151.111 kJ/mol
    CO(7) + CO[CH][O](423)                                -74.267 kJ/mol
Path reactions:
    Mfmt(1) + CO(7) <=> CO[CH]O[C]=O(170)                -462.821 kJ/mol
    CH3(10) + C(=O)O[C]=O(184) <=> CO[CH]O[C]=O(170)     -107.217 kJ/mol
    CO2(8) + CO[CH](419) <=> CO[CH]O[C]=O(170)           -109.526 kJ/mol
    CH3(10) + [CH](O[C]=O)[O](420) <=> CO[CH]O[C]=O(170)      205.413 kJ/mol
    H(30) + [CH](O[CH2])O[C]=O(421) <=> CO[CH]O[C]=O(170)      264.557 kJ/mol
    CO[CH]O[C]=O(170) <=> C(O[CH2])O[C]=O(422)            12.5529 kJ/mol
    CO(7) + CO[CH][O](423) <=> CO[CH]O[C]=O(170)          -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 #145...
Calculating density of states for isomer "CO[CH]O[C]=O(170)"
Calculating density of states for reactant channel "Mfmt(1) + CO(7)"

Using ILT method to compute k(E) for path reaction Mfmt(1) + CO(7) <=> CO[CH]O[C]=O(170).
Using ILT method to compute k(E) for path reaction CH3(10) + C(=O)O[C]=O(184) <=> CO[CH]O[C]=O(170).
Using ILT method to compute k(E) for path reaction CO2(8) + CO[CH](419) <=> CO[CH]O[C]=O(170).
Using ILT method to compute k(E) for path reaction CH3(10) + [CH](O[C]=O)[O](420) <=> CO[CH]O[C]=O(170).
Using ILT method to compute k(E) for path reaction H(30) + [CH](O[CH2])O[C]=O(421) <=> CO[CH]O[C]=O(170).
Using ILT method to compute k(E) for path reaction CO[CH]O[C]=O(170) <=> C(O[CH2])O[C]=O(422).
Using ILT method to compute k(E) for path reaction CO(7) + CO[CH][O](423) <=> CO[CH]O[C]=O(170).

Calculating phenomenological rate coefficients for Network #145...
Warning: invalid value encountered in divide
Applying modified strong collision method at 292.578 K, 0.0246349 bar...
Traceback (most recent call last):
  File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmg.py", line 740, in <module>
    execute(args)
  File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmg.py", line 359, in execute
    reactionModel.enlarge(object, database)
  File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 951, in enlarge
    self.updateUnimolecularReactionNetworks(database)
  File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 1411, 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.

This network can be found at http://rmg.mit.edu/measure/networks/e0052901b39abb259bf6b3120cfcff1a

rwest · 2011-07-19T17:30:02Z

I think this commit 72097a2 which converted all estimated KineticsData type reaction rates into a Arrhenius rates, allowing the barriers to be raised to the endothermicity of the reaction (by d35614d ), has [temporarily] made this problem go away for me. I am not entirely sure that the PES for this network would now be estimated in such a way that MEASURE works, but I think my current job does not (yet) try to explore this network.

rwest · 2011-07-28T20:08:28Z

I just got this again:

Calculating densities of states for Network #382...
Calculating phenomenological rate coefficients for Network #382...

========================================================================
Network Information (Network #387)
-------------------
Isomers:
    CH3CH2OH(42)                                         -248.823 kJ/mol
Reactant channels:
    CH3(10) + CH2OH(28)                                   109.466 kJ/mol
Product channels:
    OH(19) + C2H5(27)                                     142.111 kJ/mol
    H(17) + CH3CHOH(43)                                    142.79 kJ/mol
    H(17) + CH2CH2OH(44)                                  171.513 kJ/mol
    H(17) + CH3CH2O(31)                                   185.637 kJ/mol
    H2O(5) + C2H4(12)                                    -207.679 kJ/mol
Path reactions:
    CH3(10) + CH2OH(28) <=> CH3CH2OH(42)                  462.831 kJ/mol
    OH(19) + C2H5(27) <=> CH3CH2OH(42)                    142.111 kJ/mol
    H(17) + CH3CHOH(43) <=> CH3CH2OH(42)                  142.402 kJ/mol
    H(17) + CH2CH2OH(44) <=> CH3CH2OH(42)                 171.513 kJ/mol
    H(17) + CH3CH2O(31) <=> CH3CH2OH(42)                  185.236 kJ/mol
    H2O(5) + C2H4(12) <=> CH3CH2OH(42)                    69.1933 kJ/mol
========================================================================

Using 200 grains from -248.82 to 1570.21 kJ/mol in steps of 9.14 kJ/mol
Calculating densities of states for Network #387...
Calculating phenomenological rate coefficients for Network #387...
Warning: invalid value encountered in divide
---------------------------------------------------------------------------
ModifiedStrongCollisionError              Traceback (most recent call last)

/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmg.py in <module>()
    140     else:
    141         rmg = RMG()
--> 142         rmg.execute(args)
    143 
    144 

/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/main.pyc in execute(self, args)
    254             for spec in self.initialSpecies:
    255                 if spec.reactive:
--> 256                     self.reactionModel.enlarge(spec)
    257 
    258             # Save a restart file if desired


/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py in enlarge(self, newObject)
    595 
    596             # Recalculate k(T,P) values for modified networks

--> 597             self.updateUnimolecularReactionNetworks(database)
    598             logging.info('')
    599 

/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py in updateUnimolecularReactionNetworks(self, database)
   1058         # self = reactionModel object

   1059         for network in self.unirxnNetworks:
-> 1060             network.update(self, database, self.pressureDependence)
   1061 
   1062     def loadSeedMechanism(self, path):

/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/pdep.pyc in update(self, reactionModel, database, pdepSettings)
    411 
    412         # Calculate the rate coefficients

--> 413         K, p0 = self.calculateRateCoefficients(Tlist, Plist, Elist, method)
    414 
    415         # Generate PDepReaction objects


/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/measure/network.pyc in calculateRateCoefficients(self, Tlist, Plist, Elist, method, errorCheck)
    589                     # Apply modified strong collision method

    590                     import msc
--> 591                     K[t,p,:,:], p0[t,p,:,:,:] = msc.applyModifiedStrongCollisionMethod(T, P, Elist, densStates, collFreq, Kij, Fim, Gnj, Ereac, Nisom, Nreac, Nprod)
    592                 elif method.lower() == 'reservoir state':
    593                     # Apply reservoir state method


/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/measure/msc.so in rmgpy.measure.msc.applyModifiedStrongCollisionMethod (build/pyrex/rmgpy/measure/msc.c:2023)()
    137 
    138 
--> 139 
    140 
    141 

ModifiedStrongCollisionError: A negative steady-state concentration was encountered.
> /Users/rwest/XCodeProjects/RMGpy/RMG-Py/msc.pyx(139)rmgpy.measure.msc.applyModifiedStrongCollisionMethod (build/pyrex/rmgpy/measure/msc.c:2023)()

The network is at http://rmg.mit.edu/measure/networks/3be78458ad307d3ee4ca98e2e971f9d2

rwest · 2011-07-29T00:17:43Z

When I fix the maximumGrainSize feature and set it to 1 kcal/mol, the above network solves succesfully:

========================================================================
Network Information (Network #387)
-------------------
Isomers:
    CH3CH2OH(42)                                         -248.823 kJ/mol
Reactant channels:
    CH3(10) + CH2OH(28)                                   109.466 kJ/mol
Product channels:
    OH(19) + C2H5(27)                                     142.111 kJ/mol
    H(17) + CH3CHOH(43)                                    142.79 kJ/mol
    H(17) + CH2CH2OH(44)                                  171.513 kJ/mol
    H(17) + CH3CH2O(31)                                   185.637 kJ/mol
    H2O(5) + C2H4(12)                                    -207.679 kJ/mol
Path reactions:
    CH3(10) + CH2OH(28) <=> CH3CH2OH(42)                  462.831 kJ/mol
    OH(19) + C2H5(27) <=> CH3CH2OH(42)                    142.111 kJ/mol
    H(17) + CH3CHOH(43) <=> CH3CH2OH(42)                  142.402 kJ/mol
    H(17) + CH2CH2OH(44) <=> CH3CH2OH(42)                 171.513 kJ/mol
    H(17) + CH3CH2O(31) <=> CH3CH2OH(42)                  185.236 kJ/mol
    H2O(5) + C2H4(12) <=> CH3CH2OH(42)                    69.1933 kJ/mol
========================================================================

Using 436 grains from -248.82 to 1571.22 kJ/mol in steps of 4.18 kJ/mol
Calculating densities of states for Network #387...
Calculating phenomenological rate coefficients for Network #387...
etc...

Should we catch that type of ModifiedStrongCollisionError and try again with smaller grains?

jwallen · 2011-07-29T14:26:38Z

The solution for this particular network is to fix the estimate of the barrier height for the CH3(10) + CH2OH(28) <=> CH3CH2OH(42) reaction, which is way too high. (This is also the reason the PES drawing is cut off, as I wasn't accounting for barrier heights when determining the drawing dimensions.)

Interestingly, this network runs fine for me as a standalone MEASURE job (after fixing a small but unrelated bug in the saved input file).

… result. When getting the kinetics for a reaction generated from a template, the kinetics could conceivably come in either direction if a match for the reaction is found within a depository (e.g. training), as we don't force these to be stored in the direction the template is defined in. Therefore we need to mark the returned kinetics as being for the forward or reverse direction, so they can be properly interpreted. All kinetics estimated via group additivity or from the rules are by definition for the forward direction, as before. This fixes, among other things, the absurdly-high barrier for CH3 + CH2OH -> CH3CH2OH seen in the PES for issue #44.

faribas · 2013-06-11T20:32:28Z

RMG still sometimes produces networks that it can't solve, such as this one:

which again leads to

Traceback (most recent call last):
  File "cantherm.py", line 61, in <module>
    cantherm.execute()
  File "/Users/Fariba/Code/RMG-Py/rmgpy/cantherm/main.py", line 235, in execute
    job.execute(outputFile=outputFile, plot=self.plot)
  File "/Users/Fariba/Code/RMG-Py/rmgpy/cantherm/pdep.py", line 231, in execute
    self.K = self.network.calculateRateCoefficients(self.Tlist.value_si, self.Plist.value_si, self.method)
  File "/Users/Fariba/Code/RMG-Py/rmgpy/pdep/network.py", line 214, in calculateRateCoefficients
    self.applyModifiedStrongCollisionMethod()
  File "/Users/Fariba/Code/RMG-Py/rmgpy/pdep/network.py", line 747, in applyModifiedStrongCollisionMethod
    self.K, self.p0 = msc.applyModifiedStrongCollisionMethod(self, efficiencyModel)
  File "msc.pyx", line 52, in rmgpy.pdep.msc.applyModifiedStrongCollisionMethod (build/pyrex/rmgpy/pdep/msc.c:4064)
  File "msc.pyx", line 150, in rmgpy.pdep.msc.applyModifiedStrongCollisionMethod (build/pyrex/rmgpy/pdep/msc.c:3432)
rmgpy.pdep.msc.ModifiedStrongCollisionError: A negative steady-state concentration was encountered.

(It's only a problem at the lowest temperature (293K). If we start it at 314K it solves OK.)

Whats the best way to prevent RMG form proposing networks it can't solve, or making it more resilient?

shamelmerchant · 2013-06-11T20:49:50Z

Have you tried using the -x (or max-energy option) that will eliminate the very high energy species.

faribas · 2013-06-11T20:52:24Z

No. Option to CanTherm or RMG? Is it documented anywhere? Thanks!

github-actions · 2023-06-22T08:10:54Z

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.

ghost assigned jwallen Jul 28, 2011

rwest unassigned jwallen May 26, 2015

mjohnson541 mentioned this issue Dec 19, 2017

Max-energy option #1254

Closed

alongd added Topic: PDep Arkane labels Feb 8, 2018

github-actions bot added the stale stale issue/PR as determined by actions bot label Jun 22, 2023

github-actions bot added the abandoned abandoned issue/PR as determined by actions bot label Jul 22, 2023

github-actions bot closed this as not planned Won't fix, can't repro, duplicate, stale Jul 22, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Modified Strong Collision PDep gives negative steady-state concentration (and dies) #44

Modified Strong Collision PDep gives negative steady-state concentration (and dies) #44

rwest commented Jun 18, 2011

rwest commented Jul 19, 2011

rwest commented Jul 19, 2011

rwest commented Jul 28, 2011

rwest commented Jul 29, 2011

jwallen commented Jul 29, 2011

faribas commented Jun 11, 2013

shamelmerchant commented Jun 11, 2013

faribas commented Jun 11, 2013

github-actions bot commented Jun 22, 2023

Modified Strong Collision PDep gives negative steady-state concentration (and dies) #44

Modified Strong Collision PDep gives negative steady-state concentration (and dies) #44

Comments

rwest commented Jun 18, 2011

rwest commented Jul 19, 2011

rwest commented Jul 19, 2011

rwest commented Jul 28, 2011

rwest commented Jul 29, 2011

jwallen commented Jul 29, 2011

faribas commented Jun 11, 2013

shamelmerchant commented Jun 11, 2013

faribas commented Jun 11, 2013

github-actions bot commented Jun 22, 2023