Example Simulation

In this tutorial, we will model the reaction A + 1/2 B₂ AB using simplified placeholder values for demonstration purposes.

Tip

For anyone starting to work with MKMCXX or with microkinetic simulations in general, this tutorial is an excellent starting point which covers the basic features and settings of MKMCXX, including plotting and data interpretation.

Input file

All v3.x input files should start with the header tag # MKMCXXv3 to indicate that the input file is formatted for version 3.x. Next, input files requires five blocks to define the simulation:

&compounds: Which species are in the system
&reactions: Which elementary reaction steps link these species
&settings: Under which conditions is the microkinetic simulation executed
&runs: Which simulation scenarios should be executed
&output: How simulation output should be written to permanent storage.

Below, we will address each of these blocks.

Compounds

We assume the reaction proceeds as follows. Species A adsorbs molecularly on the catalytic surface, while B₂ adsorbs dissociatively. These adsorption steps produce the surface species A and B, respectively. The adsorbed species A and B can then combine to form AB, releasing an empty site. Finally, AB desorbs from the surface to regenerate the active site.

For this model, we therefore require three gas-phase species (A, B2, AB) and four surface species (A*, B*, AB*, and *, the latter representing an empty site).

Tip

We use the suffix * to indicate an adsorbed species. Later, when we will be plotting the results, we can use this suffix to identify adsorbed species using their labels.

For an optimal reaction, we assume that the gas phase pressures of A and B₂ are according to their stoichiometric ratio, i.e. 1 atm for A and 0.5 atm for B₂.

The &compounds block thus looks as follows

&compounds
# name ; c = concentration ; FLAGS = flags (comma-separated)
A;      c = 1.0; FLAGS = GAS, REACTANT
B2;     c = 0.5; FLAGS = GAS, REACTANT
AB;     c = 0.0; FLAGS = GAS, PRODUCT, KEY

# Adsorbed compounds
A*;         c = 0.0; FLAGS = SURFACE
B*;         c = 0.0; FLAGS = SURFACE
AB*;        c = 0.0; FLAGS = SURFACE
*;          c = 1.0; FLAGS = SURFACE
&endblock

In this block, we indicate for each compound:

Its label
Its starting concentration, indicating via a c = value keyword-value pair.
A number of FLAGS, identifying the compounds as either REACTANT or PRODUCT and whether the compound is a GAS-phase or a SURFACE-phase species. Finally, when performing sensitivity analysis, we need to identify a KEY-compound. Here, we have chosen AB to be the KEY-compound.

Note

Always define at least one KEY-compound, else any kind of sensitivity analysis will fail.
Although typically not critical if omitted, it is recommended to always indicate whether species are GAS or REACTANT type and whether they are considered to be REACTANT or PRODUCT.
Do not forget to close the block using an &endblock tag!

Reactions

To model the reactions, we will introduce four elementary reaction steps, corresponding to adsorption/desorption of the reactants and products and a surface reaction where A and B are combined. As indicated above, A and AB will adsorb/desorb molecularly and B2 will adsorb/desorb dissociatively/recombinatorily. This would yield the following reaction block.

&reactions
{A}     +   {*}         =>  {A*};                   ; model = HertzKnudsenDefault; Asite = 1e-20; m = 12; theta = 1.0; sigma = 1; S = 1; Edes = 120e3; FLAGS = DRC
{B2}    +   2{*}        =>  2{B*};                  ; model = HertzKnudsenDefault; Asite = 1e-20; m = 28; theta = 1.0; sigma = 1; S = 1; Edes = 80e3;  FLAGS = DRC
{AB}    +   {*}         =>  {AB*};                  ; model = HertzKnudsenDefault; Asite = 1e-20; m = 26; theta = 1.0; sigma = 1; S = 1; Edes = 50e3;  FLAGS = DRC
{A*}    +   {B*}        =>  {AB*} + {*};            ; model = ArrheniusDefault; Vf = 1e13; Vb = 1e13; Eaf = 120e3; Eab = 80e3; FLAGS = DRC
&endblock

In MKMCXX, each reaction is written as a reaction statement, which is the part before the first semicolon (;). The general rules are:

Species labels are enclosed in curly brackets {...}.
- Example: {A}, {B2}, {A*}, {*}
Stoichiometric factors are written in front of the label.
- Example: 2{B*} means two adsorbed B atoms.
The reaction arrow => separates reactants (left) from products (right). Note that all elementary reaction steps are assumed to be microscopically reversible, so the reactions will always be evaluated in both forward and backward direction. Thus, you only have to define them once.
Everything after the semicolon contains kinetic parameters and does not belong to the reaction statement itself.

