Some PDep reactions are generated as irreversible

[alongd]

Some PDep reactions in many RMG models are irreversible (See e.g., issue #731). From an email exchange with @jwallen, this is indeed unintentional. As a first stage or short-term solution, @whgreen suggested to add a module that runs subsequent to model completion that by default converts all => signs to = (unless the user explicitly chooses not to in the input file). This will also affect many ERC-FoundationFuelv0.9 library reactions which are defined as irreversible, e.g.:

CH3(10)+O(26)=>H(5)+H2(7)+CO(29)
O(26)+CH2(34)=>H(5)+H(5)+CO(29)
O2(4)+CH2(34)=>O(26)+CH2O(35)

Since we don't want to allow the reactions with three products to become reversible, perhaps the suggested module should prohibit this?

[rwest]

A post-processing fix won't help flux analysis (and thus model expansion decisions) while building the model, and if it's going to involve special cases for seed mechanisms and reaction libraries, it may end up being as difficult to write as fixing the problem at source. So I'd recommend spending a few more minutes looking into fixing the problem where it occurs before starting to implement a post-processing fix.

[whgreen]

Alon, Are reversible reactions with three products really a problem? Can you test in CHEMKIN and in the RMG ode solver? Bill

[alongd]

I agree. @whgreen just suggested that we should build a tool which kicks in subsequent to each PDep network, changes all => sings to =, and checks that there are no identical pre-existing reactions in the model (in which case, we should remove them). Other than that, we should look into the networks that generated the irreversible reactions, follow how they were explored, and see if we can tease out why the final reaction is still irreversible.

[whgreen]

Alon, I think we should do this fixing on irreversible library reactions as well, not just PDep. -Bill

[alongd]

So I think we should generalize the fix to every reaction added, whether PDep or library (where the user can still choose not to change reversibility of library reactions, but by default they'll be forced to be =). This way we won't ignore flux analysis.

[whgreen]

Yes. But we need a check to delete the irreversible reverse reactions if they are present Bill

Sent from my iPhone

On Aug 3, 2016, at 4:17 PM, Alon Grinberg Dana notifications@github.com<mailto:notifications@github.com> wrote:

So I think we should generalize the fix to every reaction added, whether PDep or library (where the user can still choose not to change reversibility of library reactions, but by default they'll be forced to be =). This way we won't ignore flux analysis.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/ReactionMechanismGenerator/RMG-Py/issues/736#issuecomment-237358934, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ABK_sEIVdIVcyGSDCbiM1Dw7bIWutAs8ks5qcPdRgaJpZM4JbMB2.

[mliu49]

The code to convert all pdep reactions to reversible and remove the irreversible reverse reaction already exists in RMG. The issue seems to be that it's not being called consistently.

In rmgpy.rmg.model.updateUnimolecularReactionNetworks(), RMG will

1. Merge any networks which have the same source and at least one isomer in common
2. Update invalid networks
3. Convert all net reactions in updated networks to reversible (and remove duplicate reverse rxns)

From what I can find, merging networks and adding path reactions will mark a network as invalid. The change to only iterate over updated networks happened with 0cd7425fee513c4ce4d671ce04f60e5ced55270d, which seemed to be speedup related and not expected to affect functionality.

Although it seems like every network will be updated and therefore corrected for reversibility upon any changes, this clearly isn't happening. A preliminary test shows that if step 3 is modified to act on all existing networks, then no irreversible reactions will make it into the final mechanism, with no other differences in the mechanism.

My questions for people who are more familiar with pdep (hopefully @jwallen can chime in):

• Do any code changes come to mind which have invalidated Josh's original assumption that every network would be updated following a change? If the commit referenced above introduced this issue, it seems incredible that it's been unnoticed for over two years.
• Could we be predicting incorrect rates if we're not updating some networks which should be updated? If so, then simply correcting reversibility for all networks will not be enough to fully address the issue.

[rwest]

Helpful sleuthing, @mliu49. Regarding the "if this commit introduced the issue..." question: it could be that it was working fine until some time later another change altered which networks get marked as "invalid".

[whgreen]

I think Alon found that we have published several mechanisms with irreversible Pdep reactions. They could be 2 years old.

Sent from my iPhone

[rwest]

Should rmgpy.pdep.PDepNetwork.updateConfigurations invalidate the network? I think it currently does not. It is called from rmgpy.rmg.model.CoreEdgeReactionModel.enlarge.

[mliu49]

Sorry, I realized that I made a mistake while running the test jobs. Iterating over all networks to change reversibility of net reactions does affect the model generated. For the same input file, fewer species and reactions were present in the core, and a few of the rates were different between the two models. This might be because removal of the reverse reactions affects the flux calculations.

I'm not sure if updateConfigurations should invalidate the network. It is guaranteed to be called on every network for every enlarge call with a new object (ie. reactEdge=false). If it did invalidate the network, then all of the rate constants would be updated with every iteration. It would successfully make every reaction reversible, at the cost of much (probably unneeded) computation.

[rwest]

You're right, updateConfigurations is a poor place to invalidate the network. I tried adding .invalidate() calls in a few other places that looked sensible, but it didn't help.

I now think the problem is that a reaction needs to be made reversible once it enters the core (i.e. all its reactants and products are in the core), but this may not happen on the same iteration that the network is modified in such a way as to make it marked as invalid. eg. with the 1,3-hexadiene example, HXD13 enters the core first (then CH4 and H2 from the input file), then the HXD13 dissociates into a bunch of things in PDep Network 1, which is explored. The Pdep reaction rate HXD13 (+M) => CH3 +C5H7 (+M) is evaluated and at this point it's irreversible. The path to C5H7 and CH3 is found to be fastest and is followed. First C5H7 enters the core, then in the next step, CH3. At this point, adding CH3 to the core, the PDep Network 1 is not updated in any way - not expanded or merged - and so is not marked as invalid for recalculation. However, all the products are now in the core, so it should be marked as reversible at this stage. i.e. we can't only iterate over "updated" (recently invalid) networks when looking for core reactions to make reversible. Iterating over all the networks, as you suggest, seems to fix the problem. Pull request #746 shows some debugging and a possible solution.

[alongd]

My two cents: perhaps we should add a check to see if a reaction that is added to the core (i.e., it's reactants/products were all added) is reversible, and if not then invalidate it and make reversible.

[mliu49]

I think that's exactly the problem. I'm also fairly sure now that commit 0cd7425fee513c4ce4d671ce04f60e5ced55270d introduced the issue. The old behavior was to loop through all of the core reactions, and convert any irreversible pdep reactions to reversible. The behavior post-commit was changed to loop through the net reactions of every pdep network, then check if the reaction in in the core before converting to reversible.

As @rwest said, if the reaction has not been added to the core yet, it will not (and likely never) be marked as reversible. Possible solutions:

1. Loop over all networks, which has been shown to work (method in #746)
2. Loop over all core reactions (essentially reverting 0cd7425fee513c4ce4d671ce04f60e5ced55270d)
3. Perform the conversion and search for reverse reactions upon adding a pdep reaction to the core

I'm not sure about the relative computational cost of each method. I think Method 1 involves more iterations than Method 2, since edge reactions are also included. Method 3 is similar to Method 2, but should require one less for loop.

[rwest]

I think methods 2 or 3. Probably 2 is simpler to implement and keeps things in places we'd expect to look for them, eg. keeps more of the PDep code together. But probably 3 has some benefits too.

[mliu49]

I tested method 2. It runs marginally faster than method 1, mostly due to not needing to search the list of core reactions to retrieve the reaction index. The results from method 1 and 2 are identical.

I attempted method 3, but realized that addReactionToCore is sometimes called before kinetics are generated for the reaction, so the check for thermodynamic consistency doesn't work. Moving to a point where kinetics are guaranteed to have been generated basically brings us to the current location.

I think we should go with method 2. Should I push to a new branch or to debug736?

[whgreen]

Max and Richard,

Please help me understand how this method 2 algorithm works. Is it like this:

1) Grow some P-dep networks (and non P-dep reactions) on the edge. Since they are on the edge, many of these reactions are irreversible.

2) Choose a species X to move to the core.

3) Move the non P-dep reactions involving X as a product from edge to core (iff all of the species in the reaction are now in the core)

4) Add reactions including X as a reactant. In the P-dep network that was (A+B = C, C-->Q+R, C-->X, C-->S) where A,B,C,D are in the core and Q,R,S are on the edge, this might add the elementary steps X --> Y + Z and X --> C. So now this P-dep network has been updated to (A+B = C, C-->Q+R, C=X, C-->S, X --> Y+Z).

5) I think some new P-dep networks starting from X or X + something will be initiated e.g. (X-->C, X--> Y+Z) maybe also (X+A --> D)

