Microkinetic Modeling of Nanoparticle Catalysis Using Density Functional Theory

Thesis for the degree of Licentiate of Engineering
Microkinetic Modeling of Nanoparticle Catalysis

using Density Functional Theory
Mikkel Jørgensen
Department of Physics
Chalmers University of Technology
Gothenburg, Sweden 2017
Microkinetic Modeling of Nanoparticle Catalysis using Density Functional Theory
Mikkel Jørgensen
c Mikkel Jørgensen, 2017

SE-412 96 Gothenburg
Sweden
Telephone: +46 (0)31-772 1000
Cover:
Kinetic simulations, represented by dices, close the materials gap between extended
surfaces and nanoparticles.
Printed at Chalmers Reproservice

Gothenburg, Sweden 2017
Microkinetic Modeling of Nanoparticle Catalysis using Density Functional Theory
Mikkel Jørgensen
Abstract
Heterogeneous catalysis is vitally important to modern society, and one path towards
rational catalyst design is through atomistic scale understanding. The atomistic scale
can be linked to macroscopic observables by microkinetic models based on first-principles
calculations. With the increasing accuracy of first-principles methods and growing com-
putational resources, it has become important to investigate and further develop the
methodology of microkinetic modeling, which is the theme of this thesis.
First, a procedure for mean-field microkinetic modeling of reactions over extended surfaces
is developed, where complete methane oxidation over Pd(100) and Pd(111) is studied as
an example. The model reveals how the main reaction mechanisms depend on reaction
conditions, and shows poisoning as well as promotion phenomena.
Second, the effect of entropy in microkinetic modeling is investigated, where CO oxidation
over Pt(111) is used as a model reaction. Entropy is found to affect reaction kinetics
substantially. Moreover, a method named Complete Potential Energy Sampling (CPES)
is developed as a flexible tool for estimating adsorbate-entropy.
Third, a kinetic Monte Carlo method is developed to bridge the materials gap in het-
erogeneous catalysis. The computational cost to map out the complete reaction-energy-
landscape on a nanoparticle is high, which is solved herein using generalized coordination
numbers as descriptors for reaction energies. CO oxidation over Pt is studied, and
nanoparticles are found to behave differently than the corresponding extended surfaces.
Moreover, the active site is found to vary with reaction conditions.
Finally, the reaction orders and apparent activation energies are coupled to the microscale
via the degree of rate control, which enhances the atomistic understanding of reaction
kinetics.
Keywords: Catalysis, Density Functional Theory, Microkinetic modeling, Nanoparti-

cles, Kinetic Monte Carlo, Mean-field approximation, Entropy, Methane oxidation, CO
oxidation.
i
ii
List of Publications
This thesis is based on the following appended papers:
I.
First-Principles Microkinetic Modeling of
Methane Oxidation over Pd(100) and Pd(111)
M. Jørgensen and H. Grönbeck
ACS Catalysis, 6 (2016), 6730-6738
II.
Adsorbate Entropies with Complete Potential
Energy Sampling in Microkinetic Modeling
The Journal of Physical Chemistry C, 121 (2017), 7199-7207
III.
Scaling Relations and Kinetic Monte Carlo
Simulations to Bridge the Materials Gap
in Heterogeneous Catalysis
Accepted, ACS Catalysis, (2017), doi: 10.1021/acscatal.7b01194
IV.
Connection between Macroscopic Kinetic
Measurables and the Degree of Rate Control
Submitted, Catalysis Science & Technology, (2017).
iii
iv
My contributions to the publications
Paper I
I programmed the kinetic mean-field code and performed all the calculations. I wrote the
first draft of the paper, which was finalized together with my coauthor.
Paper II
I programmed the kinetic mean-field code and performed all the calculations. I wrote the
first draft of the paper, which was finalized together with my coauthor.
Paper III
I programmed the kinetic Monte Carlo code and performed all the calculations. I wrote
the first draft of the paper, which was finalized together with my coauthor.
Paper IV
I derived the equations, and performed the numerical simulations. I wrote the first draft
of the manuscript, which was finalized together with my coauthor.
v
vi
Contents
1 Introduction 1
1.1 Thesis Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2 Calculating Rate Constants 5

2.1 Reaction Rate Constants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Density Functional Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2.1 Exchange Correlation Energy Approximations . . . . . . . . . . . . . . . . 8
2.2.2 Projector-Augmented Wave Method . . . . . . . . . . . . . . . . . . . . . . 9
2.2.3 Periodic Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Local Minimum Energies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.4 Energy Barriers with Nudged Elastic Band . . . . . . . . . . . . . . . . . . . 11
2.5 Adsorbate-Adsorbate Interactions . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.6 Modeling Entropies in Catalytic Reactions . . . . . . . . . . . . . . . . . . . . 13
2.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3 Simulating Reaction Kinetics 17

3.1 Chemical Master Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.2 Mean-Field Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.3 Kinetic Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.4 Analyzing Reaction Kinetics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.4.1 Reaction Mechanisms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.4.2 Turnover Frequencies and Coverages . . . . . . . . . . . . . . . . . . . . . . 23
3.4.3 Reaction Orders and Apparent Activation Energies . . . . . . . . . . . . . . 23
3.4.4 Degree of Rate Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4 Kinetic Models for Methane and CO Oxidation 27

4.1 Constructing a Kinetic Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.2 Complete Methane Oxidation over Pd(100) and Pd(111) . . . . . . . . . . . . 28
4.3 Entropic Effects in CO Oxidation on Pt(111) . . . . . . . . . . . . . . . . . . 31
4.4 Catalysis on Nanoparticles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5 Conclusion & Outlook 37
Acknowledgments 39
vii
viii
Chapter 1
Introduction
”What would happen if we could arrange the atoms one by one the way we want them ... I
can hardly doubt that when we have some control over the arrangement of things on a small
scale we will get an enormously greater range of possible properties that substances can
have and of different things we can do ... Atoms on a small scale behave like nothing on a
large scale, for they satisfy the laws of quantum mechanics.” Richard P. Feynman(1959) 1 .
This was how Richard Feynman envisioned the influence that nanotechnology will have
in society, nearly 60 years ago. Today, his visions of manipulating matter at the atomic
scale is close to reality. Systems can be designed atom by atom, which is a consequence
of understanding the quantum mechanical behavior of atomic systems. Understanding
the behavior of matter at the smallest scale is crucial as it enables rational design, and
manipulation, of properties in nanoparticle-based technologies.
Figure 1.1: Schematic illustration of a catalytic reaction. The catalyst provides an

alternative pathway with lower energy barriers for a chemical reaction.
One technology that benefits from atomistic understanding is heterogeneous catalysis.

Heterogeneous catalysis has a range of applications in society, such as cleaning automotive
exhaust from harmful gases, producing synthetic fuels, and manufacturing chemicals. In
the chemical industry, about 90% of the production applies catalysts 2 (p. 11) . Because
of this, small improvements can have major impacts on a global scale. A catalyst is a
substance that speeds up a chemical reaction without being consumed in the process. In
Figure 1.1, a catalytic reaction is illustrated, where reactants are converted to products
1
with lower energies as compared to the case without a catalyst. The conversion of
reactants proceeds on the surface of the catalyst, which provides an alternative pathway
for the reaction. To maximize the exposed surface area, heterogeneous catalysts are often
designed as precious metal-nanoparticles, dispersed on porous metal-oxide supports.
Since catalysis proceeds on the surface of nanoparticles, understanding bonding on
surfaces is required. The first step in a catalytic cycle is adsorption of reactants on the
surface. Next, the reactants form new chemical bonds mediated by the surface, and finally
the resulting products should desorb from the surface. For the catalytic cycle to proceed
efficiently, the adsorbate-surface bond should neither be too strong nor too weak. This
point is contained in the Sabatier principle 2 (p. 264) , which states that for a given reaction,
there is an optimal binding energy of the reactants. The fact that there is an optimal
binding energy has lead to the recognition that only a few sites might be active in a
catalytic reaction and the concept of active sites 3 . An early description to quantify the
rate of a catalytic reaction was the Arrhenius equation 2 (p. 36) , which gives an empirical
expression for the rate constant:

−Ea
k = A exp (1.1)
kB T
where Ea is the activation energy, A is a pre-exponential factor, and T is the temperature.

As Ea appears in the exponent, small variations in Ea results in large changes in the rate.
A smaller activation energy allows for a higher efficiency, since the temperature can be
kept low. The Arrhenius equation is connected to the Sabatier principle as the activation
energy usually is correlated with the bond strengths of the reactants. This is the essence
of the Brøndsted-Evans-Polanyi (BEP) principle 2(p. 267) , which states that the transition
state energy of a reaction depends on the bond strength of the adsorbates.
Besides the general and phenomenological insights contained in the Sabatier principle
and the Arrhenius equation, the microscopic picture of catalysis has been elucidated by
the computational chemistry and physics communities. The development of electronic
structure calculations; especially Density Functional Theory 4;5 , have played an important
role in studies of catalytic reactions 6 . Some of the earliest models of the catalyst surface
were small clusters that locally mimicked the active site 7 . These cluster models are
fair approximations to the catalyst surface 7 , which can explain chemical bonds between
adsorbates and surfaces qualitatively. Later, when the computational power had matured,
slab models were applied to model extended surfaces 7 , which presently is the most common
approach. Slab models have enabled quantitative descriptions of the adsorbate-surface
bond. Explicit calculations on nanoparticles still pose a special challenge due to the
system size and structural complexity. However, during the past few years, studies of
adsorption on idealized nanoparticle-models have been conducted using first-principles
calculations 8–17 .
In addition to studies of specific reactions, computational approaches have been used
to understand trends in catalytic properties. One example is the adsorption on transition
metal surfaces, where the leftmost metals in the periodic table of elements are more
reactive than the metals to the right. This trend was rationalized by the d-band model 18 ,
where the center and filling of the d-band determines the reactivity of a certain catalyst
site. Different sites on nanoparticles have different coordination numbers that correspond
2
to varying d-band centers. Thus, the d-band model suggests that the coordination number
should be a proper descriptor for binding energies 6 . To improve on the coordination
number as a descriptor, the coordination numbers of the first nearest neighbor can be
taken into account. This was recently realized with the generalized coordination number
that was shown to be a good descriptor for the reactivity on extended surfaces and
nanoparticles 19–21 .
Adsorption studies solely account for the thermodynamic aspects of catalysis. However,
to fully understand a reaction it is also central to describe the reaction kinetics. The
reaction kinetics can be elucidated using a microkinetic model, which is a mathematical
representation of a reaction as a set of intermediate steps. The benefit of constructing
a microkinetic model is that the main reaction mechanisms and kinetic bottlenecks are
accessible, which reveals the elementary steps that should be targeted in future catalyst
design. Previously, microkinetic models based on experimental measurements 22–28 and
first-principles calculations 29–43 have been formulated for various catalytic reactions.
Although this field is rapidly growing, the number of formulated microkinetic models is still
modest 44 . Models derived from first-principles calculations been formulated using several
different methods, and it is timely to investigate and develop the general methodology of
microkinetic modeling, especially to take advantage of the growing computational power.
Furthermore in the preceding literature, extended surfaces have often be used as
model-systems for nanoparticles. Using extended surfaces gives the benefit of structural
simplicity. However, it simultaneously gives rise to a, so-called, materials gap between
technical catalysis and model systems. Efforts have previously been made to bridge the
materials gap by performing experiments on surfaces and well-defined nanoparticles 45–48 .
Additionally, theoretical models have been constructed to understand the kinetics on
nanoparticles 38;40;49;50 . In these theoretical studies, the considered reaction energy
landscapes were simplified or the reactions schematic. Hence, there is a need to understand
nanoparticle catalysis using first-principles-based reaction energy landscapes.
1.1 Thesis Objectives

The purpose of this thesis is to investigate and develop the methodology of first-principles
microkinetic modeling of reactions on extended surfaces and nanoparticles. In particular,
it aims to investigate catalysis through mean-field modeling and kinetic Monte Carlo
simulations. Furthermore, methods for analysis of microkinetic simulations and comparison
with experimental results are discussed. Special emphasis is placed on the modeling of
adsorbate entropy, and how to bridge the materials gap between extended surfaces and
nanoparticles.
The necessary components for formulating, solving, and analyzing first-principles
microkinetic models are illustrated in Figure 1.2. The first step is to choose a simplified
model system of the catalyst. Simplification of the real catalyst is necessary as a technical
catalyst is complex and ill-defined. Moreover, working with a simple model system will
likely result in a clearer understanding of catalytic mechanisms. The choice of model
system can include extended surfaces as well as nanoparticles, which both will be treated
in the upcoming chapters. The next ingredient in the modeling is an energy landscape for
the reaction, which is the set of Gibbs free energy barriers for the considered elementary
3
Figure 1.2: Concept map of the thesis. The map illustrates necessary components
to model catalysis from first principles.
reactions. All Gibbs free energies in the present work are calculated using Density
Functional Theory (DFT). When the Gibbs free energy barriers have been obtained,
the rate constants can be calculated. Construction of the reaction energy landscape
and rate constants are treated in chapter 2. With access to the rate constants of all
elementary reactions, the kinetics can be simulated by solving the chemical master
equation. The master equation can in special cases be solved approximately in the mean-
field approximation, and when the mean-field approximation is not applicable, the kinetics
can be simulated with kinetic Monte Carlo. Kinetic Monte Carlo is necessary to simulate
the kinetics of reactions on nanoparticles. Simulation of reaction kinetics and analysis
of the results are treated in chapter 3. The results of the kinetic simulations should
be analyzed to understand the governing reaction mechanisms and to determine which
elementary steps that control the rate. Here, the rate controlling steps are quantified by
a degree of rate control analysis. Ultimately, the success of a first-principles microkinetic
model depends on how well it agrees with experiments. To compare with experiments, the
model can be used to simulate reaction orders and apparent activation energies. Examples
of the developed methodologies are presented in chapter 4, which describes kinetic models
of methane and CO oxidation.
4
Chapter 2
Calculating Rate Constants
To simulate reaction kinetics, it is necessary to calculate the rate constants of the

elementary steps in the microkinetic model. This chapter deals with calculations of rate
constants from first principles using Density Functional Theory.
2.1 Reaction Rate Constants

