Skip to content

Mismatched conformer degrees of freedom pass silently through TST, then surface as a bare "Singular matrix" in the MSC solve #3003

Description

@alongd

Ran into both of these while splicing QM transition states into an otherwise RMG-estimated pressure-dependent network, and they compound in a way that made the diagnosis much harder than it needed to be.

If the transition state conformer carries translation and rotation but the wells are vibration-only (or the other way round), Reaction.calculate_tst_rate_coefficient (rmgpy/reaction.py:1441) just multiplies the partition functions together — nothing anywhere in it looks at what modes either side actually has. In my case Q‡/Q_reactant was off by ~1e14 and k_TST(3200 K) came out at 1.8e26 s^-1 for a unimolecular step.

That rescales k(E), and the modified strong collision solve then dies at np.linalg.solve(a_mat, b) in rmgpy/pdep/msc.pyx with nothing but numpy.linalg.LinAlgError: Singular matrix — no grain, no J, no temperature, no pressure. I ended up bisecting by hand down to grain 346 / J=0 / 3200 K / 0.1 bar.

The bit that struck me is that the intent is already there a few lines up: b gets checked with np.isnan and logs a warning naming the grain and the isomer. a_mat has no equivalent check, and it's the argument that actually fails.

A warning in calculate_tst_rate_coefficient when the reactant and TS conformers carry different mode classes, or just the grain/J/T/P coordinates on the singular-matrix failure, would have turned a two-hour hunt into an immediate error message. Either one alone would have done it.

One thing worth flagging for anyone fixing this: the loud failure only happens when the mismatch is asymmetric. I checked afterwards, and if you put translation and rotation on the wells and the TS alike, the partition functions cancel again — the solve succeeds and the rate coefficient comes back unchanged, on a network that's still wrong. So the singular matrix isn't a general tripwire for this, which is an argument for the warning rather than relying on the crash.

Happy to open a PR for both.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions