From 6e323901c901c37f163d535efc0f20e538530745 Mon Sep 17 00:00:00 2001
From: Calvin Pieters <calvinpieters@gmail.com>
Date: Thu, 30 Jul 2026 13:58:08 +0300
Subject: [PATCH 1/8] Reject zero rate coefficients when fitting Arrhenius in
 log space

Arrhenius.fit_to_data and ArrheniusChargeTransfer.fit_to_data validated
their input rate coefficients against inf and NaN, and against mixed
signs, but not against zero. Zero passes both checks and reaches the
least-squares step, which is performed in log space: log(0) is -inf, so
numpy emits "divide by zero encountered in log" and lstsq then returns
NaN for every fitted parameter without raising.

The NaN expression is returned to the caller as if the fit had
succeeded. It surfaces much later and far from the cause -- evaluating
it raises "TypeError: Cannot convert 'complex' with non-zero imaginary
component to 'double'" from the '**' operator, which reads like a
negative temperature rather than a degenerate fit. It is also possible
for such an expression to reach the Chemkin output instead, in which
case the run completes and writes NaN kinetics with no error at all.

A rate coefficient of exactly zero typically comes from underflow at the
low end of a fitted temperature range, which is common when generating a
reverse rate coefficient for a strongly endothermic reaction. Reject it
at the same point inf and NaN are rejected, so the failure is reported
where it originates.
---
 rmgpy/kinetics/arrhenius.pyx         |  4 ++++
 test/rmgpy/kinetics/arrheniusTest.py | 16 ++++++++++++++++
 2 files changed, 20 insertions(+)

diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx
index 43b58dd8304..150cf252086 100644
--- a/rmgpy/kinetics/arrhenius.pyx
+++ b/rmgpy/kinetics/arrhenius.pyx
@@ -158,6 +158,8 @@ cdef class Arrhenius(KineticsModel):
         import scipy.stats
         if not all(np.isfinite(klist)):
             raise  ValueError("Rates must all be finite, not inf or NaN")
+        if any(klist==0):
+            raise ValueError("Rates must all be nonzero; a zero rate coefficient cannot be fit in log space.")
         if any(klist<0):
             if not all(klist<0):
                 raise ValueError("Rates must all be positive or all be negative.")
@@ -1377,6 +1379,8 @@ cdef class ArrheniusChargeTransfer(KineticsModel):
         import scipy.stats
         if not all(np.isfinite(klist)):
             raise  ValueError("Rates must all be finite, not inf or NaN")
+        if any(klist==0):
+            raise ValueError("Rates must all be nonzero; a zero rate coefficient cannot be fit in log space.")
         if any(klist<0):
             if not all(klist<0):
                 raise ValueError("Rates must all be positive or all be negative.")
diff --git a/test/rmgpy/kinetics/arrheniusTest.py b/test/rmgpy/kinetics/arrheniusTest.py
index 5913c851e91..947811c71f8 100644
--- a/test/rmgpy/kinetics/arrheniusTest.py
+++ b/test/rmgpy/kinetics/arrheniusTest.py
@@ -177,6 +177,22 @@ def test_fit_to_data(self):
         assert round(abs(arrhenius.Ea.value_si - self.arrhenius.Ea.value_si), 2) == 0
         assert round(abs(arrhenius.T0.value_si - self.arrhenius.T0.value_si), 4) == 0
 
+    def test_fit_to_data_with_zero_rate(self):
+        """
+        Test that Arrhenius.fit_to_data() rejects a rate coefficient of exactly zero.
+
+        Zero passes the finite check and the sign check, but the fit is performed in log space,
+        so log(0) = -inf enters the least-squares problem and every fitted parameter comes back
+        NaN without an exception being raised. Evaluating that expression then fails far from
+        the cause, with a TypeError about converting a complex number.
+        """
+        Tdata = np.array([300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500])
+        kdata = np.array([self.arrhenius.get_rate_coefficient(T) for T in Tdata])
+        kdata[0] = 0.0  # as if the low-temperature rate had underflowed
+
+        with pytest.raises(ValueError, match="nonzero"):
+            Arrhenius().fit_to_data(Tdata, kdata, kunits="m^3/(mol*s)")
+
     def test_fit_to_negative_data(self):
         """
         Test the Arrhenius.fit_to_data() method on negative rates

From 2a1015c5b98ad3bee388a6a61ab27554e8414b95 Mon Sep 17 00:00:00 2001
From: Calvin Pieters <calvinpieters@gmail.com>
Date: Thu, 30 Jul 2026 13:58:08 +0300
Subject: [PATCH 2/8] Keep collision limit checks alive when a rate cannot be
 evaluated

check_collision_limit_violation evaluated a rate coefficient for every
extreme T/P condition and compared it against the collision limit for
that condition. Any reaction whose rate could not be evaluated -- an
unsupported kinetics type raising ReactionError from
generate_reverse_rate_coefficient, or a degenerate reverse fit -- raised
out of check_model and aborted the run at the very end, after the model
had otherwise completed.

Evaluate each direction independently and record the rate, its collision
limit and the condition together, so that a condition which cannot be
evaluated is skipped with a warning instead of propagating. Storing the
three values as one tuple rather than in parallel lists also removes the
possibility of the comparison loop pairing a rate with another
condition's limit: the collision limit was previously appended before
the rate was evaluated, so a rate that failed left an orphan limit
behind and shifted every later comparison. Since the collision limit
grows with sqrt(T), that skew compared a high-temperature rate against a
smaller low-temperature limit and reported violations that do not exist.

In check_model, keep the broad catch as a backstop -- nothing detected
by this end-of-run diagnostic justifies discarding a completed run --
but log the traceback so a genuine bug behind the failure stays
diagnosable rather than being reduced to a one-line message.
---
 rmgpy/reaction.py          |  53 ++++++++------
 rmgpy/rmg/main.py          |  11 ++-
 test/rmgpy/reactionTest.py | 144 +++++++++++++++++++++++++++++++++++++
 test/rmgpy/rmg/mainTest.py |  59 +++++++++++++++
 4 files changed, 246 insertions(+), 21 deletions(-)

diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py
index 13817852863..d3a071c1eb5 100644
--- a/rmgpy/reaction.py
+++ b/rmgpy/reaction.py
@@ -1748,37 +1748,50 @@ def check_collision_limit_violation(self, t_min, t_max, p_min, p_max):
                 conditions.append([t_max, p_max])
         logging.debug("Checking whether reaction {0} violates the collision rate limit...".format(self))
         violator_list = []
-        kf_list = []
-        kr_list = []
-        collision_limit_f = []
-        collision_limit_r = []
+        forward_checks = []
+        reverse_checks = []
         for condition in conditions:
+            temp, pressure = condition
             if len(self.reactants) >= 2:
                 try:
-                    collision_limit_f.append(self.calculate_coll_limit(temp=condition[0], reverse=False))
+                    limit_f = self.calculate_coll_limit(temp=temp, reverse=False)
                 except ValueError:
                     continue
                 else:
-                    kf_list.append(self.get_rate_coefficient(condition[0], condition[1]))
+                    try:
+                        kf = self.get_rate_coefficient(temp, pressure)
+                    except (ReactionError, KineticsError, TypeError, ValueError,
+                            ZeroDivisionError, OverflowError) as err:
+                        logging.warning(
+                            "Skipping forward collision limit check for reaction %s at %.1f K, %.3g Pa "
+                            "because rate evaluation failed: %s",
+                            self, temp, pressure, err,
+                        )
+                    else:
+                        forward_checks.append((kf, limit_f, condition))
             if len(self.products) >= 2:
                 try:
-                    collision_limit_r.append(self.calculate_coll_limit(temp=condition[0], reverse=True))
+                    limit_r = self.calculate_coll_limit(temp=temp, reverse=True)
                 except ValueError:
                     continue
                 else:
-                    kr_list.append(self.generate_reverse_rate_coefficient().get_rate_coefficient(condition[0], condition[1]))
-        if len(self.reactants) >= 2:
-            for i, k in enumerate(kf_list):
-                if k > collision_limit_f[i]:
-                    ratio = k / collision_limit_f[i]
-                    condition = '{0} K, {1:.1f} bar'.format(conditions[i][0], conditions[i][1] / 1e5)
-                    violator_list.append([self, 'forward', ratio, condition])
-        if len(self.products) >= 2:
-            for i, k in enumerate(kr_list):
-                if k > collision_limit_r[i]:
-                    ratio = k / collision_limit_r[i]
-                    condition = '{0} K, {1:.1f} bar'.format(conditions[i][0], conditions[i][1] / 1e5)
-                    violator_list.append([self, 'reverse', ratio, condition])
+                    try:
+                        kr = self.generate_reverse_rate_coefficient().get_rate_coefficient(temp, pressure)
+                    except (ReactionError, KineticsError, TypeError, ValueError,
+                            ZeroDivisionError, OverflowError) as err:
+                        logging.warning(
+                            "Skipping reverse collision limit check for reaction %s at %.1f K, %.3g Pa "
+                            "because reverse rate evaluation failed: %s",
+                            self, temp, pressure, err,
+                        )
+                    else:
+                        reverse_checks.append((kr, limit_r, condition))
+        for direction, checks in (('forward', forward_checks), ('reverse', reverse_checks)):
+            for k, collision_limit, (temp, pressure) in checks:
+                if k > collision_limit:
+                    ratio = k / collision_limit
+                    condition = '{0} K, {1:.1f} bar'.format(temp, pressure / 1e5)
+                    violator_list.append([self, direction, ratio, condition])
         return violator_list
 
     def calculate_coll_limit(self, temp, reverse=False):
diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py
index eaa4a732eba..21815c5e535 100644
--- a/rmgpy/rmg/main.py
+++ b/rmgpy/rmg/main.py
@@ -1647,7 +1647,16 @@ def check_model(self):
             if rxn.is_surface_reaction():
                 # Don't check collision limits for surface reactions.
                 continue
-            violator_list = rxn.check_collision_limit_violation(t_min=self.Tmin, t_max=self.Tmax, p_min=self.Pmin, p_max=self.Pmax)
+            try:
+                violator_list = rxn.check_collision_limit_violation(
+                    t_min=self.Tmin, t_max=self.Tmax, p_min=self.Pmin, p_max=self.Pmax
+                )
+            except Exception:
+                logging.warning(
+                    "Skipping collision limit check for reaction %s because evaluation failed.",
+                    rxn, exc_info=True,
+                )
+                continue
             if violator_list:
                 violators.extend(violator_list)
         # Whether or not violators were found, rename 'collision_rate_violators.log' if it exists
diff --git a/test/rmgpy/reactionTest.py b/test/rmgpy/reactionTest.py
index 458de09f6a1..93f1703e938 100644
--- a/test/rmgpy/reactionTest.py
+++ b/test/rmgpy/reactionTest.py
@@ -31,6 +31,7 @@
 This module contains unit tests of the rmgpy.reaction module.
 """
 
