MicroKatc

Automated microkinetic analysis of interconnected catalytic cycles, from DFT outputs to mechanistic insight.




_{Graphical abstract © 2025 the authors, ACS Catal. 15, 4739, CC BY-NC-ND 4.0}

MicroKatc turns a set of Gaussian calculations and a list of elementary steps into a complete kinetic picture of a catalytic system. It computes thermally corrected Gibbs energies, builds and runs COPASI microkinetic models, and then answers the questions a mechanistic study needs: which steps control the rate, how the apparent activation energy responds to reaction conditions, and where the catalyst actually spends its time.

The approach is published in:

O. Abdullayev, D. Garay-Ruiz, B. Bori-Bru, C. Bo. "Microkinetic Assessment of Ligand-Exchanging Catalytic Cycles." ACS Catal. 2025, 15 (6), 4739–4745. doi:10.1021/acscatal.5c00348

Highlights

• End-to-end pipeline: raw DFT output files in, publication-ready figures out, driven by one command.
• One study file per system: steps written as chemists write them (A + 2 B <=> C via TS), energies from Gaussian output files or typed in kcal mol^-1, and every problem in the file reported with its location before anything runs.
• Apparent activation energy (E_a) from Arrhenius fits of every step's flux and every species' rate of change, swept over temperature and over the concentration of any chosen reactant.
• Degree of rate control (DRC) for every elementary step, by central finite differences parallelised across CPU cores, to identify the rate-determining and inhibiting steps.
• Multi-cycle catalyst tracking: the catalyst concentration in each competing cycle (for example 0-ligand vs. 1-ligand cycles) as conditions change.
• Conversion-time analysis: time to reach a set product yield as a function of reactant concentration.
• Result caching: simulations and fitted parameters are saved as CSV files and reused, so long parameter sweeps never repeat finished work.

How it works

Module	Role
microkatc.py	The command: run and check a study file
study.py, steps.py, study_yaml.py	Read and check the study file
thermochemistry.py, calculating_G_for_microkinetics.py	G of each species (thermochange or typed) and the forward/reverse barrier of each step
microkinetics_simulation.py	COPASI simulations with convergence checks and caching; catalyst and conversion analyses
apparent_activation_energy.py	Apparent Ea fits and DRC calculation
plotting_functions.py	All figures

Case study: Rh-catalysed hydroformylation

The example data in this repository model the hydroformylation of ethene by a homogeneous rhodium catalyst. Two cycles compete: one without (0L) and one with (1L) a coordinated PMe₃ ligand. The analysis varies the initial PMe₃ concentration (studied_species: PMe3 in the study file).

Key findings

• Adding PMe₃ shifts the catalyst from the 0L cycle into the 1L cycle. The 1L cycle has lower activation energies, becomes the main source of product, and reaches 99 % conversion much faster.
• The rate-determining step moves with the conditions: I8_0L ⇌ I9_0L controls the rate at low PMe₃, and I3_1L ⇌ I4_1L takes over at high PMe₃. The negative DRC of I3_0L ⇌ I4_0L shows that this step inhibits the 0L cycle.
• More PMe₃ also releases more CO, which poisons the catalyst at the I6 ⇌ I7 steps. This is why the E_a of the rate-determining steps rises and then plateaus.

All figures below are built by readme_figures.py from the results of the example study. Blue is the 0L cycle, pink the 1L cycle and green the product. Every run also writes the full versions, for every step and every concentration, to results/microkinetics_simulations_images/.

1. Arrhenius check

The apparent E_a comes from the slope of ln(r) against 1/T over 325–375 K, which is valid at low catalyst concentration. In each regime, product formation follows the same straight line as its rate-determining step, so the two share one E_a. Every run draws this plot for every step and keeps those with R² > 0.9.

2. Apparent activation energy

The E_a of both rate-determining steps rises with PMe₃ and then plateaus as the catalyst saturates with PMe₃. The rise comes from CO release: more PMe₃ frees more CO, which poisons the catalyst at the I6 ⇌ I7 steps. E_a is lower in the 1L cycle than in the 0L cycle, so the 1L cycle has faster kinetics and takes over product formation. The product's E_a therefore falls from 23.5 to about 21.5 kcal mol^-1, then settles at the 1L value of 22.1 kcal mol^-1.

3. Degree of rate control

I8_0L ⇌ I9_0L controls the rate in the 0L cycle, but it loses importance as the 1L cycle becomes active. The DRC of I3_1L ⇌ I4_1L grows with PMe₃, making it the controlling step at high concentration. The DRC of I3_0L ⇌ I4_0L turns negative in the transition region, which shows that this step inhibits the 0L cycle.

