Unexpectedly large impact of van der Waals interactions on the description of heterogeneous catalyzed reactions: the water gas shift reaction on Cu(321) as a case example

. The molecular mechanisms of the water gas shift reaction on Cu(321) have been chosen to investigate the effect of dispersion terms on the description of the energy profile and reaction rates. The present study based on periodic DFT calculations shows that including dispersion terms does not change the qualitative picture of the overall reaction, maintaining the rate determining step and the predominant route. However, the effect of dispersion is different for different adsorbates  reactants, intermediates or products  with a clear net effect and with no compensation of errors. Thus, in the OH+OH H 2 O+O process the dispersion effects imply up to three orders of magnitude in the calculated reaction rates; the formation of carboxyl is highly disfavoured when dispersion terms are explicitly included and finally, the reaction rate for CO 2 production (at 463 K) through cis COOH dissociation is enhanced by three orders of magnitude by including dispersion terms in the calculation of the energy barrier. Consequently, the inclusion of dispersion terms largely affects the overall potential energy profile and produces tremendous changes in the predicted reaction rates. Therefore, dispersion terms must be included when aiming at obtaining information from macroscopic simulations employing for instance microkinetic or kinetic Monte Carlo approaches, where these effects should be

Please note that technical editing may introduce minor changes to the text and/or graphics, which may alter content. The journal's standard Terms & Conditions and the Ethical guidelines still apply. In no event shall the Royal Society of Chemistry be held responsible for any errors or omissions in this Accepted Manuscript or any consequences arising from the use of any information it contains.
Accepted Manuscript www.rsc.org/pccp PCCP Introduction Density functional theory (DFT) based calculations carried out on suitable periodic surface models have enormously contributed to our understanding of heterogeneously catalyzed reactions at the molecular level to the point that they have nowadays become a rather standard tool as illustrated in recently published books. 1,2 This type of computational methodology allowed to take into account environmental effects into the equilibrium structure of surfaces exposed to gases, 3 determining rather accurate energy profiles for many heterogeneously catalyzed reactions thus unveiling the molecular mechanism behind complex processes involving many elementary steps 4 and helped to derive useful concepts as descriptors allowing for a rational design of potential new and improved catalysts. 5,6 The information extracted from the DFT based calculations often includes transition state theory (TST) reaction rate constants for the elementary steps which can be used in subsequent macroscopic simulations of complex reactions. For instance the microkinetic modeling 7 of the water gas shift reaction (WGSR) catalyzed by Cu(111) by Gokhale et al. 8 and the kinetic Monte Carlo simulations by Yang et al. 9 and Prats et al., 10 both based on DFT calculated rates, constitute an excellent example of the interpretative and predictive power of this computational approach. Moreover, the increasing use of models involving stepped surfaces [11][12][13][14][15][16] or large metallic nanoparticles 17 provides more realistic models of the catalytic active sites. Yet, one of the remaining problems in this field concerns the accuracy of the calculated total energy defining the potential energy surface. In fact, commonly used Generalized Gradient Approach 18, 19 (GGA) forms of the exchange-correlation potential such as PW91 20 or PBE 21 provide a balanced and rather accurate description of the bulk properties of the three series of transition metals whereas other broadly used functionals such as RPBE exhibit a poorer behavior and excessively stabilize surface energies. 22,23 Nevertheless, these GGA functionals do not provide accurate enough results for main group elements containing molecules 24 and, as already pointed out by Kristyan and Pulay twenty years ago, 25 neglect dispersion terms which may play a non-negligible role in the molecular picture of heterogeneously catalyzed reactions.
The first of the two shortcomings of DFT mentioned above has precisely triggered the development of new and more accurate functionals such as the widely used B3LYP hybrid functional 26 or the series of Minnesota hybrid meta-GGA functionals. [27][28][29] These Page 2 of 28 Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript are, however, seldom used in computational studies in heterogeneous catalysis due to the difficulty that these methods face when applied to metals. 30-,32 Because of their good performance, hybrid functionals are surely the most popular kind of functionals in molecular chemistry and homogeneous catalysis. 33 However, for extended systems they have large, and often excessive, computational demands as compared to GGA type functionals. This is due to the long range of the exchange interactions when making use of programs working in the real space and due to the requirement for dense Brillouin zone sampling when relying on programs using plane wave basis sets. 34 Nevertheless, it is often argued that, for chemical reactions taking place at metal surfaces, calculated relative energies are much less affected than absolute energies by the inherent errors of GGA type functionals 2,4,35,36 and this is surely one of the keys of the success of this type of calculations.
The effect of van der Waals (also known as dispersion) interactions on adsorption properties has been the focus of an intense research in the past few years, especially after the landmark contributions of Grimme and coworkers, [37][38][39] which has triggered many new theoretical developments and the appearance of a plethora of new functionals aiming to account for these terms in an accurate and non-empirical way as recently critically reviewed by Klimes and Michaelides. 40 Dispersion terms play an important role in chemical and physical processes involving biomolecules and their role in conformational related problems and in thermochemistry has been recently reviewed. 41 These terms largely affect the adsorption properties of molecules at surfaces and can even be the dominant term as in the case of aromatic molecules interaction with the basal plane of MoS 2 , 42 hydrocarbons interacting with zeolites 43 or graphene on metallic surfaces; 44 a review on the role of dispersion terms on adsorption properties has been recently published. 45 In spite of the large number of articles devoted to study the importance of dispersion terms in adsorbate-surface interactions, there is almost no information regarding the effect of dispersion terms in the energy profile of heterogeneously catalyzed reactions, especially for complex mechanisms involving several elementary steps. An important catalyzed reaction with special technological relevance 46 supported on different oxides 47 although other metals and supports have also been proposed. 48,49 The molecular mechanism for the low temperature catalyzed WGSR involves a rather large number of elementary steps and two possible routes, redox or carboxyl, are possible. 8 These have been studied in depth for the Cu(111) 8 and Cu(321) 13 surfaces; the latter one, containing different low-coordinated sites, offers a more realistic model of the catalyst. Moreover in the latter case, there is detailed information regarding the structure of the many transition state structures involved in the mechanism and thus constitutes an excellent system to check the effect of the dispersion terms on the overall energy profile and rate constants. This is precisely the goal of the present paper. We will provide compelling evidence that while the qualitative picture of the overall reaction scheme is not largely affected by the inclusion of the dispersion terms, there are significant differences in the calculated reaction rates, which have important implications in the macroscopic description of the overall process via microkinetic or kinetic Monte Carlo simulations.