In this example, we consider two types of elementary reaction models:

HertzKnudsenDefault:
Used for adsorption and desorption steps. This model applies the Hertz–Knudsen equation to calculate rates based on molecular flux, sticking probability, site area, and desorption energy. It is typically chosen for gas–surface interactions.
ArrheniusDefault:
Used for surface reactions between adsorbed species. This model follows the Arrhenius equation, defined by (temperature-independent) pre-exponential factors and activation energies for forward and backward reactions. It is typically chosen for bond-forming and bond-breaking steps on the surface.

Warning

Although elementary reaction steps are microscopically reversible, elementary reaction steps for which the HertzKnudsenDefault model is being used need to be written such that the gas-phase components are on the left-hand side of the reaction.

Hertz-Kndusen

For Hertz-Knudsen reaction, by default, the following forms for the reaction rate constants are assumed.

\[ k_{f} = \frac{pA_{\text{site}}}{\sqrt{2 \pi m k_{b} T}} S \]

and

\[ k_{b} = \frac{k_{b}T^{3}}{h^{3}} \frac{A_{\text{site}}(2 \pi m k_{b})}{\sigma \theta_{\text{rot}}} \exp \left( -\frac{\Delta E_{\text{des}}}{k_{b}T} \right) \]

Given these equations, the user needs to define, using key=value pairs, the following parameters:

Asite: Active site surface area in m^-2.
m: Mass of the species in amu. (atomic units)
theta: Rotational temperature in K
sigma: Symmetry factor (dimensionless)
S: Sticking coefficient (dimensionless)
Edes: Desorption energy in J/mol

Note

The desorption energy is the energy that needs to be invested to remove the species from the surface. As such, this should in principle be a positive value.

Arrhenius

For Arrhenius equations, we use the following form for the reaction rate constants. This functional form assumes that the pre-exponential factor is temperature-independent. This is an approximation.

\[ k_{f} = \nu_{f} \exp \left( -\frac{\Delta E_{\text{a,f}}}{RT} \right) \]

and

\[ k_{b} = \nu_{b} \exp \left( -\frac{\Delta E_{\text{b,f}}}{RT} \right) \]

Given these equations, the user needs to define, using key=value pairs, the following parameters:

Vf: Pre-exponential factor in the forward direction in s^-1.
Vb: Pre-exponential factor in the backward direction in s^-1
Eaf: Forward barrier in J/mol.
Eab: Backward barrier in J/mol.

Tip

A typical ballpark value for the forward and backward pre-exponential factors can be derived from the equation \(k_{b}T/h\) which is roughly 1e13 s^-1.

Settings

In the &settings block, we need to indicate simulation settings which will be applied to every simulation scenario. Examples include sensitivity analyses. In this simulation, we will explore reaction orders, apparent activation energy and degree of rate control. This yields the following &settings block.

&settings
ORDERS = 1
DRC = 1
EACT = 1
&endblock

Runs

Under the &runs block, we need to define which simulation scenarios are being explored. For a default "thermal" simulation. Each scenario requires 4 parameters, being:

T: The simulation temperature in K.
t: The simulation time in s.
atol: Absolute tolerance for the ODE solver.
rtol: Relative tolerance for the ODE solver.

The simulation is performed by time-integrating the system of ordinary differential equations that represent the elementary reaction steps. Initial conditions are taken from the c=value definitions in the compounds block. By default, the simulation begins with an empty catalytic surface as the initial state. Gas-phase species are held constant, corresponding to a finite-conversion simulation. Time integration is carried out using finite time steps. The solver adjusts the step size dynamically, guided by user-defined relative and absolute tolerances. These tolerances ensure that increasing the step size, improving efficiency, does not introduce errors larger than the specified thresholds.

We use the following settings where we explore the system between T=350K and T=1200K using steps of 25K. We use an absolute and relative tolerance of 1e-8 throughout and t = 1e6 seconds of simulation time.

