Computational Catalyst Screening

Grabow, Lars C.

doi:10.1039/9781839163296-00001

Skip Nav Destination

Book Chapter

Chapter 1: Computational Catalyst Screening

William A. Brookshire Department of Chemical & Biomolecular Engineering, University of Houston, Houston, Texas, 77204-4004, USA [email protected]

Search for other works by this author on:

This Site

PubMed

Google Scholar

Doi:

https://doi.org/10.1039/9781839163296-00001

Published:

20 Dec 2024
Special Collection: 2024 eBook Collection

Series: Catalysis Series

Page range:

1 - 59

Computational screening of heterogeneous catalysts based on reactivity descriptors is a very powerful method for rapid identification of promising novel catalyst candidates. This chapter outlines the overall procedure based on examples from the literature and provides step-by-step instructions with solved numerical problems for NH₃ synthesis and CO oxidation on transition metal surfaces. The theoretical foundations of the screening approach, including the d-band model, linear scaling relations, Sabatier analysis, basic microkinetic modeling and the analysis of such models, are explained at the level necessary for a novice to perform an independent screening study for other transition metal catalyzed reactions. More experienced readers may find the references for suggested additional reading and resources useful.

1.1 Introduction

Brute-force attacks are known in cryptography as (typically illegal) attempts to hack into encrypted data by systematically trying all possible key combinations of letters, digits and special characters until the correct access key or password has been found. Although a brute-force attack is guaranteed to be successful, its application is limited to very small problems because of the time required to generate and test all possible key combinations. For example, a standard 128-bit encryption key has 2¹²⁸ possible permutations. If we simply assume that a typical central processing unit (CPU) can generate 10⁹ bit flips per second (approximately 1 GHz), then the total time that is required to test all possible permutations is 2¹²⁸/10⁹ which equals 3.4 × 10²⁹ seconds or 10²² years! For obvious reasons a brute-force attack is probably going to fail for this problem and a more targeted strategy is needed.

The above example illustrates the shortcomings of a brute-force attack, but a variation of it is still one of the most widely used strategies for the development of heterogeneous catalysts in practice. By using a combinatorial chemistry approach with completely automated, high-throughput experimentation equipment, it is possible to synthesize and test enormous libraries of catalysts for their catalytic activity for a specific reaction. A good example is the search for advanced water–gas shift catalysts, in which Yaccato et al. synthesized over 50 000 catalysts and tested them in more than 250 000 experiments for low, medium and high temperature water–gas shift conditions.¹ Their effort led to a proprietary noble metal catalyst that can reduce the reactor volume by an order of magnitude without increasing the reactor cost. Although this trial-and-error approach almost always leads to an acceptable catalyst, the search space is restricted by the amount of time and resources available and many, possibly far better, candidates can be missed. The quickly evolving alternative to experimental high-throughput catalyst testing is computational catalyst screening. This approach relies on the fact that the catalyst activity for many catalytic reactions is usually determined by a small number of descriptors, which can be calculated from first-principles density functional theory (DFT) simulations and stored in a database. Populating this database with features or properties from data is the most time-consuming step in this process, but the resulting database has to be generated only once and its content may be used to engineer features that describe catalytic activity trends for other reactions. With a comprehensive database in place, it becomes a very easy task to screen thousands of database records in a short amount of time to identify catalyst candidates that possess descriptor values within the optimal range for a given reaction. Although the computational screening process can still be interpreted as a brute-force attack, the complexity of the problem has been greatly reduced. Hence, the number of materials that can be screened computationally increases drastically when compared with the experimental approach. The list of catalysts that fall into the desired range of descriptor values may be narrowed down further by using cost, stability, environmental friendliness, or any other applicable criteria.^2,3 The remaining materials can then be synthesized and experimentally tested under realistic reaction conditions. In general, not all computationally screened candidates will be good catalysts, but good catalysts will usually be included in the candidate list.

Somorjai and Li have recently reviewed the major advances in modern surface science that only became possible through the successful symbiosis of theory and surface-sensitive experimental techniques.⁴ The recent literature also contains several examples where a descriptor-based approach, both theoretically and experimentally, has led to the discovery of new catalytic materials. The following list should not be understood as an exhaustive review, but is meant to serve as inspiration to the reader and to demonstrate the wide applicability of this method. Early on, Besenbacher et al. discovered graphite resistant Ni–Au alloy catalysts for steam reforming,⁵ Jacobsen et al. found an active Co–Mo alloy for ammonia synthesis by interpolation in the periodic table,⁶ and Toulhoat and Raybaud showed that the metal–sulfur bond strength can correctly predict trends in hydro-desulfurization activity on metal–sulfide catalysts.⁷ These initial successes were followed by other prominent examples that include CO-tolerant fuel cell anodes,^8,9 Cu–Ag alloys as selective ethylene epoxidation catalysts,¹⁰ near-surface alloys for hydrogen activation¹¹ and evolution,^2,12 Ru–Pt core–shell particles for preferential CO oxidation,¹³ Ni–Zn alloys for the selective hydrogenation of acetylene,¹⁴ Sc and Y modified Pt and Pd electrodes¹⁵ and mixed-metal Pt monolayer catalysts¹⁶ for electrochemical oxygen reduction, and the rediscovery of Pt as the most active and selective catalyst for the production of hydrogen cyanide.^17,18

1.1.1 A Walk Through a Computational Catalyst Design Process: Methanation

The most comprehensive example of a success story in computational catalyst design comes from the group of Jens Nørskov, who pioneered the descriptor-based design approach and applied it to numerous reactions.^19,20 In several publications his group has studied the methanation reaction (CO + 2H₂ → CH₄ + H₂O), starting from a detailed electronic structure analysis and leading to the development of a patented technical methanation catalyst based on a Fe–Ni alloy.^21–25 In the beginning of any descriptor-based design study the first question that must be answered is: “What is the most suitable reactivity descriptor for the reaction?” This question is typically answered by thoroughly studying the underlying reaction mechanism and identifying the rate-limiting step and most abundant surface intermediates. However, intuition can sometimes replace a detailed mechanistic study and a descriptor can be found through an educated guess. In the case of the methanation reaction, CO dissociation is the most critical step in the reaction mechanism. For weakly interacting metal catalysts, the dissociation is rate limiting, whereas for strongly interacting catalysts, the surface is poisoned by adsorbed C and O atoms. This leads to the volcano curve shown in Figure 1.1(a), which shows the experimentally measured methanation activity as a function of the calculated CO dissociation energy. The top of the volcano corresponds to the maximum methanation activity and indicates the optimal value of the CO dissociation energy, which is the activity descriptor in this case. The next step in the catalyst design process is to screen a database of CO dissociation energies and search for catalysts with CO dissociation energies near the optimum. This screening may be combined with a cost estimation of the resulting material and can further be linked to a stability test. Figure 1.1(b) and (c) show Pareto plots of binary transition metal alloys for which the CO dissociation energy was estimated through a simple interpolation scheme owing to the lack of an existing database. The most active catalysts, characterized by CO dissociation energies close to the optimum, lie to the left of the graph and are connected with a solid line indicating the Pareto-optimal set. The Pareto-optimal set of Figure 1.1(b) contains the cheapest catalysts for a given value of CO dissociation energies and, similarly, Figure 1.1(c) can be used to screen for alloy stability. Only alloys with a negative alloy formation energy are stable and their stability increases as the alloy formation energy becomes more negative. Upon careful inspection, it is apparent that FeNi₃ is not only contained in both sets, but it is also located at the “knee” of the activity Pareto-optimal set, which indicates that neighboring solutions are worse with respect to either activity or cost. Clearly, FeNi₃ is a very promising catalyst candidate for the methanation reaction. This catalyst identification step concludes the theoretical design process and experimental verification of the theoretical prediction is necessary. Experimentally obtained methanation rates of Fe–Ni alloys as a function of Ni content are displayed in Figure 1.1(d) and clearly show that the computationally predicted FeNi₃ alloy is significantly more active than the alternatives. As an outcome of this tour de force in computational catalyst design, a process based on Fe–Ni alloys has been patented for the hydrogenation of carbon oxides.²⁴

Figure 1.1

View large Download slide

Computational design of a technical methanation catalyst. (a) A characteristic volcano curve is obtained when the experimentally determined methanation activity is plotted as a function of the CO dissociation energy. (b) and (c) Pareto-optimal bimetallic catalysts in terms of cost and stability. E_diss(optimal) refers to the optimal CO dissociation energy corresponding to the maximum of the activity volcano in panel (a). (d) Measured methanation activity of binary Fe–Ni alloys at T = 548 K, 2% CO in 1 bar H₂ as a function of Ni content. (a), (b) and (d) Reproduced from ref. 15 with permission from Springer Nature, Copyright 2009. (c) Reproduced from ref. 25 with permission from IOP Publishing, Copyright 2009.

In the remainder of this chapter, the background information that leads to the identification of appropriate catalyst descriptors (e.g. the d-band model, scaling relationships) is reviewed and the basic strategy for successful catalyst screening using various levels of detail (e.g. Sabatier rate vs. microkinetic model) is outlined. A step-by-step illustration of the method will be given using ammonia synthesis and CO oxidation as examples. The interested reader is encouraged to work through the examples independently at their own pace.

1.2 Starting from the Electronic Structure

1.2.1 Density Functional Theory

Computational catalyst screening would not be possible without the existence of a theory that enables calculation of the chemical properties of the catalyst and the reaction of interest. Fortunately, in the mid 1960s Hohenberg, Kohn and Sham published two seminal papers formulating two theorems, which led to the development of density functional theory (DFT).^26,27 The contributions of Walter Kohn to the development of this theory were later honored in 1998 with the Nobel Prize in Chemistry. DFT is nowadays widely used in many different areas of science and engineering, including computational chemistry, catalysis, materials science, physics, and geology. The two theorems can be summarized as:

Theorem 1. The ground state properties of a many-electron system are uniquely determined by the electron density.
Theorem 2. The total energy of a system has a minimum for the ground state electron density.

DFT provides a solution to the Schrodinger equation:

\hat{H} Ψ = E Ψ

(1.1)

where

\hat{H}

is the Hamiltonian operator, Ψ is the wave function and E is the energy. It is, in principle, an exact theory, but in practice the exact formulation of the kinetic energy term for a system of interacting electrons is unknown. In the Kohn–Sham approach, the kinetic energy is therefore approximated with the kinetic energy of a system of non-interacting electrons and a correction term, E_XC, which accounts for exchange and correlation effects in the interacting system. Although approximations for the description of the exchange–correlation energy must be made, DFT has the huge advantage over wave function-based methods that the electron density is a function of only three spatial coordinates, while the many-body wave function for N electrons depends on 3N coordinates. Thus, DFT significantly reduces the computational intensity of the problem and enables the treatment of systems of several hundred atoms. From the electron density, n( r ), all other properties of the system are determined (Theorem 1) and the total energy E is calculated using

\begin{array}{l} E [n (r)] & = & T_{KS} [n (r)] + \frac{1}{2} \iint \frac{n (r) n (r r)}{| r - r r |} d r d r r \\ + \int v_{nucl} (r) n (r) d r + E_{XC} [n (r)] \end{array}

(1.2)

The first term, T_KS[n( r )], is the kinetic energy of fictitious, non-interacting electrons and is obtained from the single-electron Kohn–Sham equation

(- \frac{ℏ^{2}}{2 m} \nabla^{2} + v_{eff}) ψ_{i} = ε_{i} ψ_{i}

(1.3)

where v_eff is the effective field defined by the nuclei and the current electron density, ℏ is the reduced Planck’s constant, m is the mass, ψ_i is the wavefunction of electronic state i, and ε_i is the Kohn–Sham eigenvalue. The second and third term in the total energy equation describe the electrostatic electron–electron interactions (Hartree energy) and the electron–nuclei interactions, respectively. The last term, E_XC, in the total energy equation depends on the unknown exchange–correlation functional, for which several approximations exist. The simplest approximation is the local density approximation (LDA), which can be derived from the case of a homogeneous electron gas and only depends on the electron density at a single point. In this case, the exchange contribution in the LDA is exact, but the correlation still has to be approximated. The LDA works remarkably well for bulk materials where the electron density varies slowly, but has insufficient accuracy for most applications in chemistry, including atoms, molecules, clusters, and surfaces.

An obvious extension to the LDA is the generalized gradient approximation (GGA), which depends not only on the local density but on the density gradient. Because the gradient correction can be implemented into a GGA functional in many different ways, there exist a variety of different GGA flavors. The most widely used GGA functionals are the Perdew–Wang 91 (PW91)^28,29 and the Perdew–Burke–Ernzerhof (PBE)³⁰ functional. Both GGA functionals have good accuracy for a wider range of problems than the LDA because they contain more physical information; however, they are not necessarily always better. The PBE functional was later revised by Hammer, Hansen and Nørskov [revised PBE (RPBE)], to improve the accuracy of chemisorption energies of atoms and small molecules on transition metal surfaces.³¹ These GGA functionals are very good general-purpose functionals and may be used as a starting point for computational catalyst design. However, the GGA still fails for problems such as the accurate prediction of band gaps in semiconductors, systems where van der Waals interactions are dominant, or for electronic structure calculations in materials with strongly correlated electrons, where self-interaction errors can be encountered. Several improvements to the GGA have been suggested [e.g. DFT + on-site Coulomb interactions (DFT + U), dispersion corrected DFT (DFT-D), meta-GGA, hybrid-GGA], but many of these functionals are problem specific or contain adjustable parameters that need to be fitted for each system. This empirical nature, along with the increased computational effort, renders these functionals generally unsuitable for computational catalyst screening. Work to improve XC functionals further is ongoing in the community, and in the next few years faster computers and new functionals will have a positive effect on the quality of DFT calculations.

Although DFT calculations are at the heart of computational catalysis, it is not strictly necessary to perform the calculations for every catalyst design project. DFT calculations for many reactions have already been published and efforts are undertaken to make these data easily available to the whole catalysis community, even on mobile devices. Yes, there is an app for that!³² However, even a non-DFT expert should understand the basic principles that underlie the theoretical results, before using them in a research project. For those readers that have a deeper interest in DFT and want to perform their own calculations, the tutorial-style book Density Functional Theory – A Practical Introduction by David S. Sholl and Janice A. Steckel is a highly recommended starting point.³³

1.2.2 The d-Band Model

Computational catalyst screening relies on the prediction of correct trends across different catalysts rather than the prediction of quantitative rates and selectivities for each catalyst. Understanding the origin of the observed trends in terms of the underlying electronic structure can therefore be very helpful during the screening process. For transition metal surfaces, trends in reactivity can be very well described and understood in terms of the d-band model (Figure 1.2) developed by Hammer and Nørskov.^34,35

Figure 1.2

View large Download slide

Schematic density of states (DOS) illustration of the d-band model. The interaction of an adsorbate state with a transition metal can be thought of as a two-step process. The interaction with the broad s-band leads to a broadening and downshift of the adsorbate states. The adsorbate states split into bonding and anti-bonding states upon interaction with the narrow transition metal d-band. Anti-bonding states that are above the Fermi level remain empty and do not weaken the chemisorption. Reproduced from ref. 20 with permission from National Academy of Sciences; and original content of figure reproduced from ref. 71 with permission from Springer Nature, Copyright 2005.

A schematic drawing of the bonding structure in a hydrogen molecule is likely to have been presented in chemistry classes. Upon bond formation, the two atomic orbitals form two new molecular orbitals, and they can be classified as a bonding and an anti-bonding orbital. In the case of a hydrogen molecule, two electrons can distribute into these orbitals and since each orbital can accommodate up to two electrons, naturally both electrons occupy the lower-lying bonding orbital. The energy that is gained by stabilizing the electrons in this process is the bond energy.

The interactions of adsorbates with transition metals are a bit more complicated but conceptually similar. A transition metal does not possess atomic orbitals but has a continuous range of available states called a “band”. In a simplified picture of the band structure in transition metals (Figure 1.2, right) it can be imagined that they have a broad s-band (turquoise) and a narrow band of localized d-electrons (red). Assuming that the coupling between the adsorbate level and the s- and d-bands of the transition metal can be separated the chemisorption energy can be written as the sum of the two interactions.

Δ E_{Chem} = Δ E_{s} + Δ E_{d}

(1.4)

