State-of-the-art and challenges in theoretical simulations of heterogeneous catalysis at the microscopic level

Abstract

Theoretical simulations of systems that represent heterogeneous catalysts constitute one of the main tools in the research for new catalytic materials. Theory plays a role in the three stages of the development ladder: characterisation, understanding and prediction. Due to the complexity of the computational methods, there is a strong need to integrate different models and cover the relevant scales in heterogeneous catalysis. This requirement constitutes an important drawback as scientists need training in several aspects of the problem including chemical, physical and engineering views of the modelling while keeping the experimental and industrial interests and needs in perspective. Here we present some of the latest developments in the field of theoretical simulations at the microscopic level while illustrating suitable examples that show how theory can shed light on several aspects of characterisation, activity, selectivity and long-term stability.

---

1 Introduction

The Holy Grail of theoretical simulations is the determination of a suitable stoichiometry and corresponding structures for a particular performance in a chemical process and the route to synthesise the active phase. The main goal of theory is to provide a deeper insight into the relation between structure and any of the three key parameters: activity, selectivity, and stability. The use of a black box set of programs containing all the theoretical models would be very practical, as that could shorten the time to commercialise new catalytic applications by designing suitable compounds and it would provide valuable suggestions for synthesis. At present, computational materials science and heterogeneous catalysis still have to overcome a large number of hurdles to reach this goal. Although the computational power accessible to theoretical simulations in heterogeneous catalysis has increased exponentially since mid-nineties and new tools and algorithms are being developed continuously, many important parts are still lacking in the present simulations. Moreover, not everything can or should be computed. Sometimes, the information stored in our detailed reaction networks may end up with a large number of elementary steps that are meaningless, or configurations with little, if any, role in the total reaction process and thus that do not contribute significantly to our overall knowledge. In this case, we need strategies to reduce the number of calculations and the computational effort required. Furthermore, the use of a single theoretical approach is not always enough and the coupling between different models is required. Particular effort needs to be devoted to the integration between different time and length scales. Therefore, the electronic structure obtained through Density Functional Theory (DFT) constitutes the keystone on which many of the longer space or time scale phenomena are built. The crucial role of DFT can be seen in Fig. 1 as adapted from the work of Vlachos.^1,2 It shows how the hierarchical structure of simulations in material science can easily be transferred to the theory in heterogeneous catalysis. The complexity of the catalytic phenomena requires the description of multiple time and length scales and thus, an intelligent way of adapting the relevant information from lower, more complete steps is crucial. In any case the detailed electronic structures need to be summarised and packed in a useful way for more qualitative approaches that, on the other hand, can reach larger time and space scales relevant in process and plant simulations.

Fig. 1 The time and length scales for different simulation tools. Adapted from Vlachos.^1,2

Recently, several reviews have been presented in the field of simulations in the lowest space and time scales. Special effort has been dedicated to describe all the properties of materials,³ to use linear-scaling relationships and screening^4,5 and lately to reach larger scales.^2,6 In the early days, the study of adsorption and the reaction energy profiles were determined for direct reaction that could convert reactants into products on a single surface. Nowadays, computational techniques can address not only the activity of different materials, but also the selectivity and the stability against reaction conditions. Desirable future goals would be (i) the implementation of filters based on descriptors to discard some of the targeted materials so that they fulfil industrial or laboratory requirements, and (ii) the ability to gather this information into usable databases for future data-mining.⁵

In the present perspective we aim at reviewing critically the latest achievements in the field, giving experimentalists a proper guide to assess the minimum requirements for the theoretical models and techniques that can be considered as reliable. In a second step, we will address the challenges theoretical simulations face.

2 State-of-the-art

Density functional theory

Density Functional Theory has achieved great success in computational catalysis. However, bridging gaps in the temperature, pressure, time and length scales accessible to DFT-calculations and those characteristic for real-world experiments still remains a challenge, even when using the most sophisticated codes. Moreover, the accuracy and efficiency of these methods depend on several technical choices,⁷ for example an exchange-correlation functional, the basis-set for the expansion of the Kohn–Sham orbitals and algorithms to solve the corresponding equations. The choice of exchange-correlation functionals and the completeness of the basis-set determine the accuracy, whereas the numerical algorithms determine the efficiency of the calculations. The minimum functional requirements for different systems and properties are summarised in Fig. 2. The hierarchy of exchange-correlation functionals that allows us to achieve an increasing accuracy of DFT results has been dubbed by John Perdew as “the Jacob's ladder of DFT”.⁸

The lowest rung of this ladder is the so-called Local Density Approximation (LDA),⁹ where the exchange-correlation energy (E_xc) for a homogeneous electron gas of the same density, as obtained from quantum Monte-Carlo simulations,¹⁰ is also applied to non-homogeneous situations. LDA has succeeded in solving many bulk¹¹ and surface problems. However, for chemical reactions occurring at surfaces, LDA usually leads to adsorbate over-binding.¹² Also, the potential-energy profiles for dissociations of diatomic molecules on metallic surfaces are badly characterised by LDA.¹³

The description is improved in the Generalised Gradient Approximation (GGA),^14,15 the second rung of the Perdew's ladder, which includes a dependence of E_xc on the local gradient of the electron density. The GGA heals the over-binding tendency of the LDA (although, with an inclination to over-correct).¹⁶ Moreover, for many chemical reactions, the GGA allows us to achieve sufficient accuracy.¹⁷ Indeed, this functional constitutes the proper description level for most of the reactions on metals, see Fig. 2. Several forms of the GGA have been proposed in the literature: PW91,¹⁸ PBE,¹⁹ revised PBE,¹⁶ PBEsol,²⁰ MA05,²¹ and WC.²² PW91 and PBE are by far the most commonly used functionals. However GGA's have two serious shortcomings. The first one is that they do not account for van der Waals (vdW) interactions that result from dynamical correlations between fluctuating charge distributions.⁷ The second problem, which arises from the approximate form of the exchange-correlation term, is the non-zero interaction of a single electron with its own density, which is known as self-interaction error (SIE). SIE is the cause of many of the failures of approximate functionals, such as too small band gaps,^23,24 wrong dissociation energies for molecules,²⁵ and a bad description of systems with localised f electrons. One approach to reduce SIE is the DFT+U method, in which Hubbard-type terms are added to account for the on-site Coulomb interactions in the localised d or f orbitals. The Hubbard parameter (U) can be fitted so that it reproduces experimental band-gaps, geometries or other properties. Unfortunately, U is not free from some arbitrariness, as fitting the parameter for one of the experimental terms does not ensure that the rest will be adequately reproduced. Some examples of oxide compounds that require this level of theoretical approximation are NiO and Ce₂O₃.²⁶ Therefore this is the minimum relevant description for reducible oxides, Fig. 2.

Fig. 2 Minimum density functional theory approximations for a set of catalytic materials and their characteristic properties.

Meta-GGA including higher-order terms of the gradient of the local kinetic energy density is the third rung on the DFT ladder.²⁷ Nevertheless, meta-GGA does not lead to a systematic improvement over the GGA, as shown by the ambiguous results for the adsorption of small molecules on transition- and noble-metal surfaces⁷ and for the stepwise hydrogenation of benzene to cyclohexene on a Ni(111) surface.²⁸

The next rung is represented by hybrid functionals that mix exactly, i.e. Hartree–Fock (HF) and DFT exchange, and describe correlation at the standard DFT level. The most popular hybrid in molecular chemistry is the B3LYP functional,^29,30 which combines LDA with HF exchange. The B3LYP functional achieves a very high accuracy for almost all properties of small molecules, but it fails when applied to metals and semiconductor solids, because the correlation part of the functional is incorrect in the homogeneous electron gas limit.^31,32 Extended systems are better represented by other hybrids such as the PBE0³³ and HSE03 functionals.³⁴ The results obtained for solids are unclear. PBE0 and HSE03 overall improve predictions for lattice parameters and bulk moduli of most solids, as well as for the band gaps in semiconductors and insulators,^31,32 moreover these hybrids offer an excellent description of insulating antiferromagnetic rare-earth and transition-metal oxides where the GGA failed to reproduce them correctly.^35,36 Atomisation energies, as well as magnetic metals, are described with a higher degree of accuracy through standard PBE.^31,32 In the latter case, hybrid functionals overestimate the exchange splitting and the magnetic moments, and broaden the d-bands,³⁷ therefore the overall description of the adsorbate–substrate complexes is not improved.³⁸ Another drawback of the hybrid-functionals is their very high computational cost.³¹ Thus, climbing up to the highest rung of DFT ladder does not guarantee a solution to all problems but hybrids reduce the SIE of pure DFT due to the mixture of a certain amount of exact exchange.¹⁶

