diff --git a/environment.yml b/environment.yml index 093a15ae75..da98a778f1 100644 --- a/environment.yml +++ b/environment.yml @@ -3,10 +3,11 @@ channels: - defaults - rmg - conda-forge + - cantera dependencies: - cairo - cairocffi - - rmg::cantera >=2.3.0 + - cantera::cantera=2.6 - conda-forge::cclib >=1.6.3 - rmg::chemprop - coolprop diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx index a02806dfb0..8c4856c560 100644 --- a/rmgpy/kinetics/arrhenius.pyx +++ b/rmgpy/kinetics/arrhenius.pyx @@ -227,11 +227,14 @@ cdef class Arrhenius(KineticsModel): """ self._A.value_si *= factor - def to_cantera_kinetics(self): + def to_cantera_kinetics(self, arrhenius_class=False): """ - Converts the Arrhenius object to a cantera Arrhenius object + Converts the RMG Arrhenius object to a cantera ArrheniusRate or + the auxiliary cantera Arrhenius class (used by falloff reactions). + Inputs for both are (A,b,E) where A is in units of m^3/kmol/s, b is dimensionless, and E is in J/kmol - Arrhenius(A,b,E) where A is in units of m^3/kmol/s, b is dimensionless, and E is in J/kmol + arrhenius_class: If ``True``, uses cantera.Arrhenius (for falloff reactions). If ``False``, uses + Cantera.ArrheniusRate """ import cantera as ct @@ -261,15 +264,18 @@ cdef class Arrhenius(KineticsModel): b = self._n.value_si E = self._Ea.value_si * 1000 # convert from J/mol to J/kmol - return ct.Arrhenius(A, b, E) + if arrhenius_class: + return ct.Arrhenius(A, b, E) + else: + return ct.ArrheniusRate(A, b, E) def set_cantera_kinetics(self, ct_reaction, species_list): """ - Passes in a cantera ElementaryReaction() object and sets its - rate to a Cantera Arrhenius() object. + Passes in a cantera Reaction() object and sets its + rate to a Cantera ArrheniusRate object. """ import cantera as ct - assert isinstance(ct_reaction, ct.ElementaryReaction), "Must be a Cantera ElementaryReaction object" + assert isinstance(ct_reaction.rate, ct.ArrheniusRate), "Must have a Cantera ArrheniusRate attribute" # Set the rate parameter to a cantera Arrhenius object ct_reaction.rate = self.to_cantera_kinetics() @@ -449,7 +455,7 @@ cdef class ArrheniusEP(KineticsModel): def set_cantera_kinetics(self, ct_reaction, species_list): """ - Sets a cantera ElementaryReaction() object with the modified Arrhenius object + Sets a cantera Reaction() object with the modified Arrhenius object converted to an Arrhenius form. """ raise NotImplementedError('set_cantera_kinetics() is not implemented for ArrheniusEP class kinetics.') @@ -703,7 +709,7 @@ cdef class ArrheniusBM(KineticsModel): def set_cantera_kinetics(self, ct_reaction, species_list): """ - Sets a cantera ElementaryReaction() object with the modified Arrhenius object + Sets a cantera Reaction() object with the modified Arrhenius object converted to an Arrhenius form. """ raise NotImplementedError('set_cantera_kinetics() is not implemented for ArrheniusBM class kinetics.') @@ -863,12 +869,13 @@ cdef class PDepArrhenius(PDepKineticsModel): """ import cantera as ct import copy - assert isinstance(ct_reaction, ct.PlogReaction), "Must be a Cantera PlogReaction object" + assert isinstance(ct_reaction.rate, ct.PlogRate), "Must have a Cantera PlogRate attribute" pressures = copy.deepcopy(self._pressures.value_si) - ctArrhenius = [arr.to_cantera_kinetics() for arr in self.arrhenius] + ctArrhenius = [arr.to_cantera_kinetics(arrhenius_class=True) for arr in self.arrhenius] - ct_reaction.rates = list(zip(pressures, ctArrhenius)) + new_rates = ct.PlogRate(list(zip(pressures, ctArrhenius))) + ct_reaction.rate = new_rates ################################################################################ diff --git a/rmgpy/kinetics/arrheniusTest.py b/rmgpy/kinetics/arrheniusTest.py index 10449215ad..f814ac49b3 100644 --- a/rmgpy/kinetics/arrheniusTest.py +++ b/rmgpy/kinetics/arrheniusTest.py @@ -233,7 +233,7 @@ def test_change_rate(self): def test_to_cantera_kinetics(self): """ Test that the Arrhenius cantera object can be set properly within - a cantera ElementaryReaction object + a cantera Reaction object """ ctArrhenius = self.arrhenius.to_cantera_kinetics() self.assertAlmostEqual(ctArrhenius.pre_exponential_factor, 1e9, 6) diff --git a/rmgpy/kinetics/chebyshev.pyx b/rmgpy/kinetics/chebyshev.pyx index aa46821489..3b36ab140a 100644 --- a/rmgpy/kinetics/chebyshev.pyx +++ b/rmgpy/kinetics/chebyshev.pyx @@ -268,7 +268,7 @@ cdef class Chebyshev(PDepKineticsModel): """ import cantera as ct import copy - assert isinstance(ct_reaction, ct.ChebyshevReaction), "Must be a Cantera ChebyshevReaction object" + assert isinstance(ct_reaction.rate, ct.ChebyshevRate), "Must have a Cantera Chebychev rate attribute" Tmin = self.Tmin.value_si Tmax = self.Tmax.value_si @@ -294,4 +294,6 @@ cdef class Chebyshev(PDepKineticsModel): except: raise Exception('Chebyshev units {0} not found among accepted units for converting to ' 'Cantera Chebyshev object.'.format(self.kunits)) - ct_reaction.set_parameters(Tmin, Tmax, Pmin, Pmax, coeffs) + + new_chebyshev = ct.ChebyshevRate(temperature_range=(Tmin, Tmax), pressure_range=(Pmin, Pmax), data=coeffs) + ct_reaction.rate = new_chebyshev diff --git a/rmgpy/kinetics/falloff.pyx b/rmgpy/kinetics/falloff.pyx index af188154e9..ad1d595ce8 100644 --- a/rmgpy/kinetics/falloff.pyx +++ b/rmgpy/kinetics/falloff.pyx @@ -220,17 +220,29 @@ cdef class Lindemann(PDepKineticsModel): self.arrheniusLow.change_rate(factor) self.arrheniusHigh.change_rate(factor) + + def set_cantera_kinetics(self, ct_reaction, species_list): """ Sets the efficiencies and kinetics for a cantera reaction. """ import cantera as ct - assert isinstance(ct_reaction, ct.FalloffReaction), "Must be a Cantera FalloffReaction object" + assert isinstance(ct_reaction.rate, ct.LindemannRate), "Must have a Cantera LindemannRate attribute" ct_reaction.efficiencies = PDepKineticsModel.get_cantera_efficiencies(self, species_list) - ct_reaction.high_rate = self.arrheniusHigh.to_cantera_kinetics() - ct_reaction.low_rate = self.arrheniusLow.to_cantera_kinetics() - ct_reaction.falloff = ct.Falloff() + ct_reaction.rate = self.to_cantera_kinetics() + + + def to_cantera_kinetics(self): + """ + Converts the Lindemann object to a cantera LindemannRate object + """ + import cantera as ct + + high_rate = self.arrheniusHigh.to_cantera_kinetics(arrhenius_class=True) + low_rate = self.arrheniusLow.to_cantera_kinetics(arrhenius_class=True) + return ct.LindemannRate(low=low_rate, high=high_rate) + ################################################################################ @@ -389,15 +401,30 @@ cdef class Troe(PDepKineticsModel): """ import cantera as ct - assert isinstance(ct_reaction, ct.FalloffReaction), "Must be a Cantera FalloffReaction object" + assert isinstance(ct_reaction.rate, ct.TroeRate), "Must have a Cantera TroeRate attribute" + ct_reaction.efficiencies = PDepKineticsModel.get_cantera_efficiencies(self, species_list) - ct_reaction.high_rate = self.arrheniusHigh.to_cantera_kinetics() - ct_reaction.low_rate = self.arrheniusLow.to_cantera_kinetics() + + ct_reaction.rate = self.to_cantera_kinetics() + + def to_cantera_kinetics(self): + """ + Converts the Troe object to a cantera Troe object + """ + import cantera as ct A = self.alpha T3 = self.T3.value_si T1 = self.T1.value_si if self.T2 is None: - ct_reaction.falloff = ct.TroeFalloff(params=[A, T3, T1]) + falloff = [A, T3, T1] else: T2 = self.T2.value_si - ct_reaction.falloff = ct.TroeFalloff(params=[A, T3, T1, T2]) + falloff = [A, T3, T1, T2] + + high = self.arrheniusHigh.to_cantera_kinetics(arrhenius_class=True) + low = self.arrheniusLow.to_cantera_kinetics(arrhenius_class=True) + return ct.TroeRate(high=high, low=low, falloff_coeffs=falloff) + + + + \ No newline at end of file diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index 0fd5188466..e5f49cb13f 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -281,39 +281,81 @@ def to_cantera(self, species_list=None, use_chemkin_identifier=False): ct_products[product_name] += 1 else: ct_products[product_name] = 1 + if self.specific_collider: # add a specific collider if exists ct_collider[self.specific_collider.to_chemkin() if use_chemkin_identifier else self.specific_collider.label] = 1 if self.kinetics: if isinstance(self.kinetics, Arrhenius): # Create an Elementary Reaction - ct_reaction = ct.ElementaryReaction(reactants=ct_reactants, products=ct_products) + ct_reaction = ct.Reaction(reactants=ct_reactants, products=ct_products, rate=ct.ArrheniusRate()) elif isinstance(self.kinetics, MultiArrhenius): # Return a list of elementary reactions which are duplicates - ct_reaction = [ct.ElementaryReaction(reactants=ct_reactants, products=ct_products) + ct_reaction = [ct.Reaction(reactants=ct_reactants, products=ct_products, rate=ct.ArrheniusRate()) for arr in self.kinetics.arrhenius] elif isinstance(self.kinetics, PDepArrhenius): - ct_reaction = ct.PlogReaction(reactants=ct_reactants, products=ct_products) + ct_reaction = ct.Reaction(reactants=ct_reactants, products=ct_products, rate=ct.PlogRate()) elif isinstance(self.kinetics, MultiPDepArrhenius): - ct_reaction = [ct.PlogReaction(reactants=ct_reactants, products=ct_products) + ct_reaction = [ct.Reaction(reactants=ct_reactants, products=ct_products, rate=ct.PlogRate()) for arr in self.kinetics.arrhenius] elif isinstance(self.kinetics, Chebyshev): - ct_reaction = ct.ChebyshevReaction(reactants=ct_reactants, products=ct_products) + ct_reaction = ct.Reaction(reactants=ct_reactants, products=ct_products, rate=ct.ChebyshevRate()) elif isinstance(self.kinetics, ThirdBody): if ct_collider is not None: - ct_reaction = ct.ThreeBodyReaction(reactants=ct_reactants, products=ct_products, tbody=ct_collider) + ct_reaction = ct.ThreeBodyReaction(reactants=ct_reactants, products=ct_products, third_body=ct_collider) else: ct_reaction = ct.ThreeBodyReaction(reactants=ct_reactants, products=ct_products) - elif isinstance(self.kinetics, Lindemann) or isinstance(self.kinetics, Troe): + elif isinstance(self.kinetics, Troe): + high_rate = self.kinetics.arrheniusHigh.to_cantera_kinetics(arrhenius_class=True) + low_rate = self.kinetics.arrheniusLow.to_cantera_kinetics(arrhenius_class=True) + A = self.kinetics.alpha + T3 = self.kinetics.T3.value_si + T1 = self.kinetics.T1.value_si + + if self.kinetics.T2 is None: + rate = ct.TroeRate( + high=high_rate, low=low_rate, falloff_coeffs=[A, T3, T1] + ) + else: + T2 = self.kinetics.T2.value_si + rate = ct.TroeRate( + high=high_rate, low=low_rate, falloff_coeffs=[A, T3, T1, T2] + ) + if ct_collider is not None: - ct_reaction = ct.FalloffReaction(reactants=ct_reactants, products=ct_products, tbody=ct_collider) + ct_reaction = ct.FalloffReaction( + reactants=ct_reactants, + products=ct_products, + tbody=ct_collider, + rate=rate, + ) else: - ct_reaction = ct.FalloffReaction(reactants=ct_reactants, products=ct_products) + ct_reaction = ct.FalloffReaction( + reactants=ct_reactants, products=ct_products, rate=rate + ) + + elif isinstance(self.kinetics, Lindemann): + high_rate = self.kinetics.arrheniusHigh.to_cantera_kinetics(arrhenius_class=True) + low_rate = self.kinetics.arrheniusLow.to_cantera_kinetics(arrhenius_class=True) + falloff = [] + rate = ct.LindemannRate(low_rate, high_rate, falloff) + if ct_collider is not None: + ct_reaction = ct.FalloffReaction( + reactants=ct_reactants, + products=ct_products, + tbody=ct_collider, + rate=rate, + ) + else: + ct_reaction = ct.FalloffReaction( + reactants=ct_reactants, products=ct_products, rate=rate + ) + else: raise NotImplementedError('Unable to set cantera kinetics for {0}'.format(self.kinetics)) diff --git a/rmgpy/reactionTest.py b/rmgpy/reactionTest.py index 64e8eb002a..a0d60fed45 100644 --- a/rmgpy/reactionTest.py +++ b/rmgpy/reactionTest.py @@ -35,6 +35,7 @@ import cantera as ct import numpy +import yaml from external.wip import work_in_progress from copy import deepcopy @@ -1505,6 +1506,161 @@ def setUp(self): self.species_list = [ch3, ethane, co2, ch4, h2o, ar, h2, h, oh, ho2, o2, co, h2o2] + + yaml_def = ''' +phases: +- name: gas + thermo: ideal-gas + kinetics: gas + elements: [O, H, C, Ar] + species: [ethane, Ar, H2(2), H(3), OH(4), HO2(5), O2(6), H2O2(7), CO(9), CH3(13), CH4(15), CO2(16), H2O(27)] + state: {T: 300.0, P: 1 atm} +species: +- name: ethane + composition: {C: 2, H: 6} + thermo: + model: NASA7 + temperature-ranges: [100.0, 954.51, 5000.0] + data: + - [3.78033, -3.24263e-03, 5.52381e-05, -6.38581e-08, 2.28637e-11, -1.16203e+04, + 5.21034] + - [4.58983, 0.0141508, -4.75962e-06, 8.60294e-10, -6.21717e-14, -1.27218e+04, + -3.61739] +- name: Ar + composition: {Ar: 1} + thermo: + model: NASA7 + temperature-ranges: [200.0, 6000.0] + data: + - [2.5, 0.0, 0.0, 0.0, 0.0, -745.375, 4.37967] +- name: H2(2) + composition: {H: 2} + thermo: + model: NASA7 + temperature-ranges: [100.0, 1959.08, 5000.0] + data: + - [3.43536, 2.1271e-04, -2.78625e-07, 3.40267e-10, -7.76031e-14, -1031.36, + -3.90842] + - [2.78816, 5.87644e-04, 1.59009e-07, -5.52736e-11, 4.34309e-15, -596.143, + 0.112747] +- name: H(3) + composition: {H: 1} + thermo: + model: NASA7 + temperature-ranges: [100.0, 4563.27, 5000.0] + data: + - [2.5, -1.91243e-12, 2.45329e-15, -1.02377e-18, 1.31369e-22, 2.54742e+04, + -0.444973] + - [2.50167, -1.43051e-06, 4.6025e-10, -6.57826e-14, 3.52412e-18, 2.54727e+04, + -0.455578] +- name: OH(4) + composition: {H: 1, O: 1} + thermo: + model: NASA7 + temperature-ranges: [100.0, 1145.75, 5000.0] + data: + - [3.51457, 2.92773e-05, -5.32163e-07, 1.01949e-09, -3.85945e-13, 3414.25, + 2.10435] + - [3.07194, 6.04016e-04, -1.39783e-08, -2.13446e-11, 2.48066e-15, 3579.39, + 4.578] +- name: HO2(5) + composition: {H: 1, O: 2} + thermo: + model: NASA7 + temperature-ranges: [100.0, 932.15, 5000.0] + data: + - [4.04594, -1.73464e-03, 1.03766e-05, -1.02202e-08, 3.34908e-12, -986.754, + 4.63581] + - [3.21024, 3.67942e-03, -1.27701e-06, 2.18045e-10, -1.46338e-14, -910.369, + 8.18291] +- name: O2(6) + composition: {O: 2} + thermo: + model: NASA7 + temperature-ranges: [100.0, 1074.55, 5000.0] + data: + - [3.53732, -1.21572e-03, 5.3162e-06, -4.89446e-09, 1.45846e-12, -1038.59, + 4.68368] + - [3.15382, 1.67804e-03, -7.69974e-07, 1.51275e-10, -1.08782e-14, -1040.82, + 6.16756] +- name: H2O2(7) + composition: {H: 2, O: 2} + thermo: + model: NASA7 + temperature-ranges: [100.0, 908.87, 5000.0] + data: + - [3.73136, 3.35071e-03, 9.35033e-06, -1.521e-08, 6.41585e-12, -1.77212e+04, + 5.45911] + - [5.41579, 2.61008e-03, -4.39892e-07, 4.91087e-11, -3.35188e-15, -1.8303e+04, + -4.02248] +- name: CO(9) + composition: {C: 1, O: 1} + thermo: + model: NASA7 + temperature-ranges: [100.0, 884.77, 5000.0] + data: + - [3.66965, -5.50953e-03, 2.00538e-05, -2.08391e-08, 7.43738e-12, 1200.77, + -12.4224] + - [2.8813, 2.31665e-03, -4.40151e-07, 4.75633e-11, -2.78282e-15, 1173.45, + -9.65831] +- name: CH3(13) + composition: {C: 1, H: 3} + thermo: + model: NASA7 + temperature-ranges: [100.0, 1337.62, 5000.0] + data: + - [3.91547, 1.84154e-03, 3.48744e-06, -3.3275e-09, 8.49964e-13, 1.62856e+04, + 0.351739] + - [3.54145, 4.76788e-03, -1.82149e-06, 3.28878e-10, -2.22547e-14, 1.6224e+04, + 1.6604] +- name: CH4(15) + composition: {C: 1, H: 4} + thermo: + model: NASA7 + temperature-ranges: [100.0, 1084.12, 5000.0] + data: + - [4.20541, -5.35556e-03, 2.51123e-05, -2.13762e-08, 5.97522e-12, -1.01619e+04, + -0.921275] + - [0.908272, 0.0114541, -4.57173e-06, 8.2919e-10, -5.66314e-14, -9719.98, + 13.9931] +- name: CO2(16) + composition: {C: 1, O: 2} + thermo: + model: NASA7 + temperature-ranges: [100.0, 988.89, 5000.0] + data: + - [3.27861, 2.74152e-03, 7.16074e-06, -1.08027e-08, 4.14282e-12, -4.84703e+04, + 5.97937] + - [4.5461, 2.91913e-03, -1.15484e-06, 2.27654e-10, -1.7091e-14, -4.89804e+04, + -1.43275] +- name: H2O(27) + composition: {H: 2, O: 1} + thermo: + model: NASA7 + temperature-ranges: [100.0, 1130.24, 5000.0] + data: + - [4.05764, -7.87933e-04, 2.90876e-06, -1.47518e-09, 2.12838e-13, -3.02816e+04, + -0.311363] + - [2.84325, 2.75108e-03, -7.8103e-07, 1.07243e-10, -5.79389e-15, -2.99586e+04, + 5.91041] +reactions: +- equation: H(3) + H(3) + M <=> H2(2) + M + type: three-body + rate-constant: + A: 1.000000e+18 + b: -1.0 + Ea: 0.0 + efficiencies: + CO2(16): 0.0 + CH4(15): 2.0 + ethane: 3.0 + H2O(27): 0.0 + H2(2): 0.0 + Ar: 0.63 + ''' + + gas = ct.Solution(yaml=yaml_def) + self.troe = Reaction(index=1, reactants=[ch3, ch3], products=[ethane], kinetics=Troe( arrheniusHigh=Arrhenius(A=(6.77e+16, 'cm^3/(mol*s)'), n=-1.18, Ea=(0.654, 'kcal/mol'), @@ -1515,40 +1671,52 @@ def setUp(self): efficiencies={Molecule(smiles="O=C=O"): 2.0, Molecule(smiles="[H][H]"): 2.0, Molecule(smiles="O"): 6.0, Molecule(smiles="[Ar]"): 0.7, Molecule(smiles="C"): 2.0, Molecule(smiles="CC"): 3.0})) - - self.ct_troe = ct.Reaction.fromCti('''falloff_reaction('CH3(13) + CH3(13) (+ M) <=> ethane (+ M)', - kf=[(6.770000e+16,'cm3/mol/s'), -1.18, (0.654,'kcal/mol')], - kf0=[(3.400000e+41,'cm6/mol2/s'), -7.03, (2.762,'kcal/mol')], - efficiencies='ethane:3.0 CO2(16):2.0 CH4(15):2.0 Ar:0.7 H2O(27):6.0 H2(2):2.0', - falloff=Troe(A=0.619, T3=73.2, T1=1180.0, T2=10000.0))''') - + + self.ct_troe = ct.Reaction.from_yaml(''' + equation: CH3(13) + CH3(13) (+M) <=> ethane (+M) # Reaction 1 + type: falloff + low-P-rate-constant: {A: 3.4e+41 cm^6/mol^2/s, b: -7.03, Ea: 2.762 kcal/mol} + high-P-rate-constant: {A: 6.77e+16 cm^3/mol/s, b: -1.18, Ea: 0.654 kcal/mol} + Troe: {A: 0.619, T3: 73.2, T1: 1180.0, T2: 1.0e+04} + efficiencies: {CH4(15): 2.0, Ar: 0.7, ethane: 3.0, H2(2): 2.0, CO2(16): 2.0, + H2O(27): 6.0}''',gas) + self.arrheniusBi = Reaction(index=2, reactants=[h, ch4], products=[h2, ch3], kinetics=Arrhenius(A=(6.6e+08, 'cm^3/(mol*s)'), n=1.62, Ea=(10.84, 'kcal/mol'), T0=(1, 'K'))) - self.ct_arrheniusBi = ct.Reaction.fromCti( - '''reaction('H(3) + CH4(15) <=> H2(2) + CH3(13)', [(6.600000e+08,'cm3/mol/s'), 1.62, (10.84,'kcal/mol')])''') - + self.ct_arrheniusBi = ct.Reaction.from_yaml(''' + equation: H(3) + CH4(15) <=> H2(2) + CH3(13) + rate-constant: {A: 6.6e+08 cm^3/mol/s, b: 1.62, Ea: 10.84 kcal/mol} + ''',gas) + self.arrheniusBi_irreversible = Reaction(index=10, reactants=[h, ch4], products=[h2, ch3], kinetics=Arrhenius(A=(6.6e+08, 'cm^3/(mol*s)'), n=1.62, Ea=(10.84, 'kcal/mol'), T0=(1, 'K')), reversible=False) - self.ct_arrheniusBi_irreversible = ct.Reaction.fromCti( - '''reaction('H(3) + CH4(15) => H2(2) + CH3(13)', [(6.600000e+08,'cm3/mol/s'), 1.62, (10.84,'kcal/mol')])''') + self.ct_arrheniusBi_irreversible = ct.Reaction.from_yaml(''' + equation: H(3) + CH4(15) => H2(2) + CH3(13) # Reaction 3 + rate-constant: {A: 6.6e+08 cm^3/mol/s, b: 1.62, Ea: 10.84 kcal/mol} + ''',gas) self.arrheniusMono = Reaction(index=15, reactants=[h2o2], products=[h2, o2], kinetics=Arrhenius(A=(6.6e+03, '1/s'), n=1.62, Ea=(10.84, 'kcal/mol'), T0=(1, 'K'))) - - self.ct_arrheniusMono = ct.Reaction.fromCti( - '''reaction('H2O2(7) <=> H2(2) + O2(6)', [(6.600000e+03,'1/s'), 1.62, (10.84,'kcal/mol')])''') + + self.ct_arrheniusMono = ct.Reaction.from_yaml(''' + equation: H2O2(7) <=> H2(2) + O2(6) # Reaction 4 + rate-constant: {A: 6600.0 1/s, b: 1.62, Ea: 10.84 kcal/mol} + ''',gas) self.arrheniusTri = Reaction(index=20, reactants=[h, h, o2], products=[h2o2], kinetics=Arrhenius(A=(6.6e+08, 'cm^6/(mol^2*s)'), n=1.62, Ea=(10.84, 'kcal/mol'), T0=(1, 'K'))) - self.ct_arrheniusTri = ct.Reaction.fromCti( - '''reaction('H(3) + H(3) + O2(6) <=> H2O2(7)', [(6.6e+08, 'cm6/mol2/s'), 1.62, (10.84,'kcal/mol')])''') + + self.ct_arrheniusTri = ct.Reaction.from_yaml(''' + equation: H(3) + H(3) + O2(6) <=> H2O2(7) # Reaction 5 + rate-constant: {A: 6.6e+08 cm^6/mol^2/s, b: 1.62, Ea: 10.84 kcal/mol} + ''',gas) self.multiArrhenius = Reaction(index=3, reactants=[oh, ho2], products=[h2o, o2], kinetics=MultiArrhenius(arrhenius=[ @@ -1556,10 +1724,17 @@ def setUp(self): T0=(1, 'K')), Arrhenius(A=(5e+15, 'cm^3/(mol*s)'), n=0, Ea=(17.33, 'kcal/mol'), T0=(1, 'K'))])) - - self.ct_multiArrhenius = [ct.Reaction.fromCti('''reaction('OH(4) + HO2(5) <=> H2O(27) + O2(6)', [(1.450000e+13,'cm3/mol/s'), 0.0, (-0.5,'kcal/mol')], - options='duplicate')'''), ct.Reaction.fromCti('''reaction('OH(4) + HO2(5) <=> H2O(27) + O2(6)', [(5.000000e+15,'cm3/mol/s'), 0.0, (17.33,'kcal/mol')], - options='duplicate')''')] + self.ct_multiArrhenius = [ + ct.Reaction.from_yaml(''' + equation: OH(4) + HO2(5) <=> H2O(27) + O2(6) # Reaction 6 + duplicate: true + rate-constant: {A: 1.45e+13 cm^3/mol/s, b: 0.0, Ea: -0.5 kcal/mol} + ''',gas), + ct.Reaction.from_yaml(''' + equation: OH(4) + HO2(5) <=> H2O(27) + O2(6) # Reaction 7 + duplicate: true + rate-constant: {A: 5.0e+15 cm^3/mol/s, b: 0.0, Ea: 17.33 kcal/mol} + ''',gas)] self.pdepArrhenius = Reaction(index=4, reactants=[ho2, ho2], products=[o2, h2o2], kinetics=PDepArrhenius( @@ -1586,11 +1761,15 @@ def setUp(self): ], ), ) - - self.ct_pdepArrhenius = ct.Reaction.fromCti('''pdep_arrhenius('HO2(5) + HO2(5) <=> O2(6) + H2O2(7)', - [(0.1, 'atm'), (8.800000e+16, 'cm3/mol/s'), -1.05, (6.461,'kcal/mol')], - [(1.0, 'atm'), (8.000000e+21,'cm3/mol/s'), -2.39, (11.18,'kcal/mol')], - [(10.0, 'atm'), (3.300000e+24,'cm3/mol/s'), -3.04, (15.61,'kcal/mol')])''') + + self.ct_pdepArrhenius = ct.Reaction.from_yaml(''' + equation: HO2(5) + HO2(5) <=> O2(6) + H2O2(7) # Reaction 8 + type: pressure-dependent-Arrhenius + rate-constants: + - {P: 0.1 atm, A: 8.8e+16 cm^3/mol/s, b: -1.05, Ea: 6.461 kcal/mol} + - {P: 1.0 atm, A: 8.0e+21 cm^3/mol/s, b: -2.39, Ea: 11.18 kcal/mol} + - {P: 10.0 atm, A: 3.3e+24 cm^3/mol/s, b: -3.04, Ea: 15.61 kcal/mol} + ''',gas) self.multiPdepArrhenius = Reaction(index=5, reactants=[ho2, ch3], products=[o2, ch4], kinetics=MultiPDepArrhenius( @@ -1620,17 +1799,25 @@ def setUp(self): ], ), ) - - self.ct_multiPdepArrhenius = [ct.Reaction.fromCti('''pdep_arrhenius('HO2(5) + CH3(13) <=> O2(6) + CH4(15)', - [(0.001, 'atm'), (9.300000e+10, 'cm3/mol/s'), 0.0, (0.0,'kcal/mol')], - [(1.0, 'atm'), (8.000000e+10, 'cm3/mol/s'), 0.0, (0.0,'kcal/mol')], - [(3.0, 'atm'), (7.000000e+10, 'cm3/mol/s'), 0.0, (0.0,'kcal/mol')], - options='duplicate')'''), - ct.Reaction.fromCti('''pdep_arrhenius('HO2(5) + CH3(13) <=> O2(6) + CH4(15)', - [(0.001, 'atm'), (7.100000e+05, 'cm3/mol/s'), 1.8, (1.133,'kcal/mol')], - [(1.0, 'atm'), (8.800000e+05, 'cm3/mol/s'), 1.77, (0.954,'kcal/mol')], - [(3.0, 'atm'), (2.900000e+05, 'cm3/mol/s'), 1.9, (0.397,'kcal/mol')], - options='duplicate')''')] + self.ct_multiPdepArrhenius = [ + ct.Reaction.from_yaml(''' + equation: HO2(5) + CH3(13) <=> O2(6) + CH4(15) # Reaction 9 + duplicate: true + type: pressure-dependent-Arrhenius + rate-constants: + - {P: 0.001 atm, A: 9.3e+10 cm^3/mol/s, b: 0.0, Ea: 0.0 kcal/mol} + - {P: 1.0 atm, A: 8.0e+10 cm^3/mol/s, b: 0.0, Ea: 0.0 kcal/mol} + - {P: 3.0 atm, A: 7.0e+10 cm^3/mol/s, b: 0.0, Ea: 0.0 kcal/mol} + ''',gas), + ct.Reaction.from_yaml(''' + equation: HO2(5) + CH3(13) <=> O2(6) + CH4(15) # Reaction 10 + duplicate: true + type: pressure-dependent-Arrhenius + rate-constants: + - {P: 0.001 atm, A: 7.1e+05 cm^3/mol/s, b: 1.8, Ea: 1.133 kcal/mol} + - {P: 1.0 atm, A: 8.8e+05 cm^3/mol/s, b: 1.77, Ea: 0.954 kcal/mol} + - {P: 3.0 atm, A: 2.9e+05 cm^3/mol/s, b: 1.9, Ea: 0.397 kcal/mol} + ''',gas)] self.chebyshev = Reaction(index=6, reactants=[h, ch3], products=[ch4], kinetics=Chebyshev( coeffs=[[12.68, 0.3961, -0.05481, -0.003606], [-0.7128, 0.731, -0.0941, -0.008587], @@ -1638,15 +1825,19 @@ def setUp(self): [-0.2403, 0.1779, 0.01946, -0.008505], [-0.1133, 0.0485, 0.03121, -0.002955]], kunits='cm^3/(mol*s)', Tmin=(300, 'K'), Tmax=(3000, 'K'), Pmin=(0.001, 'atm'), Pmax=(98.692, 'atm'))) - self.ct_chebyshev = ct.Reaction.fromCti('''chebyshev_reaction('H(3) + CH3(13) (+ M) <=> CH4(15) (+ M)', - Tmin=300.0, Tmax=3000.0, - Pmin=(0.001, 'atm'), Pmax=(98.692, 'atm'), - coeffs=[[ 9.68000e+00, 3.96100e-01, -5.48100e-02, -3.60600e-03], - [-7.12800e-01, 7.31000e-01, -9.41000e-02, -8.58700e-03], - [-5.80600e-01, 5.70000e-01, -5.53900e-02, -1.11500e-02], - [-4.07400e-01, 3.65300e-01, -1.18000e-02, -1.17100e-02], - [-2.40300e-01, 1.77900e-01, 1.94600e-02, -8.50500e-03], - [-1.13300e-01, 4.85000e-02, 3.12100e-02, -2.95500e-03]])''') + self.ct_chebyshev = ct.Reaction.from_yaml(''' + equation: H(3) + CH3(13) <=> CH4(15) # Reaction 11 + type: Chebyshev + temperature-range: [300.0, 3000.0] + pressure-range: [0.001 atm, 98.692 atm] + data: + - [9.680, 0.3961, -0.05481, -3.606e-03] + - [-0.7128, 0.731, -0.0941, -8.587e-03] + - [-0.5806, 0.57, -0.05539, -0.01115] + - [-0.4074, 0.3653, -0.0118, -0.01171] + - [-0.2403, 0.1779, 0.01946, -8.505e-03] + - [-0.1133, 0.0485, 0.03121, -2.955e-03] + ''',gas) self.thirdBody = Reaction(index=7, reactants=[h, h], products=[h2], kinetics=ThirdBody( @@ -1657,8 +1848,13 @@ def setUp(self): Molecule(smiles="[Ar]"): 0.63, Molecule(smiles="C"): 2.0, Molecule(smiles="CC"): 3.0})) - self.ct_thirdBody = ct.Reaction.fromCti('''three_body_reaction('H(3) + H(3) + M <=> H2(2) + M', [(1.000000e+18,'cm6/mol2/s'), -1.0, (0.0,'kcal/mol')], - efficiencies='CO2(16):0.0 CH4(15):2.0 ethane:3.0 H2O(27):0.0 H2(2):0.0 Ar:0.63')''') + self.ct_thirdBody = ct.Reaction.from_yaml(''' + equation: H(3) + H(3) + M <=> H2(2) + M # Reaction 12 + type: three-body + rate-constant: {A: 1.0e+18 cm^6/mol^2/s, b: -1.0, Ea: 0.0 kcal/mol} + efficiencies: {H2(2): 0.0, CH4(15): 2.0, CO2(16): 0.0, Ar: 0.63, H2O(27): 0.0, + ethane: 3.0} + ''',gas) self.lindemann = Reaction(index=8, reactants=[h, o2], products=[ho2], kinetics=Lindemann( @@ -1670,11 +1866,15 @@ def setUp(self): Molecule(smiles="O"): 6.0, Molecule(smiles="[Ar]"): 0.5, Molecule(smiles="C"): 2.0, Molecule(smiles="CC"): 3.0, Molecule(smiles="[O][O]"): 6.0})) - - self.ct_lindemann = ct.Reaction.fromCti('''falloff_reaction('H(3) + O2(6) (+ M) <=> HO2(5) (+ M)', - kf=[(1.800000e+10,'cm3/mol/s'), 0.0, (2.385,'kcal/mol')], - kf0=[(6.020000e+14,'cm6/mol2/s'), 0.0, (3.0,'kcal/mol')], - efficiencies='CO2(16):3.5 CH4(15):2.0 ethane:3.0 H2O(27):6.0 O2(6):6.0 H2(2):2.0 Ar:0.5')''') + + self.ct_lindemann = ct.Reaction.from_yaml(''' + equation: H(3) + O2(6) (+M) <=> HO2(5) (+M) # Reaction 13 + type: falloff + low-P-rate-constant: {A: 6.02e+14 cm^6/mol^2/s, b: 0.0, Ea: 3.0 kcal/mol} + high-P-rate-constant: {A: 1.8e+10 cm^3/mol/s, b: 0.0, Ea: 2.385 kcal/mol} + efficiencies: {H2O(27): 6.0, Ar: 0.5, CH4(15): 2.0, O2(6): 6.0, ethane: 3.0, + H2(2): 2.0, CO2(16): 3.5} + ''',gas) def test_arrhenius(self): """ @@ -1691,8 +1891,9 @@ def test_arrhenius(self): self.assertEqual(type(converted_obj), type(ct_obj)) # Check that the reaction string is the same self.assertEqual(repr(converted_obj), repr(ct_obj)) - # Check that the Arrhenius string is identical - self.assertEqual(str(converted_obj.rate), str(ct_obj.rate)) + # Check that the rate is the same. arrhenius string is not going to be identical + self.assertEqual(converted_obj.rate.input_data, ct_obj.rate.input_data) + def test_multi_arrhenius(self): """ @@ -1712,7 +1913,10 @@ def test_multi_arrhenius(self): # Check that the reaction string is the same self.assertEqual(repr(converted_rxn), repr(ct_rxn)) # Check that the Arrhenius rates are identical - self.assertEqual(str(converted_rxn.rate), str(ct_rxn.rate)) + self.assertAlmostEqual(converted_rxn.rate.pre_exponential_factor, ct_rxn.rate.pre_exponential_factor, places=3) + self.assertAlmostEqual(converted_rxn.rate.temperature_exponent, ct_rxn.rate.temperature_exponent) + self.assertAlmostEqual(converted_rxn.rate.activation_energy, ct_rxn.rate.activation_energy) + # self.assertEqual(converted_rxn.rate.input_data, ct_rxn.rate.input_data) def test_pdep_arrhenius(self): """ @@ -1756,41 +1960,34 @@ def test_chebyshev(self): Tests formation of cantera reactions with Chebyshev kinetics. """ ct_chebyshev = self.chebyshev.to_cantera(self.species_list, use_chemkin_identifier=True) - self.assertEqual(type(ct_chebyshev), type(self.ct_chebyshev)) - self.assertEqual(repr(ct_chebyshev), repr(self.ct_chebyshev)) - - self.assertEqual(ct_chebyshev.Tmax, self.ct_chebyshev.Tmax) - self.assertEqual(ct_chebyshev.Tmin, self.ct_chebyshev.Tmin) - self.assertEqual(ct_chebyshev.Pmax, self.ct_chebyshev.Pmax) - self.assertEqual(ct_chebyshev.Pmin, self.ct_chebyshev.Pmin) - self.assertTrue((ct_chebyshev.coeffs == self.ct_chebyshev.coeffs).all()) + self.assertEqual(type(ct_chebyshev.rate), type(self.ct_chebyshev.rate)) + self.assertEqual(ct_chebyshev.rate.temperature_range, self.ct_chebyshev.rate.temperature_range) + self.assertEqual(ct_chebyshev.rate.pressure_range, self.ct_chebyshev.rate.pressure_range) + self.assertTrue((ct_chebyshev.rate.data == self.ct_chebyshev.rate.data).all()) def test_falloff(self): """ Tests formation of cantera reactions with Falloff kinetics. """ + ct_troe = self.troe.to_cantera(self.species_list, use_chemkin_identifier=True) + self.assertEqual(type(ct_troe.rate), type(self.ct_troe.rate)) + self.assertAlmostEqual(ct_troe.rate.low_rate.pre_exponential_factor, self.ct_troe.rate.low_rate.pre_exponential_factor, places=3) + self.assertEqual(ct_troe.rate.low_rate.temperature_exponent, self.ct_troe.rate.low_rate.temperature_exponent) + self.assertEqual(ct_troe.rate.low_rate.activation_energy, self.ct_troe.rate.low_rate.activation_energy) + self.assertEqual(ct_troe.efficiencies, self.ct_troe.efficiencies) + ct_third_body = self.thirdBody.to_cantera(self.species_list, use_chemkin_identifier=True) - self.assertEqual(type(ct_third_body), type(self.ct_thirdBody)) - self.assertEqual(repr(ct_third_body), repr(self.ct_thirdBody)) - self.assertEqual(str(ct_third_body.rate), str(self.ct_thirdBody.rate)) + self.assertEqual(type(ct_third_body.rate), type(self.ct_thirdBody.rate)) + self.assertAlmostEqual(ct_third_body.rate.pre_exponential_factor, self.ct_thirdBody.rate.pre_exponential_factor, places=3) + self.assertEqual(ct_third_body.rate.temperature_exponent, self.ct_thirdBody.rate.temperature_exponent) + self.assertEqual(ct_third_body.rate.activation_energy, self.ct_thirdBody.rate.activation_energy) self.assertEqual(ct_third_body.efficiencies, self.ct_thirdBody.efficiencies) ct_lindemann = self.lindemann.to_cantera(self.species_list, use_chemkin_identifier=True) - self.assertEqual(type(ct_lindemann), type(self.ct_lindemann)) - self.assertEqual(repr(ct_lindemann), repr(self.ct_lindemann)) + self.assertEqual(type(ct_lindemann.rate), type(self.ct_lindemann.rate)) self.assertEqual(ct_lindemann.efficiencies, self.ct_lindemann.efficiencies) - self.assertEqual(str(ct_lindemann.low_rate), str(self.ct_lindemann.low_rate)) - self.assertEqual(str(ct_lindemann.high_rate), str(self.ct_lindemann.high_rate)) - self.assertEqual(str(ct_lindemann.falloff), str(self.ct_lindemann.falloff)) - - ct_troe = self.troe.to_cantera(self.species_list, use_chemkin_identifier=True) - self.assertEqual(type(ct_troe), type(self.ct_troe)) - self.assertEqual(repr(ct_troe), repr(self.ct_troe)) - self.assertEqual(ct_troe.efficiencies, self.ct_troe.efficiencies) - - self.assertEqual(str(ct_troe.low_rate), str(self.ct_troe.low_rate)) - self.assertEqual(str(ct_troe.high_rate), str(self.ct_troe.high_rate)) - self.assertEqual(str(ct_troe.falloff), str(self.ct_troe.falloff)) + self.assertEqual(str(ct_lindemann.rate.low_rate), str(self.ct_lindemann.rate.low_rate)) + self.assertEqual(str(ct_lindemann.rate.high_rate), str(self.ct_lindemann.rate.high_rate)) ################################################################################ diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index 79797d3629..4af6a43256 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -48,7 +48,7 @@ import numpy as np import psutil import yaml -from cantera import ck2cti +from cantera import ck2yaml from scipy.optimize import brute import rmgpy.util as util @@ -1098,14 +1098,14 @@ def execute(self, initialize=True, **kwargs): self.run_model_analysis() - # generate Cantera files chem.cti & chem_annotated.cti in a designated `cantera` output folder + # generate Cantera files chem.yaml & chem_annotated.yaml in a designated `cantera` output folder try: if any([s.contains_surface_site() for s in self.reaction_model.core.species]): self.generate_cantera_files(os.path.join(self.output_directory, 'chemkin', 'chem-gas.inp'), - surfaceFile=( + surface_file=( os.path.join(self.output_directory, 'chemkin', 'chem-surface.inp'))) self.generate_cantera_files(os.path.join(self.output_directory, 'chemkin', 'chem_annotated-gas.inp'), - surfaceFile=(os.path.join(self.output_directory, 'chemkin', + surface_file=(os.path.join(self.output_directory, 'chemkin', 'chem_annotated-surface.inp'))) else: # gas phase only self.generate_cantera_files(os.path.join(self.output_directory, 'chemkin', 'chem.inp')) @@ -1671,13 +1671,13 @@ def process_reactions_to_species(self, obj): def generate_cantera_files(self, chemkin_file, **kwargs): """ - Convert a chemkin mechanism chem.inp file to a cantera mechanism file chem.cti + Convert a chemkin mechanism chem.inp file to a cantera mechanism file chem.yaml and save it in the cantera directory """ transport_file = os.path.join(os.path.dirname(chemkin_file), 'tran.dat') - file_name = os.path.splitext(os.path.basename(chemkin_file))[0] + '.cti' + file_name = os.path.splitext(os.path.basename(chemkin_file))[0] + '.yaml' out_name = os.path.join(self.output_directory, 'cantera', file_name) - if 'surfaceFile' in kwargs: + if 'surface_file' in kwargs: out_name = out_name.replace('-gas.', '.') cantera_dir = os.path.dirname(out_name) try: @@ -1687,14 +1687,14 @@ def generate_cantera_files(self, chemkin_file, **kwargs): raise if os.path.exists(out_name): os.remove(out_name) - parser = ck2cti.Parser() + parser = ck2yaml.Parser() try: - parser.convertMech(chemkin_file, transportFile=transport_file, outName=out_name, quiet=True, permissive=True, + parser.convert_mech(chemkin_file, transport_file=transport_file, out_name=out_name, quiet=True, permissive=True, **kwargs) - except ck2cti.InputParseError: + except ck2yaml.InputError: logging.exception("Error converting to Cantera format.") logging.info("Trying again without transport data file.") - parser.convertMech(chemkin_file, outName=out_name, quiet=True, permissive=True, **kwargs) + parser.convert_mech(chemkin_file, out_name=out_name, quiet=True, permissive=True, **kwargs) def initialize_reaction_threshold_and_react_flags(self): num_core_species = len(self.reaction_model.core.species) diff --git a/rmgpy/rmg/mainTest.py b/rmgpy/rmg/mainTest.py index aa1f315c94..2127affb57 100644 --- a/rmgpy/rmg/mainTest.py +++ b/rmgpy/rmg/mainTest.py @@ -180,7 +180,7 @@ def test_make_cantera_input_file(self): outName = os.path.join(self.rmg.output_directory, 'cantera') files = os.listdir(outName) for f in files: - if '.cti' in f: + if '.yaml' in f: try: ct.Solution(os.path.join(outName, f)) except: @@ -411,7 +411,7 @@ def setUp(self): THERM ALL 300.000 1000.000 5000.000 -ethane(1) H 6 C 2 G100.000 5000.000 954.52 1 +ethane(1) H 6 C 2 G100.000 5000.000 954.52 1 4.58987205E+00 1.41507042E-02-4.75958084E-06 8.60284590E-10-6.21708569E-14 2 -1.27217823E+04-3.61762003E+00 3.78032308E+00-3.24248354E-03 5.52375224E-05 3 -6.38573917E-08 2.28633835E-11-1.16203404E+04 5.21037799E+00 4 @@ -492,7 +492,7 @@ def test_chemkin_to_cantera_conversion(self): Tests that good and bad chemkin files raise proper exceptions """ - from cantera.ck2cti import InputParseError + from cantera.ck2yaml import InputError for ck_input, works in self.chemkin_files.items(): os.chdir(originalPath) @@ -510,7 +510,7 @@ def test_chemkin_to_cantera_conversion(self): if works: self.rmg.generate_cantera_files(os.path.join(os.getcwd(), 'chem001.inp')) else: - with self.assertRaises(InputParseError): + with self.assertRaises(InputError): self.rmg.generate_cantera_files(os.path.join(os.getcwd(), 'chem001.inp')) # clean up diff --git a/rmgpy/tools/canteramodel.py b/rmgpy/tools/canteramodel.py index be6ea784b1..4684ce4cdb 100644 --- a/rmgpy/tools/canteramodel.py +++ b/rmgpy/tools/canteramodel.py @@ -203,7 +203,7 @@ def __init__(self, species_list=None, reaction_list=None, canteraFile='', output `species_list`: list of RMG species objects `reaction_list`: list of RMG reaction objects `reaction_map`: dict mapping the RMG reaction index within the `reaction_list` to cantera model reaction(s) indices - `canteraFile` path of the chem.cti file associated with this job + `canteraFile` path of the chem.yaml file associated with this job `conditions`: a list of `CanteraCondition` objects `sensitive_species`: a list of RMG species objects for conductng sensitivity analysis on `thermo_SA`: a boolean indicating whether or not to run thermo SA. By default, if sensitive_species is given, @@ -289,19 +289,19 @@ def refresh_model(self): def load_chemkin_model(self, chemkin_file, transport_file=None, **kwargs): """ - Convert a chemkin mechanism chem.inp file to a cantera mechanism file chem.cti + Convert a chemkin mechanism chem.inp file to a cantera mechanism file chem.yaml and save it in the output_directory Then load it into self.model """ - from cantera import ck2cti + from cantera import ck2yaml base = os.path.basename(chemkin_file) base_name = os.path.splitext(base)[0] - out_name = os.path.join(self.output_directory, base_name + ".cti") + out_name = os.path.join(self.output_directory, base_name + ".yaml") if os.path.exists(out_name): os.remove(out_name) - parser = ck2cti.Parser() - parser.convertMech(chemkin_file, transportFile=transport_file, outName=out_name, **kwargs) + parser = ck2yaml.Parser() + parser.convert_mech(chemkin_file, transport_file=transport_file, out_name=out_name, **kwargs) self.model = ct.Solution(out_name) def modify_reaction_kinetics(self, rmg_reaction_index, rmg_reaction): @@ -734,41 +734,40 @@ def check_equivalent_falloff(fall1, fall2): assert ct_rxn1.reactants == ct_rxn2.reactants, "Same reactants" assert ct_rxn1.products == ct_rxn2.products, "Same products" - if isinstance(ct_rxn1, ct.ElementaryReaction): - assert ct_rxn1.allow_negative_pre_exponential_factor == ct_rxn2.allow_negative_pre_exponential_factor, \ - "Same allow_negative_pre_exponential_factor attribute" - if ct_rxn1.rate or ct_rxn2.rate: - check_equivalent_arrhenius(ct_rxn1.rate, ct_rxn2.rate) - - elif isinstance(ct_rxn1, ct.PlogReaction): - if ct_rxn1.rates or ct_rxn2.rates: - assert len(ct_rxn1.rates) == len(ct_rxn2.rates), "Same number of rates in PLOG reaction" - - for i in range(len(ct_rxn1.rates)): - P1, arr1 = ct_rxn1.rates[i] - P2, arr2 = ct_rxn2.rates[i] - assert check_nearly_equal(P1, P2, dE), "Similar pressures for PLOG rates" - check_equivalent_arrhenius(arr1, arr2) - - elif isinstance(ct_rxn1, ct.ChebyshevReaction): - assert ct_rxn1.Pmax == ct_rxn2.Pmax, "Same Pmax for Chebyshev reaction" - assert ct_rxn1.Pmin == ct_rxn2.Pmin, "Same Pmin for Chebyshev reaction" - assert ct_rxn1.Tmax == ct_rxn2.Tmax, "Same Tmax for Chebyshev reaction" - assert ct_rxn1.Tmin == ct_rxn2.Tmin, "Same Tmin for Chebyshev reaction" - assert ct_rxn1.nPressure == ct_rxn2.nPressure, "Same number of pressure interpolations" - assert ct_rxn1.nTemperature == ct_rxn2.nTemperature, "Same number of temperature interpolations" - for i in range(ct_rxn1.coeffs.shape[0]): - for j in range(ct_rxn1.coeffs.shape[1]): - assert check_nearly_equal(ct_rxn1.coeffs[i, j], ct_rxn2.coeffs[i, j], dE), \ - "Similar Chebyshev coefficients" + if isinstance(ct_rxn1, ct.Reaction): + # may not mean it is arrhenius, need to check if it is troe, + if isinstance(ct_rxn1.rate, ct.ArrheniusRate): + assert ct_rxn1.allow_negative_pre_exponential_factor == ct_rxn2.allow_negative_pre_exponential_factor, \ + "Same allow_negative_pre_exponential_factor attribute" + if ct_rxn1.rate or ct_rxn2.rate: + check_equivalent_arrhenius(ct_rxn1.rate, ct_rxn2.rate) + elif isinstance(ct_rxn1.rate, ct.PlogRate): + if ct_rxn1.rate.rates or ct_rxn2.rate.rates: + assert len(ct_rxn1.rates) == len(ct_rxn2.rates), "Same number of rates in PLOG reaction" + + for i in range(len(ct_rxn1.rate.rates)): + P1, arr1 = ct_rxn1.rate.rates[i] + P2, arr2 = ct_rxn2.rate.rates[i] + assert check_nearly_equal(P1, P2, dE), "Similar pressures for PLOG rates" + check_equivalent_arrhenius(arr1, arr2) + + elif isinstance(ct_rxn1.rate, ct.ChebyshevRate): + assert ct_rxn1.rate.pressure_range == ct_rxn2.rate.pressure_range, "Same Prange for Chebyshev reaction" + assert ct_rxn1.rate.temperature_range == ct_rxn2.rate.temperature_range, "Same Trange for Chebyshev reaction" + assert ct_rxn1.rate.n_pressure == ct_rxn2.rate.n_pressure, "Same number of pressure interpolations" + assert ct_rxn1.rate.n_temperature == ct_rxn2.rate.n_temperature , "Same number of temperature interpolations" + for i in range(ct_rxn1.rate.data.shape[0]): + for j in range(ct_rxn1.rate.data.shape[1]): + assert check_nearly_equal(ct_rxn1.rate.data[i, j], ct_rxn2.rate.data[i, j], dE), \ + "Similar Chebyshev coefficients" elif isinstance(ct_rxn1, ct.ThreeBodyReaction): assert ct_rxn1.default_efficiency == ct_rxn2.default_efficiency, "Same default efficiency" - assert ct_rxn1.efficiencies == ct_rxn2.efficiencies, "Same efficienciess" + assert ct_rxn1.efficiencies == ct_rxn2.efficiencies, "Same efficiencies" elif isinstance(ct_rxn1, ct.FalloffReaction): assert ct_rxn1.default_efficiency == ct_rxn2.default_efficiency, "Same default efficiency" - assert ct_rxn1.efficiencies == ct_rxn2.efficiencies, "Same efficienciess" + assert ct_rxn1.efficiencies == ct_rxn2.efficiencies, "Same efficiencies" if ct_rxn1.falloff or ct_rxn2.falloff: check_equivalent_falloff(ct_rxn1.falloff, ct_rxn2.falloff) if ct_rxn1.high_rate or ct_rxn2.high_rate: diff --git a/rmgpy/tools/data/various_kinetics/chem_annotated.inp b/rmgpy/tools/data/various_kinetics/chem_annotated.inp index fb9584ca44..54acce1ae8 100644 --- a/rmgpy/tools/data/various_kinetics/chem_annotated.inp +++ b/rmgpy/tools/data/various_kinetics/chem_annotated.inp @@ -184,7 +184,6 @@ HO2(5)+CH3(13)=O2(6)+CH4(15) 1.000e+00 0.000 0.000 PLOG/ 80.000 2.500e+05 1.900 0.000 / PLOG/ 100.000 1.900e+05 1.940 0.000 / PLOG/ 650.000 7.000e+05 1.700 0.298 / - PLOG/ 2000.000 0.000e+00 0.000 0.000 / DUPLICATE