Rate constants can be estimated using Transition State Theory 51;52 (TST). The TST
rate constant for an elementary reaction is derived by dividing phase-space into a reactant
and product region by a specific dividing surface 52 . All points on the dividing surface
in phase-space define the transition state (R‡ ). Reactants (R) that cross this dividing
surface become products (P ). The main assumption of TST is that the reactant and
transition state are in equilibrium, and that the transition rate into the product region is
so low that the Boltzmann distribution remains intact 52 . The second assumption is that
the reactants cross the dividing surface only once. Thereby, the TST estimate of the rate
constant will always be larger than the actual rate constant. The reaction can be written
as:
R R‡ → P (2.1)
The TST rate constant is written as a frequency f of the molecules passing R‡ into P ,
multiplied with the ratio of the partition functions as 52 :
Z‡
k TST = f , (2.2)
ZR
where Z ‡ is the partition function in the transition state and Z R is the partition function
of the reactant state. f can be found using the Maxwell-Boltzmann speed distribution,
and assuming that R‡ has a finite width δ, such that Z ‡ = δZ ∗ . Using this, the rate
constant can be shown to take the following form 52 :
kB T Z∗
k TST = √ , (2.3)
2πµkB T Z R
where µ is an effective mass for the reaction coordinate.
An effective mass of the reaction coordinate is an abstract concept, and it is cumbersome
to work with. Therefore, it is often assumed that the reaction coordinate is a vibration
with frequency ν. Taking the classical limit of the partition function for this vibrational
mode, in the limit where kB T hν, gives a simple expression for the rate constant 51 :
kB T Z ‡
k TST = . (2.4)
h ZR
5
For a barrierless adsorption reaction, where the gas-phase reactant molecule impinges on
a substrate atom with area A, one can deduce the following expression 52
TST s0 pA
kads =√ , (2.5)
2πM kB T
where p is the pressure of the gas, A is the are of the adsorption site, M is the molecule
mass, and s0 is the sticking probability. If this expression is used for the adsorption rate
constant, the corresponding desorption rate constant kdes can be determined from the
equilibrium constant K:
−∆G
kads
K= = e kB T , (2.6)
kdes
where ∆G is the Gibbs free energy change upon adsorption. In this manner, thermody-
namic consistency is ensured.
The Gibbs free energy change is composed of an enthalpy (∆H) and an entropy (∆S):
∆G = ∆H − T ∆S (2.7)
The enthalpy change of the reaction is determined by electronic rehybridization, where
bonds are broken and new bonds are formed. Likewise, the entropy is determined by the
potential energy landscape. Therefore, it is essential to calculate the energy of atomic
systems, which can be achieved using Density Functional Theory.
2.2 Density Functional Theory

The total energy of a system of atoms is calculated using the stationary Schrödinger
equation:
Ĥ |Ψi = E |Ψi , (2.8)
where the many-body Hamiltonian Ĥ can be written:
Ĥ = ĤeZ + Ĥee + ĤZZ . (2.9)
The expression (2.9) includes contributions from electron-nuclei interactions (eZ), electron-
electron interactions (ee), and nuclei-nuclei interactions (ZZ). The corresponding
Schrödinger equation to (2.9) does not allow an analytical solution for more than one
electron, and must be solved by approximations.
One common approximation is the Born-Oppenheimer Approximation, which treats
the nuclei with classical physics and invokes the adiabatic approximation. The adiabatic
approximation states that electrons move on a much shorter time scale than the nuclei.
Therefore, the electrons instantaneously follow the nuclei during atomic motion without
undergoing transitions between stationary states. However, the stationary states vary
as the nuclei move. The adiabatic approximation holds if the ratio of the electron
to nuclei mass is small compared to the energy difference between different adiabatic
eigenstates 53 (p. 9) . The second part of the Born-Oppenheimer approximation treats
nuclei as classical, which can be rationalized due to their large thermal wavelengths (λT ).
6
λT is roughly the average de Broglie wavelength at a given temperature 53 (p.10) . Particles
with substantially larger separation than λT do not exhibit quantum phase coherence,
and the total nuclear wavefunction is simply the product of the individual wavefunctions.
When the Born-Oppenheimer approximation is invoked, it is sufficient to deal with the
electronic part of the Hamiltonian ĥe :
ĥe = T̂ + V̂ext + V̂ee = T̂e + V̂ext + V̂H + V̂XC (2.10)
Where T̂ is the kinetic energy operator for the electrons, V̂ext is the external classical
electrostatic field set up by the nuclei, and V̂ee describes the electron-electron interactions.
In the second equality V̂ee is split up into the classical Hartree energy V̂H and the
quantum mechanical exchange and correlation energy V̂XC . The Hartree energy is simply
the electrostatic energy of the field set up by the electrons, which has a non-physical
self-interaction that should be canceled by V̂XC .
The exchange energy stems from the indistinguishability principle and symmetrization
postulate of quantum mechanics, which state that a fermionic wavefunction is anti-
symmetric under particle permutation. These constraints modify the electronic energy
with the so-called exchange energy. Exchange will make electrons of identical spin less
likely to be close as two fermions cannot occupy the same quantum state. The correlation
energy is the energy related to correlation of electronic motion. In this sense, exchange
interaction is also a form of correlation.
Density Functional Theory is based on the two Hohenberg-Kohn (HK) theorems 4 . The
first HK theorem states that the ground state electronic density (ρ) uniquely determines
the external potential (V̂ext ). This implies that ρ also determines the ground state |Ψi
of the many body problem. The second HK theorem defines the variational energy as a
functional of ρ 4
Z
Ev [ρ] = hΨ[ρ]|T̂ + Ûee |Ψ[ρ]i + ρ(r)vext (r) dr (2.11)
Where vext is the external potential. It must be noted that Ψ depends on ρ and vice
versa, which implies that the problem must be solved self-consistently.
The expression for the kinetic energy of the real interacting system is unknown as
a functional of the density. This difficulty is circumvented in the Kohn-Sham scheme 5 ,
where one works with a fictitious non-interacting system that has the same density as the
real system. The benefit is that the ground state of a non-interacting system is simply
the Slater-determinant composed of the individual single particle states. In the fictitious
reference system, the Hamiltonian is 53 (p. 61)
N
X
ĤR = T̂R + vR (ri ) (2.12)
i=1
Where N is the number of electrons. The first term is the single particle non-interacting
kinetic energy operator, and vR is a reference potential which ensures that the density
coincides with the density of the interacting system. vR can be derived using variational
calculus with the constraint that ρ integrates to N ; this yields 53 (p. 63)
Z eXC
ρ(r0 ) δE
vR (r) = vext (r) + 0
dr + . (2.13)
|r − r | δρ(r)
7
The second term in (2.13) is the Hartree energy and the third term is a chemical potential,
which is formulated as a functional derivative with respect to ρ. E eXC is a modified
exchange-correlation energy that must be introduced when neglecting kinetic correlations
in the reference system. The Kohn-Sham equations are a set of equations defined by (2.12).
These equations must be solved self-consistently since vR depends on ρ and vice-versa.
Finally, the energy is evaluated by the Kohn-Sham functional 53 (p. 62)
Z ZZ
1 ρ(r)ρ(r0 ) eXC [ρ]
EKS [ρ] = TR + ρ(r)vext (r) dr + drdr0 + E (2.14)
2 |r − r0 |
The only remaining unknown quantity here is the exchange-correlation energy, which
determines the accuracy of DFT. It is interesting to note that DFT is an exact theory,
given the true exchange-correlation energy.
2.2.1 Exchange Correlation Energy Approximations

The approximation of the exchange-correlation functional is a major concern when it
comes to the accuracy of DFT. One can view the exchange-correlation energy as the
Coulomb interaction between the density and a displaced charge density, which is brought
about by exchange and correlation 53 (p. 69) . This displaced charge density is called
the exchange-correlation hole ρeXC (r, r0 ) = ρ(r0 ) [g(r, r0 ) − 1]. Here g(r, r0 ) is the pair-
correlation function, describing the probability to find an electron at r given that one is
already present at r0 .
The simplest and earliest kind of exchange-correlation approximation is the Local
Density Approximation (LDA) 5 . In LDA, the exchange-correlation is approximated by
considering the inhomogeneous electron gas as being locally homogeneous, thus applying
the exchange-correlation hole of the homogeneous electron gas. In this way, the exchange-
correlation energy is found by integration with an energy density e LDA
XC :
Z
EeXC
LDA
[ρ] = ρ(r)e LDA
XC [ρ(r)] dr . (2.15)
The exchange-correlation hole and pair-correlation function for an inhomogeneous electron

gas is however non-local. To introduce some semi-local behavior and improve on LDA, the
Generalized Gradient Approximation (GGA) can be invoked 53 (p. 86) . This approximation
improves over LDA by taking into account the gradient of the density. GGA functionals
are constructed by augmenting the LDA energy density with an exchange-correlation
enhancement factor FXC [ρ]. Thereby the exchange-correlation energy becomes
Z
EeXC
GGA
[ρ] = ρ(r)e LDA
XC [ρ(r)]FXC [ρ(r), ∇ρ(r)] dr . (2.16)
Multiple GGA functionals have been developed throughout time, and one popular choice
is the Perdew-Burke-Ernzerhof (PBE) functional 54 . Another GGA functional derived
from PBE is the revised PBE functional of Hammer et al. 55 (RPBE), which is optimized
for chemisorption. Metallic systems are generally well described using PBE, whereas
chemisorption problems might prove more difficult. For example, PBE predicts the lattice
8
constant of Pd to be 56 3.99 Å, and RPBE gives 4.02 Å. Experiments give 56 a lattice
constant of 3.876 Å, and PBE comes closer to the experimental lattice constant. PBE
predicts the chemisorption energy of CO on Pd(111) at the fcc site to be 55 -1.94 eV,
whereas RPBE yields -1.65 eV. The corresponding experimental value is 55 -1.47 eV.
Similar improvements for RPBE over PBE in the chemisorption energy exist for oxygen 55 .
This indicates that PBE performs better for the solid, whereas RPBE predicts good
chemisorption energies. In microkinetic modeling the chemisorption energies are most
important as lattice-parameters do not enter the rate constant calculations explicitly.
Therefore, the RPBE functional is applied in Paper II and Paper III to model CO oxidation
over Pt.
The exchange-correlation energy of GGA functionals can be improved further, and
different strategies exists 53 (chap. 5) such as, meta-GGAs, van der Waals functionals, and
hybrid functionals where exact exchange energies are included.
The discussion concerning DFT is straightforwardly generalized to spin-polarized
systems. This is done by decomposing ρ into a spin-up density and a spin-down density
(ρ = ρ↑ + ρ↓ ), where one must solve the Kohn-Sham equations for each spin channel
separately.
2.2.2 Projector-Augmented Wave Method

When performing DFT calculations in practice, a certain basis set is needed to represent
the Kohn-Sham orbitals. In the present work, the applied basis functions are plane waves.
Plane waves are advantageous for periodic systems as they fulfill Bloch’s theorem 53 (p. 129) .
Another advantage is that increasing the number of plane wave components in the basis
results in monotonous convergence of the energy. Other choices of basis functions include
numerical grids and atomic orbitals 53 (chap. 8) .
The wavefunctions close to the nuclei oscillate rapidly due to orthogonality, and it
would require an unreasonable plane-wave basis to describe these oscillations. To avoid
describing the oscillations explicitly, the Projector-Augmented Wave (PAW) method 57
is used. In PAW, the core electrons are usually treated as chemically frozen as they do
not directly take part in the rehybridization, and only the valence electrons are treated
explicitly. The valence wavefunctions are treated by introducing a linear transformation T
between the true wavefunction |Ψi and a pseudo wavefunction |Ψi. e Real-space is divided
up into interstitial regions between the atomic cores and augmentation regions around the
cores. The pseudo wavefunction is chosen to be a smooth function inside augmentation
regions, which needs a smaller basis set. The method requires that |Ψ fi i = |Ψi i outside
the augmentation regions, and that they match in value and derivative at the boundary.
The transformation T only acts locally in each augmentation region, and it is the identity
transformation outside.
The transformation is performed using different basis functions than plane-waves. This
basis functions can be the all-electron partial waves |φi i, which are found by integration
of the radial part of the Schrödinger Equation for the isolated atom. Here i denotes all
angular momentum quantum numbers and an index labeling the atom in question. i also
contains an index n that describes different partial waves for identical angular momenta
and atomic sites. Each |φi i is assigned a pseudo partial wave |φei i that is equal to |φi i
9
outside of the augmentation region. |φei i forms a complete set inside. Using the partial
waves as basis functions inside the augmentation regions, the transformation takes the
form 57
X
T =I+ |φi i − |φei i he
pi | (2.17)
i
e
|Ψi = T |Ψi
Where I is the identity operator and he pi | is a projector function. There is one he

pi | for
each |φei i, which fulfill the conditions:
X
|φei i he
pi | = I, pi |φej i = δij .
he (2.18)
i
It should be noted that the projectors are local operators inside each augmentation region
that project out the relevant part of the wavefunction. Hence, inside an augmentation
region the wavefunction is the augmented part of the sum in (2.17), whereas outside the
e and no augmentation takes place.
projectors are orthogonal to |Ψi
2.2.3 Periodic Systems

The atomistic model systems in the present work were chosen to be periodic systems. The
wavefunction of an electron in a periodic potential (ψk ) obeys Bloch’s theorem, which
can be stated as 53 (p. 129)
ψk (r) = eik·r uk (r), (2.19)
where k is the wavenumber and uk is a function with the periodicity of the potential.
Using the periodic zone scheme 58 (p. 223) and working in the first Brillouin zone, the
wavefunction of an infinite solid is described taking only one unit-cell into account. In
principle, all possible wavenumbers must be included in the calculation, however in practice
this is solved by Brillouin zone sampling. In Brillouin zone sampling, one represents the
density, hence the energy, as a sum over special wavenumbers in the first Brillouin zone
called k-points: X
ρ(r) = ωk |ψk (r)|2 (2.20)
k∈BZ
Where ωk are weights that depend of the symmetry of the Brillouin zone. In practice,
symmetry is applied and only points in the irreducible wedge of the Brillouin zone are
considered. The density is given as a sum over each k-point as
X Nk
X (k) (k)
ρ(r) = ωk fi |φi (r)|2 , (2.21)
k∈BZ i=1
(k) (k)
where φi (r) are the eigenstates of the Kohn-Sham equations, fi are the occupation
numbers of the orbitals, and Nk is the number of occupied electronic states with wavenum-
ber k. It is necessary to solve a set of coupled Kohn-Sham equations, one for each k-point.
The number of k-points is a quantity, which must be increased until convergence in the
10
quantity of interest is achieved. Furthermore, it should be noted that a larger supercell
will require fewer k-points for convergence as a larger unit-cell will result in a smaller
Brillouin zone.
When calculating the density it can have practical relevance to work with partial
occupation numbers, which is referred to as smearing of the Fermi surface. For metallic
systems, this might be necessary as different k-points can make new bands enter or exit
the calculation during the self-consistent solution of the Kohn-Sham equations.
2.3 Local Minimum Energies

