RRKM

Calculate unimolecular reaction rates using Rice-Rammsperger-Kassel-Marcus (RRKM) theory in Python.

The Graph Kinetic Monte Carlo simulator depends on reaction rates. The RRKM tool calculates these rates from the reactant density of states and transition-state sum of states:^[1]

\[k(E) = \frac{\sigma N^\dagger(E-E_0)}{h\rho(E)}\]

The RRKM module is supplied with the GKMC repository in the tools/rrkm/ subdirectory. There is a Python API (rrkm.py) and command-line interface (CLI; rrkm_cli.py).

Contents

• Contents
• Command-line interface
• How it works
• Public API

Command-line interface

Installation

See Additional tools > Installation for installation instructions.

Usage

Instructions for the command-line interface can be retrieved with:

uv run tools/rrkm_cli.py --help

The CLI reads quantum-chemistry output files from --input-directory and writes JSON rate files to --output-directory. Both default to the current directory. The output directory is created when it does not exist.

For example:

uv run tools/rrkm_cli.py D 1 --destination 2 --type scramble \
    --input-directory path/to/calculations \
    --output-directory path/to/rates

This writes D1to2.json to path/to/rates. For dissociation:

uv run tools/rrkm_cli.py D 1 --type dissociate \
    --input-directory path/to/calculations \
    --output-directory path/to/rates

This writes D1diss.json to the output directory. Use --reverse to also write the reverse scrambling rate and --suffix to select a suffixed dissociation transition-state file.

How it works

The RRKM API provides two functions:

• rrkm.scramble() calculates rates for isotope scrambling between two aliphatic sites.
• rrkm.dissociate() calculates rates for isotope dissociation from an aliphatic site.

Both functions return a dictionary that follows the ratefile.schema.json.

Input files

The functions read quantum-chemistry output files using cclib. Each file must contain vibrational frequencies, electronic energies, and zero-point vibrational energy. The calculation must have completed successfully.

The directory argument is the parent directory of the following layout:

For a scrambling reaction:

directory/
├── Stable/
│   ├── {isotope}{source}.log
│   └── {isotope}{destination}.log
└── TS/
    └── {isotope}{source}to{destination}.log

For a dissociation reaction:

directory/
├── Stable/
│   └── {isotope}{source}.log
├── Dissociation/
│   └── {isotope}{source}diss{suffix}.log
└── {isotope}-Dissociated.log

For scrambling, the threshold is $E_0=E_\mathrm{TS}-E_\mathrm{reactant}$ and the reaction energy is $\Delta=E_\mathrm{product}-E_\mathrm{reactant}$. For dissociation, the threshold is calculated as $E_0=E_\mathrm{product}+E_\mathrm{atom}-E_\mathrm{reactant}$, where atom_energy_eV defaults to $-13.6\ \mathrm{eV}$. The dissociation reaction energy is reported as $\Delta=0$.

The reactant must have no imaginary frequencies, and the transition state must have exactly one. The imaginary transition-state mode is excluded automatically from the state-count calculation. Frequencies are scaled by scaling_factor before they are discretized.

Public API

scramble(isotope, source, destination, ...)

Returns:

A dict containing:

The energy_eV and rate_s_inv lists have equal lengths.

Raises:

• ValueError if isotope is not "H" or "D".
• ValueError if a log file is missing required data.
• ValueError if a quantum-chemistry calculation was unsuccessful.
• ValueError if energies or frequencies are non-finite.
• ValueError if the starting structure has imaginary frequencies.
• ValueError if the transition state does not have exactly one imaginary frequency.
• ValueError if the threshold energy is negative or non-finite.
• ValueError if scaling_factor, energy_grain_eV, or energy_max_eV is invalid.
• ValueError if a frequency rounds to zero energy quanta.
• OverflowError if the state counts overflow.
• FloatingPointError if calculated rates are non-finite or negative.

dissociate(isotope, source, ...)

Returns the same rate-file dictionary as scramble(), with "reactionType": "dissociate" and "Delta": 0.0. Raises the same calculation errors as scramble().