Perspective: On the active site model in computational catalyst screening

Over the last decade, a new theme has made a staggering appearance in theoretical surface catalysis research: computational screening as a means of catalyst discovery.^1–6 The basis for such screening approaches is the realization that there are many dependencies in the adsorbate and transition state (TS) energetics, the so-called scaling relations, that largely govern the microkinetics of a given catalyst material. Brønsted-Evans-Polanyi relationships link activation barriers to thermochemical reaction energies, while thermochemical scaling relations correlate the underlying adsorption strengths of reaction intermediates to those of their constituting base elements.^6–12 Together with (sometimes drastic) assumptions on the reaction mechanism, these dependencies often allow for the description of the catalytic activity in terms of only a few simple parameters. Common examples of such parameters, generally called descriptors, are the adsorption energies of one or more reaction intermediates or constituting base elements like C, N, or O. Since these adsorption energies also follow simple trends (at least over the transition metal series^13,14), one can either deduce clear guiding principles for catalyst design or screen a large number of catalyst materials on the basis of only a limited number of first-principles calculations.

One can hardly overstate the impact that this type of work has made in recent years. Computational screening is now a widely accepted (sometimes already considered essential) strategy to guide experimental catalyst synthesis. Reports on the corresponding identification of improved catalyst materials are already piling up.^15–21 Maybe even more important though, is the conceptual framework that has come with this approach.^4,6,14,22,23 Thinking in terms of volcano plots and understanding the underlying trends that give rise to them; all this has become natural to us and thereby encompasses much, if not most, of the knowledge that had been obtained in a century of empirical catalyst research. An understanding of these trends has furthermore clarified the fundamental limitations to the maximal catalytic activity achievable for a given class of catalysts. It has provided leads on how to overcome such limitations employing, e.g., promotion, multidimensional, or multifunctional binding sites and support effects.^5,24,25

That said, it is however also clear that the approach necessarily has limitations—even though this is less emphasized in the wake of its current of success. Foremost, this concerns the reaction mechanisms utilized in these studies. Often, a single pre-determined mechanism is assumed, even when screening systems of largely varying reactivity. On the energetic level, one might wonder about inaccuracies of the employed scaling relations which sometimes are more of a loosely scattered trend rather than an exact relation. Also, it is still an open question whether scaling relations can accurately describe more complex systems containing adsorbates of extended molecular structure, high surface coverage, or long-range interactions.^6,25 And finally, one might wonder whether the typically employed low-rung density-functional theory (DFT) functionals can really compute the descriptors used in these studies adequately.²⁶ 

Addressing these issues is currently the subject of much research—be that the use of improved DFT functionals,^27–30 sensitivity analyses,^31–35 or the consideration of extended reaction networks and multiple descriptors. Even if we assume that the mechanism indeed proceeds the way we believe it does and even if we had the exact DFT functional, there is however still one crucial input which computational screening approaches in surface catalysis heavily rely on: the chosen model for the active site(s). This aspect is less frequently discussed but, in our view, possibly even more critical than any of the other issues. When writing down generic reaction mechanisms, there is little specification of the active site to which the reaction intermediate(s) bind. It is merely denoted with a “B.” Unfortunately, this does not work for DFT (or any other molecular-level computational approach required to compute the descriptors). For such calculations, we need to place atoms into the structure file. In other words, we need to know much more about the active site(s) driving the reactions. Depending on the electronic screening length of the material, we need to know at least the first shell of the neighbors around the actual substrate atoms to which the reaction intermediates coordinate. Often, this extends much further to appropriately capture the long-range electrostatics, dielectric response, or dispersive interactions involving the catalyst material. Local geometry optimization allows for the relaxation of the atomic positions within a chosen local configuration for the active site; however, this will generally not modify the configuration itself. This configuration—or atomic-scale structure—of the active site has to be chosen by us. In almost all computational studies to date (certainly all computational screening work) this is thus an input, not a result of the calculations.