The previous sections have described how the energy of an atomistic system can be obtained
given the atomic positions. However, the precise positions that correspond to local energy
minima are generally not known. Energy minima are important to find, since a system
will spend a major part of its time near the minima. The local optimization procedure
used in this thesis is the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm 59 . The
BFGS algorithm is of the Quasi-Newton type where the atomic positions are updated
according to the Hessian matrix.
Calculating adsorption energies is done in a series of steps. First, the bulk unit-cell is
used to calculate the lattice constant that corresponds to the lowest energy, for example by
fitting an equation of state. Next a surface slab is constructed, where the number of layers
required to converge the surface- or adsorption-energy is determined. Then adsorbates
are placed over sites on the surface, and the system is relaxed locally. There can be
multiple local minima on a surface, and each possible minimum should be probed to find
the lowest energy. The calculations should also be converged with respect to simulation
cell size. Furthermore, a vacuum layer is introduced perpendicular to the surface to avoid
interactions between periodic repetitions of the system. For gas-phase moelcules, a large
cell is required to avoid interactions between periodic repetitions of the system.
Vibrations of gas-phase molecules and adsorbates can also be calculated with DFT.
This is often done by assuming a harmonic potential and finding the vibrational frequencies
using finite differences. The vibrations become important when analyzing transition states
as well as correcting the bare DFT energies for zero-point motion. A zero-point energy is
present for confined matter since it defines a ∆x in Heisenberg’s uncertainty principle.
Thereby a small ∆x implies a larger ∆p, which is responsible for the zero-point energy.
2.4 Energy Barriers with Nudged Elastic Band

The energy barrier for moving a system from the initial state to the transition state plays
a central role in describing chemical reactions. The transition state energy is needed
to calculate the rate constants. Nudged Elastic Band 60 (NEB) is a method to find the
Minimum Energy Path (MEP) between an initial state and final state, and NEB is also
used to identify the transition state.
To initialize a NEB calculation, the initial and final atomic positions are specified.
Next, the atomic positions are interpolated between the initial and final state. The
resulting atomic configurations are called images, and the distance between the images
are described by a reaction coordinate r. Structural optimization is done on each image,
11
which is connected to the other images by modifying the force component parallel to
the MEP. The total force on each image is decomposed in parallel and perpendicular
components as:
F = Fk + F⊥ = −kr + F⊥ (2.22)
where k is a spring constant. Hence, each image interacts with its neighbor images by a
spring-force parallel to the MEP, which prohibits full relaxation of each image down to a
local minimum (initial or final state). The perpendicular component ensures that each
image relaxes down to the MEP at the given reaction coordinate. In this manner, the
images are stretched along the MEP.
Typically, the transition state is not found without releasing the highest energy image
from all spring forces, and reversing the true parallel force component. As a result,
the image climbs up the MEP to the true transition state. This procedure is called
the Climbing Image method 60 , which can be switched on after a coarse initial NEB
calculation. Whether the highest energy point is a saddle point on the MEP can be
verified by vibrating the highest energy image. If the image has one imaginary frequency
in the direction of the reaction coordinate, the image is a saddle point on the energy
surface and can be considered a transition state.
2.5 Adsorbate-Adsorbate Interactions

The adsorbates on a catalyst surface interact either by direct interactions or by modification
of the surface electronic structure 61;62 . Such phenomena can have large effects on the
binding energies of certain species and consequently may be important for the kinetics. A
microkinetic model should take this into account by modifying the energies as a function of
the adsorbate configurations. Paper I treats adsorbate-adsorbate interactions as functions
of the adsorbate coverages, where the energies of the most abundant species (O, C, CO)
are perturbed linearly using calculations of the differential adsorption energies. This
turns out to be important for obtaining reasonable coverages. In Paper II, an exponential
relation is fitted for the CO-CO and O-O interactions, whereas the CO-O interactions are
modified linearly. It is important to use the differential adsorption energies instead of
the average adsorption energies, as adding another adsorbate will alter the energy of the
system based on the adsorbates that are already present.
Adsorbate-adsorbate interactions can also be fitted to specific atomic arrangements.
This is achieved by so-called cluster expansion Hamiltonians 63 (p. 94) where one calculates
the energies of certain geometric occupation-patterns and fits a Hamiltonian to the system.
The expansion is truncated at a certain number of nearest neighbors, typically including
two and three-particle interactions. In Paper III, the adsorbate-adsorbate interactions are
considered between the first nearest-neighbors, and two-particle interactions are accounted
for.
Figure 2.1 shows simulated coverages plotted as a function of temperature for complete
methane oxidation over Pd(100) using the microkinetic model of Paper I. The applied
model includes O-O, C-C, and OH-OH interactions. The figure shows coverages when the
interactions are turned off and are at full strength. Without any interactions, O has a
very high coverage and there is no C on the surface. At full interaction strength, there is
12
Figure 2.1: Simulated coverages without adsorbate-adsorbate interactions (left) and
with adsorbate-adsorbate interactions (right) during complete methane oxidation
over Pd(100). Pressures: 0.61 mbar methane and 3.06 mbar oxygen.
a coverage of both C and O. Thus, the adsorbate-adsorbate interactions are of utmost

importance for obtaining reasonable coverages and kinetics.
2.6 Modeling Entropies in Catalytic Reactions

The Gibbs free energies entering the rate constants are composed of an enthalpic con-
tribution and an entropic part. Thus, to compute rate constants, the entropy difference
between the initial and transition state must be modeled.
To model entropy, it is customary to work in the canonical ensemble and thereby use
the canonical partition function Z, which is a sum over all states, available to the system.
The entropy and partition function are related by the connection between statistical
mechanics and thermodynamics 64

∂ 1 ∂Z
S=− (−kB T ln Z)V,N = kB ln Z + kB T . (2.23)
∂T Z ∂T V,N
In this way, the entropy is uniquely defined by the partition function.

The partition function of an ideal-gas molecule can be split up into translations,
rotation, and vibrations. Assuming that these degrees of freedom are independent, the
partition function is a product of the individual contributions:
Z = Ztrans Zrot Zvib . (2.24)
The translational part of the partition function of the free-gas molecule is given by 2 (p.89)
3/2
2πM kB T
Ztrans = V , (2.25)
h2
where V is the volume of the gas, and M is the mass. V is often substituted for the
pressure using the ideal gas equation of state pV = kB T .
13
The rotational partition-function can be calculated in the rigid-rotor approxima-
tion 2 (p. 92) : 3/2 p
1 8π 2 kB T
Zrot = πIA IB IC (2.26)
σ h2
where IX is the moment of inertia of the molecule around the principal axis X, and σ is
the symmetry factor of the molecule.
Finally, the vibrational partition function in the Harmonic Approximation has the
following form 2 (p. 90)
~ωi
−
Y e 2kB T
Zvib = , (2.27)
~ωi
i −
1 − e kB T
where ωi is the vibrational frequency of mode number i. The Harmonic Approximation
assumes that the molecule is subject to a parabolic potential leading to the harmonic
oscillator energies ~ωi , which can be calculated using finite differences in the energies
obtained with DFT.
The partition function of adsorbates differs from gas-phase molecules. Paper II explores
the effect of adsorbate entropy in microkinetic modeling. In this context, the method of
Complete Potential Energy Sampling (CPES) is developed to estimate adsorbate entropy
with DFT. Four different adsorbate entropy approximations were compared in Paper II:
the Harmonic Approximation, the Hindered Translator, the Free Translator, and CPES.
These approximations differ in the description of the translational partition function
parallel to the surface. All remaining degrees of freedom are described as frustrated
vibrations.
In the Harmonic Approximation, the adsorbate fixed to one site and subject to a
harmonic potential. Thereby all modes are frustrated vibrations, and (2.27) determines the
full partition function. If low-energy modes are found (< 100 cm−1 ), it can be necessary
to truncate the value, as low-energy modes are computationally uncertain.
The Hindered Translator 65;66 allows for motion of the adsorbates, depending on their
diffusion barriers. The model is simple as it ignores surface topology by only treating
the fastest path of diffusion. In the Hindered Translator, the adsorbate is subject to a
sinusoidal potential, and the two surface-parallel directions are treated with the following
partition function 66 :

πrx rx + 1 2 rx
Nsites exp − I0
Tx Tx 2Tx 2
ZTHT = 2 exp . (2.28)
1 (2 + 16rx )Tx
1 − exp −
Tx
Tx = kB T /hνx is the ratio between the thermal energy and the vibrational energy, and
rx = Wx /hνx is the ratio between the diffusion barrier Wx and the vibrational energy.
Nsites is the number of surface sites, and I0 is the zero-order modified Bessel function of the
first kind, arising from the periodic potential. The last exponential in (2.28) is accounting
for zero-point energy corrections 66 . There can be some uncertainties in working with
14
the Hindered Translator. For example, it treats surface symmetry in a rough way by
assuming that the potential energy surface for diffusion obeys V (x, y) = V (x) + V (y),
which is questionable for most surfaces. The Hindered Translator is, however, expected
to be accurate for closed-shell molecules due to the presumably flat energy landscape for
diffusion.
The Free Translator is derived by assuming completely free diffusion of the adsorbate
over the surface. The translational partition function is given by that of a 2D free-gas:
Free 2πM kB T
Ztrans =A , (2.29)
h2
where A is the area that the adsorbate is free to translate over. The free translator defines
an upper bound to the translational entropy of any adsorbate as it corresponds to a
completely flat potential energy landscape.
CPES is developed in Paper II as an alternative method for calculating entropies. The
method uses DFT calculations to map out the potential energy as a function of adsorbate
position in the unit-cell. The translational partition function is evaluated by numerical
integration of a semi-classical partition function 64 :
ZZ −V (x, y)
2πM kB T
CPES
Ztrans = e kB T dx dy, (2.30)
h2
where V (x, y) is the potential energy of the molecule at position (x, y) in the unit-cell.
V (x, y) is calculated with DFT by locally optimizing the adsorbate in the direction
perpendicular to the surface, for various (x, y). For a completely flat potential energy
surface with V (x, y) = 0, (2.30) coincides with the Free Translator partition function
(2.29). The CPES method takes the surface symmetry into account as well as the detailed
potential energy landscape. In this sense, the CPES method is as accurate as the chosen
exchange-correlation potential used in the DFT calculations. Rotations can also be
modeled using (2.30) by mapping out a potential energy surface for rotation over some
angle θ, and evaluating the integral over θ. However, a complication is that the potential
energy for rotations and translations over a surface may be mutually dependent, which
can be solved by treating translations and rotations together using (2.30) and mapping
out V (x, y, θ).
The approximations discussed in this section are only valid in the low coverage limit
as they do not include adsorbate-adsorbate interactions. However, CPES can include the
effect of adsorbate-adsorbate interactions by reformulating (2.30) as a multi-dimensional
integral over multiple atoms by mapping out a potential V (x1 , y1 , x2 , y2 , · · · ).
2.7 Summary
This chapter has described how rate constants can be calculated by DFT. To describe
the Gibbs free energy properly, modeling of adsorbate-adsorbate interactions and entropy
calculations were also discussed. The uncertainties in calculating rate constants primar-
ily come from the uncertainty in the DFT calculations, the approximate treatment of
15
adsorbate-adsorbate interactions, and approximations made to the entropy. DFT calcula-
tions and adsorbate-adsorbate interactions give a large contribution to the uncertainty,
whereas entropy has a slightly lower impact.
16
Chapter 3
Simulating Reaction Kinetics
The previous chapter described how to calculate the rate constants of catalytic reactions,
which are determined by Gibbs free energies. The present chapter deals with simulating
and analyzing reaction kinetics, where the discussion begins with the chemical master
equation.
3.1 Chemical Master Equation

The kinetics of a chemical reaction can be viewed as a set of transitions of a system into
various states. Reactions occur on a high-dimensional free energy surface with several
maxima and minima. In principle, a reaction can be followed in space and time by solving
Newton’s equations of motions with molecular dynamics. However, molecular dynamics
is unfeasible for simulations of reaction kinetics, as the time-step is several orders of
magnitude smaller than the typical time of a chemical reaction. This is due to the fact
that molecular dynamics spends most time simulating vibrations.
The chemical master equation takes advantage of the time-scale separation between
vibrations and chemical reactions by dividing phase-space into regions corresponding to
different chemical states of the system. Consider two states of the system: α and β. The
chemical master equation connected to state α is given by 63 (p. 31) :
dPα X
= [Wαβ Pβ − Wβα Pα ] , (3.1)
dt
β
where Wαβ is the transition rate from β to α, and Pα is the probability for the system to
be in α. The sum in (3.1) runs over all states β, which can bring the system into state α.
In practice, the state is often defined by the coverages on the catalyst sites. Thereby, α
can be represented as
α = {A∗ , B∗ , ∗, ∗, A∗ , ∗, A∗ , B∗ , · · ·} (3.2)
Where A and B are chemical species, and ∗ denotes an empty site. The set is simply the
species occupying the sites of the system. The chemical master equation is the key to the
kinetics of the system as it describes the time-evolution for the probability of observing
the different states. A high transition rate into a particular state and a low transition
rate out of the state, implies that its probability increases in time. Thus, the kinetics will
follow the fastest set of transitions into a set of most probable states, which defines the
equilibrium.
dPα
At equilibrium, one interesting feature of the master equation is that = 0. In
dt
this case, Wαβ Pβ (t0 ) = Wβα Pα (t0 ), where t0 is a time at equilibrium. This is known as
17
the principle of detailed balance, which reflects that each step is in equilibrium with its
own inverse.
3.2 Mean-Field Approximation

The Mean-Field Approximation (MFA) yields an approximate solution to the chemical
master equation. To derive the MFA, consider a model surface reaction A + B → AB.
The average number of species A on the surface is given by 63 (p. 105) :
 
X dhN A i X X
hNA i = Pα Nα(A) ⇒ =  [Wαβ Pβ − Wβα Pα ] Nα(A) , (3.3)
α
dt α β
(A)
where Nα is the number of species A in state α, and Pα is the probability of the system
to be in state α. The time-derivative was inserted from the master equation (3.1). The
indices of summation can be rearranged to yield 63 (p. 111) :
dhNA i X X h (A) (A)

