Diagnose mismatched conformer degrees of freedom, instead of surfacing them as a bare "Singular matrix" - #3004
Diagnose mismatched conformer degrees of freedom, instead of surfacing them as a bare "Singular matrix"#3004alongd wants to merge 2 commits into
Conversation
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
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
Regression Testing Results
Detailed regression test results.Regression test aromatics:Reference: Execution time (DD:HH:MM:SS): 00:00:00:52 aromatics Passed Core Comparison ✅Original model has 15 species. aromatics Failed Edge Comparison ❌Original model has 106 species. Non-identical thermo! ❌
thermo: Thermo group additivity estimation: group(Cs-(Cds-Cds)CsCsH) + group(Cs-(Cds-Cds)CsCsH) + group(Cs-(Cds-Cds)(Cds-Cds)CsH) + group(Cs-(Cds-Cds)CsHH) + group(Cds-CdsCsCs) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + polycyclic(s2_4_5_diene_1_5) + polycyclic(s3_4_5_ene_3) + polycyclic(s2_5_5_diene_1_5) - ring(Cyclobutene) - ring(Cyclopentene) - ring(Cyclopentene) + radical(cyclopentene-allyl) Non-identical thermo! ❌
thermo: Thermo group additivity estimation: group(Cs-(Cds-Cds)CsCsH) + group(Cs-(Cds-Cds)(Cds-Cds)CsH) + group(Cs-(Cds-Cds)(Cds-Cds)CsH) + group(Cs-CsCsHH) + group(Cds-CdsCsCs) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + Estimated bicyclic component: polycyclic(s2_3_5_ane) - ring(Cyclopentane) - ring(Cyclopropane) + ring(Cyclopentene) + ring(Cyclopropane) + polycyclic(s2_3_6_diene_0_3) + Estimated bicyclic component: polycyclic(s3_5_6_ane) - ring(Cyclohexane) - ring(Cyclopentane) + ring(1,4-Cyclohexadiene) + ring(Cyclopentene) - ring(Cyclopropane) - ring(Cyclopentene) - ring(1,4-Cyclohexadiene) + radical(cyclopentene-4) Non-identical thermo! ❌
thermo: Thermo group additivity estimation: group(Cs-CsCsCsH) + group(Cs-(Cds-Cds)CsCsH) + group(Cs-(Cds-Cds)CsCsH) + group(Cs-(Cds-Cds)CsHH) + group(Cds- Cds(Cds-Cds)Cs) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-Cds(Cds-Cds)H) + polycyclic(s2_3_5_ene_1) + polycyclic(s2_3_6_ene_1) + Estimated bicyclic component: polycyclic(s3_5_6_ane) - ring(Cyclohexane) - ring(Cyclopentane) + ring(Cyclohexene) + ring(Cyclopentene) - ring(Cyclopropane) - ring(Cyclopentene) - ring(Cyclohexene) + radical(cyclopentene-allyl) Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Errors occurred during edge comparison
ERROR conda.cli.main_run:execute(148): `conda run python scripts/checkModels.py aromatics-edge stable_regression_results/aromatics/chemkin/chem_edge_annotated.inp stable_regression_results/aromatics/chemkin/species_edge_dictionary.txt test/regression/aromatics/chemkin/chem_edge_annotated.inp test/regression/aromatics/chemkin/species_edge_dictionary.txt` failed. (See above for error)
|
| k(1bar) | 300K | 400K | 500K | 600K | 800K | 1000K | 1500K | 2000K |
|---|---|---|---|---|---|---|---|---|
| k(T): | 8.02 | 7.64 | 7.35 | 7.11 | 6.75 | 6.48 | 5.99 | 5.64 |
| k(T): | 3.54 | 4.28 | 4.73 | 5.02 | 5.39 | 5.62 | 5.91 | 6.06 |
kinetics: Arrhenius(A=(3.18266e+20,'cm^3/(mol*s)'), n=-2.694, Ea=(-0.265,'kcal/mol'), T0=(1,'K'), comment="""Estimated from node Root_Ext-5R-R_7R!H->C_N-7C-inRing in family Peroxyl_Disproportionation.""")
kinetics: Arrhenius(A=(3.2e+12,'cm^3/(mol*s)'), n=0, Ea=(4.064,'kcal/mol'), T0=(1,'K'), comment="""Estimated from node Root_Ext-5R-R_7R!H->C_N-7C-inRing_Ext-5R-R in family Peroxyl_Disproportionation.""")
kinetics: Estimated from node Root_Ext-5R-R_7R!H->C_N-7C-inRing in family Peroxyl_Disproportionation.
kinetics: Estimated from node Root_Ext-5R-R_7R!H->C_N-7C-inRing_Ext-5R-R in family Peroxyl_Disproportionation.
Errors occurred during edge comparison ⚠️
ERROR conda.cli.main_run:execute(148): `conda run python scripts/checkModels.py liquid_oxidation-edge stable_regression_results/liquid_oxidation/chemkin/chem_edge_annotated.inp stable_regression_results/liquid_oxidation/chemkin/species_edge_dictionary.txt test/regression/liquid_oxidation/chemkin/chem_edge_annotated.inp test/regression/liquid_oxidation/chemkin/species_edge_dictionary.txt` failed. (See above for error)
liquid_oxidation Passed Observable Testing ✅
Regression test nitrogen:
Reference: Execution time (DD:HH:MM:SS): 00:00:00:56
Current: Execution time (DD:HH:MM:SS): 00:00:01:02
Reference: Memory used: 909.62 MB
Current: Memory used: 908.46 MB
nitrogen Failed Core Comparison ❌
Original model has 41 species.
Test model has 41 species. ✅
Original model has 360 reactions.
Test model has 359 reactions. ❌
The original model has 1 reactions that the tested model does not have. ❌
rxn: HNO(48) + HCO(13) <=> NO(38) + CH2O(18) origin: H_Abstraction
Errors occurred during core comparison ⚠️
ERROR conda.cli.main_run:execute(148): `conda run python scripts/checkModels.py nitrogen-core stable_regression_results/nitrogen/chemkin/chem_annotated.inp stable_regression_results/nitrogen/chemkin/species_dictionary.txt test/regression/nitrogen/chemkin/chem_annotated.inp test/regression/nitrogen/chemkin/species_dictionary.txt` failed. (See above for error)
nitrogen Failed Edge Comparison ❌
Original model has 133 species.
Test model has 133 species. ✅
Original model has 983 reactions.
Test model has 981 reactions. ❌
Non-identical thermo! ❌
original: O1[C]=N1
tested: O1[C]=N1
| Hf(300K) | S(300K) | Cp(300K) | Cp(400K) | Cp(500K) | Cp(600K) | Cp(800K) | Cp(1000K) | Cp(1500K) |
|---|---|---|---|---|---|---|---|---|
| 116.46 | 53.90 | 11.62 | 12.71 | 13.49 | 13.96 | 14.14 | 13.85 | 13.58 |
| 141.64 | 58.66 | 12.26 | 12.27 | 12.09 | 11.96 | 12.26 | 12.72 | 12.15 |
thermo: Thermo group additivity estimation: group(O2s-CdN3d) + group(N3d-OCd) + group(Cd-HN3dO) + ring(Cyclopropene) + radical(CdJ-NdO)
thermo: Thermo group additivity estimation: group(O2s-CdN3d) + group(N3d-OCd) + group(Cd-HN3dO) + ring(oxirene) + radical(CdJ-NdO)
The original model has 2 reactions that the tested model does not have. ❌
rxn: HNO(48) + HCO(13) <=> NO(38) + CH2O(18) origin: H_Abstraction
rxn: HON(T)(83) + HCO(13) <=> NO(38) + CH2O(18) origin: Disproportionation
Non-identical kinetics! ❌
original:
rxn: NCO(66) <=> O1[C]=N1(126) origin: Intra_R_Add_Endocyclic
tested:
rxn: NCO(66) <=> O1[C]=N1(126) origin: Intra_R_Add_Endocyclic
| k(1bar) | 300K | 400K | 500K | 600K | 800K | 1000K | 1500K | 2000K |
|---|---|---|---|---|---|---|---|---|
| k(T): | -49.54 | -33.65 | -24.16 | -17.85 | -10.01 | -5.35 | 0.80 | 3.82 |
| k(T): | -66.25 | -46.19 | -34.19 | -26.21 | -16.28 | -10.36 | -2.54 | 1.31 |
kinetics: Arrhenius(A=(6.95187e+18,'s^-1'), n=-1.628, Ea=(88.327,'kcal/mol'), T0=(1,'K'), comment="""Estimated from node Backbone0_N-2R!H-inRing_N-1R!H-inRing_Sp-2R!H-1R!H in family Intra_R_Add_Endocyclic.""")
kinetics: Arrhenius(A=(6.95187e+18,'s^-1'), n=-1.628, Ea=(111.271,'kcal/mol'), T0=(1,'K'), comment="""Estimated from node Backbone0_N-2R!H-inRing_N-1R!H-inRing_Sp-2R!H-1R!H in family Intra_R_Add_Endocyclic.""")
Identical kinetics comments:
kinetics: Estimated from node Backbone0_N-2R!H-inRing_N-1R!H-inRing_Sp-2R!H-1R!H in family Intra_R_Add_Endocyclic.
Errors occurred during edge comparison ⚠️
ERROR conda.cli.main_run:execute(148): `conda run python scripts/checkModels.py nitrogen-edge stable_regression_results/nitrogen/chemkin/chem_edge_annotated.inp stable_regression_results/nitrogen/chemkin/species_edge_dictionary.txt test/regression/nitrogen/chemkin/chem_edge_annotated.inp test/regression/nitrogen/chemkin/species_edge_dictionary.txt` failed. (See above for error)
nitrogen Passed Observable Testing ✅
Regression test oxidation:
Reference: Execution time (DD:HH:MM:SS): 00:00:01:34
Current: Execution time (DD:HH:MM:SS): 00:00:01:34
Reference: Memory used: 791.90 MB
Current: Memory used: 792.48 MB
oxidation Passed Core Comparison ✅
Original model has 59 species.
Test model has 59 species. ✅
Original model has 694 reactions.
Test model has 694 reactions. ✅
oxidation Passed Edge Comparison ✅
Original model has 230 species.
Test model has 230 species. ✅
Original model has 1524 reactions.
Test model has 1524 reactions. ✅
oxidation Passed Observable Testing ✅
Errors occurred during observable testing ⚠️
WARNING: Initial mole fractions do not sum to one; normalizing.
Regression test sulfur:
Reference: Execution time (DD:HH:MM:SS): 00:00:00:37
Current: Execution time (DD:HH:MM:SS): 00:00:00:38
Reference: Memory used: 907.87 MB
Current: Memory used: 905.30 MB
sulfur Passed Core Comparison ✅
Original model has 27 species.
Test model has 27 species. ✅
Original model has 74 reactions.
Test model has 74 reactions. ✅
sulfur Failed Edge Comparison ❌
Original model has 89 species.
Test model has 89 species. ✅
Original model has 227 reactions.
Test model has 227 reactions. ✅
The original model has 1 reactions that the tested model does not have. ❌
rxn: O(4) + SO2(15) (+N2) <=> SO3(16) (+N2) origin: primarySulfurLibrary
The tested model has 1 reactions that the original model does not have. ❌
rxn: O(4) + SO2(15) (+N2) <=> SO3(16) (+N2) origin: primarySulfurLibrary
Errors occurred during edge comparison ⚠️
ERROR conda.cli.main_run:execute(148): `conda run python scripts/checkModels.py sulfur-edge stable_regression_results/sulfur/chemkin/chem_edge_annotated.inp stable_regression_results/sulfur/chemkin/species_edge_dictionary.txt test/regression/sulfur/chemkin/chem_edge_annotated.inp test/regression/sulfur/chemkin/species_edge_dictionary.txt` failed. (See above for error)
sulfur Passed Observable Testing ✅
Regression test superminimal:
Reference: Execution time (DD:HH:MM:SS): 00:00:00:23
Current: Execution time (DD:HH:MM:SS): 00:00:00:24
Reference: Memory used: 967.54 MB
Current: Memory used: 967.45 MB
superminimal Passed Core Comparison ✅
Original model has 13 species.
Test model has 13 species. ✅
Original model has 21 reactions.
Test model has 21 reactions. ✅
superminimal Passed Edge Comparison ✅
Original model has 18 species.
Test model has 18 species. ✅
Original model has 28 reactions.
Test model has 28 reactions. ✅
Regression test RMS_constantVIdealGasReactor_superminimal:
Reference: Execution time (DD:HH:MM:SS): 00:00:02:56
Current: Execution time (DD:HH:MM:SS): 00:00:02:20
Reference: Memory used: 2615.25 MB
Current: Memory used: 2478.56 MB
RMS_constantVIdealGasReactor_superminimal Passed Core Comparison ✅
Original model has 13 species.
Test model has 13 species. ✅
Original model has 19 reactions.
Test model has 19 reactions. ✅
RMS_constantVIdealGasReactor_superminimal Passed Edge Comparison ✅
Original model has 13 species.
Test model has 13 species. ✅
Original model has 19 reactions.
Test model has 19 reactions. ✅
RMS_constantVIdealGasReactor_superminimal Passed Observable Testing ✅
Regression test RMS_CSTR_liquid_oxidation:
Reference: Execution time (DD:HH:MM:SS): 00:00:14:29
Current: Execution time (DD:HH:MM:SS): 00:00:16:13
Reference: Memory used: 2860.23 MB
Current: Memory used: 2677.89 MB
RMS_CSTR_liquid_oxidation Failed Core Comparison ❌
Original model has 35 species.
Test model has 35 species. ✅
Original model has 189 reactions.
Test model has 189 reactions. ✅
The original model has 1 species that the tested model does not have. ❌
spc: CCH2
The tested model has 1 species that the original model does not have. ❌
spc: [CH2]CCC(C)O(93)
The original model has 1 reactions that the tested model does not have. ❌
rxn: C[CH2](6) + [CH2]CC(5) <=> pentane(2) origin: R_Recombination
The tested model has 1 reactions that the original model does not have. ❌
rxn: CCCC(C)[O](44) <=> [CH2]CCC(C)O(93) origin: intra_H_migration
Errors occurred during core comparison ⚠️
ERROR conda.cli.main_run:execute(148): `conda run python scripts/checkModels.py RMS_CSTR_liquid_oxidation-core stable_regression_results/RMS_CSTR_liquid_oxidation/chemkin/chem_annotated.inp stable_regression_results/RMS_CSTR_liquid_oxidation/chemkin/species_dictionary.txt test/regression/RMS_CSTR_liquid_oxidation/chemkin/chem_annotated.inp test/regression/RMS_CSTR_liquid_oxidation/chemkin/species_dictionary.txt` failed. (See above for error)
RMS_CSTR_liquid_oxidation Passed Edge Comparison ✅
Original model has 107 species.
Test model has 107 species. ✅
Original model has 545 reactions.
Test model has 545 reactions. ✅
RMS_CSTR_liquid_oxidation Passed Observable Testing ✅
Regression test fragment:
Reference: Execution time (DD:HH:MM:SS): 00:00:00:31
Current: Execution time (DD:HH:MM:SS): 00:00:00:33
Reference: Memory used: 753.67 MB
Current: Memory used: 753.36 MB
fragment Passed Core Comparison ✅
Original model has 10 species.
Test model has 10 species. ✅
Original model has 2 reactions.
Test model has 2 reactions. ✅
fragment Passed Edge Comparison ✅
Original model has 33 species.
Test model has 33 species. ✅
Original model has 47 reactions.
Test model has 47 reactions. ✅
fragment Passed Observable Testing ✅
Errors occurred during observable testing ⚠️
WARNING: Initial mole fractions do not sum to one; normalizing.
Regression test RMS_constantVIdealGasReactor_fragment:
Reference: Execution time (DD:HH:MM:SS): 00:00:03:26
Current: Execution time (DD:HH:MM:SS): 00:00:02:48
Reference: Memory used: 2737.40 MB
Current: Memory used: 2638.75 MB
RMS_constantVIdealGasReactor_fragment Passed Core Comparison ✅
Original model has 10 species.
Test model has 10 species. ✅
Original model has 2 reactions.
Test model has 2 reactions. ✅
RMS_constantVIdealGasReactor_fragment Passed Edge Comparison ✅
Original model has 27 species.
Test model has 27 species. ✅
Original model has 24 reactions.
Test model has 24 reactions. ✅
RMS_constantVIdealGasReactor_fragment Passed Observable Testing ✅
Errors occurred during observable testing ⚠️
WARNING: Initial mole fractions do not sum to one; normalizing.
Regression test minimal_surface:
Reference: Execution time (DD:HH:MM:SS): 00:00:00:32
Current: Execution time (DD:HH:MM:SS): 00:00:00:32
Reference: Memory used: 910.52 MB
Current: Memory used: 912.98 MB
minimal_surface Passed Core Comparison ✅
Original model has 11 species.
Test model has 11 species. ✅
Original model has 3 reactions.
Test model has 3 reactions. ✅
minimal_surface Passed Edge Comparison ✅
Original model has 38 species.
Test model has 38 species. ✅
Original model has 38 reactions.
Test model has 38 reactions. ✅
minimal_surface Passed Observable Testing ✅
beep boop this comment was written by a bot 🤖
Fixes the two diagnostics described in #3003. Both came out of splicing QM transition states into an otherwise RMG-estimated pressure-dependent network, where a mismatch silently inflated k by ~1e14 and then surfaced, thousands of lines later, as
numpy.linalg.LinAlgError: Singular matrix.calculate_tst_rate_coefficientnow warns when the transition state and a reactant don't both enumerate external modes. The TST expression assumes translation and rotation cancel betweenQ‡andQ_reactant; when one side is described vibration-only and the other isn't, they don't.It warns rather than raises — a mixed description can be deliberate, and raising would break working code. The check is deliberately a yes/no on whether external modes are present at all, rather than 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 flags differences that are perfectly correct. I had it that way first and it fired on
H + C2H4, which is what the monatomic test now pins. It warns once per reaction, since this is called once per temperature inTlist.The MSC solve now says where and why it went singular.
a_matdegenerates when the rate coefficients in a grain overwhelm the collision frequency — for two isomers it reduces to[[-k21, k12], [k21, -k12]], exactly rank-1. The error now carries the grain, J index, temperature, and the two scales with their ratio, so it's clear the input was unphysical rather than the algorithm unlucky. The intent was already there a few lines up:bis checked for NaN and logs a warning naming the grain and isomer.a_mathad nothing, and it's the argument that fails.It deliberately doesn't regularize or perturb
a_mat— it's singular because the physics handed to it was wrong, and smoothing that over turns a loud failure into a quiet wrong number.Neither change moves any number a working network computes: the first only logs, the second only runs on a path that already raised.
Testing
test/rmgpy/pdep+test/rmgpy/reactionTest.py: 86 passed, 2 skipped.test/rmgpy/statmech+test/rmgpy/kinetics: 372 passed.test/arkane: 231 passed, 7 skipped — which is the real check that the warning doesn't fire on legitimate networks.Seven new tests, and I checked each one fails without its fix:
try/exceptremoved, both singular-matrix tests go red while the well-conditioned control stays green.One thing worth mentioning for review, because it argues against a natural simplification: the singular matrix is not a general tripwire for this. I tested putting translation and rotation on the wells and the transition state uniformly — the partition functions cancel again, the solve succeeds, and the fitted rate coefficient is unchanged, on a network that's still wrong. Only the asymmetric case blows up. So the warning isn't redundant with the crash; it covers a case the crash can't see.