+import logging
 import math
 
 import cantera as ct
@@ -66,6 +67,7 @@
 from rmgpy.statmech.translation import IdealGasTranslation
 from rmgpy.statmech.vibration import HarmonicOscillator
 from rmgpy.thermo import Wilhoit, ThermoData, NASA, NASAPolynomial
+from rmgpy.transport import TransportData
 
 def order_of_magnitude(number):
     return math.floor(math.log(number, 10))
@@ -3247,3 +3249,145 @@ def test_reverse_surface_charge_transfer_rate(self):
                 kr = kr_oxidation.get_rate_coefficient(T,V)
                 K = self.rxn_oxidation.get_equilibrium_constant(T,V)
                 assert order_of_magnitude(kf/kr) == order_of_magnitude(K)
+
+
+class TestCollisionLimitViolation:
+    """
+    Contains unit tests of the Reaction.check_collision_limit_violation() method.
+    """
+
+    def setup_class(self):
+        """
+        A bimolecular reaction on both sides, so that the forward and reverse directions are
+        both checked. All four species carry transport data (required by calculate_coll_limit)
+        and thermo (required by generate_reverse_rate_coefficient).
+        """
+        self.ch3 = Species(
+            label="CH3",
+            molecule=[Molecule().from_smiles("[CH3]")],
+            transport_data=TransportData(sigma=(3.8, "angstrom"), epsilon=(144, "K")),
+            thermo=ThermoData(
+                Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"),
+                Cpdata=([9.397, 10.123, 10.856, 11.571, 12.899, 14.055, 16.195], "cal/(mol*K)"),
+                H298=(9.357, "kcal/mol"),
+                S298=(45.174, "cal/(mol*K)"),
+                Cp0=(4.0 * constants.R, "J/(mol*K)"),
+                CpInf=(10.0 * constants.R, "J/(mol*K)"),
+            ),
+        )
+        self.ch4 = Species(
+            label="CH4",
+            molecule=[Molecule().from_smiles("C")],
+            transport_data=TransportData(sigma=(3.746, "angstrom"), epsilon=(141.4, "K")),
+            thermo=ThermoData(
+                Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"),
+                Cpdata=([8.615, 9.687, 10.963, 12.301, 14.841, 16.976, 20.528], "cal/(mol*K)"),
+                H298=(-17.714, "kcal/mol"),
+                S298=(44.472, "cal/(mol*K)"),
+                Cp0=(4.0 * constants.R, "J/(mol*K)"),
+                CpInf=(13.0 * constants.R, "J/(mol*K)"),
+            ),
+        )
+        self.c2h5 = Species(
+            label="C2H5",
+            molecule=[Molecule().from_smiles("C[CH2]")],
+            transport_data=TransportData(sigma=(4.302, "angstrom"), epsilon=(252.3, "K")),
+            thermo=ThermoData(
+                Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"),
+                Cpdata=([11.635, 13.744, 16.085, 18.246, 21.885, 24.676, 29.107], "cal/(mol*K)"),
+                H298=(29.496, "kcal/mol"),
+                S298=(56.687, "cal/(mol*K)"),
+                Cp0=(4.0 * constants.R, "J/(mol*K)"),
+                CpInf=(19.0 * constants.R, "J/(mol*K)"),
+            ),
+        )
+        self.c2h6 = Species(
+            label="C2H6",
+            molecule=[Molecule().from_smiles("CC")],
+            transport_data=TransportData(sigma=(4.302, "angstrom"), epsilon=(252.3, "K")),
+            thermo=ThermoData(
+                Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"),
+                Cpdata=([12.684, 15.506, 18.326, 20.971, 25.500, 29.016, 34.595], "cal/(mol*K)"),
+                H298=(-19.521, "kcal/mol"),
+                S298=(54.799, "cal/(mol*K)"),
+                Cp0=(4.0 * constants.R, "J/(mol*K)"),
+                CpInf=(22.0 * constants.R, "J/(mol*K)"),
+            ),
+        )
+
+    def make_reaction(self, A):
+        """Return CH4 + C2H5 <=> CH3 + C2H6 with the given Arrhenius pre-exponential."""
+        return Reaction(
+            reactants=[self.ch4, self.c2h5],
+            products=[self.ch3, self.c2h6],
+            kinetics=Arrhenius(A=(A, "m^3/(mol*s)"), n=0, Ea=(0, "kcal/mol"), T0=(1, "K")),
+        )
+
+    def test_no_violation_for_physical_rate(self):
+        """A sane rate coefficient must not be reported as a collision limit violator."""
+        rxn = self.make_reaction(A=1e3)
+        assert rxn.check_collision_limit_violation(t_min=300, t_max=1500, p_min=1e5, p_max=1e6) == []
+
+    def test_violation_is_detected_and_labelled(self):
+        """An absurdly large rate coefficient must be reported in both directions."""
+        rxn = self.make_reaction(A=1e20)
+        violators = rxn.check_collision_limit_violation(t_min=300, t_max=1500, p_min=1e5, p_max=1e6)
+
+        directions = [v[1] for v in violators]
+        assert "forward" in directions
+        assert "reverse" in directions
+        for violator_rxn, direction, ratio, condition in violators:
+            assert violator_rxn is rxn
+            assert ratio > 1.0
+            # The condition string must name a temperature that was actually evaluated.
+            assert condition in ("300.0 K, 1.0 bar", "1500.0 K, 1.0 bar")
+
+    def test_rate_and_limit_stay_aligned_when_a_condition_fails(self, caplog):
+        """
+        Regression test: if the rate coefficient cannot be evaluated at one condition, the
+        surviving rates must still be compared against *their own* collision limits.
+
+        Previously the collision limit was appended before the rate was evaluated, so a failure
+        at t_min left an orphaned limit behind and every later rate was compared against the
+        wrong limit. Because the collision limit grows with sqrt(T), that skew compared the
+        high-temperature rate against the smaller low-temperature limit and manufactured a
+        spurious violation.
+        """
+        t_min, t_max, pressure = 300.0, 1500.0, 1e5
+        reference = self.make_reaction(A=1e3)
+
+        # A rate that is under the limit at t_max but would exceed the (smaller) limit at t_min.
+        limit_min = reference.calculate_coll_limit(temp=t_min, reverse=False)
+        limit_max = reference.calculate_coll_limit(temp=t_max, reverse=False)
+        assert limit_min < limit_max, "collision limit is expected to increase with temperature"
+        k_at_t_max = 0.5 * (limit_min + limit_max)
+        assert limit_min < k_at_t_max < limit_max
+
+        class RateFailsAtTmin(Reaction):
+            """
+            Reaction whose forward rate cannot be evaluated at t_min. Subclassing is used rather
+            than patching the instance because Reaction is a cdef class with read-only attributes.
+            """
+
+            def get_rate_coefficient(self, T, P=0, surface_site_density=0, potential=0):
+                if T <= t_min:
+                    raise ValueError("simulated rate evaluation failure at t_min")
+                return k_at_t_max
+
+        rxn = RateFailsAtTmin(
+            reactants=[self.ch4, self.c2h5],
+            products=[self.ch3, self.c2h6],
+            kinetics=Arrhenius(A=(1e3, "m^3/(mol*s)"), n=0, Ea=(0, "kcal/mol"), T0=(1, "K")),
+        )
+        with caplog.at_level(logging.WARNING):
+            violators = rxn.check_collision_limit_violation(
+                t_min=t_min, t_max=t_max, p_min=pressure, p_max=pressure
+            )
+
+        # Guard the premise: if the override ever stopped reaching the compiled caller, the
+        # assertion below would pass vacuously because the real rate is also under the limit.
+        assert "Skipping forward collision limit check" in caplog.text
+
+        # k_at_t_max < limit_max, so the forward direction must not be reported. If the surviving
+        # rate were compared against the orphaned t_min limit instead, it would be.
+        assert [v for v in violators if v[1] == "forward"] == []
diff --git a/test/rmgpy/rmg/mainTest.py b/test/rmgpy/rmg/mainTest.py
index 00aa547e9f8..697d000b056 100644
--- a/test/rmgpy/rmg/mainTest.py
+++ b/test/rmgpy/rmg/mainTest.py
@@ -37,8 +37,12 @@
 
 from rmgpy import get_path, settings
 from rmgpy.data.rmg import RMGDatabase