&runs
T = 350;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 375;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 400;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 425;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 450;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 475;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 500;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 525;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 550;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 575;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 600;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 625;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 650;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 675;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 700;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 725;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 750;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 775;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 800;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 825;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 850;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 875;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 900;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 925;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 950;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 975;    t = 1e6; atol = 1e-8; rtol = 1e-8
T = 1000;   t = 1e6; atol = 1e-8; rtol = 1e-8
T = 1025;   t = 1e6; atol = 1e-8; rtol = 1e-8
T = 1050;   t = 1e6; atol = 1e-8; rtol = 1e-8
T = 1075;   t = 1e6; atol = 1e-8; rtol = 1e-8
T = 1100;   t = 1e6; atol = 1e-8; rtol = 1e-8
T = 1125;   t = 1e6; atol = 1e-8; rtol = 1e-8
T = 1150;   t = 1e6; atol = 1e-8; rtol = 1e-8
T = 1175;   t = 1e6; atol = 1e-8; rtol = 1e-8
T = 1200;   t = 1e6; atol = 1e-8; rtol = 1e-8
&endblock

Tip

For the simulation time t, use a value significantly high by which you can reasonably assume that the system has reached its steady-state condition. Of course, always check it by looking into the time-profiles.
Good default values for the absolute and relative tolerances are 1e-8. Note that using values lower than ~1e-19 is not useful as this is lower than the machine precision for a 80-bit floating point number.

Output

The results of the simulation are written to a permanent medium (e.g. hard drive). In the &output-block, we can indicate which storage format should be used. If multiple storage formats are being specified, the simulation results are stored in all of these formats. In the example below, we write the simulation results in all three supported formats.

&output
FORMAT = EXCEL, TABDELIM, HDF5
&endblock

The supported formats are as follows:

HDF5: A High-performance data management and storage suite that efficiently stores heterogeneous data in a single file.
EXCEL: MS-XLSX Office Open XML Spreadsheet Format. Can be opened in Microsoft Excel or other spreadsheet software supporting this format. The format is also supported by the Python Data Analysis Library (pandas).
TABDELIM: A simple flat-text style format using tabs as the delimiting character.

Running the simulation

By placing all 5 blocks inside a single .mkx file, we can start the simulation by executing the following one-liner from a terminal application.

<path to mkmcxx executable> -i <path to input file>

Note

If you have made changes to the input file and wish to rerun the simulation, overwriting previous results, you need to add the -f instruction on the command line.

Typical output of the program is similar to the one shown below, here provided as an example.

PS D:\PROGRAMMING\MKMCXX\AB> .\mkmcxx.exe -i .\input.mkx
=================================================================================================
    _____ ______   ___  __    _____ ______   ________     ___    ___ ___    ___ ________
   |\   _ \  _   \|\  \|\  \ |\   _ \  _   \|\   ____\   |\  \  /  /|\  \  /  /|\_____  \
   \ \  \\\__\ \  \ \  \/  /|\ \  \\\__\ \  \ \  \___|   \ \  \/  / | \  \/  / ||____|\ /_
    \ \  \\|__| \  \ \   ___  \ \  \\|__| \  \ \  \       \ \    / / \ \    / /      \|\  \
     \ \  \    \ \  \ \  \\ \  \ \  \    \ \  \ \  \____   /     \/   /     \/      __\_\  \
      \ \__\    \ \__\ \__\\ \__\ \__\    \ \__\ \_______\/  /\   \  /  /\   \     |\_______\
       \|__|     \|__|\|__| \|__|\|__|     \|__|\|_______/__/ /\ __\/__/ /\ __\    \|_______|
                                                         |__|/ \|__||__|/ \|__|

=================================================================================================

MKMCXX Version: 3.0.0
Compiled with: GNU 15.1.0
Build type: Release
Compiler flags:  -O3 -DNDEBUG -flto=auto
Float type: long double
Float size: 16 bytes
Float precision: 64 bits
Build time: 2025-08-07 18:24:04
Author: Ivo Filot <<i.a.w.filot@tue.nl>>

=================================================================================================
CPU Vendor: GenuineIntel
CPU Model: 12th Gen Intel(R) Core(TM) i7-12700
Logical Threads: 20
Total RAM: 127.747 GB
OS: Windows
=================================================================================================
Reading input from: .\input.mkx
Read 7 compounds from input file.
Read 4 reactions from input file.
Done reading input file.

[Orchestrator] Launching 1155 jobs using 20 threads.
+ Running: 0   Γ£ô Completed: 1155   Γùï Queued: 0
[Orchestrator] All jobs completed successfully.


=============================================================================
                           COMPUTATION TIME REPORT
