Allow RMG to handle a specific species as a third body collider #1084

alongd · 2017-07-12T19:51:20Z

Previously RMG could not handle reactions of the sort A + B (+C) = D (+C). Moreover, reactions like H2 + Ar <=> H + H + Ar were only handled correctly if the kinetics is of Arrhenius type (like in the BurkeH2O2 libraries), but if the kinetics is in any PDep form, RMG would consider it as H2 + Ar (+M) <=> H + H + Ar (+M), wrongly accounting for an additional third body.

Usage
Simple add a specific third body collider to your library reaction as you would intuitively (i.e., no need for a SpecificCollider keyword).
For example (and this rxn is the motivation for this PR):

entry(
    index = 226,
    label = "SO2 + O (+N2) <=> SO3 (+N2)",
    degeneracy = 1,
    kinetics = Troe(
        arrheniusHigh = Arrhenius(A=(3.7e+11, 'cm^6/(mol^2*s)'), n=0, Ea=(1689, 'cal/mol'), T0=(1, 'K')),
        arrheniusLow = Arrhenius(A=(2.9e+27, 'cm^9/(mol^3*s)'), n=-3.58, Ea=(5206, 'cal/mol'), T0=(1, 'K')),
        alpha=0.43, T3=(371, 'K'), T1=(7442, 'K'), efficiencies={}),
    shortDesc = u"""[Marshall2005],[Marshall2006]""",
)

Another example w/o parenthesis, which RMG should also handle correctly now, is:

entry(
    index = 18,
    label = "HNCN + N2 <=> H + NCN +N2",
    degeneracy = 1,
    kinetics = ThirdBody(
        arrheniusLow = Arrhenius(A=(1.79e+28, 'cm^3/(mol*s)'), n=-3.44, Ea=(64502, 'cal/mol'), T0 = (1, 'K'), Tmin=(2000, 'K'), Tmax=(4000, 'K'))),
    shortDesc = u"""[Lin2000a]""",
)