Elementary steps in the water gas shift reaction
In this section we will briefly summarize the most salient features of the reaction mechanisms proposed for the WGSR. These can be grouped in two general mechanisms, The study of Fajín et al. 13 evidenced that the presence of steps increases the water adsorption energy and decreases the energy barrier of water dissociation and atomic hydrogen recombination steps which on Cu(321) are found to constitute the ratedetermining steps (rds). Interestingly, these two elementary steps are also the rds for the WGSR on Cu(111) but on the stepped Cu(321) they have similar energy barriers and reaction rates while on the flat Cu(111) surface the water dissociation has an energy barrier considerably larger than the hydrogen recombination.
In the present work, the effect of van der Waals interactions will be explicitly taken into account for all adsorption, reaction and desorption steps outlined above using Cu(321) as catalyst model as described in the next section.

Surface model and computational details
The interaction of the different reactants, intermediates and products involved in the WGSR catalysed by the Cu(321) surface has been obtained from periodic DFT calculations modelled through the usual repeated slab approach with a 2×2×1 supercell constructed using the optimum lattice parameter of 3.63 Å for the computational method chosen here and described in detail below; note that this is sufficiently close to the experimental value of 3.62 Å. 50 It is also worth pointing out that, in order to minimize lateral interactions, the unit cell for the Cu(321) slab model thus defined is larger than the one previously used by Fajín et al. 13 The 2×2×1 supercell used in the present work contains 60 Cu atoms distributed in four atomic layers as schematically shown in Figure 1 and consists of a monoclinic prism with an angle of 104.96º between the x and y axes and of 90º for the angles between x and z or y and z axes. Further, the unit cell vectors along the x, y and z directions have different lengths. The corresponding fractional coordinates of the atoms in this unit cell were obtained using the Materials Studio computer code (version 8.0). 51 The unit cell for the two-dimensional slab thus obtained was modified by adding a vacuum region of 12 Å and scaling the fractional coordinates conveniently so as to obtain a unit cell that can be replicated in three dimensions as required when using a plane-wave periodic DFT approach. The resulting slab was further modified by allowing full relaxation of the position of the uppermost 28 Cu atoms within the computational approach described below.
In order to investigate the impact of the dispersion terms in the calculated energy profile we compare results from two series of periodic DFT calculations, both carried out with the VASP code. [52][53][54] For the first series we rely on the PW91 calculation of Fajín et al. 13 whereas in the second one the effect of the van der Waals interactions has been included by adding the dispersion term obtained from the D2 method of Grimme 38 to the PBE calculated energy (PBE-D2). Note in passing by that, in spite of its semiempirical flavor, the D2 method has been shown to properly describe the physisorption and chemisorption states of graphene with Ni(111). 44 Nevertheless, to validate the present results some key calculations have been carried out with the D3 parameterization of Grimme 39 (PBE-D3) and with the more physically grounded method proposed by Tkatchenko et al. 55 (DFT-T). Note also that PW91 and PBE results for bulk properties of transition metals 22,23 and also for the description of the adsorption energy of WGSR species are very similar. 56 The valence electron density was expanded in a plane-wave basis set with a cut-off of 415 eV for the kinetic energy. The effect of core electrons in the valence electron density was taken into account through the projector augmented wave (PAW) method 57 as implemented in VASP. 58 Numerical integration in the reciprocal space was carried out by employing a 5×5×1 Monkhorst-Pack grid of special k-points. 59 The energy cut-off and k- where E slab-m and E slab refer to the total energy of the slab model representing the Cu (321) surface with and without the m adsorbate and E m corresponds to the total energy of the molecule in the gas phase computed, as usual, by placing it in a box with the same size of the unit cell for the slab. For the situations with two adsorbates above the surface unit cell the co-adsorption energy is calculated as: where E slab-m1-m2 stands for the total energy of the system formed by the two species adsorbed on the slab and E slab , E m1 and E m2 are as in Eq. (1).

