Performance of Made Simple Meta-GGA Functionals with rVV10 Nonlocal Correlation for H2 + Cu(111), D2 + Ag(111), H2 + Au(111), and D2 + Pt(111).

Accurately modeling heterogeneous catalysis requires accurate descriptions of rate-controlling elementary reactions of molecules on metal surfaces, but standard density functionals (DFs) are not accurate enough for this. The problem can be solved with the specific reaction parameter approach to density functional theory (SRP-DFT), but the transferability of SRP DFs among chemically related systems is limited. We combine the MS-PBEl, MS-B86bl, and MS-RPBEl semilocal made simple (MS) meta-generalized gradient approximation (GGA) (mGGA) DFs with rVV10 nonlocal correlation, and we evaluate their performance for the hydrogen (H2) + Cu(111), deuterium (D2) + Ag(111), H2 + Au(111), and D2 + Pt(111) gas-surface systems. The three MS mGGA DFs that have been combined with rVV10 nonlocal correlation were not fitted to reproduce particular experiments, nor has the b parameter present in rVV10 been reoptimized. Of the three DFs obtained the MS-PBEl-rVV10 DF yields an excellent description of van der Waals well geometries. The three original MS mGGA DFs gave a highly accurate description of the metals, which was comparable in quality to that obtained with the PBEsol DF. Here, we find that combining the three original MS mGGA DFs with rVV10 nonlocal correlation comes at the cost of a slightly less accurate description of the metal. However, the description of the metal obtained in this way is still better than the descriptions obtained with SRP DFs specifically optimized for individual systems. Using the Born-Oppenheimer static surface (BOSS) model, simulations of molecular beam dissociative chemisorption experiments yield chemical accuracy for the D2 + Ag(111) and D2 + Pt(111) systems. A comparison between calculated and measured E 1/2(ν, J) parameters describing associative desorption suggests chemical accuracy for the associative desorption of H2 from Au(111) as well. Our results suggest that ascending Jacob's ladder to the mGGA rung yields increasingly more accurate results for gas-surface reactions of H2 (D2) interacting with late transition metals.

Introduction

In heterogeneous catalysis, the rate-limiting step is often the dissociative chemisorption of a molecule on a surface.[1,2] The dissociation of simple hydrogen (H2) and nitrogen (N2) molecules constitutes important steps in the production of ammonia and syngas.[3−5] The dissociation of H2 is also relevant to the industrial synthesis of methanol from CO2 over a Cu/ZnO/Al2O3 catalyst, for which the dissociation of H2 is considered to be a rate-limiting step.[6−8] Calculating chemically accurate barrier heights[9] for rate-controlling reactions to obtain accurate rates of the overall reaction network[10] potentially has a large financial impact on the chemical industry since it allows theoretical screening for more efficient catalysts.[11]

Currently, density functional theory (DFT) is the only method that is computationally cheap enough to map out full potential energy surfaces (PESs) for gas-surface reactions. Development of density functionals (DFs) that can accurately describe dissociative chemisorption reactions on surfaces is important to increase the predictive power of DFT. DFs constructed using the generalized gradient approximation (GGA) that provide chemically accurate results for specific gas-surface reactions and that in some cases show transferability to chemically related systems are based on the semiempirical specific reaction parameter (SRP) approach to DFT (SRP-DFT).[12−15] However, DFs at the GGA level are always a compromise between a good description of the molecule and of the metal,[16] despite efforts to construct GGA-based DFs[17] or nonseparable gradient approximation DFs[18] that perform equally well for both solids and molecules. A good description of the metal is crucial to calculate accurate barrier heights since the barrier height might depend on the interlayer distance of the two topmost metal layers[19−21] and the amplitude of thermal motion of the metal atoms in the top layer.[22,23]

In previous work, we developed semilocal MS-PBEl, MS-B86bl, and MS-RPBEl meta-GGA (mGGA) DFs[24] based on the made simple (MS) formalism,[25,26] which yield a description of the metal that is comparable in accuracy to that obtained with PBEsol[27] DF. Additionally, the DFs provide a chemically accurate description of molecular beam experiments on the dissociative chemisorption of H2 (D2) on Cu(111)[24,28−30] and a near-chemically accurate description of similar experiments on D2 + Ag(111).[24,31] The reason behind this improved overall performance of mGGA-based DFs over GGA-based DFs is that mGGA DFs also depend on the kinetic energy density τ, which allows a DF to distinguish between regions of the electron density describing single (covalent), metallic, and weak bonds[32] via the dimensionless inhomogeneity parameter α.[25,26,32] This parameter has also been used in the construction of several other much used mGGA DFs, such as TPSS,[33] revTPSS,[34] RTPSS,[35] SCAN,[36,37] and mBEEF.[38] Several groups have now reported good simultaneous descriptions of lattice constants and adsorption energies,[38−40] or, more generally, energetics and structure,[25,26,41,42] when using mGGA DFs. The MS-RPBEl DF has also shown some success in describing the O2 + Al(111) system.[15]

In recent work, we also identified nonlocal correlation (NLC) as a key ingredient for a DF that can describe dissociative chemisorption of H2 (D2) with chemical accuracy on multiple metals[14] and not just on different crystal faces of the same metal,[43,44] which had previously only been demonstrated for reactions of CH4 with metal surfaces, i.e., Ni(111)[45] and Pt(111).[46] Here, we combine the three previously developed MS mGGA DFs with rVV10[47] nonlocal correlation to obtain the MS-PBEl-rVV10, MS-B86bl-rVV10, and MS-RPBEl-rVV10 DFs, and we will evaluate their performance for the H2 + Cu(111), Ag(111), Au(111), and Pt(111) systems. The three original MS mGGA DFs,[24] which we combine with rVV10[47] nonlocal correlation, show no van der Waals (vdW) interactions for H2 interacting with transition metals,[14] which is the best-case scenario to complement a semilocal exchange–correlation functional with (r)VV10 nonlocal correlation according to Vydrov and Van Voorhis.[48]

The PESs we computed with the three new DFs are subsequently used in quasi-classical trajectory (QCT) calculations. In the dynamics calculations, we use the Born–Oppenheimer static surface (BOSS) model, which is known to work well for activated H2 dissociation on cold metals.[49−53] Calculations that incorporate surface motion show that the impact of surface atom motion (phonons) can be neglected due to the effect on the reaction probability being small for the low-surface-temperature experiments considered here.[19,21,51,54] It is also justified to neglect the effect of electron–hole pair excitation on the reaction probability, as its effect on sticking has previously been shown to be small in calculations on H2 + Cu(111),[55−57] Ag(111),[58−60] and Ru(0001).[61] Previous research has also shown that for highly activated dissociation of H2 on cold metals, the difference between quantum dynamics (QD) and QCT calculations is marginal,[14,62,63] and there is also some evidence that the same observation holds for the nonactivated reaction of D2 + Pt(111) for all but the lowest translational energies.[13,64] Since our dynamical model is best suited to molecular beam dissociative chemisorption experiments, we will mainly compare to this kind of experiment[28−31,65,66] to assess the quality of the obtained DFs. These experiments have been performed for H2 + Cu(111),[28−30] Ag(111),[31] and Pt(111).[65−67]

Additionally, we will also compare to the associative desorption experiments that are available for the H2 (D2) + Au(111)[68] and Ag(111)[69,70] systems as in our previous work, by comparing the measured E0(ν, J) parameters characterizing the measurements with calculated E1/2(ν, J) parameters,[14] assuming detailed balance. Given that the DFs developed here are too reactive with respect to the H2 (D2) + Cu(111) system (as will be shown below), we will omit such an analysis for the recent associative desorption experiments[71] for this system here. For the H2 + Cu(111) system, it is known that the effect of surface motion cannot readily be ignored for specific observables at a high surface temperature[52] (Ts), and this may hold for the H2 + Au(111) and Ag(111) systems as well. Therefore, it is difficult to assess the quality of the developed DFs when using the BOSS model in comparison to high-surface-temperature experiments,[68,71] as will be done below. We also note that it is also possible to simulate associative desorption directly by running trajectories starting around the transition state using Metropolis sampling of the initial conditions[72−76] and that this has also been done for H2 and D2 desorbing from Cu(111). There are some limitations regarding these calculations: in earlier work,[73,74] a PES that is an approximate fit[77] to unconverged DFT calculations[78] was used. The statistical accuracy of the later work[76] is limited by the number of ab initio molecular dynamics (AIMD) trajectories that have been calculated.

