 1.1 Introduction
 1.1.1 A Walk through a Computational Catalyst Design Process: Methanation
 1.2 Starting from the Electronic Structure
 1.2.1 Density Functional Theory
 1.2.2 The dBand Model
 1.3 Identifying the Right Descriptor Set
 1.3.1 Scaling Relations for Surface Intermediates
 1.3.2 Scaling Relations for Transition States: The Brønsted–Evans–Polanyi Relationship
 1.4 The Sabatier Principle and the Volcano Curve
 1.4.1 Sabatier Analysis
 1.5 Sabatier Analysis in Practice
 1.5.1 First Example: Ammonia Synthesis
 1.5.2 Second Example: CO Oxidation
 1.6 Notes on Microkinetic Modeling
 1.6.1 Numerical Solution Strategies
 1.6.2 Entropy and Enthalpy Corrections
 1.6.3 Microkinetic Model Analysis
 1.7 CO Oxidation Catalyst Screening
 1.7.1 Numerical Microkinetic Model
 1.7.2 Degree of Rate and Catalyst Control
 1.7.3 Twodimensional CO Oxidation Volcano
 1.7.4 Effect of Lateral Interactions
 1.8 Conclusions
 Appendix
CHAPTER 1: Computational Catalyst Screening

Published:02 Dec 2013