Adsorption and co-adsorption of reactants, intermediates and products
In this subsection we discuss the effect of dispersion on the most favourable adsorption and co-adsorption configurations of the species involved in the WGSR mechanism catalysed by the Cu(321) surface. A summary of calculated results regarding adsorption energy is reported in Tables 1 and 2 whereas Figure 2 and Figure 3 report the equilibrium geometry of the adsorbed and co-adsorbed states for all involved species, respectively.

Physical Chemistry Chemical Physics Accepted Manuscript
In agreement with the previous study of Fajín et al., 13   Co-adsorption energies of the most stable configuration of reactant and product pairs for the different elementary steps in the WGSR mechanism are listed in Table 2 Table 2 also evidence that the effect of dispersion is different for different adsorbate pairs which, as we will show in the next sections, must have an influence on several energy barriers and on the resulting TST reaction rates. In fact, the effect is almost negligible (∼ 0.05 eV) for some cases such as OH a +H b , O+H, or H a +H b , it is intermediate (<0.15 eV) for some others such as OH a +OH b or CO a +O b and quite large (>0.25 eV) in a few cases such as CO a +OH b or CO 2 +H. The differences are large enough to be significant and likely to be present if other methods are used to estimate the dispersion contribution to the total energy.

Energy barriers of the elementary steps
Now we come to the most important part of the present work, namely the description of the calculated energy barriers for the different elementary steps in the WGSR on Cu(321). The energy barrier for each individual step of the reaction mechanism has been calculated as the energy difference between the transition state, located using the improved Dimer method 61 and that of the most stable adsorption (or co-adsorption) configuration for the reactant(s). For the transition state (TS) calculations, a first step involved the search of a first order saddle point with the slab structure fixed and, in a second step, the atoms in the  Table 3 where equivalent ZPE corrected PW91 values have been included for comparison. Schematic representations of the transition state geometries are given in Figure   4. For products desorption (i.e., H 2 → H 2(g) and CO 2 → CO 2(g) ) the TSs are assumed to be their final states, that is, H 2 and CO 2 in the gas phase. Thus, the energy barriers for these processes are equal to their adsorption energies in absolute value (i.e., 0.12 and 0.28 eV, respectively) and, consequently, are not included in Table 3. In the following we will discuss the effect of dispersion in the different steps by comparing to the results reported by OH dissociation. This step is endoergic by 0.48 eV, almost the same value than for Cu(111) surface, and has an energy barrier of 1.51 eV (Table 3), higher than for the flat surface (1.19 eV) 8 and, as the previous step, very similar to the value obtained without including vdW corrections (1.55 eV). 13 The high energy barrier for this reaction implies that the CO 2

Physical Chemistry Chemical Physics Accepted Manuscript
formation via the redox mechanism through direct OH dissociation in the Cu(321) surface will surely be the least frequent route among all. OH disproportionation. This is an alternative path for producing atomic oxygen, it is exoergic by -0.10 eV with an energy barrier of only 0.46 eV, becoming the main route of redox mechanism. Surprisingly, this value is increased by 0.32 eV when vdW interactions are not considered (Table 3).  (Table 2).
Carboxyl dehydrogenation. This is an exoergic step (-0.52 eV) with an energy barrier of 0.80 eV (Table 3), again significantly lower than for the flat Cu(111) surface (1.18 eV) and also lower that the PW91 value obtained (1.10 eV). 13 Since the inclusion of vdW interactions stabilizes both the transition state and the final products, CO 2 formation through carboxyl intermediate will probably play a more important role in the WGSR over Cu(321) than in Cu(111).
Carboxyl disproportionation by hydroxyl. This step involving cis-COOH is also exoergic by -0.64 eV, more than for the Cu(111) surface (-0.37 eV). The energy barrier for this process is 0.33 eV to be compared to 0.55 eV without vdW correction, which in turn needs to be compared with the result obtained for the flat surface (0.38 eV).

