Diagnose mismatched conformer degrees of freedom, instead of surfacing them as a bare "Singular matrix"

rmgpy/pdep/msc.pyx

@@ -144,7 +144,35 @@ cpdef apply_modified_strong_collision_method(network, str efficiency_model='defa
                         logging.debug(str([f_im[j, n - n_isom, r, s], dens_states[n, r, s],
                                            (2 * j_list[s] + 1), exp(-e_list[r] * beta), e_list[r], beta]))
             # Solve for steady-state population
-            x = -np.linalg.solve(a_mat, b)
+            try:
+                x = -np.linalg.solve(a_mat, b)
+            except np.linalg.LinAlgError:
+                # a_mat goes singular when the microcanonical rate coefficients in this grain
+                # overwhelm the collision frequency, so the collisional terms on the diagonal stop
+                # carrying anything the solve can use. That is almost always a sign that the input
+                # k(E) is unphysical rather than that the algorithm was unlucky -- most often
+                # because the wells and the transition states of this network are described with
+                # different degrees of freedom, so Q_TS/Q_reactant never cancels. Report where it
+                # happened and how far apart the two scales were: the bare LinAlgError carries none
+                # of this, and locating it otherwise means bisecting the grains by hand.
+                # `.max()` on an empty array raises, and a network with no reactant or product
+                # channels has a zero-sized g_nj -- so an unguarded max() here would replace the
+                # unhelpful LinAlgError with an equally unhelpful ValueError.
+                k_max = 0.0
+                if k_ij.size:
+                    k_max = max(k_max, np.abs(k_ij[:, :, r, s]).max())
+                if g_nj.size:
+                    k_max = max(k_max, np.abs(g_nj[:, :, r, s]).max())
+                freq_max = coll_freq.max() if coll_freq.size else 0.0
+                raise ModifiedStrongCollisionError(
+                    'Singular matrix encountered while solving for the steady-state population of '
+                    'network {0} at grain {1} (E = {2:g} kJ/mol), J index {3}, T = {4:g} K. The '
+                    'largest rate coefficient in this grain is {5:g} s^-1 against a collision '
+                    'frequency of {6:g} s^-1, a ratio of {7:g}. A ratio this large usually means '
+                    'k(E) is unphysical -- check that every isomer and transition state in the '
+                    'network enumerates the same translational and rotational modes.'.format(
+                        network.label, r, e_list[r] / 1000., j_list[s], temperature,
+                        k_max, freq_max, k_max / freq_max if freq_max else float('inf')))
             for n in range(n_isom + n_reac):
                 for i in range(n_isom):
                     pa[i, n, r, s] = x[i, n]

rmgpy/reaction.py

@@ -61,10 +61,32 @@
 from rmgpy.molecule.molecule import Molecule, Atom
 from rmgpy.pdep.reaction import calculate_microcanonical_rate_coefficient
 from rmgpy.species import Species
+from rmgpy.statmech.rotation import Rotation
+from rmgpy.statmech.translation import Translation
 from rmgpy.thermo import ThermoData
 
 ################################################################################
 
+# Reactions whose degrees-of-freedom mismatch has already been reported. The check runs per
+# calculate_tst_rate_coefficient() call, and that is called once per temperature, so without this
+# a single inconsistent reaction would emit one identical warning per point in Tlist.
+_DOF_MISMATCH_WARNED = set()
+
+
+def _has_external_modes(conformer):
+    """
+    Return ``True`` if `conformer` enumerates any external (translational or rotational) mode.
+
+    Deliberately a yes/no answer rather than a set of mode classes: a monatomic species
+    legitimately carries translation and no rotation, and a linear species carries a different
+    rotor from a nonlinear one, so comparing mode classes between a transition state and its
+    reactants reports differences that are entirely correct. What is never correct is one side
+    of the TST expression describing the external degrees of freedom and the other omitting them
+    altogether -- that is the mismatch which leaves Q_TS/Q_reactant uncancelled.
+    """
+    return any(isinstance(mode, (Translation, Rotation)) for mode in conformer.modes)
+
+
 # helper function for sorting
 def get_sorting_key(spc):
     # List of elements to sort by, order is intentional
@@ -1411,6 +1433,8 @@ def calculate_tst_rate_coefficient(self, T):
         is the Planck constant. :math:`\\kappa(T)` is an optional tunneling
         correction.
         """
+        self._warn_if_dof_inconsistent()
+
         # Determine TST rate constant at each temperature
         Qreac = 1.0
         E0 = 0.0
@@ -1428,6 +1452,44 @@ def calculate_tst_rate_coefficient(self, T):
 
         return k
 
+    def _warn_if_dof_inconsistent(self):
+        """
+        Warn if the transition state and the reactants carry different classes of
+        external degrees of freedom.
+
+        The TST expression is only meaningful when :math:`Q^\\ddagger` and
+        :math:`Q^\\mathrm{A} Q^\\mathrm{B}` count the same modes, so that the
+        translational and rotational contributions cancel. A transition state
+        described by a full conformer (translation + rotation + vibration) sitting
+        above wells described as vibration-only -- the usual result of splicing a
+        quantum-chemistry saddle point into an otherwise estimated network -- leaves
+        those contributions uncancelled and can inflate the rate coefficient by many
+        orders of magnitude. Nothing downstream can recover the intent, so warn here.
+
+        This only warns: a mixed description can be deliberate, and raising would
+        break working code.
+        """
+        if self.transition_state is None or self.transition_state.conformer is None:
+            return
+        ts_has_external = _has_external_modes(self.transition_state.conformer)
+        for spec in self.reactants:
+            if spec.conformer is None:
+                continue
+            if _has_external_modes(spec.conformer) == ts_has_external:
+                continue
+            key = (str(self), self.transition_state.label, spec.label, ts_has_external)
+            if key in _DOF_MISMATCH_WARNED:
+                continue
+            _DOF_MISMATCH_WARNED.add(key)
+            with_external, without_external = (
+                (self.transition_state.label, spec.label) if ts_has_external
+                else (spec.label, self.transition_state.label))
+            logging.warning(
+                '%r enumerates translational/rotational modes but %r does not. The TST expression '
+                'assumes these cancel between Q_TS and Q_reactant, so where one side omits them '
+                'k(T) can be wrong by many orders of magnitude. Describe every well and transition '
+                'state in this network the same way.', with_external, without_external)
+
     def can_tst(self):
         """
         Return ``True`` if the necessary parameters are available for using

test/rmgpy/pdep/mscTest.py

@@ -0,0 +1,133 @@
+#!/usr/bin/env python3
+
+###############################################################################
+#                                                                             #
+# RMG - Reaction Mechanism Generator                                          #
+#                                                                             #
+# Copyright (c) 2002-2026 Prof. William H. Green (whgreen@mit.edu),           #
+# Prof. Richard H. West (r.west@neu.edu) and the RMG Team (rmg_dev@mit.edu)   #
+#                                                                             #
+# Permission is hereby granted, free of charge, to any person obtaining a     #
+# copy of this software and associated documentation files (the 'Software'),  #
+# to deal in the Software without restriction, including without limitation   #
+# the rights to use, copy, modify, merge, publish, distribute, sublicense,    #
+# and/or sell copies of the Software, and to permit persons to whom the       #
+# Software is furnished to do so, subject to the following conditions:        #
+#                                                                             #
+# The above copyright notice and this permission notice shall be included in  #
+# all copies or substantial portions of the Software.                         #
+#                                                                             #
+# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR  #
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,    #
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE #
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER      #
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING     #
+# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER         #
+# DEALINGS IN THE SOFTWARE.                                                   #
+#                                                                             #
+###############################################################################
+
+"""
+Tests for the diagnostics around the modified strong collision solve.
+"""
+
+import numpy as np
+import pytest
+
+from rmgpy.exceptions import ModifiedStrongCollisionError
+from rmgpy.pdep.msc import apply_modified_strong_collision_method
+
+
+class _EnergyTransferModel:
+    def get_alpha(self, T):
+        return 1000.0
+
+
+class _Species:
+    def __init__(self):
+        self.energy_transfer_model = _EnergyTransferModel()
+
+
+class _Isomer:
+    def __init__(self):
+        self.species = [_Species()]
+
+
+class _StubNetwork:
+    """
+    The smallest network that drives ``a_mat`` singular.
+
+    Two isomers whose isomerization rate coefficients dwarf the collision frequency make the
+    collisional terms on the diagonal negligible, so the matrix reduces to
+    ``[[-k21, k12], [k21, -k12]]`` -- determinant zero, exactly rank-1. That is the shape the
+    solve fails on in practice, when k(E) has been inflated by a network whose wells and
+    transition states are described with different degrees of freedom.
+    """
+
+    label = "stub"
+
+    def __init__(self, k_isomerization=1.0e27, coll_freq=1.0e07):
+        n_isom, n_grains, n_j = 2, 4, 1
+        self.T = 1000.0
+        self.P = 1.0e05
+        self.e_list = np.linspace(0.0, 3.0e04, n_grains)
+        self.j_list = np.array([0], dtype=np.int_)
+        self.dens_states = np.ones((n_isom, n_grains, n_j))
+        self.coll_freq = np.array([coll_freq, coll_freq])
+        self.Kij = np.zeros((n_isom, n_isom, n_grains, n_j))
+        self.Kij[1, 0, :, :] = k_isomerization
+        self.Kij[0, 1, :, :] = k_isomerization
+        self.Fim = np.zeros((n_isom, 0, n_grains, n_j))
+        self.Gnj = np.zeros((0, n_isom, n_grains, n_j))
+        self.E0 = np.zeros(n_isom)
+        self.n_isom = n_isom
+        self.n_reac = 0
+        self.n_prod = 0
+        self.n_grains = n_grains
+        self.n_j = n_j
+        self.isomers = [_Isomer() for _ in range(n_isom)]
+        self.reactants = []
+        self.products = []
+
+
+class TestModifiedStrongCollisionDiagnostics:
+    def test_singular_matrix_reports_where_and_why(self):
+        """
+        A singular solve must say which grain, J index and temperature it failed at, and how far
+        k(E) had outrun the collision frequency. Bare ``numpy.linalg.LinAlgError: Singular
+        matrix`` carries none of that, and locating it otherwise means bisecting the grains by
+        hand.
+        """
+        with pytest.raises(ModifiedStrongCollisionError) as exc_info:
+            apply_modified_strong_collision_method(_StubNetwork(), efficiency_model="none")
+
+        message = str(exc_info.value)
+        assert "Singular matrix" in message
+        assert "grain 0" in message
+        assert "J index 0" in message
+        assert "T = 1000 K" in message
+        # The two scales, and the ratio between them -- the number that says the input was
+        # unphysical rather than the algorithm unlucky.
+        assert "1e+27" in message
+        assert "1e+07" in message
+        assert "1e+20" in message
+
+    def test_diagnostic_survives_a_network_with_no_reactant_channels(self):
+        """
+        ``Gnj`` is zero-sized when a network has no reactant or product channels, and calling
+        ``.max()`` on an empty array raises. The diagnostic must not trade one confusing error
+        for another.
+        """
+        network = _StubNetwork()
+        assert network.Gnj.size == 0
+
+        with pytest.raises(ModifiedStrongCollisionError):
+            apply_modified_strong_collision_method(network, efficiency_model="none")
+
+    def test_well_conditioned_network_is_untouched(self):
+        """
+        The diagnostic must only fire on the path that was already failing: a network whose rate
+        coefficients sit below the collision frequency solves as before.
+        """
+        network = _StubNetwork(k_isomerization=1.0e03, coll_freq=1.0e07)
+        apply_modified_strong_collision_method(network, efficiency_model="none")

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
@@ -1706,6 +1707,90 @@ def test_generate_reverse_rate_coefficient_troe(self):
             krevrev = reverse_reverse_kinetics.get_rate_coefficient(T, P)
             assert round(abs(korig / krevrev - 1.0), 0) == 0
 
+    def _vibration_only(self, reaction, label):
+        """Strip `label`'s conformer down to vibrations, as an estimated well would be described."""
+        from rmgpy.statmech.torsion import Torsion
+        from rmgpy.statmech.vibration import Vibration
+
+        for spec in reaction.reactants:
+            if spec.label == label:
+                spec.conformer.modes = [
+                    mode for mode in spec.conformer.modes if isinstance(mode, (Vibration, Torsion))
+                ]
+                return
+        raise AssertionError("no reactant labelled {0!r}".format(label))
+
+    def test_tst_warns_when_transition_state_has_external_modes_and_reactant_does_not(self, caplog):
+        """
+        A transition state carrying translation and rotation over a well described as
+        vibration-only should warn: the TST expression assumes those cancel between Q_TS and
+        Q_reactant, and here they do not.
+        """
+        from rmgpy.reaction import _DOF_MISMATCH_WARNED
+
+        reaction = deepcopy(self.reaction)
+        self._vibration_only(reaction, "C2H4")
+
+        _DOF_MISMATCH_WARNED.clear()
+        with caplog.at_level(logging.WARNING, logger="root"):
+            reaction.calculate_tst_rate_coefficient(1000.0)
+
+        messages = [record.getMessage() for record in caplog.records]
+        assert any("translational/rotational modes" in message for message in messages), messages
+
+    def test_tst_does_not_warn_when_degrees_of_freedom_match(self, caplog):
+        """
+        The consistent case must stay silent -- a warning on every well-formed network is noise,
+        and noise is what would stop the mismatched case from being noticed.
+        """
+        from rmgpy.reaction import _DOF_MISMATCH_WARNED
+
+        _DOF_MISMATCH_WARNED.clear()
+        with caplog.at_level(logging.WARNING, logger="root"):
+            self.reaction.calculate_tst_rate_coefficient(1000.0)
+
+        messages = [record.getMessage() for record in caplog.records]
+        assert not any("translational/rotational modes" in message for message in messages), messages
+
+    def test_tst_does_not_warn_for_a_monatomic_reactant(self, caplog):
+        """
+        An atom has translation and no rotation, and a polyatomic transition state has both.
+        That difference is physics, not a defect, so it must not warn -- this reaction has an
+        H atom as a reactant and is exactly the case a mode-class comparison gets wrong.
+        """
+        from rmgpy.reaction import _DOF_MISMATCH_WARNED
+        from rmgpy.statmech.rotation import Rotation
+
+        hydrogen = next(spec for spec in self.reaction.reactants if spec.label == "H")
+        assert not any(isinstance(mode, Rotation) for mode in hydrogen.conformer.modes)
+
+        _DOF_MISMATCH_WARNED.clear()
+        with caplog.at_level(logging.WARNING, logger="root"):
+            self.reaction.calculate_tst_rate_coefficient(1000.0)
+
+        messages = [record.getMessage() for record in caplog.records]
+        assert not any("translational/rotational modes" in message for message in messages), messages
+
+    def test_tst_dof_warning_is_not_repeated_per_temperature(self, caplog):
+        """
+        The check runs once per calculate_tst_rate_coefficient() call, and that is called once
+        per temperature -- so a mismatched reaction must still warn only once.
+        """
+        from rmgpy.reaction import _DOF_MISMATCH_WARNED
+
+        reaction = deepcopy(self.reaction)
+        self._vibration_only(reaction, "C2H4")
+
+        _DOF_MISMATCH_WARNED.clear()
+        with caplog.at_level(logging.WARNING, logger="root"):
+            for T in (500.0, 1000.0, 1500.0, 2000.0):
+                reaction.calculate_tst_rate_coefficient(T)
+
+        warnings = [
+            record for record in caplog.records if "translational/rotational modes" in record.getMessage()
+        ]
+        assert len(warnings) == 1, [w.getMessage() for w in warnings]
+
     def test_tst_calculation(self):
         """
         A test of the transition state theory k(T) calculation function,