The vdW well geometries obtained from our DFT calculations will be compared to experimental results, which are mostly obtained from the analysis of selective adsorption experiments.[79−89] In these experiments, an increase or a dip is observed in a peak associated with a rotational (rotationally mediated selective adsorption, RMSA[79]) transition or in a peak for a diffractive (corrugation mediated selective adsorption, CMSA[90,91]) transition if the translational energy passes through a value that overlaps with the energy difference between two hindered rotational or parallel translational metastable states, respectively. The H2 molecule is then trapped in the final state in the vdW well close to the surface.[81,86] The resonance energies can then be used to reconstruct the shape of the potential and thus to determine the vdW well depths and geometries. Concerning the systems investigated here, studies using experiments to analyze the vdW interaction have been performed for H2 + Cu(111),[88,89] Ag(111),[82,83,85] Au(111),[89] and Pt(111).[79−81,92]

This paper is organized as follows. Section covers the computational methods used, with Section discussing the coordinate system we use. Section details how semilocal MS mGGA DFs can be combined with rVV10[47] nonlocal correlation. In Section , we discuss the construction of the PESs, and Section details the QCT method. The way in which we calculate observables is discussed in Section , and the computational details are summarized in Section . In Section , we present and discuss our results. Section provides results regarding the description of the metal, and Section discusses static PES properties like the vdW well geometries and barrier heights. In Section , we show a comparison of our results for simulated molecular beam experiments on dissociative chemisorption with experimental sticking probabilities, and in Section , we compare our results with the outcome of associative desorption experiments. The transferability of the newly developed DFs is discussed in Section . Section concludes this paper.

Methodology

Coordinate System

In the dynamics calculations, we use the BOSS model,[12] meaning that we make the Born–Oppenheimer approximation and keep the surface atoms fixed at their ideal lattice positions. We only take into account the six degrees of freedom (DOF) of the H2 molecule (see Figure a). We use molecular coordinates in which the center of mass (COM) coordinates X,Y describe the lateral position of the molecule and Z describes the molecule–surface distance. The remaining DOFs are H2 bond length r, polar angle θ, and azimuth ϕ defining the orientation of the molecule relative to the surface (see Figure a). The geometry of the (111) face of a face-centered cubic (fcc) metal together with its high-symmetry sites is shown in relation to the coordinate system used in Figure b.

Figure 1

COM coordinate system used for the description of the H2 (D2) molecule (a). Unit cell of a (111) face of an fcc metal together with the high-symmetry sites as well as the relationship with the coordinate system chosen for H2 (D2) relative to the (111) surface (b). Origin of the COM coordinate system (X,Y,Z) = (0,0,0) is the center of an atom in the top layer of the surface. We define the polar angle and azimuth as θ = 90° and ϕ = 0°, respectively, corresponding to molecules parallel to the surface pointing along the X (or equivalent U) direction. The hexagonal close-packed (hcp) and fcc hollow sites correspond to metal atoms in the second and third layers. Note that the colored triangle marks the irreducible wedge of the (111) unit cell.

Combining Made Simple Meta-GGA Exchange–Correlation with rVV10 Nonlocal Correlation

The form of the rVV10[47] nonlocal correlation functional is similar to that of the Rutgers–Chalmers vdW-DFs[37]Here, n(r) is the electron density and Φ(r,r′) is the kernel describing the density–density interactions.[37] We note that of course the term in the above equation with β in it is in fact local, but writing this expression in this way will facilitate subsequent discussions of how semilocal functionals (SLFs) can be added to the functional defined by eq to obtain full exchange–correlation functionals with nonlocal rVV10 correlation added in. The parameter β is not present in the vdW-DF1[93] and vdW-DF2[94] nonlocal correlation functionals and is here taken to be β = 1/32(3/b)3/4 so as to ensure that Ecnon-local is zero for the homogeneous electron gas.[37] In using the full exchange–correlation functional named rVV10, the nonlocal correlation (NLC) rVV10 functional is appended to the following semilocal functional (SLF)[37,47,95,96]Here, ExrPW86 is the exchange part of a refitted version of the PW86 functional[97] and EcPBE is the PBE correlation functional.[98]Equation also defines the semilocal exchange–correlation functional to which Vydrov and Van Voorhis[48] appended their NLC VV10 functional to obtain the full exchange–correlation functional now referred to as the VV10 functional. Sabatini et al.[47] obtained the NLC rVV10 functional of eq by making a minor change to the NLC VV10 functional[48] in a clever way to make it amenable to efficient evaluation by the algorithm due to Román-Pérez and Soler[99] that can also be used to speed up the evaluation of the vdW-DF1 and vdW-DF2 density functionals of the Lundqvist–Langreth group.[93,94] To reproduce the original VV10 results as closely as possible, Sabatini et al.[47] changed one of the empirical parameters in the NLC rVV10 functional, i.e., the b parameter, from the original VV10 value of 5.9 to the rVV10 value of 6.3. Here, the b parameter can be used to control the damping of the kernel at short range, while the other empirical parameter in VV10 and rVV10, C, can be used to obtain good values for the C6 dispersion coefficients describing the long-range vdW interaction. The C parameter is taken the same[47] in the NLC rVV10 as in the NLC VV10 functional.[48]

We note that there is some ambiguity associated with the SLF Sabatini et al.[47] originally appended their NLC rVV10 functional to. In a sentence saying that they were “following the original VV10 functional definition”, they provided an equation for the full rVV10 functional in which the SLF to which the NLC rVV10 functional was appended would be given byThe equation they presented suggested that in their SLF, PBE correlation was replaced with correlation from the local density approximation (LDA) (where Sabatini et al.[47] state that they used the functional as parameterized by Perdew and Wang[100]). This SLF happens to be the same as the one used with the nonlocal vdW-DF2 functional to obtain the full vdW-DF2 functional.[94] Regarding the SLF originally used, we have been informed by one of the authors of ref (47) (i.e., De Gironcoli) that the equation provided by Sabatini et al.[47] contained a misprint and that they in fact used the expression of eq instead (private communications).

The flexibility built in to the NLC rVV10 functional through the adjustable b parameter allows it to be used in combination with a number of exchange–correlation functionals, including mGGA functionals like the SCAN functional[36] and the B97M functional incorporated into the B97M-rV functional.[96] It is in this context that we use the NLC rVV10 functional, hoping that in this way we can obtain a good description of the long-range interaction, while hopefully keeping the medium-range interaction, which we think is reasonably described with the mGGA functionals[24] we will be testing as SLFs, intact, in the spirit of Peng et al.[37] In our work, the full exchange–correlation functional then takes the following formwhere EcrevTPSS is the revTPSS[34] correlation functional that is used in the original semilocal MS mGGA DFs we developed.[24]ExMS-mGGA can be either of the three MS mGGA exchange functionals we developed previously[24] based on the MS formalism.[25] In this formalism, one interpolates between two GGAs for two extreme scenarios, namely, a single-orbital system, which describes covalent bonds (Fx0(p;c)), and one in which the bonding is metallic (Fx1(p)).[24] The exchange enhancement factor of an MS mGGA DF then becomes[25]where p = s2, with s being the reduced gradient of the electron density,[25] and Fx1(p) and Fx0(p;c) are gradient enhancement factors that depend solely on p. The numerical parameter c is optimized to exactly reproduce the exchange energy of the hydrogen atom by canceling the spurious self-interaction present in the Hartree energy in this atom.[25] For both Fx1(p) and Fx0(p;c), three expressions[24] have been used, which are PBE-like,[98] RPBE-like,[101] and B86b-like,[102] in the sense that we use the gradient enhancement expression of the PBE,[98] RPBE,[101] and B86b[102] GGA DFs but with μ = 10/81 as was done in PBEsol.[27] The difference between Fx1(p) and Fx0(p;c) is that in the case of Fx0(p;c) we replace μp by μp + c everywhere,[24] as done earlier in ref (25). The interpolation between the two extreme cases then happens through a function of the inhomogeneity parameter f(α), with the inhomogeneity parameter being defined as[25,26]Here, τW is the Von Weizsäcker kinetic energy, which is equal to the kinetic energy density associated with a single-orbital electron density,[40] and τunif is the kinetic energy of the homogeneous electron gas. Note that α will approach unity as τ ≈ τunif and τW ≪ τunif for slowly varying electron densities, while α approaches zero for densities found in covalent bonding for which τ ≈ τW.[40] The expression for f(α) can be found in refs (24) and (25).

