Influence of the van der Waals interaction in the dissociation dynamics of N 2 on W ( 110 ) from first principles

Using ab initio molecular dynamics (AIMD) calculations, we investigate the role of the van der Waals (vdW) interaction in the dissociative adsorption of N2 on W(110). Hitherto, existing classical dynamics calculations performed on six-dimensional potential energy surfaces based on density functional theory (DFT), and the semi-local PW91 and RPBE [Hammer et al. Phys. Rev. B 59, 7413 (1999)] exchange-correlation functionals were unable to fully describe the dependence of the initial sticking coefficient on the molecular beam incidence conditions as found in experiments. N2 dissociation on W(110) was shown to be very sensitive not only to short molecule-surface distances but also to large distances where the vdW interaction, not included in semilocal-DFT, should dominate. In this work, we perform a systematic study on the dissociative adsorption using a selection of existing non-local functionals that include the vdW interaction (vdW-functionals). Clearly, the inclusion of the non-local correlation term contributes in all cases to correct the unrealistic energy barriers that were identified in the RPBE at large molecule-surface distances. Among the tested vdW-functionals, the original vdW-DF by Dion et al. [Phys. Rev. Lett. 92, 246401 (2004)] and the ulterior vdW-DF2 give also an adequate description of the N2 adsorption energy and energy barrier at the transition state, i.e., of the properties that are commonly used to verify the quality of any exchange-correlation functional. However, the results of our AIMD calculations, which are performed at different incidence conditions and hence extensively probe the multi-configurational potential energy surface of the system, do not seem as satisfactory as the preliminary static analysis suggested. When comparing the obtained dissociation probabilities with existing experimental data, none of the used vdW-functionals seems to provide altogether an adequate description of the N2/W(110) interaction at short and large distances.


I. INTRODUCTION
2][3][4] H 2 diffraction experiments 5 or the possibility of quantum state preparation (detection) of the incident (scattered) molecules [6][7][8][9][10] allows one to explore in more detail the quantum nature of such interactions.Nonetheless, gas-surface experiments cannot monitor step by step the evolution of the molecule (atom) and the surface.2][13] In this respect, there are various examples showing that not only the reactivity [14][15][16] but also scattering properties such as the diffraction peaks positions, 17,18 angular distributions, [19][20][21][22][23] rotational states, [24][25][26][27][28] and energy loss 21,23,29 are reasonably described by this level of theory.Considering that in all these cases, the potential energy surfaces (PESs) have been calculated with the (hitherto standard) semilocal exchange-correlation (xc) functionals that do not correctly describe the long-range van der Waals (vdW) interaction, such an agreement is somehow unexpected at first thoughts.This is more the case when thinking about the scattering properties.As recently pointed out for O 2 scattered off Ag(111), 23 the adequacy of semilocal DFT stands in that the scattering is ultimately ruled by the regions where the rehybridization phenomena between the molecular and the surface electronic states, which are reasonably described by semilocal-DFT, occur.This finding can most likely be extended to many other systems.
Knowing that for scattering problems, the same behavior is naturally expected when treating the adsorption processes, i.e., the reactivity.Here the question that emerges is whether the vdW interaction is still negligible in those systems characterized by a transition state located at the entrance channel, where the molecule is typically far from the surface and the semilocal-DFT description might be insufficient.The dissociation of H 2 and D 2 on Ru(0001) is precisely among this group.It has recently been found that only functionals that incorporate the vdW interaction are able to reasonably well reproduce the width of the initial sticking probability curve as a function of incidence energy. 30The desorption of CO from Ru(0001) is another recent example pointing to the vdW interaction as a crucial ingredient in understanding this process. 31,32he discrepancy between theory and experiments found when comparing the sticking coefficient S 0 , i.e., the dissociation probability at zero coverage, for various incident conditions (angle Θ i and kinetic energy E i ) of N 2 on W(110) is another appealing case for which vdW may play a prominent role.The dissociation process in this system is particularly sensitive to the PES properties at large molecule-surface distances, i.e., in the region where vdW interactions can be relevant. 33,34Strictly speaking, N 2 dissociation is a barrierless process, despite the S 0 vs. E i curve at normal incidence to the surface (Θ i = 0 • ) shows the characteristic S-like profile of activated systems.At Θ i = 0 • and low E i , the dissociating molecules are first dynamically trapped around the molecular well located at about 2.6 Å from the surface before finding the pathway towards the configurations leading to dissociation.The experimental and theoretical S 0 values smaller than 3 × 10 −3 found under these incidence conditions were shown to be a consequence of the extremely narrow configurational space accessing the N 2 adsorption well. 33,34Direct dissociation through the dissociating paths, which are mainly located over bridge and hollow sites, is forbidden by the energy barriers that exist at 2-3 Å above the surface.As soon as E i increases and permits to overcome such barriers, the indirect dissociation through the N 2 well is progressively quenched in favor of the direct channel, which is found to dominate at off-normal incidence (Θ i ≥ 45 • ) 35 and at normal incidence and high E i . 34mportant to us, none of the classical molecular dynamics (MD) calculations performed on two different semilocal-DFT PESs, namely, (i) those in Refs.33 and 34 that use the PW91 xc-functional 36 and (ii) the ones in Ref. 37 that use the RPBE xc-functional, 38 fully agree with the experimental S 0 for all the incident angles and energies.On the one hand, MD simulations performed on the PW91-PES reproduce the data for all E i at 60 • incidence and for normal incidence and high E i , but fail for Θ i = 0 • and low E i .On the other hand, the results obtained with the RPBE-PES are reasonably good at Θ i = 0 • for all the experimental E i -range, but completely fail at Θ i = 60 • .By comparing the evolution of the trajectories between the two PESs, it was concluded that at large molecule-surface distances, the RPBE is being too repulsive, despite it seems to provide a better description of the molecule-surface interaction in regions close to the surface than the PW91. 37The same conclusion was later on obtained when comparing the rovibrational energy of the scattered molecules with available experimental data. 26otivated by all these previous results that pinpoint the significance of the N 2 -W(110) large distance region, our purpose here is to determine the relevance of the vdW interaction in the dissociative dynamics of N 2 on W(110).To this aim, we perform ab initio molecular dynamics (AIMD) calculations within the Born-Oppenheimer approximation (BOA) using a selection of the existing vdW-functionals.0][41][42] Among them, the vdW xc-functional developed by Dion et al. 43 and its subsequent offsprings [44][45][46][47][48][49] (henceforth denoted vdW-functionals) are particularly promising because this scheme naturally includes the vdW interaction within the DFT framework.Apart from the recent studies of Refs.30 and 50, the various vdW-functionals have been tested from a pure static point of view by looking at the binding energies of the S22 set 51 of selected complexes, 45,52 the adsorption energies of organic molecules on surfaces, 53,54 and of 2D layered materials, [55][56][57][58] as well as various properties of solids. 47These kind of studies, though essential, provide limited information, for instance, on the quality of the multidimensional PES that characterizes the interaction of molecules with surfaces.In this respect, gas-surface dynamics simulations are a powerful tool to explore this issue and, hence, to test the performance of any xc-functional in describing short-and long-range interactions.Previous gas-surface simulations of H 2 on Ru(0001), 30,59 H 2 on Cu(111), 15,16 and the above mentioned N 2 on W surfaces 37 have revealed the limitations of existing (non-empirical) semilocal xc-functionals in describing at the same foot regions of small and large density gradients.Therefore, an additional outcome of our present study is that the comparative analysis of the AIMD results will serve us to test the adequacy of different vdW-functionals in describing the multi-configurational space that rules the N 2 -W(110) interaction.
The paper is organized as follows.Section II starts with a brief overview of the MD calculations performed in Refs.33, 34, and 37 and continues by describing all the details of the AIMD calculations performed in the present study with some of the existing vdW-functionals.The static properties of the different PESs that will be relevant to understand the dissociation dynamics, as well as the results of the AIMD calculations themselves are discussed in Sec.III.The summary and conclusions of our work can be found in Sec.IV.