i
= Nα − Nβ Pβ Wαβ (3.4)
dt α β
For a reaction A + B → AB, let α correspond to a state with AB, and β corresponds
(A) (A)
to a state with A + B. Then it must hold that Nα − Nβ = −1. Additionally, given
α, then for each β a well-defined number of terms in the sum have non-zero Wαβ . This
(AB)
number is equal to the number of adjacent A+B pairs on the surface (Nβ ). With the
transition rate assumed equal for each β, one obtains 63 (p. 112) :
dhNA i X (AB)
= −Pβ W Nβ = −W hN (AB) i (3.5)
dt
β
The MFA assumes that the adsorbates are randomly distributed. In this case, the number
of pairs in state β becomes 63 (p. 112) :
(B)
(AB) (A) Nβ
Nβ = ZNβ , (3.6)
S−1
where Z is the number of nearest neighbors, S is the number of surface sites, and
(B)
Nβ /S − 1 is the probability for a site to be occupied by B. Inserting (3.6) into (3.5)
gives:
(B)
dhNA i X (A) Nβ WZ
= −Pβ W ZNβ =− hN (A) N (B) i. (3.7)
dt S−1 S−1
β
The mean of a product can be written as the product of the means augmented with the
fluctuations from the individual means:
D E D ED E Dh ih iE
N (A) N (B) = N (A) N (B) + N (A) − hN (A) i N (B) − hN (B) i . (3.8)
18
In the limit of an infinite number of sites, S − 1 ≈ S, and fluctuations can be neglected
which gives:
dhNA i W Z D (A) E D (B) E
=− N N (3.9)
dt S
Dividing by the number of sites, we find the mean-field equation for the time-evolution of
A:
dθA
= −kθA θB , k = W Z, (3.10)
dt
where θi is the coverage of species i on the surface, and k is identified as the rate constant
of the reaction. Similar equations can be derived for θB and the fraction of empty sites.
The two main assumptions of the MFA are that the surface is infinite, and that adsorbates
are randomly distributed. θA corresponds to the probability for observing species A on a
surface site. The random distribution implies that configurational entropy is accounted for
in the mean-field equations. The assumption that the surface is infinite (S 1) breaks
down for nanoparticles, where fluctuations cannot be neglected. Moreover, the MFA
breaks down for strong adsorbate-adsorbate interactions, and the situation is particularly
severe for attractive interactions that lead to pairing of adsorbates.
The advantages of the mean-field picture is that it focuses on the average behavior
of the system by describing the coverages. This means that results often are simple to
analyze, and that the computational cost of mean-field modeling is considerably smaller
than Monte Carlo simulations as the mean-field equations often can be solved on a desktop
computer within some minutes. The disadvantage of the MFA is the limitation to infinite
surfaces. Furthermore, when working with coverages spacial information is lost, which
can be central to understand a reaction mechanism.
In practice, a mean-field model is constructed by having a list of coverages (θ) for
species that may cover the surface, and a set of reactions with rates R(θ). From this
information, a system of coupled first-order ordinary differential equations is formulated:
dθi X
= cij Rj (θ) , (3.11)
dt j
where cij is an integer that describes how many i-species that are produced in reaction j.
Rj is the rate of the elementary reaction j, which is equal to the rate constant times the
relevant coverages. As an example, consider the following schematic reaction:
A2 (g) + 2∗ ↔ 2A∗ (a)

∗
B(g) + ∗ ↔ B (b)
A∗ + B ∗ → AB(g) + 2∗ (c)
The time-evolution of the coverage of A is given by:

dθA
= 2ka+ θ∗2 − 2ka− θA
2
− kc+ θA θB (3.12)
dt
where ki± is the rate constants for reaction i, and * indicates an empty site. The first term
comes from the forward-reaction of (a), which requires two sites for dissociation such that
19
θ∗ enters squared. The factor two is is present as two A are formed upon adsorption. The
second term stems from the backward reaction of (a). It is negative as A is consumed by
the elementary step. The third term arises from reaction (c), which requires one A and B.
For species B the coverage equation is:
dθB
= kb+ θ∗ − kb− θB − kc+ θA θB . (3.13)
dt
The resulting system of differential equations is often solved numerically, and in the present
work this is done using Python and SciPy 67 . Before integration, the initial conditions are
specified, which is important as different initial conditions can lead to different steady
states. The steady state is reached when all coverages are constant in time. A common
situation is that the rate constants vary by several order of magnitudes, which makes
the equations stiff. Therefore, an integration method that handles stiff problems must
be chosen, and in the present work the Backward Differentiation Formula 68 method was
found to be most effective.
3.3 Kinetic Monte Carlo

Kinetic Monte Carlo 69 (KMC) is a stochastic method to solve the chemical master
equation (3.1). KMC must be applied in situations when the adsorbate distribution is not
random or the system is finite-sized. A nanoparticles is an example of a finite system with
specific adsorbate patterns. In KMC, the system starts in a specific state, and evolves in
time to different states governed by random number generation. KMC solves the master
equation by simulating the transitions between states as a function of time. This is done
by performing a large number of reaction events and their time of observation.
The KMC method can be implemented using various algorithms, for which two
popular choices are the Variable Step Size Method (VSSM) and the First-reaction method
(FRM) 63 (Chap. 3) . VSSM and FRM are two implementations of the same method, which
is also known as the n-fold way 70 . In the present work, the shortest computational time is
obtained using the FRM algorithm. Before performing a simulation, the types of possible
reaction events are defined. When this is done, FRM can be performed by a series of
steps, which are summarized in Figure 3.1.
Step 1: The simulation is initialized by setting time t = 0, defining when the simulation
should end (tend ), choosing the reaction conditions, and initializing the site occupations.
Here an event-list is initialized to keep track of possible reaction events, their time of
occurrence, and the site where they proceed. The event-list is updated at each simulation
step.
Step 2: For all possible events, the times of occurrence are calculated, and the event-list
is populated. The time of occurrence tβα for taking the system from state α to state β is
calculated according to 63 (p. 53)
1
tβα = t − lnu (3.14)
Wβα
Where t is the current simulation time, Wβα is the rate constant of the event, and u is a
random uniform deviate on the interval ]0, 1].
20
Step 3: Determines whether the simulation has ended by checking if t ≥ tend .
Step 4: Performs the chronologically next reaction in the event-list. If the next reaction
is impossible, it is discarded and the next possible event is performed.
Step 5: The event-list is updated by adding new events made possible by the last
simulation step. This only done in a neighborhood around the site where the last event
happened, since avoiding updating globally saves a significant amount of computational
resources for large systems.
Step 6: Defines the end of the simulation, where quantities of interest are saved. To
lower the memory consumption, saving data and cleaning different lists should also be
done during the simulation.
1. Initialize Simulation : Set simulation time (t = 0), simula-

tion end-time tend , temperature, pressure, and initial site-
occupations.
2. Populate event-list and generate occurrence times tβα .
3. Is t ≥ tend ? If yes skip to 6
4. Perform next possible event in queue, and set t to tβα of the

chosen event.
5. Update list of possible events. Jump to 3
6. Finalize simulation: Save results.
Figure 3.1: First Reaction Method Algorithm.
The sites and their mutual connectivity must be represented somehow in kinetic Monte
Carlo simulations. In Paper III, neighbor lists are used to define the connection between
sites. The neighbor lists for all sites constitute a global connectivity pattern. Alternatively,
a lattice-based algorithm can be implemented, where the connected sites are adjacent on
a lattice. The implemented neighbor-list approach has an advantage over a lattice-based
approach, as incommensurate lattices and dynamical catalyst changes can be accounted
for in a simpler manner.
The time-step in the simulations is governed by the rate of the fastest reactions since
these events are executed with the highest probability. Figure 3.2 (left) shows time as a
function of simulation step during a short Monte Carlo simulation. The time increment is
small in the beginning of the simulation, and becomes larger as the simulation proceeds,
which is owing to decreasing fastest rates. Figure 3.2 (right) shows a histogram over
the time step. The shortest time-steps are mostly originating from the beginning of
the simulation, where the reaction is far from equilibrium. The histogram resembles
an exponential distribution, which is reasonable as the events essentially are Poisson
processes.
A considerable challenge in KMC for chemical reactions is that the fastest reactions
21
Figure 3.2: Time in kinetic Monte Carlo simulations of CO oxidation on a 5.2 nm
Pt nanoparticle. Left: Time as a function of simulation step. Right: A histogram
over the time step. Reaction conditions: 600 K, 20 mbar CO, and 10 mbar O2 .
The CO diffusion barrier is increased by 0.5 eV.
determine the time-step. For example, diffusion events can be several orders of magnitudes
faster than chemical reactions, which can make KMC simulations computationally unfea-
sible. To overcome this problem, one can raise the barriers of the fast events. However,
this must be done with care as the fast events need to remain equilibrated. This simple
method has previously been justified theoretically 63 (pp. 162-163) . A sophisticated solution
is to identify events that are in equilibrium during the simulation and increase these rate
constants during the simulation 71–74 . This can be useful when there are many possible
reactions and no prior knowledge of the reactions is available.
KMC enables investigations of more realistic systems than the MFA since it does not
assume complete randomness, nor an infinite number of sites. This also implies that
small nanoparticles are possible to simulate. Disadvantages of KMC simulations are that
the results are more complex to analyze, and that the computational cost is significantly
higher than in the MFA.
3.4 Analyzing Reaction Kinetics
Mircokinetic models give access to various information, which can be compared with
experimental measurements. This section discusses the information that is obtained from
a microkinetic simulation, and how to analyze the results in relation to experiments.
3.4.1 Reaction Mechanisms

One advantage of doing a microkinetic analysis is that it provides information about the
most important reaction mechanism; provided that the relevant reaction pathways are
included in the simulations. The active mechanisms are found using the net-rates, which
22
can be written in a mean-field picture as:
Y Y
Wi = ki+ θl − ki− θm . (3.15)
l m
Where l is the set of reactant species for the elementary reaction i, and m is the set
of products. In kinetic Monte Carlo, one can simply observe which reactions that are
executed most frequently. The dominant reaction mechanism will be the set of reactions
with the highest net-rates that complete a catalytic cycle.
3.4.2 Turnover Frequencies and Coverages

The Turnover Frequcency (TOF) is a measure of the number of products formed per active
site and time-unit. In a microkinetic model, the TOF is simply the net rate of product
formation. Quantitative comparison between a simulated TOF and an experimental result
is complicated, since there are many uncertainties and approximations involved in the
modeling. Hence, it is common that the simulated TOF deviates with a few orders of
magnitudes from experiments. However, the trends and variations across various reaction
conditions often resemble experimental behavior.
The coverages are a simple measure to analyze the reaction, and the most abundant
species can be compared to spectroscopic experiments that show the adsorbates present
on the surface. Furthermore, knowing the most abundant species can assist in improving
the kinetic model, for example by improving adsorbate-adsorbate interactions or entropies
for the most abundant species.
3.4.3 Reaction Orders and Apparent Activation Energies

The reaction order in a given pressure and the apparent activation energy are two
macroscopic quantities that can be obtained from a microkinetic model and compared to
experimental values.
The reaction orders are phenomenological quantities that are defined assuming a
power-law in the catalytic rate 2 (pp. 26-27) :
r = Kpnx x , (3.16)
where K is a constant with respect to pressure, px is the partial pressure in species x,

and nx is the reaction order in species x. To obtain nx from simulations, the TOF is
fitted to this power-law.
The apparent activation energy Eapp is obtained by assuming an Arrhenius relation in
the turnover frequency 2 (pp. 36-37) :

−Eapp Y nx
r = A exp px , (3.17)
kB T x
where A is a the pre-exponential factor, x is a gas-species, px is the pressure, and nx is

the reaction order. In practice, Eapp is calculated numerically by evaluating the derivative
of r versus 1/T , over a narrow temperature range.
23
Direct comparison with experimental values is not always possible since the structure
in experiments is less well-defined than in theoretical models. For the methane oxidation
reaction in Paper I, the reaction orders agree well with experiments 75;76 . However, the
apparent activation energy did not agree well with experimental values 75–77 . This is likely
due to surface-oxide formation during the experiments, which was not included in the
model. In this manner, the discrepancies between experiments and a validated microkinetic
model can provide important insights about the catalyst state during experiments.
3.4.4 Degree of Rate Control

The degree of rate control 78;79 of an elementary step (χi ) can be thought of as the
sensitivity of the catalytic rate r to the individual rate constants. χi is defined by:
!
kif ∂r
χi = , (3.18)
r ∂kif Ki
where kif is the forward rate constant of step i, and the derivative is taken with a fixed
equilibrium constant Ki . Thus, χ provides information on the slow steps that limit the
catalytic rate. In practice, the derivative is evaluated with finite differences by increasing
the rate constant of the step a small percentage, e.g. 1%. χ is important in microkinetic
modeling as it determines which reaction steps that are affecting the rate, and χ can
therefore provide useful information for catalyst design. Moreover, it is important to
model these steps with a good accuracy.
Paper IV connects the microscopic kinetics to the reaction orders and apparent
activation energies using the degree of rate control. In this way, these two macroscopic
quantities are linked to the atomistic picture. The reaction orders are shown to have a
simple relation with χ:
X ∂lnWi
nx = χi . (3.19)
i
∂lnpx
The reaction order is a weighted sum over the individual rates Wi , derived with respect to
the pressure of interest. The reaction order in the presence of a catalyst has previously been
described as a phenomenological quantity 80 , however, the connection between reaction
orders and the elementary reactions are shown in (3.19).
The relation between χ and the apparent activation energy takes the form
X X
∂Si kB T ∂Pj
Eapp = χi Ei + kB T + kB T 2 + χj Ej − + kB T 2 (3.20)
∂T j∈trans
2 ∂T
i∈vib
X ∂nx
− kB T 2 ln px
x
∂T
Where Ei is the energy barrier of elementary step i, Si is the entropic barrier, and Pj
is the sticking coefficient related to step j. nx is the reaction order in species x, and px
is the partial pressure. The first sum runs over the steps that have a vibration as the
reaction coordinate, and the second sum is over reactions with a translational reaction
24
coordinate. The third sum accounts for the pressure dependence of the total rate. It is
interesting to note that the zeroth order contribution in temperature is the largest and
arises from the elementary energy barriers. The first-order temperature terms arise from
a lost degree of freedom in the reaction coordinate, and the second order terms stem from
entropic losses and pressure dependencies.
The derived expressions provide a microscopic understanding of the reaction orders
and apparent activation energies. Furthermore, the expressions might be used to evaluate
proposed rate-determining steps and reaction mechanisms, for example by combining
DFT calculations and experimental measurements.
3.5 Summary
The present chapter has discussed how the chemical master equation describes the reaction
kinetics. The mean-field approximation was derived as an approximate solution to the
master equation, and Kinetic Monte Carlo was presented as a stochastic method to solve
the master equation. Finally, analysis of results obtained from microkinetic models was
discussed with focus on comparing with experiments.
25
26
Chapter 4
Kinetic Models for Methane and CO

