From 8deab1c219f0d4fea82d760461f6b9bba7bc3c14 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Wed, 2 Jul 2025 22:59:46 -0400
Subject: [PATCH 01/19] Minor refactor of kinetics _make_rule and rate
 averaging comments.

---
 rmgpy/data/kinetics/family.py | 13 +++++++++----
 1 file changed, 9 insertions(+), 4 deletions(-)

diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index eb7223836e..6bdd4afb02 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -4596,13 +4596,16 @@ 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 n > 0: # not needed. at this point n should always be > 0
         if isinstance(rs[0].kinetics, Arrhenius):
             arr = ArrheniusBM
         else:
@@ -4772,7 +4775,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 +4785,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 +4795,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(

From d262550d168d73af08ef1abfd328571ffaa782f9 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Fri, 3 Apr 2026 09:59:01 -0400
Subject: [PATCH 02/19] Minor refactor of kinetics _make_rule: outdent an "if
 True:" block

The previous commit refactoring had made this block of code
always run, so in this commit I remove the "if" statement and
reduce the indentation level. Looks like many lines have changed
but it's only whitespace.
---
 rmgpy/data/kinetics/family.py | 125 +++++++++++++++++-----------------
 1 file changed, 62 insertions(+), 63 deletions(-)

diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index 6bdd4afb02..0bf77cf96f 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -4605,73 +4605,72 @@ def _make_rule(rr):
         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: # not needed. at this point n should always be > 0
-        if isinstance(rs[0].kinetics, Arrhenius):
-            arr = ArrheniusBM
+
+    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)
         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)
-            else:
+            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)
+        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(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:
-                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(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
+
 
 
 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):

From 1d5455b89d281dfd30992d2e542216c8c1b4b269 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Wed, 2 Jul 2025 23:01:18 -0400
Subject: [PATCH 03/19] Don't trust Blowers-Masel rate rule with abs(n) > 5.0

These are usually a consequence of trying to fit some reactions
that don't seem to follow Blowers-Masel enthalpy-dependence.

Matt had suggested this in 226af17914369171a12f9a5e7f3dd75ff73b0d9a
but then removed it in d0126887698c3ea04d2391e3389991d1b381af54.
---
 rmgpy/data/kinetics/family.py | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index 0bf77cf96f..fb4ac56313 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -4612,12 +4612,20 @@ def _make_rule(rr):
         arr = ArrheniusChargeTransferBM
     if n > 1:
         kin = arr().fit_to_reactions(rs, recipe=recipe)
-    if n == 1 or kin.E0.value_si < 0.0:
+    if n == 1 or kin.E0.value_si < 0.0 or abs(kin.n.value_si) > 5.0:
+        # 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])
         #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)
         else:
+            if kin.E0.value_si < 0.0:
+                reason = "E0<0"
+            elif abs(kin.n.value_si) > 5.0:
+                reason = "abs(n)>5"
+            else:
+                reason = "?"
+            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 rs[list(set(range(len(rs))) - {i})]]).get_rate_coefficient(T=Tref) / rxn.get_rate_coefficient(T=Tref)

From 3e5253a3853075ce7676b6f14fe90d27bb477228 Mon Sep 17 00:00:00 2001
From: davidfarinajr <farina.d@northeastern.edu>
Date: Tue, 14 Dec 2021 16:58:43 -0500
Subject: [PATCH 04/19] ArrheniusBM.get_activation_energy can cope with
 negative E0

Negative E0 values don't make a lot of sense physically,
but might be encountered when doing an optimization.

Commit message by Richard, in 2025, while rebasing.
Original code by David in 2021.
---
 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 ea9e462d4f..4e663185ec 100644
--- a/rmgpy/kinetics/arrhenius.pyx
+++ b/rmgpy/kinetics/arrhenius.pyx
@@ -573,16 +573,22 @@ 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:
+            if dHrxn > 0:
+                Ea = dHrxn
+            else:
+                Ea = min(0.0, 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 0fa20ee7b5a01f9cb1a1b1c89d5f2c3b6deb6491 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Tue, 1 Jul 2025 21:12:24 -0400
Subject: [PATCH 05/19] Reuse dHrxn since we're in a loop.

It's also useful later.
---
 rmgpy/kinetics/arrhenius.pyx | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx
index 4e663185ec..789f1f167f 100644
--- a/rmgpy/kinetics/arrhenius.pyx
+++ b/rmgpy/kinetics/arrhenius.pyx
@@ -661,8 +661,9 @@ cdef class ArrheniusBM(KineticsModel):
             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))
 