II. AIMD CALCULATIONS DETAILS
The scattering and dissociative adsorption of N 2 on the W(110) surface were studied previously 26,33,34,37 by performing MD calculations on a (previously generated) frozen-surface (FS) six-dimensional (6D) adiabatic PES that only depends on the position of the two N atoms over the surface, since the W atoms are fixed at their positions in the relaxed W(110) surface (see below).The continuous 6D PES was constructed using the so-called corrugation reducing procedure (CRP) 60 to interpolate a set of 5610 energy values that were calculated with DFT and the generalized gradient approximation for the xc energy functional.In particular, the PW91 was used in Refs.33 and 34 and the RPBE in Ref. 37. In the following, we will refer to these PESs as PW91-CRP and RPBE-CRP, respectively.We will restrict the use of the term MD to differentiate the simulations performed on top of these precalculated PESs from the AIMD calculations done in this work.The classical trajectory is calculated in MD by solving the Hamilton-Jacobi equations of N 2 with a Bulirsch-Stoer integrator. 61For each incidence condition defined by Θ i and E i , the (dissociative) sticking coefficient S 0 was evaluated from the outcome of 5000 trajectories that sample with a Monte-Carlo procedure the azimuthal incidence angle of the beam, the N 2 orientation respect to the surface and its position over the surface unit cell.In Refs.33 and 37, the results correspond to pure classical dynamics calculations that neglect the N 2 zero point energy (ZPE).The latter was included in Ref. 34 as one of the initial conditions of the impinging N 2 (quasiclassical calculation).Note, however, that during the trajectory calculation, the evolution of this degree-of-freedom is treated classically and not quantically.The classical and quasiclassical sticking probabilities are almost identical, the latter being slightly higher than the former.
In the present work, the use of AIMD allows us to simulate a statistically accurate number of trajectories for five vdWfunctionals, which is computationally less demanding than constructing the same number of six-dimensional CRP PESs.Moreover, this strategy prevents from interpolation errors.
AIMD calculations are performed within the BOA with the  code 62 using in all cases the same supercell that was applied in Refs.34 and 37 to calculate the DFT energy grid used in the interpolation of the 6D PW91-and RPBE-CRP PESs.Briefly, the W(110) surface is modeled by means of a periodic five-layer slab with a (2 × 2) surface unit cell and a supercell vector along the surface normal of 12d, being d = 2.24 Å the theoretical interlayer distance.After surface relaxation, the periodic slabs are separated by 18.16 Å of vacuum.We have verified that, except in one case that will be discussed in Sec.III A, the lattice constants calculated with the rest of functionals differ in less than 1% respect to the PW91 value of 3.17 Å (to be compared with the nominal 3.16 Å).Thus, the use of the PW91 supercell is well justified.The electron-core interaction is treated with the plane augmented wave (PAW) approximation 63,64 in all the calculations using the vdW-functionals 65 and with the ultra-soft pseudopotentials 66 in the PW91 and RPBE calculations, as done in Refs.34 and 37, in order to quantify the interpolation accuracy with the present AIMD calculations.The electron wave functions are expanded in a plane-wave basis set with an energy cut-off of 560 eV.The fractional occupancies are determined through the broadening approach of Methfessel and Paxton with N = 1 and σ = 0.4 eV. 67The Brillouin zone integration is performed with a 4 × 4 × 1 Monkhorst-Pack grid of special k-points. 68The energy criterion for total energy self-consistency is 10 −5 eV.
The vdW interaction is included in the AIMD calculations following the vdW density functional method 43,46 that is fully based on the electron density.Aside the original vdW-DF functional by Dion et al. 43 and the subsequent improved vdW-DF2 version, 46 we use the optimized vdW functionals developed by Klimeš et al. that are also based on the original vdW-DF scheme, but use alternative less-repulsive exchange functionals. 45,47Following the notation of these authors, the optimized functionals will be denoted as optPBE-vdW, optB88-vdW, and optB86b-vdW in the following.The performance of all these functionals in describing different static properties of dimers, solids, adsorbates, etc., has been extensively discussed in the literature.Here, we briefly summarize those results of interest for our present study.The vdW-DF and still the vdW-DF2 usually overestimate the lattice constant of solids and the bond length of weakly interacting compounds as a consequence of the rather repulsive character of the original revPBE exchange interaction.In contrast, the optimized-vdW functionals provide by construction a much more accurate description of the S22 data reference set of hydrogen-and dispersion-bonded dimers.In particular, the optB88 gives the best overall description. 45When performing standard tests on solid state properties such as lattice constant, bulk moduli, and atomization energies, the three optimized vdW improve over the vdW-DF and vdW-DF2 and among them, the optB86b-vdW gives the smallest errors. 47The vdW-DF and the optPBE-vdW have, respectively, an exchange term very similar to that of RPBE and PW91.Therefore, the comparison between the results obtained with vdW-DF and RPBE or between optPBE-vdW and PW91 can provide information on the effect of the correlation correction introduced in the vdW-functionals.
The vdW-DF implemented in  is based on the Román-Pérez and Soler algorithm 69 that considerably reduces the computational effort associated to the self-consistent evaluation of the vdW energy.This property is particularly important when dealing with AIMD calculations as done here.In our case, we have verified that the computational time consumed in one ionic step is a factor of 1.5 slower with the vdW functionals than with the semilocal PW91 and RPBE.
In all AIMD calculations, we adopt the FS approximation as done in Refs.33, 34, and 37.In particular, the W surface atoms are kept fixed and equal to the relaxed positions obtained with the PW91 xc-functional.Here, we only perform pure classical trajectory calculations that neglect the ZPE of the impinging N 2 .Thus, for each incidence energy E i and polar angle Θ i , the incidence azimuthal angle and the positions of the two N atoms over the surface are randomly mapped using a conventional Monte-Carlo sampling that in practice is identical to the one used in Refs.33 and 37, where the ZPE was also neglected.All the trajectories start with the N 2 center of mass at 6.5 Å from the surface topmost layer.This distance assures negligible initial forces along the surface normal with values below 0.04 eV/Å even for the far-ranged vdW-functionals.The time steps used in the velocity-Verlet integration algorithm are 0.2 fs in the AIMD calculations performed with the optB86b-vdW and 0.5 fs in the rest.In most cases, these time steps are small enough to assure total energy conservation along the trajectory integration.However, we observed a few trajectories, in particular when using the optB88-and mostly the optB86b-vdW functionals, in which the energy conservation was rather poor.Noticeably, we verified that the values of the dissociation probabilities can extensively depend on the criteria used for energy conservation.In this respect, our convention has been to discard those trajectories with a total energy change larger than 5% of the initial incidence energy E i .All in all, the dissociation probabilities are calculated from the outcome of 100-200 trajectories.We remark that despite this number is much smaller than the one used in the MD calculations performed with any of the CRP PESs, we have verified that the same set of trajectories (i.e., same sampling of the initial conditions) used in our AIMD calculations already reproduces in the CRP PESs molecular dynamics calculations the dissociation probabilities of Refs.33 and 37.
The AIMD trajectories are counted as dissociated when the N-N internuclear distance is r ≥ 2.4 Å and as reflected when the center of mass moves away from the surface and arrives at a distance Z > 4.5 Å from the surface.In case of the vdW AIMD calculations, we have verified that at these distances, the molecules have enough energy to surmount the remaining vdW attraction.