As it has already been mentioned, dispersion (van der Waals) interactions are missing in standard DFT calculations. At present, a method that accounts for the vdW energy seamlessly and accurately is the random-phase approximation (RPA), combined with the adiabatic connection and fluctuation dissipation theorem (ACFDT).³⁹ Since this method is computationally demanding, it is mostly restricted to small systems and can serve as a benchmark for assessing the reliability of less sophisticated approaches. Lundqvist et al.^40,41 proposed a non-local correlation functional that accounts approximately for dispersion interactions. Results achieved with a nonlocal vdW-DFT correlation functional depend on the judicious combination of the local and non-local contributions to the functionals. A simpler dispersion correction is offered by the semi-empirical force fields of Grimme et al. (DFT-D2).^42,43 In this method, the dispersion contributions are calculated by pair-wise interactions from the London formula⁴⁴ leading to the ΣC₆/R⁶ term. Recently, Bučko et al.⁴⁵ applied Grimme's method to a large number of solids showing a wide range of chemical bonds such as molecular, ionic and covalent. This illustrates that Grimme's approach can be used for arbitrary systems irrespectively of the nature of dominant interactions. Moreover, a careful choice of the substrate interaction coefficient C₆^46,47 provides a good starting point to study the adsorption of molecules and films on conducting metal surfaces⁴⁸ and, potentially, on ionic or semiconducting solids.⁴⁷ Thus the Grimme or similar approaches constitute the minimum meaningful theory level when van der Waals interactions are important, i.e. in layered compounds (as PtO₂), or in physisorption, see Fig. 2.

As it was mentioned at the beginning of this section, standard calculations refer to systems under vacuum conditions. However, recent progress allows researchers to overcome these restrictions.⁴⁸ DFT combined with statistical mechanics in the grand canonical ensemble describes an adsorbate–substrate complex in equilibrium with a reactive atmosphere. Temperature effects are treated as corrections to static total-energy calculations by adding the vibrational free energy to the total energy from DFT calculations or using first principles canonical molecular dynamics simulations.⁴⁹

Models for the materials

The materials involved in heterogeneous catalysis present a wide range of properties. Metals, semiconductors and insulators are employed in many cases in mixed configurations, thus adding an extra degree of complexity; hence not in all cases the same electronic structure approximation is optimum for the different constituents, see Fig. 2. Moreover, materials can be amorphous or crystalline, or have grain boundaries or porous structures that render the systematic assessment of their properties difficult. In general, different types of materials are studied by several approaches: while the chemistry of amorphous systems relies on the very local configurations that can be achievable, the description of crystalline materials is mostly based on their periodicity. The study of such materials has advanced very significantly and as a first approximation some of the properties of amorphous structures can be simplified to the periodic ones. For this reason we will focus on them.

A problem of DFT with periodic boundary conditions is caused by the complex lattices with several potential structures and/or structural disorder due to finite temperatures. In the first case, the as-calculated energies do not contain thermal contributions to the free energy that are needed for some delicate cases when multiple minima lie close in energy. A good example of this is the re-evaluation of the structure of dawsonite (MAlCO₃(OH)₂ M = Na, K, NH₄⁺) compounds.⁵⁰ In this particular case, the experimental determination of the crystal structure is performed at finite temperatures at which the ammonium cations might rotate. The introduction of this thermal contribution is usually done by accounting for two different terms the zero-point vibrational energy and the contribution to the energy through the phonon distribution in the form of a summation.⁵⁰ This approach is also listed in Fig. 2, and constitutes the minimum description when trying to assess structures and reactivity of complex materials with close-lying structures and important phonon effects. Other catalytic solids exhibit some degree of disorder in the crystal site occupancies (this is different from amorphous disorder!). The models employed to deal with this kind of inhomogeneity can be classified into three broad groups:⁵¹ (i) methods in which a sort of average atom is defined. They allow recovering the perfect periodicity as implemented in the virtual crystal approximation (VCA), where the potential felt by electrons is generated by average atoms. (ii) Methods involving a large supercell with a more or less random distribution of ions in the sites. A useful variation of this model is the special quasi-random structure, where the ion positions in the supercell are chosen to mimic as closely as possible the most relevant near-neighbour pair and multisite correlation functions of the random substitutional alloy.⁵² (iii) Multi-configurational supercell approaches, where the site-disordered solids are described as a set of configurations in a supercell representing a piece of the solid (e.g. as in ref. 53). This kind of approaches is useful to identify favourable ordering patterns,^54,55 to evaluate the evolution of ion disorder with the temperature,⁵⁶ or to examine the segregation of impurities at solid surfaces beyond the dilute limit.⁵⁷

Once bulk structures are determined, selected cuts along the lower Miller index directions are performed and then the surface energy for the faces with different orientations is calculated. The Wulff construction,⁵⁸ a simple continuous model, determines that the lowest surface energy will be the most represented in the equilibrium structure of the metal nanoparticles in the catalysts. Extensions to materials covered by surfactant or molecular modifiers have been put forward.⁵⁹ However, the Wulff construction presents several approximations: (i) only surface energies are taken into account, i.e. the contribution from edges and other defects is overlooked and thus the representation is more accurate for medium size particles; (ii) only equilibrium structures are retrieved. Non-equilibrium structures are much more difficult to assess but, on the other hand, they are also difficult to characterise experimentally, thus this area remains somehow obscure to both experimental and theoretical worlds.

The reactivity of solid materials is usually first studied theoretically on the low Miller index facets that usually show lowest surface energies. This is at odds with the corresponding reactivity as the ability to adsorb new molecules does depend on the number of bonds of the substrate. Still several reasons justify this choice: first of all, these low energy surfaces are the most represented by the Wulff construction; moreover, the unit cells of low-index surfaces tend to be smaller than open surfaces or vicinal ones (so the number of different configurations for adsorbates is smaller and thus the number of possibilities to be investigated are lower) thus it reduces the calculation burden. This is the minimum calculation setup for a simple active phase (Case I) shown in Fig. 3.

Fig. 3 Computational models for the interaction of metals and oxides (or other supports) as a function of the type of interaction.

In some cases, the adsorption or the activation of reactants on these surfaces is not strong enough, Case II in Fig. 3. This is for instance the case of the dissociative adsorption of nitrogen on Ru surfaces that is the key step in the synthesis of ammonia in the Haber–Bosch process.¹⁷ Early in the investigation of this process, discrepancies between ultra-high vacuum measurements of adsorption and the barriers calculated by DFT were found, but after employing theoretical models containing steps the experimental value was retrieved. The reason for the highest dissociation activity of the open surface was correlated to the presence of B₅ sites that provided a transition state structure where no atoms were shared.⁶⁰ Since then, the suitable procedure when investigating a new reaction begins determining molecular and dissociative adsorption on low index surfaces. Then, if the reactivity of the surfaces seems exceedingly poor, the presence of defects is included in the study, usually starting with geometric perturbations, such as low-coordinated sites,⁶¹ atomic vacancies and/or considering the presence of impurities from a secondary metal,⁶² or fragments from some of the interacting molecules.⁶³

With the approaches presented in the previous paragraph the representation of the middle size particles is reasonable, but still two questions remain open. First of all, when modelling finite size nanoparticles, the number of atoms that can be included in the simulations is rather small, few hundred atoms at most. Secondly, the shape of the particles also has a contribution and this is somehow even more difficult to assess. Recently the convergence of the properties with the size has been analysed.⁶⁴ Even if the convergence of the metal properties is found at relatively small number of atoms in the cluster, the mobility of the surface atoms generates more dynamic properties as the cohesive energy of the surface atoms is reduced when compared to the bulk. This problem is more acute for metals like Au, for which the melting temperature is quite low.

Although we employ the cluster as a particle of 2–5 nm in diameter, the representation is far from complete. In the preparation of catalysts, there is special interest in obtaining monodisperse samples, at least in research and for characterisation. The aim would be to obtain only particles with a given diameter but this is impossible when employing chemical methods. The alternative, use of physical methods (chemical vapour deposition with charge mass selection)⁶⁵ with very detailed separation, is scientifically relevant but not viable. The reasons for that are the high cost and the low surface areas of these model catalysts. Developments are still being made in this regard, and for instance, spray deposition controls particle size much better with smaller deviations.⁶⁶ There is a hidden issue that concerns the analysis of the size and shape of the nanoparticles, in general, Transmission Electron Microscopy even in high resolution experiments cannot adequately determine the particles below 1 nm in diameter. The fact that experiments are not accurate enough poses two different problems. In principle, the activity of the catalyst, when linked to the presence of low-coordinated sites, increases when decreasing the particle size as the number of exposed atoms with respect to the total metal content is larger.⁶⁷ This is a reason for the preference of low-diameter particles. But for very small nanoparticles, activity is not a continuous function as the electronic structure of small aggregates has a discrete spectrum and their fluxionality (conformational flexibility) plays a crucial role. This is known as the non-scalable regime where each atom counts. Experiments have shown that for instance CO oxidation is possible for Au₈⁶⁸ and Au₂₀, whereas Au₅₅ is inert.⁶⁹ The analysis of the reactivity of such small clusters does then depend on their electronic structure and, as a consequence, requires a tedious examination and the results from one cannot be extrapolated to the ability of the next.⁶⁵ From a theoretical point of view this opens up millions of possibilities, but at the same time it is a never-ending story as it is almost impossible to tailor the properties for all systems.