The interaction of the adsorbate state with the broad s-band, ΔE_s, leads to a broadening of the state and a downshift in energy. Since all transition metals have broad and half-filled s-bands, the contribution ΔE_s, is approximately the same for all transition metals. In the next step, the broadened molecular state can couple to the narrow d-band of the transition metal, which causes a split into bonding and anti-bonding states, just as in the case of the molecular orbitals of a hydrogen molecule. However, in contrast to the molecular system, where a fixed number of electrons is available to occupy the new orbitals, a transition metal has a large reservoir of electrons. These electrons will fill up all states located below the Fermi level and, consequently, the more anti-bonding states are located below the Fermi level, the weaker the resulting chemisorption bond. For stronger bonds it is desirable that all the anti-bonding states are high in energy and above the Fermi level.

The location of the anti-bonding states relative to the Fermi level and, in turn, the strength of the resulting bond, is largely determined by the position of the d-band. Since the d-band spans a range of energies, a convenient way to define the term “position” needs to be found, and therefore the d-band center is introduced. The d-band center,

ϵ_{d}

⁠, is the energy-weighted average of the density of d-states ρ. It sounds complicated, but it is the same equation that is used to calculate the center of mass (just substitute

r_{COM} = ϵ_{d}

⁠, r = E and m = ρ).

ϵ_{d} = \frac{\int ρ E d E}{\int ρ d E}

(1.5)

The d-band center varies across the transition metals and according to the d-band model this will cause a variation in the interaction strength, ΔE_d. Since ΔE_s is constant, the variations in chemisorption energy across transition metals can be attributed to changes in ΔE_d and, thus the d-band center. Transition metals with higher-lying d-bands have stronger chemisorption properties. Numerous examples supporting the d-band model theory exist in the literature (ref. 20 and 22–37 therein) but, as with most rules, there are some exceptions. In systems including electronegative adsorbates with nearly filled valence shells (e.g. OH, F, Cl) on surfaces with almost fully occupied d-bands (d⁹ and d¹⁰ metals) there can be a significant repulsive interaction between the transition metal d-band and the adsorbate states, which leads to a local reversal of the d-band model trend.³⁶ For computational screening studies over large descriptor domains, however, these exceptions only affect the local fine structure and general trends are preserved.

1.3 Identifying the Right Descriptor Set

If you have ever played the party game “Taboo” you have experienced how hard it can be to describe a target word or object without using the five most closely related words listed on the card. On the other hand, if you were allowed to mention only the five forbidden words, your teammates would probably guess the target word immediately. Think of these five taboo words as the descriptors that we need to guess the performance of our catalyst. Choosing the right descriptors is the key to a successful descriptor-based catalyst design study and it needs to be done carefully.

The ideal set of descriptors needs to fulfill two conflicting requirements. First and foremost, the descriptor set has to be large enough to enable predictions that are sufficiently accurate to serve the purpose of the study. For simple reactions, a single descriptor may be sufficient to describe qualitatively the activity trends across different catalyst materials. However, product selectivity, for example, is often more sensitive to the input parameters and may require additional descriptors. If quantitative results are desired, then the number of required descriptors quickly approaches the total number of enthalpy and entropy parameters in a reaction network, which defeats the purpose of reducing the complexity of the problem by introducing the concept of descriptors. In fact, reducing the complexity is the second most important requirement that the descriptors must meet. The set of descriptors should be as small as needed to capture catalytic trends and enable fast and efficient screening for new catalyst materials.

There is no strict rule for the “correct” choice of descriptors; in fact, for most cases there is more than one set of descriptors and all of them may be equally viable. A descriptor can be any measurable intrinsic quantity of the catalyst (e.g. the d-band center), but most often descriptors are binding energies of key intermediates inferred from the detailed knowledge of the dominant reaction mechanism and the kinetically relevant steps. This information can be obtained from mechanistic studies using DFT in combination with kinetic measurements and modeling. Alternatively, the existence of scaling relations for surface intermediates and transition states can guide the descriptor selection process. These scaling relations have their physical roots in the d-band model and are discussed next.

1.3.1 Scaling Relations for Surface Intermediates

The binding energy of an adsorbed molecule is determined by the number and strength of the chemical bonds that it forms with the surface. To a first approximation these bonds can be considered to form independently of each other and it can be assumed that the bond strength only depends on the two types of atoms involved in the bond formation. In this simple picture, it would be sufficient to know through which atoms a molecule binds to the catalyst surface in order to predict its adsorption energy. Indeed, it has been shown that the adsorption energies of adsorbates within a family of similar adsorbates can be predicted using this simple idea. The resulting linear scaling relationships can be used to predict unknown adsorption energies from the adsorption energies of related surface species.^37,38 The accuracy of linear scaling relations is generally adequate to predict catalytic trends correctly, but the average errors (0.2–0.3 eV) are too large for quantitative predictions. However, in the context of computational screening, in which the relative order of different transition metal catalysts is the sole subject of interest, these errors are acceptable.

1.3.1.1 Hydrogen-containing Molecules

The discussion in this section is fully based on the scaling papers by Abild-Pedersen et al.³⁷ and Fernandez et al.,³⁸ which are highly recommended references on this topic. Each of the four panels in Figure 1.3 shows the adsorption energies, ΔE^AH_x, of the hydrogen-containing intermediates AH_x (A = C, N, O, S) plotted as a function of the adsorption energy, ΔE^A of the central atom A for a range of typical transition metal surfaces. The adsorption energies were obtained from periodic DFT calculations on close-packed terraces (black), step sites (red) and, additionally, on the face-centered cubic [FCC(100)] surface (blue) for OH_x. It can be seen that the adsorption energies ΔE^AH_x are linearly correlated with ΔE^A and given by

Δ E^{{AH}_{x}} = γ Δ E^{A} + ξ,

(1.6)

where γ is the slope and ξ is the intercept. There is certainly some scatter around the linear scaling lines, but the general trend is correctly captured and the absolute errors are within the 0.2–0.3 eV range. Upon inspection of Figure 1.3, the observant reader may have noticed that the scaling lines look different for different adsorbates and surface geometries, i.e. flat vs. stepped surfaces. But is there any physical significance to these differences?

Figure 1.3

View large Download slide

Linear scaling relationships on close-packed terraces (black), step sites (red) and the FCC(100) surface (blue) of hydrogen-containing molecules of the type AH_x with A = C, N, O, S. Reproduced from ref. 37 with permission from American Physical Society, Copyright 2007.

The first observation is that the slope of the scaling relation does not depend on the surface geometry and is constant for any given adsorbate. This is best seen in the case of the OH adsorption energy, ΔE^OH, which scales with a slope of approximately 1/2 with respect to the adsorption energy of oxygen, ΔE^O, on all three different surfaces. The same generalization applies to all the other adsorbates and can easily be rationalized using the simple binding picture, where the adsorption energy depends on the number and the type of bonds a molecule forms with the surface. For example, carbon has four valence electrons and likes to form four C–H bonds. If these H atoms are removed in a stepwise fashion, the unsatisfied carbon valencies start to interact with the catalyst surface and form C–M bonds (‘M’ denotes a metal atom on the catalyst surface) as depicted in Figure 1.4. When moving from one metal to another the C–M bond strength changes and the simple bond order method can be used to estimate how this change affects the adsorption energy of the CH_x species.

Figure 1.4

View large Download slide

Simplified bonding scheme for adsorbed CH_x species.

In general, the scaling slope γ in for hydrogen-containing adsorbates AH_x can then be predicted by

γ (x) = \frac{x_{\max} - x}{x_{\max}}

(1.7)

where x_max is the maximum number of H atoms that can bind to the central atom A. Specifically, x_max = 4 for C, x_max = 3 for N and x_max = 2 for O and S. It should be noted that the linear scaling behavior and characteristic slope of hydrogen-containing adsorbates is not limited to adsorption on metal surfaces. The same linear scaling relations have been discovered for transition metal nitrides, oxides and sulfides, which broadens the range of possible applications of scaling relations far beyond catalysis on transition metal surfaces.³⁸

Now that the origin of the slope in the linear scaling relations is understood, attention can be focused on the intercept ξ. Unfortunately, the intercept cannot be as easily predicted as the slope, but there are two important aspects that need to be mentioned. Using the bond order conservation method to estimate the scaling slope, a single reference calculation can be sufficient to predict adsorption energies on any other catalyst surface. Assume that ΔE^AH_x and ΔE^A were calculated on a reference catalyst and x_max is known to predict the slope, γ, then can be used to predict ΔE^AH_x as a function of ΔE^A by simply eliminating the intercept ξ.

Δ E^{{AH}_{x}} = Δ E_{ref}^{{AH}_{x}} + γ (x) \times (Δ E^{A} - Δ E_{ref}^{A})

(1.8)

Structure sensitivity is the other factor that affects the intercept, ξ, which is clearly demonstrated by the vertical offset between the parallel scaling lines for each adsorbate in Figure 1.3. Although structure sensitivity plays an important role in catalysis, its effect on the scaling behavior of adsorbates is not discussed here and the interested reader is referred to the detailed tutorial review by Nørskov et al. instead.³⁹

1.3.1.2 Extension to More Complex Surface Intermediates

The special case of hydrogen-containing AH_x species impressively demonstrates the relationship between the binding energies of similar adsorbates, but in the majority of catalytic reactions more complex surface intermediates are encountered. In the process of establishing linear scaling relations for complex species, it is possible to start intuitively by correlating their adsorption energies to the functional groups contained in the species or simply to the adsorption energies of the atoms through which the molecule binds. This approach is especially useful when the adsorption geometry of the involved surface species is known. Let us take the bond conservation idea introduced in the previous section for AH_x and apply it to CH_x–NH_y with x,y = 0–2 surface intermediates. These intermediates occur, for example, during the Degussa and Andrussow process for HCN synthesis from CH₄ and NH₃.^17,18,40,41 It may be assumed that the C and N atoms are connected through a single bond and that the CH_x–NH_y intermediates bind to the surface equally through the C atom and the N atom. The adsorption energy of CH_x–NH_y should then scale linearly as in

Δ E^{{CH}_{x} - {NH}_{y}} = γ \times (Δ E^{{CH}_{x}} + Δ E^{{NH}_{y}}) + ξ

(1.9)

where

Δ E^{{CH}_{x}} = \frac{4 - 1 - x}{4} Δ E^{C} + ξ_{C}

and

Δ E^{{NH}_{x}} = \frac{3 - 1 - x}{3} Δ E^{N} + ξ_{N}

⁠. Note the ‘−1’ term indicating that one valency of C and N is consumed by the C–N bond formation and is no longer available. If the expressions for ΔE^CH_x and ΔE^NH_y are substituted into and the intercepts collected, it is possible to obtain the final expression for

Δ E^{({CH}_{x} - {NH}_{y})}

⁠.

Δ E^{{CH}_{x} - {NH}_{y}} = γ_{{CH}_{x} - {NH}_{y}} (\frac{4 - 1 - x}{4} Δ E^{C} + \frac{3 - 1 - x}{3} Δ E^{N}) + ξ_{{CH}_{x} - {NH}_{y}}

(1.10)

The comparison of CH_x–NH_y adsorption energies estimated from with DFT reference calculations given in Figure 1.5(a) shows that this approach adequately reproduces the correct trend, although several outliers with errors of approximately 0.5 eV exist. A similar strategy based on bond order conservation was used by Jones et al. to derive scaling relations for C₂ hydrocarbons on transition metal surfaces.⁴²

Figure 1.5

View large Download slide

Linear scaling relations for complex adsorbates. (a) Parity plot of CH_x–NH_y adsorption energies predicted from scaling and the DFT reference calculations. (b) Linear scaling of C₂H₄ and C₂H₆ with methyl (CH₃). Reproduced from ref. 14 with permission from AAAS, Copyright 2008.

Even without deriving scaling relations on paper as done for CH_x–NH_y, it is possible to simply try different scaling relations and evaluate them on the basis of their mean absolute error and/or R² value. Figure 1.5(b) shows a very good scaling relation for C₂H₂ and C₂H₄ with methyl (CH₃), and many more such examples exist. Hence, if access to an adsorption energy database for many different surface species is available, it is possible to bypass the bond order conservation idea and just rely on statistical tools. Using this approach, it should be possible to derive a suitable descriptor for any adsorbate, but the physics behind the relationship may not be obvious.

1.3.2 Scaling Relations for Transition States: The Brønsted–Evans–Polanyi Relationship

Transition states are often considered special because they have a very short lifetime. In terms of their adsorption properties, however, they are not different from any other surface species and obey the same physical principles. It is therefore not surprising that scaling laws for transition states exist in the same way as they do for other surface intermediates. The truly remarkable aspect is that transition state scaling for a large number of reactions is nearly universal and can be described by a single linear relationship!^43–45 However, before continuing the discussion of transition state scaling it is important to clarify several technical terms.

Transition state scaling (TSS) is often used interchangeably with the term Brønsted–Evans–Polanyi or BEP relationship.^46,47 Although both terms refer to the same concept, there is a distinction between them. Transition state scaling refers to a relation between the transition state energy E_TS and the energy of either the initial state (E_IS) or final state (E_FS) of a reaction. The transition state energy, E_TS, is not equal to the activation barrier E_a of a reaction, but it can be used to calculate E_a = E_TS − E_IS. A BEP relation, on the other hand, is a relation between the activation energy E_a and the heat of a reaction ΔE, which is defined as ΔE = E_FS − E_IS. All energy quantities are indicated on a schematic potential energy surface (PES) in Figure 1.6.

Figure 1.6

View large Download slide

A schematic potential energy surface for an arbitrary surface reaction. The nomenclature is E_ref = reference energy chosen to be zero, E_IS = initial state energy, E_TS = transition state energy, E_FS = final state energy, E_a = activation energy barrier, ΔE = energy change of the reaction.

A direct comparison between a TSS relation and a BEP relation for an exhaustive data set of (de)hydrogenation reactions over transition metals is shown in Figure 1.7.⁴⁵ Since all reactions were formulated as dissociation reactions, E_diss was used to denote the final state energy E_FS; and ΔE_diss represents the energy change of the surface reaction, ΔE. Although identical raw data have been used to generate the plots, the transition state scaling in Figure 1.7(a) appears to be much better than the BEP relation in Figure 1.7(b). This initial impression is somewhat deceptive and a comparison of the mean absolute error (MAE) of the two plots indicates that the quality, in terms of absolute errors, of both relations is almost identical. The difference lies in the different energy scales used on the x- and y-axes of the plots, and the BEP graph can be loosely interpreted as a zoomed-in version of the transition state scaling graph.

Figure 1.7

View large Download slide

Transition state scaling (a) and BEP relation (b) for (de-)hydrogenation reactions. Reproduced from ref. 45 with permission from the Royal Society of Chemistry.

To demonstrate that, for most reactions, both graphs contain essentially the same information, a mathematical conversion of one graph into the other will be derived. We may assume that the initial and final state of the reaction is described by linear scaling relations that depend on the descriptors E₁ and E₂.

E_{IS} = γ_{1} E_{1} + ξ_{1}

(1.11)

E_{FS} = γ_{2} E_{2} + ξ_{2}

(1.12)

Using transition state scaling, E_TS can be found as a function of E_FS

E_{TS} = α_{1} E_{FS} + β_{1}

(1.13)

The quantity that is of practical interest in this case is not E_TS directly, but the activation energy barrier E_a, which is obtained in TSS from

\begin{array}{l} E_{a}^{TSS} & = & E_{TS} - E_{IS} = α_{1} E_{FS} + β_{1} - E_{IS} = α_{1} (γ_{2} E_{2} + ξ_{2}) + β_{1} - (γ_{1} E_{1} + ξ_{1}) \\ E_{a}^{TSS} & = & α_{1} γ_{2} E_{2} - γ_{1} E_{1} + α_{1} ξ_{2} + β_{1} - ξ_{1} . \end{array}

(1.14)

The last three terms are constants and may be summed and replaced with β^TSS; which then yields the final expression

E_{a}^{TSS} = α_{1} γ_{2} E_{2} - γ_{1} E_{1} + β^{TSS}

(1.15)

Alternatively, E_a can be obtained from the BEP relation.

\begin{array}{l} E_{a}^{BEP} & = & α_{2} Δ E + β_{2} = α_{2} (E_{FS} - E_{IS}) + β_{2} = α_{2} (γ_{2} E_{2} + ξ_{2} - γ_{1} E_{1} + ξ_{1}) + β_{2} \\ E_{a}^{BEP} & = & α_{2} γ_{2} E_{2} - α_{2} γ_{1} E_{1} + α_{2} (ξ_{2} - ξ_{1}) + β_{2} \\ E_{a}^{BEP} & = & α_{2} γ_{2} E_{2} - α_{2} γ_{1} E_{1} + β^{BEP} \end{array}