6) Compute the k(T,P)’s for the updated or new P-dep networks (i.e. the P-dep networks that originally involved species X as a unimolecular product), and update the core. For the network I wrote out above, the new or updated k(T,P) would include these core reactions:

 A+B = C   (reaction already in core, but now the rate k(T,P) is updated)

 A+B --> X  (new)

in addition there will be another P-dep network on the edge, which originally said C-->X but now is (C=X, X--> Y+Z); from that we should get a reversible reaction

C=X with a k(T,P).

There should also be a new P-dep network

 X--> C  somehow RMG should recognize that we already covered this one in the reversible reaction C=X above.

There should also be a new P-dep reaction X+A --> D (but there should also be D--> X+A on the edge, and RMG should recognize these are reverse reactions).

7) Polish the core by adding reverse reactions if necessary (e.g. change A+B --> X to A+B = X)

8) Integrate the ODEs implied by the core reaction mechanism

Is this correct? Please correct my misunderstandings. I am particularly interested in how RMG updates the k(T,P) for reactions in the core when their P-dep network expands.

-Bill

From: Max Liu [mailto:notifications@github.com] Sent: Wednesday, August 17, 2016 3:36 PM To: ReactionMechanismGenerator/RMG-Py Cc: William H Green; Mention Subject: Re: [ReactionMechanismGenerator/RMG-Py] Some PDep reactions are generated as irreversible (#736)

I think that's exactly the problem. I'm also fairly sure now that commit 0cd7425https://github.com/ReactionMechanismGenerator/RMG-Py/commit/0cd7425fee513c4ce4d671ce04f60e5ced55270d introduced the issue. The old behavior was to loop through all of the core reactions, and convert any irreversible pdep reactions to reversible. The behavior post-commit was changed to loop through the net reactions of every pdep network, then check if the reaction in in the core before converting to reversible.

As @rwesthttps://github.com/rwest said, if the reaction has not been added to the core yet, it will not (and likely never) be marked as reversible. Possible solutions:

1. Loop over all networks, which has been shown to work
2. Loop over all core reactions
3. Perform the conversion and search for reverse reactions upon adding a pdep reaction to the core

I'm not sure about the relative computational cost of each method. I think Method 1 involves more iterations than Method 2, since edge reactions are also included, but the benefit is that edge reactions are also marked as reversible for users who are interested in the edge. Method 3 is similar to Method 2, but should require one less for loop.

Do we care about irreversible reactions in the edge?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/ReactionMechanismGenerator/RMG-Py/issues/736#issuecomment-240522448, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ABK_sE4vH-wSZnXL3g6iGZmobmb-WnnJks5qg2KIgaJpZM4JbMB2.

[mliu49]

Here is my understanding:

1. RMG generates edge reactions (both pdep and non-pdep) by reacting all core species together. The pdep reactions will all be irreversible. New pdep networks will be created at this stage if needed.
2. RMG finds a species which exceeds the flux tolerance (X), which will be added to the core.
   • All edge reactions (both pdep and non-pdep) whose reactants and products are all in the core get moved to the core. (D --> X + A should be added to the core at this point, if it was already in the edge)
3. RMG alternatively finds a pdep network (A+B = C, C --> Q+R, C --> X, C --> S) which exceeds the flux tolerance, in which case it identifies the maximum leak species. (X)
   • The species is explored, meaning that all possible unimolecular reactions are generated.
     (X --> Y+Z and X --> C)
   • This is mutually exclusive with step 2 during each iteration. The path taken depends on whether a species or a pdep network was identified as reaching the flux tolerance during the simulation.
4. RMG will then go through all pdep networks to check if any isomers are now in the core, and it will explore those isomers as well. (if X entered the core, then X --> Y+Z and X --> C would be discovered here instead)
5. After obtaining new reactions from exploring isomers, RMG adds them to the pdep network the isomer was from. This step marks the network as invalid.
   Important: RMG will check if the reaction is present in either direction, so X --> C will not be added to the network, but X --> Y+Z will be added.
   The network is now (A+B = C, C --> Q+R, C --> X, C --> S, X --> Y+Z).
6. High pressure limit kinetics are generated for all new reactions. This might include both X --> Y+Z and X --> C, even though X --> C ended up not being used. I'm still a bit unclear on this point.
7. RMG updates pdep networks as needed
   • First step is to merge any networks which share a source and at least one isomer. This invalidates the network.
   • Second step is to update the rate constants, and also update the net reactions. (A+B --> X should be added here)
   • Last step is the part we've been discussing, which currently goes through the net reactions of every pdep network, checks if they're in the core, marks them as reversible, and searches for and removes the reverse reaction if it exists. Method 2 from above would change this step to search through the core reactions directly for irreversible pdep reactions.
8. If X was added to the core, then RMG will generate edge reactions again by reacting the core species. At this stage, the new reaction X+A --> D might be discovered, and a new pdep network would be created for it. However, I'm not entirely sure about this either. I don't know if there is an earlier check for whether the reaction already exists, or if it will persist until the next iteration when it is removed during step 7.
9. RMG simulates the core again, then returns to step 2 or 3 depending on the first thing to reach the flux tolerance.

I think this is accurate, but there are a lot of nuances in the code which I may have misinterpreted and many parts which I don't understand yet. I also have not looked much at the details of the actual k(T,P) calculation. Let me know if this made sense, and if you have additional questions.

[mliu49]

I created a new branch and pull request for method 2. I can also copy some of the commits from #746 if we think that some of the additional invalidate calls and asserts are necessary.

[rwest]

FYI for those following this issue, discussion continues on #752, which tries to fix it.

[whgreen]

Richard, Nice catch! -Bill