Above, we have already noted that the possibility to adapt the b parameter allows for flexibility in the combination of the NLC rVV10 functional with SLFs. In the past, several strategies have been used to arrive at a good choice of b. In perhaps the most rigorous approach, in the original papers presenting the full VV10[48] and rVV10[47] functionals, the b parameter was chosen to minimize the errors in the binding energies of weakly bonded dimers as present in the S22 database.[103] In a simplified procedure requiring fewer calculations, the b parameter has also been determined by demanding that calculations with the NLC rVV10 functional reproduce the Ar dimer energy curve determined with CCSD(T) calculations[96] as closely as possible.[37,47,95,96] In the development of functionals for specific purposes, the b parameter has also been fitted to more specific properties corresponding to these purposes. For instance, functionals have been developed that give good descriptions of layered materials by fitting the b parameter to obtain a good description of properties of these materials, after which the performance of the obtained functional is usually also tested on properties of other systems.[97,98] In the spirit of our SRP-DFT method, as described below, we follow an even more extreme approach to determine the b parameter.

The goal of the present paper is to investigate whether adding nonlocal rVV10 correlation to the MS mGGA functionals previously developed by us leads to functionals giving a better description of dissociative chemisorption of H2 on the noble-metal surfaces Cu(111), Ag(111), Au(111), and Pt(111). With this goal in mind, we investigated how closely we could reproduce the vdW interaction for the system for which the most accurate experimental results are available for this interaction (H2 + Cu(111), vdW well depths, and geometries are available from RMSA[79] or CMSA[90,91] experiments on this system). An additional reason for our choice of strategy is that most general-purpose DFs (at the GGA or mGGA level) cannot describe the interaction of H2 with transition-metal surfaces to within chemical accuracy (see, for example, Figure 6 in ref (24) or Figure 1a in ref (12)). Therefore, closely reproducing reference data for gas-phase dimers offers no guarantee that the obtained b value would be the best possible for H2 + transition-metal systems (although we will see below that this strategy would have worked for our case). However, we do check that the b parameter we adopt by considering the long-range attractive vdW interaction also yields a reasonably good description of the metal lattice constant for copper, which is a short- to medium-range interaction, to make sure that “the tail does not wag the dog”.[37]

Our tests on H2 + Cu(111) were first done with the MS-PBEl-rVV10 DF (see Figure S1). Adopting the b parameter of the original full rVV10 functional[47] (b = 6.3) yields a good description of the vdW well depth and minimum geometry, while a still reasonable lattice constant is obtained for copper (see Figure S1 and below). However, optimizing the b parameter for the MS-B86bl-rVV10 and MS-RPBEl-rVV10 functionals in this manner poses a dilemma. Using the small values of b suggested by a requirement of closely reproducing the H2 + Cu(111) vdW interaction leads to an underestimation of the copper lattice constant that we deem unacceptable (see Figures S2 and S3). This dilemma is illustrated in Figures S2 and S3 in the Supporting Information, in which the lattice constant, vdW well depth, and the position of the vdW minimum are shown as a function of b for the MS-B86bl-rVV10 and MS-RPBEl-rVV10 DFs. In these figures and Figure S1, the lattice constant has been recalculated for each value of b, after which the six-layer metal slabs are relaxed accordingly, and the vdW curve is calculated for a geometry in which H2 is parallel to the surface and above the top site. From Figures S1–S3, it is clear that reducing b yields smaller lattice constants and deeper vdW wells that are closer to the surface. Keeping these observations in mind, and noting that fitting the b parameters for the MS-PBEl-rVV10 DF to either the vdW well depth or the position of the minimum would have resulted in a value that is very similar to the original value of Sabatini et al.[47] (b = 6.3), we simply chose to adopt this value for all three functionals.

Finally, we note that the original MS mGGA exchange–correlation functionals appear to meet the same criterion as the semilocal exchange–correlation functional used by Vydrov and Van Voorhis[48] and Sabatini et al.,[47] i.e., that this functional does not yield an attractive long-range interaction (see Figure 3b in ref (14)). As Vydrov and Van Voorhis[48] point out: “it is preferable to combine VV10 with a functional that gives no significant binding in vdW complexes”. As our SLFs meet this criterion, we are not surprised that these SLFs combined with the NLC rVV10 functional yield either a good (with MS-PBEl) or still reasonable (with MS-B86bl or MS-RPBEl) description of the vdW interaction in H2 + Cu(111) with the choice of the original value of the b parameter.

Construction of the PESs

We use the corrugation reducing procedure (CRP)[104] to interpolate DFT results calculated on a grid to obtain a continuous representation of the PESs used in this work. Apart from using denser grids to improve the accuracy of the interpolated PESs, our method is analogous to the one used by Wijzenbroek et al.[105] In principle, we use the grids reported in the Supporting Information in ref (14).

In our previous work,[14] we have assessed the quality of the interpolated H2 + Cu(111) PES obtained using the B86SRP68-DF2 DF with ∼4900 randomly sampled geometries of H2 above the metal slab. Based on all of the randomly sampled points taken together, our CRP[104] interpolation had a root-mean-square (rms) error of 31 meV compared to the underlying electronic structure calculations. When only looking at the 3538 geometries that have an interaction energy of H2 with the surface lower than 4 eV, the rms error reduces to 8 meV (∼0.2 kcal/mol). Since we use the same interpolation grids as in our previous work,[14] we presume the accuracy of the obtained CRP PESs obtained in this work to be similar.

Quasi-Classical Dynamics

We compute observables using the quasi-classical trajectory (QCT) method.[106] This means that we take into account the quantum mechanical energies of the impinging H2 and D2 molecules in their initial rovibrational states. The method used is described more fully in ref (107). We integrate the equations of motion using the algorithm of Stoer and Bulirsch.[108]

To obtain reliable statistics, we propagate 200 000 trajectories per energy point when simulating a molecular beam experiment and 50 000 trajectories per energy point when calculating initial-state resolved reaction probabilities. Trajectories always start in the gas phase (Zgas = 8 Å). When r becomes bigger than some critical value (rc = 2.2 Å), the trajectory is counted as reacted. If during the propagation Z becomes bigger than Zgas, then the trajectory is counted as scattered. In all QCT calculations, we use a time step of dt = 0.001 fs. The reaction probability Pr is then calculated by dividing the number of reacted trajectories Nr by the total number of trajectories Ntotal

Computation of Observables

Molecular Beam Sticking

In the molecular beams, we simulate that the probability to find H2 with a velocity v in an interval v + dv and in a particular rovibrational state at a given nozzle temperature Tn can be described bywhere C is a normalization constant, v0 is the stream velocity, and α is the width of the velocity distribution. With eq , the reactivity of each state can be weighted according to its Boltzmann weight aswithHere, the factor gN in eq reflects the ortho/para ratio of hydrogen in the beam. For D2, gN = 2/3 (1/3) for even (odd) values of J, while for H2, gN = 1/4 (3/4) for even (odd) values of J. Z(Tn) is the partition function, kB is the Boltzmann constant, and Eν, is the energy of the rovibrational state characterized by the vibrational (ν) and rotational (J) quantum numbers. In eq , we take into account the rotational cooling of the H2 molecules due to the supersonic expansion by taking Trot = 0.8 × Tn.[30][30] Degeneracy-averaged reaction probabilities are computed from fully initial-state resolved reaction probabilities aswhere Pr(E, ν, J, m) is the fully initial-state-resolved reaction probability, with m being the magnetic rotational quantum number and E = 1/2mv2 being the translational energy. Molecular beam sticking probabilities can then be computed asAll parameters describing the molecular beams simulated in this work are listed in Table S3 in ref (14). A more exhaustive description of how molecular beam sticking probabilities can be computed can be found in ref (107). The set of initial rovibrational states taken into account in the QCT calculations is listed in Table .

Table 1

Rovibrational States Taken into Account, According to Their Boltzmann Weight, in Molecular Beam Simulations for the QCT Method for All H2 (D2) + Metal Systems

	(ν = 0) Jmax	(ν = 1) Jmax	(ν = 2) Jmax	(ν = 3) Jmax	(ν = 4) Jmax
QCT	30	30	30	30	30

Rovibrational State Populations of H2 and D2 Desorbing from Au(111)

