RMG estimates kinetics non-deterministically

[KEHANG]

The order of molecules in spc.molecule can eventually lead to same reactions but different templates matched and therefore different kinetics. This difference could be exposed as modelling diverging if newly developed features have a different molecule order than master branch.

I wrote a post using a detailed example at here. I want to take it as a chance to discuss what's the possible ways to solve this bug.

One way to make more deterministic might be redefine the equality of reaction; that is treat reactions with different templates as different reaction although all the species in the reaction are same isomorphically.

[connie]

This may be relevant: #323. This uses the idea of tracking multiple transition states in a single family. The difficulty lies in generating the reverse direction template.

I think we already allow duplicate reactions when they are present in separate reaction families.

[rwest]

We could:

1. Take the fastest rate (like we take the lowest energy isomer for thermo). But fastest at what temperature, if they have different activation energies?
2. Include all the reactions as duplicates (may lead to excessive rate estimates)
3. Include all the reactions, but weight them by the Boltzmann distribution based on the energy of the reacting isomer, so the reaction of the most stable resonance isomer might be scaled by 0.80 and the rate of reaction of the less stable resonance isomer by 0.20. This gets rather complicated, especially when you have multiple isomers on both reactants.

any other ideas?

[KEHANG]

I'm thinking in processNewReactions() to decide whether rxn1 == rxn2, we should not only consider whether reactants and products from rxn1 are same as those from rxn2, but also whether template1 from rxn1 is same as template2 from rxn2. We could still keep only rxn1 and eliminate rxn2 but make rxn1.template = [template1, template2].

[KEHANG]

And when estimating kinetics, rxn1's rate should be the sum of rate of template1 and rate of template2

[nickvandewiele]

@KEHANG your ipython notebook is very valuable as a unit test. Could you convert it into one, and add the WIP decorator onto it so that it passes the travis CI tests.

[KEHANG]

Ok, I'll make one.

[KEHANG]

@nickvandewiele this is the test I made: 5de3826

[KEHANG]

I made this PR, the test inside is decorated by @work_in_progress. Once we got a solution to the issue, we can turn the decorator off. I will merge the test into master if no more further suggestions.

[connie]

@rwest I would probably support solution 2 at least in the case of a single reacting isomer and multiple transition states. However, in the case where there there are multiple resonance forms.... it becomes more complex. Solution 3 is interesting but might be a little too hard to implement. I'm still unsure about Solution 1 and will get back to you on that.

[rwest]

We could do Solution 1, but choose "fastest at all temperatures" by (a) summing them, which is almost close if they are sufficiently different, or (b) evaluate them all at a bunch of T, find the highest at each T, and re-fit a new modified Arrhenius through these.

[goldmanm]

Thanks for the test case @KEHANG. I'm working on a fix for dealing with the kinetics caused by different resonance structures now. If anyone else is already working on a solution, let me know.

@rwest. when refitting, I want to avoid losing any data about the templates (and degeneracies of each template) used in the estimation. Max and I have three ideas of how to implement this while keeping the template and degeneracy information. Which one sounds the best to you?

1. Modify all reactions to use multiarrhenius, and store templates and degeneracies in lists.
2. Keep the default of reaction to be Arrhenius with the original types for template and degeneracy, but switch to multiarrhenius with degeneracies & templates in lists when the reaction hits multiple rate templates.
3. Do not modify the kinetics, template or degeneracy datatypes in reaction objects. Instead, switch a duplicate flag to True (similar to some of our library reactions). This will be printed out in chemkin and cantera outputs.

Any thoughts would be helpful.

[rwest]

Our MultiArrhenius class corresponds to a DUPLICATE reaction in Chemkin, so I don't think I understand option 3. Currently what Chemkin and Cantera treat as multiple reactions (marked as duplicate) A+B<=>C+D, RMG thinks of as one reaction, with MultiArrhenius kinetics.

Is proposal to allow TemplateReaction objects to store lists of templates and degeneracies instead of only a single one of each? (and families?)

[goldmanm]

I guess the idea behind item 3 was to have separate reaction objects for each non-degenerate transition state path, and not modify the types held by the TemplateReaction attributes template, degeneracy, (and possibly family too.). I was thinking something like the duplicate flag could indicate presence of multiple of these reaction objects (similar to how it's done in chemkin/cantera). The benefit of this is it might be more intuitive to understand and easier to integrate in, but also might take more memory and processing time.

If option 3 is thrown out, the proposal would be as you suggested: store lists of templates and degeneracies. This would likely require also changing Reaction, since the degeneracy parameter is held there. The difference between option 1 and 2 is whether the default type is kept (and lists are used only when multiple templates are found) or whether all reactions are made to hold lists for these parameters. Any thoughts on this?

Based on your idea, I also think making the family parameter a list would be beneficial (even though we try to design families to not overlap). It would have helped limit the nondeterministic kinetics caused when the R addition endo/exo cyclic `families overlapped a while ago.

[goldmanm]

closed with #1055