Kinetic Models

This chapter documents the kinetic models currently covered by the MKMCXX 3 manual.

Included Models

• ArrheniusDefault
• ArrheniusFreq
• ArrheniusQRatio
• HertzKnudsenDefault
• HertzKnudsenNist

ArrheniusDefault

This model is intended for elementary surface reactions with explicit forward and backward prefactors and activation barriers.

Parameters

• Vf Forward pre-exponential factor \(\nu_f\), in \(\mathrm{s^{-1}}\).
• Vb Backward pre-exponential factor \(\nu_b\), in \(\mathrm{s^{-1}}\).
• Eaf Forward activation energy \(E_{a,f}\), in \(\mathrm{J/mol}\).
• Eab Backward activation energy \(E_{a,b}\), in \(\mathrm{J/mol}\).

Equations

\[ k_f = \nu_f \exp\left(-\frac{E_{a,f}}{RT}\right) \]

with:

• \(\nu_f\): forward pre-exponential factor, in \(\mathrm{s^{-1}}\)
• \(E_{a,f}\): forward activation energy, in \(\mathrm{J/mol}\)
• \(R\): molar gas constant, in \(\mathrm{J/(mol\,K)}\)
• \(T\): absolute temperature, in \(\mathrm{K}\)

\[ k_b = \nu_b \exp\left(-\frac{E_{a,b}}{RT}\right) \]

with:

• \(\nu_b\): backward pre-exponential factor, in \(\mathrm{s^{-1}}\)
• \(E_{a,b}\): backward activation energy, in \(\mathrm{J/mol}\)
• \(R\): molar gas constant, in \(\mathrm{J/(mol\,K)}\)
• \(T\): absolute temperature, in \(\mathrm{K}\)

Example

reactions {
  {A*}+{B*}=>{AB*}+{*} @ ArrheniusDefault(
    Vf=1e13,    # forward pre-exponential factor [s^-1]
    Vb=1e13,    # backward pre-exponential factor [s^-1]
    Eaf=120e3,  # forward activation energy [J/mol]
    Eab=80e3    # backward activation energy [J/mol]
  )
}

ArrheniusQRatio

This model is intended for reactions parameterized through partition-function ratios.

Parameters

• QTSQIS Ratio \(Q_{\mathrm{TS}}/Q_{\mathrm{IS}}\) between the transition-state and initial-state partition functions, dimensionless.
• QTSQFS Ratio \(Q_{\mathrm{TS}}/Q_{\mathrm{FS}}\) between the transition-state and final-state partition functions, dimensionless.
• Eaf Forward activation energy \(E_{a,f}\), in \(\mathrm{J/mol}\).
• Eab Backward activation energy \(E_{a,b}\), in \(\mathrm{J/mol}\).

Equations

In the current implementation, the forward and backward rate constants are given by:

\[ k_f = \frac{k_b T}{h} Q_{\mathrm{TS}}/Q_{\mathrm{IS}} \exp\left(-\frac{E_{a,f}}{RT}\right) \]

with:

• \(k_b\): Boltzmann constant, in \(\mathrm{J/K}\)
• \(T\): absolute temperature, in \(\mathrm{K}\)
• \(h\): Planck constant, in \(\mathrm{J\,s}\)
• \(Q_{\mathrm{TS}}/Q_{\mathrm{IS}}\): forward partition-function ratio, dimensionless
• \(E_{a,f}\): forward activation energy, in \(\mathrm{J/mol}\)
• \(R\): molar gas constant, in \(\mathrm{J/(mol\,K)}\)

\[ k_b = \frac{k_b T}{h} Q_{\mathrm{TS}}/Q_{\mathrm{FS}} \exp\left(-\frac{E_{a,b}}{RT}\right) \]

with:

• \(k_b\): Boltzmann constant, in \(\mathrm{J/K}\)
• \(T\): absolute temperature, in \(\mathrm{K}\)
• \(h\): Planck constant, in \(\mathrm{J\,s}\)
• \(Q_{\mathrm{TS}}/Q_{\mathrm{FS}}\): backward partition-function ratio, dimensionless
• \(E_{a,b}\): backward activation energy, in \(\mathrm{J/mol}\)
• \(R\): molar gas constant, in \(\mathrm{J/(mol\,K)}\)

