Skip to Main Content
Skip Nav Destination

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 literature examples and provides step-by-step instructions with solved numerical problems for NH3 synthesis and CO oxidation on transition metal surfaces. The theoretical foundation of the screening approach, including the d-band model, linear scaling relations, Sabatier analysis, basic microkinetic modeling and the analysis of such models, is 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 literature and resources useful.

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 2128 possible permutations. If we simply assume that a typical central processing unit (CPU) can generate 109 bit flips per second (∼1 GHz), then the total time that is required to test all possible permutations is 2128/109=3.4×1029 seconds or 1022 years! For obvious reasons a brute-force attack is most likely 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, one can 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. have synthesized over 50 000 catalysts and tested them in more than 250 000 experiments for low, medium and high temperature water–gas shift conditions.1  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 large property database. Populating this property database with DFT data is the most time-consuming step in this process, but the resulting database is applicable to any reaction and only has to be generated once. 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 counterpart. 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.4  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,5  Jacobsen et al. found an active Co/Mo alloy for ammonia synthesis by interpolation in the periodic table,6  and Toulhoat and Raybaud showed that the metal–sulfur bond strength can correctly predict trends in hydro-desulfurization activity on metal–sulfide catalysts.7  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,10  near-surface alloys for hydrogen activation11  and evolution2,12 , Ru/Pt core–shell particles for preferential CO oxidation,13  Ni/Zn alloys for the selective hydrogenation of acetylene,14  Sc and Y modified Pt and Pd electrodes15  and mixed-metal Pt monolayer catalysts16  for electro-chemical oxygen reduction, and the rediscovery of Pt as the most active and selective catalyst for the production of hydrogen cyanide.17,18 

The most comprehensive example of a success story in computational catalyst design comes from the group of Jens Nørskov, who has pioneered the descriptor-based design approach and has applied it to numerous reactions.19,20  In several publications his group has studied the methanation reaction (CO+2H2 → CH4+H2O), 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 one must first answer the question: “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 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, one notices that FeNi3 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, FeNi3 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 FeNi3 alloy is significantly more active than its components. 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.24 