A. Static properties
The adsorption energy is usually one of the important static properties that serves to corroborate the quality of the xcfunctionals.Thus, we start by calculating the depth of the N 2 adsorption well with the different functionals and comparing them with the available experimental data that estimate a value of ∼450 meV. 70In this adsorption configuration, the molecule stays upright on top of a W atom.The depth of the adsorption well is calculated as where E[N 2 /W(110)] is the energy computed with N 2 in the adsorption configuration and allowing relaxation of the two topmost W layers, E(N 2 ) is the energy of the isolated molecule at its equilibrium bond length, and E[W(110)] that of the isolated five-layer slab that is relaxed along the surface normal for each xc-functional.The energy and force tolerances used in all relaxations and in the nudged elastic band (NEB) 71 calculations below are 10 −6 eV and 0.02 eV/Å, respectively.The E a values are displayed in Table I together with the distance Z of the N 2 center of mass to the W atom beneath.In all cases, we obtain a similar N-N internuclear distance of 1.13 Å that represents a slight elongation with respect to the theoretical bond length of 1.11 Å found for N 2 in gas-phase (vacuum).
For the sake of comparison, we also show in brackets the values obtained with the PW91 frozen-surface supercell that is used in all AIMD calculations.In general, we observe that the adsorption energies in brackets are about 100 meV smaller.
Such differences are more a consequence of E[N 2 /W(110)] being calculated without allowing for surface relaxation than of the slightly different z-relaxation of W(110) found for each functional with respect to the relaxed PW91 W(110) surface.Table I reveals that even if the adsorption configuration itself is only affected by small changes in the N 2 -surface distance, there are significant differences in the E a values among the different functionals.Compared with the experimental data, the RPBE, vdW-DF2, and vdW-DF seem to provide following this order a better description of the adsorption energy, while the optimized vdW-functionals (opt) would be largely overestimating the molecule-surface attractive interaction, with optB86b-vdW being the one that shows the largest discrepancy.The trend found in the adsorption energies coincides with the results found with these vdW-functionals on other systems.However, the tendency of the original vdW-DF(2) to underbinding gives, paradoxically, a better agreement with the experimental data in this system.
Following with the static analysis of the different functionals, we next focus on the PESs properties that were identified as relevant to the dissociation dynamics. 34,37As mentioned in the Introduction, at normal incidence and E i ≤ 0.5 eV, dissociation is dominated by the indirect channel, in which the dissociating molecules are first dynamically trapped in the vicinity of the molecular adsorption well.There the molecules spend few ps exploring the configurational space, until they find the minimum energy path to dissociation.At 60 • polar incidence as well as for normal incidence and E i > 0.5 eV, the indirect channel is no longer at work and dissociation occurs for those molecules that start their approach to the surface in configurations that permit dissociation after overcoming the barriers of ∼0.5 eV (with PW91, but larger with RPBE) that exist above 2 Å from the surface (direct channel).
Focusing on the indirect channel, its efficiency mostly depends on two factors: (i) the accessibility to the N 2 adsorption well and (ii) the accessibility to the minimum energy dissociation path (MEDP) from the well.Regarding point (i), the two dimensional (r, Z) cuts of the PES for each xc-functional plotted in Fig. 1 show that the vdW-functionals predict the absence of energy barriers to access the adsorption well.As shown in Ref. 37, this was also the case for the PW91, but not for the RPBE xc-functional that exhibits an energy barrier (EB) of ∼0.08 eV in line with its repulsive character above 2.5-3.0Å.However, this shallow barrier is not the only, TABLE I. Adsorption energies E a of N 2 on W(110) obtained with various xc-functionals (the zero-point energy is not considered, see Eq. ( 1)).In the adsorption configuration, the molecule is located perpendicular to the surface on top of a W atom.The Z -values are the corresponding N 2 center of mass distance from the surface.For completeness, the adsorption energies and distances obtained when no surface relaxation is allowed are written in the brackets below.The equilibrium N-N internuclear distance obtained in all the calculations is 1.13 Å.The experimental adsorption energy is estimated in ∼450 meV.neither the main, ingredient causing the success of the RPBE-CRP in predicting the S 0 (E i ) values at normal incidence.The differences between the PW91 and the RPBE PESs regarding point (ii) are very important.More precisely, the discrepancy between the experimental and the PW91-CRP initial sticking coefficient S 0 seems to be related with the incorrect theoretical imbalance between the barriers to reflection and to dissociation that strongly contributes promoting the latter process in the MD calculations performed with the PW91-CRP PES. 34Since both barriers are more similar with RPBE, the efficiency of the indirect dissociation is greatly diminished. 37In order to analyze, here, the properties of the vdW-functionals regarding point (ii), we calculate the MEDP using the NEB method as implemented in VASP. Figure 2 shows the results obtained with each functional, including the PW91 and RPBE also.The zero energy is referred to the N 2 molecule far away from the W(110) surface and neglecting the ZPE.In practice, we simply subtract in each case the energy of the isolated relaxed surface and that of the isolated gas-phase N 2 as already done to define the adsorption energy in Eq. ( 1).The MEDP is calculated for each xc-functional within the FS approximation and also allowing relaxation of the two topmost W(110) layers (nonfrozen surface approximation, NFS).In all cases, three are the number of images used to connect the initial point, which corresponds to N 2 in the adsorption well, and the final point, in which the molecular axis points along the [ 11 1] direction and each N is over the bridge site (see sketches in Fig. 2).
In general, the FS compared with the NFS results tends to overestimate the energy barriers (underestimate the adsorption wells) by ∼100 meV, but qualitatively the results are very similar.Note in passing that the correlation found between the height of the reaction barrier and the depth of the adsorption well follows the Brønsted-Evans-Polanyi relation.Focusing on the FS results obtained with the vdW-functionals, the vdW-DF curve with an energy barrier of about 0.37 eV measured from the bottom of the adsorption well is the energetically less favorable one to dissociation.Next, the optPBE-and optB88-vdW predict lower energy barriers of ∼0.22 eV and ∼0.17 eV, respectively.The lowest barrier corresponds to the optB86b-vdW with ∼0.12 eV.Compared with the vdW-DF and the optimized vdW-functionals, the vdW-DF2 appears to be so corrugated and repulsive that it makes it difficult to find with the NEB method a MEDP connecting directly the N 2 adsorption well with any of the dissociation configurations identified in Ref. 34.In fact, the values of the MEDP plotted in Fig. 2 for the NFS calculation could only be obtained by using as starting guesses for the NEB images the corresponding converged positions found for the other functionals.Despite the same procedure allowed us to calculate the MEDP for the FS case, the energy obtained for the N atoms over bridge sites (image 4, see Fig. 2) is unexpectedly high in comparison with the results found with the other functionals and with the vdW-DF2 in the NFS calculation.At this point, it is worth to note that the vdW-DF2 predicts a lattice constant for bulk W that is 2% stretched with respect to the PW91 value of 3.17 Å.Thus, the weakly attractive energy value of image 4, in which each N is over a bridge site very close to two W atoms, can simply be an artifact of the FS calculations, which forces an artificially compressed W(110) surface.For this reason, we consider useless to show the FS results for vdW-DF2.
The existence of a distinct indirect channel to dissociation was reported in temperature programmed desorption experiments performed in this system. 72The authors found that at any N 2 coverage about 10% of the pre-adsorbed N 2 dissociates on the surface during desorption scans at 128-150 K. Importantly, they also suggested that the activation energies to dissociation and to reflection have to be similar.Comparing these experimental observations with the theoretical MEDP of Fig. 2, it is clear from the NFS calculations that only the RPBE and the vdW-DF(2) MEDPs are compatible with these experiments.At this point, it is worthy to note that the similarities found with RPBE and vdW-DF seem to be a consequence of having a similar description of the exchange energy term, which dominates the short-range interaction.
This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 158.227.180.88On: Wed, 18 Feb 2015 16:09:29 Figure 3 provides information on the properties of one of the important configurations involved in the direct dissociation channel that dominates at 60 • incidence and at normal incidence and large E i .The molecule is over the hollow site with its axis oriented parallel to the surface and pointing along the [ 11 1] direction.Note that this is precisely the configuration that would directly correspond to the dissociated N 2 depicted in Fig. 2 (see sketch of image 4).The contour plots of Fig. 3 show that, independently of the vdW-functional, this configuration is characterized by an early energy barrier that appears at distances between 2.0 and 2.5 Å.After overcoming this barrier, the molecules will easily dissociate guided by the rather large adiabatic forces that tend to elongate the internuclear axis, exception made of the vdW-DF2 PES.In this case, there is still a second (late) barrier of about 90 meV.This barrier is probably not much affecting the direct dissociation because the incoming N 2 should first overcome a barrier of ca.0.67 eV, and, therefore, it should also have enough energy to pass this second barrier.Nonetheless, it is apparent that the existence of such a "double-barrier" structure will contribute to make dissociation more complicated.Furthermore, it would also explain the above mentioned difficulties we encounter in finding the MEDP with the NEB method.Comparing all the vdW-functionals, we observe similar trends to what we find for the adsorption well and the MEDP.The highest energy barriers are obtained with the vdW-DF2 and the vdW-DF in this order.The optPBE-and optB88-vdW energy barriers of ca.0.4 eV and 0.3 eV follow and the lowest value of 0.26 eV is obtained with the optB86b-vdW.The corresponding contour plots for PW91 34 and RPBE are similar to the ones of Fig. 3, except that the early barriers are around 0.51 and 0.81 eV, respectively.Thus, comparing the RPBE and the vdW-DF(2), the effect of the corrected correlation is clearly a reduction of the energy barrier located at 2.0-2.5 Å.
All the static analysis above has been guided by the results obtained from the MD calculations performed on the PW91and RPBE-CRP PESs and, therefore, it provides us a preliminary and meaningful view of the pros and cons of each vdW-FIG.3. Contour plots E(r, Z ) of the N 2 /W(110) PES for one of the configurations involved in the direct dissociation mechanism.The molecule is oriented parallel to the surface along the [ 11 1] direction with its center of mass over the hollow site (see image 4 in Fig. 2).The value of the EB is written in each plot.Contours as in Fig. 1.
functional in describing the complex dissociative dynamics found for N 2 on W(110).On the one hand, the vdW-DF and -DF2 functionals that, in common with the RPBE, describe better the experimental observations regarding the adsorption well and the balance between desorption and dissociation of the chemisorbed N 2 will probably give a better account of the indirect dissociation channel (we recall that vdW-DF and RPBE have a similar exchange term).Yet, the question that remains is if the large energy barriers that appear in those configurations involved in the direct channel are compatible with the dissociation probabilities measured at off-normal incidence angles.On the other hand, the (exchange) optimized vdW-functionals show similar properties to those found with the PW91-CRP PES, i.e., an excessively deep adsorption well from which the dissociation path compared with desorption is energetically favored.Then, the question here is if the indirect channel will be unrealistically enhanced, as it was the case for the PW91-CRP calculations.