=============================================================================
Total jobs executed:              1155
Commulative computation time:     0.1594 s
Average job time:                 0.0001 s
Min/Max job time:                 0.0001 / 0.0018 s

Type       Jobs   Total[s]     Avg[s]     Min[s]     Max[s]
-----------------------------------------------------------------------------
CORE         35     0.0227     0.0006     0.0003     0.0018  Threads: 1
ORDER       420     0.0639     0.0002     0.0001     0.0004  Threads: 1
DRC         560     0.0479     0.0001     0.0001     0.0002  Threads: 1
EACT        140     0.0249     0.0002     0.0001     0.0004  Threads: 1
-----------------------------------------------------------------------------

Writing results in HDF5 format.
Writing results in TABDELIM format.
Writing results in EXCEL format.

=================================================================================================
All done!
Total simulation (wallclock) execution time: 0.236264 seconds
=================================================================================================

The program output contains:

Program and build (compilation) information
- Version number, compiler and flags, build type, floating-point precision, and build time.
System details
- CPU vendor and model, number of threads, total RAM, and operating system.
Input summary
- Name of the input file, number of compounds and reactions read.
Execution log
- Notes on overwritten directories, job orchestration, number of jobs launched, queued, and completed.
Computation time report
- Total number of jobs executed
- Cumulative, average, minimum, and maximum job times
- Breakdown by job type (e.g. CORE, ORDER, DRC, EACT), including per-thread statistics
Output files
- Confirmation of results written in HDF5, tab-delimited, and Excel formats.
Completion summary
- Total wall-clock execution time and confirmation that all jobs finished successfully.

Results, plotting and interpretation

The results of the simulation are written to a separate results folder located in the current working directory. To visualize and interpret the results, it is recommended to make use of small Python scripts. Here, we provide examples making use of the h5py library which can efficiently read out a hdf5 output file.

Concentrations

To visualize the surface concentration as function of temperature, we use the following small Python script

import h5py
import matplotlib.pyplot as plt

with h5py.File('results/results.h5', 'r') as f:
    if 'aggregated/kinetics/concentrations' in f:
        plt.figure(dpi=144)
        dset = f['aggregated']['kinetics']['concentrations']
        data = dset[:]
        for i in range(data.shape[1]-1):
            label = dset.attrs['column_labels'][i+1]
            if '*' in label:
                plt.plot(data[:, 0], data[:, i + 1], '-o', label=label)
        plt.grid(linestyle='--')
        plt.legend()
        plt.xlabel('Temperature [K]')
        plt.ylabel('Surface concentration [-]')

This script opens the .h5 file and looks for the data entry corresponding the concentrations as function of temprature (aggregated/kinetics/concentrations). If found, a graph is constructed where we loop over the compounds and visualize only those compounds which have a * in its label, i.e., only for the surface compounds. The reason we do not visualize the gas-phase compounds is because they are kept static throughout the simulation and are thus not relevant to investigate. The graph produced by this script is found below.

Concentrations as function of temperature

From this graph, we can readily see that at low temperature, the surface is predominantly covered by A*, in line with its higher desorption energy of 120 kJ/mol versus 80 kJ/mol for B*. With increasing temperature, we observe that the surface concentration of A* decreases in favor of *. This is a common pattern. With increasing temperature, desorption and reaction become more facile and the result is that AB* can be rapidly formed and desorbed from the surface by which the surface becomes gets increasingly more free sites.

Derivatives

One core observable is how temperature then affects the turn-over-frequecy for this catalyst. To investigate this, we need to look into the derivatives, which is another entry in the results.h5 file. To visualize it, we make use of the following Python script

import h5py
import matplotlib.pyplot as plt

with h5py.File('results/results.h5', 'r') as f:
    if 'aggregated/kinetics/derivatives' in f:
        plt.figure(dpi=144)
        dset = f['aggregated']['kinetics']['derivatives']
        data = dset[:]
        for i in range(data.shape[1]-1):
            label = dset.attrs['column_labels'][i+1]
            if '*' not in label:
                plt.plot(data[:, 0], data[:, i + 1], '-o', label=label)
        plt.grid(linestyle='--')
        plt.legend()
        plt.xlabel('Temperature [K]')
        plt.ylabel('TOF [-]')

The graph produced by the above script is shown below.

Derivatives as function of temperature