(1.16)

Although and look similar, they are not identical if the descriptors E₁ and E₂ are considered to be independent variables. However, for the vast majority of reactions the descriptors are either identical, linearly correlated or at least exhibit a nearly linearly correlated behavior. Assuming a linear correlation and substituting

E_{2} = γ_{3} E_{1} + ξ_{3}

(1.17)

into and , E₂ can be eliminated to obtain

E_{a}^{TSS} = (α_{1} γ_{2} γ_{3} - γ_{1}) E_{1} + {\tilde{β}}^{TSS}

(1.18)

E_{a}^{BEP} = (α_{2} γ_{2} γ_{3} - α_{2} γ_{1}) E_{1} + {\tilde{β}}^{BEP}

(1.19)

where the constant terms have been included in

{\tilde{β}}^{TSS}

and

{\tilde{β}}^{BEP}

⁠. By equating

E_{a}^{TSS}

and

E_{a}^{BEP}

from and it can easily be seen that they are identical when

{\tilde{β}}^{TSS} = {\tilde{β}}^{BEP}

and

α_{2} = \frac{α_{1} γ_{2} γ_{3} - γ_{1}}{γ_{2} γ_{3} - γ_{1}}

⁠. This shows that, as long as the descriptors E₁ and E₂ used for the initial and final state, respectively, are linearly correlated, it does not matter whether a TSS or BEP relation is used.

Now that it has been demonstrated that the choice between TSS and BEP is mostly a personal preference, the author wants to explain his preference toward BEP over TSS relations. Assuming that a transition state can be characterized as initial state like, final state like, or anything in between. It is possible to define a similarity factor ω such that ω = 0 corresponds to a transition state that is initial state like, while ω = 1 corresponds to a final state-like transition state. The energy of the transition state, E_TS, would then scale with E_IS and E_FS weighted by the transition state similarity factor ω. With β as an arbitrary energy offset, the energy E_TS can then be expressed as

\begin{array}{l} E_{TS} & = & (1 - ω) E_{IS} + ω E_{FS} + β = ω (E_{FS} - E_{IS}) + E_{IS} + β \\ E_{TS} & = & ω Δ E + E_{IS} + β \end{array}

(1.20)

Rewriting for the calculation of an activation barrier leads to a BEP relation relating E_a with ΔE and the slope of the BEP relation, α, can be interpreted as the similarity factor, ω.

E_{a} = E_{TS} - E_{IS} = ω Δ E + β

(1.21)

In addition to the physical interpretation of the BEP slope, using a BEP relation has other practical advantages. For the calculation of rate and equilibrium constants E_a and ΔE are needed, which can be directly obtained from a BEP relation. In more advanced kinetic models, the effect of surface coverage on the stability of the surface intermediates may also be explicitly included, but including such an explicit coverage dependency on activation energies is not trivial.^48,49 Assume that DFT has been used to calculate E_a for a surface reaction at low coverage and that it is possible to predict the initial and final state energies for a high-coverage scenario. Performing additional DFT calculation for E_a at various higher surface coverages is prohibitively time consuming and an approximate method is desired to estimate the effect of coverage on E_a. In this case, the same idea used to derive ( can be applied): the constant β is eliminated by using a reference calculation, and finding E_a.

E_{a} = E_{a}^{ref} + ω (Δ E - Δ E^{ref})

(1.22)

Now that TSS and BEP relations can be fully appreciated, we discuss their most remarkable property, the fact that these relations are almost universally found for dissociation reactions.^44,45 This is referred to as the universality principle in heterogeneous catalysis and can be explained by the geometric similarities of the transition states with their respective dissociated states.⁴³ In Figure 1.7 it has already been shown that the universality principle holds for (de)hydrogenation reactions, and Figure 1.8(a) illustrates that the principle also applies to many other dissociation reactions. However, the MAE indicated in Figure 1.8(a) of 0.35 eV is larger than the typically accepted MAE (0.3 eV) and the individual TSS and/or BEP relations for smaller families of reactions listed in ref. 44 and 45 should be used if possible. For preliminary screening studies, however, it is perfectly acceptable to assume universal TSS and/or BEP behavior and refine the results later, as needed.

Figure 1.8

View large Download slide

(a) Universal TSS relation for C–C, C–O, C–N, N–O, N–N, and O–O dissociation reactions. Reproduced from ref. 44 with permission from Springer Nature, Copyright 2011. (b) Deviations from the “universal” TSS behavior are observed when the nature of the transition state switches between initial and final states alike. Reproduced from ref. 50 with permission from AIP Publishing, Copyright 2011.

At this point it should be stressed that linear TSS and BEP relations exist for many reactions on a large range of catalyst surfaces, but they are not as “universal” as they have been introduced. Even for simple di-atomic dissociation reactions (N₂, NO, O₂) on metal-substituted La-perovskites a significant deviation from the universal TSS–BEP line has been observed. The deviation can clearly be seen in Figure 1.8(b) as the plot moves to more reactive surfaces (more negative ΔE_diss) and can be explained by a continuous change from a late transition state (final state like) to an early transition state on more reactive surfaces. The same phenomenon can be seen on highly reactive metallic surfaces, but all transition metals that are typically used as catalysts fall into the linear region of the TSS–BEP line and the relation holds quite well in this range. Nevertheless, care must be taken to choose the appropriate TSS and/or BEP relation for each reaction and even more attention must be paid when these relations are used to extrapolate beyond the parameter range for which they were obtained.

In summary, scaling and BEP relations are very versatile, broadly applicable, and allow for a significant reduction of the number of energetic parameters to a small set of catalytic descriptors. Given the typical error bars of these relations the results are not meant to be taken quantitatively, but if the scaling and BEP relations are carefully chosen and correctly applied, then catalytic activity and selectivity trends over a broad range of materials can be obtained very well.

1.4 The Sabatier Principle and the Volcano Curve

For most things in life, there is an optimal concentration, amount or frequency in which they should occur; catalysis is no exception. The Nobel Prize-winning French chemist Paul Sabatier qualitatively described the adsorption properties of the optimal catalyst for a reaction as not too weak, such that the reactants don’t bind, and not too strong, which would lead to surface poisoning.⁵¹ This statement is known today as the Sabatier principle of catalysis. The principle is easily visualized in volcano-shaped curves that show the catalytic activity as a function of a catalytic descriptor. An example of such a volcano curve was already introduced in Figure 1.1(a) and a schematic representation is shown in Figure 1.9. The top of the volcano corresponds to the highest activity and the descriptor value that is “just right”. Finding the optimal value of the descriptor(s) is the central and by far the most critical point in any computational catalyst screening study. Once the optimal value is known, it requires only a simple database search to screen for materials with the desired property. The other features of the volcano, including the activity at the volcano peak and the slopes at which the activity decreases along both sides of the maximum, are of secondary importance. Since the focus is primarily on the descriptor value at the volcano peak, plus or minus a few tenths of an eV, it is not always necessary to develop a full-blown microkinetic model. Instead, a shortcut method, the Sabatier analysis, can be used to predict upper bounds for the reaction rate. Sabatier analyses have been very successfully used to identify catalytic trends for CO oxidation and NO decomposition.^52,53 As shown in Figure 1.9, the Sabatier rate reproduces the shape of the volcano curve obtained from a microkinetic model well enough to provide an estimate of the optimal descriptor value, and it is significantly easier and faster to implement.

Figure 1.9

View large Download slide

Comparison of the volcano curves obtained from a Sabatier analysis and full microkinetic models at different approaches to equilibrium γ. Reproduced from ref. 21 with permission from Elsevier, Copyright 2004.

1.4.1 Sabatier Analysis

The Sabatier analysis for surface reactions was first introduced by Bligaard et al.²¹ and uses a Sabatier rate that is obtained by following this fairly simple recipe:

Write the rate expressions for all reaction steps in the forward direction as a function of the approach to equilibrium $γ_{i} = \prod_{n} a_{n}^{ν_{n}} / K_{i}$ ⁠, where K_i is the equilibrium constant for step i, a_n is the activity of species n, and ν_n is the corresponding stoichiometric coefficient.
For each reaction step in the model, assume optimal surface coverages for the forward reaction.
From the overall approach to equilibrium, calculate the approach to equilibrium of step i under the assumption that all other steps are quasi-equilibrated.
Use TSS and/or BEP relations to calculate the forward rate constants k_i and the rate. This rate represents an upper bound on the actual rate of step i and will be denoted as $r_{i}^{\max}$ ⁠.
The Sabatier rate is the slowest rate of all maximum rates in the mechanism: $r_{Sabatier} = \min {r_{i}^{\max}}$ ⁠.

The procedure can be illustrated using a generic reaction mechanism with two steps as an example. The two reaction steps are a dissociation step, step (1), and a reactive desorption step, step (2):

$A_{2} + 2 ⁎ ⇌ 2 A^{⁎}$
$A^{⁎} + B ⇌ AB + ⁎$

where the asterisk ‘*’ denotes a free surface site and A* is the adsorbed species A. With θ denoting the surface coverage and γ the approach to equilibrium, the reaction rates in the forward direction are

r_{1} = 2 k_{1} P_{A_{2}} {θ_{⁎}}^{2} (1 - γ_{1})

(1.23)

r_{2} = k_{2} P_{B} θ_{A} (1 - γ_{2})

(1.24)

The optimal surface coverage for step (1) is a completely empty surface (θ_* = 1) and for step (2) the surface should be fully covered with A (θ_A = 1). The overall approach to equilibrium γ is closely related to the equilibrium constant and is defined as

γ = \frac{{P_{AB}}^{2}}{K_{eq} P_{A_{2}} {P_{B}}^{2}} with γ = γ_{1} {γ_{2}}^{2}

(1.25)

The allowed range for γ is 0 ≤ γ ≤ 1, where γ = 1 indicates a quasi-equilibrated reaction. From γ₁ and γ₂ can be estimated under the assumption that the other step is quasi-equilibrated and results in

γ_{1} = γ and γ_{2} = \sqrt{γ}

(1.26)

Assuming that k₁ and k₂ are known, the upper bounds for r₁ and r₂ can be written down and the Sabatier rate r_Sabatier found.

r_{1}^{\max} = 2 k_{1} P_{A_{2}} (1 - γ)

(1.27)

r_{2}^{\max} = k_{2} P_{B} (1 - \sqrt{γ})

(1.28)

r_{Sabatier} = \min {r_{1}^{\max}, r_{2}^{\max}}

(1.29)

It is straightforward to apply the same analysis to more complicated reaction mechanisms with a larger number of steps. The resulting volcano would have more than two sides, but the maximum and the corresponding descriptor value can still be found just as easily as in the case presented here. For reaction conditions where the overall approach to equilibrium approaches unity (γ → 1), the simple Sabatier analysis can significantly deviate from a full microkinetic model solution, but at low conversions far away from equilibrium where γ → 0, the agreement is generally good.

For the sake of completeness, the Sabatier–Gibbs analysis is also briefly mentioned here, although it is not particularly useful for general problems. The Sabatier–Gibbs analysis can only be applied to serial reaction mechanisms that are very simple, and it improves on the unrealistic assumption of optimal surface coverages for each step used in the Sabatier analysis. A stricter bound on the coverages can be imposed by considering thermodynamic limits in terms of the approach to equilibrium, but the procedure is somewhat cumbersome and in general it is better to use a microkinetic model if the simple Sabatier analysis fails. For additional information the reader is referred to ref. 52 and 54.

1.5 Sabatier Analysis in Practice

1.5.1 First Example: Ammonia Synthesis

Now it is finally time to combine all the concepts introduced so far and perform a first real catalyst design. The reaction used for the design study is one of the most important heterogeneously catalyzed reactions, namely ammonia (NH₃) synthesis from nitrogen (N₂) and hydrogen (H₂).

N_{2} + 3 H_{2} ⇌ 2 {NH}_{3} Δ H^{°} = - 45.9 kJ per mol ({NH}_{3})

This reaction is performed industrially using the Haber–Bosch process over iron (Fe) or ruthenium (Ru) catalysts at high temperatures (>400 °C) and high pressures (>100 bar). Using state-of-the-art surface science⁵⁵ and theoretical techniques,⁵⁶ the detailed reaction mechanism has been elucidated and a complete PES on terrace and step sites of Ru is shown in Figure 1.10.

Figure 1.10

Potential energy surface for ammonia synthesis on terrace (dashed line) and stepped (unbroken line) sites on Ru. Reproduced from ref. 56 with the permission from AAAS, Copyright 2005.

View large Download slide

Upon inspection of the PES for NH₃ synthesis the mechanism can be simplified to two main reaction steps, the dissociative adsorption of N₂ followed by its hydrogenation to NH₃. An extremely simplified reaction mechanism can therefore be written as

$N_{2} + 2 ⁎ ⇌ 2 N^{⁎} (dissociative adsorption)$
$N^{⁎} + \frac{3}{2} H_{2} ⇌ {NH}_{3} +^{⁎} (recombinative desorption)$

This mechanism neglects a number of elementary steps and does not include the possibility of adsorbed H atoms or NH_x species, but it is sufficient to illustrate the Sabatier principle. A catalyst that associates strongly with N will easily dissociate N₂, but the desorption step will be very difficult (note that in the desorption step only the energy of N* is affected by the catalyst because H₂ and NH₃ are gas phase species). In contrast, a catalyst interacting weakly with N allows for easy desorption, but the barrier for N₂ dissociation will be large. Hence, there is an optimal binding energy of N*, E_N, at which the dissociative adsorption of N₂ and NH₃ desorption are optimized simultaneously. The objective now is to find an estimate of E_N, by performing a Sabatier analysis for typical NH₃ synthesis conditions.

The simplified NH₃ synthesis mechanism is almost identical to the generic mechanism used in the introduction of the Sabatier analysis in Section 1.4. If A₂ = N₂ and B = 3/2H₂, are substituted into the equation it is possible to use the previously obtained results and apply them to the current problem. This gives the following two maximum rates for each step, where the stoichiometric factor of 3/2 for H₂ is included as the exponent to

P_{H_{2}}

in the expression for

r_{2}^{\max}

⁠.

r_{1}^{\max} = 2 k_{1} P_{N_{2}} (1 - γ)

(1.30)

r_{2}^{\max} = k_{2} {P_{H_{2}}}^{\frac{3}{2}} (1 - \sqrt{γ})

(1.31)

For the reaction conditions it may be assumed that there is a stoichiometric feed composition (H₂ : N₂ = 3 : 1), the total pressure is P = 100 bar and the temperature is T = 673 K. Using this information, the partial pressures of hydrogen and nitrogen can be calculated as

P_{H_{2}} = 75 bar

and

P_{N_{2}} = 25 bar

⁠. Furthermore, it is assumed that the reaction is run at the limit of very low conversion where γ → 0. This leaves only the two unknown rate constants k₁ and k₂. The rate constants, k_i, can be obtained from the standard Arrhenius expression, which depends on the activation energy, E_a,i, and the pre-exponential factor, v_i, for each step.

k_{i} = v_{i} \times \exp (\frac{- E_{a, i}}{k_{B} T})

(1.32)

The pre-exponential factor v_i can be calculated using transition state theory from the entropy change between the initial state and the transition state of the reaction,

Δ S_{i}^{‡}

⁠.

v_{i} = \frac{k_{B} T}{h} \times \exp (\frac{Δ S_{i}^{‡}}{k_{B}})

(1.33)

Next, the activation barrier estimates are obtained from published scaling relations, beginning with N₂ dissociation. The dissociation of diatomic molecules was the first group of heterogeneously catalyzed reactions for which a universal BEP relation was identified and the values of α = 0.87 and β = 1.34 eV were reported for stepped surfaces.⁴³

E_{a, 1} = 0.87 \times Δ E_{1} + 1.34 eV

(1.34)

The second step in the simplified mechanism lumps several hydrogenation steps into one step, therefore it seems appropriate to use the averaged BEP parameters for NH₃, NH₂, and NH dehydrogenation reactions found in Table 2 of ref. 44. Be careful here! The BEP parameters α = 0.61 and β = 1.43 eV are derived for the reverse reaction of step (2), so it is important to make sure to use them accordingly. The estimated reverse activation energy

E_{a, 2}^{rev}

is then given by

E_{a, 2}^{rev} = 0.61 \times (- Δ E_{2}) + 1.43 eV