B. Dissociative dynamics
The initial (dissociative) sticking coefficient S 0 of N 2 on W(110) was measured in Refs.73 and 74.The molecular beam experiments performed at various incidence angles Θ i and kinetic energies E i showed in all cases that S 0 is less than 3 × 10 −3 for E i ≤ 0.3 eV.Above this energy, S 0 steadily increases reaching values of ca.0.4 (0.3) for E i ∼ 2.25 eV and Θ i = 0 • (60 • ).The two experimental Θ i = 0 • data sets differ at high E i values by a factor of ∼1.3.Nevertheless, the overall qualitative agreement is good and the S 0 values at low E i agree quantitatively.In particular, the activation energy and S 0 curve increase are the characteristic features of the system and thus we consider them as valid benchmarks for our simulations.The experimental data together with the results obtained from the classical MD calculations run with the PW91-CRP 33 and RPBE-CRP 37 PESs are shown in Fig. 4. The so far discussed limitations of the theoretical PESs, which are based on semilocal-DFT, are readily apparent in this figure.Neither the PW91-CRP nor the RPBE-CRP can fully reproduce the experimental S 0 measured at these two extreme incidence polar angles.In particular, the PW91-CRP dissociation probability agrees rather well with the experimental S 0 reported at Θ i = 60 • but badly compares with the one measured at normal incidence and low incidence energies (E i < 0.5 eV).Also observed in Fig. 4, the theoretical S 0 at normal incidence exhibits a peaked structure in this energy range with a maximum value of ∼0.13 at E i ≃ 0.2 eV that is at variance with the experimental values S 0 < 3 × 10 −3 measured at the same incidence conditions.In contrast, the RPBE-CRP that provides a rather good description of the experimental dissociation probabilities at normal incidence completely fails at Θ i = 60 • .The RPBE-CRP behavior was attributed to the more repulsive character of the RPBE xc-functional that corrects the tendency of the PW91 functional for overbinding and, thus, provides better adsorption energies and a more realistic balance between dissociation and desorption barriers.In spite of it, the RPBE produces at the same time excessively large energy barriers at the entrance channel.The latter has dramatic consequences at large incidence angles because it prevents the molecules from getting close to the surface and points to the  4. Initial (dissociative) sticking coefficient S 0 as a function of the incidence energy E i for normal and 60 • incidence.The red and green lines reproduce the results obtained with the PW91-CRP 33,34 and the RPBE-CRP 37 PESs, respectively.Experiments: white squares, data at 800 K from Ref. 73 that reported error bars of about ±25% for the lowest and ±10% for the highest E i ; black squares, data at 800 K from Ref. 74.   attractive vdW interaction as a firm candidate to improve the theoretical results.
Before starting AIMD calculations that include the vdW interaction, we have first verified up to what extent the conclusions of Refs.34 and 37 that were obtained from classical dynamics calculations performed with the PW91-and RPBE-CRP PESs are also reproduced by AIMD using the same semilocal xc-functionals.Figure 4 shows that in general there is a good agreement between the previous CRP results and the present AIMD values of the sticking coefficient.Note, however that the peaked structure obtained with the PW91-CRP at normal incidence and low E i is less pronounced in the AIMD-PW91 calculations.In spite of it, the analysis of the AIMD dissociating trajectories (see below) shows that there are no qualitative changes in the dissociation dynamics, i.e., it still proceeds through the indirect channel that was already identified in Refs.33, 34, and 37. We have also checked that the percentage of molecules reaching the region Z ≤ 2.5 Å is equal (65%) in both calculations, AIMD and MD with the PW91-CRP PES. 37The discrepancy between the AIMD-PW91 and the PW91-CRP results is due to the different percentage of molecules that dissociate after reaching this region: about 18% in the PW91-CRP, 37 but 6.5% in AIMD-PW91.Hence, we attribute the difference in S 0 to the interpolation errors of the N 2 /W(110) CRP-PES that for the tilted configurations can be of 100 meV at Z < 2 Å.
Next, we proceed to investigate the role of the vdW interaction in this system.To optimize the computational effort required to calculate the dissociation probability for the different incident energies E i and polar angles Θ i used in the experiments, we start by running AIMD calculations for a selection of xc-functionals and two representative initial conditions: (i) Θ i = 0 • and E i = 0.25 eV and (ii) Θ i = 60 • and E i = 2.28 eV.Independently of the CRP PES used in the MD simulations, it was found that for the incidence conditions of case (i) dissociation is dominated by the indirect channel.Therefore, the corresponding theoretical S 0 value will depend on the quality of the PES provided by each xc-functional in the region close to the surface and accessing the adsorption well.In contrast, dissociation at Θ i = 60 • (case (ii)), being direct, is dominated by the energy barriers at the entrance channel.Therefore, condition (ii) serves as a meaningful test of the PES quality at intermediate/large distances from the surface (2-4 Å).For these reasons, fulfillment of conditions (i) and (ii) together provides a robust and stringent test for the quality of xc-functionals in the broad configurational space.
The results of the AIMD calculations for conditions (i) and (ii) that have been performed with the vdW-DF, optPBE-vdW, and optB86b-vdW are shown in Table II.These functionals have been selected as representative of the different energy landscapes relevant to the direct and the indirect dissociating dynamics that have been identified with our previous static analysis in Sec.III A. In other words, we expect that vdW-DF2 and optB88-vdW will provide results similar to those of vdW-DF and optB86b-vdW, respectively.We also recall that vdW-DF and RPBE have similar exchange terms and the same happens for optPBE-vdW and PW91.Noticeably, the AIMD results of Table II are in line with the expectations that one can infer from the limited analysis of the static PES properties.Thus, the vdW-DF that was quite promising because of its good description of the adsorption well seems to correct the PW91 pitfall for condition (i).However, it largely fails in TABLE II.Dissociation probabilities (S 0 ) and average time of the dissociation (⟨t d ⟩) and the reflection (⟨t r ⟩) events obtained from AIMD calculations performed with different xc-functionals.The standard error of S 0 and the root mean standard deviation of ⟨t (d, r ) ⟩ are written in parenthesis.Time given in ps.

