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/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.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/reaction.py b/rmgpy/reaction.py index 13817852863..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: @@ -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: diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index eaa4a732eba..ae1d350d7ef 100644 --- a/rmgpy/rmg/main.py +++ b/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.") diff --git a/rmgpy/solver/base.pyx b/rmgpy/solver/base.pyx index ad088aed9af..761705d3b31 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,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 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): """ 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) 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 diff --git a/test/rmgpy/reactionTest.py b/test/rmgpy/reactionTest.py index 458de09f6a1..a0d7b347e1a 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 @@ -56,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 @@ -66,6 +68,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 +3250,191 @@ 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 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. + """ + + 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) + 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) + + 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 KineticsError("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, 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. + 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, 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 diff --git a/test/rmgpy/rmg/mainTest.py b/test/rmgpy/rmg/mainTest.py index 00aa547e9f8..8f40563eee9 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"]], 0 + + 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() 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)