Reject zero rate coefficients and improve collision limit checks

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.")

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.")

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=?)
 

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:
@@ -1748,38 +1755,63 @@ 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 = []
+        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:
-                    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
+                    skipped += 1
                 else:
-                    kf_list.append(self.get_rate_coefficient(condition[0], condition[1]))
-            if len(self.products) >= 2:
+                    try:
+                        kf = self.get_rate_coefficient(temp, pressure)
+                    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 and reverse_kinetics is not None:
                 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
+                    skipped += 1
                 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])
-        return violator_list
+                    try:
+                        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)):
+            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, skipped
 
     def calculate_coll_limit(self, temp, reverse=False):
         """
@@ -1841,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:

rmgpy/rmg/main.py

@@ -1643,11 +1643,24 @@ 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
-            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, skipped = 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,
+                )
+                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
@@ -1682,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.")
 

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,20 @@ 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]]
+            # 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
+                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

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):
         """

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):
         """

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):
         """