Specific Reaction Parameter Density Functional Based on the Meta-Generalized Gradient Approximation: Application to H2 + Cu(111) and H2 + Ag(111).

Specific reaction parameter density functionals (SRP-DFs), which can describe dissociative chemisorption reactions on metals to within chemical accuracy, have so far been based on exchange functionals within the generalized gradient approximation (GGA) and on GGA correlation functionals or van der Waals correlation functionals. These functionals are capable of describing the molecule-metal surface interaction accurately, but they suffer from the general GGA problem that this can be done only at the cost of a rather poor description of the metal. Here, we show that it is possible also to construct SRP-DFs for H2 dissociation on Cu(111) based on meta-GGA functionals, introducing three new functionals based on the "made-simple" (MS) concept. The exchange parts of the three functionals (MS-PBEl, MS-B86bl, and MS-RPBEl) are based on the expressions for the PBE, B86b, and RPBE exchange functionals. Quasi-classical trajectory (QCT) calculations performed with potential energy surfaces (PESs) obtained with the three MS functionals reproduce molecular beam experiments on H2, D2 + Cu(111) with chemical accuracy. Therefore, these three non-empirical functionals themselves are also capable of describing H2 dissociation on Cu(111) with chemical accuracy. Similarly, QCT calculations performed on the MS-PBEl and MS-B86bl PESs reproduced molecular beam and associative desorption experiments on D2, H2 + Ag(111) more accurately than was possible with the SRP48 density functional for H2 + Cu(111). Also, the three new MS functionals describe the Cu, Ag, Au, and Pt metals more accurately than the all-purpose Perdew-Burke-Ernzerhof (PBE) functional. The only disadvantage we noted of the new MS functionals is that, as found for the example of H2 + Cu(111), the reaction barrier height obtained by taking weighted averages of the MS-PBEl and MS-RPBEl functionals is tunable over a smaller range (9 kJ/mol) than possible with the standard GGA PBE and RPBE functionals (33 kJ/mol). As a result of this restricted tunability, it is not possible to construct an SRP-DF for H2 + Ag(111) on the basis of the three examined MS meta-GGA functionals.

Introduction

Dissociative chemisorption reactions often control the rate of heterogeneously catalyzed processes,[1,2] which are of large importance to the chemical industry.[3] Important examples include dissociative chemisorption of N2 in ammonia synthesis[4] and dissociation of methane in steam reforming.[5] Accurately simulating rate-controlling reactions is critical to the calculation of accurate rates of the overall catalyzed processes.[6]

The best computational method to obtain accurate results for dissociative chemisorption reactions is currently based on the specific reaction parameter (SRP) approach to density functional theory (DFT) (SRP-DFT). In this approach, the density functional is taken as the weighted average of two functionals, using a mixing parameter that is typically fitted to obtain agreement with an experiment on dissociative chemisorption for the specific system considered. This method has now been applied successfully to three H2–metal systems[7−9] and three CH4–metal systems,[10,11] in the sense that it was possible to describe the sticking probability as a function of incidence energy with chemical accuracy (to within better than 1 kcal/mol). (Note that the terms dissociation, dissociative chemisorption, and sticking are used interchangeably in this work).

So far, the SRP density functionals (SRP-DFs) that have been developed were based[7,8,12] on exchange–correlation functionals within the generalized gradient approximation (GGA),[13−15] or they were based[9,10] on GGA exchange functionals[14−16] and Lundqvist–Langreth van der Waals correlation functionals.[17,18] Unfortunately, GGA functionals are not good at giving a simultaneously accurate description of molecule–surface interaction energies (and therefore reaction barriers) and metal surfaces (surface energies and lattice constants).[19] Specifically, GGA functionals that are good at describing adsorption underestimate metal surface energies and overestimate metal lattice constants, while GGA functionals that are good at describing metals overestimate adsorption energies,[19] in spite of efforts to design GGA functionals[20] or nonseparable gradient approximation functionals[21] that perform equally well on both properties. The problem is often, but not always, exacerbated, in the sense that lattice constants are increased further,[22] if a Lundqvist–Langreth correlation functional describing the attractive van der Waals interaction[17,18] is used with a GGA exchange functional to arrive at an SRP-DF.[9−11] A description of the metal that is fair at best may pose a problem for an SRP-DF if it is to describe sticking over a large range of surface temperatures (Ts) or sticking (often measured at low Ts)[23,24] and associative desorption (often measured at high Ts).[7,23,24] The reason is that the reaction barrier may depend on the interlayer spacing in the top two metal layers[12,25,26] as well as the amplitude of motion of the metal atoms in the top layer,[27,28] both of which are properties of the metal and depend on Ts.

A specific advantage of meta-GGA functionals is that their additional dependence on the kinetic energy density τ allows one to distinguish between regions of electron densities describing single (covalent) bonds, metallic bonds, and weak bonds.[29] A particularly elegant way is to introduce a dimensionless parameter α that is a function of τ such that α = 0 corresponds to covalent bonding, α = 1 to metallic bonding, and α ≫ 1 to weak bonding.[29] This parameter has been employed in the construction of several much used meta-GGA functionals, such as the TPSS,[30] the revTPSS,[31] the SCAN,[32] and the mBEEF[33] functionals. More recently, meta-GGA exchange functionals have been constructed on the basis of a function f(α) such that the exchange functional can be interpolated between the α = 0 and α = 1 limits and extrapolated to the α = ∞ limit. Examples of such methods include the meta-GGA made simple (mGGA-MS)[34,35] and the meta-GGA made very simple (mGGA-MVS) functionals.[36] A good simultaneous description of adsorption energies and lattice constants[33,37,38] or more generally of energetics and structure[34−36,39,40] has now been reported by several groups using meta-GGA functionals.