Experiments
PW91 RPBE vdW-DF optPBE-vdW opt86b-vdW reproducing the experimental S 0 at 60 • incidence.The vdW-DF2 is expected to also underestimate the dissociation probability at 60 • incidence in view of its overall more repulsive character already pointed.At first sight, the vdW-DF behavior seems similar to the one found with the RPBE-CRP MD.
Nonetheless, the analysis of the distances Z min at which the classical turning point occurs reveals that there are important differences.The excessive repulsive character of the RPBE at large distances from the surface is confirmed by the fact that 58% of the incoming trajectories are reflected above 2.5 Å.
Clearly the improved correlation introduced in the vdW functionals comes to correct this unrealistic behavior, since only five trajectories out of 99 are reflected above 2.5 Å with the vdW-DF functional.However, though closer to the surface, the molecules will be reflected in the end (we will come to this point below).In the other extreme, but as expected, the optB86b-vdW is facilitating in excess the dissociation of the dynamically trapped molecules, i.e., the one under condition (i).The results obtained with the optPBE-vdW are in better accordance with the experimental values of S 0 .In principle, this functional seems to reconcile the pros of the PW91 at 60 • incidence and of the RPBE at normal incidence.Still the value of 0.022 found at 0 • incidence (condition (i)) is about a factor of 10 larger that the S 0 ∼0.003 reported in experiments.
To get further insight into the changes introduced by the vdW-functionals in the N 2 /W(110) PES, we analyze in more detail different properties of the dissociation and reflection events that are predicted by each functional.Starting with the dissociation process, it is clear that the dissociation dynamics under conditions (i) and (ii) are completely different.We find that for all the studied xc-functionals, the time for dissociation is notably longer that the time for reflection in case of condition (i), while the duration of both events is rather similar for incidence condition (ii).The values of the average dissociating (⟨t d ⟩) and average reflecting (⟨t r ⟩) time shown in Table II serve to illustrate this general observation.The long dissociating times are an obvious fingerprint of the indirect dissociation channel that was identified from MD calculations performed on top of the PW91-and RPBE-CRP PESs. 34,37That this mechanism is at work for condition (i) is also confirmed by studying the evolution of the molecular orientation and position of the N 2 center of mass for the corresponding dissociating trajectories.The snapshots of Fig. 5 that show the position of the center of mass at different distances from the surface are representative of how the dissociating N 2 molecules move along the incoming part of the trajectory.Following the whole evolution, we observe that the molecules are in all cases dynamically trapped in the attractive parts of the PES, which exist between 2.0 and 3.0 Å from the surface.Except for the optB86b-vdW, the trapping region corresponds to configurations that are close to the top-vertical one associated to the N 2 adsorption well (see Fig. 1).This trapping, however, does not necessarily occur in the initial incoming part of the trajectory, but it happens in some cases after colliding few times with the surface.Remarkably, the snapshots also show that there are other regions apart from the attraction to the N 2 well that will permit dissociation in the optB86b-vdW case.In principle, a broader configurational space would facilitate dissociation and, hence, it would also explain the appreciably shorter mean dissociating time (0.69 ps) in comparison with the mean times obtained with either the PW91 (2.59 ps) or the optPBE-vdW (1.57ps).The final dissociating configuration that, for practical purposes, is defined as the one at which the internuclear N-N distance is ∼2.0 Å, is shown in the left-large panel of Fig. 5.There we observe the same configurations over hollow and bridge sites that were identified with the RPBE-CRP and PW91-CRP PESs.As anticipated from the snapshots above, the configurational space leading to dissociation is clearly much broader in the optB86b-vdW PES.The latter explains the larger dissociation probabilities obtained with this functional.
Focusing next in the reflection process, Table II shows that the average time used in a reflection event is rather similar in all cases.In particular, the average values are close to 0.3 ps for incidence condition (ii) and they are slightly larger for condition (i) due, partly, to the lower E i of the molecules in this case.Also note that for condition (i), the values of the average time for reflection are in general more spread, i.e., the RMSD (root mean standard deviation) is slightly larger than for condition (ii).The analysis of the classical turning point will show that this is in part a consequence of the different distances from the surface at which the molecules are scattered back under incidence condition (i).Furthermore, they are also indicative of the existence of molecular trapping.This point is particularly marked for optPBE-vdW and it explains the large average time and RMSD values found.We observe that with this functional, there are two different reflection regimes.While most of the molecules are rapidly reflected at an average time of about 350 fs, there is an important number of molecules that are characterized by a distinctively larger average reflection time of about 950 fs.The latter are molecules that, after bouncing back and forth over the surface few times, finally escape into vacuum instead of being dissociated.
The comparative analysis of the classical turning point Z min , i.e., the minimum distance from the surface reached by the reflecting N 2 , and of the positions over the surface where the molecules are scattered back provides meaningful information on where the repulsive parts of each PES are localized.The Z min -distributions for each xc-functional and conditions (i) and (ii) are shown in Fig. 6.The role of the vdW attraction at large Z allowing the molecules getting close to the surface is evident when comparing the distributions of any vdW-functional with those obtained with the PW91 and, in particular, with the RPBE AIMD.Thus, the artificially large barriers created by the RPBE at 2.5-3.0Å from the surface disappear with all vdWfunctionals.Comparing the distributions of the three vdWfunctionals, we observe that reflection typically occurs at larger distances with vdW-DF evidencing its less attractive character, while optB86b-vdW would be the most attractive.Focusing on the distributions for condition (i), the double-peaked structure observed in all cases, exception made of RPBE, reveals the existence of two separated Z-regions for reflection that are also identified with the two distinct reflection times discussed above.Regarding the XY positions over the surface at the turning point (not shown), there are no qualitative differences among the xc-functionals on where the repulsive parts are located, at least within our statistics.
Putting together the analysis of the dissociating and reflecting trajectories, we can extract the following conclusions regarding the limitations of the vdW-DF and the vdW-optB86b functionals in describing the N 2 /W(110) system.Clearly the optB86b-vdW functional is creating a very attractive PES that, at variance with the experimental data, facilitates in excess the reactivity of N 2 on W(110).The vdW-DF functional, which in common with RPBE gives a reasonable description of the configurational space relevant for the indirect channel, improves the N 2 -W(110) interaction at intermediate-large distances, where RPBE fails.Thus, the barriers that with RPBE appear above 3 Å from the surface are greatly reduced as a consequence of the corrected correlation.However, the unsatisfactory result obtained at 60 • for the sticking coefficient, which is similar to the RPBE S 0 value, suggests that such a reduction is insufficient at distances of 2.0-3.0Å from the surface, i.e., where the barriers ruling the direct dissociation mechanism are located.This idea agrees, for instance, with the results shown in Fig. 3 for one of the configurations involved in the direct dissociation channel.As discussed in Sec.III A, the energy barrier of about 0.81 eV found with the RPBE is reduced by 0.2 eV with the vdW-DF.However, taking into account that the PW91 is providing a good description of the sticking coefficient at 60 • , the vdW-DF barrier of 0.61 eV is still too high when compared with the PW91 barrier of 0.51 eV.
The results of Table II suggest that only the optPBE-vdW seems to improve the description of the N 2 -W(110) interaction.To further scrutinize this conclusion, we have calculated the dissociation probability for few more incidence conditions using the optPBE-vdW.All the results are plotted in Fig. 4 to facilitate the comparison with experiments and with the AIMD-PW91 and AIMD-RPBE calculations.At 60 • incidence, the dissociation probabilities with optPBE-vdW and with PW91 are very similar.Still, the threshold energy for dissociation at this incidence is not correctly described by optPBE-vdW, since it predicts a value S 0 = 0.06 for E i = 0.5 eV.This result remarks the incorrectness of the low energy barrier (0.38 eV) found for the dissociating configuration of Fig. 3.At normal incidence, the improvement over the PW91 at low E i is neither clear.While the experimental sticking coefficient takes values below 3 × 10 −3 for E i ≤ 0.5 eV that was reasonably described by the optPBE-vdW at 0.25 eV, the value obtained at 0.4 eV is even larger than the PW91 one.
In view of all these results, we have to conclude that even if the vdW interaction is crucial in the dissociative adsorption of N 2 on W(110) to have an accurate description of the long-distance region accessing the N 2 adsorption well, the presently studied vdW-functionals are not fully satisfactory to treat this problem.In this respect, it is worthy to remark that an adequate treatment of the dissociative adsorption of N 2 on W(110) is a particularly demanding problem because it equally depends on the PES properties at short and long molecule-surface distances.Similar difficulties are found in reproducing the experimental sticking coefficient of H 2 on Pd(111) at incidence energies below 125 meV, at which dynamic trapping dominates. 50Among the tested functionals, the PBE-vdW reasonably reproduces the experimental S 0 at incidence energies ≥125 meV, but fails in the low energy range.
Other reasons behind the disagreement between experiments and theory might be related with those experimental conditions not considered in our simulations, such as the rovibrational state of the incident N 2 and the creation of phonon excitations.Assuming a rovibrationally cold N 2 beam, we would not expect in view of the quasiclassical calculations of Ref. 34 that the ZPE (not included in our present AIMD) can be at the origin of the discrepancies.More likely, the role of phonons can be another factor to consider, particularly, because it might affect the indirect dissociation mechanism.In order to obtain a qualitative understanding on how the experimental surface temperature of 800 K may affect our conclusions based on the frozen surface approximation, we have used the PW91and RPBE-CRP PESs to perform MD calculations based on the generalized Langevin oscillator model (GLO), [75][76][77] which accounts for energy exchange between the molecule and the surface lattice.Using the same GLO parameters of Ref. 29 for the W(110) surface, we find that the sticking coefficient at 60 • incidence may change by a factor of 1.1-1.25,while it may increase as much as a factor of 2-2.5 for dissociation at normal incidence and low E i .These results suggest that our conclusions on the excessively attractive character of the exchangeoptimized vdW functionals will remain, since the dissociation probability at normal incidence and low E i will be even higher.Similarly, our estimated factor of 1.25 is insufficient to improve the performance of vdW-DF at 60 • incidence.