Oxidation
4.1 Constructing a Kinetic Model
Figure 4.1: Illustration of workflow for constructing a first-principles microkinetic

model.
The microkinetic models in this thesis are constructed using a workflow illustrated
in Figure 4.1. Formulation and refinement of a model is an iterative procedure, where
the cycle in Figure 4.1 is completed multiple times before settling on a final version of
the kinetic model. The workflow begins with identifying an atomistic model system of
the catalyst. This is often a crucial and difficult step as technical catalysts are complex
systems, which are ill-defined with respect to structure and composition. The next step
is to identify possible intermediates. These intermediates are connected by multiple
elementary reaction steps, which together constitute the catalytic cycle. Taking complete
methane oxidation to CO2 over Pd surfaces as an example, the atomistic model system
can be Pd(100) or Pd(111). Plausible intermediates are hydrocarbons, oxygen, carbon,
and CO. Thereby, a starting point for a possible reaction mechanism could be:
CH4 (g) + 5∗ → C∗ + 4H∗ (4.1)

O2 (g) + 2∗ → 2O∗ (4.2)
∗ ∗
C + 2O → CO2 (g) + 3∗, (4.3)
27
These steps are not elementary, and should be reformulated as elementary steps, possibly
by including alternative pathways. After the intermediates and elementary reactions have
been identified, rate constants are calculated for the considered elementary steps and the
kinetic equations are solved. Thereafter, the computational results are compared to exper-
imental data by simulating measurable quantities, such as turnover frequencies, reaction
orders, and apparent activation energies. The cycle is repeated until the model has the
desired level of detail and agrees with experimental data. A validated microkinetic model
opens up the possibility to use the simulations for predictions beyond the experimentally
investigated conditions.
4.2 Complete Methane Oxidation over Pd(100) and Pd(111)

In Paper I, the mean-field approximation is applied to model complete methane oxidation
over Pd(100) and Pd(111), which is a complex reaction with several possible intermediates
and mechanisms. The reaction is important to understand as methane is a potent
greenhouse gas, which has a larger global warming potential than CO2 81;82 . Moreover,
complete methane oxidation to CO2 is becoming an increasingly important reaction for
the automotive industry as the global natural-gas vehicle fleet is growing 83 . The overall
reaction is:
CH4 (g) + 2O2 (g) → CO2 (g) + 2H2 O(g). (4.4)
Pt and Pd are know to be active catalysts, where Pt is superior in net reducing conditions
and Pd is advantageous in net oxidizing environments 84 . Pd oxidizes readily at most
reaction conditions 82;85–87 and the high activity has previously been attributed to the
oxidized phase 85;88 , whereas other reports assigned the high activity to the metallic
state 89–91 . The metallic phase is difficult to stabilize during experiments, which makes
a first-principles microkinetic model a useful tool. Furthermore, a microkinetic model
reveals active reaction mechanisms and kinetic bottlenecks that can help improve the
reaction.
The kinetic model developed in this work extends a previous model of the reaction over
metallic Pd 36 , by taking into account a large number of reaction barriers and multiple
reaction mechanisms. One purpose of constructing a refined model is to compare with a
detailed model of the reaction over PdO(101) 35 . By comparing models of metallic Pd
and PdO, differences are revealed between the two phases, which is important in order to
determine the most active phase.
The simulations are performed in slight oxygen excess with pressures of 0.61 mbar
methane and 3.06 mbar oxygen. The model includes six gas-phase species, 16 adsorbed
species, and 32 elementary reactions, which enables multiple reaction pathways. The
included gas-phase species are: CH4 , O2 , CO2 , H2 O, CO, and H2 . Included adsorbates
are: CH3 , CH2 , CH, C, H, O2 , O, CO, OH, H2 O, CH2 OH, CH2 O, CHO, COH, OCOH,
and CHOH. The reaction is found to proceed via different mechanisms depending on the
reaction conditions. At high temperatures, methane is dehydrogenated on the surface one
hydrogen at a time to C + 4H. Thereafter, C is oxidized to CO and further to CO2 by
O. Simultaneously, H is oxidized to H2 O, mainly via OH + OH → H2 O + O. At lower
temperatures, instead C is oxidized by OH, and the mechanism involves steps such as
C + OH → COH and COH + O → OCOH.
28
Figure 4.2: Simulated kinetics of methane oxidation at 0.61 mbar methane and
3.06 mbar oxygen. (a) Turnover frequencies versus temperature. (b) Apparent
activation energies versus temperature. (c) Turnover frequency versus pressures
on Pd(100) at 598 K. (d) Degree of rate control versus temperature on Pd(100).
The coverages and TOFs are central quantities for understanding the kinetics. The
simulations predict that the most abundant adsorbates on the surfaces are O, C and CO.
O is present at all investigated temperatures on both Pd(100) and Pd(111), whereas C
and CO only are present at low temperatures. The simulated TOFs are shown in Figure
4.2 (a). The TOFs vary seven orders of magnitudes when increasing the temperature
from 400 to 1000 K, and the Pd(100) surface is found to be more active than Pd(111) at
all investigated conditions.
The apparent activation energy (Eapp ) reflects the Gibbs free energy barriers of the
rate controlling steps, as was derived in Paper IV. The simulated Eapp is shown in Figure
4.2 (b). Eapp varies between 1.2 and 0.85 eV for Pd(100) in the temperature range
400-1000 K. The variation is due to changing reaction mechanisms and rate-controlling
steps. An analytical expression is derived for Eapp in the limit of high temperatures,
which shows that Eapp is determined by the methane adsorption barrier and the oxygen
29
adsorption/desorption equilibrium as:
p
f pO2 KO2 ∆EO2
Eapp ≈ ECH − p , (4.5)
4
1+ pO2 KO2
f
where ECH 4
is the energy barrier for methane dissociative adsorption, pO2 is the oxygen
pressure, ∆EO2 is the dissociative adsorption energy of O2 , and KO2 is the equilibrium
constant for oxygen adsorption.
The reaction order reflects how the TOF responds to a change in pressure. Figure
4.2 (c) shows the logarithm of the TOF as a function of methane, oxygen, and water
pressure on Pd(100). The results for Pd(111) are similar. The analysis shows that a
higher methane pressure increases the TOF. Conversely, oxygen is inhibiting the rate by
blocking sites for methane dissociation, and water can promote the reaction by increasing
the rate of C + OH → COH, which is a rate-controlling step.
The degree of rate control analysis is a way to explain the kinetics in a compact manner.
The simulated degree of rate controls are shown in Figure 4.2 (d) for the main steps. The
main rate-controlling steps at low temperatures are C + OH → COH and C + O → CO.
At higher temperatures, dissociative methane adsorption eventually becomes the rate-
determining step. Oyxgen adsorption is slightly inhibiting over the entire temperature
range. These observations are reflected in the coverages as the most abundant species are
connected to the rate-controlling steps. The presence of C on the surface suggests that
oxidizing C has a finite rate-control, in agreement with the analysis. The reaction orders
are also explained by the degree of rate control, as is derived in Paper IV. The reaction
order in methane follows the degree of rate control for methane dissociation closely (Paper
I and Paper IV). The slightly negative oxygen order shows that oxygen adsorption is an
inhibition step, and thus has a negative rate control. The positive order in the water
pressure is readily understood by noting that increasing the OH concentration speeds up a
rate-controlling step: C + OH → COH. The apparent activation energy is also connected
to the degree of rate control (Paper IV). The high-temperature analytical expression
for the apparent activation energy (4.5), reflects that methane adsorption is the most
rate-controlling step, and oxygen adsorption is a slightly inhibiting step.
The model in Paper I for metallic Pd can be compared with methane oxidation over
PdO(101), which has been modeled in detail previously 35 . The reaction orders in the
methane pressure are high and positive for both PdO 35;75–77 and metallic Pd surfaces.
However, the PdO catalyst has qualitatively different reaction orders in the water and
oxygen pressures. Water is known to be poisoning the PdO catalyst 35;75–77;92 due to
site-blocking effects, whereas water can promote the reaction on metallic Pd. Oxygen
can promote the reaction on PdO at some conditions 35;76;77 , whereas it inhibits the
reaction on metallic Pd. These clear qualitative differences between Pd and PdO can aid
characterization of the active phase using only the reaction kinetics.
The aforementioned results apply to extended surfaces, however, technical catalysts
are nanostructured materials. The presence of low-coordinated sites on nanoparticles will
likely promote this reaction 8 as methane adsorption is the rate-determining step. Thus,
nanoparticles are expected to be more active for the reaction than extended surfaces,
given that other effects are not hindering the rate. The benefit of having a multitude of
sites in complete methane oxidation on nanoparticles remains to be investigated.
30
Figure 4.3: Schematic illustration of the potential energy landscape for translation
in (a) the Harmonic Approximation, (b) Hindered Translator, (c) Free Translator,
and (d) CPES.
4.3 Entropic Effects in CO Oxidation on Pt(111)

In Paper II, different models for adsorbate entropy are investigated. Previous first-
principles microkinetic models have typically described adsorbate-entropies in the Har-
monic Approximation 31;66;93 . Alternatively, weakly bound adsorbates have been treated
as losing one or all translational degrees of freedom upon adsorption 30;94–96 . Therefore,
detailed investigations of entropy in the context of microkinetic modeling is benefi-
cial. CO oxidation over Pt(111) is used as an archetype reaction. This reaction is
very well-studied 30;97–105 , and this enables a detailed comparison to experiments. The
reaction is modeled in the mean-field approximation using a Langmuir-Hinselwood mech-
anism 30;97;100–102;105;106 :
CO(g) + ∗ ↔ CO∗
O2 (g) + 2∗ ↔ 2O∗ (4.6)
CO∗ + O∗ → CO2 (g) + 2∗.
The entropies of CO and O are found to be qualitatively different. CO translates on a

flat potential energy surface with low diffusion barriers, which results in a high entropy.
O instead clearly prefers fcc binding and has high diffusion barriers, which results in a
low entropy. Typically, a mean-field model includes only a generic and coarse-grained
site, and it is blind to the detailed nature of the site. The coarse-grained site could entail
fcc, hcp, ontop, and bridge sites. When the adsorption energies on these sites are similar,
the translational entropy can be high and should be modeled accordingly.
Four different models for adsorbate entropy are evaluated: The Harmonic Approxima-
tion (HA), the Hindered Translator (HT), the Free Translator (FT), and the method of
Complete Potential Energy Sampling (CPES). The four models are described in detail in
section 2.6. The four methods differ in how the potential energy landscape for translation
31
Figure 4.4: Calculated adsorbate entropy versus temperature of CO (left) and O
(right) for one adsorbate in one Pt(111) fcc unit-cell.
is modeled. Figure 4.3 shows a schematic illustration of the potentials in the four models.
The Harmonic Approximation has a quadratic potential along the surface, the Hindered
Translator is sinusoidal, the Free Translator is unbound, and the CPES method uses the
potential energy landscape calculated by DFT.
The entropies predicted by the four approximations are given in Figure 4.4 for one
adsorbate in a single unit-cell. For both CO and O, the Harmonic Approximation predicts
the lowest entropies, as it fixates the adsorbate to one sitei . The Free Translator predicts
the highest entropies since it assumes barrierless adsorbate-diffusion. The Hindered
Translator and CPES predict entropies in between these two extremes. The Harmonic
Approximation gives an entropy for CO that is 0.4 meV/K lower than CPES. At 500 K,
this corresponds to a difference in adsorption energy of 0.2 eV, which can be considered a
significant difference.
The choice of entropy model affects the kinetics significantly, and the predicted TOFs
vary between the four models. The largest TOF is predicted in the Harmonic Approxima-
tion, where the adsorbate entropy is lowest, and the lowest TOF is obtained using the
Free Translator that is an upper bound to the translational entropy. The coverages are
also affected by the entropy. The lowest coverages are obtained in the Harmonic Approxi-
mation, whereas the highest coverages are found using the Free Translator. Additionally,
entropy affects the light-off temperatures, i.e. the temperature for where the turnover
frequency reaches half maximum. Figure 4.5 shows a comparison between light-off tem-
peratures obtained by the developed kinetic model and various experiments 97;99;102;105 .
The Harmonic Approximation and Free Translator under and overestimate the light-off
temperature, respectively. The light-off temperatures predicted by the Hindered Transla-
tor and CPES are in good agreement with the experimental references over the entire
pressure-range. Hence, modeling the entropy in detail moves the results significantly
closer to experimentally obtained values. The effect of entropy on the kinetics is readily
understood by noting that a high entropy implies a stronger adsorbate-surface bond. This
i For example either fcc, hcp, bridge, or ontop
32
Figure 4.5: Simulated and experimental (Exp) light-off temperatures as a function
of pressure. The experimental references are: Campbell et al. 97 (10−7 mbar),
Vogel et al. 102 (10−6 mbar), Nakao et al. 99 (10−2 mbar), and Calderón et al. 105
(10−1 -100 mbar).
effect is enhanced with increasing temperatures as ∆G = ∆E − T ∆S.

The overall conclusions of this study are that the treatment of adsorbate entropy
affects the results of the kinetic simulations, and mobile molecules should not be modeled
in the Harmonic Approximation.
4.4 Catalysis on Nanoparticles

Commercial catalysts often consist of metal nanoparticles, which implies that reactions
should be modeled using finite-sized systems. Therefore, to understand catalytic behavior
it is important to model nanoparticles. Nanoparticles are different from extended surfaces
as nanoparticles have different sites with various reactivities. As a consequence, the
kinetics is complex, and the mean-field approximation breaks down. Instead kinetic Monte
Carlo simulations should be applied.
When modeling nanoparticles, a large number of energies is required to represent the
reaction energy landscapes. Moreover, the reaction energies on extended surfaces and
nanoparticles can differ significantly 8;9;11;107 . Kinetic Monte Carlo models have previously
been formulated for nanoparticles using simplified energy landscapes 38;40 and schematic
reactions 49;50 . However, few studies have simulated a specific reaction using a detailed
reaction energy landscape.
In paper III, CO oxidation is modeled on Pt nanoparticles by developing a KMC
algorithm. The paper outlines a method to simulate nanoparticles by approximating
the reaction energy landscape. The energy landscape is approximated using generalized
coordination numbers, which lowers the computational cost significantly. Recently, the
33
Figure 4.6: Left: Generalized coordination numbers of ontop sites on a 2.8 nm
nanoparticle. Right: Activity of CO oxidation versus pressure and generalized
coordination number on a 2.8 nm Pt nanoparticle. The temperature is 850 K, and
the pressure ratio pCO /pO2 = 2.
generalized coordination number CN was shown to be a good descriptor for adsorption