Here, we test a new variant of the meta-GGA functional of the mGGA-MS type.[34,35] We explore the use of mGGA-MS functionals, in which exchange functionals are used with Perdew–Burke–Ernzerhof (PBE)-like[15] and RPBE-like[14] expressions and with a B86b-like expression.[41] We show that the three meta-GGA MS functionals provide a chemically accurate description of several molecular beam experiments for a benchmark system,[7,23,24,42−44] that is, dissociative chemisorption of H2 on Cu(111), while yielding a much more accurate description of the Cu lattice and of other metals than the previous SRP-DFs for this system. Our results imply that an SRP-DF can be constructed for H2 + Cu(111) on the basis of the two MS functions with PBE-like and RPBE-like expressions (MS-PBEl and MS-RPBEl), just like it was possible to do on the basis of the actual RPBE and PBE (or almost equivalently PW91) functionals.[7,12] The MS-PBEl functional also gives a slightly better description of experiments on sticking of D2 to Ag(111)[45] than achieved[46] with the previous SRP48 functional for H2 + Cu(111), but a chemically accurate description was not yet achieved for this system. We suggest that SRP-DFs for specific systems can be based on mixtures of the mGGA-PBE-MS and the mGGA-RPBE-MS functionals. However, reaction barriers are tunable over a range that is smaller than the range obtained with the original PBE and RPBE functionals. Therefore, it may be necessary to replace one of the two MS functionals with a more attractive or repulsive meta-GGA functional to obtain a meta-GGA SRP-DF for specific systems, such as H2 + Ag(111).

This paper is set up as follows. In Section , we describe the dynamical model used. Section describes the “made simple” mGGA functionals we used and Section how we performed DFT calculations with these functionals and the way in which the DFT data were interpolated to obtain global potential energy surfaces (PESs). Section describes the classical trajectory method we used and Section how we computed observables. In Section , we present and discuss our results. Section provides results concerning the description of the metal, Section discusses the PESs computed with the three new functionals, and in Section we present dynamics results, also comparing to experiments on H2 + Cu(111) and Ag(111). Section concludes this paper.

Method

Dynamical Model

The model used is the Born–Oppenheimer static surface model, assuming the surface atoms to be in their ideal lattice positions. In this model, the effect of electron–hole pair excitation, of surface phonon motion, and of Ts on reaction is neglected. In view of the low Ts employed in the molecular beam experiments on H2 + Cu(111),[23,47] D2 + Cu(111),[24] and D2 + Ag(111)[45] that we compare to, these approximations should be reasonable, as discussed for instance in ref (44). Note in particular that recent theoretical work shows only a modest influence of electron–hole pair excitation on sticking of H2 to Cu(111)[48] and Ag(111).[49] The six degrees of freedom in which the motion is explicitly modeled are the coordinates of H2. The center-of-mass position of H2 is described by its Cartesian coordinates X, Y, and Z, where Z is the distance to the surface. The orientation of the molecule is specified by the polar angle θ and the azimuthal angle ϕ, and r is the H–H distance. The coordinate system and the (111) surface of face-centered cubic (fcc) metals such as Cu and Ag are shown in Figure .

Figure 1

(a) Center of mass coordinate system used for the description of the H2 molecule relative to the (111) face of an fcc metal (with ABCA... stacking). (b) Surface unit cell and the sites considered for the (111) face, and the relationship with the coordinate system chosen for H2 relative to the (111) surface. The origin (X,Y,Z) = (0,0,0) of the center of mass coordinates is located in the surface plane at a top site. Polar and azimuthal angles θ and ϕ are chosen such that (θ = 90°, ϕ = 0°) corresponds to molecules parallel to the surface along the X (or equivalently U) direction. The hexagonal close packed (hcp) and fcc hollow sites correspond to metal atoms in layers 2 and 3.

Made Simple Meta-GGA Density Functionals

Aiming at a more flexible functional with the potential for simultaneously good description of lattice structure and surface adsorption, we extend beyond the GGA space used for the SRP functional to the meta-GGA functional space. In general, this can be done by making the density functional depend also on the kinetic energy density τ, in addition to the density and its gradient. We follow Sun et al.[34] in expressing the meta-GGA in the so-called “made-simple” (MS) form as an interpolation between the exchange part of two GGAs for two extreme scenarios: the uniform electron gas (UEG) (a limit that describes metallic bonding rather well) and a single-orbital system (as in covalent bonds). The exchange enhancement factor describing the increase of exchange relative to the UEG then readsFx1(p) and Fx0(p;c) are the gradient only dependent exchange-enhancement factors for the UEG and single-orbital cases, respectively. The exchange enhancement factor is used in GGA and meta-GGA functionals to obtain the exchange part of the exchange–correlation energy. The numerical parameter c will be discussed below. In eq , p = s2, where s is the reduced density gradient, which is proportional to the gradient of the electron density divided by n4/3, n being the electron density (the exact expression is given in ref (34)). The interpolation f(α) depends on the Kohn–Sham kinetic energy τ through the inhomogeneity parameter[34,35]Here, τW is the von Weizsäcker kinetic energy density, which is equal to the kinetic energy density associated with a single-orbital electron density.[38] Furthermore, τunif is the kinetic energy density of the UEG. The expressions of τW and τunif may be found in for instance ref (35). Crucial points are that for a slowly varying electron density as found in metals, α approaches 1 as τ ≈ τunif and τW ≪ τunif, whereas for a single-orbital electron density as found in covalent bonding α = 0, because τ = τW.[38] The MS functionals take advantage of this by defining a function f(α) that equals 0 for α = 1 and that equals 1 for α = 0[35]