The following expression is used to calculate state distributions of desorbing molecules[68]Here, Emax(ν, J) is the maximum kinetic energy to which the experiment was sensitive[68] in the sense that Pdeg(E, ν, J) could still be extracted reliably. These parameters have been obtained from Sven Kaufmann (private communications), who is one of the authors of ref (68), and the parameters are printed in Table S1. TS is the surface temperature. The Emax(ν, J) parameters for H2 (D2) + Au(111) are plotted in Figure S2 in ref (14). While Shuai et al.[68] integrated eq up to 5 eV, we opt to integrate only until Emax(ν, J) since the error function expressions derived in ref (68) are only reliable up to Emax(ν, J) and can yield sticking probabilities substantially bigger than 1 for high translational energies. We integrate eq by taking a right Riemann sum with ΔE = 0.2 meV. The N(ν, J) populations are normalized to the total ν = 0 population as was done in previous work.[14] The ratios of populations we calculate are solely based on the rovibrational states shown in Figure , i.e., we only go up to J = 7 for H2 and J = 9 for D2 as was done by Shuai et al.[68]

Figure 8

Rovibrational state populations of H2 and D2 desorbing from Au(111) are shown versus the rotational energy. Experimental results are shown in black,[68] and theoretical results are shown for MS-PBEl-rVV10 (red), MS-B86bl-rVV10 (blue), and MS-RPBEl-rVV10 (green) DF. The straight lines represent Boltzmann distributions for the surface temperature of the experiment.

E1/2(ν, J) Parameters

In our previous work, we listed four possible methods to obtain E1/2(ν, J) parameters, which can be used to compare to experimental E0(ν, J) parameters.[14] In this paper, we only use method B2 to compare calculated E1/2(ν, J) parameters to measured E0(ν, J) parameters for the H2 (D2) + Au(111) system. All four methods are discussed in the Supporting Information in ref[14][14], and we will only briefly discuss method B2 here.

When no measured sticking probabilities are available for the system of interest, one may choose to normalize the extracted reaction probabilities with reference to theory.[68,71] In method B1, theory is compared to experiment by extracting E1/2(ν, J) parameters usingIn other words, the E1/2(ν, J) parameter is the energy at which the degeneracy-averaged reaction probability is equal to half the saturation value, which is taken equal to the reaction probability at the maximum kinetic energy to which the experiment was sensitive.

However, for H2 (D2) + Au(111), the Emax(ν, J) parameters are not large enough to reliably extract E1/2(ν, J) parameters.[14] In method B2, the measured E0(ν, J) and Wν, values are therefore used to determine the Aν,B2 value at which the experimental reaction probability saturates according to the error function fit of the (ν, J) rovibrational state.[14,68] Effectively, in method B2, we take the Aν,B1 value and scale it accordingly[14]

Computational Details

A user-modified version 5.4.4 of the Vienna Ab initio Simulation Package[109−112] (VASP) has been used for all plane-wave periodic DFT electronic structure calculations. The modification of the computer package concerns the implementation of the mGGA DFs developed in this work.[24] In all calculations, the standard projector augmented wave (PAW) potentials[113] are used. We use the rVV10[47] nonlocal correlation functional as implemented in VASP,[37] which is based on the vdW-DF1[93,114,115] implementation by Klimeš et al.[116]

All calculations are carried out using a plane-wave cutoff energy of 600 eV together with smearing of 0.2 eV using the Methfessel–Paxton method of order 1. All slabs consist of six layers, of which the bottom two layers are fixed at their ideal bulk interlayer distance. A 2 × 2 supercell is used for calculations of the PESs with a vacuum distance of 16 Å and a (11 × 11 × 1) Γ-centered k-point grid. Lattice constants have been calculated using four-atom bulk unit cells and a (28 × 28 × 28) Monckhorst–Pack k-point grid, while slab relaxations were carried out using a (32 × 32 × 32) Γ-centered k-point grid together with a 1 × 1 supercell. For the molecule-metal surface calculations, a convergence parameter of 10–6 eV was used, and for the bulk lattice calculations, slab relaxations were used, and for the metal-atom calculations, a convergence parameter of 10–7 eV was used for the energy.

Results and Discussion

Metal Properties

Table shows the calculated lattice constants compared to zero-point energy-corrected experimental results[117] for the three MS mGGA-rVV10 DFs tested in this work as well as the original three MS mGGA DFs. For the four metals investigated here, we calculate lattice constants that are smaller than the zero-point energy-corrected experimental results, although the agreement with experiment[117] is still reasonable. The underestimation of the experimental lattice constants for the three DFs developed here is, on average, comparable to the somewhat overestimation of the lattice constants for SRP DFs designed for the reaction of H2 (D2) on transition metals at the GGA level that include nonlocal correlation.[14]

Table 2

Calculated Lattice Constants (in Å) and Relative Differences (in %) with Zero-Point Energy-Corrected Experimental Results[117]

	Cu		Ag		Au		Pt
	Å	%	Å	%	Å	%	Å	%
exp.[117]	3.596		4.062		4.062		3.913
MS-PBEl[24]	3.580	–0.4	4.090	0.7	4.084	0.5	3.906	–0.2
MS-PBEl-rVV10	3.514	–2.2	4.003	–1.4	4.034	–0.7	3.879	–0.9
MS-B86bl[24]	3.583	–0.4	4.092	0.7	4.087	0.6	3.908	–0.1
MS-B86bl-rVV10	3.518	–2.2	4.004	–1.4	4.036	–0.6	3.881	–0.8
MS-RPBEl[24]	3.590	–0.2	4.099	0.9	4.092	0.7	3.912	0
MS-RPBEl-rVV10	3.524	–2.0	4.008	–1.3	4.040	–0.5	3.884	–0.7
B86SRP68-DF2[14]	3.639	1.2	4.150	2.2	4.166	2.6	3.979	1.7
PBEα57-DF2[13]	3.656[14]	1.7	4.176[14]	2.8	4.198[14]	3.3	4.016[13]	2.6
PBEsol[27]	3.570[117]	–0.7	4.058[117]	–0.1	4.081[117]	0.5	3.932[117]	0.5

Table shows the interlayer contractions for the top two layers (in %) for Cu(111), Ag(111), Au(111), and Pt(111). When combining our three MS mGGA DFs[24] with rVV10[47] nonlocal correlation, we find that the relaxed six-layer slabs tend to expand somewhat, in contrast to the results obtained when not using nonlocal correlation.[14,24] The description of the relaxed slabs is not as good as obtained with previously developed SRP DFs,[14] and with our mGGA DFs not using nonlocal correlation[24] (see Table ).

Table 3

Relaxations of the Interlayer Distance of the Top Two Layers Relative to the Bulk Interlayer Distance in %

	Cu (%)	Ag (%)	Au (%)	Pt (%)
exp.	–1.0,[118,119] −0.7[120]	–2.5,[121] −0.5[122]	1.5[123]	1.1[124]
MS-PBEl[24]	–1.0	–0.4	1.0	1.0
MS-PBEl-rVV10	1.5	2.3	3.5	2.4
MS-B86bl[24]	–1.0	–0.5	1.0	1.0
MS-B86bl-rVV10	1.4	1.4	3.5	2.3
MS-RPBEl[24]	–1.6	–0.5	1.2	1.1
MS-RPBEl-rVV10	1.6	2.4	3.5	2.4
B86SRP68-DF2[14]	–0.4	–0.1	1.3	1.2
PBEα57-DF2[13]	–0.4[14]	0.0[14]	1.5[14]	0.8[13]
vdW-DF2[94]	0.0	0.5	2.1	1.5

The three original MS mGGA DFs[24] were developed to avoid having to compromise between a good description of the metal and a good description of the molecule–surface interaction.[16] It is therefore a somewhat disappointing result that when our three mGGA DFs are combined with nonlocal rVV10 correlation[47] this comes at the cost of a somewhat less good description of the metal. Tuning the b parameter in the implementation of rVV10 nonlocal correlation,[47] which modulates the repulsive part of the vdW description,[47] to obtain lattice constants closer to experiment unfortunately has the effect of removing the vdW wells in the PESs we calculate.

Including nonlocal correlation in a DF has a tendency to yield smaller lattice constants compared to DFs that do not include nonlocal correlation.[14,117] Our original MS mGGA DFs yield calculated lattice constants that are highly accurate.[24] Therefore, combining them with nonlocal correlation, which tends to shrink the lattice constants, leads to too small calculated lattice constants.