+from rmgpy.kinetics import Arrhenius
+from rmgpy.molecule import Molecule
+from rmgpy.reaction import Reaction
 from rmgpy.rmg.main import RMG, RMG_Memory, initialize_log, make_profile_graph
 from rmgpy.rmg.model import CoreEdgeReactionModel
+from rmgpy.species import Species
 
 originalPath = get_path()
 
@@ -593,3 +597,58 @@ def test_chemkin_to_cantera_conversion(self):
             # clean up
             os.chdir(originalPath)
             shutil.rmtree(self.dir_name)
+
+
+class TestCheckModelCollisionLimits:
+    """
+    Unit tests for the collision limit portion of RMG.check_model().
+    """
+
+    @staticmethod
+    def make_reaction(cls):
+        """Build a reaction of the given Reaction subclass with real species and kinetics."""
+        a = Species(label="A", molecule=[Molecule().from_smiles("C")])
+        b = Species(label="B", molecule=[Molecule().from_smiles("[CH3]")])
+        return cls(
+            reactants=[a, a],
+            products=[b, b],
+            kinetics=Arrhenius(A=(1e10, "m^3/(mol*s)"), n=0, Ea=(0, "kcal/mol"), T0=(1, "K")),
+        )
+
+    def make_rmg(self, tmp_path, core_reactions):
+        rmg = RMG()
+        rmg.reaction_model = CoreEdgeReactionModel()
+        rmg.reaction_model.core.species = []
+        rmg.reaction_model.edge.species = []
+        rmg.reaction_model.core.reactions = core_reactions
+        rmg.Tmin, rmg.Tmax = 300.0, 1500.0
+        rmg.Pmin, rmg.Pmax = 1.0e5, 1.0e6
+        rmg.output_directory = str(tmp_path)
+        return rmg
+
+    def test_failing_reaction_does_not_abort_the_check(self, tmp_path, caplog):
+        """
+        A reaction whose collision limit cannot be evaluated must be skipped with a warning,
+        and violators from the other core reactions must still be reported.
+        """
+
+        class ExplodingReaction(Reaction):
+            def check_collision_limit_violation(self, t_min, t_max, p_min, p_max):
+                raise ValueError("simulated collision limit failure")
+
+        class ViolatingReaction(Reaction):
+            def check_collision_limit_violation(self, t_min, t_max, p_min, p_max):
+                return [[self, "forward", 12.5, "1500.0 K, 1.0 bar"]]
+
+        exploding = self.make_reaction(ExplodingReaction)
+        violating = self.make_reaction(ViolatingReaction)
+        rmg = self.make_rmg(tmp_path, [exploding, violating])
+
+        with caplog.at_level(logging.WARNING):
+            rmg.check_model()
+
+        assert "Skipping collision limit check" in caplog.text
+        # The surviving reaction's violation must still make it into the report.
+        report = tmp_path / "collision_rate_violators.log"
+        assert report.is_file()
+        assert "Violation factor: 12.50" in report.read_text()

From 44651579222a72d1a46151fa81c80e8fded898c4 Mon Sep 17 00:00:00 2001
From: Calvin Pieters <calvinpieters@gmail.com>
Date: Thu, 30 Jul 2026 13:58:08 +0300
Subject: [PATCH 3/8] Derive the reverse rate derivative from Keq instead of
 kb/kf

compute_rate_derivative formed the reverse contribution to d(dy/dt)/dk
as (kb/kf) * [products]. Every reactor builds kb as kf/Keq, so when a
large activation energy drives kf to exactly zero, kb is zero too and
that ratio is the indeterminate 0.0/0.0. In C this does not raise: it
evaluates to NaN, which then spreads through flux, whole columns of
rate_deriv and the sensitivity right-hand side. The function only runs
under sensitivity analysis, so the integration itself was unaffected,
but any job with sensitivity enabled could produce a NaN sensitivity
matrix with no indication of where it came from.

Since kb is always kf/Keq, the ratio kb/kf is identically 1/Keq, so
divide by the stored Keq directly. This gives the same value whenever
kf > 0 and remains correct when kf underflows, where the true limit is
[products]/Keq rather than zero -- for a strongly reverse-dominated
reaction that value is large, so treating it as zero misreported the
reaction as forward-only.

Guard with `not (Keq[j] > 0)` rather than a test against zero. Keq is
inf for an irreversible reaction under SimpleReactor and zero under the
liquid and surface reactors, so both sentinels must be handled, and
writing the test this way also rejects a NaN Keq arising from bad
thermochemistry.
---
 rmgpy/solver/base.pyx           | 18 ++++----
 test/rmgpy/solver/simpleTest.py | 74 +++++++++++++++++++++++++++++++++
 2 files changed, 84 insertions(+), 8 deletions(-)