Interpolation between the exchange enhancement factor FxMS(p,α) describing the UEG and a single orbital can then be enforced by taking Fx1(p) = FxMS(p,α = 1) from a GGA functional accurately describing metallic bonding, and Fx0(p;c) = FxMS(p,α = 0)from a GGA functional accurately describing single-orbital systems, as in covalent bonds.[29] We note that in the present work we have simply taken b equals to 1, as in ref (34).

In this work, for the functional describing metallic bonding, we use the following PBE-like,[15] RPBE-like,[14] or B86b-like[41] expressions

In eqs −4c, for κ, we use the value 0.804 used in the original expressions, which corresponds to imposing the Lieb–Oxford bound on the exchange–correlation energy in eqs and 4b.[14,15] A crucial point is that to make the functionals defined in eqs −4c describe metallic bonding we always take μ = 10/81, as appropriate for metallic bonding[50] and as opposed to the actual value used in the PBE and RPBE functionals. This way, with the PBE-like expression we recover the PBEsol exchange functional,[51] which was designed to accurately describe elastic properties of metals.

To obtain the corresponding functionals describing covalent bonding, Fx0, in line with the MS strategy, we replace μp everywhere by (μp + c) in eq , with c > 0. We then optimize c for each functional to reproduce the exchange energy exactly for the hydrogen atom by cancelling the spurious self-interaction present in the Hartree energy in this atom. The strategy of hydrogen self-interaction correction is adopted in several meta-GGAs, for instance in MS functionals,[35] and in the TPSS[30] and SCAN[32] functionals. The spurious interaction arises because in DFT an electron interacts with itself through the use of a classical expression for the Coulomb interaction of electron densities, and this anomalous self-interaction even occurs for the one-electron H-atom, for which it can be computed and corrected for (in the exchange–correlation energy) exactly.[52] We choose to follow this strategy and c can be determined straightforwardly by numerical quadrature over the analytical nonrelativistic hydrogen atom density. This can be expected to lead to a reasonable GGA Fx0(p;c) for further single-orbital densities in general, importantly suppressing otherwise significant self-interaction errors in for example covalent bonds. Tuning the functional form of the made simple exchange functionals depending on the inhomogeneity of the density thus allows for more accurate general-purpose functionals than possible at the pure GGA-level, in particular with the capability of describing interactions within the metal well where α ≈ 1 as well as the inhomogeneous scenario of covalent and surface bonds.

In eqs −4c, the following values of c were obtained by numerical integration: 0.1036 (eq ), 0.07671 (eq ), and 0.08809 (eq ). For the correlation functional, we used the variant of the PBE correlation (vPBEc) used in revTPSS,[31] as was also done in the MS functionals of refs (34) and (35). We call the three MS functionals described in this way MS-PBEl, MS-RPBEl, and MS-B86bl, where the “l” stands for “like” to emphasize that we use a different value of μ and a different correlation functional than in the original PBE, RPBE, and B86b expressions. We emphasize that the expressions 4a–4c are non-empirical, like the original PBE, B86b, and RPBE expressions, and that no empirical fitting was performed for the b and κ parameters, as was done in ref (35). Finally, note that the PBE-like expression (eq ) was used before with MS functionals,[34,35] but with different values of κ, and/or c, and/or b.

DFT Calculations and Representation of PESs

All PESs used here were constructed from self-consistent, periodic DFT calculations carried out with a user modified version of the VASP program,[53−56] using the three functionals described above in Section , as well as other functionals. All calculations used projected augmented wave pseudopotentials from the VASP database.[57] All calculations used a 2 × 2 surface unit cell, a (11 × 11 × 1) Γ-centered Monkhorst–Pack k-point grid, a plane wave cut-off of 600 eV, 6 metal layers, a smearing of 0.2 eV using the Methfessel–Paxton method of order 1, and a vacuum distance between the slabs of 16 Å. Additional details and convergence tests, including details of how the metal lattice was computed and on the interlayer relaxation of the slab, are presented in the Supporting Information.

To obtain all PESs described here, the DFT data obtained with a particular functional were interpolated with the corrugation reducing procedure (CRP).[58,59] In the CRP, the six-dimensional PES V6D is written as a sum of a 6D interpolation functional I6D and two 3D potentials V3D, describing the interaction of the H-atoms with the surfaceHere, the 6D interpolation functional I6D is easier to fit than the full 6D potential because its “corrugation” has been reduced by subtracting the two 3D atom–surface potentials.[58] A similar trick is used in the interpolation of the atom–surface potentials.[58] The details of how we interpolated the PESs are mostly the same as described in ref (60) for H2 + Ru(0001) and ref (61) for H2 + Cu(111); where these details deviate, this is described in Section 2 of the Supporting Information.

Quasi-Classical Trajectory Method