Similarly. we observe that the interlayer distance between the top two layers of the relaxed six-layer slabs tends to expand somewhat when using rVV10[47] nonlocal correlation (see Table ). When not using nonlocal correlation, our three MS mGGA DFs produced interlayer distances between the top layers that were in line with experimental results.[14,24] When using rVV10[47] nonlocal correlation together with our MS mGGA DFs, our calculated interlayer distances of the top layer are still reasonable, although not as good as those obtained with GGA-based SRP DFs that use vdW-DF1[93] or vdW-DF2[94] nonlocal correlation (see Table 2 in ref (14)). We speculate that the more accurate interlayer distance calculated when using vdW-DF1[93] or vdW-DF2[94] nonlocal correlation is due to the way in which the correlation part of the full exchange–correlation functional is constructed. In the case of vdW-DF1[93] or vdW-DF2,[94] the nonlocal correlation part is combined only with fully local LDA correlation. In the case of the MS mGGA-rVV10 functionals that we test here, the NLC rVV10 functional is combined with semilocal correlation instead (see eq ). For calculating lattice constants and interlayer spacings of metals, it might be better to combine the MS meta-GGA exchange functionals we investigate with correlation functionals based on LDA correlation and a nonlocal vdW-DF1 or vdW-DF2 nonlocal correlation functional.

Static PES Properties

Figure shows vdW wells for H2 in parallel (ϕ = 0°, θ = 90°) and perpendicular (θ = 0°) orientations above a top site for Cu(111) (a), Ag(111) (b), Au(111) (c), and Pt(111) (d). All vdW well geometries and well depths computed by us are tabulated in Table , also comparing with experimental results that have been reported for H2 + Cu(111),[88,89] H2 + Ag(111),[83] H2 + Au(111),[89] and H2 + Pt(111).[80,92] Note that we use the same b value (b = 6.3) for the three DFs that use rVV10[47] nonlocal correlation. As noted in our previous work,[14] for Cu(111), the experimental well depths are in good agreement. However, the position reported by Harten et al.[89] is somewhat closer to the surface. Ambiguities in the level assignments in the study of Andersson et al.[87] are the most likely reason for the vdW well being reported somewhat closer to the surface compared to the later measurements.[88] Andersson and Persson[88] noted that their derived PES is also consistent with the earlier measurements.[87] As mentioned in our previous work,[14] we suspect that reported vdW wells for H2 + Ag(111)[83] and H2 + Au(111)[89] might possibly be too close to the surface.[87]

Figure 2

vdW wells for H2 + Cu(111) (a), Ag(111) (b), Au(111) (c), and Pt(111) (d). Solid lines represent a parallel orientation of H2 (θ = 90°, ϕ = 0°), and dotted lines a perpendicular orientation (θ = 0°), both above a top site. Experimental results are shown in black for H2 + Cu(111),[88] H2 + Ag(111),[83] H2 + Au(111),[89] and H2 + Pt(111).[80,92] In (b) and (c), the horizontal solid lines correspond to the experimental well depths. In (d), the dashed line corresponds to the result of Poelsema et al.[92] and the solid line corresponds to the result of Cowin et al.[80] Results for five DFs are shown: MS-PBEl-rVV10 (red), MS-B86bl-rVV10 (green), MS-RPBEl-rVV10 (blue), vdW-DF1[93] (magenta), and vdW-DF2[94] (light blue).

Table 4

VdW Well Depths and Positions for Cu(111), Ag(111), Au(111), and Pt(111) for H2 in Parallel Orientation (ϕ = 0°, θ = 90°) above a Top Site

Cu(111)	Z (Å)	EvdW (meV)
exp.	3.51,[88] 2.71[89]	29.5,[88] 22.2[89]
MS-PBEl-rVV10	3.66	33.1
MS-B86bl-rVV10	3.99	18.2
MS-RPBEl-rVV10	4.05	23.7
vdW-DF1[93]	3.77	52.4
vdW-DF2[94]	3.58	39.0
B86SRP68-DF2[14]	3.74	34.3

The MS-PBEl-rVV10 DF performs best with respect to the vdW well interaction for all systems investigated. Highly accurate vdW well depths are obtained for both the highly activated systems and the nonactivated H2 + Pt(111) system with this functional (see Figure and Table ). The agreement with the position of the minimum is also good for the system for which this is well known, i.e., H2 + Cu(111). The agreement with the experimental well depth obtained with the MS-B86bl-rVV10 and MS-RPBEl-rVV10 DFs is reasonable. This agreement is not as good as obtained with the MS-PBEl-rVV10 DF, but the MS-RPBEl-rVV10 results agree better with the experimental results for the well depths for H2 + Ag(111) and Au(111) than the results previously obtained with the vdW-DF1 functional (see Table and Figure b,c). As discussed above, optimization of the b parameter to better reproduce the well depth obtained with the MS-B86bl-rVV10 and MS-RPBEl-rVV10 DFs would result in unacceptably small lattice constants.

When comparing the results, it is clear that the MS-PBEl-rVV10 DF yields a better description of the H2-metal vdW wells investigated here than the vdW-DF1[93] and vdW-DF2[94] DFs, which is consistent with earlier work[47] on the binding energies of a subset of molecular configurations of the S22 dataset[103] and the argon dimer.[47] However, we note that the previously tested[14] B86SRP68-DF2 DF (which performed better than the vdW-DF1[93] and vdW-DF2[94] DFs tested in ref (14)) shows a performance that is comparable to that of the MS-PBEl-rVV10 DF (Table ). We also note that the polarizability obtained for the H2 molecule parallel and perpendicular to its molecular axis is similar for the MS-PBEl-rVV10 and vdW-DF2[94] DFs.

In principle, the b parameter in the rVV10[47] nonlocal correlation functional could be tuned to match experimental observations of the vdW geometries in future work. However, decreasing the b parameter to obtain a vdW well geometry more in line with experiment would also lead to further decreased lattice constants, thereby further worsening the agreement with experiment, and it would lead to lower dissociation barriers.

Tables –8 show barrier heights and geometries for H2 + Cu(111), Ag(111), Au(111), and Pt(111), respectively. For the activated systems, the lateness of the barriers (values of r at which the barriers occur) is not influenced by the use of rVV10[47] nonlocal correlation. However, for the bridge sites, the barrier geometries do move to slightly higher Z values. Adding rVV10[47] nonlocal correlation to our original three MS mGGA DFs yields barrier heights that are consistently lower by roughly 0.15–0.2 eV for the highly activated systems. For the barrier heights obtained with the current best SRP DFs, we refer the reader to our previous work.[14]

Table 5

Barrier Heights for H2 Reacting on Cu(111)a

	bridge			t2b			fcc
	Eb	rb	Zb	Eb	rb	Zb	Eb	rb	Zb
MS-PBEl[24]	0.629	1.002	1.198	0.847	1.350	1.390	0.988	1.339	1.267
MS-PBEl-rVV10	0.459	0.985	1.240	0.665	1.328	1.400	0.815	1.331	1.285
MS-B86bl[24]	0.683	0.997	1.205	0.895	1.351	1.391	1.048	1.343	1.267
MS-B86bl-rVV10	0.513	0.982	1.247	0.714	1.329	1.401	0.865	1.333	1.285
MS-RPBEl[24]	0.721	1.006	1.201	0.930	1.354	1.392	1.086	1.346	1.270
MS-RPBEl-rVV10	0.549	0.985	1.247	0.747	1.329	1.403	0.899	1.334	1.286

For the bridge, t2b, and hcp sites, ϕ = 0° and θ = 90°. Barrier heights are in eV, and the barrier positions are in Å.

Table 8

Barrier Heights for H2 Reacting on Pt(111)a

	t2b early			t2b late			bridge			t2h early			t2h late
	Eb	rb	Zb	Eb	rb	Zb	Eb	rb	Zb	Eb	rb	Zb	Eb	rb	Zb
MS-PBEl[24]	0.145	0.766	2.205	-0.035	1.096	1.529	0.616	0.838	1.599	0.339	0.800	1.840
MS-PBEl-rVV10	0.008	0.763	2.326	-0.211	1.087	1.538	0.445	0.837	1.634	0.180	0.828	1.809	0.217	1.195	1.525
MS-B86bl[24]	0.194	0.768	2.189	0.016	1.085	1.534	0.667	0.839	1.602	0.392	0.802	1.836
MS-B86bl-rVV10	0.056	0.763	2.313				0.493	0.836	1.633	0.235	0.820	1.846	0.263	1.205	1.525
MS-RPBEl-rVV10	0.071	0.769	2.230				0.521	0.841	1.624	0.261	0.830	1.805	0.319	1.211	1.526

For the bridge and t2b sites, ϕ = 0° and θ = 90° for the t2h site ϕ = 120°. Barrier heights are in eV, and the barrier positions are in Å.