4. Catalyst distribution

The share of the catalyst in each cycle at t = 1 h, with 5 × 10^-4 M catalyst. Beyond about 5 × 10^-4 M PMe₃, most of the catalyst is in the 1L cycle.

5. Concentration profiles

The resting states I1 and I7 of both cycles over time, at low, intermediate and high PMe₃. As PMe₃ increases, the catalyst moves from the 0L species into their 1L counterparts.

6. Time to 99 % conversion

The time for the product to reach 99 % of the maximum allowed by the limiting reactant (the threshold is adjustable). It drops from 5.7 h with the 0L cycle alone to 1.7 h once the 1L cycle takes over.

Putting it together

Why does the E_a of the 1L rate-determining step increase (figure 2) even though the 1L cycle becomes more active? As PMe₃ increases, I3_0L ⇌ I4_0L inhibits the 0L cycle and pushes the catalyst into the 1L cycle (figure 3). At the same time, the poisoning intermediates of the 0L cycle (I1_0L, I7_0L) decrease, and those of the 1L cycle (I1_1L, I7_1L) build up (above, at the low catalyst concentration of the E_a analysis). Product formation in the 1L cycle therefore gets harder, which raises its E_a.

Reproducing the paper

Running examples/hydroformylation/study.yaml takes a few minutes and reproduces the published results:

Published result	Paper	This code
Gibbs barriers of all 22 steps at 5 temperatures (SI Tables S1–S5)	0.1 kcal mol-1 precision	identical
Time to 99 % conversion, 0L only → 1L (Figure 4)	5.66 h → 1.69 h	5.67 h → 1.69 h
Apparent Ea of product formation (Figure 6)	23.5 → 21.4 → 22.1 kcal mol-1	23.5 → 21.5 → 22.1 kcal mol-1
Minimum DRC of I3_0L ⇌ I4_0L (Figure 5)	−0.22 at 3 × 10-4 M	−0.22 at 3 × 10-4 M

tests/test_paper_barriers.py checks the barriers against SI Table S3 on every pull request, and a nightly run checks the full results against the paper.

Quick start

git clone https://gitlab.com/dgarayr/thermochange.git   # thermochemical corrections, by D. Garay-Ruiz
git clone https://github.com/0rkhann/MicroKatc.git
cd MicroKatc
python3.10 -m venv .venv && source .venv/bin/activate
pip install -r requirements.txt                          # includes COPASI's Python bindings and copasi_helper
export thermochange=/path/to/thermochange
python microkatc.py run examples/hydroformylation/study.yaml   # the paper's example

To study your own reaction, describe it in one YAML study file: the steps (A + B <=> C via TS), where the energies come from (Gaussian output files or typed Gibbs energies), the conditions and the analyses. Then run python microkatc.py run my_study.yaml. The usage guide covers every key, both examples, the outputs and troubleshooting.

Modelling assumptions

1. Electronic structure: species were computed with DFT at the ωB97X-D/6-311G(d,p) level.
2. Standard state: the pressure is set from the temperature so that the ideal-gas concentration is 1 M, approximating a liquid medium: $C = \frac{P}{RT}$.
3. Barrierless steps: steps without a located transition state are treated as diffusion-controlled, with a barrier of 4 kcal/mol (Besora et al., 2018).
4. Network topology: only the listed steps between neighbouring intermediates are included. Cross-cycle or non-neighbour transformations are not modelled.
5. Apparent E_a: at low catalyst concentration, the rate constant $k_i$ in the linearised Arrhenius equation can be replaced by the step flux $r_i$ or the species rate $v_i$:

$$\ln r_i = \ln A - \frac{E_a}{RT}, \qquad \ln v_i = \ln A - \frac{E_a}{RT}$$

Citation

If you use MicroKatc in your work, please cite:

@article{Abdullayev2025,
  author  = {Abdullayev, Orkhan and Garay-Ruiz, Diego and Bori-Bru, Berta and Bo, Carles},
  title   = {Microkinetic Assessment of Ligand-Exchanging Catalytic Cycles},
  journal = {ACS Catalysis},
  year    = {2025},
  volume  = {15},
  number  = {6},
  pages   = {4739--4745},
  doi     = {10.1021/acscatal.5c00348}
}

Citation metadata is also in CITATION.cff; GitHub shows it under "Cite this repository".

Acknowledgements

MicroKatc was developed during a summer research stay at ICIQ. I sincerely thank my supervisor, Dr. Diego Garay-Ruiz, and Principal Investigator, Prof. Carles Bo, for their guidance and support. I also thank my fellow summer research colleagues, who became dear friends and made the experience truly memorable.

License

MIT