So why bother? In general, we can expect the binding of intermediates to decisively depend on the chosen model for the active site.³⁶ The choice can thus critically determine the outcome of a screening study. This would not be a problem if it was straightforward to identify the active site(s). By discussing a series of catalytic systems, the purpose of this perspective is to highlight that this is generally not the case. In fact, for many (if not most) working catalysts, our current atomic-scale characterization of the active site(s) is indeed not much deeper than the infamous “B.” We illustrate this with a few instructive examples covering transition metal (TM) nanoparticles and compound catalysts. In all cases, our present understanding of the catalytic function suffers from an incomplete picture of the active site. On the experimental side, there are dedicated efforts to sharpen this picture, for instance, through in situ or operando studies.^37–40 On the computational side, this is mirrored by detailed mechanistic studies.^41–48 

Nevertheless, for the time being, we will have to accept that there is, in general, considerable uncertainty in which active site model to employ in computational screening studies. Recognizing this is already a first important step. It dictates great care and caution when assessing the results of present-day studies, particularly for complex compound materials. As emphasized at the end of this perspective, actively embracing this uncertainty also offers intriguing prospects. We see the chance to extend computational catalyst screening beyond materials discovery towards active site engineering.

The work on scaling relations and other related theories started out for TM catalysts and “simple” reactions like the Haber-Bosch process for ammonia synthesis.^49,50 From both experimental and theoretical studies,^51,52 it is well-established that the activity for dissociative adsorption of N₂, which is believed to be the rate-limiting step in the commercially employed iron-based catalyst, is highly dependent on the geometry of the sites exposed by the metal facet. Early experimental investigations on single-crystal Fe model catalysts suggested that the bcc(111) facet contains the active site for N₂ dissociation.⁵¹ These findings were later backed up by DFT calculations identifying a favorable dissociation pathway at the so-called C7 site highlighted in Fig. 1(a). This site is characterized by a fourfold motif of under-coordinated, and therefore reactive, metal atoms.⁵³ It is thus a site that intuitively appears suitable to break the strong N₂ bond. A corresponding site allowing for high coordination of the transition state to under-coordinated metal atoms cannot be found on the low-index surfaces of the more common fcc or hcp TMs. This naturally leads the thinking towards defects such as steps on these other TMs. Indeed, (100)-type steps on hcp(0001) or fcc(111) surfaces feature a similar fourfold motif involving (at least at the upper step edge) heavily under-coordinated metal atoms, cf. Fig. 1(b). Consistently, this step site was shown to be an active site for N₂ and NO dissociation at Ru(0001) catalysts.^52,54,55

With this understanding, using (100) step sites as a representative active site motif for N₂ dissociation in ammonia synthesis appears rather compelling. The activity trends over the TM series obtained on the basis of this choice are furthermore fully consistent with our empirical knowledge and identify Fe and Ru as close to the top of the volcano.^33,49,56 So where is the catch? As apparent from this example, making the design choice for the active site model to be employed in a screening study crucially relies on detailed understanding as obtained in dedicated mechanistic studies on particular systems. When it comes to atomic structure details defining the active site, these are typical studies done on single-crystal model catalysts and often in ultra-high vacuum (UHV). This understanding can seem rather complete as in the N₂ example. More generally though, it will be rather shallow and then leaves quite some freedom in selecting what active site model we use as the basis of a screening study. In either case, the mechanistic studies are furthermore predominantly carried out for (model) catalyst materials that are empirically known to perform well for the targeted reaction. Another question is thus how far the mechanistic insight and the extracted active site models are still valid for different materials.

On the one hand, one can generally argue that we are mostly interested in the top of the activity volcano. Starting from a well-performing material that we can expect to already be rather close to the top, the screening might not have to extend to (or be reliable for) materials that are chemically very different. On the other hand, the larger the material space one can cover with a screening study, the more likely the possibility to identify a truly novel material. In addition, even within the comparatively small material space spanned by the elemental TMs, the reactivity already changes dramatically. While for less reactive metals, a maximally coordinating step site might be a good model for an efficient dissociation site, lower coordinating terrace sites may be more suitable for more reactive metals. Such variations in the dominant active site over the TM series are generally not accounted for in the bulk of present-day computational screening studies.