For the nonactivated H2 + Pt(111) system, we also find that using rVV10[47] nonlocal correlation leads to lower barriers, by about 0.15 eV. However, the picture is more complex since only three DFs show a double-barrier structure for the t2b site, namely, the MS-PBEl,[24] MS-B86bl,[24] and MS-PBEl-rVV10 DFs. The DFs without nonlocal correlation do not show a double-barrier structure for the t2h site, while the DFs that do use rVV10[47] nonlocal correlation do.

Note that observations on the vdW well depths and minimum positions extracted from RMSA[79] or CMSA[90,91] experiments usually represent averages taken over the sites in the surface unit cell. Checking for the site dependence of the vdW interaction in H2 + Cu(111), as found by Lee et al.,[125] we see essentially no dependence of the vdW interaction on the site within the unit cell (see Figure S4, which presents results for impact on three different sites obtained with the MS-PBEl-rVV10 DF). The site dependence found for the other systems and DFs treated here is similar to the results shown in Figure S4 in that it is very small.

Molecular Beam Sticking

Molecular Beam Sticking in H2 (D2) + Cu(111)

Molecular beam sticking probabilities for H2 (D2) + Cu(111) for six sets of molecular beam experiments are shown in Figure for the three MS mGGA-rVV10 DFs tested in this work and for the MS-B86bl and PBE DFs.[24] The parameters describing the molecular beam experiments[28−30] are tabulated in Table S3 in ref (14). Adding rVV10[47] nonlocal correlation to our three original mGGA DFs leads to higher sticking probabilities that are too high compared to experiment, as could be expected from its effect on the barrier heights (see Table ). MS-PBEl-rVV10, MS-B86bl-rVV10, and MS-RPBEl-rVV10 all overestimate the sticking probability and are not chemically accurate for this system. Given that the orignal three MS mGGA DFs were all chemically accurate for this system,[24] this is a somewhat disappointing result.

Figure 3

Molecular beam sticking probabilities for H2 and D2 reacting on Cu(111) for six sets of molecular beam experiments, as computed using the QCT method with the MS-B86bl[24] (purple), MS-PBEl-rVV10 (magenta), MS-B86bl-rVV10 (red), MS-RPBEl-rVV10 (blue), and PBE (green) DFs. Experimental results are shown in black.[28−30]

Molecular Beam Sticking in D2 + Ag(111)

Figure shows the sticking probabilities computed from simulations of molecular beams of D2 reacting on Ag(111) in comparison to experimental results.[31] Cottrell et al.[31] have reported molecular beam parameters that are symmetric with respect to the average collision energy. We consider these symmetric molecular beam parameters to be somewhat unphysical, as discussed in previous work from our group.[63] Therefore, we opted to use the molecular beam parameters of pure D2 reacting on Cu(111) reported by Auerbach and co-workers,[28] which likewise describe beams that are narrow in translational energy, in our simulations.

Figure 4

Molecular beam sticking probability as a function of the average incidence energy for D2 reacting on Ag(111). Experiment is shown in black.[31] QCT results are shown in blue for the following DFs: MS-PBEl-rVV10 (a), MS-B86bl-rVV10 (b), and MS-RPBEl-rVV10 (c). The values next to each data point denote the shift along the translational energy axis from the computed reaction probability to the interpolated experimental reaction probability in kJ/mol.

Mean absolute deviation (MAD) values are computed by calculating the mean distance along the incidence energy axis from the calculated sticking probability to the cubic spline interpolated experimental results. We consider DFs that yield a MAD value smaller than 1 kcal/mol (4.2 kJ/mol) to be chemically accurate.[126]Figure shows that all three DFs can be considered chemically accurate for this system, with the MS-PBEl-rVV10 DF performing best with a MAD value of 1.0 kJ/mol and the MS-RPBE-rVV10 DF performing worst with a still good MAD value of 2.0 kJ/mol. Here, we note that the distance between the computed and the measured S0 tends to increase with increasing translational energy.

This is the first time that we achieve a chemically accurate description of the D2 + Ag(111) system. GGA-based DFs with and without nonlocal correlation as well as the three original MS mGGA DFs have not been able to yield a chemically accurate description of this system.[14,24,63] The improved description of the sticking probability for this system is strictly due to the lowering of the barrier to reaction. Barrier geometries of the MS mGGA DFs that use nonlocal rVV10[47] correlation are very similar to the barrier geometries of the original MS mGGA DFs (see Table ).

Table 6

Barrier Heights for H2 Reacting on Ag(111)a

	bridge			t2b			fcc
	Eb	rb	Zb	Eb	rb	Zb	Eb	rb	Zb
MS-PBEl[24]	1.288	1.230	1.116	1.534	1.508	1.493	1.601	1.556	1.315
MS-PBEl-rVV10	1.082	1.224	1.157	1.328	1.486	1.506	1.392	1.553	1.345
MS-B86bl[24]	1.342	1.224	1.115	1.585	1.513	1.495	1.652	1.566	1.323
MS-B86bl-rVV10	1.134	1.223	1.159	1.376	1.488	1.507	1.442	1.560	1.348
MS-RPBEl-rVV10	1.171	1.226	1.161	1.410	1.489	1.508	1.479	1.560	1.349

For the bridge, t2b, and hcp sites, ϕ = 0° and θ = 90°. Barrier heights are in eV, and the barrier positions are in Å.

As we discussed in previous work from our group,[14,63] assessing the quality of the theoretical description of this system is difficult due to the lack of well-defined molecular beam parameters.[63] Additional experiments would allow us to improve the description of this system.[14]

Molecular Beam Sticking in D2 + Pt(111)

Figure shows calculations on D2 + Pt(111) for two sets of molecular beam experiments.[65,66] Note that this is a nonactivated system of which the original MS mGGA DFs gave a rather poor description.[14] Here, we find that the MS-B86bl-rVV10 DF (Figure b,e) yields the best results for both experiments, with MAD values of 2.7 and 2.0 kJ/mol, respectively, for the experiments of Luntz et al.[65] and Cao et al.[66] This may be compared to the MAD values of 1.1 kJ/mol for the experiment of Luntz et al.[65] and of 1.9 kJ/mol for the experiment of Cao et al.[66] that were obtained with the PBEα57-DF2 DF (Table ).[13]

Figure 5

Molecular beam sticking probabilities for D2 reacting on Pt(111) for the MS-PBEl-rVV10 (a, c), MS-B86bl-rVV10 (b, e), and MS-RPBEl-rVV10 (c, f) DFs. Experimental results are shown in red,[65,66] and QCT results in blue. The values next to each data point denote the shift along the translational energy axis from the computed reaction probability to the interpolated experimental reaction probability.

In general, the three MS mGGA-rVV10 DFs treated here are either in good agreement with experiment for the lower translational energies (MS-PBEl-rVV10) or for the higher translational energies (MS-B86bl-rVV10 and MS-RPBE-rVV10). The reason for this is that the MS-PBEl-rVV10 DF is the DF yielding the lowest early t2b barrier to reaction (see Table ), which allows it to describe the experiment correctly at the lowest translational energies. The other two mGGA-rVV10 DFs exhibit a higher early t2b barrier, leading to a worse description of the experiments[65,66] at low translational energies. Overall, the slope of the calculated sticking probability curve of the MS-PBEl-rVV10 DF is too steep, just right for the MS-B86bl-rVV10 DF, and somewhat too gentle for the MS-RPBEl-rVV10 DF.

Previous work from our group has indicated that the experiments of Luntz et al.[65] and Cao et al.[66] are in good agreement with each other for the lower incidence energies but somewhat diverge for the higher incidence energies.[127] The possible causes for this divergence are discussed in ref (127), where it was remarked that at high incidence energies, the reaction probabilities of Cao et al.[66] are most likely somewhat underestimated compared to the results of Luntz et al.[65] Note that a small increase of reactivity at the higher translational energies for the experiments of Cao et al.[66] could improve the agreement with experiment for the MS-B86bl-rVV10 DF. However, this system is still best described with the GGA-based SRP DF that was specifically designed for this system.[13,127]

Associative Desorption

Initial-State Resolved Reaction Probabilities in H2 (D2) + Ag(111)

Figure shows degeneracy-averaged initial-state resolved reaction probabilities for H2 and D2 reacting on Ag(111). A comparison is made to reaction probabilities extracted from associative desorption experiments assuming detailed balance.[69,70] Note that the experimental degeneracy-averaged reaction probabilities were not normalized but simply assumed to saturate at 1, which makes it hard to make a comparison with experiment. The translational energy in Figure refers to the translational energy of the desorbing molecules, which is measured by time-of-flight techniques using resonance-enhanced multiphoton ionization (REMPI) to achieve state selection.[69,70] For our purposes, and considering initial-state-selected reaction, this energy may also be considered as the collision energy in the experiment on reaction that is related to the associative desorption experiment via detailed balance.