Another of the major questions when dealing with theoretical simulations in heterogeneous catalysis is the role of the support. Few model systems have been employed and presented in the literature. The problems with these materials are related to the additional demands required by electronic structure calculations. One of them is the need to use an extremely large unit cell to account at the same time for the active nanoparticles and the substrate.⁷⁰ Also, as the metal nanoparticles are small their transferability is more compromised. Moreover, the most common carriers are significantly difficult to model. Aluminas,⁷¹ like the γ phase,⁷² present unsolved issues; silica also poses problems due to its amorphous nature that cannot be obtained from typical DFT and periodic boundary conditions or in terms of the degree of hydroxylation.

The way the substrate affects the properties of the active phase can be separated from lower to higher interaction, see Fig. 3.⁷³ For instance, even weakly interacting substrates influence the dispersion and shape of nanoparticles. In this case, not only the stoichiometry of the support surface but also mesoscopic parameters such the Brunauer–Emmett–Teller (BET) area or porosity are crucial to influence the size and shape of the metallic nanoparticles. This control is quite indirect as the dispersion of the metals on the substrate is determined by the ability of the atoms to be anchored to the surface by particularly active carrier positions. Therefore, in weak interaction to the substrate Case I and II models are the first order approximation, see Fig. 3. If the interaction is stronger, dispersion of the active phase occurs. In the case of very small particles, even charge transfer can take place although screening for these materials is very effective and thus charge contributions are accommodated for medium-size particles. As the stoichiometry induces the shape and the dispersion, there is a reasonable parameter that needs to be eliminated when comparing different catalysts. In this case, the inclusion of the carrier in the DFT model is not required as its effect can be introduced indirectly. A particular case of moderate interaction between the active-phase and the carrier takes place for some catalysts as that employed in the Deacon reaction (HCl oxidation to Cl₂ that usually is carried out by RuO₂-based materials), where such strong interactions occur and the active layer grows epitaxially on the carrier. This corresponds to Case III in Fig. 3. A single monolayer has such a large interaction with the support that becomes inactive and thus a few layers are required to eliminate the electronic perturbation induced by the support.⁷⁴

Another type of metal–support interaction can induce the formation of special sites at the interface between the active phase and the carrier, see Case IV in Fig. 3. Possible active sites that can benefit from the synergetic interaction have been put forward. CO oxidation on gold nanoparticles is a good example: CO adsorbs on the metal, and, simultaneously, oxygen can be activated by the partially reduced sample.⁷⁰ Some issues arise as the modelling of such bifunctional mechanisms requires adequate description of both the active phase and the support. In the submonolayer regime oxide growth on oxides has been observed for vanadia compounds on reducible oxides such as TiO₂ and CeO₂.^75,76 In both cases synergies between the small vanadia clusters and the oxide supports have been identified in experiments and described theoretically.

A final scenario can be envisaged where part of the reaction takes place in the carrier and spill-over to the (in principle) active phase takes place, see Case V in Fig. 3. In this case, simulations can be performed individually and the results can be combined to model the complex system. Similarly, the stability of the products shall be investigated against the acidity and basicity of the support, which can have an effect at this stage.

For very reducible oxides and/or high temperatures, a dynamical behaviour of the systems has been reported. For instance, some carriers are oxides and can lose oxygen during the reaction. Then, the active phase can change shape during the reaction and thus the number of active sites, their nature and the ability to generate mixed sites at the interface, significantly affecting the activity of the material. This is the particular case of Cu/ZnO,⁷⁷ the catalysts developed for the synthesis of methanol from CO and H₂. Still other dynamic behaviours, in particular related to the appearance of the strong-metal–support interaction,⁷⁸ can severely affect the catalytic properties and hinder the generation of a simple model for the catalyst. Obviously, these cases are the most complex model due to the long-time and long-length scale changes under true reaction conditions.

State of the catalyst

Preparation methods and pretreatments might affect the state of the catalyst. In the synthetic process a typical cleaning of the precursors is the reduction of the material. In particular, hydrogen atmospheres can change the relative stability of different surfaces. The history of the sample thus influences the ability to catalyse a given reaction. The study of the stability is usually addressed through first-principles thermodynamics which include not only the calculated energies but also the contributions of the pressure and the temperature and this constitutes the minimum computational set to address these properties, see Fig. 2. An illustrative example is given by the activity of Pd-based catalysts in the hydrogenation of alkyne–alkene mixtures in the treatment of feeds coming from oil reforming. The key issue in this case is that Pd behaves as a sponge and, in particular, when hydrogen or hydrocarbon atmospheres are present, it can generate hydrides or carbides.⁷⁹ The presence of hydrocarbons favours the formation of carbon and hydrocarbon surface deposits (coking).^81,82 Examples exist also in oxidations, like the epoxidation of olefins carried out with the help of a silver-based catalyst, Ag might convert at high oxygen conditions into the substoichiometric oxide.⁸⁰

The formation of these competitive phases and subproducts has been investigated with theoretical methods using DFT on a list of structures that include a small set of configurations containing both the main active phase and the atoms from the gas-phase, (reactants or products). Upon characterisation of such configurations by DFT, first principles thermodynamics shall be employed to account for the temperatures and the pressures of the components.⁸³ The equations that aim at obtaining the Gibbs or Helmholtz free energies are related to the chemical potentials of the gas-phase compounds and thus the corresponding surface free energies for all configurations are obtained.

A few simplifications are employed to determine the free energies. For instance, the contributions from the phonons are neglected. This approach helps to evaluate the order of the introduced error. In addition, the gas-phase energies are taken from the ideal gas approximation, and either the experimental entropies or the ideal-gas statistical mechanics approximations are introduced. Configurational entropies are completely neglected or only introduced as an added fix in an approximate manner.⁵⁶ Even if first principles thermodynamics provide a way to derive the equilibrium structure it is not free of error sources. The problem is that configurations usually are generated by hand, and therefore some of them that do not seem obvious a priori might not be introduced in the calculation pool, thus missing representativity. Genetic algorithms can be employed to avoid this issue but unfortunately their introduction in the field is not straightforward, since the codes are usually fitted to the particular problem to be investigated and require a greater deal of effort.⁸⁴

As first principles thermodynamics usually neglects a part of entropic contributions and rely on adsorption and desorption processes, the errors associated with the determination of pressures and temperatures are somehow larger than taking into account the whole reaction mechanism. Therefore, this technique is powerful and in particular qualitative when oxygen is one of the active gases and PW91 or PBE are employed to obtain the energies.

A similar situation might occur under operation conditions. Several catalysts are known to present an induction time where the as-prepared phase turns into the active phase. Moreover, some other materials that are active for a period might be inactive after being under reaction conditions for a sufficiently long time. This might occur by a change in the size and shape of the active phase, its formation or collapse, its volatilisation or the formation of carrier layers on the active part, just to name a few possibilities. The experimental detection of the changes taking place during this activation or deactivation is in general quite difficult as in several cases in situ experiments are required under reaction conditions. In some of these cases, first principles thermodynamics turns out to be a good tool to assess the state of the catalyst with respect to the reservoirs of the active gases that will be present in the reaction mixture.⁸⁵ Such approach can help filtering the structures that stay in the active phase under reaction conditions and also assessing the real surface termination. From our experience, we have seen that stability against harsh conditions might be the most important parameter for real applications. Whereas industry can cope with not so effective activities (as engineering solutions may reduce the low-activity impact), there is usually a much smaller toolbox to improve stability of the catalyst under reaction conditions.

