Converting pmutt.reaction.Reaction to modified Arrhenius equation

[googhgoo]

Describe the solution you'd like
MKM packages such as chemkin and Cantera uses modified Arrhenius equation which looks like:

k = A*T^(b)*e^(-Ea/RT)

Typically, we run isothermal reactor, so We simply use b = 0, and use A and Ea that would reproduce the transition state theory calculated k at the reactor temperature. This is okay for the isothermal reactor, but for those of who running the nonisothermal reactor needs to nail down the A, b and Ea values well for the correct simulation. It would be nice to have a function from pmutt.reaction.Reaction that fit these values correctly.

Describe alternatives you've considered
Currently, I had to write my own function to fit A, b, and Ea. I calculated Ea first by

rxn.get_E_state('transition_state','eV',include_ZPE=True) - rxn.get_E_state('reactants','eV',include_ZPE=True)

Be aware that you should not use get_E_act, here because it actually calculates get_H_act. (It doesn't make sense anyway that get_E_act needs T. Perhaps we need to differentiate it with get_U_act).

After this, I calculate a line of A values:

Ts = np.linspace(85, 800, 500)
As = [rxn.get_A(T) for T in Ts]

and I have written least-square fitting function to fit A*T^b:

def axb_fitting(x,y):
    n = len(x)
    # Finding required sum for least square methods
    sumX,sumX2,sumY,sumXY = 0,0,0,0
    for i in range(n):
        sumX = sumX + np.log(x[i])
        sumX2 = sumX2 +np.log(x[i])*np.log(x[i])
        sumY = sumY + np.log(y[i])
        sumXY = sumXY + np.log(x[i])*np.log(y[i])
    # Finding coefficients A and B
    b = (n*sumXY-sumX*sumY)/(n*sumX2-sumX*sumX)
    A = (sumY - b*sumX)/n
    # Obtaining a and b
    a = np.exp(A)
    return a,b

I hope this is implemented in the future, and also the difference between E and U should be solved.

[jonlym]

Ah, I see you've already used my solution from #187 for calculating the activation energy using electronic energy and ZPE. Also, see #187 for why we accept T for reaction.get_EoRT_act.

Our ChemkinReaction object has a beta attribute but as you pointed out we never fit it because we always ran isothermal. We'd greatly appreciate if you submit a pull request. That way we could include your fitting function into the official version!

[googhgoo]

Okay sounds good. Where should I make the "adapter" definition? I was never involved with pmutt, so I'm not sure. Should I make it a definition for the reaction class?

[jonlym]

I think a class method similar to the Reaction.from_string method would maintain consistency. Perhaps adding this to ChemkinReaction:

class ChemkinReaction:
    @classmethod
    def from_data(T, k):
        # Insert your fitting algorithm here
        return cls(A=A, beta=beta, Ea=Ea)

However, that class was specifically written for IO to Chemkin so if you find it too inflexible, you can define your own Reaction child class. If you're still interested, you can submit a pull request to the development branch and we can iterate from there.