Handling duplicate PLOGs correctly

[shamelmerchant]

I think there are few issues right now how RMG-Py handles them

1. Chemkin allows duplicate plogs having missing pressures we have strict check for it right now.

2. Sum of PLOGS - papers which have sum of Arrhenius equations must be input as
   PLOG / 10 A1 n1 Ea1/
   PLOG /10 A2 n2 Ea2/

and not
PLOG / 10 A1 n1 Ea1/
DUPLICATE

PLOG / 10 A2 n2 Ea2/
DUPLICATE

These two notations have different interpolation consequences. Unfortunately Chemkin documentation is not very clear but Ray Speth tested it out and told me we should have the same pressure twice in single PLOG reaction rather than have two duplicates.

[rwest]

Thanks for bringing this up. I don't quite understand though what should be happening (let alone what is happening) - can you add more explanation?

[rwest]

@edames has just brought this up again (and offered to fix it?!)

[enochd]

I have most of this in the form of a doc that i can send upon request, but the text below clarifies how chemkin interprets PLOG and how duplicate reactions are intended to be used.

In order to illustrate the consequences of providing multiple PLOG expressions and denoting them as 'DUPLICATE' reactions, below are two cases that show how interpolation can arrive at different rate coefficients. These cases are manually computed and compared with rates computed by Chemkin Pro:

In the following example, we want to represent a known reaction with rates that are 1 s-1 and 50 s-1 at 1 and 10 atm, respectively. We will compute the half-life of A in a batch reactor at 2atm. Note that the rate of this isomerization reaction is independent of temperature and only dependent on pressure. Although this does not affect the simulation outcome since the reaction is defined as irreversible, the thermodynamic parameters for the reactant and product have been made equivalent. A constant temperature and pressure homogeneous reactor is simulated in Chemkin-Pro. The species profiles are saved and the half-life time is linearly interpolated from the output:

t1/2 (Case1, Chemkin-Pro) = 0.21 seconds
t1/2 (Case2, Chemkin-Pro) = 0.28 seconds

Below, we manually compute the half-life for comparison with Chemkin-Pro's simulation results. Using the original PLOG formulation introduced by James Miller (X. Gou, J.A. Miller, W. Sun and Y. Ju, http://engine.princeton.edu, 2011):

for a set of pressures p = {p1, p2, . . ., pN}. N is the user specified number of pressures included in the PLOG parameterization. M is the number of duplicate pressures. Thus, the rate at a given pressure is the sum of the internally duplicate reactions. This equation is only valid exactly at one of the pressures pi specified in the PLOG statement. At any other pressure one must interpolate between the two nearest pressures. This is achieved through a linear interpolation between log(k) and log(P):

The above equation is applied to Case 1, obtaining k(2atm) = 3.25 s-1. The corresponding half-life of this species is thus ln(2)/3.25 s-1 = 0.21 seconds, in exact agreement with Chemkin-Pro's interpretation of Case 1.

In Case 2, Chemkin-Pro treats the duplicate reaction as expected (i.e., koverall,P = k1dup,P + k2dup,P). Thus, k(2atm) = 2.47 s-1, yielding a half-life of 0.28 seconds. Thus, the two rate representations are not equivalent. The representation in Case 1 is what was originally intended and allows for correctly interpolated rates coefficients. While this is an extreme example, the 2 atm rate coefficient interpolated in Case 2 yields an incorrect value. In most cases, the differences are much smaller and on the order of a few percent. Regardless, users of the PLOG formalism should adhere to the intended format illustrated by Case 1 in order to yield correctly interpolated rates.

[rwest]

Thanks for clarifying, @enochd.
Once RMG can handle these correctly, someone should fix the Glarborg seed mechanisms in the RMG-database project.

[connie]

Should be fixed now by #329

[nyee]

I have written a script that fixes old RMG mechanisms that suffered from this bug. You can find it here in the Works-In-Progress repo.