IV. SUMMARY
In summary, we have analyzed the performance of different vdW xc-functionals to describe the dissociative dynamics of N 2 on W(110).The analysis was motivated by the fact that neither the PW91 nor the RPBE semilocal xc-functionals were able to reproduce the experimental results for the dissociative adsorption in all the range of energies and angles of incidence.The question to be answered was whether the neglect of the vdW interaction is responsible for the disagreement between theory and experiments.Starting with the static analysis of the potential energy surface, we find that the non-local correlation term permits to correct the unrealistic energy barriers that appear above 3 Å from the surface with the RPBE.Our analysis based on AIMD calculations shows, however, that none of the used vdW-functionals provide a significant improvement over the results obtained with the semilocal xc-functionals, when compared with the available molecular beam data.The attractiveness of the vdW interaction at long distances affects the dynamics making the molecule get closer to the surface.In spite of that, reactivity is not significantly affected by the corrected long range behavior of the potential, at least with the presently used vdW-functionals.
The opt86b-vdW functional is too attractive in the whole configuration space and overestimates greatly the reactivity in all the incidence conditions.At the other extreme, the vdW-DF functional, though more attractive at long distances than the semilocal RPBE functional, shows roughly the same behavior concerning reactivity.As RPBE, vdW-DF reproduces the experimental low dissociation probability at low energies and normal incidence and underestimates it for 60 • incidence.This shows that this functional is still too repulsive for the configurations that lead to direct dissociation at intermediate distances.From the studied vdW functionals, the optPBE-vdW is the one that gives results closer to the experiments in all the analyzed incidence conditions.However, within our statistics, it is not clear whether it represents an improvement over the semilocal PW91 concerning reactivity, as it seems to overestimate the indirect channel to dissociation as well as the energy threshold for direct dissociation.
All in all, none of the xc-functionals analyzed so far, irrespective of including the vdW interaction or not, are able to fully characterize the multi-configurational N 2 /W(110) PES with the accuracy required to reproduce the experimental reaction probabilities in all the range of incidence conditions.Functionals such as vdW-DF and RPBE that describe properly the indirect channel to dissociation, due to a correct theoretical balance of the barriers to dissociation and reflection from the molecular adsorption well, overestimate the barriers for the configurations that lead to direct dissociation.On the contrary, functionals such as PW91 and optPBE-vdW, which seem to describe more accurately the direct dissociation channel, tend to overestimate the indirect channel due to the theoretical imbalance of the barriers to dissociation and reflection from the molecular adsorption well.Inclusion of the vdW interaction implies a more attractive PES at long distances from the surface and, therefore, a closer approach of the molecules to the surface but does not seem to improve substantially the results for the reaction probabilities.
Finally, we emphasize the interest of state-of-the-art gassurface molecular dynamics simulations as a valuable additional tool to test the quality of the developed xc-functionals, since the theoretical outputs can be compared with existing reliable data from molecular beams experiments.There is a remarkably large body of literature devoted to the performance of vdW-functionals in the study of molecule-surface interactions.However, to our knowledge only few studies, namely, Refs.30 and 50 and the present one, have considered the dynamical aspects of the interaction, which remains a scenario where the role of the vdW interaction is still to be addressed.In fact, our analysis shows that a functional can accurately describe adsorption energies and energy barriers, but at the same time can fail in providing an overall description of the multi-configurational space that ultimately rules the gassurface interaction.