To compute observables, the quasi-classical trajectory (QCT) method[62] was used. The QCT method gives a very good description of initial-state resolved reaction probabilities for activated H2–metal surface systems in general[63,64] and for H2 + Cu(111) in particular (see Figure 5 of ref (65)). The QCT method would be expected to perform even better for sticking probabilities, simulating the results of molecular beam experiments, which are highly averaged quantities, involving averages over collision energy and H2 internal states distributions. Specifically, results for the very similar, but slightly less reactive H2 + Cu(211) system[66] show that QCT sticking probability curves accurately reproduce quantum dynamical sticking probabilities down to probabilities of 0.002 (Smeets, Füchsel, and Kroes, to be published). In all calculations, we model scattering at normal incidence. Reaction and scattering probabilities are calculated by counting how many trajectories result in a particular outcome and dividing by the total number of trajectories. The H2 is initially placed at a distance from the surface where it no longer interacts with the surface (Z = 8 Å). It is considered to have dissociated when r > 2.25 Å, and the molecule is considered to have scattered when Z becomes >8 Å and the molecule is moving away from the surface. The initial conditions are simulated using standard Monte Carlo methods as described in ref (61). To obtain accurate statistics, for each incidence condition at least 100 000 trajectories were propagated. To integrate the equations of motion, the method by Stoer and Burlisch[67] was used. In the trajectories, the maximum propagation time was taken as 22 ps.

Computation of Observables

The initial (clean surface) sticking probability S0 can be computed from initial-state resolved reaction probabilities with appropriate averaging over the velocity distribution of the molecular beam and the rovibrational state distribution of the molecules in the beam.[68,69] The way this is done is described in the Supporting Information.

Results and Discussion

Description of the Metal

Equilibrium lattice constants computed with the three MS functionals are compared with the experimental values and the values computed with other functionals in Table , for the noble metals Cu, Ag, Au, and Pt. To facilitate the comparison with theory, the experimental lattice constants were corrected by subtracting a contribution due to zero-point vibrations.[70] The comparison clearly shows the advantage of the MS functionals that we already anticipated: the mean signed deviations (MSDs) from the experiment computed with the MS-PBEl (0.008 Å) and MS-B86bl (0.009 Å) functionals are considerably lower than that obtained with the all-purpose functional PBE[15] (0.015 Å), with only a GGA specifically designed for the solid state (PBEsol[51]) performing better (MSD = 0.002 Å, Table ). All MS functionals perform much better for lattice constants than the RPBE[14] functional, which consistently overestimates lattice constants, with a mean absolute deviation (MAD) equal to the MSD = 0.127 Å (Table ).

Table 1

Equilibrium Lattice Constants, the MSD, and MAD with Respect to Experiment (All in Å) Computed with the MS Functionals in This Work Are Compared to the Experimental Values and the Values Computed with Other Density Functionals[70]a

metal	Expt.	MS-PBEl	MS-B86bl	MS-RPBEl	PBE[70]	PBEsol[70]	RPBE
Cu	3.596	3.585	3.583	3.590	3.632	3.570	3.68[40]
Ag	4.062	4.091	4.092	4.099	4.152	4.058	4.23[40]
Au	4.062	4.084	4.087	4.092	4.154	4.081	4.23[71]
Pt	3.913	3.906	3.908	3.912	3.985	3.932	4.00[40]
MSD		0.008	0.009	0.015	0.015	0.002	0.127
MAD		0.017	0.018	0.019	0.019	0.017	0.127

The experimental values were corrected for zero-point energy effects in ref (70).

The performance of the PBE and RPBE functionals is relevant for the description of dissociative chemisorption: in many cases, PBE (or the very similar[15] PW91[13] functional) overestimates and RPBE underestimates the reactivity,[7,60,72] and an SRP functional or in any case an improved functional can be constructed by taking a weighted average of the PBE and RPBE functionals.[7,8]Table suggests that such GGA functionals should yield too large lattice constants (0.017 Å ≤ MSD ≤ 0.127 Å). This arises from the need to achieve a good description of the molecule–surface interaction energy: in the construction of a GGA functional, this goes at the cost of a good description of the metal lattice.[19] Finally, the good performance of the three MS functionals shown here is consistent with findings of earlier studies using a MS functional.[34]

Table shows the interlayer contractions (in %) computed for Cu(111) and Ag(111), for the interlayer distance between the first two top layers, also comparing with experiments. For Cu, especially the MS-PBEl and the MS-B86bl functionals yield good agreement with experiments, especially with the medium energy ion scattering experiments.[73] For Ag, all three MS functionals yield good agreement with the low-energy electron diffraction (LEED) experiments of Soares et al.[74] It is not clear to us what the source of discrepancy is between these experiments and the energy ion scattering experiments of Statiris et al.[75] However, we note that the LEED results of Soares et al. are in good agreement with results obtained recently[76] with the SCAN[32] functional and other functionals. It is important that a functional gives a good description of the interlayer contractions between the top two layers of a given surface, as this may have an important effect on the dissociation barrier height, as found for both H2 + Cu(111)[12,26] and Cu(100).[25]

Table 2

Relaxations of the Interlayer Lattice Spacing between the Top Two Layers Relative to the Bulk, in %, for Cu(111) and Ag(111)

metal	MS-PBEl	MS-B86bl	MS-RPBEl	Expt.
Cu	–1.0%	–1.0%	–1.6%	–1.0,[73] −0.7[77]
Ag	–0.4	–0.5	–0.5	–2.5,[75] −0.5[74]

Potential Energy Surfaces

Figure shows elbow plots of the MS-B86bl PES (i.e., plots of the dependence of the potential on r and Z for specific orientations and center-of-mass projections on the surface of H2) for four configurations in which H2 is parallel to the Cu(111) surface. Table shows the associated geometries and barrier heights, comparing to the previous values of the SRP PES, which gave dynamics results in agreement with experiments to chemical accuracy.[7] Analogous results for the MS-PBEl and MS-RPBEl functionals are given in Figures S1 and S2 and Tables S1 and S2 of the Supporting Information.

Figure 2