In MKMCXX input syntax, these partition-function ratios are supplied through:

• QTSQIS for \(Q_{\mathrm{TS}}/Q_{\mathrm{IS}}\)
• QTSQFS for \(Q_{\mathrm{TS}}/Q_{\mathrm{FS}}\)

Example

reactions {
  {CO*}+{*}=>{C*}+{O*} @ ArrheniusQRatio(
    QTSQIS=2.78e-01,  # transition-state / initial-state partition-function ratio [-]
    QTSQFS=4.43e-01,  # transition-state / final-state partition-function ratio [-]
    Eaf=68.4e3,       # forward activation energy [J/mol]
    Eab=71.7e3        # backward activation energy [J/mol]
  )
}

ArrheniusFreq

This model evaluates transition-state-theory prefactors from vibrational frequency files. Stable-state frequency files are attached to components through the freqfile attribute, while the reaction argument freqfile provides the transition-state frequencies.

The intended workflow is:

• assign stable-state vibrational data to components through freqfile
• assign the transition-state vibrational data to the reaction through freqfile
• provide Eaf and Eab as electronic forward and backward barriers in \(\mathrm{J/mol}\), excluding the zero-point-energy correction

At runtime, MKMCXX combines those inputs into harmonic transition-state-theory rate constants with an explicit vibrational partition-function ratio and an explicit ZPE correction.

Note

ArrheniusFreq expects Eaf and Eab to be electronic barriers only. The zero-point-energy contribution is added internally from the referenced frequency files at runtime. This differs from simpler Arrhenius-style parameterizations where the reported activation energy may already include the ZPE correction.

Parameters

• freqfile Transition-state vibrational frequency file for this reaction.
• Eaf Forward activation energy \(E_{a,f}\) excluding the ZPE correction, in \(\mathrm{J/mol}\).
• Eab Backward activation energy \(E_{a,b}\) excluding the ZPE correction, in \(\mathrm{J/mol}\).

Equations

For a reaction with initial state \(IS\), transition state \(TS\), and final state \(FS\), the forward and backward rate constants are evaluated as:

\[ k_f = \frac{k_b T}{h} \frac{Q^\mathrm{vib}_{TS}}{Q^\mathrm{vib}_{IS}} \exp\left(-\frac{E_{a,f} + \Delta E^\mathrm{ZPE}_f}{RT}\right) \]

with:

• \(k_b\): Boltzmann constant, in \(\mathrm{J/K}\)
• \(T\): absolute temperature, in \(\mathrm{K}\)
• \(h\): Planck constant, in \(\mathrm{J\,s}\)
• \(Q^\mathrm{vib}_{TS}/Q^\mathrm{vib}_{IS}\): forward vibrational partition-function ratio, dimensionless
• \(E_{a,f}\): forward electronic activation energy, in \(\mathrm{J/mol}\)
• \(\Delta E^\mathrm{ZPE}_f\): forward zero-point-energy correction, in \(\mathrm{J/mol}\)
• \(R\): molar gas constant, in \(\mathrm{J/(mol\,K)}\)

\[ k_b = \frac{k_b T}{h} \frac{Q^\mathrm{vib}_{TS}}{Q^\mathrm{vib}_{FS}} \exp\left(-\frac{E_{a,b} + \Delta E^\mathrm{ZPE}_b}{RT}\right) \]

with:

• \(k_b\): Boltzmann constant, in \(\mathrm{J/K}\)
• \(T\): absolute temperature, in \(\mathrm{K}\)
• \(h\): Planck constant, in \(\mathrm{J\,s}\)
• \(Q^\mathrm{vib}_{TS}/Q^\mathrm{vib}_{FS}\): backward vibrational partition-function ratio, dimensionless
• \(E_{a,b}\): backward electronic activation energy, in \(\mathrm{J/mol}\)
• \(\Delta E^\mathrm{ZPE}_b\): backward zero-point-energy correction, in \(\mathrm{J/mol}\)
• \(R\): molar gas constant, in \(\mathrm{J/(mol\,K)}\)