Figure 6

Initial-state-selected reaction probabilities Pdeg(E, ν, J) computed for H2 (D2) + Ag(111) using the MS-PBEl-rVV10 (blue), MS-B86bl-rVV10 (green), and MS-RPBEl-rVV10 (red) DFs as a function of translational energy are shown, comparing with values extracted from associative desorption experiments.[69,70] Results are shown for D2 (ν = 0, J = 2) (a), D2 (ν = 1, J = 2) (b), and H2 (ν = 0, J = 3) (c).

From Figure , it can be seen that the three MS mGGA-rVV10 DFs somewhat overestimate the degeneracy-averaged reaction probabilities for D2 for most energies (Figure a,b), but that the agreement with experiment is very good for H2 (Figure c). In previous work,[24] the MS-PBEl DF was shown to perform better than other GGA-based DFs mainly due to MS-PBEl exhibiting slightly earlier barriers. The barrier geometries of the three MS mGGA-rVV10 DFs we present here are very similar to the barrier geometries of the three original MS mGGA DFs.[24] Therefore, we can say safely that the increased reactivity obtained with the mGGA-rVV10 DFs developed here is due to their barriers to reaction being somewhat lower and not to a change in barrier geometry (see Table ).

E1/2(ν, J) Parameters Au(111)

Figure shows a comparison of measured[68]E0(ν, J) parameters to E1/2(ν, J) parameters calculated using method B2.[14]Table shows the accompanying MAD and mean signed deviations (MSD) values. We note that the experiment was performed at a surface temperature of 1063 K,[68] while the calculations have been performed using the BOSS model. Furthermore, incorporating surface motion in the dynamics calculations would lead to a broadening of the reaction probability curves.[21,50,51,54] In view of the procedure used to calculate E1/2(ν, J) parameters, an increase of reactivity at low translational energies has the potential to lower the E1/2(ν, J) parameters.[14] We also note that our calculations have been carried out employing an unreconstructed Au(111) surface. Mapping out a full PES of H2 interacting with reconstructed Au(111) is currently extremely hard to do if not impossible, due to the large unit cell size.[14] Earlier work in our group[128] has shown that dynamical barrier heights of reconstructed Au(111) are roughly 50 meV higher compared to unreconstructed Au(111), which would lead to slightly higher computed E1/2(ν, J) parameters.[14]

Figure 7

E1/2 (ν, J) parameters as a function of J obtained using method B2 for H2 and D2 reacting on Au(111). Experimental values are shown in black.[68] Red circles represent the MS-PBEl-rVV10 values, green circles represent the MS-B86bl-rVV10 values, and blue circles denote the MS-RPBEl-rVV10 values.

Table 9

Mean Absolute and Mean Signed Deviations for the Theoretical E1/2(ν, J) Parameters Compared to Experimental E0(ν, J) Values for Au(111)[68]

	MAD (eV) H2			MSD (eV) H2			MAD (eV) D2			MSD (eV) D2
Au(111)	total	ν = 0	ν = 1	total	ν = 0	ν = 1	total	ν = 0	ν = 1	total	ν = 0	ν = 1
MS-PBEl[24]	0.106	0.104	0.107	–0.106	–0.104	–0.107	0.092	0.084	0.100	–0.084	–0.056	–0.112
MS-PBEl-rVV10	0.029	0.038	0.018	–0.028	–0.038	–0.016	0.044	0.045	0.042	–0.044	–0.045	–0.042
MS-B86bl[24]	0.139	0.131	0.150	–0.139	–0.131	–0.150	0.112	0.100	0.127	–0.112	–0.100	–0.128
MS-B86b-rVV10	0.050	0.053	0.046	–0.050	–0.053	–0.046	0.061	0.059	0.065	–0.061	–0.059	–0.065
MS-RPBE-rVV10	0.062	0.061	0.062	–0.062	–0.061	–0.062	0.073	0.068	0.078	–0.073	–0.068	–0.078

Even though all three developed MS mGGA DFs overestimate the measured E0(ν, J) parameters, it is clear from Table that MS-PBEl-rVV10 achieves chemical accuracy here for H2, and is just 1 meV less than chemical accuracy (43 meV) for D2. The MAD values of all three newly developed DFs are similar to the MAD values of PBE (46 meV for H2 and 58 meV for D2).[14] Previously, we have found that the original mGGA DFs as well as various GGA-based SRP DFs that include nonlocal correlation overestimate the experimental E0(ν, J) parameters by roughly 0.1 eV.[14] Furthermore, all three developed mGGA DFs reproduce the J dependence of the E0(ν, J) parameters quite well. As discussed previously,[14] this suggests that the reactivities of the individual rovibrational states are well described relative to one another, as long as states are considered within the same vibrational level. Given the uncertainties involved in using method B2 to calculate E1/2(ν, J) parameters, we obtain excellent results using our three newly developed MS mGGA-rVV10 DFs.

Previously reported experiments implied that the recombination of H2 on Au(111) is coupled to the electronic degrees of freedom of the metal.[129−132] Currently, we cannot disentangle the effects of ehp excitation, surface motion, and surface reconstruction. In our previous work,[14] we discussed how a combined analysis of a molecular beam dissociative chemisorption experiment on a reasonably cold surface (if available) and calculations on a reconstructed Au(111) surface, together with the associative desorption experiment of Shuai et al.,[68] could in principle be used to obtain a fingerprint of ehp excitation. Additionally, if a molecular beam dissociative chemisorption experiment were to become available, this would allow us to assess if the absolute reactivity computed with the new mGGA-rVV10 DFs and shown here is accurate.[14]

Rovibrational State Populations of H2 and D2 Desorbing from Au(111)

Rovibrational state populations for H2 and D2 desorbing from Au(111) are shown in Figure . Here, we plot ln[N/gN(2J + 1)] versus the rotational energy, with N being the total population for each (ν, J) state and gN(2J + 1) being the statistical weight for rotational level J.[68] In such a plot, a Boltzmann distribution will appear as a straight line.[68] Shuai et al.[68] have used an upper integration limit of 5 eV. Since the error function fits of the experiment are only reliable below Emax(ν, J), we opt to use Emax(ν, J) as the upper integration limit, as we did in previous work.[14] Note that we use the same normalization procedure as in our previous work.[14] The solid line represents Boltzmann distribution at the surface temperature of 1063 K used in the experiment.[68]

For molecules in the ground state, it can be seen that the rotationally excited molecules lie above the line set by the Boltzmann distributions. The experimental results lie on a gentler slope than the Boltzmann distributions, indicating that rotationally excited molecules are more likely to adsorb.[68] Similarly, the results for vibrationally excited molecules lie on a line with a gentler slope than shown by the Boltzmann distributions. Additionally, the results for vibrationally excited molecules lie substantially above the line of Boltzmann distributions, thereby indicating that vibrationally excited molecules are more likely to adsorb.[68]

Table shows the ν/ν = 1:0 ratio of desorbing molecules; these ratios are calculated using the same rovibrational states as shown in Figure . Note that the difference between the experimental values shown in Table and those reported by Shuai et al.[68] is due to using Emax(ν, J) as the upper integration limit.

Table 10

Ratio of ν = 1 : ν = 0 Molecules Desorbing from Au(111) as Measured in Experiments[68] and Computed with the MS-PBEl,[24] MS-PBEl-rVV10, MS-B86bl,[24] MS-B86bl-rVV10, and MS-RPBEl-rVV10 DFs

	H2	D2
exp.[68]	0.552	0.424
MS-PBEl[24]	0.178	0.387
MS-PBEl-rVV10	0.175	0.377
MS-B86bl[24]	0.176	0.379
MS-B86bl-rVV10	0.193	0.365
MS-RPBEl-rVV10	0.180	0.366

From Figure , it is clear that the differences between all DFs shown is minimal and that the agreement between theory and experiment is best for D2. As was already reported by Shuai et al.,[68] the theoretical ratios computed with different DFs for H2 are much lower than the experimental ratio. In our previous work,[14] we speculated that this difference might be resolved by including surface motion in our dynamics calculations because the experimental time-of-flight distributions are much broader compared to the theoretical ones.[68]

It is clear that adding nonlocal correlation to the MS mGGA DFs has little effect on the ν/ν = 1:0 ratio of desorbing molecules. GGA-based DFs yielded slightly better ratios for D2 desorbing from Au(111).[14] However, also these DFs predicted desorption ratios for H2 that were much too low.