Typical catalytic formulations under real synthetic conditions are somehow modular and, besides, from the base metal a secondary metal and/or one or more dopants and/or molecular (selectivity) modifiers are added. When dealing with complex materials with more than one component, the stability of the alloy both by itself and also under reaction conditions is mandatory. These parameters related to stability can be obtained by theoretical considerations and for instance a group of descriptors including, solubility, segregation and induced segregation energies, island formation, decoration of low coordinated sites and substitution at step sites among others can be identified. Examples of these have been put forward by several groups including ours and they can be extended to oxides and other materials.^62,86 Additives can be electron donors or acceptors that affect the charge at the surface,^87,88 or alternatively they can help to stabilise particular stoichiometries or geometries. Finally, molecular modifiers might affect adsorption energies for thermodynamic selectivity purposes, or reduce the size of the ensembles that might lead to selectivity, as in the Lindlar catalyst preparation.⁸⁹ In all these cases, theoretical simulations can provide an understanding on the individual role of the different components to the activity–selectivity–stability of the catalysts under reaction conditions.

Reaction mechanisms and kinetic parameters

The study of the reaction networks relies on linking the various potential configurations of reactants, intermediates and products on the surface and with respect to the gas-phase. In this way different mechanisms including these sets are investigated. The most robust algorithms to obtain the transitions state structures that link reactants and products belong to the family of the Nudged Elastic Band (NEB) methods.⁹⁰ In general, once both the initial and final states are known, a series of images are placed between these positions and then the optimisation of the full band takes place. We need to consider, firstly, that the computational requirements for these calculations correspond to those of a static calculation multiplied by the number of images in the band. Thus, these are the heaviest calculations when modelling reactions on surfaces. Secondly, that the use of the traditional NEB does not ensure that the highest point along the reaction path is the true transition state (TS). A vibrational analysis of the structure is needed in order to assess correctly that the point fulfils the mathematical definition of a TS (i. e. that it is a minimum along all possible directions except for that of the reaction coordinate, for which it is a maximum). This issue was solved with the advent of the Climbing Image version of the NEB algorithm (CI-NEB),⁹⁰ a much more robust way to obtain the TS structure and energy. However, as it has been mentioned before, calculations are still very demanding and thus algorithms aiming at shortening calculation times and computational requirements have been proposed. Examples of those are the Adaptative Nudged Elastic Band Approach (ANEBA),⁹¹ where instead of choosing a large number of images to bracket the saddle point with high accuracy, the resolution in the neighborhood of the saddle point is increased. Other appealing methods are based on the dimer method,⁹² which is designed to search for saddle points corresponding to unknown elementary steps, where no final state is required and the TS structure is proposed; then upon calculation of the Hessian, the movement along the Potential Energy Surface and towards the TS structure is done by following this particular reaction coordinate. In any case, the search for TSs constitutes the bottleneck in the description of catalytic phenomena on surfaces, therefore further research into the algorithms will be required to make the calculations faster. As the number of databases and calculated points increase every day the inheritance of parent structures for the search of transition states in new more complex molecules shall be a way to reduce computational costs.

In many cases the reaction networks are constituted by a short list of reactions. This is the case of CO oxidation, for which two paths have been put forward and either the dissociative O₂ path or the associative one can take place, depending on the material.⁹³ This mechanism is relatively simple as only four elementary reactions need to be taken into account. When dealing with larger molecules or more complex mixtures, the situation becomes increasingly complex. For instance, in the selective hydrogenation of alkyne–alkene mixtures, the minimum set of reactions needs to take into account the following paths: (i) hydrogenation,⁹⁴ (ii) isomerisation of the intermediates on the surface,⁹⁵ (iii) carbide formation by decomposition of the hydrocarbon at low-coordinated sites,⁹⁶ and (iv) oligomerisation⁸² and/or carbon deposits formation. Several of these steps, in particular, those that reduce the activity, selectivity, and stability (as (iii) and (iv)) are likely to be common to different transformation and thus they could, in principle, be retrieved from databases (vide infra).

In other cases, the reactivity needs to be calculated for a series of potential surface structures including the corresponding carbides, hydrides and selectivity modifier molecules. This is also the case of the hydrogenations on Pd. Then the determination of the TS paths needs to be considered for each potential structure in order to derive the structure–activity (selectivity) relationships.⁸²

Still another scenario unfolds when considering sets of reactions that can be performed under different conditions. This is the case of HCN synthesis from ammonia and methane on PtRh catalysts, see Fig. 4.⁹⁷ The reaction can take place in aerobic (the one implemented in industry) and anaerobic conditions. The presence of oxygen atoms on the surface systematically reduces the energy demand of the H-dissociation steps but the price to be paid is that in the Andrussow oxidation the final C–N bond formation takes place from very strongly adsorbed atoms on the surface. Thus the main benefit found for the reaction thermodynamics when O₂ is present does not result in a large reduction in the kinetics or reaction temperatures (only 300 K lower), as there is a mechanistic switch for the C–N bond formation from partially hydrogenated molecules (in the non-oxidative case) to more bound atomic species (oxidative conditions), see Fig. 4.⁹⁸

Fig. 4 (a) Reaction barriers for the coupling of different C and N-containing fragments. (b) Schematic representation of the two different possibilities for the HCN formation under O-lean and rich conditions.

As shown above, the list of potential reaction paths and parallel elemental reactions is quite large even for simple reactions. In order to predict activities and selectivities, the amount of information stored in such calculations needs to be reduced to the most relevant parameters (descriptors) that can affect the reactivity (activity or selectivity) by modification. Early descriptors, based on the activity of oxides, were already devised by Sabatier who confronted the reactivity of different transition metals towards the decomposition of formic acid against a thermodynamic parameter, the formation energy of the oxide.⁹⁹

Over the last years, in particular in the group of Prof. Nørskov in Stanford, two different ways to describe relationships between energy parameters have been developed. The first one corresponds to the linear-scaling relationships.^100,101 In these, the binding energy of a fragment to a metal (also for oxides, nitrides, sulphides, carbides) correlates with the binding energy of the central atom to the surface. Moreover, the proportionality factor depends on the valence. Therefore, a simple explanation in terms of the density available to form different bonds and the valence of the heteroatom can be traced back in a simple yet effective manner. A second scaling relationship corresponds to the Brønsted–Evans–Polanyi rules, BEP.^102,103 BEP states that the activation barrier depends on the reaction energy for this particular elementary step. Dissociation reactions are in many cases responsible for the activity found, as they constitute the rate-determining steps.^104,105 In such situations, BEP can be further simplified to the adsorption energy of the central heteroatom. BEP relationships are an extension of the linear scaling relationships described above. With these two approaches, many mechanisms can be simplified to just one or two significant descriptors.⁹³

Microkinetic modelling

In many cases, dealing with the description of the complete reaction energy profile leaves open questions such as which are the rate-determinant or selectivity-determinant steps. This answer requires the modelling of the kinetic equations with the correct boundary conditions in order to be fully meaningful.

In general, the data obtained from Density Functional Theory with the appropriate approaches (vide supra) can be taken as starting point for the microkinetic (MK) simulations. The rate coefficients are normally obtained through the use of standard statistical thermodynamics and by using the harmonic transition state theory.^106,107 The list of kinetic elementary reactions, the site balance, and the boundary conditions introduced describe a system of Differential and Algebraic Equations that can be solved by different mathematical packages. In any case, the activity of the system can be described, the reaction orders retrieved, and the apparent activation energies and coverages of different species unravelled.⁷⁹ This completes the knowledge on the reacting conditions and, once the true path has been obtained and the relevant surface is taken into account and different reaction conditions can be easily plugged in. In some cases accuracy problems and unstable solutions might arise and thus this type of approach is not as straightforward as standard DFT packages.

Moreover, the reaction mechanism (understood as the list of reactions) does not change for the same process on a large number of catalysts. Therefore, taking into account the considerations made in the previous section, it is sometimes useful to simplify the problem by transferring all the relevant kinetic terms to the descriptors. Successful examples of such approach have been presented by the group of Prof. Nørskov for a number of processes including CO oxidation on metals and nanoparticles.⁹³ In particular, the Sabatier analysis with more than one descriptor has produced multidimensional volcano plots where the maximum of activity can be inferred. Another important point of microkinetic approaches is that it allows determining the weight of each step that composes the mechanism and the analysis of the role of the different elementary reactions.^6,108 Finally, MK modelling can assess important aspects such as the adequate analysis of the population of different species on the active substrates directly comparable to experiments. With such an approach it has been possible to determine why compensation effects occur in heterogeneous catalysis.¹⁰⁹

In some cases, the classical MK has some deficiencies. Those are related to the lack of geometric information in the equations which limits the applicability of this technique, in particular when confined systems are considered.¹⁰⁸ Kinetic Monte-Carlo methods are better suited to deal with these in homogenous lattices, as the geometric information is already available. The reader is addressed to the Perspective by Reuter and co.⁶ on the subject for a more detailed description.

3 Main challenges

Problems when addressing large/flexible molecules