Page 11 of 28
Physical Chemistry Chemical Physics

Physical Chemistry Chemical Physics Accepted Manuscript
H recombination. This step is common to the four investigated routes. According to the present result, H atoms are not provided only by water dissociation but also by carboxyl dehydrogenation. Another possible process for H production is OH dissociation which is very unfavourable even in this stepped surface. Although H 2 does not adsorb molecularly on Cu(111) it does on the Cu(321) surface by -0.12 eV (see Table 1). This reaction is endoergic by 0.31 eV, with an energy barrier of 0.78 eV, again smaller than the ZPE corrected 0.96 eV value for the Cu(111) surface.

Reaction rates of the elementary steps
From the calculated zero point corrected energy barriers and vibrational frequencies one can readily obtain the corresponding transition state theory rates at the temperature of interest. Table 3 reports the calculated rates for the elementary steps at 463 K; this is the same temperature used in previous work regarding the WGSR mechanism on the Cu(111) 8 and Cu(321) 13 surfaces, where dispersion terms were not included in the calculations.
As discussed above and in agreement with Fajín et al., 13  The dissociation of adsorbed water defines the rds in all cases and all results seem to indicate that the associative mechanism is clearly preferred. Nevertheless, the rates predicted by the present PBE-D2 calculations for the rds are one order of magnitude smaller than the values reported from PW91 not including dispersion terms. Here it is worth to mention that the low energy barrier for the reverse process leads to a reaction rate of 6.38·10 6 s -1 , six orders of magnitude larger than in the flat surface (5.25·10 0 s -1 ) and greater than the rate for the forward process. Clearly, it is not possible to extract reliable

Physical Chemistry Chemical Physics Accepted Manuscript
conclusions from the energy barrier for the rate-determining step in the forward direction only. Table 3 show that the effect of dispersion on all elementary steps is very different. In some cases including or not these effects has a variation of barely one order of magnitude in the calculated rates for the forward reactions. This is for instance the case of H 2 O → OH+H, OH → O+H, CO+O →CO 2 , cis-COOH+OH → CO 2 + H 2 O, and H+H → H 2 steps. However, in some other steps the effect of dispersion implies up to three orders of magnitude in the calculated reaction rates, as this is the case for OH+OH → H 2 O+O, which now appears as the dominant source of adsorbed O. This is especially relevant since including dispersion affects the reaction rate of the main step in CO 2 production also by three orders of magnitude but in the opposite sense. The formation of carboxyl is highly disfavoured when dispersion terms are explicitly included. The reaction rate for CO 2 production (at 463 K) through cis-COOH dissociation is enhanced by three orders of magnitude by including dispersion terms in the calculation of the energy barrier. Finally, it is worth to mention that, not surprisingly, dispersion terms largely affect the reaction rates for adsorption and desorption steps (not reported).

Results in
At this point, one may still argue that semiempirical dispersion treatments based on atom pairwise potentials may be inadequate for metallic systems and question the overall validity of the present results. The selected adsorption energy values reported above calculated using the D3 and Tatchenko methods would indicate that this is not the case.
Nevertheless, to reach a firm conclusion it is convenient to inspect energy barriers as well.
To this end, a new series of calculations has been carried out using the method recently  Table 2 to the point that the changes are within the incertitude of DFT methods. The equilibrium geometries calculated with PBE-D2 and PBE-Anderson methods are almost the same and the calculated surface reaction energy barriers differ in less than 0.05 eV. Consequently, one can firmly claim that the conclusion arising from the PBE-D2 calculations reported in the present work are physically meaningful.

Conclusions
The effect of dispersion terms on the description of the energy profile and reaction rates of a complex heterogeneously catalysed process has been studied in detail taking the water gas shift reaction on Cu(321) as a case example. This is a convenient case study because a rather large number of elementary steps and because of the existence of previous results regarding the molecular mechanism of the overall reaction on this surface 13  One must admit that the present results have been obtained from a particular choice in the method used to estimate dispersion terms. Nevertheless, calculations for the key steps have been also carried out including dispersion with two alternative methods indicating that the conclusions of the present work are sound and not biased.
qualitative description of the WGSR catalysed by Cu(321), their presence largely affects the overall potential energy profile and produces tremendous changes in the predicted reaction rates. Consequently, dispersion terms must be included when aiming at obtaining information from macroscopic simulations employing for instance microkinetic or kinetic Monte Carlo approaches, where these effects should be clearly shown. with numbers referring to the labels in Figure 1 and