(1.35)

and the forward barrier is found by using the fact that

Δ E_{2} = E_{a, 2}^{fwd} - E_{a, 2}^{rev}

⁠. The energy change of step (2), ΔE₂, depends linearly on ΔE₁ owing to the thermodynamic constraint that

\frac{1}{2} Δ E_{1} + Δ E_{2} = Δ E_{r}

⁠, where ΔE_r is the overall energy change of the reaction. So far, the electronic ground-state energy, E, at T = 0 K has been used exclusively in all equations, which is the primary output of DFT calculations, but from thermodynamics it is known that free energy changes are the correct quantities to use in the calculation of rate and equilibrium constants. It is assumed for now that the DFT-calculated electronic energy is approximately equal to the enthalpy and that the entropy change can be calculated separately. This can be justified by the fact that the electronic energy changes are significantly larger than all other contributions to the enthalpy. Thus, for the present problem it can be said that ΔE_r ≈ ΔH° = −45.9 kJ mol⁻¹ = −0.48 eV. Putting this all together yields the desired BEP relation for E_a,2.

\begin{array}{l} E_{a, 2} & = & Δ E_{2} + E_{a, 2}^{rev} = Δ E_{2} + 0.61 \times (- Δ E_{2}) + 1.43 eV \\ = & 0.39 \times Δ E_{2} + 1.43 eV = 0.39 \times (Δ E_{r} - \frac{1}{2} Δ E_{1}) + 1.43 eV \\ = & - 0.195 \times Δ E_{1} + 1.24 eV \end{array}

(1.36)

By choosing these particular BEP relations to estimate the activation barriers for steps (1) and (2), the reactivity descriptor for this reaction has been implicitly determined. Both barriers depend only on ΔE₁, which is the dissociative chemisorption energy of nitrogen, ΔE_N.

The last piece of missing information comprises the entropy changes between initial and transition states, which are needed to calculate the pre-exponential factors. The entropy values of adsorbed species and transition states depend only weakly on the transition metal surface. Therefore, it may be assumed that the entropy is independent of the catalyst. For a rough approximation, it can be further assumed that an adsorbed species or transition state has lost all degrees of freedom (it is completely locked in place on the surface), and its entropy is therefore zero. The gas-phase entropies of N₂ and H₂ at standard state are S°(N₂) = 192.77 J mol⁻¹ K⁻¹ and S°(H₂) 130.68 J mol⁻¹ K⁻¹, respectively, and the entropy changes between the initial states and the transition state can be estimated as $Δ S_{1}^{‡} = - 192.77 J {mol}^{- 1} K^{- 1}$ and $Δ S_{2}^{‡} = - 130.68 J {mol}^{- 1} K^{- 1}$ ⁠. For the first step it is assumed that N₂ loses all its gas-phase entropy upon entering the transition state and in the second step, the entropy loss of one H₂ gas-phase molecule is considered when it reacts with N* on the surface. Why was the entropy loss of only one H₂ molecule considered and not 3/2H₂ as the stoichiometric coefficient in the lumped equation would indicate? Well, the fact that step (2) is a lumped step and not an elementary reaction step makes it almost impossible to guess the correct entropy change, but it is known that it is a sequence of hydrogen adsorption and hydrogenation steps and for each H₂ molecule adsorbing onto the surface −130.68 J mol⁻¹ K⁻¹ of entropy are lost. When mechanisms with lumped reaction steps are analyzed, problems with estimating energy barriers and entropy changes will frequently occur, so it is highly recommended to formulate mechanisms that are comprised of only elementary reaction steps. Still, for this case, from eqn (1.33) and the entropy assumptions outlined above, the pre-exponential factors are calculated as v₁ = 1.20 × 10³ s⁻¹ and v₂ = 2.10 × 10⁶ s⁻¹.

At this point the reader may want to open their favorite math and graphing software (Excel would suffice), enter eqn (1.30)–(1.32), (1.34), and (1.36), and calculate the maximum N₂ dissociation rate $r_{1}^{\max}$ and maximum desorption rate $r_{2}^{\max}$ as a function of ΔE₁ = ΔE_N and plot the results. The graph should look like the one given in Figure 1.11. The optimal catalyst for ammonia synthesis is found for a descriptor value of ΔE_N of approximately −0.95 eV and has a predicted turn over frequency (TOF) of 8.7 s⁻¹. It is usually a lucky coincidence if the estimated TOF is close to an experimentally measured value in such a crude Sabatier analysis, but the location of the maximum with respect to the descriptor value in the volcano curve is typically reproduced within acceptable error bars. With a target value for the catalytic descriptor, it is possible to finally proceed to the last step of the example problem and screen for the best catalyst. For this purpose, the “database” of dissociative chemisorption energies, ΔE_N, of several transition metals listed in Table 1.1 is examined. It is found that well-known NH₃ synthesis catalysts Fe and Ru are closest to the optimal value of −0.95 eV. Considering the number of assumptions and simplifications that were required to perform the Sabatier analysis, it seems surprising that the two most active catalysts for ammonia synthesis were correctly identified. This certainly emphasizes the strength and robustness of the method. The analysis is helped by the fact that ΔE_N spans a range from ca. −4 to +6 eV and neighboring transition metals are often separated by 0.5 eV or more.

Figure 1.11

View large Download slide

Sabatier volcano for NH₃ synthesis.

Table 1.1

—

−1.27

−0.38

−0.10

2.88

−2.76

—

−0.84

−0.70

1.78

5.86

−4.33

—

−0.59

1.37

5.89

View Large

In summary, this example of computational screening for an NH₃ synthesis catalyst has demonstrated how BEP relations are applied to reduce the number of parameters to a single descriptor, ΔE_N, and how to perform a Sabatier analysis. As a result, the Sabatier volcano shown in Figure 1.12 was obtained from which two suitable catalysts for the reaction can be identified (Fe and Ru). The example also introduced approximations that can be applied to obtain a simple estimate for surface entropies, which will be improved on in Section 1.6.2.

Figure 1.12

Sabatier volcano for NH3 synthesis with the positions of transition metals labeled on the graph.

View large Download slide

Sabatier volcano for NH₃ synthesis with the positions of transition metals labeled on the graph.

1.5.2 Second Example: CO Oxidation

The Sabatier analysis for NH₃ synthesis presented in the previous section captures all the important ideas needed for a computational catalyst design project, but it is very idealized. A somewhat more complicated example is CO oxidation, for which a volcano surface can be derived as a function of two independent descriptors. This analysis was originally completed by Falsig et al. and resulted in the prediction of the most active catalysts for high (Pt, Pd) and low (Au) temperature CO oxidation.⁵³ Here, only the high temperature CO oxidation with the reaction conditions T = 600 K,

P_{O_{2}} = 0.33 bar

⁠, P_CO = 0.67 bar is discussed. The low temperature reaction can be studied analogously. For the CO oxidation reaction, the analyses start with the following reaction mechanism.

$CO + ⁎ ⇌ {CO}^{⁎} (quasi - equilibrated)$
$O_{2} + ⁎ ⇄ {O_{2}}^{⁎} (quasi - equilibrated)$
${O_{2}}^{⁎} + ⁎ ⇌ 2 O^{⁎}$
${CO}^{⁎} + O^{⁎} ⇌ {CO}_{2} + 2 ⁎$

Following Falsig et al., it is assumed that the molecular adsorption of CO and O₂ is fast and can be considered as quasi-equilibrated. The two competing steps in the reaction are then the dissociation of O₂* and the reaction between CO* and O* to form CO₂. Assuming very low conversions γ → 0 the maximum rates of steps (3) and (4) are:

r_{3}^{\max} = 2 k_{3} θ_{O_{2}} θ_{⁎}

(1.37)

r_{4}^{\max} = k_{4} θ_{CO} θ_{O}

(1.38)

The assumed optimal coverages that maximize

r_{3}^{\max}

are

θ_{O_{2}} = θ_{⁎} = 0.5

⁠, and

r_{4}^{\max}

is maximized for θ_CO = θ_O = 0.5, leading to the final expression of the Sabatier rate.

r_{Sabatier} = \min {0.5 \cdot k_{3}, 0.25 \cdot k_{4}}

(1.39)

Reaction steps (3) and (4) in the mechanism for CO oxidation are both surface reactions, i.e. all reactants and the transition state are adsorbed species, and for this type of reactions the pre-exponential factors are typically on the order of v = 10¹³ s⁻¹. The underlying assumption is that the entropy of the transition state is similar to the entropy of the initial state (ΔS^‡ = 0) and reduces to

v = \frac{k_{B} T}{h} \approx 10^{13} s^{- 1}

at room temperature. The activation barriers for steps (3) and (4) can be found from the (transition state) scaling relations identified by Falsig et al. for face centered cubic [fcc(111)] surfaces.⁵³

Δ E_{TS, 3} = 1.39 \times Δ E_{O} + 1.56 eV

(1.40)

Δ E_{TS, 4} = 0.7 \times (Δ E_{CO} + Δ E_{O}) + 0.02 eV

(1.41)

Δ E_{O_{2}} = 0.89 \times Δ E_{O} + 0.17 eV

(1.42)

The TSS equations can be rewritten to yield activation energy barriers.

Δ E_{a, 3} = Δ E_{TS, 3} - Δ E_{O_{2}} = 0.5 \times Δ E_{O} + 1.39 eV

(1.43)

Δ E_{a, 4} = Δ E_{TS, 4} - Δ E_{CO} - Δ E_{O} = - 0.3 \times (Δ E_{O} + Δ E_{CO}) + 0.02 eV

(1.44)

The Sabatier rate for CO oxidation in can now be obtained from the knowledge of ΔE_CO and ΔE_O, which have naturally evolved as the two descriptors needed to describe the CO oxidation activity. Later in this chapter, a full microkinetic model is used to construct a two-dimensional volcano for CO oxidation that depends on both descriptors, but a one-dimensional volcano slice as shown in Figure 1.13 is easier to follow. This volcano slice was prepared under the artificial assumption that all transition metals have a constant binding energy of CO, which was set to ΔE_CO = −1.20 eV. This value is within ±0.3 eV of the actual CO binding energies of Ru, Rh, Ni, Pt, and Pd. In contrast, Cu, Ag, and Au bind CO much more weakly, and their placement on the volcano in Figure 1.13 must be considered extremely approximate. CO oxidation rates on metals to the left of the maximum are limited by step (4), the removal of strongly adsorbed O*, while on the right side of the maximum the limiting step is O₂ dissociation, step (3). The optimal catalyst is the best compromise between all possible rate-determining steps, such that, at the maximum, no single step can be identified as rate-determining. It should also be mentioned that rate-limiting steps can change depending on the catalyst and reaction conditions and they should not be considered as a constant in any given reaction mechanism. A more detailed discussion of rate-limiting steps is presented later in this chapter. The important observation here is that Pt and Pd are located closest to the top of the volcano, which is exactly what was expected. Given the similarly correct identification of the optimum catalysts for NH₃ synthesis in the earlier example, the success of a Sabatier analysis can no longer be considered a coincidence!

Figure 1.13

View large Download slide

A slice through the CO oxidation volcano at T = 600 K for a constant value of ΔE_CO = −1.20 eV.

1.6 Notes on Microkinetic Modeling

Although the Sabatier analysis can be applied to reaction mechanisms with any number of elementary steps and even mechanisms with parallel reaction pathways, it is an approximation. More accurate computational screening for complex reaction mechanisms can be performed by implementing a full microkinetic model, which contains information about the dominant reaction mechanism and simultaneously predicts the activity, selectivity and surface coverages. As input to this microkinetic model only the catalytic descriptors are needed and all the missing energy information must be estimated from BEP and scaling relations, or assumed. But before continuing with more advanced computational screening examples, it is necessary to introduce some very basic concepts needed for the development of microkinetic models in the context of computational catalyst screening. Gokhale et al.,⁵⁷ Stoltze and Nørskov,⁵⁸ and the standard reference The Microkinetics of Heterogeneous Catalysis by Dumesic et al.⁵⁹ provide useful information for microkinetic modeling.

The microkinetic modeling procedure for catalyst screening will be illustrated for the CO oxidation example, which is the prototype reaction in heterogeneous catalysis.⁶⁰ The same four-step mechanism used in the Sabatier example will be used here and the net rates of the elementary steps are given by:

r_{1} = k_{1}^{+} P_{CO} θ_{⁎} - k_{1}^{-} θ_{CO}

(1.45)

r_{2} = k_{2}^{+} P_{O_{2}} θ_{⁎} - k_{2}^{-} θ_{O_{2}}

(1.46)

r_{3} = k_{3}^{+} θ_{O_{2}} θ_{⁎} - k_{3}^{-} {θ_{O}}^{2}

(1.47)

r_{4} = k_{4}^{+} θ_{CO} θ_{O} - k_{4}^{-} P_{{CO}_{2}} {θ_{⁎}}^{2}

(1.48)

The change in coverage of reaction intermediates is determined by the net rate at which they are produced or consumed in each reaction step.

\frac{d θ_{CO}}{d t} = r_{1} - r_{4}

(1.49)

\frac{d θ_{O_{2}}}{d t} = r_{2} - r_{3}

(1.50)

\frac{d θ_{O}}{d t} = 2 r_{3} - r_{4}

(1.51)

The fraction of empty sites on the surface can at any point of time be calculated from the site balance.

θ_{⁎} = 1 - θ_{O_{2}} - θ_{O} - θ_{CO}

(1.52)

represent a system of ordinary differential equations (ODE), but the analyses will be restricted to situations where the reaction runs under steady-state conditions, which means that the derivatives of the surface coverages with respect to time in are all zero.

\frac{d θ_{i}}{d t} = 0 for i = O_{2}, O, CO

(1.53)

This converts the system of ODEs into a system of non-linear algebraic equations that can be solved with standard root-finding methods. The steady-state assumption is not strictly necessary when solving a microkinetic model numerically, but the transient start-up and shut-down behavior is typically short in comparison to steady-state operation.

With eqn (1.52) and (1.53) there are four equations for the four unknown coverages $θ_{O_{2}}, θ_{O}, θ_{CO}$ ⁠, and θ*, and the system of nonlinear algebraic equations may be solved numerically. With currently available CPU speeds numerical solutions to microkinetic models for catalyst screening studies are generally preferred because they avoid the need to make any additional assumptions regarding the mechanism.

Analytical solutions are also valuable, but they almost always require additional assumptions and involve cumbersome pencil and paper math. A very common assumption is that molecular adsorption and desorption steps are fast in comparison with surface reactions and may be assumed to be quasi-equilibrated. This assumption can be made for steps (1) and (2) of the mechanism and the surface coverages of CO* and O₂* obtained as a function of the fraction of empty sites.

r_{1} \approx 0 \Rightarrow k_{1}^{+} P_{CO} θ_{⁎} = k_{1}^{-} θ_{CO} \Rightarrow θ_{CO} = K_{1} P_{CO} θ_{⁎}

(1.54)

r_{2} \approx 0 \Rightarrow k_{2}^{+} P_{O_{2}} θ_{⁎} = k_{2}^{-} θ_{O_{2}} \Rightarrow θ_{O_{2}} = K_{2} P_{O_{2}} θ_{⁎}

(1.55)

There is no reaction that can be assumed to be quasi-equilibrated to derive an equation for the coverage of O*, but from the steady-state assumption and the following is obtained.

0 = \frac{d θ_{O}}{d t} = 2 r_{3} - r_{4} = 2 (k_{3}^{+} θ_{O_{2}} θ_{⁎} - k_{3}^{-} {θ_{O}}^{2}) - k_{4}^{+} θ_{CO} θ_{O} + k_{4}^{-} P_{{CO}_{2}} {θ_{⁎}}^{2}

(1.56)

which can be solved for θ_O by substituting

θ_{O_{2}}

and

θ_{CO}

from and .

\begin{array}{l} θ_{O} & = & \frac{- k_{4}^{+} K_{1} P_{CO} + \sqrt{{(k_{4}^{+} K_{1} P_{CO})}^{2} - 8 k_{3}^{-} (k_{4}^{-} P_{{CO}_{2}} - 2 k_{3}^{+} K_{2} P_{O_{2}})}}{4 k_{3}^{-}} θ_{⁎} \\ = & W θ_{⁎} \end{array}

(1.57)

In W is introduced to denote the long fraction. For an analytical solution the site balance is used to solve for

θ_{⁎}

⁠.