This might not be so problematic for the trends in a simple reaction like ammonia synthesis. However, the limitations become clearly visible for more complex reaction networks like the reaction of synthesis gas (CO and H₂), with its many possible end products such as methane, methanol, higher alcohols, or C₂₊ Fischer-Tropsch products. With the exception of methanol formation, which can be catalyzed by a more noble metal such as Cu,^57,58 the formation of the other products requires dissociation of the CO molecule and therefore relies on more reactive metals such as Fe, Ru, Co, Rh, Ni, and Pd. Even though these six named metals have been heavily investigated and span only three groups in the TM series, there is still no consensus on the active site for the CO dissociation step on them. Several experimental and theoretical studies have argued in favor of open step and kink sites.^59–62 In analogy to the ammonia case, (100) step site models have been extensively used in screening studies of CO methanation and the reverse reaction, methane steam reforming, over the entire TM series.^{17,34,63–66} Multiple mechanisms have been suggested, though, that facilitate CO dissociation, e.g., through H-assisted pathways, at high surface coverage,^60,67–70 along hydroxyl-assisted pathways⁷¹ or in the vicinity of growing hydrocarbon chains.⁷² Several studies focusing on Fischer-Tropsch synthesis have therefore argued that already for Ru and Co, the active site for dissociation is found on the close-packed hcp(0001) surface,^68,71,72 cf. Fig. 1(c), while the lower thermodynamic stability and high affinity for reaction intermediates might actually poison the kink and step sites under reaction conditions.

Even over a range of just three TM groups, our current level of characterization and understanding does thus not allow us to conclusively decide whether an open or a close-packed active site model is the better base for a screening study or whether the change in TM reactivity might not induce a transition from one dominant site type to the other. Maybe it is precisely the interplay between different site types on materials near such a transition that makes these materials particularly good catalysts. Only the most recent studies start to address such issues, for instance, by comparing the activity trends obtained for step and terrace active site models.^70,73–75

Even within the appealing structural simplicity of elemental TMs, the discussion does not end at the level of just terrace or step sites. First of all, the low-index surfaces feature a number of different step or kink sites. In general, our present mechanistic understanding is not good enough to really discriminate between them. Second, an active site model needs more definition than just the nearest-neighbor substrate coordination. Strain and charging are just two of many important support effects that have been identified for TM nanoparticles.^76–81 Maybe we also want to (or should) include such aspects into the active site model used for the screening. Going beyond gas-phase reactions to electrocatalytic reactions occurring in solution, the definition of the active site should possibly also be extended to account for solvent effects. As outlined in a recent review, catalytic activities of metal electrocatalysts have in many instances been found to depend strongly on the electrolyte composition, in particular, alkali metal cations, anions such as $SO 4 2 −$⁠, pH, and the use of ionic liquids.⁸² 

All of these uncertainties still exist within the fundamental assumption that it is the nominal TM material (and its nanoparticle facets) that carries out the catalysis. Notwithstanding, radical (surface) morphological transitions under reaction conditions, at least for the more reactive metals, have long been suspected.^40,83–85 For Fischer-Tropsch catalysis, this pertains to the formation of (surface) carbides,^86–89 for ammonia synthesis, the formation of (surface) nitrides,^39,90,91 and for oxidation reactions, the formation of (surface) oxides.^92–95 Obviously, it is unlikely that any active site model motivated by motifs at the surface of the elemental TMs will allow one to assess the reaction energetics at such formed compound materials.