The fact that mGGA-based DFs yield somewhat lower ν/ν = 1:0 ratios than the GGA-based DFs[14] can be explained by the barriers to reaction predicted by the mGGA DFs being somewhat earlier. This allows the ν = 0 population too grow somewhat relative to the ν = 1 population, which would lower the ν/ν = 1:0 ratios.

Transferability

Previous work from our group has shown that semilocal DFs designed for the reaction of H2 and D2 dissociating on transition metals may be transferable between different crystal faces of the same metal,[43,44] but until quite recently transferability between systems in which H2 interacts with different metals had not yet been observed.[63,133] Recently, we have reported that nonlocal correlation is a key ingredient in obtaining SRP DFs for the reaction of H2 and D2 on transition metals that show this type of transferability, by showing that a DF that we designed to describe the activated reaction of H2 + Cu(111) can also describe the reaction of D2 + Pt(111) and vice versa.[14] Earlier, transferability of SRP DFs between systems in which a molecule interacts with surfaces of different metals has only been reported for CH4 reaction on Ni(111)[45] and CH4 reacting on Pt(111).[46]

In our calculations, we employ the BOSS model and thus neglect any surface temperature effects, and it is known that the BOSS model works well for activated H2 dissociation on cold metals.[19,49−51,53] Given that associative desorption experiments necessitate high surface temperatures,[68,71] it is difficult to assess the quality of the DFs we developed here for the H2 (D2) + Au(111) system, due to the absence of molecular beam sticking experiments for this system.

Here, we show that MS-PBEl-rVV10, MS-B86bl-rVV10, and MS-RPBEl-rVV10 DFs can describe molecular beam sticking experiments on D2 + Ag(111) to within chemical accuracy (see Figure ) and that the MS-B86bl-rVV10 DF can also describe the D2 + Pt(111) molecular beam sticking experiments of Luntz et al.[65] and Cao et al.[66] to within chemical accuracy (see Figure ). In the case of the H2 (D2) + Au(111) system, the MS-PBEl-rVV10 DF yields very good results with respect to the calculated E1/2(ν, J) parameters. To the best of our knowledge, this is the first time that an observable of the reaction of H2 (D2) on Au(111) that requires dynamics calculations is described with chemical accuracy. However, uncertainties remain for this system with respect to the effects of surface temperature, surface reconstruction, and ehp excitation.[14]

Thus, there now exist two groups of transferable (SRP) DFs for the reaction of H2 (D2) with transition-metal surfaces. The first group consists of GGA-based SRP DFs that use vdW-DF2[94] nonlocal correlation (B86SRP68-DF2[14] and PBEα57-DF2[13]), which can describe the H2 (D2) + Cu(111) and D2 + Pt(111) reactions to within chemical accuracy. The second group consists of the MS mGGA DFs that use rVV10[47] nonlocal correlation developed here, which can describe the D2 + Ag(111) and D2 + Pt(111) systems with chemical accuracy. Of course, there is also the nonconclusive evidence that suggests that the MS-PBEl-rVV10 DF can describe the associative desorption of H2 from Au(111) to within chemical accuracy.

At present, we cannot say which features of a PES are most important, apart from the lowest barrier to reaction. Experiments that probe different parts of a PES, like vibrationally or rotationally inelastic scattering, where the latter process depends on the anisotropy of the PES, are few and far between.[134,135] In general, we see that the MS mGGA-based DFs have somewhat earlier barriers for highly activated systems than the GGA-based SRP DFs, while for the nonactivated D2 + Pt(111) system, the barrier geometries of the MS mGGA-based DFs that include rVV10[47] nonlocal correlation are very similar to the barrier geometries of GGA-based SRP DFs that include nonlocal vdW-DF2[94] nonlocal correlation. At the moment, we cannot say which type of barrier geometry is more in line with reality.

If the suggested chemical accuracy in the description of H2 + Au(111) holds in contrast to experiment, then one could argue that the mGGA-based DFs that include rVV10[47] nonlocal correlation are an improvement over the previously developed GGA-based SRP DFs that include vdW-DF2[94] nonlocal correlation: in this case, the MS mGGA-rVV10-based DFs can describe three systems with chemical accuracy, compared to two systems for the GGA-based SRP DFs.[14] This would indicate that climbing Jacob’s ladder leads to a more universal description of the reaction of H2 on transition-metal surfaces.

Conclusions

We have combined our three previously developed MS-PBEl, MS-B86bl, and MS-RPBEl MS mGGA DFs with rVV10 nonlocal correlation to obtain the MS-PBEl-rVV10, MS-B86bl-rVV10, and MS-RPBEl-rVV10 DFs. We find that all three developed DFs can describe the molecular beam sticking experiments on dissociative chemisorption of D2 on Ag(111) with chemical accuracy. We also find that the MS-B86bl-rVV10 DF can describe two sets of molecular beam sticking experiments on dissociative chemisorption of D2 on Pt(111) with chemical accuracy. Additionally, by calculating E1/2(ν, J) parameters for the reaction of H2 on Au(111) and comparing these to experimental E0(ν, J) parameters for state-selective associative desorption, we obtain chemical accuracy with the MS-PBEl-rVV10 DF. Assessing the performance of the three developed MS mGGA-rVV10 DFs for the H2 (D2) + Au(111) system is however difficult due to the lack of well-characterized molecular beam sticking experiments of H2 (D2) on Au(111) and the lack of calculations that use a reconstructed Au(111) surface and incorporate ehp excitation.

Of the three developed MS mGGA-rVV10 DFs, MS-PBEl-rVV10 performs excellently for the known vdW well geometries. The MS-PBEl-rVV10 DF also maintains the improvements generally observed for mGGA-rVV10 DFs relative to GGA-vdW-DF2 DFs in this regard. The MS-B86bl-rVV10 and MS-RPBEl-rVV10 DFs yield vdW wells that are too shallow.

In comparison to state-selected experiments on associative desorption of H2 (D2) from Ag(111), we observe excellent agreement with experiment in the case of H2, for all three developed DFs. For H2, all three developed DFs show improvement over the three original MS mGGA DFs and over the best GGA-based SRP DFs that include vdW-DF2 nonlocal correlation. The associative desorption experiments on D2 desorbing from Ag(111) were less well described.

With respect to the molecular beam sticking probabilities of H2 (D2) + Cu(111), the three developed DFs yield sticking probabilities in line with the sticking probabilities predicted by the PBE DF, which are too high. This is in contrast to the highly accurate sticking probabilities obtained when using the original three MS mGGA DFs.

The three original MS mGGA DFs give a description of the metal that is comparable to that obtained with the PBEsol DF. Unfortunately, adding rVV10 nonlocal correlation comes at the cost of a worse description of the metal. In general, we see lattice constants that are smaller than the zero-point energy-corrected experimental results. However, in general, the underestimation of the calculated lattice constants is still smaller than the overestimation of calculated lattice constants obtained with the current best SRP DFs that include vdW-DF2 nonlocal correlation. The three developed MS mGGA-rVV10 DFs also predict interlayer distances between the top two layers that are too large compared to experimental observations.

The three MS mGGA DFs that have been combined in this work with rVV10 nonlocal correlation were not fitted to reproduce particular experiments, nor has the b parameter present in rVV10 been reoptimized. Our results show that, overall, ascending Jacob’s ladder from the GGA plus nonlocal correlation rung to the mGGA plus nonlocal correlation rung leads to somewhat more accurate results for dissociative chemisorption of H2 (D2) on noble metals, although the metals themselves are described less accurately, and the improvement does not hold for the well-studied H2 + Cu(111) system.

Table 7

Barrier Heights for H2 Reacting on Au(111)a

	bridge			t2b			fcc
	Eb	rb	Zb	Eb	rb	Zb	Eb	rb	Zb
MS-PBEl[24]	1.432	1.144	1.127	1.301	1.433	1.466	1.350	1.203	1.276
MS-PBEl-rVV10	1.251	1.148	1.159	1.139	1.425	1.475	1.172	1.216	1.299
MS-B86bl[24]	1.481	1.142	1.130	1.355	1.438	1.467	1.402	1.204	1.276
MS-B86bl-rVV10	1.302	1.147	1.162	1.192	1.427	1.476	1.224	1.216	1.299
MS-RPBEl-rVV10	1.336	1.147	1.163	1.226	1.436	1.476	1.258	1.219	1.302

For the bridge, t2b, and fcc sites, ϕ = 0° and θ = 90°. Barrier heights are in eV, and the barrier positions are in Å.