From 8bcd590f756df5de6b3827d2ac408cd2e3d7f0f8 Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Tue, 25 Aug 2026 09:10:41 +0300 Subject: [PATCH 1/2] Warn when TST partition functions count different degrees of freedom The TST expression only means anything when Q_TS and Q_reactant enumerate the same external modes, so translation and rotation cancel. Splice a quantum-chemistry transition state (translation + rotation + vibration) into a network whose wells are estimated as vibration-only and they no longer do: calculate_tst_rate_coefficient just multiplies the partition functions, with nothing looking at what modes either side carries. In the case that prompted this, Q_TS/Q_reactant was off by ~1e14 and k(3200 K) came out at 1.8e26 s^-1 for a unimolecular step. Warn rather than raise -- a mixed description can be deliberate, and raising would break working code. The check is a yes/no on whether external modes are present at all, not a comparison of mode classes: a monatomic species carries translation and no rotation, and a linear rotor differs from a nonlinear one, so comparing classes reports differences that are entirely correct. Warn once per reaction, since this is called once per temperature in Tlist. Refs #3003 --- rmgpy/reaction.py | 62 +++++++++++++++++++++++++++ test/rmgpy/reactionTest.py | 85 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 147 insertions(+) diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index 1381785286..740af735b2 100644 --- a/rmgpy/reaction.py +++ b/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 diff --git a/test/rmgpy/reactionTest.py b/test/rmgpy/reactionTest.py index 458de09f6a..1b8d06efa0 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 @@ -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, From fa7cf857f6ca2e2e50d35df7db055627521d686c Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Tue, 25 Aug 2026 09:10:41 +0300 Subject: [PATCH 2/2] Say where and why the modified strong collision solve went singular a_mat goes singular when the microcanonical rate coefficients in a grain overwhelm the collision frequency: the collisional terms on the diagonal stop carrying anything the solve can use, and for two isomers the matrix reduces to [[-k21, k12], [k21, -k12]], which is exactly rank-1. All that reaches the user is numpy.linalg.LinAlgError: Singular matrix -- no grain, no J, no temperature. Locating it means bisecting the grains by hand. Report the coordinates and the two scales that produced the failure, so the message says whether the input was unphysical rather than the algorithm unlucky. The intent was already here: b is checked for NaN a few lines above and logs a warning naming the grain and the isomer. a_mat had no equivalent, and it is the argument that fails. Deliberately does not regularize or perturb a_mat. It is singular because the physics handed to it was wrong, and smoothing that over would turn a loud failure into a quiet wrong number. Refs #3003 --- rmgpy/pdep/msc.pyx | 30 ++++++++- test/rmgpy/pdep/mscTest.py | 133 +++++++++++++++++++++++++++++++++++++ 2 files changed, 162 insertions(+), 1 deletion(-) create mode 100644 test/rmgpy/pdep/mscTest.py diff --git a/rmgpy/pdep/msc.pyx b/rmgpy/pdep/msc.pyx index 38aced47e7..8f909bbd35 100644 --- a/rmgpy/pdep/msc.pyx +++ b/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] diff --git a/test/rmgpy/pdep/mscTest.py b/test/rmgpy/pdep/mscTest.py new file mode 100644 index 0000000000..a9f23e3bef --- /dev/null +++ b/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")