diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index eb7223836e9..441e57fce13 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -3678,19 +3678,22 @@ def cross_validate(self, folds=5, template_rxn_map=None, test_rxn_inds=None, T=1 uncertainties[rxn] = self.rules.entries[entry.label][0].data.uncertainty - L = list(set(template_rxn_map[entry.label]) - set(rxns_test)) if L != []: if isinstance(L[0].kinetics, Arrhenius): kinetics = ArrheniusBM().fit_to_reactions(L, recipe=self.forward_recipe.actions) - if kinetics.E0.value_si < 0.0 or len(L) == 1: + # These conditions (eg. n>5) for rejecting the BM fit should correspond + # to those in _marke_rule, and those just below: + if len(L) == 1 or kinetics.E0.value_si < 0.0 or abs(kinetics.n.value_si) > 5.0: kinetics = average_kinetics([r.kinetics for r in L]) else: kinetics = kinetics.to_arrhenius(rxn.get_enthalpy_of_reaction(298.)) else: kinetics = ArrheniusChargeTransferBM().fit_to_reactions(L, recipe=self.forward_recipe.actions) - if kinetics.E0.value_si < 0.0 or len(L) == 1: + # These conditions (eg. n>5) for rejecting the BM fit should correspond + # to those in _marke_rule, and those just above: + if len(L) == 1 or kinetics.E0.value_si < 0.0 or abs(kinetics.n.value_si) > 5.0: kinetics = average_kinetics([r.kinetics for r in L]) else: kinetics = kinetics.to_arrhenius_charge_transfer(rxn.get_enthalpy_of_reaction(298.)) @@ -4596,79 +4599,93 @@ def _make_rule(rr): for i, rxn in enumerate(rxns): rxn.rank = ranks[i] rxns = np.array(rxns) - rs = np.array([r for r in rxns if type(r.kinetics) != KineticsModel]) + rs = np.array([r for r in rxns if type(r.kinetics) is not KineticsModel]) # KineticsModel is the base class with no data. n = len(rs) - if n > 0 and isinstance(rs[0].kinetics, Marcus): + if n == 0: + return None + + if isinstance(rs[0].kinetics, Marcus): kin = average_kinetics([r.kinetics for r in rs]) return kin data_mean = np.mean(np.log([r.kinetics.get_rate_coefficient(Tref) for r in rs])) - if n > 0: - if isinstance(rs[0].kinetics, Arrhenius): - arr = ArrheniusBM - else: - arr = ArrheniusChargeTransferBM - if n > 1: - kin = arr().fit_to_reactions(rs, recipe=recipe) - if n == 1 or kin.E0.value_si < 0.0: - kin = average_kinetics([r.kinetics for r in rs]) - #kin.comment = "Only one reaction or Arrhenius BM fit bad. Instead averaged from {} reactions.".format(n) - if n == 1: - kin.uncertainty = RateUncertainty(mu=0.0, var=(np.log(fmax) / 2.0) ** 2, N=1, Tref=Tref, data_mean=data_mean, correlation=label) + + if isinstance(rs[0].kinetics, Arrhenius): + arr = ArrheniusBM + else: + arr = ArrheniusChargeTransferBM + if n > 1: + kin = arr().fit_to_reactions(rs, recipe=recipe) + + # If you update the following conditions, also update the cross_validate function. + if n == 1 or kin.E0.value_si < 0.0 or abs(kin.n.value_si) > 5.0: + if n > 1: #we have an ArrheniusBM + if kin.E0.value_si < 0.0: + reason = "E0<0" + elif abs(kin.n.value_si) > 5.0: + reason = "abs(n)>5" else: + reason = "?" + + # still run it through the averaging function when n=1 to standardize the units and run checks + kin = average_kinetics([r.kinetics for r in rs]) + + if n == 1: + kin.uncertainty = RateUncertainty(mu=0.0, var=(np.log(fmax) / 2.0) ** 2, N=1, Tref=Tref, data_mean=data_mean, correlation=label) + else: + kin.comment = f"Blowers-Masel fit was bad ({reason}) so instead averaged from {n} reactions." + dlnks = np.array([ + np.log( + average_kinetics([r.kinetics for r in np.delete(rs, i)]).get_rate_coefficient(T=Tref) / rxn.get_rate_coefficient(T=Tref) + ) for i, rxn in enumerate(rs) + ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref + varis = (np.array([rank_accuracy_map[rxn.rank].value_si for rxn in rs]) / (2.0 * constants.R * Tref)) ** 2 + # weighted average calculations + ws = 1.0 / varis + V1 = ws.sum() + V2 = (ws ** 2).sum() + mu = np.dot(ws, dlnks) / V1 + s = np.sqrt(np.dot(ws, (dlnks - mu) ** 2) / (V1 - V2 / V1)) + kin.uncertainty = RateUncertainty(mu=mu, var=s ** 2, N=n, Tref=Tref, data_mean=data_mean, correlation=label) + else: + if n == 1: + kin.uncertainty = RateUncertainty(mu=0.0, var=(np.log(fmax) / 2.0) ** 2, N=1, Tref=Tref, data_mean=data_mean, correlation=label) + else: + if isinstance(rs[0].kinetics, Arrhenius): dlnks = np.array([ np.log( - average_kinetics([r.kinetics for r in rs[list(set(range(len(rs))) - {i})]]).get_rate_coefficient(T=Tref) / rxn.get_rate_coefficient(T=Tref) - ) for i, rxn in enumerate(rs) - ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref - varis = (np.array([rank_accuracy_map[rxn.rank].value_si for rxn in rs]) / (2.0 * 8.314 * Tref)) ** 2 - # weighted average calculations - ws = 1.0 / varis - V1 = ws.sum() - V2 = (ws ** 2).sum() - mu = np.dot(ws, dlnks) / V1 - s = np.sqrt(np.dot(ws, (dlnks - mu) ** 2) / (V1 - V2 / V1)) - kin.uncertainty = RateUncertainty(mu=mu, var=s ** 2, N=n, Tref=Tref, data_mean=data_mean, correlation=label) - else: - if n == 1: - kin.uncertainty = RateUncertainty(mu=0.0, var=(np.log(fmax) / 2.0) ** 2, N=1, Tref=Tref, data_mean=data_mean, correlation=label) + arr().fit_to_reactions(np.delete(rs, i), recipe=recipe) + .to_arrhenius(rxn.get_enthalpy_of_reaction(Tref)) + .get_rate_coefficient(T=Tref) / rxn.get_rate_coefficient(T=Tref) + ) for i, rxn in enumerate(rs) + ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref else: - if isinstance(rs[0].kinetics, Arrhenius): - dlnks = np.array([ - np.log( - arr().fit_to_reactions(rs[list(set(range(len(rs))) - {i})], recipe=recipe) - .to_arrhenius(rxn.get_enthalpy_of_reaction(Tref)) - .get_rate_coefficient(T=Tref) / rxn.get_rate_coefficient(T=Tref) - ) for i, rxn in enumerate(rs) - ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref - else: - dlnks = np.array([ - np.log( - arr().fit_to_reactions(rs[list(set(range(len(rs))) - {i})], recipe=recipe) - .to_arrhenius_charge_transfer(rxn.get_enthalpy_of_reaction(Tref)) - .get_rate_coefficient(T=Tref) / rxn.get_rate_coefficient(T=Tref) - ) for i, rxn in enumerate(rs) - ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref - varis = (np.array([rank_accuracy_map[rxn.rank].value_si for rxn in rs]) / (2.0 * 8.314 * Tref)) ** 2 - # weighted average calculations - ws = 1.0 / varis - V1 = ws.sum() - V2 = (ws ** 2).sum() - mu = np.dot(ws, dlnks) / V1 - s = np.sqrt(np.dot(ws, (dlnks - mu) ** 2) / (V1 - V2 / V1)) - kin.uncertainty = RateUncertainty(mu=mu, var=s ** 2, N=n, Tref=Tref, data_mean=data_mean, correlation=label) - - #site solute parameters - site_datas = [get_site_solute_data(rxn) for rxn in rxns] - site_datas = [sdata for sdata in site_datas if sdata is not None] - if len(site_datas) > 0: - site_data = SoluteTSData() - for sdata in site_datas: - site_data += sdata - site_data = site_data * (1.0/len(site_datas)) - kin.solute = site_data - return kin - else: - return None + dlnks = np.array([ + np.log( + arr().fit_to_reactions(np.delete(rs, i), recipe=recipe) + .to_arrhenius_charge_transfer(rxn.get_enthalpy_of_reaction(Tref)) + .get_rate_coefficient(T=Tref) / rxn.get_rate_coefficient(T=Tref) + ) for i, rxn in enumerate(rs) + ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref + varis = (np.array([rank_accuracy_map[rxn.rank].value_si for rxn in rs]) / (2.0 * constants.R * Tref)) ** 2 + # weighted average calculations + ws = 1.0 / varis + V1 = ws.sum() + V2 = (ws ** 2).sum() + mu = np.dot(ws, dlnks) / V1 + s = np.sqrt(np.dot(ws, (dlnks - mu) ** 2) / (V1 - V2 / V1)) + kin.uncertainty = RateUncertainty(mu=mu, var=s ** 2, N=n, Tref=Tref, data_mean=data_mean, correlation=label) + + #site solute parameters + site_datas = [get_site_solute_data(rxn) for rxn in rxns] + site_datas = [sdata for sdata in site_datas if sdata is not None] + if len(site_datas) > 0: + site_data = SoluteTSData() + for sdata in site_datas: + site_data += sdata + site_data = site_data * (1.0/len(site_datas)) + kin.solute = site_data + return kin + def _spawn_tree_process(family, template_rxn_map, obj, T, nprocs, depth, min_splitable_entry_num, min_rxns_to_spawn, extension_iter_max, extension_iter_item_cap): @@ -4772,7 +4789,7 @@ def average_kinetics(kinetics_list): beta=(beta,"1/m"), wr=(wr * 0.001, "kJ/mol"), wp=(wp * 0.001, "kJ/mol"), - comment="Averaged from {} reactions.".format(len(kinetics_list)), + comment=f"Averaged from {len(kinetics_list)} rate expressions.", ) elif isinstance(kinetics, SurfaceChargeTransfer): averaged_kinetics = SurfaceChargeTransfer( @@ -4782,6 +4799,7 @@ def average_kinetics(kinetics_list): alpha=alpha, V0=(V0,'V'), Ea=(Ea * 0.001, "kJ/mol"), + comment=f"Averaged from {len(kinetics_list)} rate expressions.", ) elif isinstance(kinetics, ArrheniusChargeTransfer): averaged_kinetics = ArrheniusChargeTransfer( @@ -4791,6 +4809,7 @@ def average_kinetics(kinetics_list): alpha=alpha, V0=(V0,'V'), Ea=(Ea * 0.001, "kJ/mol"), + comment=f"Averaged from {len(kinetics_list)} rate expressions.", ) else: averaged_kinetics = Arrhenius( diff --git a/rmgpy/display.py b/rmgpy/display.py index 3b6c751a10b..cf9a40d2d28 100644 --- a/rmgpy/display.py +++ b/rmgpy/display.py @@ -29,18 +29,18 @@ """ This module contains only a display() function, which, if you are running in -the IPython --pylab mode, will render things inline in pretty SVG or PNG graphics. -If you are NOT running in IPython --pylab mode, it will do nothing. +a Jupyter notebook or IPython session, will render things inline in pretty SVG or PNG graphics. +If you are NOT running in such an environment, it will do nothing. """ try: - from IPython.core.interactiveshell import InteractiveShell + from IPython import get_ipython + _ipython = get_ipython() except ImportError: + _ipython = None + +if _ipython is not None: + from IPython.display import display +else: def display(obj): pass -else: - if InteractiveShell.initialized(): - from IPython.core.display import display - else: - def display(obj): - pass diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx index ea9e462d4f8..6252b6ce950 100644 --- a/rmgpy/kinetics/arrhenius.pyx +++ b/rmgpy/kinetics/arrhenius.pyx @@ -573,16 +573,20 @@ cdef class ArrheniusBM(KineticsModel): Return the activation energy in J/mol corresponding to the given enthalpy of reaction `dHrxn` in J/mol, evaluated at 298 K. """ - 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: + # Negative E0 is unphysical, but could be encountered during optimization. + Ea = E0 + max(dHrxn, 0.0) 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) + Vp = 2 * w0 * (w0 + E0) / (w0 - E0) + Ea = (w0 + dHrxn / 2.0) * (Vp - 2 * w0 + dHrxn) * (Vp - 2 * w0 + dHrxn) / (Vp * Vp - (4 * w0 * w0) + dHrxn * dHrxn) + if (Ea < 0.0) or (dHrxn < -4 * E0): + Ea = 0.0 + elif (Ea < dHrxn) or (dHrxn > 4 * E0): + Ea = dHrxn + return Ea cpdef Arrhenius to_arrhenius(self, double dHrxn): """ @@ -613,11 +617,12 @@ 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) + self.w0 = (w0 * 0.001, 'kJ/mol') if len(rxns) == 1: rxn = rxns[0] @@ -627,8 +632,8 @@ cdef class ArrheniusBM(KineticsModel): Ea = rxn.kinetics.Ea.value_si 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) + Vp = 2 * w0 * (w0 + E0) / (w0 - E0) + out = Ea - (w0 + dHrxn / 2.0) * (Vp - 2 * w0 + dHrxn) * (Vp - 2 * w0 + dHrxn) / (Vp * Vp - (4 * w0 * w0) + dHrxn * dHrxn) return out E0 = fsolve(kfcn, w0 / 10.0)[0] @@ -642,23 +647,40 @@ cdef class ArrheniusBM(KineticsModel): 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) + Vp = 2 * w0 * (w0 + E0) / (w0 - E0) + Ea = (w0 + dHrxn / 2.0) * (Vp - 2 * w0 + dHrxn) * (Vp - 2 * w0 + dHrxn) / (Vp * Vp - (4 * w0 * 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))) + return lnA + np.log(T) * n + (-Ea / (constants.R * T)) # 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 + dHrxn = rxn.get_enthalpy_of_reaction(298.0) for T in Ts: - xdata.append([T, rxn.get_enthalpy_of_reaction(298.0)]) + xdata.append([T, dHrxn]) ydata.append(np.log(rxn.get_rate_coefficient(T))) - sigmas.append(s / (8.314 * T)) + sigmas.append(s / (constants.R * T)) + # Use BEP with alpha = 0.25 for inital guess of E0 + E0 += rxn.kinetics._Ea.value_si - 0.25 * dHrxn + lnA += np.log(rxn.kinetics.A.value_si) + n += rxn.kinetics.n.value_si + # Use the averages as intial guess + E0 /= len(rxns) + lnA /= len(rxns) + n /= len(rxns) + E0 = min(E0, w0) + w0 = max(2 * E0, w0) # Expression only works if w0>2E0, and is insensitive to w0 + self.w0 = (w0 * 0.001, 'kJ/mol') + if E0 < 0: + E0 = w0 / 100.0 xdata = np.array(xdata) ydata = np.array(ydata) @@ -667,10 +689,19 @@ cdef class ArrheniusBM(KineticsModel): keep_trying = True xtol = 1e-8 ftol = 1e-8 + attempts = 0 while keep_trying: keep_trying = False try: - params = curve_fit(kfcn, xdata, ydata, sigma=sigmas, p0=[1.0, 1.0, w0 / 10.0], xtol=xtol, ftol=ftol) + 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: + keep_trying = True + if attempts > 0: + self.w0.value_si *= 1.25 + w0 = self.w0.value_si # update local variable for use in kfcn optimization function + attempts += 1 + E0 = self.w0.value_si / 10.0 except RuntimeError: if xtol < 1.0: keep_trying = True @@ -679,7 +710,6 @@ cdef class ArrheniusBM(KineticsModel): else: raise ValueError("Could not fit BM arrhenius to reactions with xtol<1.0") - lnA, n, E0 = params[0].tolist() A = np.exp(lnA) self.Tmin = (np.min(Ts), "K") @@ -695,8 +725,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 @@ -1550,16 +1579,20 @@ cdef class ArrheniusChargeTransferBM(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: + # Negative E0 is unphysical, but could be encountered during optimization. + Ea = E0 + max(dHrxn, 0.0) 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) + Vp = 2 * w0 * (w0 + E0) / (w0 - E0) + Ea = (w0 + dHrxn / 2.0) * (Vp - 2 * w0 + dHrxn) * (Vp - 2 * w0 + dHrxn) / (Vp * Vp - (4 * w0 * w0) + dHrxn * dHrxn) + if (Ea < 0.0) or (dHrxn < -4 * E0): + Ea = 0.0 + elif (Ea < dHrxn) or (dHrxn > 4 * E0): + Ea = dHrxn + return Ea cpdef double get_rate_coefficient_from_potential(self, double T, double V, double dHrxn) except -1: """ @@ -1589,11 +1622,12 @@ cdef class ArrheniusChargeTransferBM(KineticsModel): rxn.kinetics.change_v0(0.0) 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] + Ts = np.array([300.0, 500.0, 600.0, 700.0, 800.0, 900.0, 1000.0, 1100.0, 1200.0, 1500.0]) if w0 is None: #estimate w0 w0s = get_w0s(recipe, rxns) w0 = sum(w0s) / len(w0s) + self.w0 = (w0 * 0.001, 'kJ/mol') if len(rxns) == 1: rxn = rxns[0] @@ -1603,8 +1637,8 @@ cdef class ArrheniusChargeTransferBM(KineticsModel): Ea = rxn.kinetics.Ea.value_si 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) + Vp = 2 * w0 * (w0 + E0) / (w0 - E0) + out = Ea - (w0 + dHrxn / 2.0) * (Vp - 2 * w0 + dHrxn) * (Vp - 2 * w0 + dHrxn) / (Vp * Vp - (4 * w0 * w0) + dHrxn * dHrxn) return out E0 = fsolve(kfcn, w0 / 10.0)[0] @@ -1617,23 +1651,40 @@ cdef class ArrheniusChargeTransferBM(KineticsModel): 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) + Vp = 2 * w0 * (w0 + E0) / (w0 - E0) + Ea = (w0 + dHrxn / 2.0) * (Vp - 2 * w0 + dHrxn) * (Vp - 2 * w0 + dHrxn) / (Vp * Vp - (4 * w0 * 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))) + return lnA + np.log(T) * n + (-Ea / (constants.R * T)) # 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 + dHrxn = rxn.get_enthalpy_of_reaction(298.0) for T in Ts: - xdata.append([T, rxn.get_enthalpy_of_reaction(298.0)]) + xdata.append([T, dHrxn]) ydata.append(np.log(rxn.get_rate_coefficient(T))) - sigmas.append(s / (8.314 * T)) + sigmas.append(s / (constants.R * T)) + # Use BEP with alpha = 0.25 for initial guess of E0 + E0 += rxn.kinetics._Ea.value_si - 0.25 * dHrxn + lnA += np.log(rxn.kinetics.A.value_si) + n += rxn.kinetics.n.value_si + # average E0, lnA, n over all reactions + E0 /= len(rxns) + lnA /= len(rxns) + n /= len(rxns) + E0 = min(E0, w0) + w0 = max(2 * E0, w0) # ensure w0 is at least 2*E0 + self.w0 = (w0 * 0.001, 'kJ/mol') + if E0 < 0: + E0 = w0 / 100.0 xdata = np.array(xdata) ydata = np.array(ydata) @@ -1642,10 +1693,19 @@ cdef class ArrheniusChargeTransferBM(KineticsModel): keep_trying = True xtol = 1e-8 ftol = 1e-8 + attempts = 0 while keep_trying: keep_trying = False try: - params = curve_fit(kfcn, xdata, ydata, sigma=sigmas, p0=[1.0, 1.0, w0 / 10.0], xtol=xtol, ftol=ftol) + 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.0 and attempts < 5: + keep_trying = True + if attempts > 0: + self.w0.value_si *= 1.25 + w0 = self.w0.value_si # update local variable for use in kfcn optimization function + attempts += 1 + E0 = self.w0.value_si / 10.0 except RuntimeError: if xtol < 1.0: keep_trying = True @@ -1654,7 +1714,6 @@ cdef class ArrheniusChargeTransferBM(KineticsModel): else: raise ValueError("Could not fit BM arrhenius to reactions with xtol<1.0") - lnA, n, E0 = params[0].tolist() A = np.exp(lnA) self.Tmin = (np.min(Ts), "K") @@ -1669,8 +1728,7 @@ cdef class ArrheniusChargeTransferBM(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') self._V0.value_si = 0.0 self.electrons = rxns[0].electrons