\begin{array}{l} θ_{⁎} & = & 1 - K_{2} P_{O_{2}} θ_{⁎} - W θ_{⁎} - K_{1} P_{CO} θ_{⁎} \Rightarrow θ_{⁎} \\ = & \frac{1}{1 + W + K_{1} P_{CO} + K_{2} P_{O_{2}}} \end{array}

(1.58)

can be solved for any reaction condition if the equilibrium constants, K_i, for steps (1) and (2), and the forward and reverse rate constants,

k_{i}^{+}

and

k_{i}^{-}

⁠, for steps (3) and (4) are known. Note that it is not necessary to know the rate constants of the quasi-equilibrated steps (1) and (2) and the equilibrium constants are sufficient. This analytical solution will be used later and compared with a full numerical solution of the CO oxidation model.

1.6.1 Numerical Solution Strategies

Most people will say that solving a microkinetic model numerically is easy and does not require any special solution strategies: one can simply use a standard ODE solver or root-finding method. This may be the case for microkinetic models that are parameterized for a working catalyst and predict reasonable reaction rates, but when a computational screening study is being performed the objective is to explore a large parameter space for the descriptors, which also includes unusual parameter combinations. Specifically, at the edges far away from the volcano maximum, the rate constants in the microkinetic model span many orders of magnitude (40–50 orders is not unusual) and numerical solvers will not always converge for such extremely stiff systems of ODEs. Although there is no general procedure to tackle this problem there are some guidelines.

In a computational screening study, the rate and equilibrium constants calculated using scaling and BEP relations are functions of the chosen descriptors. A volcano plot could be constructed by systematically sweeping through the parameter space for the descriptors. Root-finding methods can quickly solve the steady-state equations, but the rate equations in microkinetic models are highly non-linear and convergence is only achieved when good initial guesses are provided. Hence, it is a good idea to use the solution of the microkinetic model for one descriptor set as the initial guess for the neighboring descriptor set. A possible solution path on a grid of two descriptors is suggested in Figure 1.14. Provided that the grid points in the descriptor space are close to each other, the solutions should be sufficiently similar to obtain converged results.

Figure 1.14

View large Download slide

Suggested looping technique for the exploration of a two-dimensional descriptor space.

An alternative approach is to integrate the system of ODEs, which is computationally much more intensive, but does not depend as strongly on the initial guesses, as long as they make physical sense, i.e. no negative coverages or violation of the site balance. In order to find a steady-state solution the ODEs need to be integrated over a sufficient amount of time to ensure that steady-state has been reached. It is a good practice to verify steady-state convergence by plotting the transient behavior of the surface coverages and ensuring that they do not continue to change over time.

Another important aspect to consider is the stiffness of the ODE system, a characteristic caused by the large differences in rate constants. Most ODE solvers, even those specialized for stiff systems, will fail in those cases. If the simpler root-finding strategy also fails, it is possible to reduce the stiffness of the ODE system by artificially slowing down some of the very fast, quasi-equilibrated, reaction steps. But be cautious! If these steps are slowed down too much, they may become kinetically relevant and affect the overall rate or the reaction mechanism. This strategy is not recommended for ordinary problems, but may be attempted with the appropriate care if everything else fails.

1.6.2 Entropy and Enthalpy Corrections

In the previous examples only electronic energy changes were considered and the entropy was approximated as all or nothing. In essence, it was assumed that gas-phase species have 100% of their standard state entropy and surface species possess no entropy at all. These assumptions can certainly be improved on, and in order to construct thermodynamically consistent microkinetic models this is not just optional, but absolutely necessary. Entropy and enthalpy corrections for surface species can be calculated using statistical thermodynamics from knowledge of the vibrational frequencies, and the translational and rotational degrees of freedom (DOF). In contrast to gas-phase molecules, adsorbates cannot freely rotate and move across the surface, but the translational and rotational DOF are frustrated within the potential energy well imposed by the surface. In the harmonic limit the frustrated translational and rotational DOF can conveniently be described as vibrational modes, which in turn means that any surface adsorbate will have 3N vibrational DOFs that are all treated equally.

The temperature dependent surface enthalpy, H(T), of an adsorbate is composed of the electronic energy, E_elec, the zero-point energy, E_ZPE, and a heat capacity correction to account for the difference between the enthalpy at 0 K and at the actual temperature, T.

H (T) = E_{elec} + E_{ZPE} + \int_{0 K}^{T} C_{V} d \tilde{T}

(1.59)

The electronic energy component, E_elec, is directly obtained from either DFT calculations or from scaling relations based on DFT energies. The remaining components can be calculated in the harmonic approximation.⁶¹