Until recently, only the interaction of small molecules on metal systems was considered in theoretical simulations. Many heterogeneous catalysts deal with the activation of small molecules and their incorporation to other larger fragments (i.e. partial oxidations and hydrogenations) or the reshuffling of the internal bonds, like in the case of ammonia synthesis. Indeed, CO oxidation has been the favourite reaction for theoreticians so far.

The search for catalysts that can facilitate the implementation of new industrial processes faces here an important challenge. In the past, many of the interesting active molecules were obtained from oil. These active molecules were small and had a functional group at most. Examples are methanation, hydrogenation of alkynes in alkene mixtures, and even the Fischer–Tropsch reaction. Unfortunately, as these natural sources are being depleted fast, we will need to change gears and start employing larger molecules that present several functional groups. This poses many challenges to the modelling of the reactions that occur, the most important of which are: (i) the presence of multiple adsorption configurations due to the functionalities, (ii) the inordinate amount of potential parallel paths that need to be described with similar accuracy in order to properly obtain the selectivity, and (iii) the liquid nature of some of these compounds due to their high oxygen content. Therefore, even the investigation of the adsorption for all relevant configurations can become a challenge.¹¹⁰ Moreover, when attempting microkinetic modelling a better assessment of the number of positions occupied (i.e. improved definition of the adsorption site) would be required just to account properly for the larger size of the molecules.¹¹¹ Two examples of how to address this particular problem have been put forward by the groups of Greeley and Vlachos.^112,113 In the first case, the study of glycerol decomposition on Pt(111) was estimated by an empirical correlation scheme later redefined by DFT calculations coupled to BEP analysis for the dehydrogenation and C–C bond scission. Salciccioli and Vlachos proposed a parametrisation of the thermochemical properties of C₂H_xO₂ intermediates and transition states for different bond breaking patters C–C, C–O, and O–H in the form of a functional group approach.

Flexibility poses also an important problem when looking at large molecules adsorbed on surfaces. In homogeneous catalysis, extensive searches for the most stable conformation of the catalysts, reactants, intermediates and products are performed routinely. This is not the case of the simulation of large molecules on surfaces. Just to give a small example, the rotation of a methyl group can increase the Pauli repulsion by 0.3 eV with respect to the same methyl adsorbed on the surface but with no direct H-interaction on Ag.¹¹⁴

Solvent effects

Solvent effects can have a strong influence on the study of catalysis¹¹⁵ due to a non-homogenous electric field that perturbs the electronic clouds of the solute and therefore its geometry, changing its properties when compared to the gas-phase species.¹¹⁶ Solvent effects have been mainly characterised in homogeneous catalysis. Two main general approaches are commonly used: classical ensemble treatments and quantum mechanical continuum models.

The former include classical molecular dynamics simulations,^117–119 Monte-Carlo techniques,¹²⁰ free energy perturbation,^121–123 and Langevin dipole moments.¹²⁴ In many cases, solvent is treated explicitly as a rigid system with vdW ε/σ parameters and partial charges,^125–127 and the interaction between molecules is handled by pair-wise interactions between atoms.¹²⁸ These models give useful information about solute–solvent interactions, and describe conformational changes carefully. Their main drawback is the need for an appropriate potential function, extrapolated from other data, as well as the lack of a proper relaxation of the geometry and the dipole moment during the simulation.

The latter approach has its roots in the Onsager reaction field model.¹²⁹ According to this, the solute is placed in a cavity¹³⁰ that can be for the whole molecule or a summation of overlapping spheres centred on each atom (PCM)¹¹⁶ immersed in a continuous medium with a dielectric constant ε. A dipole in the molecule will induce a reflection dipole in the solvent, and this, in turn, (reaction) will interact with the molecular dipole, leading to a net stabilisation. Given the structureless nature of the solvent, this approach is best suited to describe apolar systems or solvents where a specific interaction with solute (such as H bonds) is missing. Unlike classical field methods, no additional information concerning potential functions and related parameters is needed. In PCM, the surface potential is calculated by numerical differentiation, and its interaction with the solvent can then be computed self-consistently. This is equivalent to carrying out the dipole expansion to infinite order, strongly improving the Onsager model. Further developments are DPCM (dielectric),¹³¹ CPCM¹³² and IEFPCM.^133–135

In heterogeneous catalysis, the focus has been on gas-phase processes and thus the introduction of solvent effects is much more recent. For instance in porous materials, like Metal–Organic Frameworks, (mixed organic inorganic compounds proposed to have tantalising properties, MOF) it has been possible to employ Born–Oppenheimer Molecular Dynamics to assess explicitly the role of water^136,137 or even Car–Parrinello calculations on the formation of the first aggregates for zeolite formation¹³⁸ (vide infra). In such cases, the solvation energies of atoms are taken properly into account for electrochemical purposes.¹³⁹ Water dynamics in the Car–Parrinello approaches has also been employed to study the relative stability of surfaces in semiconductors.¹⁴⁰ On metals, the situation is far more complex as first principle dynamics cannot be performed efficiently. In general, explicit water structures in the submonolayer regimes have been studied,¹⁴¹ and in some cases tridimensional ice layers have been employed.^142,143

Theoretical simulations have faced problems even in the description of such layers as many minima with similar energies exist and the proper evaluation of the dispersion contributions is needed.^{48,141,144,145} Still, models based on systems with different ice configurations filling the space between two metal slabs can be employed and although the model is naïve and quite rigid, some important results have been obtained.¹⁴⁶ Following the reasoning in homogeneous catalysis, it would be worth attempting a similar approach in heterogeneous systems through the definition of cavities and continuous medium. So far, few initiatives with very few results have been presented in the literature in this direction.¹⁴⁷

On the transformation of theoretical chemistry into an information technology

As computational approaches can systematically investigate different materials with a similar degree of accuracy it is possible to generate structured information in the form of databases. Over the last years the amount of information gathered has been increasing exponentially and a development is required to transform this Chemistry knowledge into a true information technology. As a consequence, several groups have started to consider the possibility of generating databases that could be mined when a new particular problem arises. Of course, several questions are open about how to structure the information to make it useful and relevant. The main challenges are: determining the main tags to be included and, for open repositories, how to assess the quality of calculations coming from different groups or subjects. However, the latter problem is also common to other collaborative information systems, and therefore it should be easier to solve.

As for the already existing databases, we would like to make special mention of those created by the group of Prof. Ceder in MIT which is devoted to the study of batteries,¹⁴⁸ and that of Prof. Nørskov's dedicated to heterogeneous catalysis.¹⁴⁹ A more general discussion on the role of these databases is presented by Prof. Lüthi and co. in ETH-Zürich.¹⁵⁰ Finally, the computational groups at ICIQ (Professors Bo's, Maseras', and Lopez') are currently developing a similar tool: SCIPIO.¹⁵¹ A snapshot of our web platform, which is compatible with different kinds of software (VASP,¹⁵² ADF,¹⁵³ and Gaussian,¹⁵⁴) is shown in Fig. 5.

Fig. 5 Snapshot of the SCIPIO storage and database tool developed at ICIQ (www.scipio.iciq.es) showing MOF-5 results.

From synthesis to structure

The hardest goal in the field is to actually evaluate preparation methods that could eventually lead to the desired molecular architectures known to be particularly active, selective, or stable. Obviously, this is a multiscale problem where the complexity in the number of different steps in the preparation of a catalyst is a long-term challenge. For instance, surfactant chemistry has generated a new route for the monodispersion of gold-based catalysts. Still elimination of these surfactants while keeping the chemical functionalities is still a challenge.^155,156 The investigation of the role of surfactants in the final structures, i.e. how they can introduce modifications in the typical Wulff structures, is still in its infancy.

Other examples of potential contributions in the field regard the synthesis of zeolites and the mesopore formation by water or by hydroxides. Experimentally, lots of evidences have been gathered concerning the position of different atoms in the ordered lattices or how the overall porosity can be increased.¹⁵⁷ However, important points still remain unclear. The remarkable study of Van Santen and co. on the initial mechanisms in the polymerisation of siloxane groups to configure zeolites, represents a true landmark in the way complex problems need to be addressed.¹³⁸ In this study, DFT calculations performed with solvent are accompanied by a Kinetic Monte-Carlo approach that extends the length- and time-scales of the study. Similar studies on the thermal stability of MOF also require at least the combination of Monte-Carlo simulations to assess the configurational contributions, and first principle molecular dynamics to obtain the true relevant parameters in the water induced lattice disruption.^136,137 The approaches above constitute a new direction where simulations can be powerful, this is the synthesis to structure challenge.

Towards a general theory of catalysis