Figure 1.1

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. Ediss(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 H2 as a function of Ni content. [(a), (b), (d) are reprinted by permission from Macmillan Publishers Ltd: Nature Chemistry, Ref. 15, 2009, (c) is from ref. 25]

Figure 1.1

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. Ediss(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 H2 as a function of Ni content. [(a), (b), (d) are reprinted by permission from Macmillan Publishers Ltd: Nature Chemistry, Ref. 15, 2009, (c) is from ref. 25]

Close modal

In the remainder of this chapter, the background information that leads to the identification of appropriate catalyst descriptors (e.g. 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 her/his own pace.

Computational catalyst screening would not be possible without the existence of a theory that enables us to calculate 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:

  1. The ground state properties of a many-electron system are uniquely determined by the electron density.

  2. The total energy of a system has a minimum for the ground state electron density.

DFT provides a solution to the Schrödinger equation:

formula
Equation 1.1

and 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, EXC, 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

formula
Equation 1.2

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

formula
Equation 1.3

where veff is the effective field defined by the nuclei and the current electron density. 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, EXC, 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)30  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 (RPBE), in order to improve the accuracy of chemisorption energies of atoms and small molecules on transition metal surfaces.31  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+U, 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 impact on the quality of DFT calculations.

Although DFT calculations are at the heart of computational catalysis, it is not strictly necessary to perform your own calculations for a 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!32  But 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.33 

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

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. (Taken from ref. 20)

Figure 1.2

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. (Taken from ref. 20)

Close modal

Many of us have likely seen a schematic drawing of the bonding structure in a hydrogen molecule in one of our previous chemistry classes. Upon bond formation, the two atomic orbitals form two new molecular orbitals and they can be distinguished 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) we can imagine that they have a broad s-band (turquoise) and a narrow band of localized d-electrons (red). Assuming that we can separate the coupling between the adsorbate level and the s- and d-bands of the transition metal, we can write the chemisorption energy as the sum of both interactions.

ΔEChemESEd
Equation 1.4

The interaction of the adsorbate state with the broad s-band, ΔES, 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 ΔES, 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, we need to find a convenient way to define the term “position”, and therefore introduce the d-band center. The d-band center, , is the energy-weighted average of the density of d-states ρ. It sounds complicated, but it is the same equation that you use to calculate the center of mass (just substitute rCOM= , r=E and m=ρ).

formula
Equation 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, ΔEd. If we recall that ΔEs is constant, the variations in chemisorption energy across transition metals can be attributed to changes in ΔEd 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 references 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 (d9 and d10 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.36  For computational screening studies, however, we can safely neglect this exception, because it only affects the local fine structure and general trends are preserved.

If you have ever played the party game “Taboo” you have experienced how hard it can be to describe a 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 hidden 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 accurate enough 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 problem complexity 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 find their physical roots in the d-band model and are discussed next.

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 we can consider these bonds to form independently of each other and assume that the bond strength only depends on the two types of atom 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 we are only interested in the relative order between different transition metal catalysts, these errors are acceptable.

The discussion in this section is fully based on the scaling papers by Abild-Pedersen et al.37  and Fernandez et al.,38  which are highly recommended references on this topic. Each of the four panels in Figure 1.3 shows the adsorption energies, ΔEAHx, of the hydrogen-containing intermediates AHx (A=C, N, O, S) plotted as a function of the adsorption energy, ΔEA 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 OHx. It can be seen that the adsorption energies ΔEAHx are linearly correlated with ΔEA and given by

ΔEAHx=γΔEA+ξ.
Equation 1.6
Figure 1.3

Linear scaling relationships on close-packed terraces (black), step sites (red) and the FCC(100) surface (blue) of hydrogen-containing molecules of the type AHx with A=C, N, O, S. Reprinted with permission from F. Abild-Pedersen, J. Greeley, F. Studt, J. Rossmeisl, T. R. Munter, P. G. Moses, E. Skúlason, T. Bligaard and J. K. Nørskov, Phys. Rev. Lett., 2007, 99, 016104–016105. Copyright 2007 by the American Physical Society. Available at: http://link.aps.org/doi/10.1103/PhysRevLett.99.016105.37 

Figure 1.3

Linear scaling relationships on close-packed terraces (black), step sites (red) and the FCC(100) surface (blue) of hydrogen-containing molecules of the type AHx with A=C, N, O, S. Reprinted with permission from F. Abild-Pedersen, J. Greeley, F. Studt, J. Rossmeisl, T. R. Munter, P. G. Moses, E. Skúlason, T. Bligaard and J. K. Nørskov, Phys. Rev. Lett., 2007, 99, 016104–016105. Copyright 2007 by the American Physical Society. Available at: http://link.aps.org/doi/10.1103/PhysRevLett.99.016105.37 

Close modal

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?

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 adsorprtion energy, ΔEOH, which scales with a slope of ∼1/2 with respect to the adsorption energy of oxygen, ΔEO, on all three different surfaces. The same generalization applies to all the other adsorbates and can easily be rationalized using our 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 valences and likes to form four C–H bonds. If we remove these H atoms stepwise, 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. Now, as we move from one metal to another, we change the C–M bond strength and we can use the simple bond order method to estimate how this change affects the adsorption energy of the CHx species.

Figure 1.4

Simplified bonding scheme for adsorbed CHx species.

Figure 1.4

Simplified bonding scheme for adsorbed CHx species.

Close modal

In general, the scaling slope γ in eqn (1.6) for hydrogen-containing adsorbates AHx can then be predicted by

formula
Equation 1.7

where xmax is the maximum number of H atoms that can bind to the central atom A. In particular, xmax=4 for C, xmax=3 for N and xmax=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.38 

Now that we understand the origin of the slope in the linear scaling relations, we can focus our attention 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 we calculated ΔEAHx and ΔEA on a reference catalyst and we know xmax to predict the slope, γ. Then we can use eqn (1.6) to predict ΔEAHx as a function of ΔEA by simply eliminating the intercept ξ.

formula
Equation 1.8

Structure sensitivity is the other factor that affects the intercept, ξ, which is clearly demonstrated by the 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 Prof. Nørskov instead.39 

The special case of hydrogen-containing AHx 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, we may intuitively start 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 AHx and apply it to CHx–NHy with x,y=0–2 surface intermediates. These intermediates occur, for example, during the Degussa and Andrussow process for HCN synthesis from CH4 and NH3.17,18,40,41  We may assume that the C and N atoms are connected through a single bond and that the CHx–NHy intermediates bind to the surface equally through the C atom and the N atom. The adsorption energy of CHx–NHy should then scale linearly as in

ΔECHx−NHy=γ·(ΔECHxENHy)+ξ
Equation 1.9

where and . 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 we now substitute the expressions for ΔECHx and ΔENHy into eqn (1.9) and collect the intercepts, we arrive at our final expression for ΔE(CHx−NHy).

formula
Equation 1.10

The comparison of CHx–NHy adsorption energies estimated from eqn (1.10) 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 ca. 0.5 eV exist. A similar strategy based on bond order conservation was used by Jones et al. to derive scaling relations for C2 hydrocarbons on transition metal surfaces.42 

Figure 1.5

Linear scaling relations for complex adsorbates. (a) Parity plot of CHx–NHy adsorption energies predicted from scaling and the DFT reference calculations. (b) Linear scaling of C2H4 and C2H6 with methyl (CH3) from F. Studt et al., Science, 320, 1320–1322 reprinted with permission from AAAS.14 

Figure 1.5

Linear scaling relations for complex adsorbates. (a) Parity plot of CHx–NHy adsorption energies predicted from scaling and the DFT reference calculations. (b) Linear scaling of C2H4 and C2H6 with methyl (CH3) from F. Studt et al., Science, 320, 1320–1322 reprinted with permission from AAAS.14 

Close modal

Even without deriving scaling relations on paper as done for CHx–NHy, one can simply try different scaling relations and evaluate them on the basis of their mean absolute error and/or R2 value. Figure 1.5(b) shows a very good scaling relation for C2H2 and C2H4 with methyl (CH3), and many more such examples exist. Hence, if one has access to an adsorption energy database for many different surface species, 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.

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 universal and can be described by a single linear relationship!43–45  However, before we continue to discuss the beauty of transition state scaling it is important to clarify several technical terms used in the following discussion.

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 ETS and the energy of either the initial state (EIS) or final state (EFS) of a reaction. The transition state energy, ETS, is not equal to the activation barrier Ea of a reaction, but it can be used to calculate Ea=ETS−EIS. A BEP relation, on the other hand, is a relation between the activation energy Ea and the heat of a reaction ΔE, which is defined as ΔE=EFS−EIS. All energy quantities are indicated on a schematic potential energy surface (PES) in Figure 1.6.

Figure 1.6

A schematic potential energy surface for an arbitrary surface reaction. The nomenclature is Eref=reference energy chosen to be zero, EIS=initial state energy, ETS=transition state energy, EFS=final state energy, Ea=activation energy barrier, ΔE=energy change of the reaction.

Figure 1.6

A schematic potential energy surface for an arbitrary surface reaction. The nomenclature is Eref=reference energy chosen to be zero, EIS=initial state energy, ETS=transition state energy, EFS=final state energy, Ea=activation energy barrier, ΔE=energy change of the reaction.

Close modal

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.45  Since all reactions were formulated as dissociation reactions, Ediss was used to denote the final state energy EFS; and ΔEdiss 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) in both plots tells us 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

Transition state scaling (a) and BEP relation (b) for (de-)hydrogenation reactions. (Taken from ref. 45)

Figure 1.7

Transition state scaling (a) and BEP relation (b) for (de-)hydrogenation reactions. (Taken from ref. 45)

Close modal

To convince ourselves that, for most reactions, both graphs contain essentially the same information, we will derive a mathematical conversion of one graph into the other. We start by postulating that the initial and final state of the reaction is described by linear scaling relations that depend on the descriptors E1 and E2.

EIS1E11
Equation 1.11
EFS2E22
Equation 1.12

Using transition state scaling we can find ETS as a function of EFS.

ETS1EFS1
Equation 1.13

The quantity that is of practical interest for us is not ETS directly, but the activation energy barrier Ea, which is obtained in TSS from

graphic

formula
Equation 1.14

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

formula
Equation 1.15

Alternatively, Ea can be obtained from the BEP relation.

formula
Equation 1.16

Although eqns (1.15) and (1.16) look similar, they are not identical if the descriptors E1 and E2 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. If we assume a linear correlation and substitute

E23E13
Equation 1.17

into eqns (1.15) and (1.16), we can eliminate E2 and get

formula
Equation 1.18
formula
Equation 1.19

where the constant terms have been included in and . By equating and from eqns (1.18) and (1.19) it can easily be seen that they are identical when and . This shows that, as long as the descriptors E1 and E2 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 we have demonstrated that the choice between TSS and BEP is mostly a personal preference, the author wants to motivate his preference toward BEP over TSS relations. Let us assume that a transition state can be characterized as initial state like, final state like, or anything in between. We define ω 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, ETS, would then scale with EIS and EFS weighted by the transition state character ω. With β as an arbitrary energy offset, the energy ETS can then be expressed as

formula
Equation 1.20

Rewriting eqn (1.20) for the calculation of an activation barrier leads to a BEP relation relating Ea with ΔE and the slope of the BEP relation, α, can be interpreted as the transition state character, ω.

Ea=ETSEIS=ωΔE
Equation 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 Ea 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 we have used DFT to calculate Ea for a surface reaction at low coverage and we are able to predict the initial and final state energies for a high-coverage scenario. Performing additional DFT calculation for Ea at various higher surface coverages is prohibitively time consuming and an approximate method is desired to estimate the effect of coverage on Ea. In this case, we can apply the same idea used to derive eqn (1.8): we eliminate the constant β by using a reference calculation, and we can find Ea.

formula
Equation 1.22

Now that we can fully appreciate the meaning of TSS and BEP relations, we can 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.43  In Figure 1.7 it has already been shown that the universality principle holds for (de)hydrogenation reactions, and Figure 1.8(a) demonstrates 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/BEP relations for smaller families of reactions listed in refs. 44,45 should be used if possible. For preliminary screening studies, however, it is perfectly acceptable to assume universal TSS/BEP behavior and refine the results later, if needed.

Figure 1.8

(a) Universal TSS relation for C–C, C–O, C–N, N–O, N–N, and O–O dissociation reactions. With kind permission from Springer Science and Business Media,44  S. Wang et al., Cat. Lett., 2011, 141, 370–373. (b) Deviations from the “universal” TSS behavior are observed when the nature of the transition state switches between initial and final states alike. From Ref. 50 with permission from A. Vojvodic et al., J. Chem. Phys., 2011, 134, 244509. Copyright 2011, AIP Publishing LLC.

Figure 1.8

(a) Universal TSS relation for C–C, C–O, C–N, N–O, N–N, and O–O dissociation reactions. With kind permission from Springer Science and Business Media,44  S. Wang et al., Cat. Lett., 2011, 141, 370–373. (b) Deviations from the “universal” TSS behavior are observed when the nature of the transition state switches between initial and final states alike. From Ref. 50 with permission from A. Vojvodic et al., J. Chem. Phys., 2011, 134, 244509. Copyright 2011, AIP Publishing LLC.

Close modal

At this point it should be stressed that linear TSS/BEP relations exist for many reactions on a large range of catalyst surfaces, but they are not as “universal” as they have been introduced so far. Even for simple di-atomic dissociation reactions (N2, NO, O2) on metal-substituted La-perovskites a significant deviation from the universal TSS/BEP line has been observed.50  The deviation can clearly be seen in Figure 1.8(b) as one moves to more reactive surfaces (more negative ΔEdiss) 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 phenomena 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/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 very well be obtained.

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.51  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 we care primarily about 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

Comparison of the volcano curves obtained from a Sabatier analysis and full microkinetic models at different approaches to equilibrium γ. Reprinted from T. Bligaard, et al., J. Catal., 2004, 224, 206–217 with permission from Elsevier.21 

Figure 1.9

Comparison of the volcano curves obtained from a Sabatier analysis and full microkinetic models at different approaches to equilibrium γ. Reprinted from T. Bligaard, et al., J. Catal., 2004, 224, 206–217 with permission from Elsevier.21 

Close modal

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

  1. Write the corresponding rate expressions for all reaction steps in the forward direction as a function of the approach to equilibrium , where Ki is the equilibrium constant for step i, an is the activity of species n and vn is the corresponding stoichiometric coefficient.

  2. For each reaction step in the model, assume optimal surface coverages for the forward reaction.

  3. From the overall approach to equilibrium, calculate the approach to equilibrium of step i under the assumption that all other steps are quasi-equilibrated.

  4. Use TSS/BEP relations to calculate the forward rate constants ki and the rate. This rate represents an upper bound on the actual rate of step i and will be denoted as .

  5. The Sabatier rate is the slowest rate of all maximum rates in the mechanism: .

Let us illustrate the procedure 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):