diff --git a/rmgpy/solver/base.pyx b/rmgpy/solver/base.pyx
index ad088aed9af..dca598bb3df 100644
--- a/rmgpy/solver/base.pyx
+++ b/rmgpy/solver/base.pyx
@@ -1323,7 +1323,7 @@ cdef class ReactionSystem(DASx):
         k_j is the rate parameter for the jth core reaction.
         """
         cdef np.ndarray[np.int_t, ndim=2] ir, ip
-        cdef np.ndarray[np.float64_t, ndim=1] kf, kr, C, deriv
+        cdef np.ndarray[np.float64_t, ndim=1] kf, Keq, C, deriv
         cdef np.ndarray[np.float64_t, ndim=2] rate_deriv
         cdef double fderiv, rderiv, flux, V
         cdef int j, num_core_reactions, num_core_species
@@ -1334,7 +1334,7 @@ cdef class ReactionSystem(DASx):
         ip = self.product_indices
 
         kf = self.kf
-        kr = self.kb
+        Keq = self.Keq
 
         num_core_reactions = len(self.core_reaction_rates)
         num_core_species = len(self.core_species_concentrations)
@@ -1355,12 +1355,14 @@ cdef class ReactionSystem(DASx):
             else:  # three reactants
                 fderiv = C[ir[j, 0]] * C[ir[j, 1]] * C[ir[j, 2]]
 
-            if ip[j, 1] == -1:  # only one reactant
-                rderiv = kr[j] / kf[j] * C[ip[j, 0]]
-            elif ip[j, 2] == -1:  # only two reactants
-                rderiv = kr[j] / kf[j] * C[ip[j, 0]] * C[ip[j, 1]]
-            else:  # three reactants
-                rderiv = kr[j] / kf[j] * C[ip[j, 0]] * C[ip[j, 1]] * C[ip[j, 2]]
+            if not (Keq[j] > 0.0):
+                rderiv = 0.0
+            elif ip[j, 1] == -1:  # only one product
+                rderiv = C[ip[j, 0]] / Keq[j]
+            elif ip[j, 2] == -1:  # only two products
+                rderiv = C[ip[j, 0]] * C[ip[j, 1]] / Keq[j]
+            else:  # three products
+                rderiv = C[ip[j, 0]] * C[ip[j, 1]] * C[ip[j, 2]] / Keq[j]
 
             flux = fderiv - rderiv
             gderiv = rderiv * kf[j] * RT_inverse
diff --git a/test/rmgpy/solver/simpleTest.py b/test/rmgpy/solver/simpleTest.py
index a6646402016..a63b3ff322c 100644
--- a/test/rmgpy/solver/simpleTest.py
+++ b/test/rmgpy/solver/simpleTest.py
@@ -785,3 +785,77 @@ def test_get_const_spc_indices(self):
 
         # Only "CH4" should be marked constant
         assert rxn_system.const_spc_indices == [core_species.index(a)]
+
+    def test_compute_rate_derivative_kf_underflow(self):
+        """
+        Test that compute_rate_derivative() stays finite when the forward rate coefficient
+        underflows to exactly zero.
+
+        The reverse rate coefficient is built as kb = kf / Keq, so the reverse contribution to
+        the rate derivative was computed as (kb / kf) * [products]. When a large activation
+        energy drives kf to exactly 0.0, that ratio is the indeterminate 0.0/0.0, which
+        evaluates to NaN in C without raising, and the NaN then spreads through the whole
+        sensitivity matrix. Dividing by Keq directly avoids the indeterminate form.
+        """
+        ch3 = Species(
+            molecule=[Molecule().from_smiles("[CH3]")],
+            thermo=ThermoData(
+                Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"),
+                Cpdata=(
+                    [9.397, 10.123, 10.856, 11.571, 12.899, 14.055, 16.195],
+                    "cal/(mol*K)",
+                ),
+                H298=(9.357, "kcal/mol"),
+                S298=(45.174, "cal/(mol*K)"),
+            ),
+        )
+        c2h6 = Species(
+            molecule=[Molecule().from_smiles("CC")],
+            thermo=ThermoData(
+                Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"),
+                Cpdata=(
+                    [12.684, 15.506, 18.326, 20.971, 25.500, 29.016, 34.595],
+                    "cal/(mol*K)",
+                ),
+                H298=(-19.521, "kcal/mol"),
+                S298=(54.799, "cal/(mol*K)"),
+            ),
+        )
+
+        # An activation energy this large gives exp(-Ea/RT) ~ exp(-1600) at 300 K, which is far
+        # below the smallest representable double and therefore underflows to exactly 0.0.
+        rxn = Reaction(
+            reactants=[c2h6],
+            products=[ch3, ch3],
+            kinetics=Arrhenius(A=(686.375e6, "1/s"), n=0, Ea=(4000, "kJ/mol"), T0=(1, "K")),
+        )
+
+        T = 300
+        P = 1.0e5
+        core_species = [c2h6, ch3]
+        core_reactions = [rxn]
+
+        rxn_system = SimpleReactor(
+            T,
+            P,
+            initial_mole_fractions={c2h6: 0.5, ch3: 0.5},
+            n_sims=1,
+            termination=[],
+        )
+        rxn_system.initialize_model(core_species, core_reactions, [], [])
+
+        assert rxn_system.kf[0] == 0.0, "the test requires the forward rate to underflow to zero"
+
+        rate_deriv = rxn_system.compute_rate_derivative()
+        assert np.isfinite(rate_deriv).all(), "compute_rate_derivative() produced NaN or inf"
+
+        # The reverse contribution is [products] / Keq, not zero: check the derivative of the
+        # net rate with respect to kf against that analytic value.
+        keq = rxn.get_equilibrium_constant(T)
+        c_c2h6 = rxn_system.core_species_concentrations[core_species.index(c2h6)]
+        c_ch3 = rxn_system.core_species_concentrations[core_species.index(ch3)]
+        expected_flux = c_c2h6 - c_ch3 * c_ch3 / keq
+        assert expected_flux < 0, "the reverse contribution is expected to dominate, not vanish"
+        # compute_rate_derivative() scales its result by the reactor volume before returning.
+        expected = -rxn_system.V * expected_flux
+        assert abs(rate_deriv[core_species.index(c2h6), 0] - expected) <= abs(1e-6 * expected)

From 15cafa053eaad7f515aac097e21681191759bc18 Mon Sep 17 00:00:00 2001
From: Calvin Pieters <calvinpieters@gmail.com>
Date: Thu, 30 Jul 2026 13:58:08 +0300
Subject: [PATCH 4/8] Check the reverse collision limit when reactants lack
 transport data

calculate_coll_limit raises ValueError when a species involved in the
requested direction has no transport data. In the forward direction that
was handled with `continue`, which skips the rest of the loop body for
that condition and therefore also skips the reverse-direction check --
even though the reverse check depends only on the products, which may
have perfectly good transport data.

Fall through instead, so each direction is evaluated on its own merits.
A reaction whose reactants lack transport data is now still checked in
the reverse direction rather than being silently exempted from the
collision limit report.
---
 rmgpy/reaction.py          |  4 ++--
 test/rmgpy/reactionTest.py | 17 +++++++++++++++++
 2 files changed, 19 insertions(+), 2 deletions(-)

diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py
index d3a071c1eb5..3e74cdf38f5 100644
--- a/rmgpy/reaction.py
+++ b/rmgpy/reaction.py
@@ -1756,7 +1756,7 @@ def check_collision_limit_violation(self, t_min, t_max, p_min, p_max):
                 try:
                     limit_f = self.calculate_coll_limit(temp=temp, reverse=False)
                 except ValueError:
-                    continue
+                    pass
                 else:
                     try:
                         kf = self.get_rate_coefficient(temp, pressure)
@@ -1773,7 +1773,7 @@ def check_collision_limit_violation(self, t_min, t_max, p_min, p_max):
                 try:
                     limit_r = self.calculate_coll_limit(temp=temp, reverse=True)
                 except ValueError:
-                    continue
+                    pass
                 else:
                     try:
                         kr = self.generate_reverse_rate_coefficient().get_rate_coefficient(temp, pressure)
diff --git a/test/rmgpy/reactionTest.py b/test/rmgpy/reactionTest.py
index 93f1703e938..ce67bfc43f9 100644
--- a/test/rmgpy/reactionTest.py
+++ b/test/rmgpy/reactionTest.py
@@ -3391,3 +3391,20 @@ def get_rate_coefficient(self, T, P=0, surface_site_density=0, potential=0):
         # k_at_t_max < limit_max, so the forward direction must not be reported. If the surviving
         # rate were compared against the orphaned t_min limit instead, it would be.
         assert [v for v in violators if v[1] == "forward"] == []
+
+    def test_missing_reactant_transport_still_checks_reverse(self):
+        """
+        Missing transport data on a reactant must not suppress the reverse-direction check,
+        which depends only on the products.
+        """
+        rxn = self.make_reaction(A=1e20)
+        original = rxn.reactants[0].transport_data
+        rxn.reactants[0].transport_data = TransportData(sigma=(0, "angstrom"), epsilon=(0, "K"))
+        try:
+            violators = rxn.check_collision_limit_violation(t_min=300, t_max=1500, p_min=1e5, p_max=1e6)
+        finally:
+            rxn.reactants[0].transport_data = original
+
+        directions = [v[1] for v in violators]
+        assert "forward" not in directions
+        assert "reverse" in directions

From a6ff22fc523467b0a3cd7205fb756c392ab1efac Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Thu, 30 Jul 2026 13:58:08 +0300
Subject: [PATCH 5/8] Fix argument documentation and table markup in
 canteramodel

The docstrings for generate_cantera_conditions and
Cantera.generate_conditions documented arguments named T0List, P0List
and V0List, which do not exist; the parameters are Tlist, Plist and
Vlist. They also used a mixture of backticks and single quotes for the
same argument names, and listed the supported reactor types as
unmarked continuation lines.

Correct the names, mark the reactor types up as a list, and give both
tables the closing border row that reStructuredText simple tables
require -- without it docutils rejects the table with "Malformed table.
No bottom table border found." Separate the second table from the
sentence above it with a blank line, since a table has to begin a new
block to be recognised at all.
---
 rmgpy/tools/canteramodel.py | 75 +++++++++++++++++--------------------
 1 file changed, 35 insertions(+), 40 deletions(-)

diff --git a/rmgpy/tools/canteramodel.py b/rmgpy/tools/canteramodel.py
index e488917f658..3c7eae4282e 100644
--- a/rmgpy/tools/canteramodel.py
+++ b/rmgpy/tools/canteramodel.py
@@ -150,29 +150,24 @@ def __str__(self):
 
 def generate_cantera_conditions(reactor_type_list, reaction_time_list, mol_frac_list, surface_mol_frac_list=None, Tlist=None, Plist=None, Vlist=None):
     """