\begin{matrix} E_{ZPE} & = & \frac{1}{2} \sum_{i = 0}^{# DOF} h v_{i} \\ \int_{0 K}^{T} C_{v} d \tilde{T} & = & \sum_{i = 0}^{# DOF} \frac{h v_{i}}{e^{h v_{i}} / k_{B} T - 1} \\ S (T) & = & k_{B} \sum_{i = 0}^{# DOF} [\frac{h v_{i}}{k_{B} T (e^{h v_{i} / k_{B} T} - 1)} - \ln (1 - e^{- h v_{i} / k_{B} T})] \end{matrix}

(1.60)

The complete set of vibrational modes for an adsorbate is not easily obtained and usually it is necessary to perform an additional vibrational analysis using DFT. In this vibrational analysis it is important to pay attention to the low frequency modes, corresponding to the frustrated translation, for two reasons: (i) DFT is notoriously inaccurate in the calculation of low frequency modes, and (ii) low frequency modes contribute the most to the surface entropy. It is advised to check the entropy contributions of the translational modes separately, and to limit them to the entropy of a two-dimensional surface gas, which can be derived from the partition function of a particle that can move freely across the surface.

S_{trans}^{2 D} = k_{B} [\ln (\frac{2 π m k_{B} T}{C_{S} h^{2}}) + 2]

(1.61)

In this equation m denotes the mass of the particle and C_S is the number of sites per unit area, which is roughly 10¹⁵ sites cm⁻² for transition metal surfaces. Experience has shown that the entropy and enthalpy corrections do not vary appreciably across different transition metal surfaces and for computational screening it is commonly assumed that these contributions are constant.

1.6.3 Microkinetic Model Analysis

Microkinetic models can provide much more information than “just” surface coverages and reaction rates of individual steps. It is always a good idea to perform a sensitivity analysis on the model input parameters, to obtain additional valuable information. First, more attention should to be paid to parameters to which the model shows a high sensitivity, and these parameters need to be carefully measured or calculated. Non-sensitive parameters, on the other hand, may be roughly estimated or even guessed. Second, the sensitivity of a microkinetic model towards certain parameters can be interpreted as a measure of the degree to which a step is rate controlling (degree of rate control). This information can provide direction for useful catalyst modifications which could optimize performance (degree of catalyst control).

1.6.3.1 Degree of Rate Control

The degree of rate control, X_RC, was originally proposed by Campbell as a quantitative measure to identify rate-controlling steps.^62,63 It was later generalized to quantify also the effects of adsorbate binding (poisoning) on the reaction rate.⁶⁴ The degree of rate control of step i is formally defined as the normalized partial derivative of the overall rate with respect to the rate constant, k_i, while keeping the equilibrium constant, K_i, and the rate constants, k_j, of all other steps constant.

X_{RC, i} = \frac{k_{i}}{r} {(\frac{\partial r}{\partial k_{i}})}_{K_{i}, k_{j}} = {(\frac{\partial \ln r}{\partial \ln k_{i}})}_{K_{i}, k_{j}}

(1.62)

An X_RC,i value of zero indicates a quasi-equilibrated step, whose rate constant has no effect on the overall rate, while X_RC,i = 1 indicates a single rate-controlling step in a serial reaction mechanism. For serial reaction mechanisms with only one product, it has also been found that

\sum_{i = 0}^{# steps} X_{RC, i} = 1

⁠, but in reaction mechanisms with parallel pathways or more than one possible product an analogous rule has not been identified. In the context of microkinetic modeling, is not easily applicable, but it can be rewritten into a finite-difference equation that can be readily evaluated in a microkinetic model.

X_{RC, i} = \frac{k_{i}}{r} {(\frac{\partial r}{\partial k_{i}})}_{K_{i}, k_{j}} \approx \frac{k_{i}^{0}}{r^{0}} \times \frac{r - r^{0}}{k_{i} - k_{i}^{0}} = \frac{r - r^{0}}{x \times r^{0}} with k_{i} = (1 + x) k_{i}^{0}

(1.63)

In practice X_RC,i is then calculated by changing the forward and reverse rate constants of step i simultaneously by a small amount, typically 10% (x = 0.1), and monitoring the change of the overall reaction rate with respect to the reference rate, r⁰, for the original values of k_i. Never forget to apply this change to both rate constants, otherwise K_i and the overall equilibrium constant, K_eq, will be altered as well!

Fundamentally, the definition of the degree of rate control in implies a change in the Gibbs free energy of the transition state of step i, as indicated for the transition state of CH₄ activation in Figure 1.15, while holding the energy of all other states in the potential energy surface constant. In more mathematical terms, this interpretation can be derived from the Arrhenius expression, which can be written as

\ln (\frac{k_{i}}{c}) = \frac{- G_{a, i}}{k_{B} T}

(1.64)

where c is a constant (⁠

c = \frac{k_{B} T}{h}

in transition state theory) and G_a,i is the free energy of activation for step i. Upon substituting into the definition of X_RC,i in degree of rate control is obtained in terms of the free energy of the transition state of step i, (G_TS,i).

X_{RC, i} = {(\frac{\partial \ln r}{\partial (\frac{- G_{TS, i}}{k_{B} T})})}_{G_{TS, j \neq i,} G_{n}}

(1.65)

Equivalently, it is possible to change the free energy G_n of an intermediate n and calculate its degree of rate control. This is schematically shown for the free energy of C* in Figure 1.15 and define the thermodynamic degree of rate control, X_TRC,_n, for an intermediate n can be defined as

X_{TRC, n} = {(\frac{\partial \ln r}{\partial (\frac{- G_{n}}{k_{B} T})})}_{G_{TS, i,} G_{m \neq n}}

(1.66)

and are, in principle identical and are collectively referred to as the generalized degree of rate control. In typical scenarios, the degree of rate control for transition states is positive, indicating that lowering the transition state energy increases the reaction rate, and the degree of rate control for surface intermediates is negative, which means that stabilizing the intermediate will poison the surface and result in a slower reaction rate. This can be interpreted more generally for heterogeneously catalyzed reactions. On the optimal catalyst the overall reaction proceeds along a “smooth” PES from the reactants to the products. Large potential energy barriers or deep wells result in a corrugated PES and slow reaction rates.

Figure 1.15

Illustration of the generalized degree of rate control XRC and degree of catalyst control XCC (dashed PES) for CH4 steam reforming. XRC(CH4-TS) = 0.8, XRC(C) = −0.26, XCC(C) = 0.11. Reproduced from ref. 65 with permission from AAAS, Copyright 2009.

View large Download slide

Illustration of the generalized degree of rate control X_RC and degree of catalyst control X_CC (dashed PES) for CH₄ steam reforming. X_RC(CH4-TS) = 0.8, X_RC(C) = −0.26, X_CC(C) = 0.11. Reproduced from ref. 65 with permission from AAAS, Copyright 2009.

1.6.3.2 Degree of Catalyst Control

Campbell’s degree of rate control is an extraordinarily useful concept for the analysis of reaction mechanisms, the identification of rate-limiting steps and the poisoning effect of adsorbates, but it implies that the energies of intermediates and transition states can be changed independently, one at a time. While this is the common procedure for a sensitivity analysis, we already know that the stabilities of surface intermediates and transition states are intimately linked through scaling and BEP relations, which find their physical foundation in the d-band model. It is therefore unrealistic to assume that catalysts can be found that, for example, reduce the activation energy barrier of a rate-limiting step without simultaneously affecting other parts of the PES. This leads to the idea of the degree of catalyst control X_CC, which is a constrained sensitivity analysis of a microkinetic model with respect to a reactivity descriptor E_i.⁶⁵ The definition of X_CC,i is analogous to the generalized degree of rate control and may be expressed as

X_{CC, i} = {(\frac{\partial \ln r}{\partial (\frac{- E_{i}}{k_{B} T})})}_{E_{j \neq i, BEP, Scaling}}

(1.67)

In it is indicated that all other descriptor values E_j≠i and the implemented scaling and BEP relations are kept constant when taking the partial derivative. The rationale for this definition is that the catalytic descriptors are the only independent variables and X_CC will indicate the search direction for an improved catalyst. The optimal catalyst is characterized by X_CC,i = 0 for all descriptors E_i. X_CC is proportional to the slope of the volcano curve and following it in the uphill direction will get you to the top!

The difference between X_RC and X_CC is shown graphically in Figure 1.15 for the sensitivity towards the binding energy of carbon, ΔE_C. Carbon adsorbs quite strongly and is located in a well of the PES, which explains the negative value (−0.26) for X_RC(C). Decreasing the depth of the well would smooth out the PES and increase the reaction rate. However, decreasing the stability of C* will also affect the stability of CH*, CH₂*, CH₃* surface intermediates and the associated transition states between them through scaling and BEP relations (dashed line in Figure 1.15). In particular, the transition state energy of CH₄ activation will shift up in energy. The large value of X_RC(CH₄-TS) = 0.8 indicates that CH₄ activation is a rate-determining step and increasing the barrier will have a large negative effect on the overall reaction rate. The overall effect of changing the descriptor ΔE_C on the PES is given by X_cc(C) = 0.11, indicating that a catalyst with stronger ΔE_C will have higher activity. The explanation is that the increased poisoning effect of C* is more than compensated for by lowering the barrier for the CH₄ activation step.

The degrees of rate and catalyst control both have their merits and do not conflict with each other. Determining which to use is dependent on the problem at hand. X_RC should be used to analyze the reaction mechanism and understand the importance of individual steps and adsorbates, whereas X_CC should be used to guide catalyst design.

1.7 CO Oxidation Catalyst Screening

A preliminary Sabatier analysis of the CO oxidation reaction was performed in Section 1.5.2, and an analytical solution derived under the assumption that the adsorption of CO and O₂ are quasi-equilibrated in Section 1.6. Now a numerical solution to the complete microkinetic model as a function of the descriptors ΔE_CO and ΔE_O will be formulated. The reaction mechanism will be analyzed in terms of rate and catalyst control, and at the end of this section, the effect of high surface coverages on the volcano curve will also be briefly addressed.

1.7.1 Numerical Microkinetic Model

Microkinetic models are typically custom-made and can be written in any programming or scripting language that the reader is most familiar with. It is recommended to choose a simulation environment that already provides lots of the functionality needed to solve the model, i.e. a fast ODE solver for stiff problems, a root-finder, and n-dimensional curve fitting routines, if parameter optimization is also required. Here Python is used, which is an object-oriented scripting language with good performance, and it is freely available for all popular operating systems (Windows, MacOS, Linux; http://python.org). The Python syntax should be self-explanatory as long as you are aware that comments are declared with ‘#’ and that list and array indices start with ‘0’, not with ‘1’. Readers who are interested in additional information are referred to the online tutorial available for Python newbies (http://docs.python.org/tutorial). Besides the standard Python interpreter, the NumPy (http://numpy.scipy.org), SciPy (http://www.scipy.org), mpmath (http://code.google.com/p/mpmath), and matplotlib (http://matplotlib.sourceforge.net) modules are also needed as extensions. These modules provide a variety of math routines and plotting functions, very similar to the Matlab environment. All examples here have been tested with Python 2.6, but other versions of Python should work as well. A complete Python script reproducing all results of this section can be found in the Appendix.

1.7.1.1 Constants and Conditions

It is a good idea to define the reaction conditions and useful physical constants globally at the beginning of the code, so that they are available to all functions. This is also the place to import additional modules, such as numpy.

import numpy as np

# Reaction conditions

T	=	600	#	K
PCO	=	0.33	#	bar
PO2	=	0.67	#	bar
PCO2	=	1	#	bar

# Physical constants and conversion factors

J2eV	=	6.24150974E18	#	eV/J
Na	=	6.0221415E23	#	mol−1
h	=	6.626068E-34 * J2eV	#	in eV*s
kb	=	1.3806503E-23 * J2eV	#	in eV/K
kbT	=	kb * T	#	in eV

1.7.1.2 From Descriptors to Rate Constants

The next important step is to implement the scaling and BEP relations together with any assumptions regarding entropies of transition states and surface states in order to obtain the necessary rate constants. Since this procedure will be needed multiple times for each combination of descriptors, it is defined as a function.

def get_rate_constants(EO,ECO):

# This function applies scaling and BEP to determine all rate constants

# depending on the two descriptors, EO and ECO, provided as input.

# Gas phase energies referenced to CO and O2

ECOg	=	EO2g = 0
ECO2g	=	−3.627	#	in eV calculated from DFT

# Scaling for EO2

EO2

0.89 * EO + 0.17

eq. (1.42)

# Gas phase entropies converted to eV/K

SCOg	=	197.66 * J2eV/Na	#	eV/K
SO2g	=	205.0 * J2eV/Na
SCO2g	=	213.74 * J2eV/Na

# Surface entropies are assumed to be zero

SCO = SO2 = SO = 0

# Reaction energies

dE	=	np.zeros(4)	#	array initialization
dE [0]	=	ECO − ECOg
dE [1]	=	EO2 − EO2g
dE [2]	=	2*EO − EO2
dE [3]	=	ECO2g − ECO − EO

# Entropy changes

dS	=	np.zeros(4)	#	array initialization
dS [0]	=	SCO − SCOg
dS [1]	=	SO2 − SO2g
dS [2]	=	2*SO − SO2
dS [3]	=	SCO2g − SCO − SO

# Activation energy barriers from BEP

Ea	=	np.zeros(4)	#	array initialization
Ea [2]	=	0.5 * EO + 1.39	#	eq. (1.43)
Ea [3]	=	−0.3 * (EO + ECO) + 0.02	#	eq. (1.44)

# Entropy changes to the transition state

STS	=	np.zeros(4)	#	array initialization
STS [0]	=	−0.25*SCOg	#	loss of CO entropy assumed
STS [1]	=	−0.25*SO2g	#	loss of O2 entropy assumed
STS [2]	=	0	#	surface reaction
STS [3]	=	0	#	surface reaction

# Calculate equilibrium and rate constants

K	=	[0]*4	#	equilibrium constants
kf	=	[0]*4	#	forward rate constants
kr	=	[0]*4	#	reverse rate constants

for i in range(4):

dG	=	dE [i] − T*dS [i]
K [i]	=	np.exp(−dG/kbT)

# Enforce Ea > 0, and Ea > dE, independent of what the descriptors are

Ea [i]	=	max([0,dE [i],Ea [i]])
kf [i]	=	kbT/h * np.exp(STS [i]/kb) * np.exp(−Ea [i]/kbT)
kr [i]	=	kf [i]/K [i]	#	enforce thermodynamic consistency

return (kf,kr)

In this function the scaling/BEP relations, eqn (1.42)–(1.44) have been used, and it has been arbitrarily assumed that CO and O₂ lose 25% of their gas-phase entropy when entering the transition state for adsorption. Using collision theory for the adsorption rate constants⁶⁶ of these two steps is a more elegant alternative and the reader is encouraged to see how it affects the result of the model. There are two more noteworthy points. First, the model ensures that the activation barriers, E_a,i, are all positive, even for descriptor pairs that would otherwise predict negative values for E_a,i. Second, thermodynamic consistency is enforced throughout the model by calculating the reverse rate constants from the equilibrium constants for each step.

1.7.1.3 Rate Equations

This step is simple. Rate equations for the individual steps are needed. The rates depend on the coverages, θ_i, and the rate constants. All parameters are passed as an array to the function that returns the rates.

def get_rates(theta,kf,kr):

# returns the rates depending on the current coverages theta

# Extract elements of theta and assign them

# to more meaningful variables

tCO	=	theta [0]	#	theta of CO
tO2	=	theta [1]	#	theta of O2
tO	=	theta [2]	#	theta of O
tstar	=	1.0 − tCO − tO2 − tO	#	site balance for tstar

# Calculate the rates: eqn (1.45)–(1.48)

rate	=	[0]*4	#	array with 4 zeros
rate [0]	=	kf [0] * PCO * tstar − kr [0] * tCO
rate [1]	=	kf [1] * PO2 * tstar − kr [1] * tO2
rate [2]	=	kf [2] * tO2 * tstar − kr [2] * tO * tO
rate [3]	=	kf[3] * tCO * tO − kr [3] * PCO2 * tstar * tstar

return rate

1.7.1.4 System of ODEs vs. Steady-state Equations

As discussed in Section 1.6.1 the microkinetic model may be solved as a system of ODEs or non-linear algebraic equations using the steady-state assumption. It turns out that, regardless of which approach is used, the function that must be passed to an ODE solver or numerical root-finding method is the same! Here, the more general case of the ODE system is chosen. Note that the previously defined function was named get_rates().

def get_odes(theta,t,kf,kr):

# returns the system of ODEs d(theta)/dt, calculated at the current value of theta (and time t)

rate

get_rates(theta,kf,kr)

calculate the current rates

# Time derivatives of theta

dt	=	[0] * 3
dt [0]	=	rate [0] − rate [3]	#	d(tCO)/dt
dt [1]	=	rate [1] − rate [2]	#	d(tO2)/dt
dt [2]	=	2 * rate [2] − rate [3]	#	d(tO).dt

return dt

1.7.1.5 Solving the Model

All that is left at this point is to set the values of the descriptors and solve the model. For single point calculations using only one descriptor set, it is preferable to use an ODE solver and use a completely empty surface as the initial guess for the coverages. The integration time required to reach steady-state depends on the problem and can vary from 1 × 10³ to 1 × 10⁸ seconds, or sometimes longer. Some might wonder why a reaction could possibly take 1 × 10⁸ seconds, or greater than 3 years, to equilibrate, but keep in mind that the goal is to scan a wide descriptor range, which will result in some very small rate constants. For this example, the Pt values for the descriptors and the odeint ODE solver provided by scipy were used. Note that the default tolerance and step size had to be tweaked in order for odeint to work for this stiff problem. For more complex problems the reader is advised to identify a better performing solver optimized for stiff systems of ODEs.

# Solve the model for Pt with EO = −1.25, ECO = −1.22

ECO	=	−1.22	#	CO oxidation descriptors
EO	=	−1.25
(kf,kr)	=	get_rate_constants (EO,ECO)	#	get the rate constants for the given descriptors

# Use scipy’s odeint to solve the system of ODEs from scipy.integrate import odeint

# As initial guess we assume an empty surface theta0 = (0., 0., 0.)

# Integrate the ODEs for 1E6 sec (enough to reach steady-state)

theta	=	odeint(get_odes,	#	system of ODEs
		theta0,	#	initial guess
		[0,1E6],	#	time span
		args = (kf,kr),	#	additional arguments to get_odes()
		h0 = 1E-36,	#	initial time step
		mxstep = 90000,	#	maximum number of steps
		rtol = 1E-12,	#	relative tolerance
		atol = 1E-15)	#	absolute tolerance

print_output(EO,ECO,theta [−1, :])

In the last line the helper function print_output (given below, but in the actual Python script it needs to be defined before it is called) is called with the descriptor values and the final coverages. The variable theta is a (t × c) matrix with the coverages c for each time step t. Usually, only the last row is of interest, conveniently accessed by using ‘−1’ as an index, but the time evolution is available as well for transient studies or verification that steady-state has been reached.

def print_output(EO,ECO,theta):

# Prints the solution of the model

(kf,kr)	=	get_rate_constants(EO,ECO)
rates	=	get_rates(theta,kf,kr)

print "For the descriptors EO =", EO, "and ECO =", ECO,"the result is:"

for r,rate in enumerate(rates):

print "Step", r,": rate =", rate,", kf = ", kf [r], ", kr = ", kr [r]

print "The coverages for CO*, O2*, and O* are:"

for t in theta:

print t

As a result of this exercise, the steady-state coverages on Pt of θ_CO = 0.215 ML, θ_O₂ = 8.4 × 10⁻⁴, and θ_O = 5.6 × 10⁻³ and a steady-state CO₂ formation rate of 6.16 × 10³ s⁻¹ are found.

1.7.2 Degree of Rate and Catalyst Control

Now that it is possible to predict a reaction rate for a single catalyst, only a few lines of extra code are needed to calculate the degree of rate control and the degree of catalyst control. However, the microkinetic model needs to be solved repeatedly for different parameters and it is more convenient if first a function, solve_ode, is written that takes the rate constants as input and returns the solution of the model. This function may look like this.

def solve_ode(kf, kr, theta0=(0.,0.,0.)):

# Solve the system of ODEs using scipy.integrate.odeint

# Assumes an empty surface as initial guess if nothing else is provided

from scipy.integrate import odeint

# Integrate the ODEs for 1E6 sec (enough to reach steady-state)

theta	=	odeint(get_odes,	#	system of ODEs
		theta0,	#	initial guess
		[0,1E6],	#	time span
		args = (kf,kr),	#	arguments to get_odes()
		h0 = 1E-36,	#	initial time step
		mxstep = 90000,	#	maximum number of steps
		rtol = 1E-12,	#	relative tolerance
		atol = 1E-15)	#	absolute tolerance

return theta [−1,:]

The degree of rate control, X_RC, is then calculated by repeatedly calling solve_ode with rate constants that are systematically changed for each step.

As input to calculate_Xrc the unaltered rate constants and the corresponding rate are provided.

def calculate_Xrc(r0,kf0,kr0):

# Calculates Xrc by systematically changing the

# rate constants of each step by 10% around the reference value

delta = 0.1	#	change of 10%
Xrc_rates = np.zeros(4)	#	array for storing rates
for s in range(4):	#	loop over all steps

# initialize rate constants with reference values

kf	=	kf0 [:]
kr	=	kr0 [:]
kf [s]	=	(1 + delta) * kf0 [s]
kr [s]	=	(1 + delta) * kr0 [s]

theta =		solve_ode(kf,kr)	#	Solve with the modified k’s
rates =		get_rates(theta, kf, kr)	#	Get the new rates
Xrc_rates [s]	=	rates [3]	#	Save the rate of CO2 production

# And calculate Xrc for all steps

Xrc = (Xrc_rates-r0)/(delta*r0)

return Xrc

In the example for CO oxidation on Pt with ΔE_O = −1.25 eV and ΔE_CO = −1.22 eV the degrees of rate control for steps (1) and (2) are zero (these steps are quasi-equilibrated), step (3) is almost exclusively rate-controlling with X_RC,3 = 0.99, and step (4) has an X_RC,4 = 0.01. In this example $\sum_{i} X_{RC, i} = 1$ ⁠, which is expected for a serial reaction mechanism.

Next, the degree of catalyst control X_cc is calculated by varying the values of ΔE_O and ΔE_CO by 0.05 eV. This step size should generally work well except in close proximity to the volcano peak, where smaller step sizes are preferable. Rather than providing the rate constants directly, the original values of ΔE_O and ΔE_CO are passed to the function calculate_Xcc and the corresponding rate constants are calculated on the fly.

def calculate_Xcc(r0,EO0,ECO0):

# Calculate Xcc by varying EO and ECO around their

# reference values by 0.05 eV

delta	=	0.05	#	change of 0.05 eV
Xcc_rates	=	np.zeros(2)	#	array for storing rates

# Start by modifying EO0

(kf,kr) = get rate constants(EO0 + delta,ECO0)

theta = solve ode(kf,kr)	#	Solve with the modified k’s
rates = get_rates(theta,kf,kr)	#	Get the new rates
Xcc_rates [0] = rates [3]	#	Save the rate of CO2 production

# Repeat for ECO0

(kf,kr) = get rate constants(EO0,ECO0 + delta)

theta = solve ode(kf,kr)	#	Solve with the modified k’s
rates = get_rates(theta,kf,kr)	#	Get the new rates
Xcc_rates [1] = rates [3]	#	Save the rate of CO2 production

Xcc = (np.log(Xcc_rates)–np.log(r0))/(–delta/kbT)

return Xcc

If this analysis is performed for CO oxidation on Pt, it yields X_cc(ΔE_O) = 1.38 and X_cc(ΔE_CO) = −0.28, which indicates that by stabilizing O* and destabilizing CO* the reaction rate can be increased. Indeed, on the two-dimensional volcano plot in Figure 1.16 it can clearly be seen that the maximum is located to the left of Pt in the direction of stronger O* binding. The resolution of the volcano plot is however not sufficient to judge whether stronger or weaker CO* binding would increase the activity. It appears as though ΔE_CO for Pt is already near to its optimal value.

Figure 1.16

CO oxidation activity, log10(TOF), as a function of the two catalytic descriptors ΔEO and ΔECO. (a) Numerical solution. (b) Analytical solution derived in Section 1.7, using eqn (1.54)–(1.58). The positions of the closed-packed surfaces of Pt, Pd, Cu, and Rh are indicated.

View large Download slide

CO oxidation activity, log₁₀(TOF), as a function of the two catalytic descriptors ΔE_O and ΔE_CO. (a) Numerical solution. (b) Analytical solution derived in Section 1.7, using eqn (1.54)–(1.58). The positions of the closed-packed surfaces of Pt, Pd, Cu, and Rh are indicated.

1.7.3 Two-dimensional CO Oxidation Volcano

After this thorough analysis of the microkinetic model results around a single point (Pt), the only remaining task is to calculate the CO oxidation rate over a larger descriptor range and visualize the results as a two-dimensional volcano plot as shown in Figure 1.16. As a comparison to the numerical solution in Figure 1.16(a), the analytical solution from Section 1.6 [eqn (1.54)–(1.58)] is given in Figure 1.16(b). If it is recalled that the only assumption that was made in the derivation of the analytical solution was that the CO and O₂ adsorption steps are quasi-equilibrated and that the numerically calculated degree of rate control of both steps was zero, it is not surprising that the plots are in perfect agreement.

Figure 1.16 can ideally be generated by simply looping over the entire parameter range, but in practice that will often result in convergence errors due to numerical precision problems with the largely varying rate constants in certain regions of the parameter space. For faster results it is also desirable to use a root-finding method (e.g. mpmath.findroot), which requires good initial guesses for convergence. Here, the alternating looping technique visualized in Figure 1.14 was used, and the solution of the previous grid point served as an initial guess for the next one. However, even this approach is not failure proof and a sanity check (e.g. are all coverages between 0 and 1?) of the solutions in each grid point must be performed before moving on to the next grid point. If the solution is not physical, then to a slower, but more stable, solver must be used in order to generate a new solution and then it may be possible to continue again with the fast root-finding method. An example of how this may be implemented is provided below.

# Generate the data for 2D volcano plots

# and save it in a matrix

gridpoints = 20

num_data = np.zeros([gridpoints,gridpoints])

# Generate the initial guess for the first point

(kf,kr)	=	get_rate_constants(−2.5,−2.5)
theta0	=	solve_ode(kf,kr)

# Start the scan through the parameter space

# Alternate directions for alternating columns

EO_range0	=	np.linspace(−2.5,0.0,gridpoints)
ECO_range0	=	np.linspace(−2.5,0.0,gridpoints)

for i, EO in enumerate(EO_range0):

if i%2:	#	true if odd
ECO_range = ECO_range0.copy()[::−1]	#	[::−1]

reverses the list

j_range = range(gridpoints)[::−1]

else:

ECO_range = ECO_range0.copy()

j_range = range(gridpoints)

for j,ECO in enumerate(ECO_range):

(kf, kr) = get_rate_constants(EO,ECO)

theta = solve_findroot(kf,kr,theta0)

# Check if solution is physical

if (theta > 1.0).any() or (theta < 0.0).any():

# Generate a new initial guess using odeint and solve again

theta0 = solve_ode(kf,kr)

theta = solve_findroot(kf,kr,theta0)

rates = get_rates(theta,kf,kr)

# Store the solutions

num_data [j_range [j],i] = rates [3]

# Store the numerical solution as initial guess

# for the next grid point

theta0 = theta.copy()

The numerical rate data collected in num_data can then be plotted using the contourf function from matplotlib. The full example code for solving the microkinetic model including the generation of the two-dimensional volcano graph can be found in the appendix.

1.7.4 Effect of Lateral Interactions

Throughout this chapter it has been assumed that binding energies and activation energy barriers are constant, but at higher surface coverages these quantities can be strongly influenced by adsorbate–adsorbate interactions. An accurate description of these interactions, specifically when they induce local ordering and lead to island formation or other characteristic patterns, would require the use of kinetic Monte Carlo (kMC) simulations. Since computational catalyst screening is largely based on the use of rather approximate scaling and BEP relations with mean absolute errors of approximately 0.2 eV, the use of highly sophisticated kMC methods does not seem appropriate. Adsorbate–adsorbate interactions can also be described within the mean-field approximation, which is less detailed but suitable for inclusion in microkinetic models.^48,67 In the simplest case of interacting atomic adsorbates the lateral interactions have been understood in terms of d-band structure changes, which, in turn, allows for the prediction of these interactions.⁶⁸ Just as lateral interactions alter the binding energy of adsorbates, they also affect the stability of transition states and thereby change the activation energy barriers for surface reactions. William Schneider’s group has studied this phenomenon in great detail by calculating coverage-dependent activation barriers for NO oxidation, which were later coupled with first principles-based cluster expansions for improved accuracy.^69,70 Current efforts to systematically improve the description of lateral interactions in microkinetic models are slowly closing the gap between the mean-field approximation and spatially resolved kMC simulations.⁶⁷

It is well known that adsorbate–adsorbate interactions can significantly change the surface energetics and overall catalytic activity, but the pressing question in computational catalyst screening is: Do lateral interactions alter the predicted reactivity trends? This question has been answered for the CO oxidation reaction and the short answer is ‘No’. The long answer can be found in ref. 49 and is briefly summarized here. Taking repulsive O–O, CO–CO and CO–O interactions directly into account for adsorbates and indirectly for activation barriers through BEP relations, a coverage-dependent microkinetic model for CO oxidation can be constructed and used for screening purposes. Because the lateral interactions vary across catalysts and no correlation between the interactions and the activity descriptors ΔE_O and ΔE_CO could be identified, the volcano in Figure 1.17(a) was calculated using the assumption that the interaction energies for Pt(111) are constant for the whole range of descriptors. The result is that in the low coverage regime there is no difference between coverage corrected and uncorrected volcanoes, and at high coverage the effect of surface poisoning is reduced by the destabilization of surface species. The reduced poisoning effect leads to a higher activity and a broadening of the volcano peak, as clearly seen in Figure 1.17(a), but the position of the peak and the relative order of metals remain constant [not shown here, but a contour plot with finer resolution around the peak is given as Figure 3(c) of ref. 49]. Using catalyst-independent lateral interaction energies is instructive for showing the effect of lateral interactions on the shape of the volcano curve, but the preferred visualization technique is to indicate the shift of the metal positions from their low coverage descriptor values to their high coverage values on the uncorrected volcano plot. The result of this technique is shown in Figure 1.17(b). The corrected metal positions all fall onto the line that indicates the transition from the high to the low coverage regime. As long as surface oxide formation can be excluded, it is possible to think of these strongly adsorbing catalysts as operating at their saturation coverage under reaction conditions. The effective (or coverage adjusted) descriptor values can then be evaluated at the corresponding surface coverage and used to indicate the coverage-corrected position. While the long answer is fundamentally interesting, the important take-home message is that lateral interactions usually do not alter the position of the volcano peak nor change the relative order of the activity of metal catalysts.

Figure 1.17

View large Download slide

The effect of adsorbate–adsorbate interactions on the CO oxidation volcano. (a) Adsorbate interactions characteristic for Pt were included in the generation of the volcano plot. (b) The “standard” activity volcano is shown, but the positions of the metals have been corrected to account for the effect of surface coverage. Reproduced from ref. 49 with permission from Springer Nature, Copyright 2010.

1.8 Conclusions

This chapter has provided a tutorial style overview of the rapidly emerging field of computational catalyst screening and the reader should feel adequately informed and prepared to attempt a catalyst screening project on their own. Although the field is still developing, it has already reached such a high level of complexity that several topics could only be mentioned in passing and references for further study have had to be provided. The descriptor-based approach relies heavily on assumptions and simplifications, but its success is well founded in the numerous examples of new catalyst discoveries in recent years. It is this success that emphasizes the applicability, strength, and robustness of the method for the search for novel and improved heterogeneous catalysts.

This chapter is concluded by stating that computational catalyst screening works in some sense just like a professional dating or matchmaking service. How so? Well, a client (= reaction) who is seeking the help of a professional matchmaker (= us) in order to find their perfect partner (= catalyst) will first have to go through an in-depth interview. During the interview, the matchmaker thoroughly analyzes the client’s character and personal situation and describes them with certain adjectives (= descriptors), such as attractive (= affinity to adsorbate ‘X’), financially stable (= low cost material), emotionally stable (=poisoning resistant), desire to get married (= long-term stability), and so on. Once the matchmaker has developed a full profile of the client, they will enter that information into a database and run a search algorithm that systematically screens a database for potential partners. The screening process will yield a list of, perhaps 10, matches that all have a certain compatibility index (= deviation from the optimal descriptor value), and on the following dates (= experimental testing) the client will need to find out whether the love of their life (= optimal catalyst) was amongst the identified candidates. The matchmaker’s key to success is the creation of an accurate client profile (= identification of the right descriptors), the quality of the matching algorithm (= scaling, BEP, and kinetic model) and the size of the client (= catalyst material) database.

As you can see, there are quite a few parallels between professional matchmaking services and computational catalyst screening. The next time you find yourself in a situation where you need to explain to a non-expert what you do for a living, try starting with: “I am a professional matchmaker and my clients are slow chemical reactions…”!

Appendix

Python Code for the numerical solution of the microkinetic model for CO oxidation and the generation of a two-dimensional volcano plot.

# !/usr/bin/env python

import numpy as np

# Reaction conditions

T	=	600	#	K
PCO	=	0.33	#	bar
PO2	=	0.67	#	bar
PCO2	=	1	#	bar

# Physical constants and conversion factors

J2eV	=	6.24150974E18	#	eV/J
Na	=	6.0221415E23	#	mol−1
h	=	6.626068E-34 * J2eV	#	in eV*s
kb	=	1.3806503E-23 * J2eV	#	in eV/K
kbT	=	kb*T	#	in eV

def get_rate_constants(EO,ECO):

# This function applies scaling and BEP to determine all rate constants

# depending the two descriptors EO, ECO provided as input.

# Gas phase energies with reference to CO and O2

ECOg	=	EO2g = 0
ECO2g	=	−3.627	#	in eV calculated from DFT

# Scaling for EO2

EO2 = 0.89 * EO + 0.17 # eq. (1.42)

# Gas phase entropies converted to eV/K

SCOg	=	197.66 * J2eV/Na	#	eV/K
SO2g	=	205.0 * J2eV/Na
SCO2g	=	213.74 * J2eV/Na

# Surface entropies are assumed to be zero

SCO = SO2 = SO = 0

# Reaction energies

dE	=	np.zeros(4)	#	array initialization
dE [0]	=	ECO − ECOg
dE [1]	=	EO2 − EO2g
dE [2]	=	2*EO − EO2
dE [3]	=	ECO2g − ECO − EO

# Entropy changes

dS	=	np.zeros(4)	#	array initialization
dS [0]	=	SCO − SCOg
dS [1]	=	SO2 − SO2g
dS [2]	=	2*SO − SO2
dS [3]	=	SCO2g − SCO − SO

# Activation energy barriers from BEP

Ea	=	np.zeros(4)	#	array initialization
Ea [2]	=	0.5*EO + 1.39	#	eq. (1.43)
Ea [3]	=	−0.3*(EO + ECO) + 0.02	#	eq. (1.44)

# Entropy changes to the transition state

STS	=	np.zeros(4)	#	array initialization
STS [0]	=	−0.25*SCOg	#	loss of CO entropy assumed
STS [1]	=	−0.25*SO2g	#	loss of O2 entropy assumed
STS [2]	=	0	#	surface reaction
STS [3]	=	0	#	surface reaction

# Calculate equilibrium and rate constants

K	=	[0]*4	#	equilibrium constants
kf	=	[0]*4	#	forward rate constants
kr	=	[0] *4	#	reverse rate constants

for i in range(4):

dG = dE [i] − T*dS [i]

K [i] = np.exp(–dG/kbT)

# enforce Ea > 0, and Ea > dE

# independent of what the descriptors are

Ea [i]	=	max([0,dE [i], Ea [i]])
kf [i]	=	kbT/h * np.exp(STS [i]/kb) * np.exp (−Ea [i]/kbT)
kr [i]	=	kf [i]/K [i] # enforce thermodynamic consistency

return (kf,kr)

def get_rates(theta,kf,kr):

# returns the rates depending on the current coverages theta

# Extract elements of theta and assign them

# to more meaningful variables

tCO	=	theta [0]	#	theta of CO
tO2	=	theta [1]	#	theta of O2
tO	=	theta [2]	#	theta of O
tstar	=	1.0 − tCO − tO2 − tO	#	site balance for tstar

# Calculate the rates: eqns (1.45)–(1.48)

rate	=	[0]*4	# array with 4 zeros
rate [0]	=	kf [0] * PCO * tstar − kr [0] * tCO
rate [1]	=	kf [1] * PO2 * tstar − kr [1] * tO2
rate [2]	=	kf [2] * tO2 * tstar − kr [2] * tO * tO
rate [3]	=	kf [3] * tCO * tO − kr [3] * PCO2 * tstar * tstar

return rate

def get_odes(theta,t,kf,kr):

# returns the system of ODEs d(theta)/dt

# calculated at the current value of theta (and time t, not used)

rate = get_rates(theta,kf,kr) # calculate the current rates

# Time derivatives of theta

dt	=	[0]*3
dt [0]	=	rate [0] − rate [3]	#	d(tCO)/dt
dt [1]	=	rate [1] − rate [2]	#	d(tO2)/dt
dt [2]	=	2 * rate [2] − rate [3]	#	d(tO)/dt

return dt

def print_output(EO,ECO,theta):

# Prints the solution of the model

(kf,kr)	=	get_rate_constants(EO,ECO)
rates	=	get_rates(theta,kf,kr)

print "For the descriptors EO =",EO,"and ECO =",

ECO,"the result is:"

for r,rate in enumerate(rates):

print "Step",r,": rate =",rate,", kf =",kf [r],", kr=",kr [r]

print "The coverages for CO*, O2*, and O* are:"

for t in theta:

print t

def solve_ode(kf,kr,theta0=(0.,0.,0.)):

# Solve the system of ODEs using scipy.integrate.odeint

# Assumes an empty surface as initial guess if nothing else is provided

from scipy.integrate import odeint

# Integrate the ODEs for 1E6 sec (enough to reach steady-state)

theta	=	odeint(get odes,	#	system of ODEs
		theta0,	#	initial guess
		[0,1E6],	#	time span
		args = (kf,kr),	#	arguments to get odes()
		h0 = 1E-36,	#	initial time step
		mxstep = 90000,	#	maximum number of steps
		rtol = 1E-12,	#	relative tolerance
		atol = 1E-15)	#	absolute tolerance

return theta [−1,:]

def solve_findroot(kf,kr,theta0):

# Use mpmath’s findroot to solve the model

from mpmath import mp, findroot

mp.dps = 25

mp.pretty = True

def get_findroot_eqns(*args):

return get_odes(args,0,kf,kr)

theta = findroot(get_findroot_eqns, tuple(theta0), multidimensional=True)

return np.array(theta)

def solve_analytically(kf,kr):

# Using the analytical solution from section 7, eq. (48)–(52)

# The equation for W was rewritten to avoid numerical precision errors

x3	=	kf [2] * kf [1]/kr [1] * PO2
y3	=	kr [2]
x4	=	kf [3] * kf [0]/kr [0] * PCO
y4	=	kr [3] * PCO2

numerator = 8*y3*(2*x3 + y4)

denominator = x4**2

if (numerator < 1.0e-4*denominator):

W = (x4/(4*y3))*(0.5*(numerator/denominator))

else:

W = (x4/(4*y3))*(−1 + np.sqrt(1 + numerator/denominator))

tstar = 1./(1 + W + kf [0]/kr [0] * PCO + kf [1]/kr [1] * PO2)

tco	=	kf [0]/kr [0] * PCO * tstar
tO2	=	kf [1]/kr [1] * PO2 * tstar
tO	=	W * tstar

return (tCO,tO2,tO)

def calculate_Xrc(r0,kf0,kr0):

# Calculates Xrc by systematically changing the

# rate constants of each step by 10% around the reference value

delta = 0.1	#	changeof10%
Xrc_rates = np.zeros(4)	#	array for storing rates
for s in range(4):	#	loop over all steps

# initialize rate constants with reference values

kf	=	kf0 [:]
kr	=	kr0 [:]
kf [s]	=	(1 + delta) * kf0 [s]
kr [s]	=	(1 + delta) * kr0 [s]

# Solve the ODEs with the modified k’s

theta = solve_ode(kf,kr)
rates = get_rates (theta,kf,kr)	#	Get the new rates
Xrc_rates [s] = rates [3]	#	Save the rate of CO2 production

# And calculate Xrc for all steps

Xrc = (Xrc_rates-r0)/(delta*r0)

return Xrc

def calculate_Xcc(r0,EO0,ECO0):

# Calculate Xcc by varying EO and ECO around their

# reference values by 0.05 eV

delta = 0.05	change of 0.05 eV
Xcc_rates = np.zeros(2)	array for modified rates

# Start by modifying EO0

(kf,kr) = get_rate_constants(EO0 + delta,ECO0)
theta = solve_ode(kf,kr)	#	Solve with the modified k’s
rates = get_rates(theta,kf,kr)	#	Get the new rates
Xcc_rates [0] = rates [3]	#	Save the rate of CO2 production

# Repeat for ECO0

(kf,kr) = get_rate_constants(EO0,ECO0 + delta)
theta = solve ode(kf,kr)	#	Solve with the modified k’s
rates = get_rates(theta,kf,kr)	#	Get the new rates
Xcc_rates [1] = rates [3]	#	Save the rate of CO2 production
Xcc = (np.log(Xcc_rates)-np.log(r0))/(−delta/kbT)

return Xcc

# Solve the model for Pt with EO = ECO = −1.22

ECO	=	−1.22	#	CO oxidation descriptors
EO	=	−1.25

(kf0,kr0) = get_rate_constants (EO,ECO)

get the rate constants for the given descriptor

theta0 = solve_ode(kf0,kr0)	#	Solve the model for the reference values
print_output(EO,ECO,theta0)	#	Print the output

# Sensitivity Analysis

rates = get_rates(theta0,kf0,kr0)

r0 = rates [3]	#	CO2 production rate as reference
Xrc = calculate_Xrc(r0,kf0,kr0)	#	Call the Xrc function

print "\nThe Degrees of Rate control are:"

print Xrc

Xcc = calculate_Xcc(r0,EO,ECO) # Call the Xcc function print "\nThe Degrees of Catalyst control for EO and ECO are:"

print Xcc

print "\nStarting to generate 2D data…"

# Generate the data for 2D volcano plots

# and save it in a matrix

# numerical, analytical at the same time gridpoints = 20

num_data	=	np.zeros([gridpoints,gridpoints])
ana_data	=	np.zeros([gridpoints,gridpoints])

# Generate the initial guess for the first point

(kf,kr)	=	get_rate_constants(−2.5,−2.5)
theta0	=	solve_ode(kf,kr)

# Start the scan through the parameter space

# Alternate directions for alternating columns

EO_range0	=	np.linspace(−2.5,0.0,gridpoints)
ECO_range0	=	np.linspace(−2.5,0.0,gridpoints)

for i,EO in enumerate(EO_range0):

if i%2:	#	true if odd
ECO_range = ECO_range0.copy()[::−1]	#	reverses the list
j_range = range(gridpoints)[::−1]

else:

ECO_range = ECO_range0.copy()

j_range = range(gridpoints)

for j,ECO in enumerate(ECO_range):

(kf,kr) = get_rate_constants(EO,ECO)

theta = solve_findroot(kf,kr,theta0)

# Check if solution is physical

if (theta > 1.0).any() or (theta o 0.0).any():

# Generate a new initial guess using odeint and solve again

theta0 = solve_ode(kf,kr)

theta = solve_findroot(kf,kr,theta0)

rates = get_rates(theta,kf,kr)

# Store the solutions

num_data [j_range [j],i] = rates [3] # the numerical solution

theta_ana = solve_analytically(kf,kr)

ana_data [j_range [j],i] = get_rates(theta_ana, kf,kr)[3]

# Store the numerical solution as

initial guess for the next grid point

theta0 = theta.copy()

# Generate a 2D volcano plot using matplotlib

from matplotlib.pyplot import *

for d,data in enumerate([num_data, ana_data]):

figure()

levels = np.linspace(−30,6,19)

graph = contourf(EO_range0,ECO_range0,np.log10(data), levels,cmap = cm.jet)

xlabel(r’$\Delta$E$_O$/eV’)

ylabel(r’$\Delta$E$_{CO}$/eV’)

colorbar()

scatter([−1.67,−1.25,−1.30,−2.25],

[−0.34,−1.22,−1.37,−1.58],color = ‘k’)

dx = 0.04

dy = 0.04

text(−1.67 + dx,−0.34 + dy,’Cu’)

text(−1.25 + dx,−1.22 + dy,’Pt’)

text(−1.30 + dx,−1.37,’Pd’)

text(−2.25 + dx,−1.58 + dy,’Rh’)

axis([−2.5,0.,−2.5,0.])

savefig(’COox_Volcano_’+str(d)+’.png’)

Acknowledgements

The author would like to thank Dr Hakan Demir for his input and proofreading the chapter.

References

Yaccato

Carhart

Hagemeyer

Herrmann

Lesik

Strasser

Volpe

Turner

Weinberg

Grasselli

R. K.

Brooks

C. J.

Pigos

J. M.

Comb. Chem. High Throughput Screening

2010

, vol.

(pg.

318

330

)

Google Scholar

Greeley

Jaramillo

T. F.

Bonde

Chorkendorff

I. B.

Nørskov

J. K.

Nat. Mater.

2006

, vol.

(pg.

909

913

)

Google Scholar

PubMed

Mavrikakis

Nat. Mater.

2006

, vol.

(pg.

847

848

)

Google Scholar

PubMed

Somorjai

G. A.

Top. Catal.

2010

, vol.

(pg.

311

325

)

Google Scholar

Besenbacher

Chorkendorff

Clausen

B. S.

Hammer

Molenbroek

A. M.

Nørskov

J. K.

Stensgaard

Science

1998

, vol.

279

pg.

1913

Google Scholar

PubMed

Jacobsen

C. J.

Dahl

Clausen

B. S.

Bahn

Logadottir

Nørskov

J. K.

J. Am. Chem. Soc.

2001

, vol.

123

(pg.

8404

8405

)

Google Scholar

PubMed

Toulhoat

Raybaud

J. Catal.

2003

, vol.

216

(pg.

)

Google Scholar

Strasser

Fan

Devenney

Weinberg

W. H.

Liu

Nørskov

J. K.

J. Phys. Chem. B

2003

, vol.

107

(pg.

11013

11021

)

Google Scholar

Nilekar

A. U.

Sasaki

Farberow

C. A.

Adzic

R. R.

Mavrikakis

J. Am. Chem. Soc.

2011

, vol.

133

(pg.

18574

18576

)

Google Scholar

PubMed

Linic

Jankowiak

Barteau

M. A.

J. Catal.

2004

, vol.

224

(pg.

489

493

)

Google Scholar

Greeley

Mavrikakis

Nat. Mater.

2004

, vol.

(pg.

810

815

)

Google Scholar

PubMed

Jaramillo

T. F.

Jørgensen

K. P.

Bonde

Nielsen

J. H.

Horch

Chorkendorff

Science

2007

, vol.

317

(pg.

100

102

)

Google Scholar

PubMed

Alayoglu

Nilekar

A. U.

Mavrikakis

Eichhorn

Nat. Mater.

2008

, vol.

(pg.

333

338

)

Google Scholar

PubMed

Studt

Abild-Pedersen

Bligaard

Sorensen

R. Z.

Christensen

C. H.

Nørskov

J. K.

Science

2008

, vol.

320

(pg.

1320

1322

)

Google Scholar

PubMed

Greeley

Stephens

I. E. L.

Bondarenko

A. S.

Johansson

T. P.

Hansen

H. A.

Jaramillo

T. F.

Rossmeisl

Chorkendorff

Nørskov

J. K.

Nat. Chem.

2009

, vol.

(pg.

552

556

)

Google Scholar

PubMed

Zhang

Vukmirovic

M. B.

Sasaki

Nilekar

A. U.

Mavrikakis

Adzic

R. R.

J. Am. Chem. Soc.

2005

, vol.

127

(pg.

12480

12481

)

Google Scholar

PubMed

Grabow

L. C.

Studt

Abild-Pedersen

Petzold

Kleis

Bligaard

Nørskov

J. K.

Angew. Chem., Int. Ed.

2011

, vol.

(pg.

4601

4605

)

Google Scholar

Grabow

L. C.

Studt

Abild-Pedersen

Petzold

Kleis

Bligaard

Nørskov

J. K.

Angew. Chem.

2011

, vol.

123

(pg.

4697

4701

)

Google Scholar

Nørskov

J. K.

Bligaard

Rossmeisl

Christensen

C. H.

Nat. Chem.

2009

, vol.

(pg.

)

Google Scholar

PubMed

Nørskov

J. K.

Abild-Pedersen

Studt

Bligaard

Proc. Natl. Acad. Sci. U. S. A.

2011

, vol.

108

(pg.

937

943

)

Google Scholar

PubMed

Bligaard

Nørskov

J. K.

Dahl

Matthiesen

Christensen

C. H.

Sehested

J. Catal.

2004

, vol.

224

(pg.

206

217

)

Google Scholar

Andersson

M. P.

Bligaard

Kustov

Larsen

K. E.

Greeley

Johannessen

Christensen

C. H.

Nørskov

J. K.

J. Catal.

2006

, vol.

239

(pg.

501

506

)

Google Scholar

Sehested

Larsen

Kustov

Frey

Johannessen

Bligaard

Andersson

Nørskov

Christensen

Top. Catal.

2007

, vol.

(pg.

)

Google Scholar

Andersson

Bligaard

Christensen

C. H.

Johannessen

Kustov

Larsen

K. E.

Nørskov

J. K.

and

Sehested

, Process and Catalyst For Hydrogenation of Carbon Oxides, US Pat., US7790776B2,

2010

Munter

T. R.

Landis

D. D.

Abild-Pedersen

Jones

Wang

Bligaard

Comput. Sci. Discovery

2009

, vol.

pg.

015006

Google Scholar

Kohn

Sham

L. J.

Phys. Rev.

1965

, vol.

140

pg.

A1133

Google Scholar

Hohenberg

Kohn

Phys. Rev.

1964

, vol.

136

pg.

B864

Google Scholar

Perdew

J. P.

Chevary

J. A.

Vosko

S. H.

Jackson

K. A.

Pederson

M. R.

Singh

D. J.

Fiolhais

Phys. Rev. B: Condens. Matter Mater. Phys.

1992

, vol.

pg.

6671

Google Scholar

Perdew

J. P.

Chevary

J. A.

Vosko

S. H.

Jackson

K. A.

Pederson

M. R.

Singh

D. J.

Fiolhais

Phys. Rev. B: Condens. Matter Mater. Phys.

1993

, vol.

pg.

4978

Google Scholar

Perdew

J. P.

Burke

Ernzerhof

Phys. Rev. Lett.

1996

, vol.

(pg.

3865

3868

)

Google Scholar

PubMed

Hammer

Hansen

Nørskov

Phys. Rev. B: Condens. Matter Mater. Phys.

1999

, vol.

(pg.

7413

7421

)

Google Scholar

Hummelshøj

J. S.

Abild-Pedersen

Studt

Bligaard

Nørskov

J. K.

Angew. Chem., Int. Ed.

2011

, vol.

(pg.

272

274

)

Google Scholar

Sholl

D. S.

and

Steckel

J. A.

Density Functional Theory: A Practical Introduction

Wiley-VCH, Verlag

2009

, vol.

135

Google Scholar

Hammer

Nørskov

J. K.

Adv. Catal.

2000

, vol.

(pg.

129

)

Google Scholar

Hammer

Nørskov

J. K.

Surf. Sci.

1995

, vol.

343

(pg.

211

220

)

Google Scholar

Xin

Linic

J. Chem. Phys.

2010

, vol.

132

(pg.

221101

221104

)

Google Scholar

Abild-Pedersen

Greeley

Studt

Rossmeisl

Munter

T. R.

Moses

P. G.

Skúlason

Bligaard

Nørskov

J. K.

Phys. Rev. Lett.

2007

, vol.

(pg.

016104

016105

)

Google Scholar

Fernandez

E. M.

Moses

P. G.

Toftelund

Hansen

H. A.

Martinez

J. I.

Abild-Pedersen

Kleis

Hinnemann

Rossmeisl

Bligaard

Nørskov

J. K.

Angew. Chem., Int. Ed.

2008

, vol.

(pg.

4683

4686

)

Google Scholar

Nørskov

J. K.

Bligaard

Hvolbaek

Abild-Pedersen

Chorkendorff

Christensen

C. H.

Chem. Soc. Rev.

2008

, vol.

(pg.

2163

2171

)

Google Scholar

Gómez-Díaz

López

J. Phys. Chem. C

2011

, vol.

115

(pg.

5667

5674

)

Google Scholar

Gómez-Díaz

Vargas-Fuentes

López

J. Chem. Phys.

2011

, vol.

135

pg.

124707

Google Scholar

Jones

Studt

Abild-Pedersen

Nørskov

J. K.

Bligaard

Chem. Eng. Sci.

2011

, vol.

(pg.

6318

6323

)

Google Scholar

Nørskov

J. K.

Bligaard

Logadottir

Bahn

Hansen

L. B.

Bollinger

Bengaard

Hammer

Sljivancanin

Mavrikakis

Dahl

Jacobsen

C. J. H.

J. Catal.

2002

, vol.

209

(pg.

275

278

)

Google Scholar

Wang

Temel

Shen

Jones

Grabow

L. C.

Studt

Bligaard

Abild-Pedersen

Christensen

C. H.

Nørskov

J. K.

Catal. Lett.

2011

, vol.

141

(pg.

370

373

)

Google Scholar

Wang

Petzold

Tripkovic

Kleis

Howalt

J. G.

Skúlason

Fernández

E. M.

Hvolbæk

Jones

Toftelund

Falsig

Björketun

Studt

Abild-Pedersen

Rossmeisl

Nørskov

J. K.

Bligaard

Phys. Chem. Chem. Phys.

2011

, vol.

(pg.

20760

20765

)

Google Scholar

Brønsted

J. N.

Chem. Rev.

1928

, vol.

(pg.

231

338

)

Google Scholar

Evans

M. G.

Polanyi

Trans. Faraday Soc.

1938

, vol.

(pg.

)

Google Scholar

Grabow

L. C.

Gokhale

A. A.

Evans

S. T.

Dumesic

J. A.

Mavrikakis

J. Phys. Chem. C

2008

, vol.

112

(pg.

4608

4617

)

Google Scholar

Grabow

L. C.

Hvolbæk

Nørskov

J. K.

Top. Catal.

2010

, vol.

(pg.

298

310

)

Google Scholar

Vojvodic

Calle-Vallejo

Guo

Wang

Toftelund

Studt

Martínez

J. I.

Shen

Man

I. C.

Rossmeisl

Bligaard

Noørskov

J. K.

Abild-Pedersen

J. Chem. Phys.

2011

, vol.

134

pg.

244509

Google Scholar

Sabatier

Ber. Dtsch. Chem. Ges.

1911

, vol.

(pg.

1984

2001

)

Google Scholar

Falsig

Bligaard

Rass-Hansen

Kustov

A. L.

Christensen

C. H.

Nørskov

J. K.

Top. Catal.

2007

, vol.

(pg.

117

120

)

Google Scholar

Falsig

Hvolbæk

Kristensen

I. S.

Jiang

Bligaard

Christensen

C. H.

Nørskov

J. K.

Angew. Chem., Int. Ed.

2008

, vol.

(pg.

4835

4839

)

Google Scholar

Falsig

, Understanding Catalytic Activity Trends for NO Decomposition and CO Oxidation using Density Functional Theory and Microkinetic Modeling, PhD Thesis, Technical University of Denmark, Kongens Lyngby, Denmark, 2010.

Ertl

Angew. Chem., Int. Ed.

2008

, vol.

(pg.

3524

3535

)

Google Scholar

Honkala

Hellman

Remediakis

I. N.

Logadottir

Carlsson

Dahl

Christensen

C. H.

Nørskov

J. K.

Science

2005

, vol.

307

(pg.

555

558

)

Google Scholar

Gokhale

A. A.

Kandoi

Greeley

J. P.

Mavrikakis

Dumesic

J. A.

Chem. Eng. Sci.

2004

, vol.

(pg.

4679

4691

)

Google Scholar

Stoltze

and

Nørskov

J. K.

, Theoretical Modelling of Catalytic Reactions, in

Handbook of Heterogeneous Catalysis

, ed.

Ertl

Knözinger

Schüth

and

Weitkamp

Wiley-VCH Verlag GmbH & Co. KGaA

Weinheim, Germany

, 2nd edn,

2008

Google Scholar

Dumesic

J. A.

Rudd

D. F.

Aparicio

L. L.

Rekoske

J. E.

and

Treviño

A. A.

The Microkinetics of Heterogeneous Catalysis

American Chemical Society

Washington, DC

1993

Google Scholar

Freund

H.-J.

Meijer

Scheffler

Schlogl

Wolf

Angew. Chem., Int. Ed.

2011

, vol.

(pg.

10064

10094

)

Google Scholar

Cramer

C. J.

Essentials of Computational Chemistry: Theories and Models

John Wiley & Sons Ltd

Chichester, England

, 2nd edn,

2004

Google Scholar

Campbell

C. T.

Top. Catal.

1994

, vol.

(pg.

353

366

)

Google Scholar

Campbell

C. T.

J. Catal.

2001

, vol.

204

(pg.

520

524

)

Google Scholar

Stegelmann

Andreasen

Campbell

C. T.

J. Am. Chem. Soc.

2009

, vol.

131

(pg.

8077

8082

)

Google Scholar

Nørskov

J. K.

Bligaard

Kleis

Science

2009

, vol.

324

(pg.

1655

1656

)

Google Scholar

Cortright

R. D.

Dumesic

J. A.

Adv. Catal.

2001

, vol.

pg.

161

Google Scholar

Hellman

Honkala

J. Chem. Phys.

2007

, vol.

127

pg.

194704

Google Scholar

İnoĝlu

Kitchin

Phys. Rev. B: Condens. Matter Mater. Phys.

2010

, vol.

(pg.

)

Google Scholar

Getman

R. B.

Schneider

W. F.

ChemCatChem

2010

, vol.

(pg.

1450

1460

)

Google Scholar

Schmidt

D. J.

Wolverton

Schneider

W. F.

J. Catal.

2012

, vol.

286

(pg.

)

Google Scholar

Nilsson

Pettersson

L. G. M.

Hammer

Bligaard

Christensen

C. H.

Nørskov

J. K.

Catal. Lett.

2005

, vol.

100

(pg.

111

114

)

Google Scholar

Log in

Chapter 1: Computational Catalyst Screening

1.1 Introduction

1.1.1 A Walk Through a Computational Catalyst Design Process: Methanation

1.2 Starting from the Electronic Structure

1.2.1 Density Functional Theory

1.2.2 The d-Band Model

1.3 Identifying the Right Descriptor Set

1.3.1 Scaling Relations for Surface Intermediates

1.3.1.1 Hydrogen-containing Molecules

1.3.1.2 Extension to More Complex Surface Intermediates

1.3.2 Scaling Relations for Transition States: The Brønsted–Evans–Polanyi Relationship

1.4 The Sabatier Principle and the Volcano Curve

1.4.1 Sabatier Analysis

1.5 Sabatier Analysis in Practice

1.5.1 First Example: Ammonia Synthesis

1.5.2 Second Example: CO Oxidation

1.6 Notes on Microkinetic Modeling

1.6.1 Numerical Solution Strategies

1.6.2 Entropy and Enthalpy Corrections

1.6.3 Microkinetic Model Analysis

1.6.3.1 Degree of Rate Control

1.6.3.2 Degree of Catalyst Control

1.7 CO Oxidation Catalyst Screening

1.7.1 Numerical Microkinetic Model

1.7.1.1 Constants and Conditions

1.7.1.2 From Descriptors to Rate Constants

1.7.1.3 Rate Equations

1.7.1.4 System of ODEs vs. Steady-state Equations

1.7.1.5 Solving the Model

1.7.2 Degree of Rate and Catalyst Control

1.7.3 Two-dimensional CO Oxidation Volcano

1.7.4 Effect of Lateral Interactions

1.8 Conclusions

Appendix

Acknowledgements

References

This Feature Is Available To Subscribers Only