energies on nanoparticles and surfaces 19–21 . CN is an extension of the conventional
coordination number, which accounts for the coordination numbers of the nearest neighbors
as:
X CNi
CN = , (4.7)
CNmax
i∈{N N }
where the sum runs over the nearest neighbors to the site in question, CNi is the
conventional coordination number of the nearest neighbor i, and CNmax is the maximally
possible coordination number of neighbors in the bulk. In a fcc crystal, the ontop site
has CNmax = 12, a bridge has CNmax = 18, a three-fold hollow has CNmax = 22, and
for a four-fold hollow site CNmax = 26. Figure 4.6 (left) shows CN for ontop sites on a
nanoparticle, where each color corresponds to one CN and hence one adsorption energy.
Furthermore, in Paper III the reaction energy barriers are obtained by a Brøndsted-Evans-
Polanyi (BEP) relation, which relates the adsorption energies to the transition state.
Using BEP relations, coverage dependencies in the reaction barriers are accounted for
automatically. In this manner, CN fully describes the reaction energy landscape with
limited computational effort.
CO oxidation over Pt is taken as a model reaction to study the effects of modeling
nanoparticles on the reaction kinetics. By comparing simulations on nanoparticles with
the extended Pt(111) surface, conclusions can be drawn concerning the materials gap.
The model predicts that nanoparticles have a higher TOF than the Pt(111) surface. One
reason is that the rate is not limited by adsorption as opposed to Pt(111). This is owing to
the low-coordinated sites on the nanoparticle that can contribute to the reaction at high
temperatures. Thus, the active site varies with reaction conditions, and as a consequence
the TOF scales with particle size. The most active sites at 850 K are illustrated in Figure
4.6 (right). Here, the catalytic activity is divided by the highest activity at each pressure,
34
Figure 4.7: Simulated CO coverage in CO oxidation over a 5.2 nm Pt nanoparticle.
Left: Temperature of 600 K. Right: Temperature of 1200 K. Pressures: 20 mbar
CO and 10 mbar O2 . The CO diffusion barrier is increased by 0.5 eV.
and it is normalized to the number of sites with the given CN. At low pressures, edges
and corners (CN< 5.66) are active, which implies that the TOF declines with increasing
particle size. At high pressures, the edges and facets (CN > 4.25) are most active, and
the TOF grows with particle size. Similarly, temperature affects the most active site. At
high temperatures, the corners and edges are most active, whereas low temperatures make
the edges and facets most active. These results demonstrate that nanoparticles cannot be
represented fully by extended surfaces. Furthermore, the results clearly indicate that the
mean-field approximation breaks down on nanoparticles as a mean-field model cannot
account for the site-distribution, which is important to capture the kinetics.
Paper III shows the average kinetic behavior, which is most relevant in studying
catalysts that operate over long time scales. However, the time evolution of a simulation
is also instructive to consider. Figure 4.7 shows the CO coverage on a nanoparticle
during a short simulation. At 600 K (left), the facets contain free sites, whereas the
edges and corners are completely filled. At 1200 K (right), the coverages are lower on
the edges, whereas the corners and facets are empty. One interesting observation is
that there are large fluctuations in the coverages, which are most pronounced for the
low-coordinated sites as they are less abundant. The magnitude of the fluctuations
increases with temperature and becomes considerable for smaller systems. A mean-field
model assumes that such fluctuations can be neglected, which is clearly not valid on
nanoparticles.
To understand the relation between structure and catalytic activity, it is instructive
to sequentially increase the complexity of the model system. An example is to start
by understanding extended surfaces, then stepped surfaces, nanoparticles, supported
nanoparticles, etc. However, the level of detail can reach a point where the results are too
complex and unclear to comprehend. Moreover, as DFT does not yield chemical accuracy,
at some point the results might be inaccurate. The present work has taken a step from
extended surfaces to truncated octahedral nanoparticles by using scaling relations derived
from DFT. This approach might seem rough at a first glance. However, calculating all
35
energies explicitly is less meaningful as the octahedral nanoparticle is to be considered a
model system, and due to the limited accuracy of DFT.
From the developed model, structure-activity relationships can be investigated, which
can be useful in catalyst design. The results indicate that a search for one single active site
might be a too simplistic approach on nanoparticles as various sites likely collaborate to
enable the reaction to follow a path of least resistance. With the developed method, kinetic
modeling of actual reactions can be performed routinely with limited computational effort,
and doing this will be a significant step towards bridging the materials gap.
36
Chapter 5
Conclusion & Outlook
The presented work has aimed to investigate and develop the methodology for first-
principles microkinetic modeling of catalytic reactions. This was realized in a series of
papers that addressed various topics.
Complete methane oxidation over extended palladium surfaces was modeled in detail,
which revealed the rate-controlling steps and detailed reaction mechanisms. The mecha-
nisms on metallic Pd were found to be fundamentally different from those on PdO, which
leads to interesting variations between the kinetics of the two phases.
The effect of adsorbate entropy was studied by developing a mean-field model of CO
oxidation over Pt(111). This revealed how adsorbate entropy affects the kinetic behavior.
Furthermore, Complete Potential Energy Sampling was developed to have a systematic
method for estimating entropies of adsorbates.
Additionally, this thesis has addressed the materials gap between extended surfaces
and nanoparticles in heterogeneous catalysis, where a microkinetic model was developed
for CO oxidation on Pt nanoparticles and Pt(111). The model was based on scaling
relations in the generalized coordination numbers, and the kinetics was simulated by
developing a kinetic Monte Carlo algorithm. Generalized coordination numbers provide
a computationally cheap procedure to account for the reaction energy landscape on
nanoparticles.
Finally, as a method to link theoretical simulations and experimental studies, the degree
of rate control was shown to be connected to reaction orders and apparent activation
energies. In this way, the present work has provided a microscopic understanding of these
two phenomenological quantities.
Compared to previous literature, this thesis has developed the methodology for
modeling reactions in heterogeneous catalysis. Previously, methane oxidation over Pd
surfaces has been modeled with a limited number of mechanisms. By inclusion of several
alternative reactions, a more reliable reaction mechanism was obtained in the present
work. Previously, reactions have mostly been simulated on extended surfaces. This thesis
has developed a scheme to facilitate simulations on nanoparticles by using the generalized
coordination numbers. The entropy of adsorbates has previously been of secondary concern
in microkinetic models. The present work has investigated the importance of entropy in
kinetic modeling, and it has developed the CPES method to calculate the entropy from
first principles. Reaction orders and apparent activation energies have previously not been
connected to the microscale. This work has enabled such an understanding by coupling
the reaction orders and apparent activation energies to the degree of rate control.
To further improve microkinetic modeling, there are several interesting aspects that can
be studied to extend this thesis. For example, the method of modeling adsorbate-adsorbate
interactions can be studied in greater detail, both for adsorption energies and entropies.
37
The developed kinetic Monte Carlo method can be applied for improved catalyst screening
studies, since it explicitly simulates nanoparticles. One limitation of the present work
is that it treats the catalyst as a static entity without including morphological changes.
Such dynamical changes can possibly be simulated by combining the developed Monte
Carlo algorithm with molecular dynamics. Furthermore, an oxide-support can be included
in the simulations to study possible spill-over phenomena.
Future catalyst design can benefit from understanding the active reaction mechanisms,
thermodynamics, and kinetic behavior in catalysis. In the future, when the computational
power has matured sufficiently, catalyst design can possibly be realized using computational
methods to a larger degree. An important part of this is the development of methodologies
for accurate microkinetic modeling.
As Richard Feynman envisioned it in his talk 1 almost 60 years ago, presently we have
a much greater range of possibilities for manipulating and designing nanoscaled systems.
When designing nanoscale technology, it is beneficial to understand the quantum mechan-
ical behavior of systems of atoms. This thesis has taken steps towards understanding the
detailed behavior of atomic systems.
38
Acknowledgments
The research presented herein was carried out at the Division of Chemical Physics and
Competence Centre for Catalysis at Chalmers University of Technology, Göteborg, Sweden
in the period August 2015 to July 2017.
The research is funded by Chalmers Area of Advance Nanoscience and Nanotechnology,

the Competence Centre for Catalysis, and the Swedish Research Council. Computational
time has been granted by SNIC at C3SE (Göteborg) and PDC (Stockholm).
The Competence Centre for Catalysis is hosted by Chalmers University of Technology

and financially supported by the Swedish Energy Agency and the member companies AB
Volvo, ECAPS AB, Haldor Topsøe A/S, Scania CV AB, Volvo Car Corporation AB and
Wärtsilä Finland Oy.
In addition I would like to acknowledge a non-exhaustive list of people:
My main supervisor, Henrik Grönbeck. Thank you for investing a lot of time and effort
in supervising this project and for our many scientific discussions. It has been very
developing to learn from your analytical skills and scientific insights.
My co-supervisor, Anders Hellman. You are acknowledged for serious as well as informal
discussions about theory in practice.
All my colleagues at Chemical Physics are acknowledged. Thank you for creating a
pleasant and nice working environment. I especially would like to thank Adam, Anna,
Alvaro, Baochang, Chris, Lin, Matthias, Maxime, Michael, Mikael, Oskar, and Unni.
I would also like to acknowledge the members of the Competence Centre for Catalysis.
In particular I would like to thank Carl-Robert, David, Emma, Johan, Peter, Ting, and
Torben.
At last but not least, I want to thank my family and Dorthe for making this possible by
providing me with extraordinary support and patience. Without you this thesis would
not have been.
Finally, thank you all for entering the same many-body wavefunction as I. This might
mean more than we realize.
39
40
Bibliography
[1] Feynman, R. P. There’s Plenty of Room at the Bottom. J. Microelectromech. Syst.,

1 (1992), 60–66. doi: 10.1109/84.128057.
[2] Chorkendorff, I. and Niemantsverdriet, J. W. Concepts of Modern Catalysis and
Kinetics. WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim, second, revised and
enlarged edition, 2007.
[3] Taylor, H. S. A Theory of the Catalytic Surface. Proc. R. Soc. Lond. A, 108 (1925),
105–111. doi: 10.1098/rspa.1925.0061.
[4] Hohenberg, P. and Kohn, W. Inhomogeneous Electron Gas. Phys. Rev., 136 (1964),
B864–B871. doi: 10.1103/PhysRev.136.B864.
[5] Kohn, W. and Sham, L. J. Self-Consistent Equations Including Exchange and
Correlation Effects. Phys. Rev., 140 (1965), A1133–A1138. doi: 10.1103/Phys-
Rev.140.A1133.
[6] Nørskov, J. K.; Abild-Pedersen, F.; Studt, F.; and Bligaard, T. Density Functional
Theory in Surface Chemistry and Catalysis. PNAS, 108 (2011), 937–943. doi:
10.1073/pnas.1006652108.
[7] Van Santen, R. A. and Neurock, M. Concepts in Theoretical Heteroge-
neous Catalytic Reactivity. Catal. Rev. Sci. Eng., 37 (1995), 557–698. doi:
10.1080/01614949508006451.
[8] Kozlov, S. M. and Neyman, K. M. Insights from Methane Decomposition on Nanos-
tructured Palladium. J. Catal., 337 (2016), 111–121. doi: 10.1016/j.jcat.2016.02.010.
[9] Fajı́n, J. L. C.; Bruix, A.; Cordeiro, M. N. D. S.; Gomes, J. R. B.; and Illas,
F. Density Functional Theory Model Study of Size and Structure Effects on
Water Dissociation by Platinum Nanoparticles. J. Chem. Phys., 137 (2012). doi:
10.1063/1.4733984.
[10] Jennings, P. C.; Aleksandrov, H. A.; Neyman, K. M.; and Johnston, R. L. A
DFT study of Oxygen Dissociation on Platinum based Nanoparticles. Nanoscale, 6
(2014), 1153–1165. doi: 10.1039/c3nr04750d.
[11] Viñes, F.; Lykhach, Y.; Staudt, T.; Lorenz, M. P. A.; Papp, C.; Steinrück, H.-P.;
Libuda, J.; Neyman, K. M.; and Görling, A. Methane Activation by Platinum:
Critical Role of Edge and Corner Sites of Metal Nanoparticles. Chem. Eur. J., 16
(2010), 6530–6539. doi: 10.1002/chem.201000296.
[12] Lopez-Acevedo, O.; Kacprzak, K. A.; Akola, J.; and Häkkinen, H. Quantum Size
Effects in Ambient CO oxidation Catalysed by Ligand-Protected Gold Clusters.
Nat. Chem., 2 (2010), 329–334. doi: 10.1038/NCHEM.589.
41
[13] Zheng, Z. and Greeley, J. Theoretical study of CO Adsorption on Au Catalysts
under Environmental Catalytic Conditions. Catal. Commun., 52 (2014), 78–83. doi:
10.1016/j.catcom.2014.03.016.
[14] Yudanov, I. V.; Sahnoun, R.; Neyman, K. M.; Rösch, N.; Hoffmann, J.; Schauermann,
S.; Johánek, V.; Unterhalt, H.; Rupprechter, G.; Libuda, J.; and Freund, H.-J. CO
Adsorption on Pd Nanoparticles: Density Functional and Vibrational Spectroscopy
Studies. J. Phys. Chem. B, 107 (2003), 255–264. doi: 10.1021/jp022052b.
[15] Posada-Borbón, A.; Heard, C. J.; and Grönbeck, H. Cluster Size Effects in Ethylene
Hydrogenation over Palladium. J. Phys. Chem. C, 121 (2017), 10870–10875. doi:
10.1021/acs.jpcc.6b12072.
[16] Kleis, J.; Greeley, J.; Romero, N. A.; Morozov, V. A.; Falsig, H.; Larsen, A. H.; Lu,
J.; Mortensen, J. J.; Dulak, M.; Thygesen, K. S.; Nørskov, J. K.; and Jacobsen,
K. W. Finite Size Effects in Chemical Bonding: From Small Clusters to Solids.
Catal. Lett., 141 (2011), 1067–1071. doi: 10.1007/s10562-011-0632-0.
[17] Li, L.; Larsen, A. H.; Romero, N. A.; Morozov, V. A.; Glinsvad, C.; Abild-Pedersen,
F.; Greeley, J.; Jacobsen, K. W.; and Nørskov, J. K. Investigation of Catalytic
Finite-Size-Effects of Platinum Metal Clusters. J. Phys. Chem. Lett., 4 (2013),
222–226. doi: 10.1021/jz3018286.
[18] Hammer, B. and Nørskov, J. K. Why Gold is the Noblest of All the Metals. Nature,
376 (2002), 238–240. doi: 10.1038/376238a0.
[19] Calle-Vallejo, F.; Martı́nez, J. I.; Garcı́a-Lastra, J. M.; Sautet, P.; and Loffreda,
D. Fast Prediction of Adsorption Properties for Platinum Nanocatalysts with
Generalized Coordination Numbers. Angew. Chem., Int. Ed., 53 (2014), 8316–8319.
doi: 10.1002/anie.201402958.
[20] Calle-Vallejo, F.; Loffreda, D.; Koper, M. T. M.; and Sautet, P. Introducing Struc-
tural Sensitivity into Adsorption–Energy Scaling Relations by means of Coordination
Numbers. Nat. Chem., 7 (2015), 403–410. doi: 10.1038/nchem.2226.
[21] Calle-Vallejo, F.; Tymoczko, J.; Colic, V.; Vu, Q. H.; Pohl, M. D.; Morgenstern, K.;
Loffreda, D.; Sautet, P.; Schuhmann, W.; and Bandarenka, A. S. Finding Optimal
Surface Sites on Heterogeneous Catalysts by Counting Nearest Neighbors. Science,
350 (2015), 185–189. doi: 10.1126/science.aab3501.
[22] Stoltze, P. and Nørskov, J. K. Bridging the ”Pressure Gap” between Ultrahigh-
Vacuum Surface Physics and High-Pressure Catalysis. Phys. Rev. Lett., 55 (1985),
2502–2505. doi: 10.1103/PhysRevLett.55.2502.
[23] Aparicio, L. M. and Dumesic, J. A. Ammonia Synthesis Kinetics: Surface Chemistry,

