diff --git a/rmgpy/pdep/msc.pyx b/rmgpy/pdep/msc.pyx index 38aced47e7b..8f909bbd35a 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/rmgpy/reaction.py b/rmgpy/reaction.py index 13817852863..740af735b22 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/pdep/mscTest.py b/test/rmgpy/pdep/mscTest.py new file mode 100644 index 00000000000..a9f23e3beff --- /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") diff --git a/test/rmgpy/reactionTest.py b/test/rmgpy/reactionTest.py index 458de09f6a1..1b8d06efa08 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,