(1)

(2)

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

formula
Equation 1.23
r2=k2PBθA(1−γ2)
Equation 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

formula
Equation 1.25

The allowed range for γ is 0≤γ≤1, where γ=1 indicates a quasi-equilibrated reaction. From eqn (1.25) we can estimate γ1 and γ2 under the assumption that the other step is quasi-equilibrated and obtain

formula
Equation 1.26

Assuming that we know k1 and k2 we can now write down the upper bounds for r1 and r2 and find the Sabatier rate rSabatier.

formula
Equation 1.27
formula
Equation 1.28
formula
Equation 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 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 one is better off using a microkinetic model if the simple Sabatier analysis fails. For additional information the reader is referred to refs. 52,54.

Now it is finally time to combine all the concepts introduced so far and perform our first real catalyst design. The reaction we will use for the design study is one of the most important heterogeneously catalyzed reactions, namely ammonia (NH3) synthesis from nitrogen (N2) and hydrogen (H2).

graphic

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 science55  and theoretical techniques,56  the detailed reaction mechanism has been elucidated and a complete potential energy surface (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) and stepped (solid) sites on Ru. K. Honkala et al., Science, 2005, 307, 555–558 reprinted with the permission of AAAS.56 

Figure 1.10

Potential energy surface for ammonia synthesis on terrace (dashed) and stepped (solid) sites on Ru. K. Honkala et al., Science, 2005, 307, 555–558 reprinted with the permission of AAAS.56 

Close modal

Upon inspection of the PES for NH3 synthesis one can simplify the mechanism to two main reaction steps, the dissociative adsorption of N2 followed by its hydrogenation to NH3. An extremely simplified reaction mechanism can therefore be written as

(1)

(2)

This mechanism neglects a number of elementary steps and does not include the possibility of adsorbed H atoms or NHx species, but it is sufficient to illustrate the Sabatier principle. A catalyst that associates strongly with N will easily dissociate N2, 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 H2 and NH3 are gas phase species). In contrast, a catalyst interacting weakly with N allows for easy desorption, but the barrier for N2 dissociation will be large. Hence, there is an optimal binding energy of N*, EN, at which the dissociative adsorption of N2 and NH3 desorption are optimized simultaneously. We will now try to find an estimate of EN, by performing a Sabatier analysis for typical NH3 synthesis conditions.

The simplified NH3 synthesis mechanism is almost identical to the generic mechanism used in the introduction of the Sabatier analysis in Section 1.4. If we substitute A2=N2 and B=3/2H2, we can directly use the previously obtained results and apply them to the current problem. This gives us the following two maximum rates for each step, where the stoichiometric factor of 3/2 for H2 is included as the exponent to in the expression for .

formula
Equation 1.30
formula
Equation 1.31