Elbow plots (i.e. V(Z,r)) resulting from the H2 + Cu(111) PES computed with the MS-B86bl functional and interpolated with the CRP method for four high symmetry configurations with the molecular axis parallel to the surface (θ = 90°) as depicted by the insets, for (a) top site and ϕ = 0°, (b) bridge site and ϕ = 90° (the bridge-to-hollow global minimum barrier geometry), (c) fcc site and ϕ = 0°, and (d) t2f site and ϕ = 120°. Barrier geometries are indicated with white circles, and the corresponding barrier heights are given in Table .

Table 3

H–H Distance at the Barrier (rb, in Å), the Molecule–Surface Distance at the Barrier (Zb, in Å), and the Barrier Height (Eb, in kJ/mol) as Computed for H2 + Cu(111) with the MS-B86bl Functionala

dissociation route	Eb	rb	Zb
bth	65.9 (60.6)	1.00 (1.03)	1.21 (1.16)
ttb	86.4 (86.0)	1.35 (1.40)	1.39 (1.39)
t2f (ϕ = 120°)	78.1 (74.4)	1.22 (1.27)	1.28 (1.27)
fcc (ϕ = 0°)	101.1 (97.7)	1.34 (1.59)	1.27 (1.27)

Values in parentheses are the SRP values from ref (7). The values are given for four different dissociation geometries (see Figure ). In all cases, H2 is parallel to the surface.

The barrier heights Eb computed with the MS-B86bl functional are in good agreement with the previous SRP results, overestimating the SRP barriers by 0.4–5.3 kJ/mol (Table ). This already suggests that the MS-B86bl functional should give a quite good description of dissociative chemisorption of H2 on Cu(111): the molecular beam sticking probabilities computed with the SRP functional slightly overestimated the experimental values, although agreement was achieved to within chemical accuracy. The barrier geometries obtained with the MS-B86bl functional are in quite good agreement with the SRP barrier geometries (Table ), except perhaps for the high barrier fcc geometry considered, for which the PES is quite flat in r around the barrier geometry (see Figure c). Note that the barriers tend to be a bit earlier (i.e., they occur at a smaller H–H distance) for the MS-B86bl PES than for the SRP PES. Finally, the MS-PBEl barrier heights are in even better agreement with the SRP results (see Table S1, comparing to Table ).

Finally, an important issue for the construction of semiempirical functionals is the tunability of the barrier height that can be achieved with them. With an SRP functional that is a weighted average of the GGA functionals PW91 (which is very similar to PBE[15]) and RPBE, the minimum barrier height for H2 + Cu(111) can be tuned between 46.8 and 78.9 kJ/mol (see Table ), that is, over quite a large range of 33 kJ/mol. However, with a trial SRP functional based on the MS-PBEl and MS-RPBEl functionals (which of the three MS functionals tested yield the lowest respectively the highest barriers, comparing Tables , S1, and S2) the minimum barrier height can only be tuned between 60.7 and 69.6 kJ/mol (Table ), a range of only 9 kJ/mol. The present results thus suggest that semiempirical functionals based on the made simple mGGA approach with eq can yield a much better description of the metal lattice, but at the cost of a reduced tunability of the reaction barrier height. Our finding for meta-GGA functions based on PBE-like and RPBE-like expressions is analogous to results of Garza et al. for the TPSS[30] (incorporating a PBE-like expression for the exchange enhancement factor) and their RTPSS functional (incorporating an RPBE-like expression).[78] To obtain a better agreement with molecular chemisorption energies (and greater tunability between TPSS and RTPSS), they relaxed the constraint that their RTPSS functional should reproduce the exact energy of the H-atom (i.e., that it should correct for self-interaction of this atom exactly).[78]

Table 4

Minimum Barrier Height (Eb, in kJ/mol) as Computed for H2 + Cu(111) with the Three Different MS Functionals, and with the PW91 and RPBE Functionals, for Bridge-To-Hollow Dissociationa

functional	Eb
PW91	46.8[7]
MS-PBEl	60.7
MS-B86bl	65.9
MS-RPBEl	69.6
RPBE	78.9[7]

In all cases, H2 is parallel to the surface.

Figure shows elbow plots of the MS-B86bl PES computed for four configurations in which H2 is parallel to the Ag(111) surface. Table shows the associated geometries and barrier heights, comparing to the previous values of the SRP48 PES,[46] which gave sticking probabilities that were shifted to higher incidence energies by 6.6–7.6 kJ/mol with respect to results of molecular beam experiments[45] (note that the SRP48 functional is an SRP functional for H2 + Cu(111)[12] but not for H2 + Ag(111)[46]). The results for the MS-PBEl functional are given in Figure S3 and Table S3.

Figure 3

Elbow plots (i.e. V(Z,r)) resulting from the H2 + Ag(111) PES computed with the MS-B86bl functional and interpolated with the CRP method for four high symmetry configurations with the molecular axis parallel to the surface (θ = 90°) as depicted by the insets, for (a) top site and ϕ = 0°, (b) bridge site and ϕ = 90° (the bridge-to-hollow global minimum barrier geometry), (c) fcc site and ϕ = 0°, and (d) t2f site and ϕ = 120°. Barrier geometries are indicated with white circles, and the corresponding barrier heights are given in Table .

Table 5

H–H Distance at the Barrier (rb, in Å), the Molecule–Surface Distance at the Barrier (Zb, in Å), and the Barrier Height (Eb, in kJ/mol) as Computed for H2 + Ag(111) with the MS-B86bl Functionala