Series: Catalysis Series
L. C. Grabow, in Computational Catalysis, ed. A. Asthagiri and M. J. Janik, The Royal Society of Chemistry, 2013, pp. 158.
Download citation file:
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 stepbystep instructions with solved numerical problems for NH_{3} synthesis and CO oxidation on transition metal surfaces. The theoretical foundation of the screening approach, including the dband 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.
1.1 Introduction
Bruteforce 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 bruteforce 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 128bit encryption key has 2^{128} possible permutations. If we simply assume that a typical central processing unit (CPU) can generate 10^{9} bit flips per second (∼1 GHz), then the total time that is required to test all possible permutations is 2^{128}/10^{9}=3.4×10^{29} seconds or 10^{22} years! For obvious reasons a bruteforce 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 bruteforce 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, highthroughput 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 trialanderror 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 highthroughput 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 firstprinciples density functional theory (DFT) simulations and stored in a large property database. Populating this property database with DFT data is the most timeconsuming 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 bruteforce 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 descriptorbased 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 hydrodesulfurization activity on metal–sulfide catalysts.^{7 } These initial successes were followed by other prominent examples that include COtolerant fuel cell anodes,^{8,9 } Cu/Ag alloys as selective ethylene epoxidation catalysts,^{10 } nearsurface alloys for hydrogen activation^{11 } and evolution^{2,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 electrodes^{15 } and mixedmetal Pt monolayer catalysts^{16 } for electrochemical oxygen reduction, and the rediscovery of Pt as the most active and selective catalyst for the production of hydrogen cyanide.^{17,18 }
1.1.1 A Walk through a Computational Catalyst Design Process: Methanation
The most comprehensive example of a success story in computational catalyst design comes from the group of Jens Nørskov, who has pioneered the descriptorbased design approach and has applied it to numerous reactions.^{19,20 } In several publications his group has studied the methanation reaction (CO+2H_{2} → CH_{4}+H_{2}O), starting from a detailed electronic structure analysis and leading to the development of a patented technical methanation catalyst based on a Fe/Ni alloy.^{21–25 } In the beginning of any descriptorbased 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 ratelimiting 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 paretooptimal set. The paretooptimal 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 FeNi_{3} is not only contained in both sets, but it is also located at the “knee” of the activity paretooptimal set, which indicates that neighboring solutions are worse with respect to either activity or cost. Clearly, FeNi_{3} is a very promising catalyst candidate for the methanation reaction. This catalyst identification step concludes the theoretical design process and experimental verification of the theoretical prediction is necessary. Experimentally obtained methanation rates of Fe/Ni alloys as a function of Ni content are displayed in Figure 1.1(d) and clearly show that the computationally predicted FeNi_{3} 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 }
In the remainder of this chapter, the background information that leads to the identification of appropriate catalyst descriptors (e.g. dband 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 stepbystep 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.
1.2 Starting from the Electronic Structure
1.2.1 Density Functional Theory
Computational catalyst screening would not be possible without the existence of a theory that enables 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:
The ground state properties of a manyelectron system are uniquely determined by the electron density.
The total energy of a system has a minimum for the ground state electron density.
DFT provides a solution to the Schrödinger equation:
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 noninteracting electrons and a correction term, E_{XC}, which accounts for exchange and correlation effects in the interacting system. Although approximations for the description of the exchangecorrelation 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 manybody 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
The first term, T_{KS}[n(r)], is the kinetic energy of fictitious, noninteracting electrons and is obtained from the singleelectron Kohn–Sham equations
where v_{eff} 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, E_{XC}, in the total energy equation depends on the unknown exchangecorrelation 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–BurkeErnzerhof (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 generalpurpose 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 selfinteraction errors can be encountered. Several improvements to the GGA have been suggested (e.g. DFT+U, DFTD, metaGGA, hybridGGA), 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 nonDFT 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 tutorialstyle book Density Functional Theory – A Practical Introduction by David S. Sholl and Janice A. Steckel is a highly recommended starting point.^{33 }
1.2.2 The dBand Model
Computational catalyst screening relies on the prediction of correct trends across different catalysts rather than the prediction of quantitative rates and selectivities for each catalyst. Understanding the origin of the observed trends in terms of the underlying electronic structure can therefore be very helpful during the screening process. For transition metal surfaces, trends in reactivity can be very well described and understood in terms of the dband model (Figure 1.2) developed by Hammer and Nørskov.^{34,35 }
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 antibonding 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 lowerlying 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 sband (turquoise) and a narrow band of localized delectrons (red). Assuming that we can separate the coupling between the adsorbate level and the s and dbands of the transition metal, we can write the chemisorption energy as the sum of both interactions.
The interaction of the adsorbate state with the broad sband, ΔE_{S}, leads to a broadening of the state and a downshift in energy. Since all transition metals have broad and halffilled sbands, the contribution ΔE_{S}, is approximately the same for all transition metals. In the next step, the broadened molecular state can couple to the narrow dband of the transition metal, which causes a split into bonding and antibonding 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 antibonding states are located below the Fermi level, the weaker the resulting chemisorption bond. For stronger bonds it is desirable that all the antibonding states are high in energy and above the Fermi level.
The location of the antibonding states relative to the Fermi level and, in turn, the strength of the resulting bond, is largely determined by the position of the dband. Since the dband spans a range of energies, we need to find a convenient way to define the term “position”, and therefore introduce the dband center. The dband center, , is the energyweighted average of the density of dstates ρ. It sounds complicated, but it is the same equation that you use to calculate the center of mass (just substitute r_{COM}= , r=E and m=ρ).
The dband center varies across the transition metals and according to the dband model this will cause a variation in the interaction strength, ΔE_{d}. If we recall that ΔE_{s} is constant, the variations in chemisorption energy across transition metals can be attributed to changes in ΔE_{d} and, thus the dband center. Transition metals with higherlying dbands have stronger chemisorption properties. Numerous examples supporting the dband 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 dbands (d^{9} and d^{10} metals) there can be a significant repulsive interaction between the transition metal dband and the adsorbate states, which leads to a local reversal of the dband 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.
1.3 Identifying the Right Descriptor Set
If you have ever played the party game “Taboo” you have experienced how hard it can be to describe a 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 descriptorbased 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 dband 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 dband model and are discussed next.
1.3.1 Scaling Relations for Surface Intermediates
The binding energy of an adsorbed molecule is determined by the number and strength of the chemical bonds that it forms with the surface. To a first approximation 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.
1.3.1.1 Hydrogencontaining Molecules
The discussion in this section is fully based on the scaling papers by AbildPedersen 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, ΔE^{AHx}, of the hydrogencontaining intermediates AH_{x} (A=C, N, O, S) plotted as a function of the adsorption energy, ΔE^{A} of the central atom A for a range of typical transition metal surfaces. The adsorption energies were obtained from periodic DFT calculations on closepacked terraces (black), step sites (red) and, additionally, on the facecentered cubic [FCC(100)] surface (blue) for OH_{x}. It can be seen that the adsorption energies ΔE^{AHx} are linearly correlated with ΔE^{A} and given by
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, ΔE^{OH}, which scales with a slope of ∼1/2 with respect to the adsorption energy of oxygen, ΔE^{O}, on all three different surfaces. The same generalization applies to all the other adsorbates and can easily be rationalized using 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 CH_{x} species.
In general, the scaling slope γ in eqn (1.6) for hydrogencontaining adsorbates AH_{x} can then be predicted by
where x_{max} is the maximum number of H atoms that can bind to the central atom A. In particular, x_{max}=4 for C, x_{max}=3 for N and x_{max}=2 for O and S. It should be noted that the linear scaling behavior and characteristic slope of hydrogencontaining 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 ΔE^{AHx} and ΔE^{A} on a reference catalyst and we know x_{max} to predict the slope, γ. Then we can use eqn (1.6) to predict ΔE^{AHx} as a function of ΔE^{A} by simply eliminating the intercept ξ.
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 }
1.3.1.2 Extension to More Complex Surface Intermediates
The special case of hydrogencontaining AH_{x} species impressively demonstrates the relationship between the binding energies of similar adsorbates, but in the majority of catalytic reactions more complex surface intermediates are encountered. In the process of establishing linear scaling relations for complex species, 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 AH_{x} and apply it to CH_{x}–NH_{y} with x,y=0–2 surface intermediates. These intermediates occur, for example, during the Degussa and Andrussow process for HCN synthesis from CH_{4} and NH_{3}.^{17,18,40,41 } We may assume that the C and N atoms are connected through a single bond and that the CH_{x}–NH_{y} intermediates bind to the surface equally through the C atom and the N atom. The adsorption energy of CH_{x}–NH_{y} should then scale linearly as in
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 ΔE^{CHx} and ΔE^{NHy} into eqn (1.9) and collect the intercepts, we arrive at our final expression for ΔE^{(CHx−NHy)}.
The comparison of CH_{x}–NH_{y} 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 C_{2} hydrocarbons on transition metal surfaces.^{42 }
Even without deriving scaling relations on paper as done for CH_{x}–NH_{y}, one can simply try different scaling relations and evaluate them on the basis of their mean absolute error and/or R^{2} value. Figure 1.5(b) shows a very good scaling relation for C_{2}H_{2} and C_{2}H_{4} with methyl (CH_{3}), 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.
1.3.2 Scaling Relations for Transition States: The Brønsted–Evans–Polanyi Relationship
Transition states are often considered special because they have a very short lifetime. In terms of their adsorption properties, however, they are not different from any other surface species and obey the same physical principles. It is therefore not surprising that scaling laws for transition states exist in the same way as they do for other surface intermediates. The truly remarkable aspect is that transition state scaling for a large number of reactions is 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 E_{TS} and the energy of either the initial state (E_{IS}) or final state (E_{FS}) of a reaction. The transition state energy, E_{TS}, is not equal to the activation barrier E_{a} of a reaction, but it can be used to calculate E_{a}=E_{TS}−E_{IS}. A BEP relation, on the other hand, is a relation between the activation energy E_{a} and the heat of a reaction ΔE, which is defined as ΔE=E_{FS}−E_{IS}. All energy quantities are indicated on a schematic potential energy surface (PES) in Figure 1.6.
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, E_{diss} was used to denote the final state energy E_{FS}; and ΔE_{diss} represents the energy change of the surface reaction, ΔE. Although identical raw data have been used to generate the plots, the transition state scaling in Figure 1.7(a) appears to be much better than the BEP relation in Figure 1.7(b). This initial impression is somewhat deceptive and a comparison of the mean absolute error (MAE) 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 yaxes of the plots, and the BEP graph can be loosely interpreted as a zoomedin version of the transition state scaling graph.
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 E_{1} and E_{2}.
Using transition state scaling we can find E_{TS} as a function of E_{FS.}
The quantity that is of practical interest for us is not E_{TS} directly, but the activation energy barrier E_{a}, which is obtained in TSS from
The last three terms are constants and may be summed and replaced with β^{TSS}, which then yields the final expression
Alternatively, E_{a} can be obtained from the BEP relation.
Although eqns (1.15) and (1.16) look similar, they are not identical if the descriptors E_{1} and E_{2} 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
into eqns (1.15) and (1.16), we can eliminate E_{2} and get
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 E_{1} and E_{2} 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 statelike transition state. The energy of the transition state, E_{TS}, would then scale with E_{IS} and E_{FS} weighted by the transition state character ω. With β as an arbitrary energy offset, the energy E_{TS} can then be expressed as
Rewriting eqn (1.20) for the calculation of an activation barrier leads to a BEP relation relating E_{a} with ΔE and the slope of the BEP relation, α, can be interpreted as the transition state character, ω.
In addition to the physical interpretation of the BEP slope, using a BEP relation has other practical advantages. For the calculation of rate and equilibrium constants E_{a} and ΔE are needed, which can be directly obtained from a BEP relation. In more advanced kinetic models, the effect of surface coverage on the stability of the surface intermediates may also be explicitly included, but including such an explicit coverage dependency on activation energies is not trivial.^{48,49 } Assume that we have used DFT to calculate E_{a} for a surface reaction at low coverage and we are able to predict the initial and final state energies for a highcoverage scenario. Performing additional DFT calculation for E_{a} at various higher surface coverages is prohibitively time consuming and an approximate method is desired to estimate the effect of coverage on E_{a}. In this case, 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 E_{a}.
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.
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 diatomic dissociation reactions (N_{2}, NO, O_{2}) on metalsubstituted Laperovskites 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 ΔE_{diss}) and can be explained by a continuous change from a late transition state (final state like) to an early transition state on more reactive surfaces. The same 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.
1.4 The Sabatier Principle and the Volcano Curve
For most things in life, there is an optimal concentration, amount or frequency in which they should occur; catalysis is no exception. The Nobel Prizewinning 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 volcanoshaped 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 fullblown 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.
1.4.1 Sabatier Analysis
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:
Write the corresponding rate expressions for all reaction steps in the forward direction as a function of the approach to equilibrium , where K_{i} is the equilibrium constant for step i, a_{n} is the activity of species n and v_{n} is the corresponding stoichiometric coefficient.
For each reaction step in the model, assume optimal surface coverages for the forward reaction.
From the overall approach to equilibrium, calculate the approach to equilibrium of step i under the assumption that all other steps are quasiequilibrated.
Use TSS/BEP relations to calculate the forward rate constants k_{i} and the rate. This rate represents an upper bound on the actual rate of step i and will be denoted as .
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
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
The allowed range for γ is 0≤γ≤1, where γ=1 indicates a quasiequilibrated reaction. From eqn (1.25) we can estimate γ_{1} and γ_{2} under the assumption that the other step is quasiequilibrated and obtain
Assuming that we know k_{1} and k_{2} we can now write down the upper bounds for r_{1} and r_{2} and find the Sabatier rate r_{Sabatier}.
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.
1.5 Sabatier Analysis in Practice
1.5.1 First Example: Ammonia Synthesis
Now it is finally time to combine all the concepts introduced so far and perform 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 (NH_{3}) synthesis from nitrogen (N_{2}) and hydrogen (H_{2}).
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 stateoftheart surface science^{55 } 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.
Upon inspection of the PES for NH_{3} synthesis one can simplify the mechanism to two main reaction steps, the dissociative adsorption of N_{2} followed by its hydrogenation to NH_{3}. 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 NH_{x} species, but it is sufficient to illustrate the Sabatier principle. A catalyst that associates strongly with N will easily dissociate N_{2}, but the desorption step will be very difficult (note that in the desorption step only the energy of N* is affected by the catalyst because H_{2} and NH_{3} are gas phase species). In contrast, a catalyst interacting weakly with N allows for easy desorption, but the barrier for N_{2} dissociation will be large. Hence, there is an optimal binding energy of N*, E_{N}, at which the dissociative adsorption of N_{2} and NH_{3} desorption are optimized simultaneously. We will now try to find an estimate of E_{N}, by performing a Sabatier analysis for typical NH_{3} synthesis conditions.
The simplified NH_{3} synthesis mechanism is almost identical to the generic mechanism used in the introduction of the Sabatier analysis in Section 1.4. If we substitute A_{2}=N_{2} and B=3/2H_{2}, 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 H_{2} is included as the exponent to in the expression for .
For the reaction conditions we may assume that we have a stoichiometric feed composition (H_{2}:N_{2}=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 P_{N2}=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 k_{1} and k_{2}. We can obtain the rate constants, k_{i}, from the standard Arrhenius expression, which depends on the activation energy, E_{a,i}, and the preexponential factor, ν_{i}, for each step.
The preexponential factor ν_{i} can be calculated using transition state theory from the entropy change between the initial state and the transition state of the reaction, .
Next, we obtain the activation barrier estimates from published scaling relations, beginning with N_{2} 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 }
The second step in the simplified mechanism lumps several hydrogenation steps into one step, therefore it seems appropriate to use the averaged BEP parameters for NH_{3}, NH_{2}, 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
and we find the forward barrier by using the fact that . The energy change of step (2), ΔE_{2}, depends linearly on ΔE_{1} owing to the thermodynamic constraint that ½ΔE_{1}+ΔE_{2}=ΔE_{r}, where ΔE_{r} is the overall energy change of the reaction. So far we have exclusively used the electronic groundstate 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 DFTcalculated 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 ΔE_{r} ≈ ΔH°=−45.9 kJ mol^{−1}=–0.48 eV. Putting this all together we get the desired BEP relation for E_{a,2}.
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 ΔE_{1}, which is the dissociative chemisorption energy of nitrogen, ΔE_{N}.
The last piece of missing information comprises the entropy changes between initial and transition states, which are needed to calculate the preexponential 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 gasphase entropies of N_{2} and H_{2} at standard state are S°(N_{2})=192.77 J mol^{−1} K^{−1} and S°(H_{2}) 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 N_{2} loses all its gasphase entropy upon entering the transition state and in the second step, the entropy loss of one H_{2} gasphase molecule is considered when it reacts with N* on the surface. Why did we consider the entropy loss of only one H_{2} molecule and not 3/2H_{2} 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 H_{2} molecule adsorbing onto the surface –130.68 J mol^{−1} K^{−1}of entropy are lost. When mechanisms with lumped reaction steps are analyzed, problems with estimating energy barriers and entropy changes will frequently occur, so it is highly recommended to formulate mechanisms that are comprised of only elementary reaction steps! Still, for this case, from eqn (1.33) and the entropy assumptions outlined above, the preexponential factors are calculated as ν_{1}=1.20×10^{3} s^{−1} and ν_{2}=2.10×10^{6} 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 N_{2} dissociation rate and maximum desorption rate as a function of ΔE_{1}=ΔE_{N} and plot the results. The graph should look like the one given in Figure 1.11. The optimal catalyst for ammonia synthesis is found for a descriptor value of ΔE_{N}≈−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, ΔE_{N}, 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 ΔE_{N} spans a range from ca. –4 to +6 eV and neighboring transition metals are often separated by 0.5 eV or more.
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 
W  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 
W  Re  Os  Ir  Pt  Au 
−4.33  –  –  −0.59  1.37  5.89 
In summary, this example of computational screening for an NH_{3} synthesis catalyst has demonstrated how BEP relations are applied to reduce the number of parameters to a single descriptor, ΔE_{N}, and how to perform a Sabatier analysis. As a result 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.
1.5.2 Second Example: CO Oxidation
The Sabatier analysis for NH_{3} 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, P_{O2}=0.33 bar, P_{CO}=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 O_{2} is fast and can be considered as quasiequilibrated. The two competing steps in the reaction are then the dissociation of O_{2}* and the reaction between CO* and O* to form CO_{2}. Assuming very low conversions (γ→0) the maximum rates of steps (3) and (4) are:
The assumed optimal coverages that maximize are , and is maximized for θ_{CO}=θ_{O}=0.5, leading to our final expression of the Sabatier rate.
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 preexponential factors are typically on the order of ν=10^{13} 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 }
The TSS equations can be rewritten to yield activation energy barriers.
The Sabatier rate for CO oxidation in eqn (1.39) can now be obtained from the knowledge of ΔE_{CO} and ΔE_{O}, 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 twodimensional volcano for CO oxidation that depends on both descriptors, but studying a onedimensional volcano slice as shown in Figure 1.13 is easier to follow. This volcano slice was prepared under the artificial assumption that all transition metals have a constant binding energy of CO, which was set to ΔE_{CO}=–1.20 eV. This value is within ±0.3 eV of the actual CO binding energies 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 O_{2} dissociation, step (3). The optimal catalyst is the best compromise between all possible ratedetermining steps, such that, at the maximum, no single step can be identified as ratedetermining. It should also be mentioned that ratelimiting 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 ratelimiting 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 NH_{3} synthesis in the earlier example, the success of a Sabatier analysis can no longer be considered a coincidence!
1.6 Notes on Microkinetic Modeling
Although the Sabatier analysis can be applied to reaction mechanisms with any number of elementary steps and even mechanisms with parallel reaction pathways, it is an approximation. More accurate computational screening for complex reaction mechanisms can be performed by implementing a full microkinetic model, which contains information about the dominant reaction mechanism and simultaneously predicts the activity, selectivity and surface coverages. As input to this microkinetic model only the catalytic descriptors are needed and all the missing energy information must be estimated from BEP and scaling relations or assumed. But before 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 fourstep mechanism used in the Sabatier example will be used here and the net rates of the elementary steps are given by:
The change in coverage of reaction intermediates is determined by the net rate at which they are produced or consumed in each reaction step.
And the fraction of empty sites on the surface can at any point of time be calculated from the site balance.
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 steadystate conditions, which means that the derivatives of the surface coverages with respect to time in eqns (1.49)–(1.51) are all zero.
This converts the system of ODEs into a system of nonlinear algebraic equations that can be solved with standard rootfinding methods. The steadystate assumption is not strictly necessary when solving a microkinetic model numerically, but the transient startup and shutdown behavior is typically short in comparison to steadystate 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 quasiequilibrated. We can make this assumption for steps (1) and (2) of our mechanism and then find the surface coverages of CO* and O_{2}* as a function of the fraction of empty sites.
There is not a reaction that we can assume to be quasiequilibrated to derive an equation for the coverage of O*, but from the steadystate assumption and eqn (1.51) we find
which can be solved for θ_{O} by substituting and θ_{CO} from (1.54) and (1.55).
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 θ_{*}.
Equation (1.58) can be solved for any reaction condition if the equilibrium constants, K_{i}, for steps (1) and (2), and the forward and reverse rate constants, and , for steps (3) and (4) are known. Note that we do not need to know the rate constants of the quasiequilibrated 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.
1.6.1 Numerical Solution Strategies
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 rootfinding 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. Rootfinding methods can quickly solve the steadystate equations, but the rate equations in microkinetic models are highly nonlinear 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.
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 steadystate solution the ODEs need to be integrated over a sufficient amount of time to ensure that steadystate has been reached. It is a good practice to verify steadystate 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 rootfinding strategy also fails, it is possible to reduce the stiffness of the ODE system by artificially slowing down some of the very fast, quasiequilibrated, 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.
1.6.2 Entropy and Enthalpy Corrections
In the previous examples we only considered electronic energy changes and approximated the entropy as all or nothing. In essence, we assumed that gasphase 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 gasphase molecules, adsorbates cannot freely rotate and move across the surface, but the translational and rotational DOF are frustrated within the potential energy well imposed by the surface. In the harmonic limit the frustrated translational and rotational DOF can conveniently be described as vibrational modes, which in turn means that any surface adsorbate will have 3N vibrational DOFs that are all treated equally.
The temperature dependent surface enthalpy, H(T), of an adsorbate is composed of the electronic energy, E_{elec}, the zeropoint energy, E_{ZPE}, and a heat capacity correction to account for the difference between the enthalpy at 0 K and at the actual temperature, T.
The electronic energy component, E_{elec}, is directly obtained from either DFT calculations or from scaling relations based on DFT energies. The remaining components can be calculated in the harmonic approximation.^{61 }
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 twodimensional surface gas, which can be derived from the partition function of a particle that can move freely across the surface.
In this equation m denotes the mass of the particle and C_{S} is the number of sites per unit area, which is roughly 10^{15} 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.
1.6.3 Microkinetic Model Analysis
Microkinetic models can provide much more information than “just” surface coverages and reaction rates of individual steps. It is always a good idea to perform a sensitivity analysis on the model input parameters, to obtain additional valuable information. First, more attention should to be paid to parameters to which the model shows a high sensitivity, and these parameters need to be carefully measured or calculated. Nonsensitive parameters, on the other hand, may be roughly estimated or even guessed. Second, the sensitivity of a microkinetic model towards certain parameters can be interpreted as a measure of the degree to which a step is rate controlling (degree of rate control). This information can provide direction for useful catalyst modifications which could optimize performance (degree of catalyst control).
1.6.3.1 Degree of Rate Control
The degree of rate control, X_{RC}, was originally proposed by Campbell as a quantitative measure to identify ratecontrolling 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, k_{i}, while keeping the equilibrium constant, K_{i}, and the rate constants, k_{j}, of all other steps constant.
An X_{RC,i} value of zero indicates a quasiequilibrated step, whose rate constant has no effect on the overall rate, while X_{RC,i}=1 indicates a single ratecontrolling 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 finitedifference equation that can be readily evaluated in a microkinetic model.
In practice X_{RC,i} is then calculated by changing the forward and reverse rate constants of step i simultaneously by a small amount, typically 10% (x=0.1), and monitoring the change of the overall reaction rate with respect to the reference rate, r^{0}, for the original values of k_{i}. Never forget to apply this change to both rate constants, otherwise K_{i} and the overall equilibrium constant, K_{eq}, will be altered as well!
Fundamentally, the definition of the degree of rate control in 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 CH_{4} 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
where c is a constant ( in transition state theory) and G_{a,i} is the free energy of activation for step i. Upon substituting eqn (1.64) into the definition of X_{RC,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, (G_{TS,i}).
Equivalently, one can change the free energy G_{n} of an intermediate n and calculate its degree of rate control. This is schematically shown for the free energy of C* in Figure 1.15 and one can define the thermodynamic degree of rate control, X_{TRC,n}, for an intermediate n as
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.
1.6.3.2 Degree of Catalyst Control
Campbell's degree of rate control is an extraordinarily useful concept for the analysis of reaction mechanisms, the identification of ratelimiting 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 dband model. It is therefore unrealistic to assume that we can find catalysts that, for example, reduce the activation energy barrier of a ratelimiting step without simultaneously affecting other parts of the PES. This leads to the idea of the degree of catalyst control X_{CC}, which is a constrained sensitivity analysis of a microkinetic model with respect to a reactivity descriptor E_{i}.^{65 } The definition of X_{CC,i} is analogous to the generalized degree of rate control and may be expressed as
In eqn (1.67) it is indicated that all other descriptor values E_{j≠i} and the implemented scaling and BEP relations are kept constant when taking the partial derivative. The rationale for this definition is that the catalytic descriptors are the only independent variables and X_{CC} will indicate the direction one must look to identify an improved catalyst. The optimal catalyst is characterized by X_{CC,i}=0 for all descriptors E_{i}. X_{CC} is proportional to the slope of the volcano curve and following it in the uphill direction will get you to the top!
The difference between X_{RC} and X_{CC} is shown graphically in Figure 1.15 for the sensitivity towards the binding energy of carbon, ΔE_{C}. Carbon adsorbs quite strongly and is located in a well of the PES, which explains the negative value (–0.26) for X_{RC}(C). Decreasing the depth of the well would smooth out the PES and increase the reaction rate. However, decreasing the stability of C* will also affect the stability of CH*, CH_{2}*, CH_{3}* surface intermediates and the associated transition states between them through scaling and BEP relations (dashed line in Figure 1.15). In particular, the transition state energy of CH_{4} activation will shift up in energy. The large value of X_{RC}(CH_{4}TS)=0.8 indicates that CH_{4} activation is a ratedetermining step and increasing the barrier will have a large negative impact on the overall reaction rate. The overall effect of changing the descriptor ΔE_{C} on the PES is given by X_{cc}(C)=0.11, indicating that a catalyst with stronger ΔE_{C} will have higher activity. The explanation is that the increased poisoning effect of C* is more than compensated by lowering the barrier for the CH_{4} activation step.
The degrees of rate and catalyst control both have their merits and do not conflict with each other. Determining which to use is dependent on the problem at hand. X_{RC} should be used to analyze the reaction mechanism and understand the importance of individual steps and adsorbates, whereas X_{CC} should be used to guide catalyst design.
1.7 CO Oxidation Catalyst Screening
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 O_{2} are quasiequilibrated in Section 1.6. Now we will formulate a numerical solution to the complete microkinetic model as a function of the descriptors ΔE_{CO} and ΔE_{O}. 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.
1.7.1 Numerical Microkinetic Model
Microkinetic models are typically custommade 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 rootfinder, and ndimensional curve fitting routines, if parameter optimization is also required. Here we use python, which is an objectoriented 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 selfexplanatory 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.
1.7.1.1 Constants and Conditions
It is a good idea to define the reaction conditions and useful physical constants globally at the beginning of the code, so that they are available to all functions. This is also the place to import additional modules, such as numpy.
1.7.1.2 From Descriptors to Rate Constants
The next important step is to implement the scaling and BEP relations together with any assumptions regarding entropies of transition states and surface states in order to obtain the necessary rate constants. Since we will need this procedure multiple times for each combination of our descriptors, we define it as a function.
In this function we have used the scaling/BEP relations, eqns (1.42)–(1.44), and we have arbitrarily assumed that CO and O_{2} lose 25% of their gasphase entropy when entering the transition state for adsorption. Using collision theory for the adsorption rate constants^{66 } of these two steps is a more elegant alternative and the reader is encouraged to see how it affects the result of the model. There are two more noteworthy points. First, the model ensures that the activation barriers, E_{a,i}, are all positive, even for descriptor pairs that would otherwise predict negative values for E_{a,i}. Second, thermodynamic consistency is enforced throughout the model by calculating the reverse rate constants from the equilibrium constants for each step.
1.7.1.3 Rate Equations
This step is simple. 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.
1.7.1.4 System of ODEs vs. SteadyState Equations
As discussed in Section 1.6.1 the microkinetic model may be solved as a system of ODEs or nonlinear algebraic equations using the steadystate 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 rootfinding 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().
1.7.1.5 Solving the Model
All that is left at this point is to set the values of the descriptors and solve the model. For single point calculations using only one descriptor set, it is preferable to use an ODE solver and use a completely empty surface as the initial guess for the coverages. The integration time required to reach steadystate depends on the problem and can vary from 1×10^{3} to 1×10^{8} seconds. Some might wonder why a reaction could possibly take 1×10^{8} 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.
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 steadystate has been reached.
As a result of this exercise, we find the steadystate coverages on Pt of θ_{CO}=0.215 ML, θ_{O2}=8.4×10^{−4}, and θ_{O}=5.6×10^{−3} and a steadystate CO_{2} formation rate of 6.16×10^{3} s^{−1}.
1.7.2 Degree of Rate and Catalyst Control
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.
The degree of rate control, X_{RC}, is then calculated by repeatedly calling solve_ode with rate constants that are systematically changed for each step. As input to calculate_Xrc we provide the unaltered rate constants and the corresponding rate.
In our example for CO oxidation on Pt with ΔE_{O}=–1.25 eV and ΔE_{CO}=−1.22 eV the degrees of rate control for steps (1) and (2) are zero (these steps are quasiequilibrated), step (3) is almost exclusively ratecontrolling with X_{RC,3}=0.99, and step (4) has an X_{RC,4}=0.01. In this example , which is expected for a serial reaction mechanism.
Next we calculate the degree of catalyst control X_{cc} by varying the values of ΔE_{O} and ΔE_{CO} by 0.05 eV. This step size should generally work well except in close proximity to the volcano peak, where smaller step sizes are preferable. Rather than providing the rate constants directly, we pass the original values of ΔE_{O} and ΔE_{CO} to the function calculate_Xcc and the corresponding rate constants are calculated on the fly.
If we perform this analysis for CO oxidation on Pt, we get X_{cc}(ΔE_{O})=1.38 and X_{cc}(ΔE_{CO})=–0.28, which indicates that by stabilizing O* and destabilizing CO* the reaction rate can be increased. Indeed, on the twodimensional volcano plot in Figure 1.16 it can clearly be seen that the maximum is located to the left of Pt in the direction of stronger O* binding. The resolution of the volcano plot is however not sufficient to judge whether stronger or weaker CO* binding would increase the activity. It appears as though ΔE_{CO} for Pt is already near to its optimal value.
1.7.3 Twodimensional CO Oxidation Volcano
After this thorough analysis of the microkinetic model results around a single point (Pt), the only remaining task is to calculate the CO oxidation rate over a larger descriptor range and visualize the results as a twodimensional 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 O_{2} adsorption steps are quasiequilibrated 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 rootfinding 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 rootfinding method. An example of how this may be implemented is provided below.
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 twodimensional volcano graph can be found in the Appendix.
1.7.4 Effect of Lateral Interactions
Throughout this chapter it has been assumed that binding energies and activation energy barriers are constant, but at higher surface coverages these quantities can be strongly influenced by adsorbate–adsorbate interactions. An accurate description of these interactions, specifically when they induce local ordering and lead to island formation or other characteristic patterns, would require the use of kinetic Monte Carlo (kMC) simulations. Since computational catalyst screening is largely based on the use of rather approximate scaling and BEP relations with mean absolute errors of ca. 0.2 eV, the use of highly sophisticated kMC methods does not seem appropriate. Adsorbate–adsorbate interactions can also be described within the meanfield 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 dband 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 coveragedependent activation barriers for NO oxidation, which were later coupled with first principlesbased 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 meanfield 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 coveragedependent microkinetic model for CO oxidation can be constructed and used for screening purposes. Because the lateral interactions vary across catalysts and no correlation between the interactions and the activity descriptors ΔE_{O} and ΔE_{CO} could be identified, the volcano in Figure 1.17(a) was calculated 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 catalystindependent 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 coveragecorrected position. While the long answer is fundamentally interesting, the important takehome message is that lateral interactions neither alter the position of the volcano peak nor change the relative order of the activity of metal catalysts.
1.8 Conclusions
This chapter has provided a tutorial style overview of the rapidly emerging field of computational catalyst screening and the reader should feel adequately informed and prepared to attempt a catalyst screening project on 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 descriptorbased 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 indepth 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 (= longterm 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 nonexpert what you do for a living, try starting with: “I am a professional matchmaker and my clients are slow chemical reactions…”!
Appendix
Python Code for the numerical solution of the microkinetic model for CO oxidation and the generation of a twodimensional volcano plot.