From dcbd55a80fa0ffff6f1eeb1461f0ed247ad139e6 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Mon, 17 Jan 2022 13:48:15 -0500 Subject: [PATCH 1/3] revised arrbm `get_activation_energy` method We previously did not account for E0 < 0. Here, for E0 < 0, we set Ea to E0 for exothermic rxns and dHrxn for endothermic. We also incorporate aspects of Reaction `fix_barrier_height` method to set Ea to 0 for exothermic rxns and Ea < 0, and set Ea to dHrxn if Ea < dHrxn for endothermic rxns. --- rmgpy/kinetics/arrhenius.pyx | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx index 8c4856c560..f2bd4bfc22 100644 --- a/rmgpy/kinetics/arrhenius.pyx +++ b/rmgpy/kinetics/arrhenius.pyx @@ -562,16 +562,22 @@ cdef class ArrheniusBM(KineticsModel): Return the activation energy in J/mol corresponding to the given enthalpy of reaction `dHrxn` in J/mol. """ - cdef double w0, E0 + cdef double w0, E0, Ea E0 = self._E0.value_si - if dHrxn < -4 * self._E0.value_si: - return 0.0 - elif dHrxn > 4 * self._E0.value_si: - return dHrxn + w0 = self._w0.value_si + if E0 < 0: + if dHrxn > 0: + Ea = dHrxn + else: + Ea = E0 else: - w0 = self._w0.value_si Vp = 2 * w0 * (2 * w0 + 2 * E0) / (2 * w0 - 2 * E0) - return (w0 + dHrxn / 2.0) * (Vp - 2 * w0 + dHrxn) ** 2 / (Vp ** 2 - (2 * w0) ** 2 + dHrxn ** 2) + Ea = (w0 + dHrxn / 2.0) * (Vp - 2 * w0 + dHrxn) ** 2 / (Vp ** 2 - (2 * w0) ** 2 + dHrxn ** 2) + if (dHrxn < 0.0 and Ea < 0.0) or (dHrxn < -4 * E0): + Ea = 0.0 + elif (dHrxn > 0.0 and Ea < dHrxn) or (dHrxn > 4 * E0): + Ea = dHrxn + return Ea cpdef Arrhenius to_arrhenius(self, double dHrxn): """ From 833b9735de9af859b1c7861e834c73495ee61d29 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Mon, 17 Jan 2022 13:55:32 -0500 Subject: [PATCH 2/3] revised arrbm fitting to reactions procedure Use the `get_activation_energy` method in the objective function. For intial guess of params, use average of A and n, and use BEP for guess of Ea. --- rmgpy/kinetics/arrhenius.pyx | 137 +++++++++++++++++------------------ 1 file changed, 66 insertions(+), 71 deletions(-) diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx index f2bd4bfc22..b2f3dff8b4 100644 --- a/rmgpy/kinetics/arrhenius.pyx +++ b/rmgpy/kinetics/arrhenius.pyx @@ -603,82 +603,78 @@ cdef class ArrheniusBM(KineticsModel): assert w0 is not None or recipe is not None, 'either w0 or recipe must be specified' if Ts is None: - Ts = [300.0, 500.0, 600.0, 700.0, 800.0, 900.0, 1000.0, 1100.0, 1200.0, 1500.0, 2000.0] + Ts = np.array([300.0, 500.0, 600.0, 700.0, 800.0, 900.0, 1000.0, 1100.0, 1200.0, 1500.0, 2000.0]) if w0 is None: #estimate w0 w0s = get_w0s(recipe, rxns) w0 = sum(w0s) / len(w0s) - - if len(rxns) == 1: - T = 1000.0 - rxn = rxns[0] - dHrxn = rxn.get_enthalpy_of_reaction(T) - A = rxn.kinetics.A.value_si - n = rxn.kinetics.n.value_si - Ea = rxn.kinetics.Ea.value_si + self.w0 = (w0 * 0.001, 'kJ/mol') + + def kfcn(xs, lnA, n, E0): + T = xs[:,0] + dHrxn = xs[:,1] + self.E0 = (E0, 'J/mol') + Ea = np.array([self.get_activation_energy(dHrxn[i]) for i in range(len(dHrxn))]) + return lnA + np.log(T ** n * np.exp(-Ea / (constants.R * T))) - def kfcn(E0): - Vp = 2 * w0 * (2 * w0 + 2 * E0) / (2 * w0 - 2 * E0) - out = Ea - (w0 + dHrxn / 2.0) * (Vp - 2 * w0 + dHrxn) * (Vp - 2 * w0 + dHrxn) / (Vp * Vp - (2 * w0) * (2 * w0) + dHrxn * dHrxn) - return out - - if abs(dHrxn) > 4 * w0 / 10.0: - E0 = w0 / 10.0 - else: - E0 = fsolve(kfcn, w0 / 10.0)[0] - - self.Tmin = rxn.kinetics.Tmin - self.Tmax = rxn.kinetics.Tmax - self.comment = 'Fitted to {0} reaction at temperature: {1} K'.format(len(rxns), T) - else: - # define optimization function - def kfcn(xs, lnA, n, E0): - T = xs[:,0] - dHrxn = xs[:,1] - Vp = 2 * w0 * (2 * w0 + 2 * E0) / (2 * w0 - 2 * E0) - Ea = (w0 + dHrxn / 2.0) * (Vp - 2 * w0 + dHrxn) * (Vp - 2 * w0 + dHrxn) / (Vp * Vp - (2 * w0) * (2 * w0) + dHrxn * dHrxn) - Ea = np.where(dHrxn< -4.0*E0, 0.0, Ea) - Ea = np.where(dHrxn > 4.0*E0, dHrxn, Ea) - return lnA + np.log(T ** n * np.exp(-Ea / (8.314 * T))) - - # get (T,dHrxn(T)) -> (Ln(k) mappings - xdata = [] - ydata = [] - sigmas = [] - for rxn in rxns: - # approximately correct the overall uncertainties to std deviations - s = rank_accuracy_map[rxn.rank].value_si/2.0 - for T in Ts: - xdata.append([T, rxn.get_enthalpy_of_reaction(T)]) - ydata.append(np.log(rxn.get_rate_coefficient(T))) - - sigmas.append(s / (8.314 * T)) - - xdata = np.array(xdata) - ydata = np.array(ydata) - - # fit parameters - boo = True - xtol = 1e-8 - ftol = 1e-8 - while boo: - boo = False - try: - params = curve_fit(kfcn, xdata, ydata, sigma=sigmas, p0=[1.0, 1.0, w0 / 10.0], xtol=xtol, ftol=ftol) - except RuntimeError: - if xtol < 1.0: - boo = True - xtol *= 10.0 - ftol *= 10.0 - else: - raise ValueError("Could not fit BM arrhenius to reactions with xtol<1.0") + # get (T,dHrxn(T)) -> (Ln(k) mappings + xdata = [] + ydata = [] + sigmas = [] + E0 = 0.0 + lnA = 0.0 + n = 0.0 + for rxn in rxns: + # approximately correct the overall uncertainties to std deviations + s = rank_accuracy_map[rxn.rank].value_si/2.0 + # Use BEP with alpha = 0.25 for inital guess of E0 + E0 += rxn.kinetics._Ea.value_si - 0.25 * rxn.get_enthalpy_of_reaction(298) + lnA += np.log(rxn.kinetics.A.value_si) + n += rxn.kinetics.n.value_si + for T in Ts: + xdata.append([T, rxn.get_enthalpy_of_reaction(298)]) + ydata.append(np.log(rxn.get_rate_coefficient(T))) + sigmas.append(s / (constants.R * T)) + # Use the average of the E0s as intial guess + E0 /= len(rxns) + lnA /= len(rxns) + n /= len(rxns) + E0 = min(E0, w0) + if E0 < 0: + E0 = w0 / 100.0 + + xdata = np.array(xdata) + ydata = np.array(ydata) + + # fit parameters + boo = True + xtol = 1e-8 + ftol = 1e-8 + attempts = 0 + while boo: + boo = False + try: + params = curve_fit(kfcn, xdata, ydata, sigma=sigmas, p0=[lnA, n, E0], xtol=xtol, ftol=ftol) + lnA, n, E0 = params[0].tolist() + if abs(E0/self.w0.value_si) > 1 and attempts < 5: + boo = True + if attempts > 0: + self.w0.value_si *= 1.25 + attempts += 1 + E0 = self.w0.value_si / 10.0 + except RuntimeError: + if xtol < 1.0: + boo = True + xtol *= 10.0 + ftol *= 10.0 + else: + raise ValueError("Could not fit BM arrhenius to reactions with xtol<1.0") - lnA, n, E0 = params[0].tolist() - A = np.exp(lnA) + A = np.exp(lnA) - self.Tmin = (np.min(Ts), "K") - self.Tmax = (np.max(Ts), "K") - self.comment = 'Fitted to {0} reactions at temperatures: {1}'.format(len(rxns), Ts) + self.Tmin = (np.min(Ts), "K") + self.Tmax = (np.max(Ts), "K") + self.comment = 'Fitted to {0} reactions at temperatures: {1}'.format(len(rxns), Ts) # fill in parameters A_units = ['', 's^-1', 'm^3/(mol*s)', 'm^6/(mol^2*s)'] @@ -686,8 +682,7 @@ cdef class ArrheniusBM(KineticsModel): self.A = (A, A_units[order]) self.n = n - self.w0 = (w0, 'J/mol') - self.E0 = (E0, 'J/mol') + self.E0 = (E0 * 0.001, 'kJ/mol') return self From e80a8ff11b4adb293603afd449bb8185dc8682e0 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Tue, 18 Jan 2022 11:39:04 -0500 Subject: [PATCH 3/3] modifed arrbm `fit_to_data` unit test Test now compares the fitted rate to the rate it was trained on to make sure they agree. --- rmgpy/kinetics/arrheniusTest.py | 59 ++++++++++++++++++++++++++++++--- 1 file changed, 54 insertions(+), 5 deletions(-) diff --git a/rmgpy/kinetics/arrheniusTest.py b/rmgpy/kinetics/arrheniusTest.py index f814ac49b3..031d080d93 100644 --- a/rmgpy/kinetics/arrheniusTest.py +++ b/rmgpy/kinetics/arrheniusTest.py @@ -41,7 +41,7 @@ from rmgpy.molecule.molecule import Molecule from rmgpy.reaction import Reaction from rmgpy.species import Species -from rmgpy.thermo import NASA, NASAPolynomial +from rmgpy.thermo import NASA, NASAPolynomial, ThermoData ################################################################################ @@ -460,6 +460,42 @@ def setUp(self): comment="""Thermo library: Spiekermann_refining_elementary_reactions""" ) + CF2 = Species().from_adjacency_list( + """ + 1 F u0 p3 c0 {2,S} + 2 C u0 p1 c0 {1,S} {3,S} + 3 F u0 p3 c0 {2,S} + """ + ) + CF2.thermo = NASA( + polynomials=[ + NASAPolynomial(coeffs=[2.28591,0.0107608,-1.05382e-05,4.89881e-09,-8.86384e-13,-24340.7,13.1348], Tmin=(298,'K'), Tmax=(1300,'K')), + NASAPolynomial(coeffs=[5.33121,0.00197748,-9.60248e-07,2.10704e-10,-1.5954e-14,-25190.9,-2.56367], Tmin=(1300,'K'), Tmax=(3000,'K')) + ], + Tmin=(298,'K'), Tmax=(3000,'K'), Cp0=(33.2579,'J/mol/K'), CpInf=(58.2013,'J/mol/K'), + comment="""Thermo library: halogens""" + ) + C2H6 = Species(smiles="CC") + C2H6.thermo = ThermoData( + Tdata = ([300,400,500,600,800,1000,1500],'K'), + Cpdata = ([12.565,15.512,18.421,21.059,25.487,28.964,34.591],'cal/(mol*K)','+|-',[0.8,1.1,1.3,1.4,1.5,1.5,1.2]), + H298 = (-20.028,'kcal/mol','+|-',0.1), + S298 = (54.726,'cal/(mol*K)','+|-',0.6), + comment="""Thermo library: DFT_QCI_thermo""" + ) + CH3CF2CH3 = Species(smiles="CC(F)(F)C") + CH3CF2CH3.thermo = NASA( + polynomials = [ + NASAPolynomial(coeffs=[3.89769,0.00706735,0.000140168,-3.37628e-07,2.51812e-10,-68682.1,8.74321], Tmin=(10,'K'), Tmax=(436.522,'K')), + NASAPolynomial(coeffs=[2.78849,0.0356982,-2.16715e-05,6.45057e-09,-7.47989e-13,-68761.2,11.1597], Tmin=(436.522,'K'), Tmax=(3000,'K')), + ], + Tmin = (10,'K'), Tmax = (3000,'K'), Cp0 = (33.2579,'J/(mol*K)'), CpInf = (249.434,'J/(mol*K)'), + comment="""Thermo library: CHOF_G4""" + ) + kinetics = Arrhenius(A=(0.222791,'cm^3/(mol*s)'), n=3.59921, Ea=(320.496,'kJ/mol'), T0=(1,'K'), Tmin=(298,'K'), Tmax=(2500,'K'), comment="""Training Rxn 54 for 1,2_Insertion_carbene""") + self.reaction = Reaction(reactants=[CF2,C2H6],products=[CH3CF2CH3],kinetics=kinetics) + self.reaction_w0 = 519000 # J/mol + def test_a_factor(self): """ Test that the ArrheniusBM A property was properly set. @@ -510,17 +546,25 @@ def test_fit_to_data(self): """ Test the ArrheniusBM.fit_to_data() method. """ + Tdata = np.array([300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500]) + reactant = Molecule(smiles=self.rsmi) product = Molecule(smiles=self.psmi) reaction = Reaction(reactants=[Species(molecule=[reactant], thermo=self.r_thermo,)], products=[Species(molecule=[product], thermo=self.p_thermo)], kinetics=self.arrhenius, ) - + kdata = np.array([reaction.kinetics.get_rate_coefficient(T) for T in Tdata]) arrhenius_bm = ArrheniusBM().fit_to_reactions([reaction], w0=self.w0) - self.assertAlmostEqual(arrhenius_bm.A.value_si, self.arrhenius_bm.A.value_si, delta=1.5e1) - self.assertAlmostEqual(arrhenius_bm.n.value_si, self.arrhenius_bm.n.value_si, 1, 4) - self.assertAlmostEqual(arrhenius_bm.E0.value_si, self.arrhenius_bm.E0.value_si, 1) + arrhenius = arrhenius_bm.to_arrhenius(reaction.get_enthalpy_of_reaction(298)) + for T, k in zip(Tdata, kdata): + self.assertAlmostEqual(k, arrhenius.get_rate_coefficient(T), delta=1e-6 * k) + + arrhenius_bm = ArrheniusBM().fit_to_reactions([self.reaction], w0=self.reaction_w0) + arrhenius = arrhenius_bm.to_arrhenius(self.reaction.get_enthalpy_of_reaction(298)) + kdata = np.array([self.reaction.kinetics.get_rate_coefficient(T) for T in Tdata]) + for T, k in zip(Tdata, kdata): + self.assertAlmostEqual(k, arrhenius.get_rate_coefficient(T), delta=1e-6 * k) def test_get_activation_energy(self): """ @@ -1195,3 +1239,8 @@ def test_generate_reverse_rate_coefficient(self): Arrhenius(A=(-2.1e+11,"cm^3/(mol*s)"), n=0.618, Ea=(30251,"cal/mol"), T0=(1,"K"))]) ]), duplicate=True) test_reaction.generate_reverse_rate_coefficient() + +################################################################################ + +if __name__ == '__main__': + unittest.main(testRunner=unittest.TextTestRunner(verbosity=2))