dissociation route	Eb	rb	Zb
bth	129.5 (133.1)	1.22 (1.27)	1.12 (1.10)
ttb	152.9 (163.1)	1.51 (1.57)	1.50 (1.51)
t2f (ϕ = 120°)	145.9 (152.4)	1.40 (1.45)	1.33 (1.34)
fcc (ϕ = 0°)	159.4 (164.0)	1.57 (1.67)	1.32 (1.34)

Values in parentheses are the SRP48 values from ref (46). The values are given for four different dissociation geometries (see Figure ). In all cases, H2 is parallel to the surface.

The barrier heights Eb computed with the MS-B86bl functional are lower than the previous SRP48 results, underestimating the SRP48 barriers by 3.6–10.2 kJ/mol. This might be taken to suggest that the MS-B86bl functional should give a quite good description of dissociative chemisorption of D2 on Ag(111), as the SRP48 functional gave sticking probabilities that were shifted to higher incidence energies by 6.6–7.6 kJ/mol[46] with respect to experiment.[45] However, although the barrier geometries obtained with the MS-B86bl functional are in quite good agreement with the previous SRP48 barrier geometries (Table ), as for H2 + Cu(111) (see Table ), the barriers for H2 + Ag(111) tend to be a bit earlier (occur at smaller H–H distance) for the MS-B86bl PES. As discussed below, this is relevant also to the reaction dynamics, and dynamics calculations are required to see whether the MS-B86bl PES leads to higher sticking probabilities for D2 + Ag(111) than the SRP48 PES used earlier,[46] as would be needed for better agreement with the experiment.

Dynamics Results: Molecular Beam Sticking Probabilities

The S0 values computed with the three new MS functionals and with the SRP48 functional are compared to experimental values for H2 + Cu(111)[23,47] and D2 + Cu(111)[24] in Figure . As can be seen, the S0 values computed with the three new MS functionals are in excellent agreement with all experiments shown. The best agreement with the experiments of Auerbach and Rettner and co-workers[23,24] (Figure a–e) is obtained with the MS-B86bl functional, except perhaps for the lowest incidence energy. At the lowest incidence energies, the MS-PBEl functional would seem to give better results, but this may be an artifact of the use of the QCT method (reaction probabilities computed with the MS-B86bl PES being smaller than 0.002 in these cases), as the QCT method ignores tunneling contributions. The best agreement with the experiment of Rendulic and co-workers (Figure f) is obtained with the MS-PBEl functional, of which the overall performance is very similar to the performance of the SRP48 functional. However, these are details, and the main message is that the three MS functionals all yield excellent agreement with the molecular beam experiments, while also giving a very good description of the metal lattice.

Figure 4

S0 values computed with the three new MS functionals and with the SRP48 functional are compared to experimental values for H2 + Cu(111) measured in ref (23) (panels a–c), D2 + Cu(111) measured in ref (24) (panels d,e), and H2 + Cu(111) measured in ref (47) (panel f). Experimental results are presented in black, and computational results in red (MS-PBEl), green (B86bl), blue (MS-RPBEl), and purple (SRP48).

To obtain a measure of the quality of the functionals for H2 + Cu(111), one can compute the mean distance along the incidence energy axis from the computed S0 to the interpolated experimental values. In Figure , we show the MS-B86bl results comparing to those experiments for which enough data were available to perform cubic spline interpolation of the experimental data points. Computing the MAD (i.e., the mean of the calculated distances), we obtain MAD values of 0.3 kJ/mol for the experiment of Rettner et al.[23] using pure H2 beams (Figure a), 1.7 kJ/mol for the experiment using seeded D2 beams for a nozzle temperature of 2100 K of Michelsen et al.,[24] and 2.0 kJ/mol for the pure H2 beam experiment of Rendulic and co-workers.[47] This illustrates that the MS-B86b yields a chemically accurate description (MAD-values less than 4.2 kJ/mol) of the sticking of H2 and D2 on Cu(111). The same is true for the MS-PBEl functional (MAD values of 3.2, 4.1, and 0.3 kJ/mol, see Figure S4) and for the MS-RPBEl functional (MAD values of 1.2, 3.2, and 3.0 kJ/mol, see Figure S5), although the MS-B86bl functional gives a slightly better overall performance. Therefore, all three non-empirical MS functionals used individually yield agreement with these H2 + Cu(111) experiments to within chemical accuracy. These results also imply that it is possible to construct an SRP-DF on the basis of these functionals, for instance, an SRP-DF that is a weighted average of the MS-PBEl (lowest barriers) and the MS-RPBEl functional (highest barriers).

Figure 5

S0 values computed with the MS-B86bl functional (blue symbols) are compared to experimental values (red symbols) for H2 + Cu(111) measured in ref (23) (panel a), D2 + Cu(111) measured in ref (24) (panel b), and for H2 + Cu(111) measured in ref (47) (panel c). Horizontal lines indicate the distances along the energy axis between the computed S0 and the cubic-spline interpolated experimental curve, and the MAD is the mean value of these distances.