The search for a new complete theory of catalysis has been under discussion for decades/years. In 2006, Nature¹⁵⁸ presented a list of what chemists wanted to know and, taking the sentence by Berthelot “Chemistry creates its object”, it discussed the main issues to be developed in the field. Indeed the first question regarded the design of molecules/materials with specific properties and, therefore, chemistry was identified as “a science of particulars” including that “it would be ludicrous to look for a general theory of catalysis that applies to all enzymes, materials surfaces and so on”. The methods employed by theoreticians to study all the potential materials are very similar in all cases and thus the formulation of parallelisms is possible. The main questions are still related to the ability to understand the typical jargon used by scientists with different backgrounds: such as engineering, physics, chemistry or organic chemistry.

Despite these critical assessments there have been several attempts to prove that the answer is not so straightforward. For instance, during the gold rush that started in 2000,¹⁵⁹ several processes that were proved to take place on gold nanoparticles were transferred to organometallic compounds and conversely. In turn, ruthenium is known to generate epoxides when prepared as organometallic and nanoclusters (known as polyoxometalates, POM) but these properties are not retrieved for ruthenium oxide. Experimental comparison between the activity of organometallic and Ru surfaces for polymerisation has also been put forward.¹⁶⁰ Last but not least, some sulphides that recall the structures present in the enzymes turn to be active in the electrochemical conversion of protons to hydrogen.¹⁶¹ Thus, the number of proofs is important, but the degree of transferability of the chemical properties from isolated to tridimensional structures is not straightforward.

Analysing the similarities and differences in these materials is possible by employing an equivalent level of theory to enzymes, organometallic systems and surfaces. We have studied a couple of these systems. For instance, the nature of gold catalysis is based on different properties.¹⁶² The differential adsorption of alkynes leads to selectivity, but organometallic alkenes are more likely to be coordinated to gold cations. Thus the similar chemistry, i.e. alkyne activation, has two different origins either a thermodynamic (heterogeneous) or a kinetic (homogeneous) one. As for the epoxidation of alkenes it occurs selectively on many Ru based homogeneous catalysts¹⁶³ and in nanoclusters¹⁶⁴ but the reaction is completely unselective for the surface oxide.¹⁶⁵ The reason for the selectivity is linked to the tridimensional structure and the curvature. One of the early stages of this reaction is given by the coordination of oxygen. On RuO₂, oxygen adsorbs dissociatively and the chemical potential of oxygen is high enough to prevent the formation of the intermediate. Curved surfaces as in POM compounds result in site isolation as in the organometallic compounds, thus also keeping the inability to split oxygen which improves the selectivity.

4 Conclusions

In the present perspective we have reviewed in a critical way the strengths and drawbacks of state-of-the-art theoretical methods when used to describe heterogeneous catalysis. The increasing capacity of computers and algorithms has opened a new field with great opportunities but several bottlenecks for the spread use of such integrated methods are still present. We also consider that training of both experimentally- and theoretically-oriented scientists in the abilities and limits of the present techniques will improve the synergies between both, enhancing the progress in the near future. The increased value that industries are giving to theoretical simulations is a clear demonstration of the future role of these techniques.

As for the challenges that the theoretical simulations of catalytic systems face, the integration of different computational methods and the ability to introduce larger molecules, solvent effects and transfer the practical results coming from several methods is a must. Finally, the creation of common databases that condense the results gathered by several groups and the identification of synthesis to structure rules would certainly help to generate a more general theory of catalysis get closer to application.

Acknowledgements

We thank the MICINN for projects CTQ2009-07753/BQU, CSD2006-0003, ERC-Starting Grant Bio2chem-d 2010-StG-258406, and BSC-RES for providing generous computational resources. We would like to thank Prof. J. Pérez-Ramírez, Dr Grau-Crespo and O. A. Salawu for useful discussions.

Notes and references