For the reaction conditions we may assume that we have a stoichiometric feed composition (H2:N2=3:1), the total pressure is P=100 bar and the temperature is T=673 K. Using this information we can calculate the partial pressures of hydrogen and nitrogen as =75 bar and PN2=25 bar. Further, we assume that the reaction is run at the limit of very low conversion where (γ→0). This leaves us with only the two unknown rate constants k1 and k2. We can obtain the rate constants, ki, from the standard Arrhenius expression, which depends on the activation energy, Ea,i, and the pre-exponential factor, νi, for each step.

formula
Equation 1.32

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

formula
Equation 1.33

Next, we obtain the activation barrier estimates from published scaling relations, beginning with N2 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.43 

Ea,1=0.87·ΔE1+1.34 eV
Equation 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 NH3, NH2, 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 our step (2), so we will have to make sure to use them accordingly. The estimated reverse activation energy is then given by

formula
Equation 1.35

and we find the forward barrier by using the fact that . The energy change of step (2), ΔE2, depends linearly on ΔE1 owing to the thermodynamic constraint that ½ΔE1+ΔE2=ΔEr, where ΔEr is the overall energy change of the reaction. So far we have exclusively used the electronic ground-state energy, E, at T=0 K in all equations, which is the primary output of DFT calculations, but from thermodynamics we know that free energy changes are the correct quantities to use in the calculation of rate and equilibrium constants. We assume 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 we can say that ΔEr ≈ ΔH°=−45.9 kJ mol−1=–0.48 eV. Putting this all together we get the desired BEP relation for Ea,2.

formula
Equation 1.36

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

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, we may assume the entropy to be independent of the catalyst. For a rough approximation, we can further assume that an adsorbed species or transition state has lost all degrees of freedom (it is completely locked in place on the surface), and the entropy is therefore zero. The gas-phase entropies of N2 and H2 at standard state are S°(N2)=192.77 J mol−1 K−1 and S°(H2) 130.68 J mol−1 K−1, respectively, and we can now estimate the entropy changes between the initial states and the transition state as and . For the first step it is assumed that N2 loses all its gas-phase entropy upon entering the transition state and in the second step, the entropy loss of one H2 gas-phase molecule is considered when it reacts with N* on the surface. Why did we consider the entropy loss of only one H2 molecule and not 3/2H2 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 we know that it is a sequence of hydrogen adsorption and hydrogenation steps and for each H2 molecule adsorbing onto the surface –130.68 J mol−1 K−1of 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 ν1=1.20×103 s−1 and ν2=2.10×106 s−1.

At this point the reader may want to open her/his favorite math and graphing software (Excel would suffice), enter the eqns (1.30)–(1.32), (1.34), and (1.36), and calculate the maximum N2 dissociation rate and maximum desorption rate as a function of ΔE1=ΔEN 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 ΔEN≈−0.95 eV and has a predicted turn over frequency (TOF) of 8.7 s−1. 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 we can finally proceed to the last step of the example problem and screen for the best catalyst. For this purpose, let us examine the “database” of dissociative chemisorption energies, ΔEN, of several transition metals listed in Table 1.1, and we find that 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 ΔEN spans a range from ca. –4 to +6 eV and neighboring transition metals are often separated by 0.5 eV or more.

Table 1.1

Excerpt of the periodic table with dissociative chemisorption energies of N2EN) given in eV. Reprinted from T. Bligaard, et al., J. Catal., 2004, 224, 206–217 with permission from Elsevier.21 

Cr Mn Fe Co Ni Cu 
– – −1.27 −0.38 −0.10 2.88 
Mo Tc Ru Rh Pd Ag 
−2.76 – −0.84 −0.70 1.78 5.86 
Re Os Ir Pt Au 
−4.33 – – −0.59 1.37 5.89 
Cr Mn Fe Co Ni Cu 
– – −1.27 −0.38 −0.10 2.88 
Mo Tc Ru Rh Pd Ag 
−2.76 – −0.84 −0.70 1.78 5.86 
Re Os Ir Pt Au 
−4.33 – – −0.59 1.37 5.89 
Figure 1.11

Sabatier volcano for NH3 synthesis.

Figure 1.11

Sabatier volcano for NH3 synthesis.

Close modal

In summary, this example of computational screening for an NH3 synthesis catalyst has demonstrated how BEP relations are applied to reduce the number of parameters to a single descriptor, ΔEN, and how to perform a Sabatier analysis. As a result we obtained the Sabatier volcano shown in Figure 1.12 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 in Section 1.6.2.

Figure 1.12

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

Figure 1.12

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

Close modal

The Sabatier analysis for NH3 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 we can derive a volcano surface 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.53  Here, we discuss only the high temperature CO oxidation with the reaction conditions T=600 K, PO2=0.33 bar, PCO=0.67 bar. The low temperature reaction can be studied analogously. For the CO oxidation reaction we start with the following reaction mechanism.

(1)

(2)

(3)

(4)

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

formula
Equation 1.37
formula
Equation 1.38

The assumed optimal coverages that maximize are , and is maximized for θCOO=0.5, leading to our final expression of the Sabatier rate.

rSabatier=min{0.5·k3,0.25·k4}
Equation 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 ν=1013 s−1. The underlying assumption is that the entropy of the transition state is similar to the entropy of the initial state (ΔS‡=0) and eqn (1.33) reduces to 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 fcc(111) surfaces.53 

ΔETS,3=1.39·ΔEO+1.56 eV
Equation 1.40
ΔETS,4=0.7·(ΔECOEO)+0.02 eV
Equation 1.41
formula
Equation 1.42

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

formula
Equation 1.43
ΔEa,4ETS,4−ΔEO−ΔEO=−0.3·(ΔEOECO)+0.02 eV
Equation 1.44

The Sabatier rate for CO oxidation in eqn (1.39) can now be obtained from the knowledge of ΔECO and ΔEO, which have naturally evolved again as the two descriptors that we need 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 studying 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 ΔECO=–1.20 eV. This value is within ±0.3 eV of the actual CO binding energies on 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 O2 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 we expected. Given the similarly correct identification of the optimum catalysts for NH3 synthesis in the earlier example, the success of a Sabatier analysis can no longer be considered a coincidence!

Figure 1.13

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

Figure 1.13

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

Close modal

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 we continue 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.,57  Stoltze and Nørskov,58  and the standard reference The Microkinetics of Heterogeneous Catalysis by Dumesic et al.59  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.60  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:

formula
Equation 1.45
formula
Equation 1.46
formula
Equation 1.47
formula
Equation 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.

formula
Equation 1.49
formula
Equation 1.50
formula
Equation 1.51

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

formula
Equation 1.52

Eqns (1.49)–(1.51) represent a system of ordinary differential equations (ODE), but we will restrict ourselves to situations where the reaction runs under steady-state conditions, which means that the derivatives of the surface coverages with respect to time in eqns (1.49)–(1.51) are all zero.