The discussion in Sec. II underscored that the design choice of the active site model in a screening study already faces considerable uncertainty even for apparently “simple” TM catalysts. Yet, this is nothing compared to the situation we are confronted with when studying the much more structurally complex compound materials. This is particularly cumbersome since these materials are increasingly focused on in modern energy and catalysis research, and the vast material space spanned by them so much as calls for an exploration through computational screening. To make things worse, such materials like oxides, (oxy-)nitrides, or carbides already challenge us in many other aspects important for a first-principles based screening. This includes the necessity for more advanced DFT treatments (often at least at the GGA + U or hybrid functional level) to appropriately describe the bonding^27,96,97 as well as energy trends that are not necessarily as simple and inevitably established as at the d-band filling dictated TMs.^6,98

In comparison to TMs, we also have to live with a much less elaborate atomic-scale understanding, even for the few fruit fly reactions carried out by compound catalysts which the large majority of mechanistic studies has focused on to date. One of these reactions, somewhat comparable to ammonia synthesis in its simplicity and popularity, is the water oxidation or oxygen evolution reaction (OER) at the prototypical photo(-electro)catalyst TiO₂.^99,100 Work pertaining to this reaction on single-crystals has largely focused on rutile-structured TiO₂ and its most stable (110) facet.¹⁰¹ This surface is among the structurally most simple oxide surfaces one can imagine. Its catalytic activity is generally ascribed to the under-coordinated surface metal atoms exposed by the stoichiometric termination,^{94,95,101–103} cf. Fig. 2. For the OER, detailed computational studies established the oxygen binding strength to this so-called coordinatively unsaturated (cus) site as a suitable descriptor for the catalytic activity.^104,105 The volcano-shaped activity trends obtained with this descriptor provide an intuitive picture of the overpotential-determining reaction steps, and similar to the ammonia synthesis example, the predicted overall activity ordering of rutile (and a number of other) oxides seems generally consistent with available experimental data.^104,106

This conclusive picture has recently been questioned by a recalculation of the binding of the underlying reaction intermediates on TiO₂(110) at the more accurate hybrid functional level.¹⁰⁷ At this level, the initial potential-determining proton-coupled electron transfer step towards adsorbed OHB already exhibits a thermodynamic barrier that straddles the TiO₂ band gap. Particularly when considering an additional kinetic barrier for this step,^{103,108–110} this suggests that rutile TiO₂ would not be able to photocatalytically split water—something that is rather difficult to reconcile with this material being the fruit fly among water splitting catalysts.^99,100 An obvious explanation for this discrepancy would be that the studied cus site on pristine TiO₂(110) is not the active site responsible for this photocatalytic function of TiO₂. If so, the question is how much we can trust the computational screening results if they are based on an active site model that is not even correct for one of the very few oxides for which we have gained a somewhat deeper characterization at all.

Recognizing their importance for oxide properties in general,¹¹² point vacancies appear as a straightforward alternative for the active site. In fact, studies of the latter on TiO₂(110) indeed show that they can act as attractors for localized charge carriers^113,114 which, in turn, could drive the OER.¹¹⁵ However, point vacancies are just one example of how not only the structure but also the composition can be changed at the surface of compound materials like oxides. Both structure and composition can deviate significantly from a mere stoichiometric truncation of the bulk lattice and are furthermore a sensitive function of the surrounding environment (reaction conditions). This is illustrated in Fig. 2, which shows the change of the most stable surface termination of rutile RuO₂(110) under increasingly O-rich conditions. Specifically, near-ambient O atmospheres up to elevated temperatures stabilize an O-rich termination in which the cus site is only available in the form of a point vacancy surrounded by terminal O atoms (or hydroxyl groups in humid environments).^111,116,117 This simple example demonstrates nicely that the surface termination (i.e., structure and composition) of the compound surface under reaction conditions does not necessarily have to correspond to the surface termination characterized, for instance, under UHV conditions. This is rather consequential, since most of the experimental atomic-scale characterization we have available to date was precisely obtained under UHV. Great care must therefore be exerted when using such information to establish an active site model for a screening study.