1. D. G. Vlachos, AIChE J., 2012, 58, 1314–1325.
2. M. Salciccioli, M. Stamatakis, S. Caratzoulas and D. G. Vlachos, Chem. Eng. Sci., 2011, 66, 4319–4355.
3. E. A. Carter, Science, 2008, 321, 800–803.
4. J. K. Nørskov, T. Bligaard, B. Hvolbaek, F. Abild-Pedersen, I. Chorkendorff and C. H. Christensen, Chem. Soc. Rev., 2008, 37, 2163–2171.
5. J. K. Nørskov, T. Bligaard, J. Rossmeisl and C. H. Christensen, Nat. Chem., 2009, 1, 37–46.
6. M. K. Sabbe, M.-F. Reyniers and K. Reuter, Catal. Sci. Technol., 2012 10.1039/c2cy20261a.
7. J. Hafner, J. Comput. Chem., 2008, 29, 2044–2078.
8. V. E. Van Doren, C. Van Alsenoy and P. Geerlings, Density functional theory and its application to materials : Antwerp, Belgium, 8–10 June 2000, American Institute of Physics, Melville, N.Y., 2001.
9. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, 1133–1135.
10. D. M. Ceperley and B. J. Alder, Phys. Rev. Lett., 1980, 45, 566–569.
11. P. Haas, F. Tran and P. Blaha, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 085104.
12. J. L. F. Da Silva, C. Stampfl and M. Scheffler, Phys. Rev. Lett., 2003, 90, 066104.
13. G. R. Darling and S. Holloway, Rep. Prog. Phys., 1995, 58, 1595–1672.
14. D. C. Langreth and M. J. Mehl, Phys. Rev. Lett., 1981, 47, 446–450.
15. D. C. Langreth and M. J. Mehl, Phys. Rev. B: Condens. Matter Mater. Phys., 1983, 28, 1809–1834.
16. B. Hammer, L. B. Hansen and J. K. Nørskov, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 7413–7421.
17. K. Honkala, A. Hellman, I. N. Remediakis, A. Logadottir, A. Carlsson, S. Dahl, C. H. Christensen and J. K. Nørskov, Science, 2005, 307, 555–558.
18. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh and C. Fiolhais, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 6671–6687.
19. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
20. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. L. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
21. R. Armiento and A. E. Mattsson, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 72, 085108.
22. Z. G. Wu and R. E. Cohen, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 235116.
23. R. W. Godby, M. Schluter and L. J. Sham, Phys. Rev. Lett., 1986, 56, 2415–2418.
24. P. Rinke, A. Qteish, J. Neugebauer, C. Freysoldt and M. Scheffler, New J. Phys., 2005, 7, 126.
25. A. K. Kelkkanen, B. I. Lundqvist and J. K. Nørskov, J. Chem. Phys., 2009, 131, 046102.
26. J. L. F. Da Silva, M. V. Ganduglia-Pirovano, J. Sauer, V. Bayer and G. Kresse, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 75, 045121.
27. J. P. Perdew, S. Kurth, A. Zupan and P. Blaha, Phys. Rev. Lett., 1999, 82, 2544–2547.
28. E. G. Moroni, G. Kresse, J. Hafner and J. Furthmuller, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 56, 15629–15646.
29. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
30. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627.
31. J. Paier, R. Hirschl, M. Marsman and G. Kresse, J. Chem. Phys., 2005, 122, 234102.
32. J. Paier, M. Marsman and G. Kresse, J. Chem. Phys., 2007, 127, 024103.
33. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
34. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207–8215.
35. C. Franchini, V. Bayer, R. Podloucky, J. Paier and G. Kresse, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 72, 045132.
36. C. Franchini, R. Podloucky, J. Paier, M. Marsman and G. Kresse, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 75, 195128.
37. A. Kiejna, G. Kresse, J. Rogal, A. De Sarkar, K. Reuter and M. Scheffler, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 035404.
38. A. Stroppa and G. Kresse, New J. Phys., 2008, 10, 063020.
39. J. Harl and G. Kresse, Phys. Rev. Lett., 2009, 103, 056401.
40. M. Dion, H. Rydberg, E. Schroder, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett., 2004, 92, 246401.
41. K. Lee, E. D. Murray, L. Z. Kong, B. I. Lundqvist and D. C. Langreth, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 081101.
42. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
43. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
44. F. London, Trans. Faraday Soc., 1937, 33, 8b–26.
45. T. Bučko, J. Hafner, S. Lebegue and J. G. Angyan, J. Phys. Chem. A, 2010, 114, 11814–11824.
46. V. G. Ruiz, W. Liu, E. Zojer, M. Scheffler and A. Tkatchenko, Phys. Rev. Lett., 2012, 108, 146103.
47. G. X. Zhang, A. Tkatchenko, J. Paier, H. Appel and M. Scheffler, Phys. Rev. Lett., 2011, 107, 245501.
48. P. Błoński and N. Lopez, J. Phys. Chem. C, 2012, 116, 15484–15492.
49. J. Hafner, J. Phys.: Condens. Matter, 2010, 22, 384205.
50. Z. Łodziana, G. Stoica and J. Pérez-Ramírez, Inorg. Chem., 2011, 50, 2590–2598.
51. R. Grau-Crespo and U. V. Waghmare, Molecular Modeling for the Design of Novel Performance Chemicals and Materials, CRC Press, 2012.
52. A. Zunger, S. H. Wei, L. G. Ferreira and J. E. Bernard, Phys. Rev. Lett., 1990, 65, 353.
53. R. Grau-Crespo, S. Hamad, C. R. A. Catlow and N. H. de Leeuw, J. Phys.: Condens. Matter, 2007, 19, 256201.
54. R. Grau-Crespo, N. H. de Leeuw and C. R. A. Catlow, Chem. Mater., 2004, 16, 1954–1960.
55. J. Carrasco, N. Lopez and F. Illas, Phys. Rev. Lett., 2004, 93, 225502.
56. R. Grau-Crespo, A. Y. Al-Baitai, I. Saadoune and N. H. De Leeuw, J. Phys.: Condens. Matter, 2010, 22, 255401.
57. R. Grau-Crespo, N. H. De Leeuw, S. Hamad and U. V. Waghmare, Proc. R. Soc. A, 2011, 467, 1925–1938.
58. W. G. Wulff, Z. Kryst. Miner., 1901, 34, 449–530.
59. I. N. Remediakis and et. al., personal communication.
60. S. Dahl, A. Logadottir, R. C. Egeberg, J. H. Larsen, I. Chorkendorff, E. Tornqvist and J. K. Nørskov, Phys. Rev. Lett., 1999, 83, 1814–1817.
61. N. Lopez, T. V. W. Janssens, B. S. Clausen, Y. Xu, M. Mavrikakis, T. Bligaard and J. K. Nørskov, J. Catal., 2004, 223, 232–235.
62. N. Lopez and C. Vargas-Fuentes, Chem. Commun., 2012, 48, 1379–1391.
63. R. T. Vang, K. Honkala, S. Dahl, E. K. Vestergaard, J. Schnadt, E. Laegsgaard, B. S. Clausen, J. K. Nørskov and F. Besenbacher, Nat. Mater., 2005, 4, 160–162.
64. J. Kleis, J. Greeley, N. A. Romero, V. A. Morozov, H. Falsig, A. H. Larsen, J. Lu, J. J. Mortensen, M. Dulak, K. S. Thygesen, J. K. Nørskov and K. W. Jacobsen, Catal. Lett., 2011, 141, 1067–1071.
65. U. Heiz, A. Sanchez, S. Abbet and W. D. Schneider, J. Am. Chem. Soc., 1999, 121, 3214–3217.
66. W. J. Stark, S. E. Pratsinis and A. Baiker, Chimia, 2002, 56, 485–489.
67. M. Mavrikakis, P. Stoltze and J. K. Nørskov, Catal. Lett., 2000, 64, 101–106.
68. B. Yoon, H. Hakkinen, U. Landman, A. S. Worz, J. M. Antonietti, S. Abbet, K. Judai and U. Heiz, Science, 2005, 307, 403–407.
69. H.-G. Boyen, G. Kastle, F. Weigl, B. Koslowski, C. Dietrich, P. Ziemann, J. P. Spatz, S. Riethmüller, C. Hartmann, M. Moller, G. Schmid, M. G. Garnier and P. Oelhafen, Science, 2002, 297, 1533–1536.
70. I. N. Remediakis, N. Lopez and J. K. Nørskov, Angew. Chem., Int. Ed., 2005, 44, 1824–1826.
71. Z. Łodziana, N. Y. Topsoe and J. K. Nørskov, Nat. Mater., 2004, 3, 289–293.
72. M. Digne, P. Sautet, P. Raybaud, P. Euzen and H. Toulhoat, J. Catal., 2002, 211, 1–5.
73. J. Venables, Introduction to Surface and Thin Film Processes, Cambridge University Press, Cambridge, 2000.
74. D. Teschner, R. Farra, L. D. Yao, R. Schlogl, H. Soerijanto, R. Schomacker, T. Schmidt, L. Szentmiklosi, A. P. Amrute, C. Mondelli, J. Pérez-Ramírez, G. Novell-Leruth and N. Lopez, J. Catal., 2012, 285, 273–284.
75. M. Calatayud and C. Minot, J. Phys. Chem. B, 2004, 108, 15679–15685.
76. M. V. Ganduglia-Pirovano, C. Popa, J. Sauer, H. Abbott, A. Uhl, M. Baron, D. Stacchiola, O. Bondarchuk, S. Shaikhutdinov and H.-J. Freund, J. Am. Chem. Soc., 2010, 132, 2345–2349.
77. P. L. Hansen, J. B. Wagner, S. Helveg, J. R. Rostrup-Nielsen, B. S. Clausen and H. Topsoe, Science, 2002, 295, 2053–2055.
78. S. J. Tauster, S. C. Fung, R. T. K. Baker and J. A. Horsley, Science, 1981, 211, 1121–1125.
79. D. Teschner, J. Borsodi, A. Wootsch, Z. Revay, M. Havecker, A. Knop-Gericke, S. D. Jackson and R. Schlogl, Science, 2008, 320, 86–89.
80. M. Schmid, A. Reicho, A. Stierle, I. Costina, J. Klikovits, P. Kostelnik, O. Dubay, G. Kresse, J. Gustafson, E. Lundgren, J. N. Andersen, H. Dosch and P. Varga, Phys. Rev. Lett., 2006, 96, 146102.
81. S. Helveg, C. Lopez-Cartes, J. Sehested, P. L. Hansen, B. S. Clausen, J. R. Rostrup-Nielsen, F. Abild-Pedersen and J. K. Nørskov, Nature, 2004, 427, 426–429.
82. M. Garcia-Mota, B. Bridier, J. Pérez-Ramírez and N. Lopez, J. Catal., 2010, 273, 92–102.
83. X. G. Wang, A. Chaka and M. Scheffler, Phys. Rev. Lett., 2000, 84, 3650–3653.
84. U. Martinez, L. B. Vilhelmsen, H. H. Kristoffersen, J. Stausholm-Moller and B. Hammer, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 205434.
85. N. Lopez, J. Gomez-Segura, R. P. Marin and J. Pérez-Ramírez, J. Catal., 2008, 255, 29–39.
86. G. Novell-Leruth, G. Carchini and N. Lopez, submitted.
87. J. J. Mortensen, B. Hammer and J. K. Nørskov, Phys. Rev. Lett., 1998, 80, 4333–4336.
88. J. Gomez-Diaz, K. Honkala and N. Lopez, Surf. Sci., 2010, 604, 1552–1557.
89. M. Garcia-Mota, J. Gomez-Diaz, G. Novell-Leruth, C. Vargas-Fuentes, L. Bellarosa, B. Bridier, J. Pérez-Ramírez and N. Lopez, Theor. Chem. Acc., 2011, 128, 663–673.
90. G. Henkelman and H. Jonsson, J. Chem. Phys., 2000, 113, 9978–9985.
91. P. Maragakis, S. A. Andreev, Y. Brumer, D. R. Reichman and E. Kaxiras, J. Chem. Phys., 2002, 117, 4651–4658.
92. G. Henkelman and H. Jonsson, J. Chem. Phys., 1999, 111, 7010–7022.
93. H. Falsig, B. Hvolbaek, I. S. Kristensen, T. Jiang, T. Bligaard, C. H. Christensen and J. K. Nørskov, Angew. Chem., Int. Ed., 2008, 47, 4835–4839.
94. D. Mei, P. A. Sheth, M. Neurock and C. M. Smith, J. Catal., 2006, 242, 1–15.
95. J. Andersin, N. Lopez and K. Honkala, J. Phys. Chem. C, 2009, 113, 8278–8286.
96. L. Nykanen, J. Andersin and K. Honkala, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 075417.
97. L. Andrussow, Angew. Chem., 1935, 48, 593–595.
98. J. Gomez-Diaz and N. Lopez, J. Phys. Chem. C, 2011, 115, 5667–5674.
99. I. Chorkendorff and J. W. Niemantsverdriet, Concepts of modern catalysis and kinetics, Wiley-VCH, Weinheim, 2003.
100. F. Abild-Pedersen, J. Greeley, F. Studt, J. Rossmeisl, T. R. Munter, P. G. Moses, E. Skulason, T. Bligaard and J. K. Nørskov, Phys. Rev. Lett., 2007, 99, 016105.
101. A. Vojvodic, F. Calle-Vallejo, W. Guo, S. Wang, A. Toftelund, F. Studt, J. I. Martinez, J. Shen, I. C. Man, J. Rossmeisl, T. Bligaard, J. K. Nørskov and F. Abild-Pedersen, J. Chem. Phys., 2011, 134, 244509.
102. J. N. Brønsted, Chem. Rev., 1928, 5, 231–338.
103. M. G. Evans and M. Polanyi, Trans. Faraday Soc., 1938, 34, 11–24.
104. T. Bligaard, J. K. Nørskov, S. Dahl, J. Matthiesen, C. H. Christensen and J. Sehested, J. Catal., 2004, 224, 206–217.
105. J. K. Nørskov, T. Bligaard, A. Logadottir, S. Bahn, L. B. Hansen, M. Bollinger, H. Bengaard, B. Hammer, Z. Sljivancanin, M. Mavrikakis, Y. Xu, S. Dahl and C. J. H. Jacobsen, J. Catal., 2002, 209, 275–278.
106. P. W. Atkins and J. de Paula, Atkins' Physical Chemistry, Oxford University Press, London, 8th edn, 2006.
107. K. J. Laidler, Chemical kinetics, Pearson Education, 2004.
108. K. Reuter and M. Scheffler, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 045433.
109. D. Teschner, G. Novell-Leruth, R. Farra, A. Knop-Gericke, R. Schlögl, L. Szentmiklósi, M. González Hevia, H. Soerijanto, R. Schomäcker, J. Pérez-Ramírez and N. López, Nat. Chem., 2012 DOI:10.1038/nchem.1411.
110. B. Bridier, D. Karhanek, J. Pérez-Ramírez and N. Lopez, ChemCatChem, 2012 DOI:10.1002/cctc.201200021.
111. A. Toftelund, PhD Thesis, DTU, 2012.
112. B. Liu and J. Greeley, J. Phys. Chem. C, 2011, 115, 19702–19709.
113. M. Salciccioli and D. G. Vlachos, J. Phys. Chem. A, 2012, 116, 4621–4628.
114. D. Torres, N. Lopez, F. Illas and R. M. Lambert, Angew. Chem., Int. Ed., 2007, 46, 2055–2058.
115. W. M. C. Sameera and F. Maseras, Comput. Mol. Sci., 2012, 2, 375–385.
116. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999–3093.
117. J. A. McCammon and M. Karplus, Annu. Rev. Phys. Chem., 1980, 31, 29–45.
118. A. Warshel, S. T. Russell and A. K. Churg, Proc. Natl. Acad. Sci. U. S. A., 1984, 81, 4785–4789.
119. C. F. Wong and J. A. McCammon, J. Am. Chem. Soc., 1986, 108, 3830–3832.
120. W. L. Jorgensen and C. Ravimohan, J. Chem. Phys., 1985, 83, 3050–3054.
121. W. L. Jorgensen and J. K. Buckner, J. Phys. Chem., 1987, 91, 6083–6085.
122. D. L. Beveridge and F. M. Dicapua, Annu. Rev. Biophys. Biophys. Chem., 1989, 18, 431–492.
123. P. A. Bash, U. C. Singh, F. K. Brown, R. Langridge and P. A. Kollman, Science, 1987, 235, 574–576.
124. A. Warshel and M. Levitt, J. Mol. Biol., 1976, 103, 227–249.
125. M. V. Bobetic and J. A. Barker, Phys. Rev. B: Condens. Matter Mater. Phys., 1970, 2, 4169–4175.
126. J. A. Barker, R. A. Fisher and R. O. Watts, Mol. Phys., 1971, 21, 657–673.
127. G. C. Maitland and E. B. Smith, Mol. Phys., 1971, 22, 861–868.
128. O. Tapia, J. Math. Chem., 1992, 10, 139–181.
129. L. Onsager, J. Am. Chem. Soc., 1936, 58, 1486–1493.
130. J. G. Kirkwood, J. Chem. Phys., 1934, 2, 351–361.
131. M. Cossi and V. Barone, J. Chem. Phys., 1998, 109, 6246–6254.
132. A. Klamt and G. Schuurmann, J. Chem. Soc., Perkin Trans., 1993, 2, 799–805.
133. E. Cances, B. Mennucci and J. Tomasi, J. Chem. Phys., 1997, 107, 3032–3041.
134. B. Mennucci, E. Cances and J. Tomasi, J. Phys. Chem. B, 1997, 101, 10506–10517.
135. E. Cances and B. Mennucci, J. Math. Chem., 1998, 23, 309–326.
136. L. Bellarosa, S. Calero and N. Lopez, Phys. Chem. Chem. Phys., 2012, 14, 7240–7245.
137. L. Bellarosa, J. M. Castillo-Sánchez, T. J. H. Vlugt, S. Calero and N. Lopez, Chem.–Eur. J., 2012 , in press.
138. X.-Q. Zhang, T. T. Trinh, R. A. van Santen and A. P. J. Jansen, J. Am. Chem. Soc., 2011, 133, 6613–6625.
139. J. K. Nørskov, J. Rossmeisl, A. Logadottir, L. Lindqvist, J. R. Kitchin, T. Bligaard and H. Jonsson, J. Phys. Chem. B, 2004, 108, 17886–17892.
140. X.-Q. Gong, A. Selloni, M. Batzill and U. Diebold, Nat. Mater., 2006, 5, 665–670.
141. J. Carrasco, A. Michaelides, M. Forster, S. Haq, R. Raval and A. Hodgson, Nat. Mater., 2009, 8, 427–431.
142. J. Andersin, P. Parkkinen and K. Honkala, J. Catal., 2012, 290, 118–125.
143. B. N. Zope, D. D. Hibbitts, M. Neurock and R. J. Davis, Science, 2010, 330, 74–78.
144. P. J. Feibelman, Science, 2002, 295, 99–102.
145. J. Carrasco, B. Santra, J. Klimes and A. Michaelides, Phys. Rev. Lett., 2011, 106, 026101.
146. K. Tonigold and A. Gross, J. Comput. Chem., 2012, 33, 695–701.
147. J. S. Hummelshoj, F. Abild-Pedersen, F. Studt, T. Bligaard and J. K. Nørskov, Angew. Chem., Int. Ed., 2012, 51, 272–274.
148. http://ceder.mit.edu/ .
149. D. Landis, J. Hummelshoj, S. Nestorov, J. Greeley, M. Dulak, T. Bligaard, J. K. Nørskov and K. Jacobsen, Comput. Sci. Eng., 2012, 1.
150. A. Elsener, C. C. M. Samson, M. P. Braendle, P. Buehlmann and H. P. Luethi, Chimia, 2007, 61, 165–168.
151. J. J. Iglesias, C. Bo, F. Maseras and N. Lopez, www.scipio.iciq.es.
152. VASP, www.vasp.at.
153. ADF, www.scm.com.
154. Gaussian, www.gaussian.com.
155. O. Lopez-Acevedo, K. A. Kacprzak, J. Akola and H. Hakkinen, Nat. Chem., 2010, 2, 329–334.
156. J. A. Lopez-Sanchez, N. Dimitratos, C. Hammond, G. L. Brett, L. Kesavan, S. White, P. Miedziak, R. Tiruvalam, R. L. Jenkins, A. F. Carley, D. Knight, C. J. Kiely and G. J. Hutchings, Nat. Chem., 2011, 3, 551–556.
157. J. Pérez-Ramírez, Nat. Chem., 2012, 4, 250–251.
158. P. Ball, Nature, 2006, 442, 500–502.
159. M. Haruta, Nature, 2005, 437, 1098–1099.
160. N. M. Esfandiari and S. A. Blum, J. Am. Chem. Soc., 2011, 133, 18145–18147.
161. Y. Hou, B. L. Abrams, P. C. K. Vesborg, M. E. Björketun, K. Herbst, L. Bech, A. M. Setti, C. D. Damsgaard, T. Pedersen, O. Hansen, J. Rossmeisl, S. Dahl, J. K. Nørskov and I. Chorkendorff, Nat. Mater., 2011, 10, 434–438.
162. M. Garcia-Mota, N. Cabello, F. Maseras, A. M. Echavarren, J. Pérez-Ramírez and N. Lopez, ChemPhysChem, 2008, 9, 1624–1629.
163. I. Serrano, M. I. Lopez, I. Ferrer, A. Poater, T. Parella, X. Fontrodona, M. Sola, A. Llobet, M. Rodriguez and I. Romero, Inorg. Chem., 2011, 50, 6044–6054.
164. R. Neumann and M. Dahan, Nature, 1997, 388, 353–355.
165. N. Lopez and G. Novell-Leruth, Phys. Chem. Chem. Phys., 2010, 12, 12217–12222.

---

Footnote

† This article is dedicated to the memory of Dr Jaime Gómez-Díaz who did his PhD in our group (2007–2011).

---