Currently we only allow one specific collider per reaction (As far as I can see, I don't think more than one would be necessary, and probably also isn't implemented in Chemkin/Cantera, although I haven't checked). E.g., A + B (+N2) (+Ar) <=> C (+N2) (+Ar) isn't valid in RMG, and will throw an appropriate error.

The implementation of this feature in the solver (simple.pyx) is through the Peff calculation.

Two tests were added:

testReadSpecificCollider in chemkinTest.py
testSpecificColliderModel in simpleTest.py

Lastly, this PR also allows RMG to consider specificColliders species with parenthesis in their name such as N2(4) or CH2(S) or CH2(S)(15), as was recently fixed in Cantera and will be added to the next Chemkin release.

Addresses #1070

alongd · 2017-07-12T20:04:25Z

This specificCollider PR currently fails on the new testSpecificColliderModel test.
I'm still looking into it, but would love another pair of eyes helping me figure out whether an additional change to the solver should be made 😯 💭

rwest · 2017-07-12T20:14:07Z

I don't have time to look properly at this this week, but I did have a question about your second example above:

entry(
    index = 18,
    label = "HNCN + N2 <=> H + NCN +N2",
    degeneracy = 1,
    kinetics = ThirdBody(
        arrheniusLow = Arrhenius(A=(1.79e+28, 'cm^3/(mol*s)'), n=-3.44, Ea=(64502, 'cal/mol'), T0 = (1, 'K'), Tmin=(2000, 'K'), Tmax=(4000, 'K'))),
    shortDesc = u"""[Lin2000a]""",
)

What is this meant to represent?

If it's just r˲ = k[HNCN][N₂] then wouldn't
Arrhenius(A=(1.79e+28, 'cm^3/(mol*s)'), n=-3.44, Ea=(64502, 'cal/mol'))
fully describe it, rather than
ThirdBody(arrheniusLow=Arrhenius(A=(1.79e+28, 'cm^3/(mol*s)'), n=-3.44, Ea=(64502, 'cal/mol'))) ?

Perhaps you could say how it would look in CHEMKIN syntax? Or how would you describe the kinetics?

alongd · 2017-07-13T18:27:15Z

I agree, there's no advantage in writing a reaction with a specificCollider in a ThirdBody format, i.e.:

entry(
    index = 1,
    label = "H + O2 + CO2 <=> O + OH + CO2",
    degeneracy = 1,
    kinetics = ThirdBody(
        arrheniusLow = Arrhenius(A=(1.79e+28, 'cm^6/(mol^2*s)'), n=-3.44, Ea=(64502, 'cal/mol'), T0 = (1, 'K'), Tmin=(2000, 'K'), Tmax=(4000, 'K'))),
)

and it's indeed much simpler (and recommended?) to write it as a normal Arrhenius, as pointed out by @rwest.

However, if a user happens to use the ThirdBody kinetics format for a library reaction with a specific collider, we would like RMG to handle that correctly. Before this fix, the reaction in the example above would have the following wrong Chemkin representation (with two colliders):

! Reaction index: Chemkin #2; RMG #2
! Library reaction: FFCM1(-)
! Flux pairs: CO2(6), CO2(6); O2(1), O(7); H(3), OH(4); 
H(3)+O2(1)+CO2(6)+M=OH(4)+O(7)+CO2(6)+M             1.790e+34 -3.440    64.502

But now RMG will pick that up, add a comment that a specificCollider was identified, and output the kinetics as if a normal Arrhenius kinetics reaction was inputted:

! Reaction index: Chemkin #2; RMG #2
! Library reaction: FFCM1(-)
! Specific third body collider: CO2
! Flux pairs: O2(1), O(6); H(3), OH(4); 
H(3)+O2(1)+CO2(7)=OH(4)+O(6)+CO2(7)                 1.790e+28 -3.440    64.502

alongd · 2017-07-14T02:12:54Z

Added kinetics documentation. Also addresses RMG-database#191

alongd · 2017-07-14T02:20:15Z

Couldn't figure out how to add relative references to the Classes in modification.rst pointing outside the users folder in the documentation (so added absolute).
Any ideas?

rwest · 2017-07-14T03:52:54Z

If the sole aim is to protect against user error, I would suggest stopping with a helpful error message, rather than trying to guess what they intended. If they made a mistake, then we may guess incorrectly.

(Also, rather unlikely, but if someone actually did want to input "H + O2 + CO2 +M <=> O + OH + CO2 +M" or "A + CO2 + M <=> B + CO2 + M" there would be no way to do so in the new format?)

alongd · 2017-07-14T12:15:38Z

@rwest: Will do. Thanks for that.
Let's see that I got it right: a reaction with a third body collider in parenthesis, whether specific or not, should only be applied to kinetics formats that specify High and Low rates with different corresponding orders (units).
In other words, we should only allow the parenthesis form (+M) or (+specificSpecies) in Troe and Lindemann, right?

alongd · 2017-07-14T13:09:24Z

I see that Chebyshev also gets a (+M), meaning a specificCollider should be allowed there as well?

mjohnson541

Look mostly good! Just a few small details.

mjohnson541 · 2017-07-14T17:53:17Z

rmgpy/reaction.py

        """
        arrow = ' <=> '
        if not self.reversible: arrow = ' => '
-        return arrow.join([' + '.join([str(s) for s in self.reactants]), ' + '.join([str(s) for s in self.products])])
+        if self.specificCollider:
+            return arrow.join([' + '.join([str(s) for s in self.reactants])+' '+str(self.specificCollider), ' + '.join([str(s) for s in self.products])+' '+str(self.specificCollider)])


should the parentheses around the specificCollider be added here?

mjohnson541 · 2017-07-14T18:30:39Z

rmgpy/data/kinetics/library.py

+            collider = re.search('\(\+[^\)]+\)',reactants)
+            if collider is not None:
+                collider = collider.group(0) # save string value rather than the object
+                assert collider == re.search('\(\+[^\)]+\)',products).group(0), "Third body colliders in reaction {0} in kinetics library {1} are not identical!".format(rxn_string, self.label)


Would it make sense to changes this to "missing from products or are not identical"?

mjohnson541 · 2017-07-14T18:46:42Z

rmgpy/data/kinetics/family.py

@@ -29,7 +29,7 @@
 ################################################################################

 """
-This module contains functionality for working with kinetics libraries.
+This module contains functionality for working with kinetics families.


Can you squish this one with the one before?

mjohnson541 · 2017-07-14T19:23:10Z

rmgpy/solver/simple.pyx

+
+    """
+    specificColliderSpecies:
+    a list that contains Species classes of species which are specific third body colliders


I know this makes sense in general, but it reads kind of strangely to me can you change "Species classes of" to "object references to"?

mjohnson541 · 2017-07-14T19:44:41Z

rmgpy/solver/simple.pyx

+                    j = self.reactionIndex[rxn]
+                    pdepSpecificColliderReactionIndices.append(j)
+                    self.pdepSpecificColliderKinetics.append(rxn.kinetics)
+                    for spc in itertools.chain(coreSpecies, edgeSpecies):


Why can't you just use:
self.specificColliderSpecies.append(rxn.specificCollider)
instead of this loop structure? Is there something missing in the rxn.specificCollider object? Can you add a comment about this?

mjohnson541 · 2017-07-14T19:55:58Z

rmgpy/reaction.py

-        else:
-            logging.debug('did not modify A factor when modifying degeneracy since ' \
-                          'the reaction rate was not set')
+        #else:


If we're not going to let it print the error message I think we might as well delete these lines and clean up the code.

mjohnson541 · 2017-07-14T20:01:19Z

documentation/source/users/rmg/database/modification.rst

+
+**MultiArrhenius**
+(see `MultiArrhenius Class <http://reactionmechanismgenerator.github.io/RMG-Py/reference/kinetics/multiarrhenius.html>`_ for details)::
+


Can you specify that MultiArrhenius is the sum of multiple Arrhenius expressions (as opposed to being one Arrhenius over one Trange and one Arrhenius over another Trange)

mjohnson541 · 2017-07-14T20:06:47Z

documentation/source/users/rmg/database/modification.rst

+	        kunits='cm^3/(mol*s)', Tmin=(300, 'K'), Tmax=(3000, 'K'), Pmin=(0.0013156, 'atm'), Pmax=(131.56, 'atm')))
+
+
+The ``Troe``, ``Lindemann``, and ``Chebyshev`` pressure dependent formats could also be defined with a specific collider if needed, for example::


Can you add a title so this doesn't get passed over as part of the Chebyshev section?

rwest · 2017-07-14T20:24:41Z

I don't know how chemkin would interpret (+N2) (meaning a specificCollider) on a CHEB or a PLOG, both of which are defined in terms of pressure (Atm) rather than species concentrations. Perhaps it would figure out the partial pressure of N2. I have never seen such an expression "in the wild".

Now that we parse things like (+N2) in our RMG reaction definitions, I wonder if we should parse (ad require) (+M) to signify pressure-dependent reactions, and +M to signify third-body reactions (low pressure limit), like Chemkin and Cantera. Consequences would have to be thought through, but it might reduce confusion for people translating between RMG and Chemkin or Cantera formats.

At the moment the chemkin reactions A (+M) = B + C (+M) and A = B + C and A +M = B + C +M are all written as A = B + C in RMG (but with different kinetics attributes).

Perhaps not for this pull request, but might be worth thinking about.

alongd · 2017-07-14T21:54:18Z

While Chemkin's pre-processor failed on a PLOG reaction with a specific collider:

  1. CH3(17)+CH2NH2(91)(+Ar)=NCC(1)(+Ar)            1.00E+00    0.0        0.0
      Error...PLOG and other pressure-depency mutually exclusive...
      Rate coefficients at P=4.93E-01(atm)           9.381E+54  -12.3       22.7
      Rate coefficients at P=7.72E-01(atm)           4.176E+52  -11.6       21.8
      Rate coefficients at P=1.21E+00(atm)           1.011E+50  -10.8       20.6
      Rate coefficients at P=1.89E+00(atm)           1.393E+47   -9.9       19.2
      Rate coefficients at P=2.96E+00(atm)           1.186E+44   -9.0       17.7
      Error...no pres-dependency option defined...

it did allow it in a CHEB reaction (also, RMG outputs CHEB with (+M), but not PLOG).
I'm not sure why CHEB gets a different treatment in this respect (but at least it's consistent).

The current state of this branch is to allow a specific collider only for reactions that get a (+M): Troe, Lindemann, and Chebyshev.

The Chebyshev case is indeed strange. Generally, if a specific collider is present, RMG calculates the rate at the partial pressure of the collider species alone. For Chebyshev this could be a problem if it isn't defined for this specific (low) pressure.

I tried running a CHEB format with a specific collider again in Chemkin, this time using a minor species (not Ar). Chemkin failed right away since the intermediate species' concentration was 0 in the beginning of the simulation.

I think we shouldn't allow specific colliders in Chebyshev, although using one of the bath gases as a specific collider could make some sense.

alongd · 2017-07-14T22:06:18Z

Thanks @mjohnson541!
I addressed all comments, the fixes for library.py and family.py were squashed.
We now also throw an error if a specificCollider is given for Chebyshev, following @rwest's and my comments above (unless there are any objections?)
Let me know if everything's OK, and I'll organize and squash.

alongd · 2017-07-14T22:06:49Z

BTW, what's up with Travis lately?

specificCollider is a Species, unlike reactants and products which are lists of Species

A specificCollider species could appear in one of the folowing forms: 'A + B (+S) = C + D (+S)`, or `A + B + S = C + D + S`, depending on the kenetics type of the PDep reaction. RMG first finds a (+S) type specificCollider, then also looks for a species that appears both in the reactants and products (regardless of order). If found, verify that this is the only scpecificCollider (as we currently only support one specificCollider per reaction), and remove it from the reactants and products lists. It is saved as a specificCollider attribute of the Reaction.

Generalizes reading the third body collider whether it is (+M), (+m) or a specific species, i.e., (+N2). Gives RMG the ability to read the collider even if the species name has paranthesis, e.g., N2(5) or CH2(T)(15). This is still not implemented in Chemkin and Cantera, but both have open issues to implement it in the future. Also adds a comment to verbose Chemkin file that a specific third body collider was identified.

Test that a Chemkin reaction with a specific species as a third body collider can be properly read, even if the species name contains parenthesis such as (+CH(T))

- family.py: Added the attribute, though currently not fully tested in families. - pdep.py - LibraryReaction in database.py - depository.py

- Tests that a specific collider is correctly handeled by the solver - Also tests that RMG can reac a specific collider with parenthesis in its name

Other than adding - pdepSpecificColliderKinetics - specificColliderSpecies - pdepSpecificColliderReactionIndices The main change is in the Peff calculation

Due to the numerous calls to __setDegeneracy(), this debug message appears repeatedly while debuging (or running tests with -d). In fact, it could get printed so many times that it fills the CMD stack.

- a detailed list of all available libraries and their description - examples of all libraries kinetic formats - a note regarding the implementation of a specificCollider

alongd mentioned this pull request Jul 12, 2017

Implementing kinetics for a specific third body collider #1070

Closed

alongd force-pushed the collider branch from 6386237 to 37f97d6 Compare July 12, 2017 19:56

mjohnson541 self-requested a review July 12, 2017 21:13

alongd force-pushed the collider branch 2 times, most recently from a4585f7 to 61e84bb Compare July 13, 2017 20:33

alongd added Type: Feature Topic: PDep Status: Ready for Review PR is complete and ready to be reviewed labels Jul 13, 2017

alongd force-pushed the collider branch from d1bee42 to f54e687 Compare July 14, 2017 02:15

alongd force-pushed the collider branch from a72a594 to 36b8cac Compare July 14, 2017 14:16

mjohnson541 requested changes Jul 14, 2017

View reviewed changes

alongd force-pushed the collider branch from 36b8cac to a8633b4 Compare July 14, 2017 21:57

mjohnson541 approved these changes Jul 15, 2017

View reviewed changes

alongd force-pushed the collider branch from d775b8e to ee625bb Compare July 17, 2017 13:09

alongd added Status: Accepted This PR is approved and ready to be merged and removed Status: Ready for Review PR is complete and ready to be reviewed labels Jul 17, 2017

alongd added 2 commits July 17, 2017 16:35

Added a specificCollider attribute to the Reaction class

aa8e238

specificCollider is a Species, unlike reactants and products which are lists of Species

alongd added 8 commits July 17, 2017 16:35

Test that a specificCollider can be read

4e0fb0a

Test that a Chemkin reaction with a specific species as a third body collider can be properly read, even if the species name contains parenthesis such as (+CH(T))

Handeling specificCollider attribute

491cbb4

- family.py: Added the attribute, though currently not fully tested in families. - pdep.py - LibraryReaction in database.py - depository.py

Enable model.py to handel specificCollider as a model Species

e100a7e

Added a testSpecificColliderModel

ea26608

- Tests that a specific collider is correctly handeled by the solver - Also tests that RMG can reac a specific collider with parenthesis in its name

Modified simple.pyx solver to account for specificCollider

3e4f2bf

Other than adding - pdepSpecificColliderKinetics - specificColliderSpecies - pdepSpecificColliderReactionIndices The main change is in the Peff calculation

Removed logging.debug message in reaction.py __setDegeneracy()

f7f1fee

Due to the numerous calls to __setDegeneracy(), this debug message appears repeatedly while debuging (or running tests with -d). In fact, it could get printed so many times that it fills the CMD stack.

Added kinetics documentation

619c2cf

- a detailed list of all available libraries and their description - examples of all libraries kinetic formats - a note regarding the implementation of a specificCollider

alongd force-pushed the collider branch from ee625bb to 619c2cf Compare July 17, 2017 20:35

alongd merged commit 8fbf885 into master Jul 17, 2017

mliu49 deleted the collider branch August 1, 2017 19:34

alongd mentioned this pull request Dec 7, 2017

Add specificCollider to areIdenticalSpeciesReferences() #1246

Merged

rwest mentioned this pull request Feb 14, 2022

ck2yaml crashes when parsing a third body reaction with PLOG Cantera/cantera#1165

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Allow RMG to handle a specific species as a third body collider #1084

Allow RMG to handle a specific species as a third body collider #1084

alongd commented Jul 12, 2017 •

edited

Loading

alongd commented Jul 12, 2017 •

edited

Loading

rwest commented Jul 12, 2017

alongd commented Jul 13, 2017 •

edited

Loading

alongd commented Jul 14, 2017

alongd commented Jul 14, 2017

rwest commented Jul 14, 2017

alongd commented Jul 14, 2017

alongd commented Jul 14, 2017

mjohnson541 left a comment

mjohnson541 Jul 14, 2017

mjohnson541 Jul 14, 2017

mjohnson541 Jul 14, 2017

mjohnson541 Jul 14, 2017

mjohnson541 Jul 14, 2017

mjohnson541 Jul 14, 2017

mjohnson541 Jul 14, 2017

mjohnson541 Jul 14, 2017

rwest commented Jul 14, 2017

alongd commented Jul 14, 2017

alongd commented Jul 14, 2017 •

edited

Loading

alongd commented Jul 14, 2017


		MultiArrhenius
		(see `MultiArrhenius Class <http://reactionmechanismgenerator.github.io/RMG-Py/reference/kinetics/multiarrhenius.html>`_ for details)::

		kunits='cm^3/(mol*s)', Tmin=(300, 'K'), Tmax=(3000, 'K'), Pmin=(0.0013156, 'atm'), Pmax=(131.56, 'atm')))


		The ``Troe``, ``Lindemann``, and ``Chebyshev`` pressure dependent formats could also be defined with a specific collider if needed, for example::

Allow RMG to handle a specific species as a third body collider #1084

Allow RMG to handle a specific species as a third body collider #1084

Conversation

alongd commented Jul 12, 2017 • edited Loading

alongd commented Jul 12, 2017 • edited Loading

rwest commented Jul 12, 2017

alongd commented Jul 13, 2017 • edited Loading

alongd commented Jul 14, 2017

alongd commented Jul 14, 2017

rwest commented Jul 14, 2017

alongd commented Jul 14, 2017

alongd commented Jul 14, 2017

mjohnson541 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rwest commented Jul 14, 2017

alongd commented Jul 14, 2017

alongd commented Jul 14, 2017 • edited Loading

alongd commented Jul 14, 2017

alongd commented Jul 12, 2017 •

edited

Loading

alongd commented Jul 12, 2017 •

edited

Loading

alongd commented Jul 13, 2017 •

edited

Loading

alongd commented Jul 14, 2017 •

edited

Loading