To put the performance of the three MS meta-GGA functionals in perspective, in Figure we compare the seeded beam D2 + Cu(111) experimental results of Michelsen et al.[24] to the S0 computed with two standard GGA functionals (PBE[15] and RPBE[14]) and three standard meta-GGA functionals (TPSS,[30] revTPSS,[31] and SCAN[32]). The PBE and RPBE results straddle the experimental S0 over a rather large energy interval, in agreement with the large tunability of SRP functionals taken as weighted averages of these functionals. As previously observed for H2 + Ru(0001),[60] the revTPSS functional improves over the performance of PBE, which is consistent with the design purpose of the former being to function well for both condensed matter physics (metals) and quantum chemistry (molecules). We also see an improvement going from PBE to TPSS, but not as much as with revTPSS. The SCAN functional gives the worst performance of all. This functional obeys a maximum number of known exact constraints and performs better than PBE for thermochemical data and reaction barriers of molecules and for lattice constants of solids.[32] Our finding that it performs more poorly than PBE for dissociative chemisorption may seem surprising but it is consistent with studies that find that SCAN overbinds more than PBE for chemisorption on metals.[78,79] The reasons for this are presently not fully understood. Garza et al. have speculated that the result that the most constrained non-empirical meta-GGA (i.e., SCAN) performs poorly at describing chemisorption on metals is due to inherent limitations of the form of semilocal functionals.[78] They seem to have based their suggestion partly on their finding that a meta-GGA functional (RTPSS) that performs quite well at describing molecular chemisorption on metals can be obtained by relaxing a constraint, that is, the constraint that the H-atom should be described exactly by correcting for the self-interaction.[78] Patra et al. have speculated that the SCAN functional should overbind CO on transition metal surfaces primarily due to density driven errors in the self-consistent SCAN energy.[79]

Figure 6

S0 values measured in seeded beam experiments on D2 + Cu(111)[24] (black symbols and lines) are shown as a function of Ei, comparing to the S0 values computed for these experiments with the PBE[15] (red), RPBE[14] (green), TPSS[30] (purple), rev-TPSS[31] (orange), and SCAN[32] (blue symbols and lines) functionals.

The revTPSS functional shows a MAD value with the D2 + Cu(111) experiment of Figure of 6.8 kJ/mol (see Figure S6). While this presents good agreement for a standard semilocal functional, it does not yet correspond to chemical accuracy. We also note that it is not possible to construct an SRP density functional on the basis of the three standard meta-GGA functionals tested: none of these functionals underestimates the reaction probability. However, it is probably possible to construct an SRP functional on the basis of one of these three meta-GGA functionals (with the best choice probably being the revTPSS functional) and the MS-RPBEl functional, which consistently underestimates the reactivity of H2 on Cu(111) (see Figure ), or the RTTPS meta-GGA functional, which has a performance on molecular chemisorption on metals that is comparable to that of the RPBE functional.[78]

The S0 values computed with the MS-PBEl and MS-B86bl functionals are compared to the experimental values for D2 + Ag(111)[45] in Figure . These two functionals perform slightly better than the SRP48 functional, which gave a MAD of 7 kJ/mol,[46] while MAD values of 4.5 and 5.5 kJ/mol were obtained with the MS-PBEl and MS-B86bl functionals, respectively. The improvement of the performance is not as large as one might assume based on barrier heights only (which for MS-B86bl were lower than the SRP48 values by 3.6–10.2 kJ/mol, see above). However, the sticking of D2 on Ag(111) is dominated by reaction of high vibrational states (the dominant contribution comes from v = 3).[46] The barriers are earlier on the MS-B86bl surface than on the SRP48 PES (see Table ), and as summarized by the Polanyi rules,[80,81] the later the barrier is, the more the reaction of vibrationally excited molecules is promoted. Thus, with the B86bl PES, the reaction of the high vibrational states is slightly less promoted, and this to some extent cancels the effect of the lower barriers. However, the main point is that the new MS functionals perform rather well for D2 + Ag(111) and in fact slightly better than the SRP48 functional while also yielding a very accurate description of the Ag(111) lattice.

Figure 7

S0 values computed with the MS-PBEl (blue symbols, upper panel) and MS-B86bl (purple symbols, lower panel) functionals are shown as a function of Ei, comparing to the values measured (red symbols) in molecular beam experiments,[45] for D2 + Ag(111). Also indicated are the distances along the energy axis of the computed points to the cubic spline interpolated experimental sticking probability curve. The MAD is the mean value of these distances.

Like the SRP48 functional, the MS-PBEl and MS-B86bl functionals, which may in principle be counted as SRP functionals for H2 + Cu(111) (i.e., with zero mixing coefficient of the other functional the SRP functional would be based on), are not yet transferable to H2 + Ag(111): they do not give a chemically accurate description of the existing molecular beam experiments for this system (such transferability was observed for the SRP density functional for methane + Ni(111),[10] to methane + Pt(111)[11]). As the MS-PBEl functional yields the highest reactivity of the three MS functionals tested here for H2 + Ag(111), our results suggest that it should not be possible to base an SRP-DF for this system on a combination of two of these three MS functionals. We attribute this result to the rather limited tunability of the MS functionals tested here. However, our present results for H2 + Ag(111) (Figure ) and for H2 + Cu(111) (Figure ) suggest that an SRP-DF for H2 + Ag(111) might be constructed on the basis of the MS-PBEl and one of the three meta-GGA functionals (SCAN, TPSS, or revTPSS). For further discussion on the agreement between theory and experiment for sticking of D2 on Ag(111), the reader is referred to ref (46).

Finally, initial-state selected reaction probabilities Pdeg(Ei,v,j) computed with the MS-PBEl and MS-B86bl functionals are compared with the values extracted from associative desorption experiments[82,83] and computed with the SRP48 functional[46] for H2 and D2 + Ag(111) in Figure . Interestingly, changing to the new MS functionals from the previously used SRP48 functional now leads to more distinct improvement than observed for the sticking (Figure and discussion above). The reason is that the associative desorption experiments were performed for low vibrational states of H2 and D2 (v = 0 and 1), so that the effect of the earlier barriers in the MS PESs is less pronounced than for the sticking, which is dominated by the reaction of v = 3 D2 as noted before.

Figure 8