-    Creates a list of cantera conditions from from the arguments provided.
-
-    ======================= ====================================================
-    Argument                Description
-    ======================= ====================================================
-    `reactor_type_list`        A list of strings of the cantera reactor type. List of supported types below:
-        IdealGasReactor: A constant volume, zero-dimensional reactor for ideal gas mixtures
-        IdealGasConstPressureReactor: A homogeneous, constant pressure, zero-dimensional reactor for ideal gas mixtures
-        IdealGasConstPressureTemperatureReactor: A homogenous, constant pressure and constant temperature, zero-dimensional reactor
-                            for ideal gas mixtures (the same as RMG's SimpleReactor)
-
-    `reaction_time_list`      A tuple object giving the ([list of reaction times], units)
-    `mol_frac_list`           A list of molfrac dictionaries with species object keys
-                           and mole fraction values
-    `surface_mol_frac_list`   A list of molfrac dictionaries with surface species object keys
-                           and mole fraction values
-    To specify the system for an ideal gas, you must define 2 of the following 3 parameters:
-    `T0List`                A tuple giving the ([list of initial temperatures], units)
-    'P0List'                A tuple giving the ([list of initial pressures], units)
-    'V0List'                A tuple giving the ([list of initial specific volumes], units)
-
-
-    This saves all the reaction conditions into the Cantera class.
+    Creates a list of cantera conditions from the arguments provided.
+
+    ======================== ====================================================
+    Argument                 Description
+    ======================== ====================================================
+    `reactor_type_list`      A list of strings of the cantera reactor type. List of supported types below:
+                               - IdealGasReactor: A constant volume, zero-dimensional reactor for ideal gas mixtures
+                               - IdealGasConstPressureReactor: A homogeneous, constant pressure, zero-dimensional reactor for ideal gas mixtures
+                               - IdealGasConstPressureTemperatureReactor: A homogenous, constant pressure and constant temperature, zero-dimensional reactor for ideal gas mixtures (same as RMG's SimpleReactor)
+    `reaction_time_list`     A tuple object giving the ([list of reaction times], units)
+    `mol_frac_list`          A list of molfrac dictionaries with species object keys and mole fraction values
+    `surface_mol_frac_list`  A list of molfrac dictionaries with surface species object keys and mole fraction values
+    `Tlist`                  A tuple giving the ([list of initial temperatures], units)
+    `Plist`                  A tuple giving the ([list of initial pressures], units)
+    `Vlist`                  A tuple giving the ([list of initial specific volumes], units)
+    ======================== ====================================================
+
+    Note: To specify the system for an ideal gas, you must define 2 of the following 3 parameters: `Tlist`, `Plist`, `Vlist`
     """
 
     def convert_to_quantity_list(input_list):
@@ -286,23 +281,23 @@ def __init__(self, species_list=None, reaction_list=None, canteraFile='', output
     def generate_conditions(self, reactor_type_list, reaction_time_list, mol_frac_list, surface_mol_frac_list=None, Tlist=None, Plist=None, Vlist=None):
         """
         This saves all the reaction conditions into the Cantera class.
-        ======================= ====================================================
-        Argument                Description
-        ======================= ====================================================
-        `reactor_type_list`        A list of strings of the cantera reactor type. List of supported types below:
-            IdealGasReactor: A constant volume, zero-dimensional reactor for ideal gas mixtures
-            IdealGasConstPressureReactor: A homogeneous, constant pressure, zero-dimensional reactor for ideal gas mixtures
-            IdealGasConstPressureTemperatureReactor: A homogenous, constant pressure and constant temperature, zero-dimensional reactor 
-                                for ideal gas mixtures (the same as RMG's SimpleReactor)
-
-        `reaction_time_list`      A tuple object giving the ([list of reaction times], units)
-        `mol_frac_list`           A list of molfrac dictionaries with species object keys
-                               and mole fraction values
-        `surface_mol_frac_list`   A list of molfrac dictionaries with surface species object keys and mole fraction values
-        To specify the system for an ideal gas, you must define 2 of the following 3 parameters:
-        `T0List`                A tuple giving the ([list of initial temperatures], units)
-        'P0List'                A tuple giving the ([list of initial pressures], units)
-        'V0List'                A tuple giving the ([list of initial specific volumes], units)
+
+        ======================== ====================================================
+        Argument                 Description
+        ======================== ====================================================
+        `reactor_type_list`      A list of strings of the cantera reactor type. List of supported types below:
+                                   - IdealGasReactor: A constant volume, zero-dimensional reactor for ideal gas mixtures
+                                   - IdealGasConstPressureReactor: A homogeneous, constant pressure, zero-dimensional reactor for ideal gas mixtures
+                                   - IdealGasConstPressureTemperatureReactor: A homogenous, constant pressure and constant temperature, zero-dimensional reactor for ideal gas mixtures (same as RMG's SimpleReactor)
+        `reaction_time_list`     A tuple object giving the ([list of reaction times], units)
+        `mol_frac_list`          A list of molfrac dictionaries with species object keys and mole fraction values
+        `surface_mol_frac_list`  A list of molfrac dictionaries with surface species object keys and mole fraction values
+        `Tlist`                  A tuple giving the ([list of initial temperatures], units)
+        `Plist`                  A tuple giving the ([list of initial pressures], units)
+        `Vlist`                  A tuple giving the ([list of initial specific volumes], units)
+        ======================== ====================================================
+
+        Note: To specify the system for an ideal gas, you must define 2 of the following 3 parameters: `Tlist`, `Plist`, `Vlist`
         """
 
         self.conditions = generate_cantera_conditions(reactor_type_list, reaction_time_list, mol_frac_list, surface_mol_frac_list, Tlist, Plist, Vlist)

From 8dff1f3b1df46f39b483b5d5483405e633a7b58f Mon Sep 17 00:00:00 2001
From: Calvin Pieters <calvinpieters@gmail.com>
Date: Fri, 31 Jul 2026 18:01:22 +0300
Subject: [PATCH 6/8] Drop underflowed samples before fitting a reverse rate
 coefficient

Raising on a zero rate coefficient reports the problem at its source, but
it also fails a fit that could have succeeded: a reverse rate coefficient
that underflows at the cold end of the fitting range still has a perfectly
good fit over the rest of the range, and the underflowed samples carry no
information for a fit performed in log space.

Filter those samples out where klist is built, in the five reverse_*_rate
helpers, so the common case produces a usable fit instead of an error. The
validation in fit_to_data stays as a backstop for callers that construct
klist themselves, and is extended to StickingCoefficient.fit_to_data and
SurfaceChargeTransfer.fit_to_data, which perform the same log-space fit and
had the same gap.

Addresses review comments on rmgpy/kinetics/arrhenius.pyx.
---
 rmgpy/kinetics/surface.pyx |  5 +++
 rmgpy/reaction.py          | 64 +++++++++++++++++++++++++++++---------
 test/rmgpy/reactionTest.py | 42 +++++++++++++++++++++----
 3 files changed, 90 insertions(+), 21 deletions(-)

diff --git a/rmgpy/kinetics/surface.pyx b/rmgpy/kinetics/surface.pyx
index 5e1eca4b80b..9802eedde8b 100644
--- a/rmgpy/kinetics/surface.pyx
+++ b/rmgpy/kinetics/surface.pyx
@@ -169,6 +169,9 @@ cdef class StickingCoefficient(KineticsModel):
         """
         import scipy.stats
 
+        if any(klist==0):
+            raise ValueError("Rates must all be nonzero; a zero rate coefficient cannot be fit in log space.")
+
         assert len(Tlist) == len(klist), "length of temperatures and rates must be the same"
         if len(Tlist) < 3 + three_params:
             raise KineticsError('Not enough degrees of freedom to fit this Arrhenius expression')
@@ -973,6 +976,8 @@ cdef class SurfaceChargeTransfer(KineticsModel):
         import scipy.stats
         if not all(np.isfinite(klist)):
             raise  ValueError("Rates must all be finite, not inf or NaN")
+        if any(klist==0):
+            raise ValueError("Rates must all be nonzero; a zero rate coefficient cannot be fit in log space.")
         if any(klist<0):
             if not all(klist<0):
                 raise ValueError("Rates must all be positive or all be negative.")
diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py
index 3e74cdf38f5..95476419cda 100644
--- a/rmgpy/reaction.py
+++ b/rmgpy/reaction.py
@@ -1142,8 +1142,9 @@ def reverse_arrhenius_rate(self, k_forward, reverse_units, Tmin=None, Tmax=None)
         klist = np.zeros_like(Tlist)
         for i in range(len(Tlist)):
             klist[i] = kf.get_rate_coefficient(Tlist[i]) / self.get_equilibrium_constant(Tlist[i])
+        Tfit, kfit = _drop_zero_rate_samples(Tlist, klist)
         kr = Arrhenius()
-        kr.fit_to_data(Tlist, klist, reverse_units, kf.T0.value_si)
+        kr.fit_to_data(Tfit, kfit, reverse_units, kf.T0.value_si)
         kr.solute = kf.solute
         return kr
 
@@ -1166,8 +1167,9 @@ def reverse_surface_arrhenius_rate(self, k_forward, reverse_units, Tmin=None, Tm
         klist = np.zeros_like(Tlist)
         for i in range(len(Tlist)):
             klist[i] = kf.get_rate_coefficient(Tlist[i]) / self.get_equilibrium_constant(Tlist[i])
+        Tfit, kfit = _drop_zero_rate_samples(Tlist, klist)
         kr = SurfaceArrhenius()
-        kr.fit_to_data(Tlist, klist, reverse_units, kf.T0.value_si)
+        kr.fit_to_data(Tfit, kfit, reverse_units, kf.T0.value_si)
         kr.solute = kf.solute
         return kr
 
@@ -1193,8 +1195,9 @@ def reverse_sticking_coeff_rate(self, k_forward, reverse_units, surface_site_den
             klist[i] = \
                 self.get_surface_rate_coefficient(Tlist[i], surface_site_density=surface_site_density) / \
                 self.get_equilibrium_constant(Tlist[i], surface_site_density=surface_site_density)
+        Tfit, kfit = _drop_zero_rate_samples(Tlist, klist)
         kr = SurfaceArrhenius()
-        kr.fit_to_data(Tlist, klist, reverse_units, kf.T0.value_si)
+        kr.fit_to_data(Tfit, kfit, reverse_units, kf.T0.value_si)
         kr.solute = kf.solute
         return kr
 
@@ -1219,8 +1222,9 @@ def reverse_surface_charge_transfer_rate(self, k_forward, reverse_units, Tmin=No
         klist = np.zeros_like(Tlist)
         for i in range(len(Tlist)):
             klist[i] = kf.get_rate_coefficient(Tlist[i],V0) / self.get_equilibrium_constant(Tlist[i],V0)
+        Tfit, kfit = _drop_zero_rate_samples(Tlist, klist)
         kr = SurfaceChargeTransfer(alpha=kf.alpha.value, electrons=-1*self.electrons, V0=(V0,'V'))
-        kr.fit_to_data(Tlist, klist, reverse_units, kf.T0.value_si)
+        kr.fit_to_data(Tfit, kfit, reverse_units, kf.T0.value_si)
         kr.solute = kf.solute
         return kr
 
@@ -1243,8 +1247,9 @@ def reverse_arrhenius_charge_transfer_rate(self, k_forward, reverse_units, Tmin=
         klist = np.zeros_like(Tlist)
         for i in range(len(Tlist)):
             klist[i] = kf.get_rate_coefficient(Tlist[i],V0) / self.get_equilibrium_constant(Tlist[i],V0)
+        Tfit, kfit = _drop_zero_rate_samples(Tlist, klist)
         kr = ArrheniusChargeTransfer(alpha=kf.alpha.value, electrons=-1*self.electrons, V0=(V0,'V'))
-        kr.fit_to_data(Tlist, klist, reverse_units, kf.T0.value_si)
+        kr.fit_to_data(Tfit, kfit, reverse_units, kf.T0.value_si)
         kr.solute = kf.solute
         return kr
 
@@ -1737,7 +1742,9 @@ def check_collision_limit_violation(self, t_min, t_max, p_min, p_max):
         """
         Warn if a core reaction violates the collision limit rate in either the forward or reverse direction
         at the relevant extreme T/P conditions. Assuming a monotonic behaviour of the kinetics.
-        Returns a list with the reaction object and the direction in which the violation was detected.
+        Returns ``(violator_list, skipped)``, where `violator_list` holds the reaction object and the
+        direction in which each violation was detected, and `skipped` counts the direction/condition
+        pairs that could not be evaluated.
         """
         conditions = [[t_min, p_min]]
         if t_min != t_max:
@@ -1750,40 +1757,52 @@ def check_collision_limit_violation(self, t_min, t_max, p_min, p_max):
         violator_list = []
         forward_checks = []
         reverse_checks = []
+        skipped = 0
+        reverse_kinetics = None
+        if len(self.products) >= 2:
+            try:
+                reverse_kinetics = self.generate_reverse_rate_coefficient()
+            except (ReactionError, KineticsError, ZeroDivisionError, OverflowError) as err:
+                logging.warning(
+                    "Skipping reverse collision limit check for reaction %s because the reverse rate "
+                    "coefficient could not be generated: %s",
+                    self, err,
+                )
+                skipped += len(conditions)
         for condition in conditions:
             temp, pressure = condition
             if len(self.reactants) >= 2:
                 try:
                     limit_f = self.calculate_coll_limit(temp=temp, reverse=False)
                 except ValueError:
-                    pass
+                    skipped += 1
                 else:
                     try:
                         kf = self.get_rate_coefficient(temp, pressure)
-                    except (ReactionError, KineticsError, TypeError, ValueError,
-                            ZeroDivisionError, OverflowError) as err:
+                    except (ReactionError, KineticsError, ZeroDivisionError, OverflowError) as err:
                         logging.warning(
                             "Skipping forward collision limit check for reaction %s at %.1f K, %.3g Pa "
                             "because rate evaluation failed: %s",
                             self, temp, pressure, err,
                         )
+                        skipped += 1
                     else:
                         forward_checks.append((kf, limit_f, condition))
-            if len(self.products) >= 2:
+            if len(self.products) >= 2 and reverse_kinetics is not None:
                 try:
                     limit_r = self.calculate_coll_limit(temp=temp, reverse=True)
                 except ValueError:
-                    pass
+                    skipped += 1
                 else:
                     try:
-                        kr = self.generate_reverse_rate_coefficient().get_rate_coefficient(temp, pressure)
-                    except (ReactionError, KineticsError, TypeError, ValueError,
-                            ZeroDivisionError, OverflowError) as err:
+                        kr = reverse_kinetics.get_rate_coefficient(temp, pressure)
+                    except (ReactionError, KineticsError, ZeroDivisionError, OverflowError) as err:
                         logging.warning(
                             "Skipping reverse collision limit check for reaction %s at %.1f K, %.3g Pa "
                             "because reverse rate evaluation failed: %s",
                             self, temp, pressure, err,
                         )
+                        skipped += 1
                     else:
                         reverse_checks.append((kr, limit_r, condition))
         for direction, checks in (('forward', forward_checks), ('reverse', reverse_checks)):
@@ -1792,7 +1811,7 @@ def check_collision_limit_violation(self, t_min, t_max, p_min, p_max):
                     ratio = k / collision_limit
                     condition = '{0} K, {1:.1f} bar'.format(temp, pressure / 1e5)
                     violator_list.append([self, direction, ratio, condition])
-        return violator_list
+        return violator_list, skipped
 
     def calculate_coll_limit(self, temp, reverse=False):
         """
@@ -1854,6 +1873,21 @@ def generate_high_p_limit_kinetics(self):
         """
         raise NotImplementedError("generate_high_p_limit_kinetics is not implemented for all Reaction subclasses.")
 
+def _drop_zero_rate_samples(Tlist, klist):
+    """
+    Return the subset of `Tlist` and `klist` for which the rate coefficient is nonzero.
+
+    Reverse rate coefficients underflow to exactly zero at the cold end of the fitting range for
+    strongly endothermic reactions. Such a sample carries no information for a fit performed in
+    log space, and log(0) would make every fitted parameter NaN, so it is dropped instead.
+    """
+    nonzero = klist != 0
+    if not nonzero.all():
+        logging.debug("Dropping %d of %d rate coefficient samples that underflowed to zero "
+                      "before fitting.", (~nonzero).sum(), len(klist))
+    return Tlist[nonzero], klist[nonzero]
+
+
 def _same_object(object1, object2, _check_identical=False, _only_check_label=False,
              _generate_initial_map=False, _strict=True, _save_order=False):
     if _only_check_label:
diff --git a/test/rmgpy/reactionTest.py b/test/rmgpy/reactionTest.py
index ce67bfc43f9..a0d7b347e1a 100644
--- a/test/rmgpy/reactionTest.py
+++ b/test/rmgpy/reactionTest.py
@@ -57,9 +57,10 @@
     StickingCoefficient,
     SurfaceChargeTransfer,
 )
+from rmgpy.exceptions import KineticsError
 from rmgpy.molecule import Molecule
 from rmgpy.quantity import Quantity
-from rmgpy.reaction import Reaction
+from rmgpy.reaction import Reaction, _drop_zero_rate_samples
 from rmgpy.species import Species, TransitionState
 from rmgpy.statmech.conformer import Conformer
 from rmgpy.statmech.rotation import NonlinearRotor
@@ -3251,6 +3252,31 @@ def test_reverse_surface_charge_transfer_rate(self):
                 assert order_of_magnitude(kf/kr) == order_of_magnitude(K)
 
 
+class TestReverseRateFitWithUnderflow:
+    """
+    Tests that a reverse rate coefficient can still be fitted when it underflows at low T.
+    """
+
+    def test_zero_samples_are_dropped_before_fitting(self):
+        """
+        A reverse rate coefficient that underflows to zero at the cold end of the fitting range
+        must not prevent the fit; the affected samples are dropped instead.
+        """
+        Tlist = np.array([300.0, 500.0, 700.0, 900.0, 1200.0, 1500.0])
+        klist = np.array([0.0, 1e-8, 1e-4, 1e-2, 1.0, 10.0])
+
+        Tfit, kfit = _drop_zero_rate_samples(Tlist, klist)
+        assert len(Tfit) == 5
+        assert 300.0 not in Tfit
+        assert (kfit > 0).all()
+
+        # The surviving samples fit cleanly, where the full list would give NaN parameters.
+        arrhenius = Arrhenius().fit_to_data(Tfit, kfit, kunits="m^3/(mol*s)")
+        assert np.isfinite(arrhenius.A.value_si)
+        assert np.isfinite(arrhenius.n.value_si)
+        assert np.isfinite(arrhenius.Ea.value_si)
+
+
 class TestCollisionLimitViolation:
     """
     Contains unit tests of the Reaction.check_collision_limit_violation() method.
@@ -3326,12 +3352,14 @@ def make_reaction(self, A):
     def test_no_violation_for_physical_rate(self):
         """A sane rate coefficient must not be reported as a collision limit violator."""
         rxn = self.make_reaction(A=1e3)
-        assert rxn.check_collision_limit_violation(t_min=300, t_max=1500, p_min=1e5, p_max=1e6) == []
+        violators, skipped = rxn.check_collision_limit_violation(t_min=300, t_max=1500, p_min=1e5, p_max=1e6)
+        assert violators == []
+        assert skipped == 0
 
     def test_violation_is_detected_and_labelled(self):
         """An absurdly large rate coefficient must be reported in both directions."""
         rxn = self.make_reaction(A=1e20)
-        violators = rxn.check_collision_limit_violation(t_min=300, t_max=1500, p_min=1e5, p_max=1e6)
+        violators, _ = rxn.check_collision_limit_violation(t_min=300, t_max=1500, p_min=1e5, p_max=1e6)
 
         directions = [v[1] for v in violators]
         assert "forward" in directions
@@ -3371,7 +3399,7 @@ class RateFailsAtTmin(Reaction):
 
             def get_rate_coefficient(self, T, P=0, surface_site_density=0, potential=0):
                 if T <= t_min:
-                    raise ValueError("simulated rate evaluation failure at t_min")
+                    raise KineticsError("simulated rate evaluation failure at t_min")
                 return k_at_t_max
 
         rxn = RateFailsAtTmin(
@@ -3380,13 +3408,14 @@ def get_rate_coefficient(self, T, P=0, surface_site_density=0, potential=0):
             kinetics=Arrhenius(A=(1e3, "m^3/(mol*s)"), n=0, Ea=(0, "kcal/mol"), T0=(1, "K")),
         )
         with caplog.at_level(logging.WARNING):
-            violators = rxn.check_collision_limit_violation(
+            violators, skipped = rxn.check_collision_limit_violation(
                 t_min=t_min, t_max=t_max, p_min=pressure, p_max=pressure
             )
 
         # Guard the premise: if the override ever stopped reaching the compiled caller, the
         # assertion below would pass vacuously because the real rate is also under the limit.
         assert "Skipping forward collision limit check" in caplog.text
+        assert skipped == 1
 
         # k_at_t_max < limit_max, so the forward direction must not be reported. If the surviving
         # rate were compared against the orphaned t_min limit instead, it would be.
@@ -3401,10 +3430,11 @@ def test_missing_reactant_transport_still_checks_reverse(self):
         original = rxn.reactants[0].transport_data
         rxn.reactants[0].transport_data = TransportData(sigma=(0, "angstrom"), epsilon=(0, "K"))
         try:
-            violators = rxn.check_collision_limit_violation(t_min=300, t_max=1500, p_min=1e5, p_max=1e6)
+            violators, skipped = rxn.check_collision_limit_violation(t_min=300, t_max=1500, p_min=1e5, p_max=1e6)
         finally:
             rxn.reactants[0].transport_data = original
 
         directions = [v[1] for v in violators]
         assert "forward" not in directions
         assert "reverse" in directions
+        assert skipped == 2  # forward direction skipped at both conditions

From 5a02e7a22f74db0791c93af4cca4eea7611bc405 Mon Sep 17 00:00:00 2001
From: Calvin Pieters <calvinpieters@gmail.com>
Date: Fri, 31 Jul 2026 18:01:22 +0300
Subject: [PATCH 7/8] Mark irreversible reactions with Keq = inf in every
 reactor

SimpleReactor marks an irreversible reaction with kb = 0 and Keq = inf,
but the liquid, surface and mass-flow-batch reactors left Keq at the 0.0
from np.zeros_like. Consumers therefore had to accept either sentinel:
compute_rate_derivative's `not (Keq[j] > 0)` guard was written to catch
the zero case for that reason.

Set kb = 0 and Keq = inf explicitly in the other three reactors so the
sentinel agrees by design. The guard is unchanged and still correct, but
now the only thing that can trip it is NaN thermochemistry -- a reversible
reaction cannot have Keq == 0, since get_equilibrium_constant raises in
that case. Say so in the comment.

Addresses review comment on rmgpy/solver/base.pyx.
---
 rmgpy/solver/base.pyx      | 6 ++++++
 rmgpy/solver/liquid.pyx    | 3 +++
 rmgpy/solver/mbSampled.pyx | 3 +++
 rmgpy/solver/surface.pyx   | 3 +++
 4 files changed, 15 insertions(+)

diff --git a/rmgpy/solver/base.pyx b/rmgpy/solver/base.pyx
index dca598bb3df..761705d3b31 100644
--- a/rmgpy/solver/base.pyx
+++ b/rmgpy/solver/base.pyx
@@ -1355,6 +1355,12 @@ cdef class ReactionSystem(DASx):
             else:  # three reactants
                 fderiv = C[ir[j, 0]] * C[ir[j, 1]] * C[ir[j, 2]]
 
+            # kb is always built as kf/Keq, so kb/kf is identically 1/Keq. Dividing by Keq avoids
+            # the 0.0/0.0 that kb/kf becomes when kf underflows, which is NaN in C and silently
+            # poisons the sensitivity matrix. Every reactor marks an irreversible reaction with
+            # Keq = inf, which passes this test and correctly gives C/inf = 0. A reversible
+            # reaction cannot have Keq == 0 (get_equilibrium_constant raises), so the only thing
+            # that trips this branch is NaN thermochemistry, where zero is the safe answer.
             if not (Keq[j] > 0.0):
                 rderiv = 0.0
             elif ip[j, 1] == -1:  # only one product
diff --git a/rmgpy/solver/liquid.pyx b/rmgpy/solver/liquid.pyx
index aefd38e030d..4095861738f 100644
--- a/rmgpy/solver/liquid.pyx
+++ b/rmgpy/solver/liquid.pyx
@@ -166,6 +166,9 @@ cdef class LiquidReactor(ReactionSystem):
             if rxn.reversible:
                 self.Keq[j] = rxn.get_equilibrium_constant(self.T.value_si)
                 self.kb[j] = self.kf[j] / self.Keq[j]
+            else:
+                self.kb[j] = 0.0
+                self.Keq[j] = np.inf
 
     def get_threshold_rate_constants(self, model_settings):
         """
diff --git a/rmgpy/solver/mbSampled.pyx b/rmgpy/solver/mbSampled.pyx
index 18d76381400..572dd130531 100644
--- a/rmgpy/solver/mbSampled.pyx
+++ b/rmgpy/solver/mbSampled.pyx
@@ -237,6 +237,9 @@ cdef class MBSampledReactor(ReactionSystem):
             if rxn.reversible:
                 self.Keq[j] = rxn.get_equilibrium_constant(self.T.value_si)
                 self.kb[j] = self.kf[j] / self.Keq[j]
+            else:
+                self.kb[j] = 0.0
+                self.Keq[j] = np.inf
 
     def set_colliders(self, core_reactions, edge_reactions, core_species):
         """
diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx
index 8fa1fb7ec7c..98c9f7008e4 100644
--- a/rmgpy/solver/surface.pyx
+++ b/rmgpy/solver/surface.pyx
@@ -313,6 +313,9 @@ cdef class SurfaceReactor(ReactionSystem):
                 # which applies the coverage-dependent correction to Keq at runtime.
                 self.Keq[j] = rxn.get_equilibrium_constant(self.T.value_si)
                 self.kb[j] = self.kf[j] / self.Keq[j]
+            else:
+                self.kb[j] = 0.0
+                self.Keq[j] = np.inf
 
     def log_initial_conditions(self, number=None):
         """

From 1f3f0c7c450483c56bf3f0887e7016d2ea5b79f3 Mon Sep 17 00:00:00 2001
From: Calvin Pieters <calvinpieters@gmail.com>
Date: Fri, 31 Jul 2026 18:01:22 +0300
Subject: [PATCH 8/8] Report how many collision limit checks were skipped

Three changes to check_collision_limit_violation, all from review:

Narrow the caught exceptions to ReactionError, KineticsError,
ZeroDivisionError and OverflowError. TypeError and ValueError indicate
programmatic bugs rather than a reaction whose kinetics cannot be
evaluated, and swallowing them here would hide real defects.

Generate the reverse rate coefficient once per reaction rather than once
per condition. It does not depend on T or P, so it was being rebuilt up to
four times, and a reaction whose reverse coefficient cannot be generated
was reported once per condition instead of once.

Return the number of skipped direction/condition pairs alongside the
violator list, and count skipped reactions in check_model, so that a model
where nothing could be evaluated no longer reports "No collision rate
violators found in the model's core" as though it had been fully checked.
---
 rmgpy/reaction.pxd         |  2 +-
 rmgpy/rmg/main.py          | 13 ++++++++++++-
 test/rmgpy/rmg/mainTest.py |  2 +-
 3 files changed, 14 insertions(+), 3 deletions(-)

diff --git a/rmgpy/reaction.pxd b/rmgpy/reaction.pxd
index e17b74c2e89..a2fddc96271 100644
--- a/rmgpy/reaction.pxd
+++ b/rmgpy/reaction.pxd
@@ -147,7 +147,7 @@ cdef class Reaction:
 
     cpdef ensure_species(self, bint reactant_resonance=?, bint product_resonance=?, bint save_order=?)
 
-    cpdef list check_collision_limit_violation(self, float t_min, float t_max, float p_min, float p_max)
+    cpdef tuple check_collision_limit_violation(self, float t_min, float t_max, float p_min, float p_max)
 
     cpdef calculate_coll_limit(self, float temp, bint reverse=?)
 
diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py
index 21815c5e535..ae1d350d7ef 100644
--- a/rmgpy/rmg/main.py
+++ b/rmgpy/rmg/main.py
@@ -1643,12 +1643,14 @@ def check_model(self):
 
         # Check all core reactions (in both directions) for collision limit violation
         violators = []
+        skipped_checks = 0
+        skipped_reactions = 0
         for rxn in self.reaction_model.core.reactions:
             if rxn.is_surface_reaction():
                 # Don't check collision limits for surface reactions.
                 continue
             try:
-                violator_list = rxn.check_collision_limit_violation(
+                violator_list, skipped = rxn.check_collision_limit_violation(
                     t_min=self.Tmin, t_max=self.Tmax, p_min=self.Pmin, p_max=self.Pmax
                 )
             except Exception:
@@ -1656,7 +1658,9 @@ def check_model(self):
                     "Skipping collision limit check for reaction %s because evaluation failed.",
                     rxn, exc_info=True,
                 )
+                skipped_reactions += 1
                 continue
+            skipped_checks += skipped
             if violator_list:
                 violators.extend(violator_list)
         # Whether or not violators were found, rename 'collision_rate_violators.log' if it exists
@@ -1691,6 +1695,13 @@ def check_model(self):
                     violators_f.write(
                         f"{rxn_string}\n" f"Direction: {direction}\n" f"Violation factor: {ratio:.2f}\n" f"Violation condition: {condition}\n\n\n"
                     )
+        elif skipped_checks or skipped_reactions:
+            logging.info(
+                "No collision rate violators found among the checks that could be evaluated. "
+                "%d individual checks and %d whole reactions were skipped because their rate "
+                "coefficients could not be evaluated; see the warnings above.",
+                skipped_checks, skipped_reactions,
+            )
         else:
             logging.info("No collision rate violators found in the model's core.")
 
diff --git a/test/rmgpy/rmg/mainTest.py b/test/rmgpy/rmg/mainTest.py
index 697d000b056..8f40563eee9 100644
--- a/test/rmgpy/rmg/mainTest.py
+++ b/test/rmgpy/rmg/mainTest.py
@@ -638,7 +638,7 @@ def check_collision_limit_violation(self, t_min, t_max, p_min, p_max):
 
         class ViolatingReaction(Reaction):
             def check_collision_limit_violation(self, t_min, t_max, p_min, p_max):
-                return [[self, "forward", 12.5, "1500.0 K, 1.0 bar"]]
+                return [[self, "forward", 12.5, "1500.0 K, 1.0 bar"]], 0
 
         exploding = self.make_reaction(ExplodingReaction)
         violating = self.make_reaction(ViolatingReaction)