From 7817c92bea4510027050da8afc0ebc8d82607da7 Mon Sep 17 00:00:00 2001
From: davidfarinajr <farina.d@northeastern.edu>
Date: Tue, 14 Dec 2021 16:58:43 -0500
Subject: [PATCH 06/19] Better starting guesses for optimization of lnA, n, E0.

When doing a least squares optimization of the Blowers Masel parameters,
we were using initial guesses of lnA = 1, n = 1, and E0 = w0/10.
We now use values derived from the Arrhenius parameters of the reactions
in question.

Commit message by Richard when rebasing in 2025.
Code changes by David in 2021.
---
 rmgpy/kinetics/arrhenius.pyx | 16 +++++++++++++++-
 1 file changed, 15 insertions(+), 1 deletion(-)

diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx
index 789f1f167f..76a206dd15 100644
--- a/rmgpy/kinetics/arrhenius.pyx
+++ b/rmgpy/kinetics/arrhenius.pyx
@@ -658,6 +658,9 @@ cdef class ArrheniusBM(KineticsModel):
             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
@@ -666,6 +669,17 @@ cdef class ArrheniusBM(KineticsModel):
                     xdata.append([T, dHrxn])
                     ydata.append(np.log(rxn.get_rate_coefficient(T)))
                     sigmas.append(s / (8.314 * 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)
+            if E0 < 0:
+                E0 = w0 / 100.0
 
             xdata = np.array(xdata)
             ydata = np.array(ydata)
@@ -677,7 +691,7 @@ cdef class ArrheniusBM(KineticsModel):
             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)
                 except RuntimeError:
                     if xtol < 1.0:
                         keep_trying = True

From ae48f6ec482192722f93fdbfb7bcf607aa631a24 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Tue, 1 Jul 2025 21:52:04 -0400
Subject: [PATCH 07/19] Simplify and hopefully optimize exponential in fitting.

---
 rmgpy/kinetics/arrhenius.pyx | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx
index 76a206dd15..a350ae8783 100644
--- a/rmgpy/kinetics/arrhenius.pyx
+++ b/rmgpy/kinetics/arrhenius.pyx
@@ -652,7 +652,7 @@ cdef class ArrheniusBM(KineticsModel):
                 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)))
+                return lnA + np.log(T) * n + (-Ea / (8.314472 * T))
 
             # get (T,dHrxn(T)) -> (Ln(k) mappings
             xdata = []

From dd64341abe986e9c642aca6a41b808e170034193 Mon Sep 17 00:00:00 2001
From: davidfarinajr <farina.d@northeastern.edu>
Date: Tue, 14 Dec 2021 16:58:43 -0500
Subject: [PATCH 08/19] Revised Blowers-Masel ArrheniusBM fitting.

If fitted E0 exceeds the w0 then try again,
and gradually increase the w0 up to 5 times.

Use kJ/mol for w0 and E0.

Commit message by Richard when rebasing in 2025.
Original changes by David Farina in 2021.
---
 rmgpy/kinetics/arrhenius.pyx | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx
