Skip to content

ArrheniusBM fitting procedure #2264

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
davidfarinajr wants to merge 3 commits into from
Closed

ArrheniusBM fitting procedure #2264

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
dcbd55a
revised arrbm `get_activation_energy` method
davidfarinajr Jan 17, 2022
833b973
revised arrbm fitting to reactions procedure
davidfarinajr Jan 17, 2022
e80a8ff
modifed arrbm `fit_to_data` unit test
davidfarinajr Jan 18, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
157 changes: 79 additions & 78 deletions rmgpy/kinetics/arrhenius.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
"""
Expand All @@ -597,91 +603,86 @@ 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)']
order = len(rxns[0].reactants)
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

Expand Down
59 changes: 54 additions & 5 deletions rmgpy/kinetics/arrheniusTest.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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


################################################################################
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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):
"""
Expand Down Expand Up @@ -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))