In this situation, ab initio thermodynamic studies of the type underlying Fig. 2 are increasingly used to identify the most stable surface terminations under applied gas-phase conditions in heterogeneous catalysis⁸⁵ or under applied potential and pH in electrocatalysis.^{104–106,118–121} Further combining this with Wulff constructions provides information about nanoparticle shapes and the dominant facets that are exposed.^122–126 Both the knowledge on the dominant facets and the surface termination may then be used for the generation of suitable active site models for the targeted reaction conditions. Unfortunately, this is also not without flaws. A most fundamental limitation derives from the employed methodology itself. The predictive power of present-day ab initio thermodynamic approaches is restricted to the pool of candidate structures explicitly tested. In other words, the method only yields the most stable structure among the ones tested. If the truly most stable structure, say a complex surface reconstruction, was not inside the pool, it will not be identified. Instead, another surface termination will erroneously be pre-dicted to be stable and may mislead us in the generation of our active site model.

Another limitation is that the predominant active site does not necessarily have to be located on the most stable facet, nor does it have to be the most stable, and thus most prevalent, site at this facet. This aspect was recently highlighted for the OER on cobalt oxides, which are of significant interest as earth-abundant anode materials for practical photocatalytic water splitting devices.¹²⁷ Computational studies in the literature have targeted the mixed oxide phase Co₃O₄,¹⁰⁶ as well as the more oxidized CoOOH phase.^120,121 A screening study of OER catalysts that included Co₃O₄¹⁰⁶ assumed that the active site is a single Co center on the (001) surface, which is the most stable surface termination under OER conditions. With the same thermodynamic stability argument, a site on the $( 101 4 ¯ )$ surface of CoOOH was examined.^120,121 Despite the fact that these sites are likely the most prevalent sites under experimental conditions, a recent experimental study by Frei and co-workers on Co₃O₄¹²⁸ suggested that less prevalent defect sites account for most of the observed activity. This was based on the fact that the activity measured for a single active site using a transient isotopic FTIR technique was 100 times higher than the average activity of all exposed Co centers.

This situation was further analyzed by a first-principles microkinetic study that explicitly examined the OER activity of a number of sites on different Co₃O₄ facets.^129,130 This comprised single (dual) sites involving one (two) redox-active Co centers on the most stable (001) and second most dominant (110) facet, as well as on a step site on the (001) facet as featured at the (311) vicinal. As summarized in Fig. 3, the single center site on the (001) surface that was focused on in previous studies exhibits an insignificant activity. All other sites largely surpass this activity, with turnover frequencies (TOFs) that are up to eight orders of magnitude larger. Intriguingly, the site giving the highest TOF also changes as the applied potential is increased. Sites containing more highly-coordinated O species are more active at lower applied potentials, while sites containing less highly-coordinated O species dominate at higher applied potentials. Close to the applied potential corresponding to the experimental work of Frei and co-workers,¹²⁸ the dual Co step edge site is several orders of magnitude more active than any other site, cf. Fig. 3. This is consistent with the experimental observation that the OER is predominantly driven by a minority defect site.

Both the results from the experiment and the first-principles microkinetic study support the conclusion that the thermodynamically most stable single Co center site at the most stable (001) facet plays only an insignificant role for the electrocatalytic function. This example thus nicely demonstrates the danger of basing a computational screening study on such an active site model as determined through ab initio thermodynamics. However, ab initio thermodynamic studies are presently among the few means we have at all to motivate a specific geometric structural motif for the active site under reaction conditions. Until new methodology allows us to overcome this situation, the conclusion is that we obviously will have to live with a considerable uncertainty with respect to which active site model we use for catalyst screening. This does not necessarily invalidate the trends and understanding we derive from the screening. Yet, the minimum we have to ask is how robust this understanding and the possibly favored catalyst materials are with respect to this uncertainty.