FIG. 2 .
FIG. 2. MEDP calculated with the nudged elastic band method using different xc-functionals as indicated by the different symbols.Results obtained with the frozen surface approximation (FS-calculation, top-left figure) and allowing relaxation of the two topmost layers (NFS-calculation, bottom-left figure).Lines to guide the eye.The molecule is initially (image 0) in the adsorption well that exists atop the W atom and ends up (image 4) parallel to the surface with the center of mass over the hollow site and the molecular axis pointing along the [ 11 1] direction.The calculated configurations along the MEDP (images 0-4) are depicted in the right panels.

FIG. 5 .
FIG.5.Center of mass position over the (2 × 2) surface unit cell of the dissociating N 2 as obtained from AIMD performed with different xc-functionals for the incidence condition Θ i = 0 • and E i = 0.25 eV.Large grey circles represent the W atoms at the surface topmost layer.Left figure: XY-positions at which the internuclear N-N distance gets the value of 2 Å.Right panels: XY-positions when first reaching the distance Z = 6.5, 3.0, 2.6, and 1.5 Å from the surface.

FIG. 6 .
FIG. 6. Classical turning point distributions Z min of N 2 scattered off W(110) as obtained from the AIMD calculations performed with the semilocal PW91 and RPBE and with different vdW-functionals.The incidence conditions are Θ i = 0 • and E i = 0.25 eV (top panel) and Θ i = 60 • and E i = 2.28 eV (bottom panel).
. 1. Contour plots E(r, Z ) of the N 2 /W(110) PES for the molecular adsorption configuration calculated with different vdW-functionals.Black solid (white dashed) contour lines, separated by 0.2 eV, correspond to positive (negative) potential energy values.The zero potential energy, which corresponds to the N 2 gas-phase minimum, is shown by white solid lines.