index a350ae8783..8bae422e94 100644
--- a/rmgpy/kinetics/arrhenius.pyx
+++ b/rmgpy/kinetics/arrhenius.pyx
@@ -619,11 +619,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]
@@ -688,10 +689,18 @@ 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=[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
+                        attempts += 1
+                        E0 = self.w0.value_si / 10.0
                 except RuntimeError:
                     if xtol < 1.0:
                         keep_trying = True
@@ -700,7 +709,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")
@@ -716,8 +724,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 d49b5697a5317eae66ad1bf004c16e6d465efd6b Mon Sep 17 00:00:00 2001
From: Nora Khalil <khalil.nor@FACF2XQ43G76D.tld>
Date: Tue, 28 May 2024 15:25:09 -0400
Subject: [PATCH 09/19] modifying BM fitting so w0 is always at least 2 E0

---
 rmgpy/kinetics/arrhenius.pyx | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx
index 8bae422e94..ccd69ea58b 100644
--- a/rmgpy/kinetics/arrhenius.pyx
+++ b/rmgpy/kinetics/arrhenius.pyx
@@ -679,6 +679,8 @@ cdef class ArrheniusBM(KineticsModel):
             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
 

From 2fc779ac021d10b26d5be286b565b13fd96f2cca Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Wed, 2 Jul 2025 23:21:22 -0400
Subject: [PATCH 10/19] ArrheniuschargeTransferBM.get_activation_energy can
 cope with negative E0

Like the commit
ArrheniusBM.get_activation_energy can cope with negative E0
but for ArrheniuschargeTransferBM because we have horrible
code duplication issues.
---
 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 ccd69ea58b..74a6760555 100644
--- a/rmgpy/kinetics/arrhenius.pyx
+++ b/rmgpy/kinetics/arrhenius.pyx
@@ -1580,16 +1580,22 @@ 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:
+            if dHrxn > 0:
+                Ea = dHrxn
+            else:
+                Ea = min(0.0, 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) * (Vp - 2 * w0 + dHrxn) / (Vp * Vp - (2 * w0) * (2 * w0) + dHrxn * dHrxn)
+            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 double get_rate_coefficient_from_potential(self, double T, double V, double dHrxn) except -1:
         """

From b6ddb4d35fe424f11b66b6d9955c3896a23fd3c4 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Wed, 2 Jul 2025 23:34:59 -0400
Subject: [PATCH 11/19] Update ArrheniusChargeTransferBM.fit_to_reactions

These changes mirror those made to the ArrheniusBM.fit_to_reactions
so look there for detailed commit messages explaining the changes.
This, unfortunately, is a bunch of duplicate code.
---
 rmgpy/kinetics/arrhenius.pyx | 38 +++++++++++++++++++++++++++++-------
 1 file changed, 31 insertions(+), 7 deletions(-)

diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx
index 74a6760555..b9c694f884 100644
--- a/rmgpy/kinetics/arrhenius.pyx
+++ b/rmgpy/kinetics/arrhenius.pyx
@@ -1625,11 +1625,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]
@@ -1657,19 +1658,36 @@ cdef class ArrheniusChargeTransferBM(KineticsModel):
                 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)))
+                return lnA + np.log(T) * n + (-Ea / (8.314 * 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))
+                # 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)
@@ -1678,10 +1696,18 @@ 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
+                        attempts += 1
+                        E0 = self.w0.value_si / 10.0
                 except RuntimeError:
                     if xtol < 1.0:
                         keep_trying = True
@@ -1690,7 +1716,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")
@@ -1705,8 +1730,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
 

From 98ac807437671f753e9dc22a72cc17c71b612e04 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Mon, 6 Apr 2026 10:40:04 -0400
Subject: [PATCH 12/19] Handle negative E0 in ArrheniusBM and
 ArrheniusChargeTransferBM by adjusting activation energy calculation.

A negative intrinsic barrier height is not very physical
given Blowers and Masels original derivation
but it could be encountered in some cases of optimization etc.
This choice of what to do when E0 < 0 seems as good as any.
It is continuous and monotonic when gradually decreasing E0 past E0 = 0. And it allows E0 < 0 values to actually change the resulting Ea values, and for the latter to go negative.
(but not more negative that E0).

Beware code duplication here - need to do things twice.
---
 rmgpy/kinetics/arrhenius.pyx | 12 ++++--------
 1 file changed, 4 insertions(+), 8 deletions(-)

diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx
index b9c694f884..bd86c33632 100644
--- a/rmgpy/kinetics/arrhenius.pyx
+++ b/rmgpy/kinetics/arrhenius.pyx
@@ -577,10 +577,8 @@ cdef class ArrheniusBM(KineticsModel):
         E0 = self._E0.value_si
         w0 = self._w0.value_si
         if E0 < 0:
-            if dHrxn > 0:
-                Ea = dHrxn
-            else:
-                Ea = min(0.0, E0)
+            # Negative E0 is unphysical, but could be encountered during optimization.
+            Ea = E0 + max(dHrxn, 0.0)
         else:
             Vp = 2 * w0 * (2 * w0 + 2 * E0) / (2 * w0 - 2 * E0)
             Ea = (w0 + dHrxn / 2.0) * (Vp - 2 * w0 + dHrxn) ** 2 / (Vp ** 2 - (2 * w0) ** 2 + dHrxn ** 2)
@@ -1584,10 +1582,8 @@ cdef class ArrheniusChargeTransferBM(KineticsModel):
         E0 = self._E0.value_si
         w0 = self._w0.value_si
         if E0 < 0:
-            if dHrxn > 0:
-                Ea = dHrxn
-            else:
-                Ea = min(0.0, E0)
+            # Negative E0 is unphysical, but could be encountered during optimization.
+            Ea = E0 + max(dHrxn, 0.0)
         else:
             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)

From 3b1461a96994f2244d3c8125bd53cc0d13bef67e Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Mon, 6 Apr 2026 11:00:33 -0400
Subject: [PATCH 13/19] Refactor ArrheniusBM and ArrheniusChargeTransferBM

To improve activation energy calculations for negative E0 scenarios.
These changes shouldn't actually change the result.
Just some slight code optimization (eg. x*x instead of x**2)
and simplification (eg. do (x/y) instead of (2*x)/(2*y))

Beware of code duplication.
---
 rmgpy/kinetics/arrhenius.pyx | 32 ++++++++++++++++----------------
 1 file changed, 16 insertions(+), 16 deletions(-)

diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx
index bd86c33632..0f64aa0c11 100644
--- a/rmgpy/kinetics/arrhenius.pyx
+++ b/rmgpy/kinetics/arrhenius.pyx
@@ -580,11 +580,11 @@ cdef class ArrheniusBM(KineticsModel):
             # Negative E0 is unphysical, but could be encountered during optimization.
             Ea = E0 + max(dHrxn, 0.0)
         else:
-            Vp = 2 * w0 * (2 * w0 + 2 * E0) / (2 * w0 - 2 * E0)
-            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):
+            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 (dHrxn > 0.0 and Ea < dHrxn) or (dHrxn > 4 * E0):
+            elif (Ea < dHrxn) or (dHrxn > 4 * E0):
                 Ea = dHrxn
         return Ea
 
@@ -632,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]
@@ -647,8 +647,8 @@ 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 + (-Ea / (8.314472 * T))
@@ -1585,11 +1585,11 @@ cdef class ArrheniusChargeTransferBM(KineticsModel):
             # Negative E0 is unphysical, but could be encountered during optimization.
             Ea = E0 + max(dHrxn, 0.0)
         else:
-            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)
-            if (dHrxn < 0.0 and Ea < 0.0) or (dHrxn < -4 * E0):
+            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 (dHrxn > 0.0 and Ea < dHrxn) or (dHrxn > 4 * E0):
+            elif (Ea < dHrxn) or (dHrxn > 4 * E0):
                 Ea = dHrxn
         return Ea
 
@@ -1636,8 +1636,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]
@@ -1650,8 +1650,8 @@ 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 + (-Ea / (8.314 * T))

From 3422d0ba27a02f36d50ca7d11ccefa9a2bd37abe Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Mon, 6 Apr 2026 14:35:02 -0400
Subject: [PATCH 14/19] Use constants.R instead of 8.314 or 8.314472

---
 rmgpy/data/kinetics/family.py | 4 ++--
 rmgpy/kinetics/arrhenius.pyx  | 8 ++++----
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index fb4ac56313..b2eab92084 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -4631,7 +4631,7 @@ def _make_rule(rr):
                         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
+            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()
@@ -4659,7 +4659,7 @@ def _make_rule(rr):
                         .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
+            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()
diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx
index 0f64aa0c11..d46a82c5e1 100644
--- a/rmgpy/kinetics/arrhenius.pyx
+++ b/rmgpy/kinetics/arrhenius.pyx
@@ -651,7 +651,7 @@ cdef class ArrheniusBM(KineticsModel):
                 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 + (-Ea / (8.314472 * T))
+                return lnA + np.log(T) * n + (-Ea / (constants.R * T))
 
             # get (T,dHrxn(T)) -> (Ln(k) mappings
             xdata = []
@@ -667,7 +667,7 @@ cdef class ArrheniusBM(KineticsModel):
                 for T in Ts:
                     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)
@@ -1654,7 +1654,7 @@ cdef class ArrheniusChargeTransferBM(KineticsModel):
                 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 + (-Ea / (8.314 * T))
+                return lnA + np.log(T) * n + (-Ea / (constants.R * T))
 
             # get (T,dHrxn(T)) -> (Ln(k) mappings
             xdata = []
@@ -1670,7 +1670,7 @@ cdef class ArrheniusChargeTransferBM(KineticsModel):
                 for T in Ts:
                     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)

From 42718a82544370e5081ac262a8c5fc98dca9a591 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Mon, 6 Apr 2026 14:45:02 -0400
Subject: [PATCH 15/19] Fix w0 bug in Blowers-Masel fitting.

We changed the w0, but not the instance of it that's used inside the
objective function inside the fitting loop.

(Beware code duplication in ArrheniusBM and ArrheniusChargeTransferBM)
---
 rmgpy/kinetics/arrhenius.pyx | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx
index d46a82c5e1..6252b6ce95 100644
--- a/rmgpy/kinetics/arrhenius.pyx
+++ b/rmgpy/kinetics/arrhenius.pyx
@@ -699,6 +699,7 @@ cdef class ArrheniusBM(KineticsModel):
                         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:
@@ -1702,6 +1703,7 @@ cdef class ArrheniusChargeTransferBM(KineticsModel):
                         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:

From 159069914b2a9ff535278cb912d3b9a1836eee05 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Mon, 6 Apr 2026 14:46:38 -0400
Subject: [PATCH 16/19] Simplify some array slicing for leave-one-out
 selection.

For the leave-one-out selection of reactions from rs,
replace
    rs[list(set(range(len(rs))) - {i})]
with
    np.delete(rs, i)
(which returns a copy)
---
 rmgpy/data/kinetics/family.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index b2eab92084..e319017a88 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -4628,7 +4628,7 @@ def _make_rule(rr):
             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 rs[list(set(range(len(rs))) - {i})]]).get_rate_coefficient(T=Tref) / rxn.get_rate_coefficient(T=Tref)
+                        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
@@ -4646,7 +4646,7 @@ def _make_rule(rr):
             if isinstance(rs[0].kinetics, Arrhenius):
                 dlnks = np.array([
                     np.log(
-                        arr().fit_to_reactions(rs[list(set(range(len(rs))) - {i})], recipe=recipe)
+                        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)
@@ -4654,7 +4654,7 @@ def _make_rule(rr):
             else:
                 dlnks = np.array([
                     np.log(
-                        arr().fit_to_reactions(rs[list(set(range(len(rs))) - {i})], recipe=recipe)
+                        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)

From 3d0b9b1c0004e6e95ac5714a4a4acb6acb4e3d3c Mon Sep 17 00:00:00 2001
From: Nora Khalil <khalil.nor@northeastern.edu>
Date: Wed, 8 Apr 2026 11:28:56 -0400
Subject: [PATCH 17/19] Refactoring to avoid AttributeError in _make_rule

---
 rmgpy/data/kinetics/family.py | 17 +++++++++++------
 1 file changed, 11 insertions(+), 6 deletions(-)

diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index e319017a88..764c6a3c7c 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -4613,18 +4613,23 @@ def _make_rule(rr):
     if n > 1:
         kin = arr().fit_to_reactions(rs, recipe=recipe)
     if n == 1 or kin.E0.value_si < 0.0 or abs(kin.n.value_si) > 5.0:
-        # 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])
-        #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)
-        else:
+        
+        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(

From 5cbdf83098758c4d2cb926fd7cd1f0d8a6cb3375 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Thu, 9 Apr 2026 22:45:21 -0400
Subject: [PATCH 18/19] Update the cross_validate method to match the
 _make_rule method, for fitting kinetics.

I tried using a helper function that could be re-used, but it started to get too messy
and complicated. This unfortunately requires continued maintenance keeping three
different bits of code in sync doing the same thing, eg. rejecting BM fits when abs(n)>5
or E0<0.
---
 rmgpy/data/kinetics/family.py | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index 764c6a3c7c..441e57fce1 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.))
@@ -4612,8 +4615,9 @@ def _make_rule(rr):
         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"
@@ -4621,15 +4625,13 @@ def _make_rule(rr):
                 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:
 
+        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(

From d9ab313c024b9f0a9cde3c2f7d8dbbcb55f1441b Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Thu, 9 Apr 2026 23:08:06 -0400
Subject: [PATCH 19/19] Modernize display.py IPython detection

Replace deprecated InteractiveShell.initialized() with get_ipython(),
the canonical way to detect an active IPython/Jupyter session.
Also update the docstring to drop the legacy --pylab reference.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 rmgpy/display.py | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/rmgpy/display.py b/rmgpy/display.py
index 3b6c751a10..cf9a40d2d2 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