The quantification of trends and dependencies in the reaction energetics at solid surfaces has paved the way for powerful first-principles computational screening studies in heterogeneous catalysis. The discussion in this perspective underscores that the employed active site model is a most crucial design choice in such studies. The validity of the obtained activity trends and of the concomitantly favored catalyst materials stands and falls with the degree to which the used model is indeed representative of the true active site(s). In principle, this is self-evident. In practice, and dazzled by the undoubted accomplishments of the seminal screening studies, it might be good to take a step back and remind us of this simple truth.

Initial screening studies focused on transition metal catalysts. In the absence of reaction-induced morphological transitions, the structural complexity of this material class seems limited. Still, even if we only consider the most obvious active site types, there are at least the different terrace and step sites of all low-index facets to choose from. For structurally more complex compound materials, the number of conceivable site types explodes. At present, our lack of knowledge on which one to best base the screening is presumably the biggest roadblock to an extended and effective computational exploration of the corresponding vast material spaces.

Recognizing the considerable uncertainty about the nature of the active site for many (or most) catalytic systems, one may wonder why or to what extent screening studies work. Unfortunately, the question why a given screening study works or fails in a specific case is very difficult to answer. There could be at least three reasons why a study successfully identifies the most active catalyst: (i) it assumed the correct active site and mechanism, (ii) the results are insensitive to these assumptions, and (iii) the authors simply “got lucky” (e.g., fortuitous cancellation of errors). As to screening studies that fail to identify the most active catalyst, there are unfortunately not enough (if any) examples available in the literature to allow to decidedly judge when and why this might happen. This may be due to the fact that these failures are usually not published. Part of the focus of this perspective is to encourage detailed experimental and theoretical studies that allow us to narrow down the candidate active sites and mechanisms and thereby help to answer the above question. In this respect, detailed mechanistic studies focusing on one catalyst material are a valuable and necessary complementary companion to screening approaches.

Computational screening neither will, nor strives to, provide as accurate and comprehensive a description as a detailed mechanistic study. Particularly for complex compound materials, it is unlikely that any one employed active site model will really be a faithful representation of the predominant active sites over the entire range of material space investigated. A useful way to think about the choice-of-active-site problem might instead be to view the space of possible site models in the same way as the space of possible materials: as something to explore. Straightforwardly, one can thus perform the materials screening separately for a number of different sites and see how robust the obtained insights are with respect to the “optimum” material. Several very recent studies have started to go in this direction.^70,73,74 Conversely, any variations in the trend over the explored sites serve to identify how the catalyst performance depends on the different structural motifs.⁷⁵ Screening studies considering several sites can naturally encompass changes of the predominant active site, for instance, due to strong changes in the reactivity of the material or due to changes in the reaction conditions (temperature, applied potential, etc.).

Systematically extending the screening over many candidate sites, an intriguing perspective is to establish “material–site-type–activity” maps. Just as much as the present “material–activity” volcanos help to identify promising materials, such extended maps will identify geometric motifs that would further optimize the performance of a given material—and thereby guide synthesis endeavors aiming to stabilize such motifs. In light of this, an interesting direction for methodological work would be to employ global geometry optimization approaches to generate (diverse) pools of candidate sites. Accounting for a simultaneous presence of different active site types within extended microkinetic models would in turn even allow one to capture (and engineer) a possible synergistic interplay. In this respect, the evident uncertainty in the choice of active site today offers exciting prospects for new research directions in the future. Eventually, actively embracing the space of active site models in screening studies may thus even provide leads towards the targeted design of optimized sites in future catalysts. On route to such grand goals, we can in fact start tomorrow and with a very modest step—simply by ceasing the annoying tendency to shift any information on the employed active site model to the supporting information or to not include it in a screening publication at all.

M.A. would like to thank the Alexander von Humboldt Foundation for a post-doctoral fellowship. Funding through the Deutsche Forschungsgemeinschaft is acknowledged within Project No. RE1509/27-1.