Rate Expressions, and Kinetic Analysis. Top. Catal., 1 (1994), 233–252. doi:
10.1007/BF01492278.
42
[24] Aparicio, L. M.; Rossini, S. A.; Sanfilippo, D. G.; Rekoske, J. E.; Treviño, A. A.;
and Dumesic, J. A. Microkinetic Analysis of Methane Dimerization Reaction. Ind.
Eng. Chem. Res., 30 (1991), 2114–2123. doi: 10.1021/ie00057a009.
[25] Dahl, S.; Sehested, J.; Jacobsen, C. J. H.; Törnqvist, E.; and Chorkendorff, I.
Surface Science Based Microkinetic Analysis of Ammonia Synthesis over Ruthenium
Catalysts. J. Catal., 192 (2000), 391–399. doi: 10.1006/jcat.2000.2857.
[26] Dumesic, J. A. and Trevino, A. A. Kinetic Simulation of Ammonia Synthesis

Catalysis. J. Catal., 116 (1989), 119–129. doi: 10.1016/0021-9517(89)90080-8.
[27] Frank, B.; Jentoft, F. C.; Soerijanto, H.; Kröhnert, J.; Schlögl, R.; and Schomäcker,
R. Steam Reforming of Methanol over Copper-Containing Catalysts: Influence
of Support Material on Microkinetics. J. Catal., 246 (2007), 177–192. doi:
10.1016/j.jcat.2006.11.031.
[28] Prasad, V.; Karim, A. M.; Arya, A.; and Vlachos, D. G. Assessment of Overall
Rate Expressions and Multiscale, Microkinetic Model Uniqueness via Experimental
Data Injection: Ammonia Decomposition on Ru/γ-Al2 O3 for Hydrogen Production.
Ind. Eng. Chem. Res., 48 (2009), 5255–5265. doi: 10.1021/ie900144x.
[29] Stoltze, P. and Nørskov, J. K. An Interpretation of the High-Pressure Kinetics of

Ammonia Synthesis Based on a Microscopic Model. J. Catal., 110 (1988), 1–10.
doi: 10.1016/0021-9517(88)90291-6.
[30] Falsig, H.; Hvolbæk, B.; Kristensen, I. S.; Jiang, T.; Bligaard, T.; Christensen, C. H.;
and Nørskov, J. K. Trends in the Catalytic CO Oxidation Activity of Nanoparticles.
Angew. Chem. Int. Ed., 47 (2008), 4835–4839. doi: 10.1002/anie.200801479.
[31] Honkala, K.; Hellman, A.; Remediakis, I. N.; Logadottir, A.; Carlsson, A.; Dahl, S.;
Christensen, C. H.; and Nørskov, J. K. Ammonia Synthesis from First-Principles
Calculations. Science, 307 (2005), 555–558. doi: 10.1126/science.1106435.
[32] Hansgen, D. A.; Vlachos, D. G.; and Chen, J. G. Using First Principles to Predict
Bimetallic Catalysts for the Ammonia Decomposition Reaction. Nat. Chem., 2
(2010), 484–489. doi: 10.1038/nchem.626.
[33] Hansen, E. W. and Neurock, M. Modeling Surface Kinetics with First-Principles-

Based Molecular Simulation. Chem. Eng. Sci., 54 (1999), 3411–3421. doi:
10.1016/S0009-2509(98)00489-8.
[34] Gokhale, A. A.; Dumesic, J. A.; and Mavrikakis, M. On the Mechanism of Low-
Temperature Water Gas Shift Reaction on Copper. J. Am. Chem. Soc., 130 (2008),
1402–1414. doi: 10.1021/ja0768237.
[35] Van den Bossche, M. and Grönbeck, H. Methane Oxidation over PdO(101) Revealed
by First-Principles Kinetic Modeling. J. Am. Chem. Soc., 137 (2015), 12035–12044.
doi: 10.1021/jacs.5b06069.
43
[36] Trinchero, A.; Hellman, A.; and Grönbeck, H. Methane Oxidation over Pd and
Pt studied by DFT and Kinetic Modeling. Surf. Sci., 616 (2013), 206–213. doi:
10.1016/j.susc.2013.06.014.
[37] Heard, C. J.; Hu, C.; Skoglundh, M.; Creaser, D.; and Grönbeck, H. Kinetic Regimes
in Ethylene Hydrogenation over Transition-Metal Surfaces. ACS Catal., 6 (2015),
3277–3286. doi: 10.1021/acscatal.5b02708.
[38] Mei, D.; Du, J.; and Neurock, M. First-Principles-Based Kinetic Monte Carlo
Simulation of Nitric Oxide Reduction over Platinum Nanoparticles under Lean-Burn
Conditions. Ind. Eng. Chem. Res., 49 (2010), 10364–10373. doi: 10.1021/ie100999e.
[39] Loffreda, D.; Delbecq, F.; Vigné, F.; and Sautet, P. Catalytic Hydrogenation of
Unsaturated Aldehydes on Pt(111): Understanding the Selectivity from First-
Principles Calculations. Angew. Chem. Int. Ed., 44 (2005), 5279–5282. doi:
10.1002/anie.200500913.
[40] Yang, L.; Karim, A.; and Muckerman, J. T. Density Functional Kinetic Monte
Carlo Simulation of Water-Gas Shift Reaction on Cu/ZnO. J. Phys. Chem. C., 117
(2013), 3414–3425. doi: 10.1021/jp3114286.
[41] Rogal, J.; Reuter, K.; and Scheffler, M. CO Oxidation on Pd(100) at Technologically
Relevant Pressure Conditions: First-Principles Kinetic Monte Carlo study. Phys.
Rev. B, 77 (2008), 155410. doi: 10.1103/PhysRevB.77.155410.
[42] Piccinin, S. and Stamatakis, M. CO Oxidation on Pd(111): A First-Principles-
Based Kinetic Monte Carlo Study. ACS Catal., 4 (2014), 2143–2152. doi:
10.1021/cs500377j.
[43] Kandoi, S.; Greeley, J.; Sanchez-Castillo, M. A.; Evans, S. T.; Gokhale, A. A.;
Dumesic, J. A.; and Mavrikakis, M. Prediction of Experimental Methanol Decom-
position Rates on Platinum from First Principles. Top. Catal., 37 (2006), 17–28.
doi: 10.1007/s11244-006-0001-1.
[44] Schlögl, R. Heterogeneous Catalysis. Angew. Chem. Int. Ed., 54 (2015), 3465–3520.
doi: 10.1002/anie.201410738.
[45] Österlund, L.; Grant, A. W.; and Kasemo, B. Lithographic Techniques in Nanocatal-
ysis, chapter Lithographic Techniques in Nanocatalysis, pages 269–341. Springer
Berlin Heidelberg, 2007. doi: 10.1007/978-3-540-32646-5.
[46] Schauermann, S. and Freund, H.-J. Model Approach in Heterogeneous Catalysis:
Kinetics and Thermodynamics of Surface Reactions. Acc. Chem. Res., 48 (2015),
2775–2782. doi: 10.1021/acs.accounts.5b00237.
[47] Papp, C. From Flat Surfaces to Nanoparticles: In Situ Studies of the Reactivity of
Model Catalysts. Catal. Lett., 147 (2017), 2–19. doi: 10.1007/s10562-016-1925-0.
[48] Freund, H.-J. Model Systems in Heterogeneous Catalysis: Selectivity Studies at the
Atomic Level. Top. Catal., 48 (2008), 137–144. doi: 10.1007/s11244-008-9052-9.
44
[49] Zhdanov, V. P. and Kasemo, B. Monte Carlo Simulation of the Kinetics of Rapid
Reactions on Nanometer Catalyst Particles. Surf. Sci,, 405 (1998), 27–37. doi:
10.1016/S0039-6028(97)01078-9.
[50] Persson, H.; Thormählen, P.; Zhdanov, V. P.; and Kasemo, B. Monte Carlo
Simulations of the Kinetics of Catalytic Reactions on Nanometer-Sized Particles,
with Diffusion over Facet Boundaries. J. Vac. Sci. Technol. A, 17 (1999), 1721–1726.
doi: 10.1116/1.581880.
[51] Eyring, H. The Activated Complex and the Absolute Rate of Chemical Reactions.
Chem. Rev., 17 (1935), 65–77. doi: 10.1021/cr60056a006.
[52] Mills, G.; Jónsson, H.; and Schenter, G. K. Reversible work Transition State Theory:
Application to Dissociative Adsorption of Hydrogen. Surf. Sci., 324 (1995), 305–337.
doi: 10.1016/0039-6028(94)00731-4.
[53] Kohanoff, J. Electronic Structure Calculations for Solids and Molecules: Theory
and Computational Methods. Cambridge University Press, New York (USA), 2006.
[54] Perdew, J. P.; Burke, K.; and Ernzerhof, M. Generalized Gradient Approxima-
tion Made Simple. Phys. Rev. Lett., 77 (1996), 3865–3868. doi: 10.1103/Phys-
RevLett.77.3865.
[55] Hammer, B.; Hansen, L. B.; and Nørskov, J. K. Improved Adsorption Energetics
within Density-Functional Theory using Revised Perdew-Burke-Ernzerhof Function-
als. Phys. Rev. B, 59 (1999), 7413–7421. doi: 10.1103/PhysRevB.59.7413.
[56] Haas, P.; Tran, F.; and Blaha, P. Calculation of the Lattice Constant of Solids
with Semilocal Functionals. Phys. Rev. B, 79 (2009), 085104. doi: 10.1103/Phys-
RevB.79.085104.
[57] Blöchl, P. Projector Augmented-Wave Method. Phys. Rev. B., 50 (1994), 17953–
17979. doi: 10.1103/PhysRevB.50.17953.
[58] Kittel, C. Introduction to Solid State Physics. John Wiley & Sons, Inc., Hoboken
(USA), 8th edition edition, 2005.
[59] Head, J. D. and Zerner, M. C. A Broyden–Fletcher–Goldfarb–Shanno Optimization

Procedure for Molecular Geometries. Chem. Phys. Lett., 122 (1985), 264–270. doi:
10.1016/0009-2614(85)80574-1.
[60] Henkelman, G.; Uberuaga, B. P.; and Jónsson, H. A Climbing Image Nudged Elastic
Band Method for Finding Saddle Points and Minimum Energy Paths. J. Chem.
Phys., 113 (2000), 9901–9904. doi: 10.1063/1.1329672.
[61] Hoffmann, M. J.; Medford, A. J.; and Bligaard, T. Framework for Scalable Adsorbate–
Adsorbate Interaction Models. J. Phys. Chem. C, 120 (2016), 13087–13094. doi:
10.1021/acs.jpcc.6b03375.
45
[62] Van den Bossche, M. and Grönbeck, H. Adsorbate Pairing on Oxide Surfaces: Influ-
ence on Reactivity and Dependence on Oxide, Adsorbate Pair, and Density Func-
tional. J. Phys. Chem. C, 121 (2017), 8390–8398. doi: 10.1021/acs.jpcc.6b12789.
[63] Jansen, A. P. J. An Introduction to Kinetic Monte Carlo Simulations of Surface

Reactions. Springer Berlin Heidelberg, 2012. doi: 10.1007/978-3-642-29488-4.
[64] Gould, H. and Tobochnik, J. Statistical and Thermal Physics With Computer
Applications. Princeton University Press, Woodstock, Oxfordshire, 2010. P 318.
[65] Hill, T. L. Introduction to Statistical Thermodynamics. Dover Publications Inc.,

New York, 1986. Chapter 9.
[66] Sprowl, L. H.; Campbell, C. T.; and Árnadóttir, L. Hindered Translator and
Hindered Rotor Models for Adsorbates: Partition Functions and Entropies. J. Phys.
Chem. C, 120 (2016), 9719–9731. doi: 10.1021/acs.jpcc.5b11616.
[67] Jones, E.; Oliphant, T.; Peterson, P.; and et al. SciPy: Open source scientific
tools for Python, 2001. URL https://docs.scipy.org/doc/scipy/reference/
tutorial/integrate.html, [Online; accessed 2017-05-10].
[68] Iserles, A. A First Course in the Numerical Analysis of Differ-

ential Equations. Cambridge University Press, second edition edition,
2009. Chapter 2 - Multistep Methods, [Books24x7 version] Available from
http://common.books24x7.com.proxy.lib.chalmers.se/toc.aspx?bookid=30922.
[69] Gillespie, D. T. A General Method for Numerically Simulating the Stochastic Time
Evolution of Coupled Chemical Reactions. J. Comput. Phys., 22 (1976), 403–434.
doi: 10.1016/0021-9991(76)90041-3.
[70] Chatterjee, A. and Vlachos, D. G. An Overview of Spatial Microscopic and Accel-

erated Kinetic Monte Carlo Methods. J. Comput.-Aided Mater. Des., 14 (2007),
253–308. doi: 10.1007/s10820-006-9042-9.
[71] Chatterjee, A. and Voter, A. F. Accurate Acceleration of Kinetic Monte Carlo