Initial-state selected reaction probabilities Pdeg(Ei,v,j) computed with the MS-PBEl and MS-B86bl functionals are shown as a function of Ei, comparing with values extracted from associative desorption experiments[82,83] and computed with the SRP48 functional[46] for H2 and D2 + Ag(111). Results are shown for (v = 0, j = 2) D2 (a), (v = 1, j = 2) D2 (b), and (v = 0, j = 3) H2 (c).

Conclusions

The main goal of this study was to determine whether, with a meta-GGA functional constructed within the “made simple” approach, it would be possible to get a chemically accurate description of the dissociative chemisorption of H2 on Cu(111), while at the same time obtaining a better description of the Cu lattice than possible with previous SRP functionals based on the GGA. A second goal was to determine whether with the meta-GGA “made simple” functionals constructed here it should be possible to also get a more accurate description of the dissociative chemisorption of H2 on and its associative desorption from Ag(111) than was previously possible with the SRP48 GGA functional for H2 + Cu(111).

To determine the answer to these questions, we computed bulk lattice constants for Cu, Ag, Au, and Pt, interlayer lattice spacing relaxations for Cu(111) and Ag(111) and PESs for H2 + Cu(111) and Ag(111). We did this for three meta-GGA functionals based on the MS concept. In this approach, a function of the kinetic energy density is defined that effectively allows one to vary the exchange functional according to whether the binding in a certain region is metallic or covalent. The exchange parts of the three functionals (MS-PBEl, MS-B86bl, and MS-RPBEl) are based on the expressions for the PBE, B86b, and RPBE exchange functionals, respectively. The three new MS functionals yield metal bulk lattice constants with an accuracy intermediate between the all-purpose PBE GGA functional and the PBEsol functional, a GGA functional designed with the specific goal of accurately reproducing observables for the solid state. Likewise, the interlayer lattice spacing relaxations for the top two layers of Cu(111) and Ag(111) are in good agreement with experimental values.

The barrier heights and geometries obtained for H2 + Cu(111) were in good agreement with those obtained earlier with the original SRP functional for H2 + Cu(111). More importantly, the sticking probability curves computed with the three MS functionals and the QCT method agreed with experiments of Rettner and Auerbach and co-workers and of Rendulic and co-workers to within chemical accuracy. Similarly, the sticking probability curves computed with the MS-PBEl and MS-B86bl functionals for D2 + Ag(111) agree slightly better with the molecular beam experiments of Hodgson and co-workers than dynamics calculations based on the SRP48 GGA functional designed for H2 + Cu(111). The same is true for initial-state selected reaction probabilities computed for H2 and D2 + Ag(111) and the initial-state selected reaction probabilities extracted from associative desorption experiments of Hodgson and co-workers on these systems.

The main conclusions from our work are therefore that, considering the two systems investigated, (i) it is possible to construct non-empirical meta-GGA “made simple” functionals for these two H2–metal systems that describe the dissociative chemisorption reaction as accurately as previous semiempirical functionals based on GGA functionals, while simultaneously giving a more accurate description of the metal lattice, and (ii) on the basis of these MS functionals (in particular, MS-PBEl and MS-RPBEl), an SRP-DF can be constructed for H2 + Cu(111), but not for H2 + Ag(111). This limitation for H2 + Ag(111) is due to a potential disadvantage of the SRP approach based on the MS meta-GGA functionals tested: results for H2 + Cu(111) suggest that these candidate SRP functionals are less “tunable” than analogous semiempirical GGA functionals for barrier heights. Specifically, for H2 + Cu(111), the minimum barrier height varied by 9 kJ/mol going from MS-PBEl to MS-RPBEl, while it varied by 33 kJ/mol going from PBE to RPBE.

We also investigated the behavior of the standard meta-GGA functionals SCAN, TPSS, and rev-TPPS and noted that they all overestimated the sticking probability for D2 + Cu(111), so that they cannot be combined with each other to obtain an SRP functional for this system. However, the rev-TPSS functional provided a rather good description of the sticking of D2 on Cu(111), with a MAD of the computed and measured sticking probability curves of only 6.8 kJ/mol. Also, it might be possible to construct an SRP-DF for H2 + Ag(111) on the basis of the revTPSS and MS-PBEl functionals, which would probably also give a good description of the Ag lattice. More generally, a good strategy for constructing an SRP-DF based on meta-GGA functionals might be to start with a weighted average of MS-PBEl, and MS-RPBEl, if necessary replacing MS-PBEl with TPSS or revTPSS, or MS-RPBEl by RTPSS,[78] depending on whether a more attractive or more repulsive component functional is required. An alternative strategy for “casting the net wider” might be to relax slightly the condition that either the PBE-like or the RPBE-like MS functional for covalent bonding should exactly correct for the self-interaction correction of the H-atom, by allowing c to vary slightly in eq or 4b, as done in the RTPSS functional.[78]

The next step would be to apply the MS meta-GGA functionals to the other molecule–metal surface systems for which chemically accurate reaction barrier heights are now available [H2 + Pt(111) and methane + Ni(111), Pt(111), and Pt(211)]. We anticipate that for these systems it will be necessary to accurately model the van der Waals attractive interaction of the molecule with the metal surface. Preliminary results we obtained with the MS functionals investigated here for H2 + Pt(111) suggest that these functionals overestimate the reaction barrier height for this system. This is in line with previous findings that modeling the van der Waals interaction is necessary for an accurate description of weakly activated dissociative chemisorption of H2 on metal surfaces such as Ru(0001)[60] and of dissociative chemisorption of methane on Ni(111).[10,84]