formula
Equation 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 eqns (1.52) and (1.53) there are four equations for the four unknown coverages , θ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/desorption steps are fast in comparison to surface reactions and may be assumed to be quasi-equilibrated. We can make this assumption for steps (1) and (2) of our mechanism and then find the surface coverages of CO* and O2* as a function of the fraction of empty sites.

formula
Equation 1.54
formula
Equation 1.55

There is not a reaction that we can assume to be quasi-equilibrated to derive an equation for the coverage of O*, but from the steady-state assumption and eqn (1.51) we find

formula
Equation 1.56

which can be solved for θO by substituting and θCO from (1.54) and (1.55).

formula
Equation 1.57

In eqn (1.57) we introduced W to denote the long fraction. For an analytical solution one now simply uses the site balance eqn (1.52) to solve for θ*.

formula
Equation 1.58

Equation (1.58) can be solved for any reaction condition if the equilibrium constants, Ki, for steps (1) and (2), and the forward and reverse rate constants, and , for steps (3) and (4) are known. Note that we do not need to know the rate constants of the quasi-equilibrated steps (1) and (2) and the equilibrium constants are sufficient. We will use this analytical solution later and compare it to a full numerical solution of the CO oxidation model.

Most people will tell you 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 we do a computational screening study we want to explore a large parameter space for our 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. 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 to 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

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

Figure 1.14

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

Close modal

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 ensure 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 pay attention! If you slow down these steps 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.

In the previous examples we only considered electronic energy changes and approximated the entropy as all or nothing. In essence, we 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 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, Eelec, the zero-point energy, EZPE, and a heat capacity correction to account for the difference between the enthalpy at 0 K and at the actual temperature, T.

formula
Equation 1.59

The electronic energy component, Eelec, 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.61 

graphic
graphic

formula
Equation 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.

formula
Equation 1.61

In this equation m denotes the mass of the particle and CS is the number of sites per unit area, which is roughly 1015 sites cm−2 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.

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).

The degree of rate control, XRC, was originally proposed by Campbell as a quantitative measure to identify rate-controlling steps.62,63  It was later generalized to quantify also the impact of adsorbate binding (poisoning) on the reaction rate.64  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, ki, while keeping the equilibrium constant, Ki, and the rate constants, kj, of all other steps constant.

formula
Equation 1.62

An XRC,i value of zero indicates a quasi-equilibrated step, whose rate constant has no effect on the overall rate, while XRC,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 , 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, eqn (1.62) is not easily applicable, but it can be rewritten into a finite-difference equation that can be readily evaluated in a microkinetic model.

formula
Equation 1.63

In practice XRC,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, r0, for the original values of ki. Never forget to apply this change to both rate constants, otherwise Ki and the overall equilibrium constant, Keq, will be altered as well!

Fundamentally, the definition of the degree of rate control in eqn (1.62) implies a change in the Gibbs free energy of the transition state of step i, as indicated for the transition state of CH4 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

formula
Equation 1.64
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 with the kind permission of Science, ref. 65)

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 with the kind permission of Science, ref. 65)

Close modal

where c is a constant ( in transition state theory) and Ga,i is the free energy of activation for step i. Upon substituting eqn (1.64) into the definition of XRC,i in eqn (1.62) we obtain the degree of rate control in terms of the free energy of the transition state of step i, (GTS,i).

formula
Equation 1.65

Equivalently, one can change the free energy Gn 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 one can define the thermodynamic degree of rate control, XTRC,n, for an intermediate n as

formula
Equation 1.66

Eqns (1.65) and (1.66) 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.

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 we can find catalysts 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 XCC, which is a constrained sensitivity analysis of a microkinetic model with respect to a reactivity descriptor Ei.65  The definition of XCC,i is analogous to the generalized degree of rate control and may be expressed as

formula
Equation 1.67