with

\[ \Delta E^\mathrm{ZPE}_f = E^\mathrm{ZPE}_{TS} - E^\mathrm{ZPE}_{IS} \]

\[ \Delta E^\mathrm{ZPE}_b = E^\mathrm{ZPE}_{TS} - E^\mathrm{ZPE}_{FS} \]

The vibrational partition functions of the stable states are assembled from the reaction stoichiometry. If the initial state contains components \(i\) with stoichiometric coefficients \(\nu_i^\mathrm{IS}\), then

\[ Q^\mathrm{vib}_{IS} = \prod_i \left(q_i^\mathrm{vib}\right)^{\nu_i^\mathrm{IS}} \]

and similarly for the final state:

\[ Q^\mathrm{vib}_{FS} = \prod_j \left(q_j^\mathrm{vib}\right)^{\nu_j^\mathrm{FS}} \]

The transition-state contribution comes directly from the reaction-level frequency file:

\[ Q^\mathrm{vib}_{TS} = q^\mathrm{vib}_{TS} \]

For one vibrational frequency file with real vibrational modes \(\tilde{\nu}_m\) in \(\mathrm{cm}^{-1}\), MKMCXX evaluates:

\[ q^\mathrm{vib}(T) = \prod_m \frac{1}{1-\exp\left(-\frac{h c \tilde{\nu}_m}{k_b T}\right)} \]

\[ E^\mathrm{ZPE} = \frac{1}{2} \sum_m h c \tilde{\nu}_m N_A \]

where \(N_A\) is Avogadro's constant so that the ZPE contribution is expressed in \(\mathrm{J/mol}\), consistent with \(E_{a,f}\), \(E_{a,b}\), and \(RT\).

Cutoff Treatment

To avoid unphysically large low-frequency vibrational contributions, MKMCXX can apply a simulation-level vibrational cutoff:

\[ \tilde{\nu}_m^\mathrm{eff} = \max(\tilde{\nu}_m, \tilde{\nu}_\mathrm{min}) \]

where \(\tilde{\nu}_\mathrm{min}\) is provided through vib_cutoff in \(\mathrm{cm}^{-1}\). The same effective frequency is used consistently in both the partition function and ZPE evaluation.

Notes

• A component without a freqfile contributes no vibrational correction in this model, which is equivalent to \(q^\mathrm{vib}=1\) and \(E^\mathrm{ZPE}=0\).
• Empty surface sites can therefore be left without a vibrational file.
• Eaf and Eab should be provided without a ZPE correction, because that correction is inferred from the referenced frequency files at runtime.
• The current implementation uses the real frequencies listed under vibrations.frequencies_cm-1 in the YAML files.

Example

components {
  CO* {phase=surface, init=0.0, freqfile="ISFS/co.yaml"}  # stable-state frequencies for CO*
  C*  {phase=surface, init=0.0, freqfile="ISFS/c.yaml"}   # stable-state frequencies for C*
  O*  {phase=surface, init=0.0, freqfile="ISFS/o.yaml"}   # stable-state frequencies for O*
  *   {phase=surface, init=1.0}                           # vacant surface site
}

reactions {
  {CO*}+{*}=>{C*}+{O*} @ ArrheniusFreq(
    freqfile="TS/co_diss.yaml",  # transition-state vibrational frequencies
    Eaf=68.4e3,                  # forward electronic activation energy [J/mol]
    Eab=71.7e3                  # backward electronic activation energy [J/mol]
  )
}

simulation {
  baseline @ transient_solve(
    model=thermal,  # temperature-controlled simulation mode
    vib_cutoff=100, # minimum vibrational frequency used in partition functions [cm^-1]
    solver={method=cvode_bdf, atol=1e-12, rtol=1e-10, t_end=1e4},
    conditions=[{T=600}]  # simulation temperature [K]
  )
}

In this example, MKMCXX automatically builds:

• \(Q^\mathrm{vib}_{IS}\) from the reactant-side component frequency files
• \(Q^\mathrm{vib}_{FS}\) from the product-side component frequency files
• \(Q^\mathrm{vib}_{TS}\) from TS/co_diss.yaml

and then evaluates the forward and backward rate constants from the equations above.

HertzKnudsenDefault

This model is intended for adsorption and desorption steps using the Hertz-Knudsen formulation.

Parameters

• Asite Active site area per adsorption site \(A_{\text{site}}\), in \(\mathrm{m^2}\).
• m Mass of the adsorbing species \(m\), in \(\mathrm{amu}\).
• theta Rotational temperature \(\theta_{\text{rot}}\) used in the desorption expression, in \(\mathrm{K}\).
• sigma Symmetry number \(\sigma\) used in the desorption expression, dimensionless.
• S Sticking coefficient \(S\) for adsorption, dimensionless.
• Edes Desorption energy \(E_{\text{des}}\), in \(\mathrm{J/mol}\).

Equations

In the current implementation, the forward adsorption rate constant is:

\[ k_f = \frac{A_{\text{site}} S}{\sqrt{2 \pi m \, \mathrm{amu} \, k_b T}} \; p_{\mathrm{ref}} \]

with:

• \(A_{\text{site}}\) in \(\mathrm{m^2}\)
• \(m\) in \(\mathrm{amu}\)
• \(k_b\) in \(\mathrm{J/K}\)
• \(T\) in \(\mathrm{K}\)
• \(S\) dimensionless
• \(p_{\mathrm{ref}} = 101325 \; \mathrm{Pa/atm}\)

The factor \(p_{\mathrm{ref}}\) is the explicit pressure-conversion factor used by the code. This is relevant because gas pressures are provided in atm in MKMCXX inputs, while the Hertz-Knudsen flux expression itself is written in SI units.

\[ k_b = \frac{k_bT^3}{h^3}\frac{A_{\text{site}}(2\pi m \, \mathrm{amu} \, k_b)}{\sigma\theta_{\text{rot}}} \exp\left(-\frac{E_{\text{des}}}{RT}\right) \]

with:

• \(A_{\text{site}}\) in \(\mathrm{m^2}\)
• \(m\) in \(\mathrm{amu}\)
• \(k_b\) in \(\mathrm{J/K}\)
• \(T\) in \(\mathrm{K}\)
• \(h\) in \(\mathrm{J\,s}\)
• \(\sigma\) dimensionless
• \(\theta_{\text{rot}}\) in \(\mathrm{K}\)
• \(E_{\text{des}}\) in \(\mathrm{J/mol}\)
• \(R\) in \(\mathrm{J/(mol\,K)}\)

Example

reactions {
  {A}+{*}=>{A*} @ HertzKnudsenDefault(
    Asite=1e-20,  # active site area [m^2]
    m=12,         # gas-phase molecular mass [amu]
    theta=1.0,    # rotational temperature [K]
    sigma=1,      # rotational symmetry number [-]
    S=1,          # sticking coefficient [-]
    Edes=120e3    # desorption energy [J/mol]
  )
}

HertzKnudsenNist

This model extends the Hertz-Knudsen approach with thermochemical data from the NIST-style parameterization.

Parameters

• Asite Active site area per adsorption site \(A_{\text{site}}\), in \(\mathrm{m^2}\).
• m Mass of the adsorbing species \(m\), in \(\mathrm{amu}\).
• S Sticking coefficient \(S\) for adsorption, dimensionless.
• Qads Adsorbate partition-function contribution \(Q_{\mathrm{ads}}\) used in the desorption term, dimensionless.
• Edes Desorption energy \(E_{\mathrm{des}}\), in \(\mathrm{J/mol}\).
• A Shomate coefficient \(A\), in \(\mathrm{J/(mol\,K)}\).
• B Shomate coefficient \(B\), in \(\mathrm{J/(mol\,K\,kK)}\).
• C Shomate coefficient \(C\), in \(\mathrm{J/(mol\,K\,kK^2)}\).
• D Shomate coefficient \(D\), in \(\mathrm{J/(mol\,K\,kK^3)}\).
• E Shomate coefficient \(E\), in \(\mathrm{J\,kK^2/(mol\,K)}\).
• F Shomate coefficient \(F\), in \(\mathrm{kJ/mol}\).
• G Shomate coefficient \(G\), in \(\mathrm{J/(mol\,K)}\).
• H Shomate coefficient \(H\), in \(\mathrm{kJ/mol}\).