The graph shows that reactants are consumed and products are formed according to their stoichiometric ratios, with A : B₂ : AB = 1 : 1/2 : 1. The turnover frequencies (TOF) exhibit a maximum as a function of temperature, which is a typical result in heterogeneous catalysis and reflects the Sabatier principle. At low temperatures, the available thermal energy is insufficient to overcome reaction barriers. At high temperatures, adsorbates desorb rapidly instead of remaining on the surface to react. As a result, the highest TOF is achieved only at intermediate temperatures.

Note

Comparing this result to the concentrations profile, observe that this optimum in reactivity corresponds to the situation where A* and * are roughly in equal ratios.

Reaction orders

Both experimentally and through simulation, an important kinetic property that can be determined is the reaction order. The reaction order describes how the overall reaction rate responds to changes in the partial pressure of a gas-phase species. In other words, it measures the sensitivity of the forward rate to variations in the concentration of a particular reactant.

Mathematically, the reaction order with respect to species x is defined as

\[ n_{x} = p_{x} \, \frac{\partial \ln r^{+}}{\partial p_{x}} \]

where

\(p_{x}\) is the partial pressure of species x
\(r^{+}\) is the forward reaction rate

A positive reaction order indicates that increasing the pressure of species x accelerates the reaction, while a negative value indicates an inhibitory effect. Reaction orders are often non-integer values in heterogeneous catalysis, reflecting the complex interplay between adsorption, surface reactions, and desorption.

import h5py
import matplotlib.pyplot as plt

with h5py.File('results/results.h5', 'r') as f:
    if 'aggregated/kinetics/derivatives' in f:
        plt.figure(dpi=144)
        dset = f['aggregated']['sensitivity']['orders']
        data = dset[:]
        for i in range(data.shape[1]-1):
            label = dset.attrs['column_labels'][i+1]
            if '*' not in label:
                plt.plot(data[:, 0], data[:, i + 1], '-o', label=label)
        plt.grid(linestyle='--')
        plt.legend()
        plt.xlabel('Temperature [K]')
        plt.ylabel('Reaction order [-]')

The graph produced by the script above looks as follows

Reaction order as function of temperature

The graph shows that at low temperatures the reaction order with respect to A is negative (–1). This means that increasing the partial pressure of A actually inhibits the overall reaction rate, which can be understood as a consequence of surface poisoning: excess adsorption of A blocks active sites that would otherwise participate in the formation of products. As the temperature increases, the reaction order in A rises and eventually reaches +1. In this regime, higher partial pressures of A promote the reaction, consistent with A becoming the rate-controlling reactant.

For B₂, the reaction order remains constant at 0.5 across the entire temperature range. This fractional order suggests that the rate depends on the square root of the B₂ pressure, which is typical for dissociative adsorption processes where two surface sites are required. For AB, the reaction order is 0.0, indicating that the rate is independent of the partial pressure of AB. This behavior is expected for a product, as its presence in the gas phase does not directly influence the forward rate of reaction.

Apparent activation energy

Another key kinetic property that can be obtained from both experiment and simulation is the apparent activation energy. The apparent activation energy characterizes how the overall reaction rate changes with temperature. Unlike the intrinsic activation energy of an individual elementary step, the apparent activation energy reflects the combined influence of all relevant processes, including adsorption, surface reactions, and desorption. As such, it provides valuable insight into which steps are rate-controlling under given operating conditions.

Mathematically, the apparent activation energy is defined as

\[ E_\text{app} = RT \, \frac{\partial \ln r^{+}}{\partial T} \]

where

\(R\) is the universal gas constant
\(T\) is the temperature
\(r^{+}\) is the forward reaction rate

A higher apparent activation energy indicates stronger temperature dependence of the reaction rate, while a lower value suggests that other effects, such as adsorption equilibria, are compensating for the intrinsic barrier. In heterogeneous catalysis, the apparent activation energy can therefore differ significantly from the activation energy of the elementary surface reaction, and its variation with conditions often reveals mechanistic details about the catalytic cycle.

Via the Python script as shown below, we can visualize the apparent activation energy

import h5py
import matplotlib.pyplot as plt

with h5py.File('results/results.h5', 'r') as f:
    if 'aggregated/kinetics/derivatives' in f:
        plt.figure(dpi=144)
        dset = f['aggregated']['sensitivity']['eapp']
        data = dset[:]
        for i in range(data.shape[1]-1):
            label = dset.attrs['column_labels'][i+1]
            plt.plot(data[:, 0], data[:, i + 1] / 1000, '-o')
        plt.grid(linestyle='--')
        plt.xlabel('Temperature [K]')
        plt.ylabel('Apparent activation energy [kJ/mol]')