Simulations through the Modification of Rate Constants. J. Chem. Phys., 132
(2010), 194101. doi: 10.1063/1.3409606.
[72] Fichthorn, K. A. and Lin, Y. A Local Superbasin Kinetic Monte Carlo Method. J.
Chem. Phys., 138 (2013), 164104. doi: 10.1063/1.4801869.
[73] Puchala, B.; Falk, M. L.; and Garikipati, K. An Energy Basin Finding Algorithm
for Kinetic Monte Carlo Acceleration. J. Chem. Phys., 132 (2010), 134104. doi:
10.1063/1.3369627.
[74] Dybeck, E. C.; Plaisance, C. P.; and Neurock, M. Generalized Temporal Acceleration
Scheme for Kinetic Monte Carlo Simulations of Surface Catalytic Processes by
Scaling the Rates of Fast Reactions. J. Chem. Theory Comput., 13 (2017), 1525–
1538. doi: 10.1021/acs.jctc.6b00859.
46
[75] Zhu, G.; Han, J.; Zemlyanov, D. Y.; and Ribeiro, F. H. The Turnover Rate for the
Catalytic Combustion of Methane over Palladium Is Not Sensitive to the Structure
of the Catalyst. J. Am. Chem. Soc., 126 (2004), 9896–9897. doi: 10.1021/ja049406s.
[76] Zhu, G.; Han, J.; Zemlyanov, D. Y.; and Ribeiro, F. H. Temperature Dependence of
the Kinetics for the Complete Oxidation of Methane on Palladium and Palladium
Oxide. J. Phys. Chem. B, 109 (2005), 2331–2337. doi: 10.1021/jp0488665.
[77] Monteiro, R.; Zemlyanov, D.; Storey, J. M.; and Ribeiro, F. H. Turnover Rate and
Reaction Orders for the Complete Oxidation of Methane on a Palladium Foil in
Excess Dioxygen. J. Catal., 199 (2001), 291–301. doi: 10.1006/jcat.2001.3176.
[78] Campbell, C. T. Future Directions and Industrial Perspectives Micro- and macro-
kinetics: their Relationship in Heterogeneous Catalysis. Topics in Catalysis, 1
(1994), 353–366. doi: 10.1007/BF01492288.
[79] Stegelmann, C.; Andreasen, A.; and Campbell, C. T. Degree of Rate Control: How
Much the Energies of Intermediates and Transition States Control Rates. J. Am.
Chem. Soc., 131 (2009), 8077–8082. doi: 10.1021/ja9000097.
[80] McNaught, A. D. and Wilkinson, D. IUPAC. Compendium of Chemical Terminology.

Blackwell Scientific Publications, Oxford, 2nd ed. (the ”gold book”) edition, 1997.
[Online accessed 6th of April 2017 - https://doi.org/10.1351/goldbook.O04322 ].
[81] Gélin, P. and Primet, M. Complete Oxidation of Methane at Low Temperature

over Noble Metal Based Catalysts: a Review. Appl. Catal. B-Environ., 39 (2002),
1–37. doi: 10.1016/S0926-3373(02)00076-0.
[82] Hoflund, G. B.; Li, Z.; Epling, W. S.; Göbel, T.; Schneider, P.; and Hahn, H.
Catalytic Methane Oxidation Over Pd Supported on Nanocrystalline and Polycrys-
talline TiO2 Mn3O4, CeO2 and ZrO2. React. Kinet. Catal. Lett., 70 (2000), 97–103.
doi: 10.1023/A:1010362632223.
[83] Hao, H.; Liu, Z.; Zhao, F.; and Li, W. Natural Gas as Vehicle Fuel in
China: A Review. Renewable Sustainable Energy Rev., 62 (2016), 521–533. doi:
10.1016/j.rser.2016.05.015.
[84] Burch, R.; Loader, P.; and Urbano, F. Some Aspects of Hydrocarbon Activation on
Platinum Group Metal Combustion Catalysts. Catal. Today, 27 (1996), 243–248.
doi: 10.1016/0920-5861(95)00194-8.
[85] Farrauto, R. J.; Hobson, M. C.; Kennelly, T.; and Waterman, E. M. Catalytic
Chemistry of Supported Palladium for Combustion of Methane. Appl. Catal. A-Gen.,
81 (1992), 227–237. doi: 10.1016/0926-860X(92)80095-T.
[86] Lundgren, E.; Gustafson, J.; Mikkelsen, A.; Andersen, J. N.; Stierle, A.; Dosch, H.;
Todorova, M.; Rogal, J.; Reuter, K.; and Scheffler, M. Kinetic Hindrance during
the Initial Oxidation of Pd(100) at Ambient Pressures. Phys. Rev. Lett., 92 (2004),
046101. doi: 10.1103/PhysRevLett.92.046101.
47
[87] Ketteler, G.; Ogletree, D. F.; Bluhm, H.; Liu, H.; Hebenstreit, E. L. D.; and
Salmeron, M. In Situ Spectroscopic Study of the Oxidation and Reduction of
Pd(111). J. Am. Chem. Soc., 127 (2005), 18269–18273. doi: 10.1021/ja055754y.
[88] Baldwin, T. R. and Burch, R. Catalytic Combustion of Methane over Supported
Palladium Catalysts: II. Support and Possible Morphological Effects. Appl. Catal.,
66 (1990), 359–381. doi: 10.1016/S0166-9834(00)81649-8.
[89] Hicks, R. F.; Qi, H.; Young, M. L.; and Lee, R. G. Structure Sensitivity of Methane
Oxidation over Platinum and Palladium. J. Catal., 122 (1990), 280–294. doi:
10.1016/0021-9517(90)90282-O.
[90] Hicks, R. F.; Qi, H.; Young, M. L.; and Lee, R. G. Effect of Catalyst Structure on
Methane Oxidation over Palladium on Alumina. J. Catal., 122 (1990), 295–306.
doi: 10.1016/0021-9517(90)90283-P.
[91] Oh, S. H.; Mitchell, P. J.; and Siewert, R. M. Methane Oxidation over Alumina-
Supported Noble Metal Catalysts With and Without Cerium Additives. J. Catal.,
132 (1991), 287–301. doi: 10.1016/0021-9517(91)90149-X.
[92] Ciuparu, D.; Lyubovsky, M. R.; Altman, E.; Pfefferle, L. D.; and Datye, A. Catalytic
Combustion of Methane over Palladium-based Catalysts. Cat. Rev. - Sci. Eng., 44
(2002), 593–649. doi: 10.1081/CR-120015482.
[93] Qi, W.; Ran, J.; Wang, R.; Du, X.; Shi, J.; and Ran, M. Kinetic Mechanism
of Effects of Hydrogen Addition on Methane Catalytic Combustion over Pt(111)
Surface: A DFT Study with Cluster Modeling. Comp. Mater. Sci., 111 (2016),
430–442. doi: 10.1016/j.commatsci.2015.09.002.
[94] Zhu, Y.-A.; Chen, D.; Zhou, X.-G.; and Yuan, W.-K. DFT Studies of Dry Re-
forming of Methane on Ni Catalyst. Catal. Today., 148 (2009), 260–267. doi:
10.1016/j.cattod.2009.08.022.
[95] Bligaard, T.; Nørskov, J. K.; Dahl, S.; Matthiesen, J.; Christensen, C. H.; and
Sehested, J. The Brønsted–Evans–Polanyi Relation and the Volcano Curve in Hetero-
geneous Catalysis. J. Catal., 224 (2004), 206–217. doi: 10.1016/j.jcat.2004.02.034.
[96] Riedel, J. N.; Rötzer, M. D.; Jørgensen, M.; Vej-Hansen, U. G.; Pedersen, T.; Sebök,
B.; Schweinberger, F. F.; Vesborg, P. C. K.; Hansen, O.; Schiøtz, J.; Heiz, U.;
and Chorkendorff, I. H2/D2 Exchange Reaction on Mono-Disperse Pt Clusters:
Enhanced Activity from Minute O2 Concentrations. Catal. Sci. Technol., 6 (2016),
6893–6900. doi: 10.1039/C6CY00756B.
[97] Campbell, C. T.; Ertl, G.; Kuipers, H.; and Segner, J. A Molecular Beam Study of
the Catalytic Oxidation of CO on a Pt(111) surface. J. Chem. Phys., 73 (1980),
5862–5873. doi: 10.1063/1.440029.
[98] Park, J. Y.; Zhang, Y.; Grass, M.; Zhang, T.; and Somorjai, G. A. Tuning of Cat-
alytic CO Oxidation by Changing Composition of Rh-Pt Bimetallic Nanoparticles.
Nano Lett., 8 (2008), 673–677. doi: 10.1021/nl073195i.
48
[99] Nakao, K.; Watanabe, O.; Sasaki, T.; Ito, S.-I.; Tomishige, K.; and Kunimori,
K. CO Oxidation on Pd(111), Pt(111), and Rh(111) Surfaces studied by In-
frared Chemiluminescence Spectroscopy. Surf. Sci., 601 (2007), 3796–3800. doi:
10.1016/j.susc.2007.04.015.
[100] Palmer, R. L. and Smith Jr., J. N. Molecular Beam Study of CO Oxidation on a (111)
Platinum Surface. J. Chem. Phys., 60 (1974), 1453–1463. doi: 10.1063/1.1681219.
[101] Gerrard, A. L. and Weaver, J. F. Kinetics of CO Oxidation on High-Concentration

Phases of Atomic Oxygen on Pt(111). J. Chem. Phys., 123 (2005), 224703. doi:
10.1063/1.2126667.
[102] Vogel, D.; Spiel, C.; Suchorski, Y.; Trinchero, A.; and Schlögl, R. Local Catalytic
Ignition during CO Oxidation on Low-Index Pt and Pd Surfaces: A Combined
PEEM, MS, and DFT Study. Angew. Chem., Int. Ed., 51 (2012), 10041–10044.
doi: 10.1002/anie.201204031.
[103] Wrobel, R. J.; Becker, S.; and Weiss, H. Second/Additional Bistability in a CO
Oxidation Reaction on Pt(111): An Extension and Compilation. J. Phys. Chem. C,
116 (2012), 22287–22292. doi: 10.1021/jp302270n.
[104] Völkening, S. and Wintterlin, J. CO Oxidation on Pt(111) - Scanning Tunneling

Microscopy Experiments and Monte Carlo Simulations. J. Chem. Phys., 114 (2001),
6382–6395. doi: 10.1063/1.1343836.
[105] Calderón, S. K.; Grabau, M.; Óvári, L.; Kress, B.; Steinrück, H.-P.; and Papp, C.
CO Oxidation on Pt(111) at Near Ambient Pressures. J. Chem. Phys., 144 (2016),
044706. doi: 10.1063/1.4940318.
[106] Alavi, A.; Hu, P.; Deutsch, T.; Silvestrelli, P. L.; and Hutter, J. CO Oxidation
on Pt(111): An Ab Initio Density Functional Theory Study. Phys. Rev. Lett., 80
(1998), 3650–3653. doi: 10.1103/PhysRevLett.80.3650.
[107] Cheng, Z.; Fine, N. A.; and Lo, C. S. Platinum Nanoclusters Exhibit Enhanced
Catalytic Activity for Methane Dehydrogenation. Top. Catal., 55 (2012), 345–352.
doi: 10.1007/s11244-012-9803-5.
49
50

Microkinetic Modeling of Nanoparticle Catalysis Using Density Functional Theory

Uploaded by

Copyright:

Available Formats

Microkinetic Modeling of Nanoparticle Catalysis Using Density Functional Theory

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Microkinetic Modeling of Nanoparticle Catalysis Using Density Functional Theory

Uploaded by

Copyright:

Available Formats

Thesis for the degree of Licentiate of Engineering

Microkinetic Modeling of Nanoparticle Catalysis

c Mikkel Jørgensen, 2017

Printed at Chalmers Reproservice

Keywords: Catalysis, Density Functional Theory, Microkinetic modeling, Nanoparti-

This thesis is based on the following appended papers:

2 Calculating Rate Constants 5

3 Simulating Reaction Kinetics 17

4 Kinetic Models for Methane and CO Oxidation 27

5 Conclusion & Outlook 37

Figure 1.1: Schematic illustration of a catalytic reaction. The catalyst provides an

One technology that benefits from atomistic understanding is heterogeneous catalysis.

where Ea is the activation energy, A is a pre-exponential factor, and T is the temperature.

1.1 Thesis Objectives

Calculating Rate Constants

To simulate reaction kinetics, it is necessary to calculate the rate constants of the

2.1 Reaction Rate Constants

2.2 Density Functional Theory

Ĥ = ĤeZ + Ĥee + ĤZZ . (2.9)

2.2.1 Exchange Correlation Energy Approximations

The exchange-correlation hole and pair-correlation function for an inhomogeneous electron

2.2.2 Projector-Augmented Wave Method

Where I is the identity operator and he pi | is a projector function. There is one he

2.2.3 Periodic Systems

2.3 Local Minimum Energies

2.4 Energy Barriers with Nudged Elastic Band

2.5 Adsorbate-Adsorbate Interactions

a coverage of both C and O. Thus, the adsorbate-adsorbate interactions are of utmost

2.6 Modeling Entropies in Catalytic Reactions

In this way, the entropy is uniquely defined by the partition function.

Z = Ztrans Zrot Zvib . (2.24)

Simulating Reaction Kinetics

3.1 Chemical Master Equation

3.2 Mean-Field Approximation

dhNA i X X h (A) (A)

A2 (g) + 2∗ ↔ 2A∗ (a)

The time-evolution of the coverage of A is given by:

3.3 Kinetic Monte Carlo

1. Initialize Simulation : Set simulation time (t = 0), simula-

4. Perform next possible event in queue, and set t to tβα of the

Figure 3.1: First Reaction Method Algorithm.

3.4 Analyzing Reaction Kinetics

3.4.1 Reaction Mechanisms

3.4.2 Turnover Frequencies and Coverages

3.4.3 Reaction Orders and Apparent Activation Energies

where K is a constant with respect to pressure, px is the partial pressure in species x,

where A is a the pre-exponential factor, x is a gas-species, px is the pressure, and nx is

3.4.4 Degree of Rate Control

Kinetic Models for Methane and CO

4.1 Constructing a Kinetic Model

Figure 4.1: Illustration of workflow for constructing a first-principles microkinetic

CH4 (g) + 5∗ → C∗ + 4H∗ (4.1)

4.2 Complete Methane Oxidation over Pd(100) and Pd(111)

4.3 Entropic Effects in CO Oxidation on Pt(111)

The entropies of CO and O are found to be qualitatively different. CO translates on a

i For example either fcc, hcp, bridge, or ontop

effect is enhanced with increasing temperatures as ∆G = ∆E − T ∆S.

4.4 Catalysis on Nanoparticles

generalized coordination number CN was shown to be a good descriptor for adsorption