Equations

In the current implementation, the forward adsorption rate constant is:

\[ k_f = \frac{A_{\text{site}} S}{\sqrt{2 \pi m \, \mathrm{amu} \, k_b T}} \; p_{\mathrm{ref}} \]

with:

• \(A_{\text{site}}\) in \(\mathrm{m^2}\)
• \(m\) in \(\mathrm{amu}\)
• \(k_b\) in \(\mathrm{J/K}\)
• \(T\) in \(\mathrm{K}\)
• \(S\) dimensionless
• \(p_{\mathrm{ref}} = 101325 \; \mathrm{Pa/atm}\)

The factor \(p_{\mathrm{ref}}\) is the explicit pressure-conversion factor used by the code. This is relevant because gas pressures are provided in atm in MKMCXX inputs, while the Hertz-Knudsen flux expression itself is written in SI units.

For the backward rate constant, MKMCXX first evaluates the gas-phase entropy and enthalpy from the Shomate-style parameters A through H using \(t = T / 1000\):

\[ S_{\mathrm{gas}} = A \ln t + Bt + \frac{Ct^2}{2} + \frac{Dt^3}{3} - \frac{E}{2t^2} + G \]

where \(S_{\mathrm{gas}}\) is in \(\mathrm{J/(mol\,K)}\).

\[ H_{\mathrm{gas}} = At + \frac{Bt^2}{2} + \frac{Ct^3}{3} + \frac{Dt^4}{4} - \frac{E}{t} + F - H \]

where \(H_{\mathrm{gas}}\) is in \(\mathrm{kJ/mol}\) and \(t\) is the reduced temperature variable in \(\mathrm{kK}\).

The backward desorption rate constant is then computed as:

\[ k_b = \frac{k_f}{Q_{\mathrm{ads}}} \exp\left(\frac{S_{\mathrm{gas}}}{R}\right) \exp\left(-\frac{E_{\mathrm{des}} + 1000 H_{\mathrm{gas}}}{RT}\right) \]

with \(R\) in \(\mathrm{J/(mol\,K)}\), \(E_{\mathrm{des}}\) in \(\mathrm{J/mol}\), and the factor \(1000\) converting \(H_{\mathrm{gas}}\) from \(\mathrm{kJ/mol}\) to \(\mathrm{J/mol}\) before it is combined with \(E_{\mathrm{des}}\).

Notes

• Qads corresponds to the adsorbate partition-function contribution \(Q_{\mathrm{ads}}\) used in the desorption expression.
• The Shomate contribution is evaluated from the supplied coefficients \(A\) through \(H\).
• This documentation follows the implementation in src/kinetics/reactions/hertz_knudsen_nist.cpp.

Example

reactions {
  {CO}+{*}=>{CO*} @ HertzKnudsenNist(
    Asite=1e-19,  # active site area [m^2]
    m=28,         # gas-phase molecular mass [amu]
    S=1,          # sticking coefficient [-]
    A=25.56759,   # Shomate coefficient A [J/(mol K)]
    B=6.09613,    # Shomate coefficient B [J/(mol K kK)]
    C=4.054656,   # Shomate coefficient C [J/(mol K kK^2)]
    D=-2.671301,  # Shomate coefficient D [J/(mol K kK^3)]
    E=0.131021,   # Shomate coefficient E [J kK^2/(mol K)]
    F=-118.0089,  # Shomate coefficient F [kJ/mol]
    G=227.3665,   # Shomate coefficient G [J/(mol K)]
    H=-110.5271,  # Shomate coefficient H [kJ/mol]
    Qads=1,       # adsorbate partition-function contribution [-]
    Edes=80e3     # desorption energy [J/mol]
  )
}