which yields the following graph

Apparent activation energy as function of temperature

The graph shows that at low temperature the apparent activation energy is approximately 200 kJ/mol. This value can be rationalized by considering the elementary steps required for the reaction to proceed under these conditions. Since the surface is almost fully covered by A, one A must first be removed, which requires about 120 kJ/mol. The subsequent adsorption of half a B₂ molecule stabilizes the system, releasing roughly 40 kJ/mol. Together with the intrinsic reaction barrier of 120 kJ/mol, this leads to an effective apparent activation energy close to 200 kJ/mol.

As the temperature increases, the surface becomes less saturated with A*, and more vacant sites are available. This lowers the apparent activation energy, reflecting a kinetically more favorable regime where surface crowding is no longer the primary limitation.

At the optimal reaction temperature, the apparent activation energy approaches zero. This corresponds to the maximum in the turnover frequency (TOF), where the rate is locally insensitive to changes in temperature—analogous to probing the extremum of a function.

At even higher temperatures, the apparent activation energy becomes negative. In this regime, increasing the temperature reduces the overall reaction rate because adsorbates desorb too rapidly to participate effectively in surface reactions. Consequently, lowering the temperature actually enhances the rate. This inversion is a hallmark of the Sabatier principle, which predicts that maximum catalytic activity occurs only within an intermediate temperature (or binding energy) window.

Degree of rate control

Another valuable kinetic descriptor is the Degree of Rate Control (DRC).
The DRC quantifies how sensitive the overall reaction rate is to changes in the rate constant of an individual elementary step. In other words, it measures how much "control" each step exerts over the overall kinetics of the catalytic cycle.

Mathematically, the DRC of step i is defined as

\[ \chi_i = \frac{k_i}{r} \, \frac{\partial r}{\partial k_i} = \frac{\partial \ln r}{\partial \ln k_{i}} \]

where

\(k_i\) is the rate constant (both forward and backward) of elementary step i
\(r\) is the overall reaction rate

The DRC takes values between 0 and 1 in many practical cases. A value close to 1 indicates that the step is strongly rate-controlling: small changes in its rate constant directly scale the overall rate. A value close to 0 means the step has little influence on the kinetics under the given conditions. Interestingly, negative values can also occur, implying that accelerating that step would actually reduce the overall reaction rate (e.g. by shifting equilibria or surface coverages unfavorably).

The DRC thus provides a powerful tool for mechanistic analysis: it identifies which elementary steps are kinetically relevant, guides catalyst design by pointing to critical barriers, and helps explain why apparent kinetic parameters vary with conditions.

The following Python script extracts the DRC results from results.h5 and visualizes it.

import h5py
import matplotlib.pyplot as plt

with h5py.File('results/results.h5', 'r') as f:
    if 'aggregated/kinetics/derivatives' in f:
        plt.figure(dpi=144)
        dset = f['aggregated']['sensitivity']['drc']
        data = dset[:]
        for i in range(data.shape[1]-1):
            label = dset.attrs['column_labels'][i+1]
            plt.plot(data[:, 0], data[:, i + 1], '-o', label=label)
        plt.grid(linestyle='--')
        plt.legend()
        plt.xlabel('Temperature [K]')
        plt.ylabel('Degree of rate control [-]')

DRC as function of temperature

From the graph above, we can see that only one elementary step—the surface reaction—has a significant impact on the overall kinetics. Its DRC value is equal to 1, which means that this step fully controls the reaction rate under the given conditions. In practical terms, any change to the rate constant of this step (for example by lowering its activation barrier) would lead to a directly proportional change in the overall rate.

All other steps, such as adsorption and desorption, have DRC values close to 0. This indicates that variations in their rate constants do not appreciably influence the overall kinetics. These steps are equilibrated under the current conditions and therefore play no decisive role in limiting the reaction rate.

The identification of a single step with a DRC of 1 highlights it as the rate-determining step. This is consistent with the classical interpretation of catalytic cycles: while many processes occur simultaneously, only the slowest, kinetically dominant step governs the observable rate. From a catalyst design perspective, this also means that efforts to improve performance should be targeted at reducing the barrier of this particular elementary step.