In eqn (1.67) it is indicated that all other descriptor values Ej≠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 XCC will indicate the direction one must look to identify an improved catalyst. The optimal catalyst is characterized by XCC,i=0 for all descriptors Ei. XCC 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 XRC and XCC is shown graphically in Figure 1.15 for the sensitivity towards the binding energy of carbon, ΔEC. Carbon adsorbs quite strongly and is located in a well of the PES, which explains the negative value (–0.26) for XRC(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*, CH2*, CH3* 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 CH4 activation will shift up in energy. The large value of XRC(CH4-TS)=0.8 indicates that CH4 activation is a rate-determining step and increasing the barrier will have a large negative impact on the overall reaction rate. The overall effect of changing the descriptor ΔEC on the PES is given by Xcc(C)=0.11, indicating that a catalyst with stronger ΔEC will have higher activity. The explanation is that the increased poisoning effect of C* is more than compensated by lowering the barrier for the CH4 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. XRC should be used to analyze the reaction mechanism and understand the importance of individual steps and adsorbates, whereas XCC should be used to guide catalyst design.

We have already performed a preliminary Sabatier analysis of the CO oxidation reaction in Section 1.5.2, and derived an analytic solution under the assumption that the adsorption of CO and O2 are quasi-equilibrated in Section 1.6. Now we will formulate a numerical solution to the complete microkinetic model as a function of the descriptors ΔECO and ΔEO. We will analyze the reaction mechanism 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.

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. an ODE solver, a root-finder, and n-dimensional curve fitting routines, if parameter optimization is also required. Here we use python, which is an object-oriented scripting language with very 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, we also need 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 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.

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.

graphic

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 we will need this procedure multiple times for each combination of our descriptors, we define it as a function.

graphic

In this function we have used the scaling/BEP relations, eqns (1.42)–(1.44), and we have arbitrarily assumed that CO and O2 lose 25% of their gas-phase entropy when entering the transition state for adsorption. Using collision theory for the adsorption rate constants66  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, Ea,i, are all positive, even for descriptor pairs that would otherwise predict negative values for Ea,i. Second, thermodynamic consistency is enforced throughout the model by calculating the reverse rate constants from the equilibrium constants for each step.

This step is simple. We need rate equations for the individual steps. 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.

graphic

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 you want to use, 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 we named the previously defined function get_rates().

graphic

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×103 to 1×108 seconds. Some might wonder why a reaction could possibly take 1×108 seconds, or greater than 3 years, to equilibrate, but keep in mind that our goal is to scan a wide descriptor range, which will result in some very small rate constants. For this example, we use the Pt values for the descriptors and the odeint ODE solver provided by scipy. Note that the default tolerance and step size had to be tweaked in order for odeint to work for this stiff problem.

graphic

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, we are only interested in the last row, conveniently accessed by using ‘−1’ as index, but the time evolution is available as well for transient studies or verification that steady-state has been reached.

graphic

As a result of this exercise, we find the steady-state coverages on Pt of θCO=0.215 ML, θO2=8.4×10−4, and θO=5.6×10−3 and a steady-state CO2 formation rate of 6.16×103 s−1.

Now that we are able 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, we need to solve the microkinetic model repeatedly for different parameters and it is more convenient if we first write a function, solve_ode, that takes the rate constants as input and returns the solution of the model. This function may look like this.

graphic

The degree of rate control, XRC, is then calculated by repeatedly calling solve_ode with rate constants that are systematically changed for each step. As input to calculate_Xrc we provide the unaltered rate constants and the corresponding rate.

graphic

In our example for CO oxidation on Pt with ΔEO=–1.25 eV and ΔECO=−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 XRC,3=0.99, and step (4) has an XRC,4=0.01. In this example , which is expected for a serial reaction mechanism.

Next we calculate the degree of catalyst control Xcc by varying the values of ΔEO and ΔECO 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, we pass the original values of ΔEO and ΔECO to the function calculate_Xcc and the corresponding rate constants are calculated on the fly.

graphic

If we perform this analysis for CO oxidation on Pt, we get Xcc(ΔEO)=1.38 and Xcc(ΔECO)=–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 ΔECO 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, eqns (1.54)–(1.58). The positions of the closed-packed surfaces of Pt, Pd, Cu, and Rh are indicated.

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, eqns (1.54)–(1.58). The positions of the closed-packed surfaces of Pt, Pd, Cu, and Rh are indicated.

Close modal

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 [eqns (1.54)–(1.58)] is given in Figure 1.16(b). If we recall that the only assumption that was made in the derivation of the analytical solution was that the CO and O2 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 to the next one. However, even this approach is not fail 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 one must resort to a slower, but more stable, solver in order to generate a new solution and then possibly continue again with the fast root-finding method. An example of how this may be implemented is provided below.

graphic

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.

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 ca. 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.68  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 pursuits 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.67 

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 ΔEO and ΔECO could be identified, the volcano in Figure 1.17(a) was calculated under 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, you may 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 neither alter the position of the volcano peak nor change the relative order of the activity of metal catalysts.

Figure 1.17

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. L. C. Grabow, et al., ‘Understanding Trends in Catalytic Activity: The Effect of Adsorbate Adsorbate Interactions for CO Oxidation Over Transition Metals’, Topic. Catal., 2010, 53, 298–310 is given with kind permission from Springer Science and Business Media.

Figure 1.17

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. L. C. Grabow, et al., ‘Understanding Trends in Catalytic Activity: The Effect of Adsorbate Adsorbate Interactions for CO Oxidation Over Transition Metals’, Topic. Catal., 2010, 53, 298–310 is given with kind permission from Springer Science and Business Media.

Close modal

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 his/her own. Although the field is still relatively new, it has already reached such a high level of complexity that several topics could only be mentioned in passing and references for further study 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.

Let us conclude this chapter 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 her/his perfect partner (= catalyst) will first have to go through an in-depth interview. During the interview, the matchmaker thoroughly analyzes the client's character/personal situation and describes it 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, he or she 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 his/her 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…”!

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

graphic

1.
Yaccato
 
K.
Carhart
 
R.
Hagemeyer
 
A.
Herrmann
 
M.
Lesik
 
A.
Strasser
 
P.
Volpe
 
A.
Turner
 
H.
Weinberg
 
H.
Grasselli
 
R. K.
Brooks
 
C. J.
Pigos
 
J. M.
Combin. Chem. High Throughput Screen.
2010
, vol. 
13
 (pg. 
318
-
330
)
2.
Greeley
 
J.
Jaramillo
 
T. F.
Bonde
 
J.
Chorkendorff
 
I. B.
Nørskov
 
J. K.
Nature Mat.
2006
, vol. 
5
 (pg. 
909
-
913
)
3.
Mavrikakis
 
M.
Nature Mat.
2006
, vol. 
5
 (pg. 
847
-
848
)
4.
Somorjai
 
G. A.
Li
 
Y.
Top. Catal.
2010
, vol. 
53
 (pg. 
311
-
325
)
5.
Besenbacher
 
F.
Chorkendorff
 
I.
Clausen
 
B. S.
Hammer
 
B.
Molenbroek
 
A. M.
Nørskov
 
J. K.
Stensgaard
 
I.
Science
1998
, vol. 
279
 pg. 
1913
 
6.
Jacobsen
 
C. J.
Dahl
 
S.
Clausen
 
B. S.
Bahn
 
S.
Logadottir
 
A.
Nørskov
 
J. K.
J. Am. Chem. Soc.
2001
, vol. 
123
 (pg. 
8404
-
8405
)
7.
Toulhoat
 
H.
J. Catal.
2003
, vol. 
216
 (pg. 
63
-
72
)
8.
Strasser
 
P.
Fan
 
Q.
Devenney
 
M.
Weinberg
 
W. H.
Liu
 
P.
Nørskov
 
J. K.
J Phys. Chem. B
2003
, vol. 
107
 (pg. 
11013
-
11021
)
9.
Nilekar
 
A. U.
Sasaki
 
K.
Farberow
 
C. A.
Adzic
 
R. R.
Mavrikakis
 
M.
J. Am. Chem. Soc.
2011
, vol. 
133
 (pg. 
18574
-
18576
)
10.
Linic
 
S.
Jankowiak
 
J.
Barteau
 
M. A.
J. Catal.
2004
, vol. 
224
 (pg. 
489
-
493
)
11.
Greeley
 
J.
Mavrikakis
 
M.
Nature Mat.
2004
, vol. 
3
 (pg. 
810
-
815
)
12.
Jaramillo
 
T. F.
Jørgensen
 
K. P.
Bonde
 
J.
Nielsen
 
J. H.
Horch
 
S.
Chorkendorff
 
I.
Science
2007
, vol. 
317
 (pg. 
100
-
102
)
13.
Alayoglu
 
S.
Nilekar
 
A. U.
Mavrikakis
 
M.
Eichhorn
 
B.
Nature Mat.
2008
, vol. 
7
 (pg. 
333
-
338
)
14.
Studt
 
F.
Abild-Pedersen
 
F.
Bligaard
 
T.
Sorensen
 
R. Z.
Christensenand
 
C. H.
Nørskov
 
J. K.
Science
2008
, vol. 
320
 (pg. 
1320
-
1322
)
15.
Greeley
 
J.
Stephens
 
I. E. L.
Bondarenko
 
A. S.
Johansson
 
T. P.
Hansen
 
H. A.
Jaramillo
 
T. F.
Rossmeisl
 
J.
Chorkendorff
 
I.
Nørskov
 
J. K.
Nature Chem.
2009
, vol. 
1
 (pg. 
552
-
556
)
16.
Zhang
 
J.
Vukmirovic
 
M. B.
Sasaki
 
K.
Nilekar
 
A. U.
Mavrikakis
 
M.
Adzic
 
R. R.
J. Am. Chem. Soc.
2005
, vol. 
127
 (pg. 
12480
-
12481
)
17.
Grabow
 
L. C.
Studt
 
F.
Abild-Pedersen
 
F.
Petzold
 
V.
Kleis
 
J.
Bligaard
 
T.
Nørskov
 
J. K.
Ange. Chem. Int. Ed.
2011
, vol. 
50
 (pg. 
4601
-
4605
)
18.
Grabow
 
L. C.
Studt
 
F.
Abild-Pedersen
 
F.
Petzold
 
V.
Kleis
 
J.
Bligaard
 
T.
Nørskov
 
J. K.
Ange. Chem.
2011
, vol. 
123
 (pg. 
4697
-
4701
)
19.
Nørskov
 
J. K.
Bligaard
 
T.
Rossmeisl
 
J.
Christensen
 
C. H.
Nature Chem.
2009
, vol. 
1
 (pg. 
37
-
46
)
20.
Nørskov
 
J. K.
Abild-Pedersen
 
F.
Studt F
 
F.
Bligaard
 
T.
Proc. Natl. Acad. Sci. USA
2011
, vol. 
108
 (pg. 
937
-
943
)
21.
Bligaard
 
T.
Nørskov
 
J. K.
Dahl
 
S.
Matthiesen
 
J.
Christensen
 
C. H.
Sehested
 
J.
J. Catal.
2004
, vol. 
224
 (pg. 
206
-
217
)
22.
Andersson
 
M. P.
Bligaard
 
T.
Kustov
 
A.
Larsen
 
K. E.
Greeley
 
J.
Johannessen
 
T.
Christensen
 
C. H.
Nørskov
 
J. K.
J. Catal.
2006
, vol. 
239
 (pg. 
501
-
506
)
23.
Sehested
 
J.
Larsen
 
K.
Kustov
 
A.
Frey
 
A.
Johannessen
 
T.
Bligaard
 
T.
Andersson
 
M.
Nørskov
 
J.
Christensen
 
C.
Topic. Catal.
2007
, vol. 
45
 (pg. 
9
-
13
)
24.
M.
Andersson
,
T.
Bligaard
,
C. H.
Christensen
,
T.
Johannessen
,
A.
Kustov
,
K. E.
Larsen
,
J. K.
Nørskov
and
J.
Sehested
, “Process and Catalyst For Hydrogenation of Carbon Oxides”, US Patent US7790776 B2 (
2010
)
25.
Munter
 
T. R.
Landis
 
D. D.
Abild-Pedersen
 
F.
Jones
 
G.
Wang
 
S.
Bligaard
 
T.
Comp. Sci. Discov.
2009
, vol. 
2
 pg. 
015006
 
26.
Kohn
 
W.
Sham
 
L. J.
Phys. Rev.
1965
, vol. 
140
 pg. 
A1133
 
27.
Hohenberg
 
P.
Kohn
 
W.
Phys. Rev.
1964
, vol. 
136
 pg. 
B864
 
28.
Perdew
 
J. P.
Chevary
 
J. A.
Vosko
 
S. H.
Jackson
 
K. A.
Pederson
 
M. R.
Singh
 
D. J.
Fiolhais
 
C.
Phys. Rev. B
1992
, vol. 
46
 pg. 
6671
 
29.
Perdew
 
J. P.
Chevary
 
J. A.
Vosko
 
S. H.
Jackson
 
K. A.
Pederson
 
M. R.
Singh
 
D. J.
Fiolhais
 
C.
Phys. Rev. B
1993
, vol. 
48
 pg. 
4978
 
30.
Perdew
 
J. P.
Burke
 
K.
Ernzerhof, M.
 
M.
Phys. Rev. Lett.
1996
, vol. 
77
 (pg. 
3865
-
3868
)
31.
Hammer
 
B.
Hansen
 
L.
Nørskov
 
J.
J. Phys. Rev. B
1999
, vol. 
59
 (pg. 
7413
-
7421
)
32.
Hummelshøj
 
J. S.
Abild-Pedersen
 
F.
Studt
 
F.
Bligaard
 
T.
Nørskov
 
J. K.
Ange. Chem. Int. Ed. Engl.
2011
, vol. 
51
 (pg. 
272
-
274
)
33.
D. S.
Sholl
and
J. A.
Steckel
,
Density Functional Theory: A Practical Introduction
,
Wiley-VCH
,
Verlag
,
2009
, vol. 135
34.
Hammer
 
B.
Nørskov
 
J. K.
Adv. Catal.
2000
, vol. 
45
 (pg. 
71
-
129
)
35.
Hammer
 
B.
Nørskov
 
J. K.
Surf. Sci.
1995
, vol. 
343
 (pg. 
211
-
220
)
36.
Xin
 
H.
Linic
 
S.
J. Chem.Phys.
2010
, vol. 
132
 (pg. 
221101
-
221104
)
37.
Abild-Pedersen
 
F.
Greeley
 
J.
Studt
 
F.
Rossmeisl
 
J.
Munter
 
T. R.
Moses
 
P. G.
Skúlason
 
E.
Bligaard
 
T.
Nørskov
 
J. K.
Phys. Rev. Lett.
2007
, vol. 
99
 (pg. 
016104
-
016105
)
38.
Fernandez
 
E. M.
Moses
 
P. G.
Toftelund
 
A.
Hansen
 
H. A.
Martinez
 
J. I.
Abild-Pedersen
 
F.
Kleis
 
J.
Hinnemann
 
B.
Rossmeisl
 
J.
Bligaard
 
T.
Nørskov
 
J. K.
Ange. Chem.-Int. Ed.
2008
, vol. 
47
 (pg. 
4683
-
4686
)
39.
Nørskov
 
J. K.
Bligaard
 
T.
Hvolbaek
 
B.
Abild-Pedersen
 
F.
Chorkendorff
 
I.
Christensen
 
C. H.
Chem. Soc. Rev.
2008
, vol. 
37
 (pg. 
2163
-
2171
)
40.
Gómez-Díaz
 
J.
López
 
N.
J. Phys. Chem. C
2011
, vol. 
115
 (pg. 
5667
-
5674
)
41.
Gómez-Díaz
 
J.
Vargas-Fuentes
 
C.
López
 
N.
J. Chem. Phys.
2011
, vol. 
135
 pg. 
124707
 
42.
Jones
 
G.
Studt
 
F.
Abild-Pedersen
 
F.
Nørskov
 
J. K.
Bligaard
 
T.
Chem. Engin. Sci.
2011
, vol. 
66
 (pg. 
6318
-
6323
)
43.
Nørskov
 
J. K.
Bligaard
 
T.
Logadottir
 
A.
Bahn
 
S.
Hansen
 
L. B.
Bollinger
 
M.
Bengaard
 
H.
Hammer
 
B.
Sljivancanin
 
Z.
Mavrikakis
 
M.
Xu
 
Y.
Dahl
 
S.
Jacobsen
 
C. J. H.
J. Catal.
2002
, vol. 
209
 (pg. 
275
-
278
)
44.
Wang
 
S.
Temel
 
B.
Shen
 
J.
Jones
 
G.
Grabow
 
L. C.
Studt
 
F.
Bligaard
 
T.
Abild-Pedersen
 
F.
Christensen
 
C. H.
Nørskov
 
J. K.
Catalysis Letters
2011
, vol. 
141
 (pg. 
370
-
373
)
45.
Wang
 
S.
Petzold
 
V.
Tripkovic
 
V.
Kleis
 
J.
Howalt
 
J. G.
Skúlason
 
E.
Fernández
 
E. M.
Hvolbæig;k
 
B.
Jones
 
G.
Toftelund
 
A.
Falsig
 
H.
Björketun
 
M.
Studt
 
F.
Abild-Pedersen
 
F.
Rossmeisl
 
J.
Nørskov
 
J. K.
Bligaard
 
T.
Phys. Chem. Chem. Phys.
2011
, vol. 
13
 (pg. 
20760
-
20765
)
46.
Brønsted
 
J. N.
Chem. Rev.
1928
, vol. 
5
 (pg. 
231
-
338
)
47.
Evans
 
M. G.
M.
 
Polanyi
Trans. Faraday Soc.
1938
, vol. 
34
 (pg. 
11
-
24
)
48.
Grabow
 
L. C.
Gokhale
 
A. A.
Evans
 
S. T.
Dumesic
 
J. A.
Mavrikakis
 
M.
J. Phys. Chem. C
2008
, vol. 
112
 (pg. 
4608
-
4617
)
49.
Grabow
 
L. C.
Hvolbæig;k
 
B.
Nørskov
 
J. K.
Topic. Catal.
2010
, vol. 
53
 (pg. 
298
-
310
)
50.
Vojvodic
 
A.
Calle-Vallejo
 
F.
Guo
 
W.
Wang
 
S.
Toftelund
 
A.
Studt
 
F.
Martínez
 
J. I.
Shen
 
J.
Man
 
I. C.
Rossmeisl
 
J.
Bligaard
 
T.
Noørskov
 
J. K.
Abild-Pedersen
 
F.
J. Chem. Phys.
2011
, vol. 
134
 pg. 
244509
 
51.
Sabatier
 
P.
Berichte Deutsch. Chem. Gesells.
1911
, vol. 
44
 (pg. 
1984
-
2001
)
52.
Falsig
 
H.
Bligaard
 
T.
Rass-Hansen
 
J.
Kustov
 
A. L.
Christensen
 
C. H.
Nørskov
 
J. K.
Topic. Catal.
2007
, vol. 
45
 (pg. 
117
-
120
)
53.
Falsig
 
H.
Hvolbæig;k
 
B.
Kristensen
 
I. S.
Jiang
 
T.
Bligaard
 
T.
Christensen
 
C. H.
Nørskov
 
J. K.
Ange. Chem.-Int. Ed.
2008
, vol. 
47
 (pg. 
4835
-
4839
)
54.
H.
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
55.
Ertl
 
G.
Ange. Chem. Int. Ed.
2008
, vol. 
47
 (pg. 
3524
-
3535
)
56.
Honkala
 
K.
Hellman
 
A.
Remediakis
 
I. N.
Logadottir
 
A.
Carlsson
 
A.
Dahl
 
S.
Christensen
 
C. H.
Nørskov
 
J. K.
Science
2005
, vol. 
307
 (pg. 
555
-
558
)
57.
Gokhale
 
A. A.
Kandoi
 
S.
Greeley
 
J. P.
Mavrikakis
 
M.
Dumesic
 
J. A.
Chem. Engin. Sci.
2004
, vol. 
59
 (pg. 
4679
-
4691
)
58.
P.
Stoltze
and
J. K.
Nørskov
,
“Theoretical Modelling of Catalytic Reactions”
in
Handbook of Heterogeneous Catalysis
, 2nd Edition, Eds.: G. Ertl, H. Knözinger, F. Schüth and J. Weitkamp,
Wiley-VCH Verlag GmbH & Co. KGaA
,
Weinheim, Germany
,
2008
59.
J. A.
Dumesic
,
D. F.
Rudd
,
L. L.
Aparicio
,
J. E.
Rekoske
and
A. A.
Treviño
,
The Microkinetics of Heterogeneous Catalysis
,
American Chemical Society
,
Washington, DC
,
1993
60.
Freund
 
H.-J.
Meijer
 
G.
Scheffler
 
M.
Schlögl
 
R.
Wolf
 
M.
Ange. Chem. Int. Ed.
2011
, vol. 
50
 (pg. 
10064
-
10094
)
61.
C. J.
Cramer
,
Essentials of Computational Chemistry: Theories and Models
, 2nd edn,
John Wiley & Sons Ltd.
,
Chichester, England
,
2004
62.
Campbell
 
C. T.
Topic. Catal.
1994
, vol. 
1
 (pg. 
353
-
366
)
63.
Campbell
 
C. T.
J. Catal.
2001
, vol. 
204
 (pg. 
520
-
524
)
64.
Stegelmann
 
C.
Andreasen
 
A.
Campbell
 
C. T.
J. Am. Chem. Soc.
2009
, vol. 
131
 (pg. 
8077
-
8082
)
65.
Nørskov
 
J. K.
Bligaard
 
T.
Kleis
 
J.
Science
2009
, vol. 
324
 (pg. 
1655
-
1656
)
66.
Cortright
 
R. D.
Dumesic
 
J. A.
Adv. Catal.
2001
, vol. 
46
 pg. 
161
 
67.
Hellman
 
A.
Honkala
 
K.
J. Chem. Phys.
2007
, vol. 
127
 pg. 
194704
 
68.
İnoĝlu
 
N.
Kitchin
 
J.
Phys. Rev. B
2010
, vol. 
82
 (pg. 
1
-
5
)
69.
Getman
 
R. B.
Schneider
 
W. F.
ChemCatChem
2010
, vol. 
2
 (pg. 
1450
-
1460
)
70.
Wu
 
C.
Schmidt
 
D. J.
Wolverton
 
C.
Schneider
 
W. F.
J. Catal.
2012
, vol. 
286
 (pg. 
88
-
94
)
Close Modal

or Create an Account

Close Modal
Close Modal