In [20]:

    

from rmgpy.chemkin import loadChemkinFile
from rmgpy.kinetics import KineticsData, ArrheniusEP, ThirdBody, Lindemann, Troe, Chebyshev, PDepArrhenius, MultiArrhenius, MultiPDepArrhenius, getRateCoefficientUnitsFromReactionOrder
from rmgpy.reaction import Reaction
import numpy


    

In [2]:

    

path = '/home/mjliu/Documents/ReverseReactions/kinetics.inp'
dict = '/home/mjliu/Documents/ReverseReactions/dictionary.txt'
speciesList, reactionList = loadChemkinFile(path)


    

In [17]:

    

reactionList[0].kinetics


    

    
Out[17]:





PDepArrhenius(pressures=([0.03947,1,10,100],'atm'), arrhenius=[MultiArrhenius(arrhenius=[Arrhenius(A=(3.62e+84,'cm^3/(mol*s)'), n=-21.48, Ea=(42190,'cal/mol'), T0=(1,'K')), Arrhenius(A=(2.16e+43,'cm^3/(mol*s)'), n=-10.07, Ea=(12890,'cal/mol'), T0=(1,'K'))]), MultiArrhenius(arrhenius=[Arrhenius(A=(2.63e+69,'cm^3/(mol*s)'), n=-16.52, Ea=(40770,'cal/mol'), T0=(1,'K')), Arrhenius(A=(2.89e+26,'cm^3/(mol*s)'), n=-4.68, Ea=(7584,'cal/mol'), T0=(1,'K'))]), MultiArrhenius(arrhenius=[Arrhenius(A=(7.59e+49,'cm^3/(mol*s)'), n=-10.65, Ea=(32900,'cal/mol'), T0=(1,'K')), Arrhenius(A=(1.04e+17,'cm^3/(mol*s)'), n=-1.72, Ea=(4254,'cal/mol'), T0=(1,'K'))]), MultiArrhenius(arrhenius=[Arrhenius(A=(1.35e+27,'cm^3/(mol*s)'), n=-4.11, Ea=(18320,'cal/mol'), T0=(1,'K')), Arrhenius(A=(1.27e+20,'cm^3/(mol*s)'), n=-2.72, Ea=(5193,'cal/mol'), T0=(1,'K'))])])

In [4]:

    

len(reactionList)


    

    
Out[4]:





4

In [5]:

    

def writePLOG(kinetics):
    string = ""
    numReactants = 1
    
    if hasattr(kinetics,'highPlimit') and kinetics.highPlimit is not None:
        arrhenius = kinetics.highPlimit
        string += '{0:<9.3e} {1:<9.3f} {2:<9.3f}'.format(
            arrhenius.A.value_si / (arrhenius.T0.value_si ** arrhenius.n.value_si) * 1.0e6 ** (numReactants - 1),
            arrhenius.n.value_si,
            arrhenius.Ea.value_si / 4184.
            )
    else:
        # Print dummy values that Chemkin parses but ignores
        string += '{0:<9.3e} {1:<9.3f} {2:<9.3f}'.format(1, 0, 0)
        
    string += '\n'
    
    for P, arrhenius in zip(kinetics.pressures.value_si, kinetics.arrhenius):
        if isinstance(arrhenius, MultiArrhenius):
            for arrh in arrhenius.arrhenius:
                string += '    PLOG/ {0:<9.3f} {1:<9.3e} {2:<9.3f} {3:<9.3f}/\n'.format(P / 101325.,
                arrh.A.value_si / (arrh.T0.value_si ** arrh.n.value_si) * 1.0e6 ** (numReactants - 1),
                arrh.n.value_si,
                arrh.Ea.value_si / 4184.
                )
        else:
            string += '    PLOG/ {0:<9.3f} {1:<9.3e} {2:<9.3f} {3:<9.3f}/\n'.format(P / 101325.,
                arrhenius.A.value_si / (arrhenius.T0.value_si ** arrhenius.n.value_si) * 1.0e6 ** (numReactants - 1),
                arrhenius.n.value_si,
                arrhenius.Ea.value_si / 4184.
            )
            
    return string


    

In [26]:

    

def reverseRateOld(reaction):
    kunits = getRateCoefficientUnitsFromReactionOrder(len(reaction.products))
            
    kf = reaction.kinetics

    if isinstance(kf, PDepArrhenius):
        if kf.Tmin is not None and kf.Tmax is not None:
            Tlist = 1.0/numpy.linspace(1.0/kf.Tmax.value, 1.0/kf.Tmin.value, 50)
        else:
            Tlist = 1.0/numpy.arange(0.0005, 0.0034, 0.0001)
        Plist = kf.pressures.value_si
        K = numpy.zeros((len(Tlist), len(Plist)), numpy.float64)
        for Tindex, T in enumerate(Tlist):
            for Pindex, P in enumerate(Plist):
                K[Tindex, Pindex] = kf.getRateCoefficient(T, P) / reaction.getEquilibriumConstant(T)
        kr = PDepArrhenius()
        kr.fitToData(Tlist, Plist, K, kunits, kf.arrhenius[0].T0.value)
        return kr


    

In [33]:

    

string1 =  writePLOG(reverseRateOld(reactionList[2]))


    

In [29]:

    

def reverseRateNew(reaction):
    kunits = getRateCoefficientUnitsFromReactionOrder(len(reaction.products))
            
    kf = reaction.kinetics

    if isinstance(kf, PDepArrhenius):
        kr = PDepArrhenius()
        kr.pressures = kf.pressures
        kr.arrhenius = []
        rxn = Reaction(reactants = reaction.reactants, products = reaction.products)            
        
        for kinetics in kf.arrhenius:
            rxn.kinetics = kinetics
            kr.arrhenius.append(rxn.generateReverseRateCoefficient())
            
        return kr


    

In [30]:

    

string2 =  writePLOG(reverseRateNew(reactionList[2]))


    

In [34]:

    

print string1
print string2
print string1 == string2


    

    




1.000e+00 0.000     0.000    
    PLOG/ 0.005     2.416e+11 -0.639    65.688   /
    PLOG/ 0.022     9.317e+11 -0.624    65.688   /
    PLOG/ 0.099     2.046e+12 -0.555    65.689   /

1.000e+00 0.000     0.000    
    PLOG/ 0.005     2.416e+11 -0.639    65.688   /
    PLOG/ 0.022     9.317e+11 -0.624    65.688   /
    PLOG/ 0.099     2.046e+12 -0.555    65.689   /

True

In [35]:

    

print writePLOG(reverseRateNew(reactionList[0]))


    

    




1.000e+00 0.000     0.000    
    PLOG/ 0.039     1.199e+92 -23.015   92.578   /
    PLOG/ 0.039     7.152e+50 -11.605   63.278   /
    PLOG/ 1.000     8.708e+76 -18.055   91.158   /
    PLOG/ 1.000     9.569e+33 -6.215    57.972   /
    PLOG/ 10.000    2.513e+57 -12.185   83.288   /
    PLOG/ 10.000    3.443e+24 -3.255    54.642   /
    PLOG/ 100.000   4.470e+34 -5.645    68.708   /
    PLOG/ 100.000   4.205e+27 -4.255    55.581   /

In [ ]: