Designing new SRP density functionals including non-local vdW-DF2 correlation for H₂ + Cu(111) and their transferability to H₂ + Ag(111), Au(111) and Pt(111)

Abstract

Specific reaction parameter density functionals (SRP-DFs) that can describe dissociative chemisorption molecular beam experiments of hydrogen (H₂) on cold transition metal surfaces with chemical accuracy have so far been shown to be only transferable among different facets of the same metal, but not among different metals. We design new SRP-DFs that include non-local vdW-DF2 correlation for the H₂ + Cu(111) system, and evaluate their transferability to the highly activated H₂ + Ag(111) and H₂ + Au(111) systems and the non-activated H₂ + Pt(111) system. We design our functionals for the H₂ + Cu(111) system since it is the best studied system both theoretically and experimentally. Here we demonstrate that a SRP-DF fitted to reproduce molecular beam sticking experiments for H₂ + Cu(111) with chemical accuracy can also describe such experiments for H₂ + Pt(111) with chemical accuracy, and vice versa. Chemically accurate functionals have been obtained that perform very well with respect to reported van der Waals well geometries, and which improve the description of the metal over current generalized gradient approximation (GGA) based SRP-DFs. From a systematic comparison of our new SRP-DFs that include non-local correlation to previously developed SRP-DFs, for both activated and non-activated systems, we identify non-local correlation as a key ingredient in the construction of transferable SRP-DFs for H₂ interacting with transition metals. Our results are in excellent agreement with experiment when accurately measured observables are available. It is however clear from our analysis that, except for the H₂ + Cu(111) system, there is a need for more, more varied, and more accurately described experiments in order to further improve the design of SRP-DFs. Additionally, we confirm that, when including non-local correlation, the sticking of H₂ on Cu(111) is still well described quasi-classically.

---

1 Introduction

Hydrogen (H₂) dissociation on various transition metal surfaces is an example of an intensely studied activated elementary surface reaction within surface science. Chemically accurate computation of rate-controlling states is essential in order to accurately describe the complex overall processes that take place during heterogeneous catalysis under real world conditions.^1–5 Industrially H₂ dissociation is an important step in the production of methanol from CO₂ over a Cu/ZnO/Al₂O₃ catalyst.^6–8 Additionally H₂ dissociation is an important process in the production of syngas and ammonia.⁴ Increasing the predictive power of theoretical models of heterogeneous catalysis potentially has an important financial impact on the chemical industry.⁹

An important step to increasing the predictive power of theoretical models is to create density functionals (DFs) that are chemically accurate for specific systems,^1,10–13i.e. DFs that can describe reaction barrier heights to within 1 kcal mol⁻¹.⁵ A next step is to investigate what ingredients of a DF that is chemically accurate for one system might make it transferable to another system without loss of accuracy. Presently, specific reaction parameter (SRP) density functional theory (DFT) is the only method that can describe the interaction of H₂ with metal surfaces with demonstrated chemical accuracy, while simultaneously being computationally cheap enough to make large comparative studies feasible. Therefore the design of accurate DFs is highly important to the field. The availability of transferable specific reaction parameter density functionals (SRP-DFs) has the potential to greatly speed up theoretical heterogeneous catalysis research. It does so by avoiding the need to design a new DF for each system of interest. The availability of transferable SRP-DFs would greatly improve the predictive power of theory for systems for which only sparse experimental results have been published.

The fitting of SRP-DFs is meticulous work and presently¹⁴ requires experimental data as reference data.^1,10–13,15 Transferability of a DF among systems in which one specific molecule interacts with surfaces of different metals has so far only been reported for the DF designed for CH₄ dissociation on Ni(111),¹³ which could also describe the dissociation of CH₄ on Pt(111) with chemical accuracy,¹⁶ where Ni and Pt belong to the same group. So far, for SRP-DFs fitted to reproduce molecular beam adsorption experiments for H₂ interacting with transition metals transferability was shown among systems in which H₂ interacts with different faces of the same metal,^17,18 but not among systems where the interaction is with surfaces of different metals.^12,15,19,20

The transferability of a SRP-DF that was fitted for the CH₄ + Ni(111) system to the CH₄ + Pt(111) system suggests that non-local correlation may be an important ingredient for a transferable DF, as this SRP-DF contains non-local correlation. For this reason here we investigate the design of new SRP-DFs that include non-local correlation for H₂ + transition metal systems. In this work we present two new SRP-DFs featuring GGA exchange but using non-local correlation. These DFs were fitted to experiments on the H₂ + Cu(111) system, since theoretically^{10,12,21–44} and experimentally^45–61 this is the best studied system. For this system we can have the most confidence that discrepancies between theory and experiment can be attributed to either shortcomings in DF design or the limitations of using the Born–Oppenheimer static surface (BOSS) model. It is well known that the BOSS model works well for activated H₂ dissociation on cold metals.^{23,27–29,62} Additionally we evaluate the performance of a SRP-DF that was fitted to the H₂ + Pt(111) system¹¹ for the H₂ + Cu(111) system.

In undertaking this study we have two aims. The first is to identify features of the newly constructed SRP-DFs that increase their transferability to other systems. Since we use the BOSS model, a direct assessment of the quality of a given DF is really only possible for molecular beam dissociative chemisorption experiments on reasonably cold metal surfaces. Apart from the H₂ + Cu(111) system we will consider such experiments for the H₂ + Pt(111) system^63,64 and the H₂ + Ag(111) system.⁵⁴ For the latter two systems there exists some uncertainty about the validity of the molecular beam parameters describing the experiments.^19,65 We however feel that, although there are some uncertainties in the parameters describing the experiments, nevertheless valuable insights on transferability can be derived by analysing the predicted reactivity of the SRP-DFs considered here.

Our second aim is to analyse the limits of our dynamical model to the extent that this is possible. We hope that a detailed analysis of the H₂ + Cu(111) system's large body of experimental work^45–61 will indicate how to proceed with improving the theoretical description of this system. To this end we will analyse both associative desorption experiments^47,51 and dissociative chemisorption experiments.^45,46 Naturally, our primary motive is to achieve chemical accuracy. We have also carried out a full quantum mechanical molecular beam simulation by carrying out a large number of fully initial-state resolved quantum dynamical (QD)^66,67 calculations for the H₂ + Cu(111) system. This is important because the inclusion of a van der Waals well in the PES might lead to discrepancies between quasi-classical trajectory (QCT)⁶⁸ and QD results compared to the good agreement that was obtained for these two methods for the H₂ + Cu(211) system³⁶ using the SRP48 DF,^10,62 which does not employ non-local correlation.

For the H₂ + Cu(111) system it is known that the effect of surface motion cannot readily be ignored for specific observables at high surface temperature⁶² (T_s). Analysing associative desorption and dissociative chemisorption experiments as linked through detailed balance,⁶⁹ might allow us to disentangle the effects of surface motion and the non-adiabatic contributions of electron–hole pair (ehp) excitations, a methodology that was suggested by Shuai et al.⁷⁰ for the H₂ + Au(111) system. If detailed balance is applicable then an analysis of both associative desorption and dissociative chemisorption experiments should yield the same result. A detailed analysis might therefore allow us to identify which dynamical effects not included in the BOSS model may have to be included in future work. Here we will make a direct comparison to experimental effective barrier (E₀(ν, J)) parameters.⁴⁷ Even though a complementary molecular beam dissociative chemisorption experiment is not available for the associative desorption experiment of Shuai et al.⁷⁰ on H₂ + Au(111), we extend our analysis also to this system. We 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,^71–75 and that this has also been done for H₂ and D₂ desorbing from Cu(111). However, the calculations published so far do have limitations. The early work^72,73 used a PES that is an approximate fit³⁹ to unconverged DFT calculations⁴⁰ using the PW91 DF,⁷⁶ and the statistical accuracy of the results of the later work⁷⁵ was limited by the number of ab initio molecular dynamics (AIMD) trajectories that could be run. However, an interesting aspect of the later work⁷⁵ is that the effects of surface atom motion and electron–hole pair excitation could also be investigated.

Furthermore we will treat vibrationally and rotationally inelastic scattering for the H₂ + Cu(111) system, since the opinion has been voiced that these properties might be extra sensitive to the van der Waals well, which is present in potential energy surfaces (PESs) computed here with the use of non-local correlation.²⁴

For the H₂ + Ag(111) system the only molecular beam dissociative chemisorption experiments we are aware of are those of Hodgson and coworkers.⁵⁴ We will also make a comparison to initial-state resolved associative desorption experiments for this system,^77,78 as was recently done by Jiang and Guo⁷⁹ using quantum dynamics calculations on a permutation invariant polynomial (PIP) neural network potential⁸⁰ and in work done in our group.^12,19 Earlier experimental work on the H₂ + Ag(111) system suggested that H₂ prefers to physisorb on silver surfaces,^81–83 and that the dissociation of H₂ on silver is endothermic,⁸⁴ exhibiting a relatively low barrier for associative desorption of H₂.^84–86 For this reason some earlier experimental studies have focused on scattering at low translational energies of H₂ from Ag(111)^87–90 and Ag(110).⁹¹ Recent theoretical studies that addressed the effects of electron–hole pair excitation have shown very interesting effects at low translational energies with respect to inelastic scattering and dynamical steering.^92–94 The non-adiabatic energy loss during the dynamics was however shown to be small.⁹⁴ This is in agreement with work on H₂ + Cu(111), which also showed little effect of electron–hole pair excitation on sticking.^22,95 Therefore we presume the BOSS model to be accurate enough for our first aim, which is to identify the features of SRP-DFs that contribute to their transferability.

Our DFT calculations using van der Waals correlation functionals^96,97 also yield results regarding the geometry and the depth of the van der Waals wells for the systems investigated, and we will compare these results to the experimental results for the systems investigated here. In many cases the experimental results come from an analysis of experiments on selective adsorption.^{59,61,89,90,98–104} In these experiments, an increase or a dip is observed in a peak for a diffractive (corrugation mediated selective adsorption, CMSA^105,106) or a rotational (rotationally mediated selective adsorption, RMSA⁹⁸) transition if the translational energy goes through a value that coincides with the energy between two parallel translational or hindered rotational metastable states, respectively. In the final state, the H₂ molecule is trapped in the van der Waals well close to the surface.^100,103 Information about the resonance energies that are present can be used to reconstruct the shape of the potential, and thus the van der Waals well geometries and well depths. The H₂-(111) metal–surface systems investigated here exhibit little corrugation. For this reason experiments using RMSA of HD, in which rotational excitation is used to probe the bound levels of the gas-surface potential, have been particularly important for gathering information concerning van der Waals interactions in these systems. The off-center position of the center of mass of HD results in very pronounced resonances when using the RMSA technique.^61,98–101 Experimentalists have also been able to carry out RMSA measurements using H₂ or D₂ instead of HD.^89,90,102 van der Waals well depths can also be obtained from the temperature dependence of the Debye–Waller attenuation of peaks for (rotationally) inelastic diffraction,^107,108 from potential inversion using calculations on (rotationally inelastic) diffractive scattering using the eikonal approximation,¹⁰⁷ and from potential inversion using measurements on phonon-assisted RMSA (also called rotationally mediated focused inelastic resonances, RMFIR¹⁰⁹). Concerning the systems investigated here, studies using experiments to analyze the van der Waals interaction have been performed for H₂ + Cu(111),^59,61 H₂ + Ag(111),^89,90,101 H₂ + Au(111),⁶¹ and H₂ + Pt(111).^98–100,110

This paper is organized as follows. Section 2 covers the computational methods used, with Section 2.1 discussing the coordinate system we use. Section 2.2 details the construction of SRP-DFs, and Section 2.3 the interpolation of the PESs based on SRP DFT data. The QCT and QD methods are described in Sections 2.4 and 2.5, respectively. The way in which we calculate observables is discussed in Section 2.6, and the computational details are summarized in Section 2.7.

Section 3 is the results section. Section 3.1 concerns the electronic structure of the investigated systems and the fitted PESs for these systems. Section 3.2 concerns simulations of molecular beam sticking experiments for three H₂ (D₂) + metal systems. Initial-state resolved reaction probabilities are presented in Section 3.3, and Section 3.4 presents parameters, E_1/2(ν, J), that can be used to compare with associative desorption experiments. Rotational quadrupole alignment parameters, which characterize the dependence of reaction on the alignment of H₂ relative to the surface (parallel or end-on) are presented in Section 3.5 for the D₂ + Cu(111) system. Rotationally and vibrationally inelastic scattering of H₂ from Cu(111) is considered in Section 3.6. State distributions of molecules desorbing from Au(111) are presented in Section 3.7.

The discussion is presented in Section 4. Calculated metal properties are discussed in Section 4.1, and static properties of the calculated PESs in Section 4.2. Molecular beam sticking results are discussed in Section 4.3. In Section 4.4 we analyze initial-state resolved reaction probabilities, and compare extracted E_1/2(ν, J) parameters with values extracted from associative desorption experiments on H₂ + Cu(111) and Au(111). We also compare initial-state resolved reaction probabilities computed with theory with experimental data for H₂ + Ag(111), and the state distributions of molecules desorbing from Au(111). In Section 4.4 we also discuss theoretical and experimental values of rotational quadrupole alignment parameters for D₂ + Cu(111). The inelastic scattering of H₂ from Cu(111) is discussed in Section 4.5 The quality of the QCT results is evaluated by comparison with QD results for H₂ + Cu(111) in Section 4.6. Section 4.7 contains an overview of the comparison theory-experiment for all the H₂–metal surface systems we treat in this work. The transferability of the created SRP-DFs is discussed in Section 4.8. In Section 4.9 we discuss how and when qualitative differences in the adiabatic descriptions of molecular beam sticking and associative desorption can be interpreted as a possible signature of the role of ehp excitation. Conclusions are presented in Section 5.

2 Methodology

2.1 Coordinate system

All calculations in this work are carried out using the BOSS model.¹⁰ This means that the atoms of the metal slab are fixed to their ideal lateral positions, and that we only take into account the six molecular degrees of freedom (DOF) of the impinging H₂ (D₂) molecule (see Fig. 1a). Three of the DOFs taken into account are the centre of mass (COM) coordinates X, Y and Z, where (X,Y) describes the lateral position of the molecule and Z describes the molecule–surface distance. The other DOFs are the H₂ bond length r, the polar orientation angle with respect to the surface normal θ and the azimuthal angle ϕ. The geometry of the (111) face of an fcc metal together with its high symmetry sites is shown in relation to the coordinate system used in Fig. 1b.

Fig. 1 The COM coordinate system used for the description of the H₂ (D₂) molecule (a). The unit cell of a (111) face of a fcc metal together with the high symmetry sites as well as the relationship with the coordinate system chosen for H₂ (D₂) relative to the (111) surface (b). The origin of the COM coordinate system (X,Y,Z) = (0,0,0) is at an atom in the top surface layer (a top site). We define the polar angle and azimuth such that (θ = 90°,ϕ = 0°) corresponds to molecules parallel to the surface pointing along the X (or equivalent U) direction. The hcp and fcc hollow sites correspond to metal atoms in the second and third layer, respectively.

2.2 SRP DFT

We use periodic DFT calculations to construct PESs, testing DFs at the GGA, the meta-GGA (mGGA), and GGA + non-local correlation level, where non-local correlation refers to either vdW-DF1 non-local correlation⁹⁶ or vdW-DF2 non-local correlation⁹⁷ (see the types of DF defined in Table 1). In many cases we test a DF that is a SRP-DF for at least one of the four systems considered. In the present context a SRP-DF is constructed by taking a weighted average of a DF that overestimates the sticking probability, and one DF that underestimates the sticking probability for the system of interest.¹⁰ More specifically, the mixing often occurs for the exchange part of the DF, as was done for many of the previously developed SRP-DFs.^11,13,16,26

Table 1 The exchange–correlation DFs used in this work are presented in alphabetical order. Also shown is the type of each functional as well as the constituent exchange and correlation parts

Name	Type	Exchange	Correlation
B86SRP68-DF2	vdW-DF	0.68 B86r^111 + 0.32 RPBE^112	vdW-DF2^97
MS-B86bl^12	mGGA	MS-B86bl^12	revTPSS^113
MS-PBEl^12	mGGA	MS-PBEl^12	revTPSS^113
optPBE-DF1^114	vdW-DF	optPBE^114	vdW-DF1^96
PBE^115	GGA	PBE^115	PBE^115
PBEα57-DF2^11	vdW-DF	PBEα = 0.57^116	vdW-DF2^97
PBEsol^117	GGA	PBEsol^117	PBE^115
RPBE^112	GGA	RPBE^112	PBE^115
SRP48^62	GGA	0.52 PBE^115 + 0.48 RPBE^112	PBE^117
SRPsol63-DF2	vdW-DF	0.63 PBEsol^117 + 0.37 RPBE^112	vdW-DF2^97
vdW-DF1^96	vdW-DF	revPBE^118	vdW-DF1^96
vdW-DF2^97	vdW-DF	rPW86^119	vdW-DF2^97

---

The SRP-DFs developed in this work are all DFs in which the exchange part is taken at the GGA level, and the exchange correlation functional takes the following form:

Here α is the SRP mixing parameter, ρ is the three dimensional electron density and ∇ρ is the gradient of the electron density. E¹_X(ρ,∇ρ) and E²_X(ρ,∇ρ) are the two DFs that are to be mixed into the exchange part of the SRP-DF. The non-local correlation part here can correspond to the non-local correlation used in the vdW-DF1 or vdW-DF2 DFs.^96,97 The DFs used in this work (i.e., the B86SRP68-DF2 and SRPsol63-DF2 DFs), and the other DFs considered in this work are shown in Table 1. Table 1 also shows the type of each DF and the exchange and correlation components contained in each function (see eqn (1)).

2.3 Construction of the PESs

A continuous representation of the PESs is obtained by the interpolation of DFT results calculated on a grid using the corrugation reducing procedure (CRP).¹²⁰ The method we use is analogous to that used by Wijzenbroek et al.,⁴¹ but we used denser grids to represent the full six dimensional molecule–surface interaction potential and the three dimensional atom–surface interaction potential to further increase the accuracy of the resulting CRP¹²⁰ PESs with respect to the underlying DFT calculations. Details are presented in the ESI.†

2.4 Quasi-classical dynamics

The QCT calculations presented in this work are carried out on six dimensional PESs and assume quasi-classical initial conditions.⁶⁸ This means that we take into account the quantum mechanical energies of the impinging H₂ and D₂ molecules in their initial rovibrational states. The method used is described more fully in ref. 121. The equations of motion are integrated using the method of Stoer and Bulirsch.¹²²

When simulating a molecular beam experiment 200 000 trajectories are propagated per energy point, and when calculating initial-state resolved reaction probabilities 50 000 trajectories are propagated per energy point. All trajectories start in the gas phase, at Z_gas = 8 Å. For all QCT calculations we use a time step of dt = 0.001 fs. Trajectories are assumed to result in reaction if r becomes bigger then some critical value r_c (2.2 Å) and in scattering if Z becomes bigger then Z_gas which is also the starting point of all trajectories.

2.5 Quantum dynamics

Six dimensional quantum dynamics (QD) calculations are performed by solving the time-dependent Schrödinger equation,using the time-dependent wave packet (TDWP) method^123,124 with our in-house computer package.^66,67 Here, Ψ(Q⃑,t) denotes the nuclear wave function of H₂ at time t with Q⃑ being the position vector. Furthermore we employ the following Hamiltonian in order to take into account the six degrees of freedom of H₂:Here, M and μ are the mass and reduced mass of H₂, Ĵ²(θ,ϕ) is the angular momentum operator and V(Q⃑) is the six dimensional PES. The scattered wave packet is analysed using the scattering matrix formalism,¹²⁵ yielding fully initial-state resolved S-matrix elements for vibrationally, rotationally, and diffractionally inelastic scattering. From the S-matrix elements the corresponding state-to-state probabilities for scattering at the incident energy E can be obtained.^66,67 Subsequently the sticking probability can be computed^66,67 by subtracting the sum of the scattering probabilities from one.

Further information on how we construct the initial wave packet can be found in Section S2 of the ESI.† Table S1 (ESI†) presents parameters describing all the initial wave packets used, and Table S2 (ESI†) shows the rovibrational states taken into account in the molecular beam simulations of sticking with the QCT and QD methods.

2.6 Computation of observables

2.6.1 Simulating molecular beam sticking. In order to compute molecular beam sticking probabilities the translational energy and rovibrational state distributions need to be taken into account according to the nozzle temperature, T_n. The sticking probability S₀ is computed usingHere, 〈E〉 is the average translational energy, and the probability for a molecule present in the beam to be in a rovibrational state described by the vibrational quantum number ν and the angular momentum quantum number J and to have a velocity between v and v + dv is denoted by:

P(v,ν, J,T_n)dv = P_flux(v;T_n)dv × P_int(ν, J,T_n). (5)

The flux-weighted velocity distribution P_flux is a function of T_n and is determined by the width parameter α and the stream velocity v₀ according to⁴⁴

P_flux(v;T_n)dv = Cv³e^{−(v−v₀)2/α2}dv (6)

with C being a normalization constant. Through eqn (4) and (5) the reactivity of each rovibrational state is weighted according to its Boltzmann weight as follows:with

f(ν, J,T_n) = (2J + 1) × e^{(−(E_ν,0−E_0,0)/k_BT_vib)} × e^{(−(E_ν,J−E_ν,0)/k_BT_rot)}. (8)

Here, k_B is the Boltzmann constant and E_ν,J is the energy of the quantum state characterized by ν and J. The first and second Boltzmann factor describe vibrational and rotational state populations, respectively. Here we take into account the effect of rotational cooling during the supersonic expansion by taking the rotational temperature to be T_rot = 0.8 × T_n,⁵¹ while the vibrational temperature, T_vib, is taken to be equal to T_n. The factor g_N in eqn (7) reflects the ortho/para ratio of hydrogen in the beam. For H₂, g_N is 1/4 (3/4) for even (odd) values of J, and for D₂, g_N = 2/3 (1/3) for even (odd) values of J.

In the QCT calculations presented in this work the probability distribution P(v, ν, J, T_n) is randomly sampled as described in ref. 121. All parameters describing molecular beam experiments used for the calculations presented here can be found in Table S3 (ESI†). The reaction probability is then computed by dividing the number of adsorbed trajectories, N_ads, by the total number of calculated trajectories, N, i.e. P_r = N_ads/N.

We compute initial-state selected but degeneracy averaged reaction probabilities, P_deg(E, ν, J), as:

where P_r(E, ν, J, m_J) is the fully initial-state resolved reaction probability, m_J is the magnetic rotational quantum number, and E is the translational energy. For the QD calculations it is not possible to directly sample P(v, ν, J, T_n). Molecular beam reaction probabilities for the QD method are instead calculated from initial-state resolved reaction probabilities in the same manner as discussed in our previous work.³⁶

2.6.2 Comparing to experimental E₀(ν, J) parameters. Experimentally, for H₂–metal surface reactions the initial state-selected reaction probabilities are usually obtained^45,47,51,70 from associative desorption measurements using the principle of detailed balance.⁶⁹ Experiments on H₂ associatively desorbing from metals typically measure the (unnormalized) state-resolved translational energy distributions of molecules from the surface using resonance-enhanced multi-photon ionization (REMPI).^18,45,51 These distributions, P_des(E, ν, J), may be related to the degeneracy averaged initial-state resolved reaction probability, using:The extracted reaction probabilities are usually fitted to a sigmoid function, e.g. the function involving the error function:Here, the A_ν,J values are the saturation values of the extracted degeneracy averaged reaction probabilities, and the effective barrier height (E₀(ν, J)) is the incidence energy at which P_deg(E, ν, J) first becomes equal to . Using eqn (11), E₀(ν, J) also is the inflexion point of the reaction probability curve if the saturation value A_ν,J corresponds to the absolute saturation value. W_ν,J represents the width of the reaction probability curve.

If the proportionality factor implicit in eqn (10) would also be measured in the experiment (for instance, because the exact state-selective flux would have been measured), it should be possible to directly extract absolute values of the initial-state selected reaction probabilities from the experiment by assuming detailed balance. In this case, the directly extracted A_ν,J value would be the true saturation value of the reaction probability. These values could then be directly compared to computed degeneracy averaged initial-state resolved reaction probabilities (eqn (9)).

In practice, even if the prefactor in eqn (10) is not measured in associative desorption, it is still possible to extract normalized values from a wholly experimental procedure, if measured sticking probabilities are also available. This has been done for H₂ + Cu(111)^47,51 and D₂ + Cu(111).^45,47,62 In this method, which we call Method A1, essentially the sticking probability is described in terms of the initial-state selected probabilities extracted from the associative desorption experiments, thereby obtaining the saturation values describing the latter. Method A1 is described more fully in the ESI.†

If no sticking experiments are available, as for the H₂ + Au(111) system also studied here, the experimentalists may chose not to normalize the extracted reaction probabilities in an absolute sense. However, the extracted reaction probabilities may still be normalized relative to one another. This was done in recent experiments on H₂ and D₂ + Au(111), in which A_ν,J was set to one for (ν = 0,J = 6) H₂ and for (ν = 0,J = 0) D₂, and the A_ν,J values for different (ν, J) states of H₂ and D₂, respectively, reflected the values of the reaction probabilities relative to these reference states.⁷⁰ We will call this Method A2, where the A in A2 emphasizes that this method is also wholly experimental.

If no measured sticking probabilities are available for the system of interest, one may still choose to normalize reaction probabilities extracted from associative desorption experiments, but now with reference to theory.^47,70 We label such methods with “B” to emphasize that the normalization is done with reference to theory. In the methods we are aware of, the experimentalists define a translational energy E_max(ν, J), which is the maximum translational energy for which the not yet normalized value can still be accurately extracted using eqn (10).^47,70 At higher E this becomes difficult because the desorption flux becomes small due to the exponential factor in eqn (10), leading to too much noise in the determined P_deg(E,ν, J). Parameters can then be described in two ways (Methods B1 and B2). Briefly, in Method B1 the saturation parameters are determined by setting them equal to the theoretical reaction probability computed for E = E_max(ν, J). The E_1/2(ν, J) parameters is then determined as the energy at which the computed reaction probability equals half the saturation value. Method B2 aims to improve upon this. Method B1 was previously followed in experimental papers to enable comparison with theory for H₂, D₂ and HD desorbing from Cu surfaces,⁴⁷ and for H₂ and D₂ desorbing from Au(111).⁷⁰ Further details of Methods B1 and B2 are presented in the ESI.†

2.6.3 Rotational quadrupole alignment parameters. The rotational quadrupole alignment parameter, A⁽²⁾₀(E,ν, J), is computed from initial-state resolved reaction probabilities as follows:¹²⁶The rotational quadrupole alignment parameter is a measure of the dependence of the reaction on the alignment of H₂ relative to the surface.

2.6.4 Rovibrational state populations of H₂ and D₂ desorbing from Au(111). State distributions of desorbing molecules are calculated in the following manner:⁷⁰Here T_S is the surface temperature, and E_max(ν, J) is the maximum kinetic energy sensitivity of the experiment,⁷⁰ which is plotted as a function of J in Fig. S2 (ESI†). To make a comparison between theory and experiment possible, the experimental P_deg(E,ν, J) are replaced by the error function expressions of ref. 70. In order to make this comparison valid, we only integrate eqn (13) up to E_max(ν,J). The error function fits derived in ref. 70 are only valid below E_max(ν,J) and can yield sticking probabilities substantially bigger than one for higher energies. Integration of eqn (13) is done by taking a right Riemann sum with a dE of 0.2 meV. The N(ν, J) populations are normalized to the total ν = 0 population according to:The ratio ν = 1 : ν = 0 can then calculated as:The upper limits to J used in eqn (14) and (15) are discussed below.

2.7 Computational details

All the electronic structure calculations were carried out by performing plane wave periodic DFT calculations using a user modified version of the Vienna Ab Initio Simulation Package^127–130 (VASP). The modification of the computer package concerns an interface to VASP with the LibXC density functional library.¹³¹ The standard VASP projector augmented wave (PAW) potentials¹³² and vdW-DF correlation^96,97 as implemented in VASP^133,134 were used for all calculations at the GGA level except for the SRP48 calculations on Pt(111), for which the standard VASP ultrasoft pseudopotentials¹³⁵ and PBE correlation¹¹⁵ were used. Calculations done using a mGGA use mGGA correlation (see Table 1).

All calculations at the GGA level presented in this work have been carried out using a plane wave cutoff energy of 450 eV together with smearing of 0.2 eV using the Methfessel–Paxton method of order 1. The input parameters for calculations with a mGGA DF can be found in ref. 12. Lattice constants have been calculated using a four atom bulk unit cell and a 28 × 28 × 28 Monckhorst–Pack k-point grid. All metal slabs consist of six layers of which the bottom two layers were fixed at the ideal bulk interlayer distance. Slab relaxation has been carried out using a 1 × 1 supercell, a 32 × 32 × 1 Γ-centered k-point grid and a vacuum distance of 16 Å. PES calculations have been carried out using a 3 × 3 supercell, a 11 × 11 × 1 Γ-centered k-point grid and a vacuum distance of 16 Å.

3 Results

3.1 Electronic structure

3.1.1 Description of the metal. Table S4 (ESI†) shows calculated lattice constants for different DFs, comparing with zero-point energy corrected experimental results.¹³⁶Table 2 shows the percentage change of the distance between the top two layers in the metal slab relative to the calculated bulk interlayer distance, also comparing to experimental results.^{56–58,137–141}

Table 2 Relaxations of the interlayer distance between the top two layers relative to the bulk interlayer distance in %

	Cu	Ag	Au	Ni	Pd	Pt
Exp.	−1.0%,^56,58 −0.7%^57	−2.5%,^137 −0.5%^138	1.5%^139	−0.07%^56	1.3%^140	1.1%^141
SRP48^10	−1.3%^10	−0.6%^19	0.6%^142	−1.0%^20	0.0%	1.0%
VdW-DF1^96	−0.2%	0.1%	1.6%	−1.1%	0.7%	1.3%
vdW-DF2^97	0.0%	0.5%	2.1%	−1.1%	0.8%	1.5%
B86SRP68-DF2	−0.4%	−0.1%	1.3%	−1.1%^20	0.7%	1.2%
SRPsol63-DF2	−0.4%	−0.2%	1.4%	−1.1%	0.6%	1.0%
PBEα57-DF2^11	−0.4%	0.0%	1.5%	−0.8%^20	0.6%	0.8%
optPBE-DF1^114	−0.8%^41	−0.2%	1.0%		0.2%	0.8%
MS-B86bl^12	−1.0%	−0.5%	1.0%		0.3%	1.0%
PBE^115	−0.3%^143	−0.2%^143	1.0%^143		−0.5%^143	0.9%^143
PBEsol^117	−0.4%^143	−0.1%^143	0.8%^143		0.6%^143	0.8%^143

---

3.1.2 H₂ + metal surface PESs. Barrier heights and geometries for H₂ + Cu(111) for high symmetry geometries are shown in Table 3. The energetic corrugation ξ, which is the highest minus the lowest barrier height, is shown as well. Elbow plots for four geometries are shown in Fig. 2 for the B86SRP68-DF2 SRP-DF.

Table 3 Barrier heights for H₂ reacting on Cu(111) for high symmetry geometries. For the bridge site ϕ = 90° and θ = 90°, for the t2b and hcp geometries ϕ = 0° and θ = 90°. The high symmetry locations are shown in Fig. 1b, the t2b geometry refers to the COM of the molecule being placed on a top site, with the H-atoms pointing to the nearest bridge sites. Barrier heights are in eV, and the barrier positions in Å. Additionally the energetic corrugation, ξ, is also shown in eV

	Bridge			t2b			hcp			ξ
	E b	r b	Z b	E b	r b	Z b	E b	r b	Z b
SRP48^10,41	0.636	1.030	1.172	0.887	1.396	1.394	1.047	1.539	1.269	0.411
B86SRP68-DF2	0.725	1.045	1.169	0.936	1.380	1.392	1.056	1.379	1.273	0.331
SRPsol63-DF2	0.712	1.045	1.167	0.925	1.379	1.392	1.038	1.373	1.270	0.326
PBEα57-DF2	0.720	1.063	1.142	0.934	1.385	1.390	1.035	1.373	1.268	0.315
optPBE-DF1^41	0.712	1.053	1.165	0.915	1.382	1.396	1.070	1.427	1.271	0.358
MS-B86bl^12	0.683	0.997	1.205	0.895	1.351	1.391	1.079	1.369	1.266	0.396

---

Fig. 2 Elbow plots, i.e. V(Z,r) resulting from the H₂ + Cu(111) PES computed using the B86SRP68-DF2 DF and interpolated using the CRP method for four high symmetry geometries in which the molecular axis is parallel to the surface (θ = 90°) as depicted by the insets for (a) the top site and ϕ = 0°, (b) the bridge site and ϕ = 90°, (c) the fcc site and ϕ = 0°, and (d) the t2f site and ϕ = 120°. Barrier positions are indicated with white circles.

Barrier heights and positions for H₂ + Ag(111), Au(111) and Pt(111) are shown in Tables S5, S6 and S7 (ESI†) respectively. With respect to the H₂ + Pt(111) system the most striking result is that only the PBEα57-DF2¹¹ and MS-PBEl¹² DFs exhibit a double barrier structure for the top-to-bridge (t2b) geometrie whereas the other DFs tested do not. The PBEα57-DF2¹¹ SRP-DF is the only DF that predicts a negative early barrier to reaction for this reaction.

We have checked the fit accuracy of our CRP¹²⁰ PES for the B86SRP68-DF2 DF for H₂ + Cu(111) using ∼4900 randomly sampled geometries. Based on all the randomly sampled points taken together our CRP¹²⁰ fit has a root mean square (rms) error of 31 meV. When only looking at the 3538 geometries that have an interaction energy of H₂ with the surface lower then 4 eV the rms error reduces to 8 meV (∼0.2 kcal mol⁻¹). Our CRP¹²⁰ PES is thus highly accurate with respect to the underlying electronic structure calculations. Since the other PESs calculated for this paper have been constructed in the same manner, we presume their accuracy with respect to the electronic structure calculations to be similar, and high enough for our purposes.

3.1.3 van der Waals wells. Fig. 3a and b show van der Waals wells for H₂ in a parallel orientation (ϕ = 0°,θ = 90°) above a top atom of Cu(111) for different DFs, comparing with experimental results.⁵⁹ Panel a shows calculated van der Waals wells for the non-standard DFs with non-local van der Waals correlation investigated in this work. Panel b shows van der Waals well depths for SRP-DFs developed previously in our group,^10,12 as well as for the two standard vdW-DF1⁹⁶ and vdW-DF2⁹⁷ DFs. Agreement with experiment is best for the B86SRP68-DF2 DF.

Fig. 3 The computed interaction energy of H₂ parallel to the Cu(111) surface above a top site (ϕ = 0°,θ = 90°) is compared with experimental results⁵⁹ (black). Panel (a) shows the calculated van der Waals well for the B86SRP68-DF2 (red), PBEα57-DF2, SRPsol63-DF2 and optPBE-DF1⁴¹ DFs, and panel (b) for the SRP48,¹⁰ MS-B86bl,¹² vdW-DF1⁹⁶ and vdW-DF2⁹⁷ DFs.

All van der Waals well depths and geometries for the systems and DFs investigated in this work are tabulated in Table S8 (ESI†). With respect to the van der Waals well depths for H₂ + Ag(111), H₂ + Au(111), and H₂ + Pt(111) we find depths that are in good agreement with experimental work^61,90,99,110 for the B86SRP68-DF2 DF.

3.2 Molecular beam sticking probabilities

Molecular beam sticking probabilities computed with five DFs for H₂ and D₂ reacting on Cu(111) are shown in Fig. 4, comparing to experimental results of Auerbach and coworkers^45,51 and Rendulic and coworkers.⁴⁶ Fig. S5 (ESI†) shows the comparison with an additional experiment for the five DFs discussed here, namely for the pure D₂ molecular beams of Auerbach and coworkers.⁴⁵ The difference between theory and experiment is assessed by determining how far the theoretical result needs to be shifted along the incidence energy axis to be superimposed on a spline interpolated curve going through the experimental results. Values of the mean absolute deviation (MAD) are calculated as the mean of the absolute number of these shifts for a particular set of molecular beam experiments. From the MAD values it can be seen that all five DFs considered describe the experiments on H₂ + Cu(111) shown in Fig. 4 with chemical accuracy. Fig. 5 shows comparisons to two additional sets of molecular beam experiments of Auerbach and coworkers^45,51 for H₂ reacting on Cu(111) for a more limited set of DFs.

Fig. 4 Molecular beam sticking probabilities for H₂ and D₂ reacting on Cu(111) for three sets of molecular beam experiments, as computed with five SRP-DFs. Experimental results are shown in red,^45,46,51 QCT results are shown in blue, and QD results are shown in green in two panels (panels e and f). 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⁻¹. The MAD values are the mean value of the energy differences for each experiment.

Fig. 5 Molecular beam sticking probabilities for H₂ reacting on Cu(111). Experimental values are shown in black,⁵¹ computed reaction probabilities are shown for the B86SRP68-DF2 (red), PBEα57-DF2 (blue), and the SRP48 (green) DFs. The solid lines correspond to QCT calculations, the dashed lines to QD calculations.

For the B86SRP68-DF2 DF QD results are also shown for the experiments concerning H₂ in Fig. 4e, f and in Fig. 5. Note that in these figures the QCT results are based on more rovibrational states than those included in the QD calculations (see Table S2, ESI†). In Fig. 6 we also explicitly compare QCT and QD results obtained while averaging over the same rovibrational states, and compare those to the QCT results shown in Fig. 4e, f and in Fig. 5. Fig. S6 (ESI†) shows the same data as shown in Fig. 6, but with the reaction probability on a log-scale. Table S9 (ESI†) displays the computed MADs for all experiments considered.

Fig. 6 Molecular beam sticking probabilities for H₂ reacting on Cu(111). Experimental values are shown in black,⁵¹ B86SRP68-DF2 QCT results are shown in green, B86SRP68-DF2 QCT results based on the same rovibrational states as taken into account in the QD calculations in blue, and B86SRP68-DF2 QD results in red.

Fig. 7 shows molecular beam sticking experiments for D₂ reacting on Ag(111). Experimental results of Hodgson and coworkers¹⁴⁴ are also shown. The calculated results are obtained using the pure D₂ molecular beam parameters of Auerbach and coworkers⁴⁵ obtained from experiments on Cu(111). The DFs treated in this work, as well as the MS-PBEl mGGA,¹² reproduce the experiment to almost within chemical accuracy. The SRP48 DF^19,62 yields the worst and the MS-PBEl¹² the best performance.

Fig. 7 Molecular beam sticking probability as a function of the average incidence energy for D₂ reacting on Ag(111). Experiment is shown in black.¹⁴⁴ QCT results are shown for the following DFs: SRP48 (purple),¹⁹ MS-PBEl (red),¹² B86SRP68-DF2 (blue), and PBEα57-DF2 (green). The dotted line represents the experimental curve shifted by 4.2 kJ mol⁻¹, denoting the limit to chemical accuracy.

In Fig. 8 molecular beam sticking probabilities for D₂ reacting on Pt(111) for three DFs are shown, comparing to the molecular beam experiments of Luntz et al.⁶⁴ A comparison to the experimental results of Cao et al.⁶⁰ is shown in Fig. S7 (ESI†). For the comparison to the experiment of Luntz et al.⁶⁴ the molecular beam parameters of Groot et al.¹⁴⁵ are used, while Cao et al.⁶⁰ have actually reported their molecular beam parameters. Fig. 8a shows that the B86SRP68-DF2 DF describes the experiment⁶⁴ with overall chemical accuracy, albeit that the energy shifts are larger than 4.2 kJ mol⁻¹ at the lowest energies. However, Fig. S7a (ESI†) shows that the experiments of Cao et al.⁶⁰ are not described to within chemical accuracy using the B86SRP68-DF2 DF. The SRP48 and MS-PBEl DFss are not able to describe either experiment to within chemical accuracy.

Fig. 8 Molecular beam sticking probabilities for D₂ reacting on Pt(111) for (a) the PBEα57-DF2, (b) B86SRP68-DF2, (c) SRP48 and (d) MS-PBEl DFs. Experimental results are shown in red,⁶⁴ 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.

The parameters describing the beams used in the experiment can be found in Table S3 (ESI†). MADs and mean signed deviations (MSDs) for all simulated molecular beam experiments considered here can be found in Table S9 (ESI†).

3.3 Initial-state resolved reaction probabilities

Initial-state resolved reaction probabilities will be presented for H₂ reacting on both Cu(111) and Ag(111). Fig. 9 shows fully initial-state resolved reaction probabilities for H₂ reacting on Cu(111), as obtained with the QD and QCT methods. At the level of degeneracy averaged reaction probabilities the agreement is very good (Fig. 9a and c), but the QD method predicts a slightly larger orientational dependence of the reaction probability (Fig. 9b). Fig. S8 (ESI†) shows a comparison of degeneracy averaged reaction probabilities obtained using the QD and QCT method for the 20 rovibrational states included in the QD calculations. The agreement between the QD and QCT method is very good for all states, though there are some differences for the J < 3 rovibrational states.

Fig. 9 Reaction probabilities computed for H₂ + Cu(111) with QD calculations (solid lines) and QCT calculations (dashed lines) for normal incidence using the B86SRP68-DF2 DF. Degeneracy averaged reaction probabilities computed for J = 1 for both the ground state and the first vibrationally excited state (a). Fully initial-state resolved reaction probabilities are shown for he m_J = 0 and 1 states belonging to J = 1 for both the ground vibrational and the first excited vibrationally state (b). Reaction probabilities computed for the J = 0 state for both the ground vibrational and the first excited vibrational state (c).

Fig. 10 shows degeneracy averaged initial-state resolved reaction probabilities for H₂ and D₂ reacting on Ag(111). A comparison is made to reaction probabilities extracted from associative desorption experiments assuming detailed balance.^77,78 The agreement with experiment seems best for D₂ and when using the MS-PBEl mGGA.¹² Note, however, that the P_deg(E, ν, J) extracted from experiments were not normalized, but simply assumed to saturate at unity, making it hard to perform an accurate comparison with experiment.

Fig. 10 Initial-state selected reaction probabilities P_deg(E, ν, J) computed for H₂ (D₂) + Ag(111) using the MS-PBEl (purple), SRP48¹⁹ (red), B86SRP68-DF2 (green) and the PBEα57-DF2 (blue) DFs are shown as a function of translational energy, comparing with values with values extracted from associative desorption experiments.^77,78 Results are shown for D₂ (ν = 0,J = 2) (a), D₂ (ν = 1,J = 2) (b), and H₂ (ν = 0,J = 3) (c).

3.4 E _1/2(ν, J) parameters

E _1/2(ν, J) parameters calculated for H₂ (D₂) + Cu(111) using Methods A1 and B1 are shown in Fig. 11 and 12, respectively, also comparing with experiment.⁴⁷ Using Method B1, the DFs that include non-local correlation qualitatively reproduce the rotational hindering observed experimentally for (ν = 0), i.e. the increase of E_1/2(ν, J) parameters with increasing J for low J before decreasing with increasing J. A third degree polynomial has been fitted to the calculated E_1/2(ν, J) parameters, which describes the dependence of the E_1/2(ν, J) parameters on J. The polynomials are shown, without the energy axis offset resulting from the fit to a third degree polynomial, in Fig. S9 (ESI†).

Fig. 11 E _1/2(ν, J) parameters as a function of J obtained using Method A1 for H₂ and D₂ reacting on Cu(111). Red circles represent the SRP48 values,⁶² green circles the MS-B86bl values,¹² blue circles represent the B86SRP68-DF2 values with the solid blue circles corresponding to QD calculations, magenta circles represent the PBEα57-DF2 values, and solid black circles represent experimental results.⁴⁷

Fig. 12 E _1/2(ν, J) parameters as a function of J obtained using Method B1 for H₂ and D₂ reacting on Cu(111). Red circles represent the SRP48 values,⁶² green circles the MS-B86bl values,¹² blue circles the B86SRP68-DF2 values with the solid blue circles corresponding to QD calculations, magenta circles represent the PBEα57-DF2 values, and solid black circles represent experimental results.⁴⁷

E _1/2(ν, J) parameters calculated for H₂ (D₂) + Au(111) using Method B2 are shown in Fig. 13, also comparing to experiment.⁷⁰

Fig. 13 E _1/2(ν, J) parameters as a function of J obtained using Method B2 for H₂ and D₂ reacting on Au(111), experimental values are shown in black.⁷⁰ Red circles represent the SRP48 values,¹⁴² green circles the B86SRP68-DF2 values, blue circles the PBEα57-DF2¹¹ values, purple circles the MS-PBEl¹² values and red triangles the PBE values.¹⁴²

Accompanying MAD and mean signed deviations (MSD) values of the computed E_1/2(ν, J) parameters from the experimental values are presented in Table S10 (ESI†) for both H₂ (D₂) + Cu(111) and H₂ (D₂) + Au(111).

3.5 Rotational quadrupole alignment parameters Cu(111)

In Fig. 14 we compare calculated rotational quadrupole alignment parameters using the QCT method to experimental ones for D₂ desorbing from Cu(111).⁴⁹ Note that a positive A⁽²⁾₀(ν, J) indicates a preference for parallel reaction, a negative value a preference for perpendicular reaction, and zero means the reaction proceeds independent of D₂ orientation. We observe only a monotonous increase of the rotational quadrupole alignment parameters with decreasing translational energy, indicating that at translational energies close to the reaction threshold molecules prefer to react in a parallel orientation.

Fig. 14 Panel (a) shows rotational quadrupole alignment parameter, A⁽²⁾₀(ν = 0,J = 11), and panel (b) shows the rotational quadrupole alignment parameter A⁽²⁾₀(ν = 1,J = 6) for D₂ reacting with Cu(111). Experimental results are shown in black.⁴⁹ Theoretical results using the QCT method are shown for the B86SRP68-DF2 (red), SRP48¹⁰ (green), PBEα57-DF2¹¹ (blue), and the MS-B86bl¹² (magenta) DFs.

3.6 Inelastic scattering of H₂ from Cu(111)

Vibrationally inelastic scattering probabilities, P(ν = 0,J→ν = 1,J = 3) for H₂ scattering from Cu(111) are shown in Fig. 15 for scattering from the initial J = 1, 3 and 5 states. Results were obtained using the QD method. The onset of the vibrational inelastic scattering probabilities is correlated with the onset of reactivity for each particular state. The DFs yield similar results for the initial (ν = 0,J = 5) rovibrational state, except that the SRP48 DF yields smaller vibrational excitation probabilities for E > 0.8 eV (panel c). For the (ν = 0,J = 1) and (ν = 0,J = 3) initial rovibrational states the differences are larger (panels a and b).

Fig. 15 Vibrationally inelastic scattering probabilities for P(ν = 0,J → ν = 1,J = 3). Shown are results for J = 1 (a), J = 3 (b), and J = 5 (c). QD results for the following DFs are shown: B86SRP68-DF2 (black), SRPsol63-DF2 (red), PBEα57-DF2 (green), optPBE-DF1 (blue), MS-B86bl (magenta), and SRP48²⁴ (orange).

Fig. 16 shows the ratio of rotationally elastic and inelastic scattering probabilities P(ν = 1,J = 0 → ν = 1,J = 2)/P(ν = 1,J = 0 → ν = 1,J = 0) computed with two DFs and comparing with experiment.⁵⁴ Note that both curves need to be shifted by 40 meV in order to overlap with the onset of the experimental curves measured for a surface temperature of 300 K.⁵⁴

Fig. 16 A comparison of the ratios of theoretical and experimental rotationally inelastic scattering probabilities P(ν = 1,J = 0 → ν = 1,J = 2)/P(ν = 1,J = 0 → ν = 1,J = 0) for H₂ impinging on Cu(111). Experimental results are shown in black⁵⁴ with the error bars representing maximum deviations in repeated measurements constituting estimated standard deviations. QD results for the SRP48 DF are shown in red¹⁰ and QD results for the B86SRP68-DF2 DF are shown in blue. Dashed lines constitute the calculated ratio's of rotationally inelastic scattering shifted by 40 meV.

3.7 Rovibrational state populations of H₂ and D₂ desorbing from Au(111)

Fig. 17 shows rovibrational state populations of H₂ and D₂ desorbing from Au(111). Here we plot ln[N/g_N(2J + 1)] versus the rotational energy, with N being the total population for each (ν, J) state (see eqn (13)) and g_N(2J + 1) being the statistical weight for rotational level J.⁷⁰ For D₂, g_N = 2 for even J and 1 for odd J; for H₂, g_N = 1 for even J and 3 for odd J. In such a plot a Boltzmann distribution will appear as a straight line.⁷⁰ Here we integrate eqn (13) up to E_max(ν, J).

Fig. 17 Rovibrational state populations of H₂ and D₂ desorbing from Au(111) are shown versus the rotational energy. Experimental results are shown in black,⁷⁰ theoretical results are shown in red for the SRP48 DF¹⁴² and in blue for the B86SRP68-DF2 DF. The straight lines represent Boltzmann distributions for the surface temperature of the experiment.

Fig. S10 (ESI†) shows the rovibrational state populations of H₂ desorbing from Au(111) as reported in by Shuai et al.⁷⁰ together with the values we calculate using eqn (13) with an upper integration limit of 5 eV to be in line with the procedure used by Shuai et al.⁷⁰ as outlined in a private communication.¹⁴⁶

Table 4 shows the ratio of the fluxes of molecules desorbing in the first excited and in the vibrational ground state for both H₂ and D₂ desorbing from Au(111). The ratios are calculated using eqn (15) and are solely based on the rovibrational states for which results are shown in Fig. 17. Again, here we integrate eqn (13) up to E_max(ν, J).

Table 4 The ratio of ν = 1 : ν = 0 molecules desorbing from Au(111) as measured in experiments⁷⁰ and computed with the SRP48 and B86SRP68-DF2 DFs

	H2	D2
Exp.^70	0.552	0.424
SRP48	0.250	0.473
B86SRP68-DF2	0.249	0.522

---

4 Discussion

Our aim has been to develop new SRP-DFs that include non-local correlation for the H₂ + Cu(111) system, and afterwards assess the transferability of these DFs to other systems. The reason for taking this approach, instead of fitting the DF to best reproduce experiments on all systems shown in the Results section, is that numerous experimental results are available for the reaction of H₂ and D₂ on Cu(111).^45–55 For the non-copper systems discussed in this work experimental results are sparse,^17,70 and there is discussion about the validity of the available parameters describing the molecular beams used in the experiments on H₂ (D₂) + Ag(111) and Pt(111).^19,60,64,65

The good agreement between different dissociative chemisorption experiments^45,46,51 on the reaction of H₂ (D₂) + Cu(111), and their resultant constraints for a to be developed SRP-DF, provides an opportunity to design the best performing SRP-DF for this system yet reported. The DF that best describes sticking probabilities obtained from dissociative chemisorption molecular beam experiments for H₂ (D₂) + Cu(111), will be considered the best performing DF. We choose this definition since our calculations are carried out using the BOSS model. From the literature it is known that the BOSS model works rather well for activated H₂ dissociation on cold metals.^{23,27–29,62,147}

Comparisons to experimental results obtained from associative desorption experiments will not be included in the assessment of the quality of the newly constructed SRP-DFs for two reasons. The first reason is that associative desorption experiments are carried out at high surface temperatures. Since we have carried out calculations using the BOSS model we neglect surface temperature effects. The second reason is that in obtaining state-specific information from associative desorption experiments requires the assumption of detailed balance, which is strictly speaking not applicable if an electronically non-adiabatic process is involved and energy exchange with the surface is allowed. Since neither process can be ruled out we feel it safer to base our judgement on the sticking experiments. We still discuss the valuable experimental results on associative desorption since they do provide insight into the reaction dynamics. However, as we will discuss, it is fraught with difficulty to make a direct quantitative comparison between calculations on dissociative chemisorption and associative desorption experiments without improving our dynamical model by incorporating phonons and ehp excitations, which is challenging to do.^148–150 Recently, this has been done for H₂ + Cu(111)⁷⁵ using ab initio molecular dynamics with electronic friction (AIMDEF) calculations employing the PBE¹¹⁵ DF. It is hard to draw firm conclusions from this work on the effect of electron–hole pair excitation as the statistical accuracy of the AIMDEF calculations is limited through the small number of AIMDEF trajectories.

4.1 Metal properties

Using a GGA DF for a theoretical description of gas surface dynamics is in most cases a compromise between a good description of the metal slab and a good description of the interaction of a molecule with the metal slab.¹⁵¹ Table S4 (ESI†) shows the calculated lattice constants for all DFs investigated in this work. Highly accurate results are only achieved using PBEsol¹¹⁷ and our previously developed mGGA MS-B86bl.¹² All our candidate SRP-DFs yield improved results over the vdW-DF1,⁹⁶ vdW-DF2,⁹⁷ and SRP48¹⁰ DFs. We generally see a large improvement in the calculated lattice constant if a DF has at least one component with a μ value closer to the second gradient expansion , as used in the PBEsol,¹¹⁷ MS-B86bl,¹² SRPsol63-DF2 and the B86SRP68-DF2 DFs.

When looking at the relaxations of the interlayer distance of the two top most layers relative to the bulk interlayer distance (Table 2) no clear trend can be discerned. The best performing DF is again MS-B86bl.¹² We do note that the SRP48 DF appears to produce top interlayer distances that on the whole are closer to the experimental values than DFs obtained combining GGA exchange DFs with vdW-DF2 non-local correlation. The reason for this is unclear.

4.2 Static PES properties

The experimental van der Waals well depth that has been measured for H₂ + Cu(111)⁵⁹ can provide us with a constraint when fitting a new SRP-DF for this system. The new SRP-DFs we present here are B86SRP68-DF2 and SRPsol63-DF2 (Table 1). The two original van der Waals DFs, namely vdW-DF1⁹⁶ and vdW-DF2,⁹⁷ yield wells that are too deep, especially for vdW-DF1,⁹⁶ as also found in earlier work¹⁵² (Fig. 3b, see also Table S8, ESI†). The two previous SRP-DFs for this system produce a very tiny (SRP48¹⁰) or no (MS-B86bl¹²) van der Waals well at all (Fig. 3b). Of the four DFs considered in Fig. 3b, the vdW-DF2⁹⁷ DF produces the best results. Other exploratory calculations carried out by us showed that using vdW-DF1⁹⁶ correlation yields van der Waals wells that are much to deep and too close to the surface, although they are not shown here. The better performance of vdW-DF2 correlation is most likely due to the large-N asymptote correction used in this DF.⁹⁷

Fig. 3a shows the same van der Waals well depth for the SRP-DFs that include vdW-DF2 correlation, and one SRP-DF that includes vdW-DF1 correlation. These four DFs are all chemically accurate with respect to the reactivity of H₂ on Cu(111), as can be seen from the MAD values in Fig. 4. The PBEα57-DF2 DF has originally been fitted to reproduce experiments for H₂ (D₂) + Pt(111),¹¹ and the optPBE-DF1¹¹⁴ DF has previously been shown to be chemically accurate H₂ (D₂) + Cu(111).⁴¹ The exchange part of the optPBE-DF1 DF was optimized to reduce intermediate range effects to avoid double counting when combining it with non-local vdW-DF1 correlation.¹¹⁴ It is clear that the choice of exchange functional greatly impacts the depth and position of the van der Waals well. The difference between the depths of the deepest and the shallowest van der Waals well obtained with the non-standard DFs using vdW-DF2 non-local correlation⁹⁷ (Fig. 3a, 23.4 meV) is greater than the difference between the depths of the vdW-DF1 and vdW-DF2 wells (Fig. 3b, 13.4 meV, see also Table S8, ESI†). It can also been seen that going from PBEα57-DF2 to SRPsol63-DF2 the calculated van der Waals well more closely resembles experiment.⁵⁹ The closest agreement with experiment is achieved using the B86SRP68-DF2 DF.

All van der Waals wells computed by us are tabulated in Table S8 (ESI†), also comparing with experimental van der Waals wells that have been reported for H₂ + Cu(111),^59,61 H₂ + Ag(111),⁹⁰ H₂ + Au(111)⁶¹ and H₂ + Pt(111).^99,110 With respect to the Cu(111) well depth the experimental results are in reasonable agreement with each other. However, the position reported by Harten et al.⁶¹ is somewhat closer to the surface. This difference in the Z dependence of the van der Waals well can be attributed to ambiguities in the bound state level assignments from the Feshbach resonances in the earlier experiment,¹⁰⁴ and it has been remarked that the results obtained later⁵⁹ are in fact also consistent with the earlier measurements.¹⁰⁴ Additionally we suspect that the reported van der Waals wells for H₂ + Ag(111)⁹⁰ and H₂ + Au(111)⁶¹ might possibly be too close to the surface¹⁰⁴ for the same reason. This drawback of using the RMSA technique⁹⁸ might be alleviated by redoing the potential inversion on the basis of the original data on Feshbach resonances with more advanced theoretical models,¹⁵³e.g., using an analysis in which the molecule–surface potential is not laterally averaged. In yet a different approach, instead of using the RMSA approach⁹⁸ Poelsema et al.¹¹⁰ presented a combined thermal energy atom scattering/thermal desorption spectroscopy (TEAS/TDS) study of the H₂ + Pt(111) system, obtaining van der Waals well geometries that were subsequently accurately reproduced by theory.¹¹ For the H₂ + metal(111) systems studied here it would certainly be advantageous if additional experimental data were to become available addressing the van der Waals interaction, using either a sophisticated analysis of results for RMSA studies or through combined TEAS/TDS studies. New experiments would allow for a better comparison between theory and experiment with respect to the predictions obtained by the inclusion of non-local correlation in DFT calculations on the systems addressed here.

We can however say that the B86SRP68-DF2 DF yields van der Waals well depths that agree well with experiments for the highly activated systems H₂+ Cu(111), H₂ + Ag(111) and H₂ + Au(111). For the weakly activated H₂ + Pt(111) system the B86SRP68-DF2 appears to somewhat underestimate the depth of the van der Waals well, while the SRP DF for this system (PBEα57-DF2) yields a value in the range of the experimental values.

4.3 Molecular beam sticking

4.3.1 Molecular beam sticking of H₂ (D₂) + Cu(111): QCT results. The fitting of the candidate SRP-DFs was done by reproducing six different sets of molecular beam experiments^45,46,51 for which the molecular beam parameters were taken from (the ESI† of) ref. 10 and 28. Fig. 4 compares results of three sets of molecular beam sticking probabilities, S₀, with results computed with the SRP48,¹⁰ B86SRP68-DF2, SRPsol63-DF2, PBEα57-DF2¹¹ and optPBE-DF1^41,117 DFs using the QCT method. We focus on these three sets of molecular beam experiments because for these enough experimental data points are available to perform a cubic spline interpolation. The quality of the DFs is assessed by computing MAD values, i.e. the mean distance along the incidence energy axis from the computed S₀ data point to the interpolated experimental data point with the same S₀ value.

It is clear from the total MAD values that all DFs evaluated in Fig. 4 are chemically accurate with respect to the three sets of molecular beam experiments for H₂ (D₂) + Cu(111) shown in Fig. 4, and that the agreement between theory and experiment is good for all molecular beam conditions. The theoretical results shown here were obtained using the BOSS model. The experiments of Michelsen et al.⁴⁵ and Rettner et al.⁵³ considered here employed a low surface temperature of 120 K, the experiment of Berger et al.⁴⁶ was reportedly done on a ‘cold’ surface, and from the literature it is known that the BOSS model works well for activated H₂ dissociation on cold metals.^23,26–29 Another advantage of fitting a SRP-DF to these sets of molecular beam experiments is that they cover both H₂ and D₂ for very different experimental conditions with respect to the nozzle temperature, the average collision energy, and the width of the velocity distributions.¹⁰

The B86SRP68-DF2 DF exhibits a MAD of 1.4 kJ mol⁻¹ for the experiments shown in Fig. 4, which is the lowest value obtained with the five DFs discussed. Of the DFs that include non-local correlation the B86SRP68-DF2 DF also performs best with respect to the calculated lattice constants, as can be seen in Table S4 (ESI†). In addition it performs best with respect to the shape and depth of the van der Waals well. We therefore select the B86SRP68-DF2 DF as the new, and most accurate, SRP-DF for the H₂ + Cu(111) system. From this point onward we will mainly focus on the results obtained with the newly selected B86SRP68-DF2 DF and the PBEα57-DF2¹¹ DF, and on how the performance of these two DFs compares to that of the original SRP48¹⁰ DF and of our previously developed mGGA DFs.¹²

In the ESI† one additional comparison to experiment is made in Fig. S5, concerning pure D₂ molecular beams.⁴⁵ This comparison highlights a limitation of assessing the quality of a candidate SRP-DF by computing the mean distance along the incidence energy axis from the computed S₀ value to interpolated experimental values. In this experiment the average translational energy does not monotonically increase with increasing nozzle temperature. Due to the high sensitivity of the sticking probability to the width of the velocity distribution of the molecular beam, the sticking probability also does not monotonically increase with the average translational energy. In Fig. S5 (ESI†) the application of our quality assessment strategy leads, in some cases, to large deviations with respect to the interpolated experimental results (note, however, that only for one data point and for one DF chemical accuracy was not achieved, see Fig. S5c, ESI†). Our quality assessment strategy works best if the reactivity increases monotonically with increasing average translational energy.

4.3.2 Molecular beam sticking of H₂ + Cu(111): QD results. Due to the computational expense of calculating molecular beam sticking probabilities using the QD method we only carried out QD molecular beam simulations for H₂ + Cu(111) for the B86SRP68-DF2 DF (Fig. 4e, f and 5). We have used the same methodology to carry out QD molecular beam simulations as in our previous work on Cu(211).³⁶ Overall, for the four experiments considered with both methods, the MAD value obtained with the B86SRP68-DF2 DF increases from a QCT value of 1.3 kJ mol⁻¹ to 1.6 kJ mol⁻¹ for QD. In the QCT calculations more rovibrational states are taken into account compared to the QD calculations (see Table S2, ESI†). In Fig. 6 we explicitly compare QD and QCT results calculated from the same set of rovibrational states, and compare those results to the QCT results shown in Fig. 4e, f and 5. Overall the agreement between QD and QCT sticking probabilities is very good when both are calculated from the same set of initial rovibrational states. There are however small differences when looking at narrow low average translational energy molecular beams (see e.g.Fig. 6a).

Fig. 6a and b show that for narrow low average translational energy molecular beams⁵¹ the QD method predicts slightly larger sticking probabilities compared to QCT sticking probabilities calculated from the same set of rovibrational states. This small gap between QD and QCT sticking probabilities based on the same set of rovibrational states is slightly bigger then what was obtained for H₂ + Cu(211) when using the SRP48 DF.³⁶ Though this suggests that quantum effects might play a small role in the dynamics, it is also possible that the slightly higher sticking probability predicted by QD is due to the underlying reaction probability curves for specific included rovibrational states showing more structure than for H₂ + Cu(211).³⁶ Since the molecular beam sticking probabilities are very small in Fig. 6a and b, they could be very sensitive to increased structure (and noise) in the underlying reaction probability curves. We shall further discuss the differences between fully initial-state resolved QD and QCT reaction probabilities in Section 4.6.

Fig. 6c and d show very good agreement between QD and QCT results, also for the QCT results that were based on more initial rovibrational states. Only for the highest nozzle temperature points in Fig. 6c (see Table S3, ESI†) taking into account more initial rovibrational states in the QD than those listed in Table S2 (ESI†) might be advisable. This is also the reason that QD now predicts a somewhat lower sticking probability than QCT in Fig. 4e. The QD and QCT results based on the same set of initial rovibrational states in Fig. 6c and d are however in good agreement at these high nozzle temperature points, indicating that for all but the lowest average translational energies H₂ + Cu(111) is well described quasi-classically.

Note that for a nozzle temperature of 2000 K we take into account all rovibrational states that have a Boltzmann weight >0.001. Highly excited rovibrational states, either with high ν, high J, or both, yield high reaction probabilities at low translational energies, therefor the effect of not taking into account these rovibrational states might be larger than expected from their Boltzmann weight. It is however computationally very expensive to take into account all initial rovibrational states that have a Boltzmann weight >0.001 at a nozzle temperature of 2300 K in the QD calculations, i.e. the nozzle temperature for the point in Fig. 6c where the QCT results and the QCT results based on the same, smaller, set of initial rovibrational states as the QD calculations diverge most. Doing this would nearly double the amount of wave packet calculations, and these additional calculations would also require larger basis sets.

4.3.3 Molecular beam sticking of D₂ + Pt(111). Fig. 8 shows a comparison of calculations on D₂ + Pt(111) to the molecular beam experiments of Luntz et al.⁶⁴ The PBEα57-DF2 DF was originally fitted¹¹ in order to reproduce these experiments. In Fig. 8b we see that the B86SRP68-DF2 DF can also describe these experiments with overall chemical accuracy (MAD = 3.1 kJ mol⁻¹), although the agreement with experiment at low incidence energies is just shy of being chemically accurate, which would adversely affect the extraction of the minimum barrier height for this system. The SRP48 DF also describes the experiments with overall chemical accuracy, but at low E the agreement with experiment is really poor (Fig. 8c). The MS-PBEl DF does not agree with experiments to within overall chemical accuracy (Fig. 8d).

The discrepancy between theory and experiment at low incidence energies for the SRP48 and B86SRP68-DF2 DFs most likely arises because these DFs exhibit a too high early barrier to reaction at the top site minimum barrier geometry (see Table S7, ESI†). These two DFs also do not posses the double barrier structure for the t2b site that the PBEα57-DF2 DF predicts. The MS-PBEl functions does predict a double barrier structure for the t2b site. The early t2b barrier predicted by the MS-PBEl DF is however very high when compared to results obtained with the other DFs, which is most likely the root cause of its poor performance for this system.

The molecular beam experiments of Cao et al.⁶⁰ are not as well described by the B86SRP68-DF2 DF, as can be seen in Fig. S7 (ESI†). However we note that the increase of the MAD value in going from Fig. 8a to Fig. S7a (ESI†) is similar in size to what was reported for the PBEα57-DF2 DF⁶⁵ (see Table S9, ESI†). Earlier work from our group has shown that the experimental results of Luntz et al.⁶⁴ and Cao et al.⁶⁰ are in good agreement with each other for the lower incidence energies but somewhat diverge for the higher incidence energies.⁶⁵ The possible origins of the discrepancy between these two sets of experimental data are discussed in ref. 65. Here it was surmised that at high average incidence energies the reaction probabilities of Cao et al.⁶⁰ are most likely somewhat underestimated compared to the results of Luntz et al.⁶⁴ because the average incidence energies are somewhat underestimated in the experiment of Cao et al.⁶⁰

4.3.4 Molecular beam sticking of D₂ + Ag(111). We now make a comparison to molecular beam experiments on D₂ + Ag(111).¹⁴⁴ Even though silver is only one row below copper in the periodic table, the SRP48 DF that was fitted to reproduce experiments on Cu(111) was not able to describe experiments on Ag(111) with chemical accuracy.¹⁹Fig. 7 shows the computed S₀ for D₂ + Ag(111). Hodgson and coworkers¹⁴⁴ have reported translational energy distributions that were symmetric in the energy domain. In our view, and as discussed in previous work from our group,¹⁹ the symmetric translational energy distributions are somewhat unphysical. Therefore we opted to use the molecular beam parameters of pure D₂ reacting on Cu(111) reported by Auerbach and coworkers,⁴⁵ which likewise describe beams that are narrow in translational energy. Here we see that the PBEα57-DF2 and B86SRP68-DF2 DF are similar in accuracy to our previously developed MS-PBEl mGGA DF,¹² and that these three DFs predict a reactivity just shy of chemical accuracy when compared to the molecular beam experiment of Hodgson and coworkers (see Table S9, ESI†). Note that the PBEα57-DF2 DF seems to perform worse than the other three DFs at the lowest translational energy. All three SRP-DFs yield an overall performance that is better than that of SRP48.

Although we are not yet able to describe the molecular beam experiment with chemical accuracy, the improvement of the DFs that include non-local correlation over the SRP48 DF again suggests that non-local correlation is an important ingredient for constructing SRP-DFs describing H₂ + metal systems that incorporate GGA exchange. The MS-PBEl mGGA¹² does not include non-local correlation but performs similarly well as the GGA based SRP-DFs that do include non-local correlation. Perhaps further improvement is possible by adding non-local correlation to the MS-PBEl DF.

4.4 Associative desorption

4.4.1 Comparing to experimental E₀(ν, J) parameters. When comparing to experimental E₀(ν, J) parameters by calculating E_1/2(ν, J) parameters from calculated degeneracy averaged reaction probabilities we effectively try to model an associative desorption experiment as an initial-state resolved dissociative chemisorption experiment. Additionally our QCT results are obtained using the BOSS model, while the associative desorption experiments necessitate high surface temperatures. The J dependence of the calculated E_1/2(ν, J) parameters is a measure of how accurate the reactivity of the individual rovibrational states is described relative to each other. The trend is here to be understood as the dependence of the E₀(ν, J) parameters on J, which we can visualize using a third degree polynomial fit to the calculated results (see e.g. Fig. S9, ESI†).

In this work we have not included MDEF calculations in which the effect of ehp excitations is modeled as a classical friction force. In our previous work³⁶ the MDEF method shifted the E_1/2(ν, J) parameters to slightly higher values since, again, we model an associative desorption experiment using a dissociative chemisorption calculations. Including the effect of ehp excitations in the dynamics here then has a similar effect in both cases, and shifts the E_1/2(ν, J) parameters to slightly higher energies. However, as mentioned in our previous work,³⁶ if ehp excitations are important then assuming detailed balance to extract degeneracy averaged reaction probability curves is, strictly speaking, not correct. More specifically, we would expect that if we applied electronic friction to the simulation of an associative desorption experiment in the manner discussed here, the predicted translational energy distributions would shift to higher energies as opposed to lower energies as expected in a direct simulation of associative desorption. In other words: extracting E_1/2(ν, J) parameters from dissociative chemisorption calculations applying electronic friction would shift our E_1/2(ν, J) parameters to higher energies instead of the expected lower energies, because the effective barrier would go up. However, we note that the MDEF calculations in our previous work³⁶ only shifted the trend in E_1/2(ν, J) parameters to slightly higher values on the energy axis and did not influence the observed trend in their J-dependence. This is in accordance with very recent direct simulations of associative desorption of H₂ from Cu(111)⁷⁵ using AIMD and AIMDEF calculations and the PBE DF, which find little effect of ehp excitations modeled at the local density friction approximation (LDFA)¹⁵⁴ level.

In Section 4.10 we will further discuss how a combined analysis of dissociative chemisorption and associative desorption experiments might be used in the future to determine a possible fingerprint of non-adiabatic effects.

4.4.1.1 E _1/2(ν, J) parameters Cu(111). The comparison of the measured E₀(ν, J) parameters values with the E_1/2(ν, J) parameters computed using Method B1 is shown in Fig. 12. The E_max(ν, J) parameters needed in this method (see eqn (S4), ESI†) were taken from Tables S4 and S6 of ref. 47. The comparison of the experimental E₀(ν, J) parameters with theoretical E_1/2(ν, J) parameters extracted using Method A1 is presented in Fig. 11. In this method we used A = 0.325 for H₂ and A = 0.513 for D₂, as obtained in ref. 47.

Table S10 (ESI†) presents MAD and MSD values obtained with both methods. The following conclusions can be drawn: for almost all DFs Method A1 (based fully on experiments) yields the lowest MAD and MSD values. The only exception occurs for H₂ + Cu(111) when using the SRP48 DF, for which Method B1 gives a lower MAD value (41 meV), although the difference with Method A1 (43 meV) is quite small. From this point of view, Method A1 works better.

The use of Method A1 for H₂ would seem to yield conclusions that are more consistent with the conclusions from the comparisons of the sticking probabilities. Specifically, the mGGA DF and the DFs containing non-local van der Waals correlation all perform better than SRP48 DF for H₂. Nonetheless, SRP48 performs best for D₂, regardless of whether Method A1 or B1 is used. We also note that the better behavior of the other DFs in procedure A1 is to some extent suspect due to the rather low A value employed for H₂ (0.325). This A value is much lower than the A value extracted for D₂ in Method A1 (0.513). This may well be a simple artifact resulting from the method followed: the A value determined for D₂ is likely to be more accurate because the sticking experiments were done for a kinetic energy up to 0.83 eV,⁴⁵ whereas the sticking experiments for H₂ only went up to about 0.5 eV.⁵¹ This suggest that the A-value for D₂ (Fig. S1c, ESI†) is much more accurate than for H₂ (Fig. S1a, ESI†), and it is not clear why the A value for H₂ should differ much from it. Furthermore, the A value established for H₂ in Method A1 is much lower than the A-values established for H₂ in Method B1 (see also Fig. S1a and c, ESI†), while for D₂ the A values extracted with the two methods resemble each other, and the A values extracted with Method B1 for H₂, much more (see Fig. S1b–d, ESI†).

Finally, we note at this stage the E_1/2(ν, J) parameters computed using Method A1 do, in general, not reproduce the subtle trend found experimentally that the E₀(ν, J) parameters first increase somewhat with J (see Fig. 11) (attributed to rotational hindering^45,47,51).

The E_1/2(ν, J) parameters calculated using Method B1 (especially the ones calculated with DFs incorporating non-local van der Waals correlation) better reproduce this subtle rotational hindering effect (Fig. 12). The reasonable performance of Method B1 is also clear from Fig. S9 (ESI†), which presents a third degree polynomial fit to the E_1/2(ν, J) parameters as a function of J obtained using the B1 method. The polynomial fits are shown without the energy axis offset. DFs that include non-local correlation reproduce the subtle rotational hindering effect, DFs that do not include non-local correlation do not or hardly show rotational hindering. Additionally, for H₂, the agreement with the experimental dependence of the E₀(ν, J) parameters on J improves when using the QD method for vibrationally exited molecules.

That DFs including non-local correlation better reproduce the subtle rotational hindering effect with the use of Method B1 is wholly due to the rovibrational state dependence of the E_max(ν, J) parameters. The extracted degeneracy averaged reaction probabilities in fact monotonically increase with increasing J (Fig. S8, ESI†), and this is true for all DFs used in this work. Reproducing rotational hindering based on these degeneracy averaged reaction probability curves is therefore not possible when selecting the same A value for all rovibrational states, as done with Method A1.

The E_1/2(ν, J) parameters calculated using both Methods A1 and B1 do reproduce the clear trend (Fig. 11 and 12) that at high J the E₀(ν, J) parameters decrease with J (attributed to energy transfer from rotation to the reaction coordinate as the rotational constant decreases when the molecule stretches to reach the dissociation barrier^25,36,126).

We are aware of one single PES that does reproduce the rotational hindering effect as observed in the experiment, namely the LEPS PES³⁹ used by Dai and Light³⁷ for six-dimensional QD calculations. As discussed in the ESI,† we have investigated whether the rotational hindering observed by Dai and Light³⁷ could be due to their QD calculations being unconverged. We find that this is not the case, but that the observed difference with our calculations using the best SRP-DF (B86SRP68-DF2) only arises for J < 3, suggesting that the rotational hindering effect is very subtle (see Fig. S4, ESI†). Further research is needed to check whether the difference between the calculations could be due to the LEPS fit³⁹ being inaccurate, the underlying electronic structure calculations being unconverged,^39,40 or both. See the ESI† for further details.

Reported density functional molecular dynamics (DFMD) calculations that include surface motion have shown that at low translational energies surface temperature effects somewhat increase the reactivity of D₂ reacting on Cu(111).^23,62 Increased reactivity at low translational energies might lower calculated E_1/2(ν, J) parameters for both Methods A1 and B1. We have extracted E_1/2(ν, J) parameters using both Methods A1 and B1 from Nattino et al.⁶² The results are shown in Fig. 18. The reported data was obtained using the SRP48 DF.⁶² There were only three rovibrational states for which a data point was available at a translational energy higher than the maximum kinetic energy sensitivity of the experiment.⁴⁷ To validate the obtained E_1/2(ν, J) parameters we have also extracted E_1/2(ν, J) parameters using the same methodology for the reported logistics function fits to BOSS direct dynamics data.⁶² The good agreement for both Methods A1 and B1 between our SRP48 E_1/2(ν, J) parameters and those obtained from logistics function fits to BOSS direct dynamics data validates this comparison. For both the A1 and B1 method the single point for (ν = 0,J = 11) suggests that a small decrease can be expected from allowing the surface to move, for the (ν = 1,J = 6) point this is not so clear. There is however a dramatic decrease of the E_1/2(ν, J) parameter for (ν = 1,J = 4). This is a clear indication that, at least for low J, taking into account surface motion leads to lower E_1/2(ν, J) parameters, and this might thus be partly responsible for the observed rotational hindering. The decrease of the E_1/2(ν = 1,J = 4) value with the introduction of surface motion is less pronounced when using the A1 method.

Fig. 18 E _1/2(ν, J) parameters as a function of J for D₂ reacting on Cu(111). Closed symbols pertain to E_1/2(ν, J) parameters calculated using the A1 method, open symbols refer to the B1 method. Experimental results are shown in black,⁴⁷ SRP48 results calculated in this work are shown in red. DFMD (green) and BOSS (blue) results obtained using the SRP48 DF have been extracted from ref. 62.

Note however that the E₀(ν, J) parameters are extremely sensitive to the quality of the logistics function fits. The agreement between our SRP48 E_1/2(ν, J) parameters and those obtained from logistics function fits to BOSS direct dynamics for which no data point existed at a translational energy higher then the maximum kinetic energy sensitivity to which the experiment was sensitive was not so good.

The above observations warrant the following tentative conclusions: within the BOSS approximation, the mGGA DF and the DFs containing non-local correlation perform best for sticking. However, assuming Method B1 to be best, the SRP48 performs best for associative desorption. For associative desorption and again within the BOSS approximation, the two different methods (A1 and B1) for extracting the E_1/2(ν, J) parameters describing the reaction probabilities extracted from associative desorption experiments yield rather different results. It follows that, at this stage, the associative desorption experiments are not as useful as the sticking experiments for assessing the accuracy of theory. Hopefully this can be changed in future by taking into account surface atom motion and ehp excitation in the theory, and by computing associative desorption fluxes directly from theory, as was recently done using the PBE DF by Galparsoro et al.⁷⁵

4.4.1.2 E _1/2(ν, J) parameters Au(111). The comparison of the measured E₀(ν, J) values with the E_1/2(ν, J) parameters computed using Method B2 is shown in Fig. 13. The E_max(ν, J) parameters needed in this method (see eqn (S4), ESI†) have been obtained from a private communication,¹⁴⁶ and have been plotted in Fig. S2 (ESI†) against the E₀(ν, J) parameters extracted from the measurements.⁷⁰ We note that the E_max(ν, J) parameters typically equal for the ν = 0 states and for the ν = 1 states, where W = 0.31 eV for ν = 0 and 0.29 eV for ν = 1⁷⁰ (see Fig. S2, ESI†). This means that Method B1 will yield unreliable values for the E_1/2(ν, J) parameters. It also means that Method B2 relies heavily on extrapolation in the method to determine the A values by anchoring the measured reaction probabilities to the reaction probabilities computed at E_max(ν, J). One should therefore exercise extreme caution when comparing the theory to experiment.

Table S10 (ESI†) shows the MAD end MSD values obtained with Method B2 for H₂ + Au(111). The following tentative conclusions can be drawn: the E_1/2(ν, J) parameters computed with the mGGA DF tested here and with all DFs employing van der Waals correlation substantially overestimate the measured E₀(ν, J) parameters, with MAD values of approximately 0.1 eV for H₂. The PBE DF would appear to perform best, with a MAD value of 46 meV for H₂. This conclusions agrees with the observation that PBE reaction probabilities¹⁴² allowed better fits of the measured time-of-flight spectra of H₂ and D₂ desorbing from Au(111)⁷⁰ than the SRP48 reaction probabilities.

A caveat here is that the experiment was performed with a surface temperature of 1063 K, while all calculations were performed using the BOSS model. Allowing surface motion during the dynamics would lead to broadening of the sticking probability curves.^{23,29,155,156} A higher sticking probability at lower translational energies could potentially lower the theoretical E_1/2(ν, J) parameters. Additionally, we have performed our calculations on an unreconstructed Au(111) surface because the surface unit cell of reconstructed Au(111) is at present too big to map out a full PES using DFT calculations. Earlier work in our group indicated that the barriers to H₂ dissociation are somewhat higher on the reconstructed Au(111) surface and that using an unreconstructed surface might lead to the underestimation of dynamical barrier heights by about 50 meV (∼1 kcal mol⁻¹),¹⁴² which would lead to slightly higher E_1/2(ν, J) parameters and therefore to increased disagreement with the measured values.

We are not able to resolve the contradiction posed by Shuai et al.⁷⁰ that the vibrational efficacies computed with the SRP48 DF are in good agreement with experiment (as may be derived from Fig. 13) but that the ratio of desorbed molecules in the vibrational ground state versus the vibrationally excited state is not (see Table 4). The good reproduction of the J dependence of the E₀(ν, J) parameters by theory (Fig. 13) suggests that the reactivity of the individual rovibrational states relative to each other is accurately described by theory as long as states are considered with the same vibrational level. Previously reported experiments implied that the recombination of H₂ on Au(111) is coupled to the electronic degrees of freedom of the metal.^157–160 In line with Shuai et al.⁷⁰ we think that non-adiabatic effects together with surface motion effects and surface reconstruction represent the most likely causes for the lower translational energy distributions of the desorbing H₂ (ν = 0) molecules compared to theory. If molecular beam dissociative chemisorption experiments on a reasonably cold surface would become available for this system (like for H₂ + Ag(111)) this would allow for a more direct comparison to experiment of QCT and MDEF calculations. Molecular beam adsorption experiments would also allow us to check if the absolute reactivity predicted by the DFs shown here is in agreement with experiment. Therefore, at present, we cannot corroborate or refute the conclusion reached by Shuai et al.,⁷⁰ namely that the experimentally observed lower translational energy distributions compared to theoretical predictions (see Fig. 1 of Shuai et al.⁷⁰) is most likely due to ehp excitations in the desorption dynamics.

4.4.2 Rovibrational state populations of H₂ and D₂ desorbing from Au(111). Fig. 17 shows the rovibrational state populations of H₂ and D₂ desorbing from Au(111). Note that we have consistently applied the normalization procedure outlined by eqn (14) to the objects shown in Fig. 17. In the case of (ν = 0) the populations deviate from the slope set by the Boltzmann distributions at the surface temperature of 1063 K indicating that rotationally excited molecules are more likely to adsorb.⁷⁰ The populations of vibrationally excited molecules also lie on gentler slope than implied by the Boltzmann distribution and are consistently higher than would be obtained with the Boltzmann distribution indicating that vibrationally excited molecules are more likely to adsorb. Both these observations are in line with the Polanyi rules for a late barrier system like H₂ + Au(111).^161,162

Fig. S10 (ESI†) shows the state distributions of molecules desorbing from Au(111) as reported by Shuai et al.⁷⁰ in their Fig. 2 together with experimental results calculated by us using eqn (13), an upper integration limit of 5 eV, and our normalization procedure, which were the boundary conditions and integration parameters suggested to us by Shuai et al.^70,146 As can be seen from Fig. S10 (ESI†) our “experimental results” calculated using eqn (13) and the error function fits reported in ref. 70 are in good agreement with the experimental results reported by Shuai et al.,⁷⁰i.e. the calculated curves have the correct shape and can, to a good approximation, be superimposed on one another, provided that they are shifted by a constant value as explained in the caption of Fig. S10 (ESI†). However, we are not able to reproduce the normalization strategy employed by Shuai et al.⁷⁰ Also when we shift the SRP48 (ν = 0) and (ν = 1) curves by the same value we cannot exactly superimpose our computed SRP48 results with their computed SRP48 results for both vibrational states simultaneously. It is not clear to us how this discrepancy arises.

The relative populations for the (ν = 0) and (ν = 1) rovibrational states is not affected by this discrepancy. The difference between our work and Shuai et al.⁷⁰ with respect to the ν = 1 : ν = 0 ratios arises because we use E_max(ν, J) as the upper integration limit in eqn (13). Shuai et al.⁷⁰ have used 5 eV as the upper integration limit.¹⁴⁶ The only reason we choose E_max(ν, J) as the upper integration boundary is because the reported error function fits are only reliable below E_max(ν, J), for some rovibrational states the error function fits can yield sticking probabilities much larger than unity for high kinetic energies. The ratios we calculate are shown in Table 4. We note that when we use the upper integration limit of 5 eV, we reproduce the ν = 1 : ν = 0 ratios reported by Shuai et al.⁷⁰ In our view only integrating up to E_max(ν, J) is a more fair way of calculating the N(ν, J) populations, though on the scale of Fig. 17 the difference between integrating to E_max(ν, J) or 5 eV would not be visible. Note that the overwhelming majority of the area under the Gaussian fits to the time-of-flight curves, as reported in Tables S1–S4 of ref. 70, lies well below E_max(ν, J). Note also that we calculate the ν = 1 : ν = 0 ratios only using the rovibrational states shown in Fig. 17, which is the same set of rovibrational states used by Shuai et al.⁷⁰

In Fig. 17 the difference between the desorbing populations computed with the SRP48 and B86SRP68-DF2 SRP-DFs is minimal. The agreement between theory and experiment is best for D₂, although the qualitative agreement between theory and experiment is reasonable for both H₂ and D₂. It can be seen in Table 4 that there is only a reasonable agreement for D₂ with respect to the ν = 1 : ν = 0 population ratios. The theoretical ν = 1 : ν = 0 population ratio for H₂ is however too low, a result similar to what was reported by Shuai et al.⁷⁰ The difference between theory experiment can perhaps be explained by the experimental time-of-flight distributions being much broader than the theoretical ones, see Fig. 1 of ref. 70. Taking into account surface motion in the theoretical calculations might well improve the agreement with experiment with respect to both the rovibrational state distributions of desorbing molecules as well as the ν = 1 : ν = 0 population ratio.

4.4.3 Initial-state resolved reaction probabilities Ag(111). In the case of molecular beam sticking results for D₂ + Ag(111) the MS-PBEl function performed similarly well as the GGA-exchange based SRP-DFs that include non-local correlation. This is most likely due to the slightly earlier barriers to reaction predicted by the MS-PBEl DF, leading to less promotion of reaction by vibrationally excited H₂,¹² as expected from the Polanyi rules.^161,162 From this argument it follows that although the molecular beam reaction probabilities are similar the reactivity of individual rovibrational states should be different.

In Fig. 10 we see that, especially for D₂, the MS-PBEl mGGA DF has the best agreement with the initial-state selected reaction probabilities extracted from the associative desorption experiments. The good performance of the MS-PBEl DF can be explained by the slightly earlier barriers, as discussed in our previous work.¹² The DFs that include non-local correlation do not show such a large improvement over the SRP48 DF, while they do for sticking.

Without new experimental work for this system, especially molecular beam experiments covering a wide range of translational energies and nozzle temperatures, it will be difficult to further improve the theoretical description of this system. Additional experiments (e.g. a molecular beam sticking experiment on D₂ seeded in H₂ and going op to a translational energy of 0.8 eV as done for H₂ + Cu(111)⁴⁵) would also allow us to assess more accurately if the dynamics predicted by the MS-PBEl mGGA or the dynamics predicted by the GGA based SRP-DFs that do and do not include non-local correlation are more in line with experimental observations.

4.4.4 Rotational quadrupole alignment parameters: H₂ + Cu(111). In Fig. 14 we compare calculated rotational quadrupole alignment parameters for D₂ to experimental ones measured for D₂ desorbing from Cu(111).⁴⁹ We observe only a monotonic increase of the rotational quadrupole alignment parameters with decreasing translational energy, indicating that at translational energies close to the reaction threshold molecules prefer to react in a parallel orientation. This is in line with what has been reported in the literature for H₂ and D₂ associatively desorbing from Cu(111)^26,49 and Cu(100),^18,25 and can be explained by invoking a static effect of orientational hindering in which rotating molecules scatter when their initial orientation does not conform to the lowest barrier geometry.²⁵ With increasing translational energy the rotational quadrupole alignment parameter approaches zero since all molecules, irrespective of their orientation, will have enough energy to react.⁴⁹ The experimental trend is reproduced by all DFs shown in Fig. 14a and b, though the calculated values are higher than the experimental values. The theoretical results presented here have been obtained within the BOSS approximation. francescorqa have shown that incorporating surface motion in the dynamics using the DFMD technique leads to better agreement with experiment, for D₂ + Cu(111).

Note that the B86SRP68-DF2, PBEα57-DF2¹¹ and MS-B86bl¹² DFs are in good agreement with each other for both the (ν = 0,J = 11) state (Fig. 14a) and the (ν = 1,J = 6) state (Fig. 14b), but that these three DFs predict slightly higher rotational quadrupole alignment parameters than the SRP48 DF.²⁶ Given that the SRP48 rotational quadrupole alignment parameters were decreased, but still somewhat too large when surface atom motion was introduced,²⁶ the present results suggest that the SRP48 DF yields the best description of this observable.

4.5 Inelastic scattering of H₂ from Cu(111)

In this section we will discuss inelastic scattering results for H₂, obtained with QD. We start with the vibrationally inelastic scattering results for H₂ shown in Fig. 15. We specifically show QD results since the previously voiced expectation that vibrationally inelastic scattering should be well described using the QCT method for translational energies above the lowest barrier to reaction¹⁶³ has been shown not to hold.²⁴ Here we discuss the inelastic scattering probability P(ν = 0, J → ν = 1, J = 3) for three different initial J states. Panels a and b suggest that for (J = 1) and (J = 3) and for the DFs that use non-local correlation the vibrationally inelastic scattering probability is correlated with the depth of the van der Waals well (see Table S5, ESI†). A deeper van der Waals well is correlated with higher vibrationally inelastic scattering probabilities. Fig. 15c shows that for (J = 5) all DFs yield vibrational excitation probabilities in reasonable agreement with each other. Bringing TOF spectra for vibrational excitation from (ν = 0,J) to (ν = 1,J = 3) in better agreement with experimental results would require a substantial increase of the vibrational excitation probabilities computed with the SRP48 DF (by a factor of 2–3),^24,163 which is obtained with none of the DFs tested. We conclude that better agreement with experiment probably requires a different dynamical model, as suggested also by earlier work.^24,163

From the computed ratio of rotationally inelastic scattering probabilities P(ν = 1, J = 0 → ν = 1, J = 2)/P(ν = 1, J = 0 → ν = 1, J = 0) shown in Fig. 16 it is clear that the B86SRP68-DF2 DF performs not as good as the SRP48 DF.¹⁰ The shifted SRP48 curve follows the experiment more closely. Both curves need to be shifted by 40 meV in order to better overlap with the experiment performed using a surface temperature of 300 K.⁵⁴ The overlap with experiment of the shifted curves only holds until 0.14 eV, but the experimentalists noted that at higher energies the measurements became more difficult.⁵⁴ It is rather surprising that both computed ratios need to be shifted by roughly the same amount in order to overlap with experiment, since the SRP48 DF overestimates the initial sticking probability in molecular beam experiment while the B86SRP68-DF2 DF does not. From the literature it is also known that including surface motion during the dynamics might lead to broadening, and an earlier onset, of inelastic scattering probabilities.^23,29 We speculate that allowing surface motion and ehp excitation during the dynamics might obviate the need for the shift in order to superimpose the calculated curves with experiment.

4.6 QD vs. QCT for H₂ + Cu(111)

In Fig. 9 initial state-resolved reaction probabilities are shown calculated using the B86SRP68-DF2 DF. Here the (J = 0) and (J = 1) state for both the vibrational ground state and the first vibrationally excited state are shown because the differences between the QD and QCT method are most prevalent for the low lying rotational states. The QD reaction probability curves show more structure than was shown for H₂ + Cu(211),³⁶ but the agreement between the QD and QCT method for degeneracy averaged reaction probabilities (Fig. 9a and c) is still very good. From the comparison between QD and QCT degeneracy averaged reaction probabilities shown in Fig. S8 (ESI†) it can be seen that the differences between the QD and QCT method get smaller with increasing J for J > 3, though small differences remain even for high J states.

The biggest difference between the QD and QCT method are observed in Fig. 9b. Here we show fully initial-state resolved reaction probabilities, thereby distinguishing between a 'cartwheeling' molecules rotating in a plane parallel to the surface normal (m_J = 0) and a 'helicoptering' molecules rotating in a plane perpendicular to the surface normal. In line with our previous work, we observe that QD predicts a slightly larger preference for molecules reacting parallel to the surface. The rovibrational states shown in Fig. 9 are the same rovibrational states shown in Fig. 2 of our previous work for H₂ reacting on Cu(211),³⁶ in which the agreement for degeneracy averaged reaction probabilities was nearly perfect.

In recent experimental work kaufmannCu211 have reported a previously unobserved “slow reaction channel” for H₂ associatively desorbing from Cu(111) and Cu(211). In this channel, the reaction could be facilitated by trapping on the surface and distortion of the surface due to thermal motion forming a reactive site.⁴⁷ Even though our PES now contains a van der Waals well that might facilitate trapping during the reaction dynamics, we do not yet see evidence of the recently reported slow reaction channel for H₂ + Cu(111).⁴⁷ The translational energy range used in our calculations overlaps with the translational energies at which the slow channel reactivity was observed.⁴⁷ We can therefore rule out quantum effects (like tunneling, see ref. 36) during the dynamics as the origin of this slow reaction channel for H₂ + Cu(111), as we did before for H₂ + Cu(211).³⁶ We therefore propose, as done earlier for H₂ + Cu(211),³⁶ that the slow reaction channel reported by Kaufmann et al.⁴⁷ originates from the very high surface temperature of 923 K used in the associative desorption experiments. Presently it is not possible to take surface motion explicitly into account in QD calculations, and it is challenging to do so in QCT calculations.^148–150 Galparsoro et al.⁷⁵ likewise did not yet find evidence for the slow reaction channel in their AIMD calculations.

4.7 Overall description of systems

When looking at the H₂ + Ag(111) and H₂ + Au(111) systems considered in this work together, one stark realization is that further development of chemically accurate DFs for H₂ reacting on transition metal surfaces is still heavily stymied by a lack of experimental data. This is bad news as presently semi-empirical DFT seems to be the only path to extracting chemically accurate information on barriers to reaction. Relying on non-empirical constraints on DF design is not yet feasible, as illustrated by the poor performance of the SCAN¹⁶⁴ DF for H₂ + Cu(111).¹² Additionally, taking another step upwards on Jacob's ladder from a GGA or mGGA towards hybrid DFs is computationally very expensive, if not prohibitively so.

For both the activated and non-activated reactions of H₂ on transitions metals there is now only a single well studied system, namely H₂ + Cu(111) (and maybe H₂ + Pt(111)⁶⁵). What we mean by well studied is that there should be different kinds of well described experiments. For example, a combination of an associative desorption experiment and a dissociative chemisorption experiment should be available, or sticking probabilities for normal and off-normal incidence. It is also critical that the experimental conditions are described accurately.^{15,19,20,65,142}

Without new and detailed experiments on, at least, the related H₂ + Pd(111) and H₂ + Ag(111) or H₂ + Au(111) systems it is not possible to grasp the overarching trends in reactivity imposed by the position of these metals in the periodic table. In many aspects we are dancing in the dark with respect to DF design. The consequence of this is that, presently, theory can only provide models with limited predictive power.

4.7.1 H₂ (D₂) + Cu(111). The H₂ + Cu(111) system is the best described system of the ones treated here. Low surface temperature molecular beam sticking experiments are very accurately described using the BOSS model. The associative desorption experiments are however less well described by the new SRP-DFs that have been designed to reproduce low surface temperature molecular beam experiments with calculations using the BOSS model. E_1/2(ν, J) parameters obtained from reported DFMD data⁶² suggest that, at least for this system, better agreement with experiment can be attained by including the surface degrees of freedom in the dynamics. On the other hand it appears that the agreement with E₀(ν, J) parameters measured in an associative desorption experiment is also increased for vibrationally excited molecules when using the QD method.

The large amount of experimental^{45–61,165,166} and theoretical work^{10,12,18,21–44,155,167,168} available have allowed us to get the best description of this system so far. We can however not yet point to one DF that is clearly the best DF for this system. Currently two DFs compete for being the best DF for this system, namely B86SRP68-DF2 and MS-B86bl.¹² The latter has a better description of the metal and might therefor be better when looking at diffraction probabilities. The MS-B86bl DF however misses any description of van der Waals forces. In all our simulations for the H₂ + Cu(111) system the B86SRP68-DF2 and MS-B86bl DFs perform similarly well. With the information available now one might argue that the MS-B86bl DF is the best DF for this system since its description of the metal is much better than provided by the B86SRP68-DF2 DF, and because the effect of including non-local correlation is only apparent when calculating E_1/2(ν, J) parameters for this system. The MS-B86bl DF is however not transferable to weakly activated systems. In our view, the DFs that are more generally applicable, i.e. the B86SRP68-DF2 and PBEα57-DF2 DFs, are currently the best DFs. A good next step could be to use non-local correlation together with the MS-B86bl DF.

4.7.2 D₂ + Ag(111). For the D₂ + Ag(111) system it is more difficult to assess the quality of our theoretical description due to the lack of well defined molecular beam parameters.¹⁹ DFs that use GGA-exchange and non-local correlation, and the MS-PBEl DF predict roughly similar molecular beam sticking probabilities. The comparison to the initial-state resolved reaction probabilities suggests that the MS-PBEl DF performs best due to its slightly lower and earlier barriers.¹² Since the MS-PBEl DF has a better description of the metal, a better description of the initial-state resolved reaction probabilities, and performs similar to the other candidate SRP-DFs concerning molecular beam sticking, one might argue that the MS-PBEl DF is currently the best DF for this system. As said before, the best DF should also exhibit transferability. Therefore we suggest, with some hesitation, that the B86SRP68-DF2 DF is the best DF at the moment for the H₂ + Ag(111) system. More and better defined experiments will allow us to refine our description of this system.

4.7.3 H₂ (D₂) + Au(111). With respect to the H₂ + Au(111) system we cannot with certainty asses the quality of any of the DFs tested here. Based on the good reproduction of the experimental trend in E_1/2(ν, J) parameters and the reasonable agreement between theory and experiment with respect to the state distributions of desorbing molecules, we can infer that the reactivity of the rovibrational states relative to each other is described reasonably well. We cannot say anything about the accuracy of the barrier without additional experiments or improvements in our dynamical model that will allow us to disentangle the effects of surface temperature, surface reconstruction and ehp excitations (see Section 4.10).

Shuai et al.⁷⁰ suggested that the PBE DF is better then the SRP48 DF because calculated time-of-flight distributions correlated slightly less worse with experimental observations for the PBE DF. This assertion implicitly assumes the validity of detailed balance for this system. The main conclusion of the experimentalists was however that the bad agreement between theory and experiment points toward strong non-adiabatic effects, which, if true, would invalidate the assumption of detailed balance.⁷⁰

4.7.4 D₂ + Pt(111). Ghassemi et al.¹¹ have previously designed a SRP-DF for the D₂ + Pt(111) system. Although the B86SRP68-DF2 DF could describe the experiments of Luntz et al.⁶⁴ to within chemical accuracy, the description of the experiments of Cao et al.⁶⁰ was not as good. As was the case for the H₂ + Ag(111) system there is some discussion about molecular beam parameters describing different experiments, but the experiments are in reasonably good agreement with each other.⁶⁵ Overall this system is best described by the PBEα57-DF2 DF that was specifically designed for this system.¹¹

4.8 Transferability

So far SRP-DFs fitted to reproduce molecular beam adsorption experiments for H₂ and D₂ where shown to be transferable among crystal faces of the same metal,^17,18 but not among different metals.^19,20 Here we show examples in which a SRP-DF that was fitted to reproduce molecular beam experiments of the activated late barrier system of H₂ reacting on Cu(111) can also describe the non-activated early barrier system of D₂ reacting on Pt(111) with chemical accuracy, and vice versa. Transferability to a different substrate of a SRP-DF fitted to reproduce gas-surface experiments has only been reported for CH₄ dissociation on Ni(111)¹³ to CH₄ dissociation on Pt(111).¹⁶

The SRP48 DF for H₂ + Cu(111) is not transferable to the H₂ + Ag(111) system¹⁹ or to the H₂ + Pt(111) system. We have previously shown that a SRP-DF based on the mGGA that does not include non-local correlation (MS-PBEl¹²) greatly improves the transferability from H₂ + Cu(111) to H₂ + Ag(111),¹² but Fig. 8 shows that this DF is not transferable to the weakly activated H₂ + Pt(111) system.

The only group of DFs that might be considered transferable between both highly activated and weakly activated systems are DFs that include non-local correlation. We demonstrated that a SRP-DF fitted to H₂ + Pt(111), PBEα57-DF2,¹¹ can describe H₂ + Cu(111) with overall chemical accuracy and that a new SRP-DF fitted to H₂ + Cu(111), B86SRP68-DF2, can describe the D₂ + Pt(111) experiments of Luntz et al.⁶⁴ with chemical accuracy. We speculate that the transferability between the Cu(111) and Pt(111) systems might be improved by taking into account relativistic corrections in our DFT calculations beyond those already included in the PAW potentials,¹³² which at this accuracy level might be important.^169,170

Both the B86SRP68-DF2 and PBEα57-DF2 DFs more or less predict the same reactivity for the H₂ + Ag(111) system. It is not possible however to call the DFs transferable to this system, yet. The lack of well described dissociative chemisorption experiments for this system does not yet allow us to make a broad statement about the accuracy of the theoretical description of this system. At present our description appears to be just shy of chemical accuracy.

4.9 Adiabatic description of S₀ and E_1/2(ν, J), a possible fingerprint for ehp excitations

The assumption of detailed balance entails that associative desorption is the inverse of dissociative chemisorption. In an adiabatic picture there is just one reason for a possible divergence of the obtained reaction probabilities. This is based on surface temperature, which is usually much higher in the associative desorption experiments than in the sticking experiments. This might lead to a breakdown of the detailed balance assumption that is usually involved when modeling associative desorption experiments with calculations on dissociative chemisorption. We note that it might be possible to model the associative desorption experiment directly,^75,171–173 thereby negating the need to invoke the principle of detailed balance and investigate if associative desorption is indeed the inverse of dissociative chemisorption.

Including ehp excitations in dissociative chemisorption calculations would lower the reactivity thereby increasing the effective barrier.²² Including ehp excitations in hypothetical associative desorption calculations, where molecules would start at the transition state and then desorb, would probably shift the translational energy distributions of desorbing molecules to lower energies and lead to lower effective barriers. When accounting for the effect of surface temperature, the difference in predicted reactivity as embodied by the E_1/2(ν, J) parameters obtained by including ehp excitations in associative desorption and dissociative chemisorption calculations, and their differences with results from adiabatic calculations, might then be taken to be a fingerprints for the effect of ehp excitation.

For H₂ + Cu(111) we know from DFMD calculations⁶² and other approaches^{23,29,155,156} that the effect of surface motion on the reactivity in dissociative chemisorption is small, even for high surface temperatures. Fig. 18 shows some evidence that the broadening of reaction probability curves might affect the calculated E_1/2(ν, J) parameters for low J. However, in an adiabatic picture, assuming detailed balance there should be no difference between calculations on dissociative adsorption and associative desorption.

We believe that this suggests a reason for only the SRP48 DF being chemically accurate for H₂ + Cu(111) for both dissociative chemisorption and associative desorption, since it overestimates the former and underestimates the latter predicted reactivity. The new SRP-DFs that are very accurate for dissociative chemisorption on cold surfaces, for which the BOSS model is valid,^{23,27–29,62,147} must underestimate the reactivity obtained from associative desorption by at least the extent to which surface temperature would increase the reactivity in dissociative chemisorption. Any remaining discrepancy can be safely attributed to the effect of ehp excitations.

The analysis of the H₂ + Au(111) system is more complicated due to the lack of dissociative chemisorption experiments and the current inability to take into account surface reconstruction (and thereby surface motion). Without at least either a dissociative chemisorption experiment or calculations using a reconstructed surface using the BOSS model, it is not yet possible to disentangle the contributions of the surface temperature, surface reconstruction and ehp excitations to the reactivity. Additionally, more detailed associative desorption and dissociative chemisorption experiments for the H₂ + Ag(111) system would allow for a systematic investigation into the effect of ehp excitations on reactivity for highly activated late barrier reactions.

5 Conclusions

We have constructed new SRP-DFs that include non-local correlation for the H₂ (D₂) + Cu(111) system and assessed the transferability of these DFs to the H₂ (D₂) + Ag(111), H₂ (D₂) + Au(111) and H₂ (D₂) + Pt(111) systems. All newly tested and developed DFs are based on GGA-exchange and use non-local correlation to describe dissociative chemisorption of H₂ (D₂) on Cu(111) within chemical accuracy, and, to the extent that it can be assessed, improve the transferability to the other systems discussed in this work over the previously reported SRP48 and MS-B86bl SRP-DFs.

The new SRP-DFs improve the description of the metal over the previously available SRP-DFs based on mixing GGA exchange, especially concerning calculated lattice constants. In general, the new SRP-DFs with non-local correlation exhibit higher and later barriers to reaction in combination with a slightly lower energetic corrugation. We also find that vdW-DF2 non-local correlation performs better than vdW-DF1 correlation for all tried combinations with different exchange functionals, except when the exchange part of a functional was specifically optimized for use with vdW-DF1 correlation. The B86SRP68-DF2 functional best reproduces the measured van der Waals well depths for H₂ + Cu(111), H₂ + Ag(111) and H₂ + Au(111).

SRP-DFs that include non-local correlation, namely B86SRP68-DF2 and PBEα57-DF2, are transferable from the highly activated late barrier H₂ + Cu(111) system to the weakly activated earlier barrier H₂ + Pt(111) system and vice versa. This feat could not be demonstrated with GGA and mGGA SRP-DFs that do not include non-local correlation. Assessing the transferability of the tested and developed SRP-DFs to H₂ + Ag(111) and H₂ + Au(111) is difficult due to the lack of well characterized molecular beam experiments. The SRP-DFs that include non-local correlation predict similar results for molecular beam sticking of D₂ + Ag(111), which are just shy of chemical accuracy. However it should be noted that there is some discussion about the validity of the beam parameters describing this particular molecular beam experiment.

A detailed analysis of associative desorption experiments on Cu(111) suggest that accurate calculation of E_1/2(ν, J) parameters requires an improvement of our dynamical model. Describing the surface degrees of freedom might close the gap between the excellent description of dissociative chemisorption and the good description of associative desorption, for molecules in the vibrational ground state. Any discrepancy in predicted reactivity between simulated associative desorption and dissociative chemisorption remaining after taking into account the effect of surface atom motion can then most likely be attributed to ehp excitation. Lack of additional experiments for the H₂ + Au(111) system, specifically a well described dissociative chemisorption experiment, presently keeps us from disentangling the effects of surface reconstruction, surface temperature and ehp excitation for this system.

We have carried out a full molecular beam simulation for the H₂ + Cu(111) system using the QD method and the B86SRP68-DF2 DF for sticking in this system, which is the best performing DF for this system, and which includes non-local correlation. Overall H₂ + Cu(111) is very well described quasi-classically when looking at molecular beam reaction probabilities or degeneracy averaged reaction probabilities. At the level of molecular beam sticking, and degeneracy averaged reaction probabilities, the differences between the QD and QCT method are very small. The QD method predicts slightly higher reaction probabilities for molecular beam sticking for very narrow low average translational energy molecular beams when comparing to QCT results based on the same set of initial rovibrational states. When looking at initial-state resolved reaction probabilities the QD method predicts a somewhat larger orientational dependence of the reaction, in favor of molecules reacting in a parallel orientation. With respect to vibrationally and rotationally inelastic scattering of H₂ from Cu(111) the B86SRP68-DF2 DF performs almost as well as the previous best SRP-DF for this system, namely the SRP48 DF.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We would like to thank Dr M. F. Somers for the fruitful discussions about the different dynamical models employed in this work, as well as Dr E. N. Ghassemi. Without her previous work on most of the systems treated here this paper would not have been possible. Lastly we would like to thank our student B. Hagedoorn BSc, whose bachelor thesis constituted the start of the H₂ + Au(111) project. We would also like to thank Dr S. Kaufmann for useful discussions and for providing us with the E_max(ν, J) parameters from the experiments of H₂ (D₂) + Au(111). This work was supported financially through a NWO/CW TOP grant (No. 715.017.001), and by a grant of supercomputer time from NWO Exacte en Natuurwetenschappen (NWO-ENW, grant number 2019.015).

References

1. G.-J. Kroes, J. Phys. Chem. Lett., 2015, 6, 4106–4114.
2. C. A. Wolcott, A. J. Medford, F. Studt and C. T. Campbell, J. Catal., 2015, 330, 197–207.
3. M. K. Sabbe, M.-F. Reyniers and K. Reuter, Catal. Sci. Technol., 2012, 2, 2010–2024.
4. G. Ertl, Angew. Chem., Int. Ed., 2008, 47, 3524–3535.
5. G.-J. Kroes, Science, 2008, 321, 794–797.
6. K. Waugh, Catal. Today, 1992, 15, 51–75.
7. L. Grabow and M. Mavrikakis, ACS Catal., 2011, 1, 365–384.
8. M. Behrens, F. Studt, I. Kasatkin, S. Kühl, M. Hävecker, F. Abild-Pedersen, S. Zander, F. Girgsdies, P. Kurr and B.-L. Kniep, et al. , Science, 2012, 336, 893–897.
9. G. B. Park, T. N. Kitsopoulos, D. Borodin, K. Golibrzuch, J. Neugebohren, D. J. Auerbach, C. T. Campbell and A. M. Wodtke, Nat. Rev. Chem., 2019, 1–10.
10. C. Daz, E. Pijper, R. Olsen, H. Busnengo, D. Auerbach and G. Kroes, Science, 2009, 326, 832–834.
11. E. N. Ghassemi, M. Wijzenbroek, M. F. Somers and G.-J. Kroes, Chem. Phys. Lett., 2017, 683, 329–335.
12. E. W. F. Smeets, J. Voss and G.-J. Kroes, J. Phys. Chem. A, 2019, 123, 5395–5406.
13. F. Nattino, D. Migliorini, G.-J. Kroes, E. Dombrowski, E. A. High, D. R. Killelea and A. L. Utz, J. Phys. Chem. Lett., 2016, 7, 2402–2406.
14. K. Doblhoff-Dier, J. Meyer, P. E. Hoggan and G.-J. Kroes, J. Chem. Theory Comput., 2017, 13, 3208–3219.
15. J. Boereboom, M. Wijzenbroek, M. Somers and G. Kroes, J. Chem. Phys., 2013, 139, 244707.
16. D. Migliorini, H. Chadwick, F. Nattino, A. Gutiérrez-González, E. Dombrowski, E. A. High, H. Guo, A. L. Utz, B. Jackson and R. D. Beck, et al. , J. Phys. Chem. Lett., 2017, 8, 4177–4182.
17. E. N. Ghassemi, E. W. F. Smeets, M. F. Somers, G.-J. Kroes, I. M. Groot, L. B. Juurlink and G. Füchsel, J. Phys. Chem. C, 2019, 123, 2973–2986.
18. L. Sementa, M. Wijzenbroek, B. Van Kolck, M. Somers, A. Al-Halabi, H. F. Busnengo, R. Olsen, G.-J. Kroes, M. Rutkowski and C. Thewes, et al. , J. Chem. Phys., 2013, 138, 044708.
19. E. Nour Ghassemi, M. Somers and G.-J. Kroes, J. Phys. Chem. C, 2018, 122, 22939–22952.
20. T. Tchakoua, E. W. F. Smeets, M. Somers and G.-J. Kroes, J. Phys. Chem. C, 2019, 123, 20420–20433.
21. Ž. Šljivančanin and B. Hammer, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 65, 085414.
22. P. Spiering and J. Meyer, J. Phys. Chem. Lett., 2018, 9, 1803–1808.
23. P. Spiering, M. Wijzenbroek and M. Somers, J. Chem. Phys., 2018, 149, 234702.
24. G.-J. Kroes, J. Juaristi and M. Alducin, J. Phys. Chem. C, 2017, 121, 13617–13633.
25. M. Somers, D. McCormack, G.-J. Kroes, R. Olsen, E. Baerends and R. Mowrey, J. Chem. Phys., 2002, 117, 6673–6687.
26. F. Nattino, C. Díaz, B. Jackson and G.-J. Kroes, Phys. Rev. Lett., 2012, 108, 236104.
27. G.-J. Kroes and C. Daz, Chem. Soc. Rev., 2016, 45, 3658–3700.
28. C. Díaz, R. A. Olsen, D. J. Auerbach and G.-J. Kroes, Phys. Chem. Chem. Phys., 2010, 12, 6499–6519.
29. M. Wijzenbroek and M. F. Somers, J. Chem. Phys., 2012, 137, 054703.
30. D. A. McCormack, G.-J. Kroes, R. A. Olsen, J. A. Groeneveld, J. N. van Stralen, E. J. Baerends and R. C. Mowrey, Faraday Discuss., 2000, 117, 109–132.
31. D. A. McCormack, G.-J. Kroes, R. A. Olsen, J. A. Groeneveld, J. N. van Stralen, E. J. Baerends and R. C. Mowrey, Chem. Phys. Lett., 2000, 328, 317–324.
32. J. Chen, X. Zhou and B. Jiang, J. Chem. Phys., 2019, 150, 061101.
33. G. Darling and S. Holloway, J. Chem. Phys., 1994, 101, 3268–3281.
34. S. Sakong and A. Groß, Surf. Sci., 2003, 525, 107–118.
35. G.-J. Kroes, E. Pijper and A. Salin, J. Chem. Phys., 2007, 127, 164722.
36. E. W. F. Smeets, G. Füchsel and G.-J. Kroes, J. Phys. Chem. C, 2019, 123, 23049–23063.
37. J. Dai and J. C. Light, J. Chem. Phys., 1998, 108, 7816–7820.
38. J. Dai and J. C. Light, J. Chem. Phys., 1997, 107, 1676–1679.
39. J. Dai and J. Z. Zhang, J. Chem. Phys., 1995, 102, 6280–6289.
40. B. Hammer, M. Scheffler, K. W. Jacobsen and J. K. Nørskov, Phys. Rev. Lett., 1994, 73, 1400.
41. M. Wijzenbroek, D. M. Klein, B. Smits, M. F. Somers and G.-J. Kroes, J. Phys. Chem. A, 2015, 119, 12146–12158.
42. G. Füchsel, K. Cao, S. Er, E. W. F. Smeets, A. W. Kleyn, L. B. F. Juurlink and G.-J. Kroes, J. Phys. Chem. Lett., 2018, 9, 170–175.
43. G.-J. Kroes, M. Pavanello, M. Blanco-Rey, M. Alducin and D. J. Auerbach, J. Chem. Phys., 2014, 141, 054705.
44. H. A. Michelsen and D. J. Auerbach, J. Chem. Phys., 1991, 94, 7502–7520.
45. H. Michelsen, C. Rettner, D. Auerbach and R. Zare, J. Chem. Phys., 1993, 98, 8294–8307.
46. H. Berger, M. Leisch, A. Winkler and K. Rendulic, Chem. Phys. Lett., 1990, 175, 425–428.
47. S. Kaufmann, Q. Shuai, D. J. Auerbach, D. Schwarzer and A. M. Wodtke, J. Chem. Phys., 2018, 148, 194703.
48. G. Anger, A. Winkler and K. Rendulic, Surf. Sci., 1989, 220, 1–17.
49. H. Hou, S. Gulding, C. Rettner, A. Wodtke and D. Auerbach, Science, 1997, 277, 80–82.
50. G. Comsa and R. David, Surf. Sci., 1982, 117, 77–84.
51. C. Rettner, H. Michelsen and D. Auerbach, J. Chem. Phys., 1995, 102, 4625–4641.
52. C. Rettner, D. Auerbach and H. Michelsen, J. Vac. Sci. Technol., A, 1992, 10, 2282–2286.
53. C. Rettner, H. Michelsen and D. Auerbach, Chem. Phys., 1993, 175, 157–169.
54. A. Hodgson, P. Samson, A. Wight and C. Cottrell, Phys. Rev. Lett., 1997, 78, 963.
55. H. Michelsen, C. Rettner and D. Auerbach, Surf. Sci., 1992, 272, 65–72.
56. T. Okazawa, F. Takeuchi and Y. Kido, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 72, 075408.
57. S. Lindgren, L. Walldén, J. Rundgren and P. Westrin, Phys. Rev. B: Condens. Matter Mater. Phys., 1984, 29, 576.
58. K. Chae, H. Lu and T. Gustafsson, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 14082.
59. S. Andersson and M. Persson, Phys. Rev. Lett., 1993, 70, 202.
60. K. Cao, G. Füchsel, A. W. Kleyn and L. B. Juurlink, Phys. Chem. Chem. Phys., 2018, 20, 22477–22488.
61. U. Harten, J. P. Toennies and C. Wöll, J. Chem. Phys., 1986, 85, 2249–2258.
62. F. Nattino, A. Genova, M. Guijt, A. S. Muzas, C. Daz, D. J. Auerbach and G.-J. Kroes, J. Chem. Phys., 2014, 141, 124705.
63. K. Cao, R. van Lent, A. Kleyn and L. Juurlink, Chem. Phys. Lett., 2018, 706, 680–683.
64. A. Luntz, J. Brown and M. Williams, J. Chem. Phys., 1990, 93, 5240–5246.
65. E. N. Ghassemi, M. F. Somers and G.-J. Kroes, J. Phys. Chem. C, 2019, 123, 10406–10418.
66. E. Pijper, G.-J. Kroes, R. A. Olsen and E. J. Baerends, J. Chem. Phys., 2002, 117, 5885–5898.
67. G.-J. Kroes and M. F. Somers, J. Theor. Comput. Chem., 2005, 04, 493–581.
68. L. M. Raff and M. Karplus, J. Chem. Phys., 1966, 44, 1212–1229.
69. M. Cardillo, M. Balooch and R. Stickney, Surf. Sci., 1975, 50, 263–278.
70. Q. Shuai, S. Kaufmann, D. J. Auerbach, D. Schwarzer and A. M. Wodtke, J. Phys. Chem. Lett., 2017, 8, 1657–1663.
71. A. Perrier, L. Bonnet, D. Liotard and J.-C. Rayez, Surf. Sci., 2005, 581, 189–198.
72. A. Perrier, L. Bonnet and J.-C. Rayez, J. Phys. Chem. A, 2006, 110, 1608–1617.
73. A. Perrier, L. Bonnet and J.-C. Rayez, J. Chem. Phys., 2006, 124, 194701.
74. C. Daz, A. Perrier and G. Kroes, Chem. Phys. Lett., 2007, 434, 231–236.
75. O. Galparsoro, S. Kaufmann, D. J. Auerbach, A. Kandratsenka and A. M. Wodtke, Phys. Chem. Chem. Phys., 2020, 22, 17532–17539.
76. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh and C. Fiolhais, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 6671.
77. M. Murphy and A. Hodgson, Surf. Sci., 1997, 390, 29–34.
78. M. Murphy and A. Hodgson, Phys. Rev. Lett., 1997, 78, 4458.
79. B. Jiang and H. Guo, Phys. Chem. Chem. Phys., 2014, 16, 24704–24715.
80. B. Jiang and H. Guo, J. Chem. Phys., 2014, 141, 034109.
81. G. Lee and E. Plummer, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 51, 7250.
82. P. Sprunger and E. Plummer, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 48, 14436.
83. P. Avouris, D. Schmeisser and J. Demuth, Phys. Rev. Lett., 1982, 48, 199.
84. F. Healey, R. Carter, G. Worthy and A. Hodgson, Chem. Phys. Lett., 1995, 243, 133–139.
85. D. H. Parker, M. E. Jones and B. E. Koel, Surf. Sci., 1990, 233, 65–73.
86. G. Lee, P. Sprunger, M. Okada, D. Poker, D. Zehner and E. Plummer, J. Vac. Sci. Technol., A, 1994, 12, 2119–2123.
87. H. Asada, Surf. Sci., 1979, 81, 386–408.
88. J. M. Horne, S. C. Yerkes and D. R. Miller, Surf. Sci., 1980, 93, 47–63.
89. C.-f. Yu, K. B. Whaley, C. S. Hogg and S. J. Sibener, Phys. Rev. Lett., 1983, 51, 2210.
90. C. Yu, K. B. Whaley, C. S. Hogg and S. J. Sibener, J. Chem. Phys., 1985, 83, 4217–4234.
91. L. Mattera, R. Musenich, C. Salvo and S. Terreni, Faraday Discuss., 1985, 80, 115–126.
92. R. J. Maurer, Y. Zhang, H. Guo and B. Jiang, Faraday Discuss., 2019, 214, 105–121.
93. Y. Zhang, R. J. Maurer, H. Guo and B. Jiang, Chem. Sci., 2019, 10, 1089–1097.
94. R. J. Maurer, B. Jiang, H. Guo and J. C. Tully, Phys. Rev. Lett., 2017, 118, 256001.
95. A. S. Muzas, J. I. Juaristi, M. Alducin, R. D. Muiño, G.-J. Kroes and C. Daz, J. Chem. Phys., 2012, 137, 064707.
96. M. Dion, H. Rydberg, E. Schröder, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett., 2004, 92, 246401.
97. K. Lee, É. D. Murray, L. Kong, B. I. Lundqvist and D. C. Langreth, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 081101.
98. J. P. Cowin, C.-F. Yu, S. J. Sibener and J. E. Hurst, J. Chem. Phys., 1981, 75, 1033–1034.
99. J. P. Cowin, C.-F. Yu, S. J. Sibener and L. Wharton, J. Chem. Phys., 1983, 79, 3537–3549.
100. J. P. Cowin, C.-F. Yu and L. Wharton, Surf. Sci., 1985, 161, 221–233.
101. C.-F. Yu, C. S. Hogg, J. P. Cowin, K. B. Whaley, J. C. Light and S. J. Sibener, Isr. J. Chem., 1982, 22, 305–314.
102. J. Perreau and J. Lapujoulade, Surf. Sci., 1982, 122, 341–354.
103. A. Kaufhold and J. P. Toennies, Surf. Sci., 1986, 173, 320–336.
104. S. Andersson, L. Wilzén and M. Persson, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 38, 2967.
105. J. Lennard-Jones and A. Devonshire, Nature, 1936, 137, 1069–1070.
106. H. Hoinkes and H. Wilsch, Helium Atom Scattering from Surfaces, Springer, 1992, pp. 113–172.
107. M. F. Bertino, F. Hofmann and J. P. Toennies, J. Chem. Phys., 1997, 106, 4327–4338.
108. P. Nieto, E. Pijper, D. Barredo, G. Laurent, R. A. Olsen, E.-J. Baerends, G.-J. Kroes and D. Faras, Science, 2006, 312, 86–89.
109. M. Bertino, S. Miret-Artés, J. Toennies and G. Benedek, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 56, 9964.
110. B. Poelsema, K. Lenz and G. Comsa, J. Phys.: Condens. Matter, 2010, 22, 304006.
111. I. Hamada, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 121103.
112. B. Hammer, L. B. Hansen and J. K. Nørskov, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 7413.
113. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, L. A. Constantin and J. Sun, Phys. Rev. Lett., 2009, 103, 026403.
114. J. Klimeš, D. R. Bowler and A. Michaelides, J. Phys.: Condens. Matter, 2009, 22, 022201.
115. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
116. G. K. Madsen, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 75, 195108.
117. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
118. Y. Zhang and W. Yang, Phys. Rev. Lett., 1998, 80, 890.
119. J. P. Perdew and W. Yue, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 33, 8800.
120. H. Busnengo, A. Salin and W. Dong, J. Chem. Phys., 2000, 112, 7641–7651.
121. G. Füchsel, M. del Cueto, C. Díaz and G.-J. Kroes, J. Phys. Chem. C, 2016, 120, 25760–25779.
122. J. Stoer and R. Bulirsch, Introduction to numerical analysis, Springer Science & Business Media, 2013, vol. 12.
123. R. Kosloff, Time-Dependent Quantum Molecular Dynamics, Springer, 1992, pp. 97–116.
124. D. Kosloff and R. Kosloff, J. Comput. Phys., 1983, 52, 35–53.
125. G. G. Balint-Kurti, R. N. Dixon and C. C. Marston, Int. Rev. Phys. Chem., 1992, 11, 317–344.
126. R. N. Zare and W. G. Harter, Phys. Today, 1989, 42, 68.
127. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
128. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
129. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
130. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
131. S. Lehtola, C. Steigemann, M. J. Oliveira and M. A. Marques, SoftwareX, 2018, 7, 1–5.
132. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
133. J. c. v. Klimeš, D. R. Bowler and A. Michaelides, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 195131.
134. G. Román-Pérez and J. M. Soler, Phys. Rev. Lett., 2009, 103, 096102.
135. D. Vanderbilt, Phys. Rev. B: Condens. Matter Mater. Phys., 1990, 41, 7892.
136. P. Haas, F. Tran and P. Blaha, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 085104.
137. P. Statiris, H. Lu and T. Gustafsson, Phys. Rev. Lett., 1994, 72, 3574.
138. E. Soares, G. Leatherman, R. Diehl and M. Van Hove, Surf. Sci., 2000, 468, 129–136.
139. R. Nichols, T. Nouar, C. Lucas, W. Haiss and W. Hofer, Surf. Sci., 2002, 513, 263–271.
140. H. Ohtani, M. Van Hove and G. Somorjai, Surf. Sci., 1987, 187, 372–386.
141. D. L. Adams, H. Nielsen and M. A. Van Hove, Phys. Rev. B: Condens. Matter Mater. Phys., 1979, 20, 4789.
142. M. Wijzenbroek, D. Helstone, J. Meyer and G.-J. Kroes, J. Chem. Phys., 2016, 145, 144701.
143. A. Patra, J. E. Bates, J. Sun and J. P. Perdew, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, E9188–E9196.
144. C. Cottrell, R. Carter, A. Nesbitt, P. Samson and A. Hodgson, J. Chem. Phys., 1997, 106, 4714–4722.
145. I. Groot, H. Ueta, M. Van der Niet, A. Kleyn and L. Juurlink, J. Chem. Phys., 2007, 127, 244701.
146. S. Kaufmann, personal communication.
147. C. Y. Wei, S. P. Lewis, E. J. Mele and A. M. Rappe, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 10062–10068.
148. K. Shakouri, J. Behler, J. Meyer and G.-J. Kroes, J. Phys. Chem. Lett., 2017, 8, 2131–2136.
149. N. Gerrits, H. Chadwick and G.-J. Kroes, J. Phys. Chem. C, 2019, 123, 24013–24023.
150. N. Gerrits, K. Shakouri, J. Behler and G.-J. Kroes, J. Phys. Chem. Lett., 2019, 10, 1763–1768.
151. L. Schimka, J. Harl, A. Stroppa, A. Grüneis, M. Marsman, F. Mittendorfer and G. Kresse, Nat. Mater., 2010, 9, 741–744.
152. K. Lee, K. Berland, M. Yoon, S. Andersson, E. Schröder, P. Hyldgaard and B. I. Lundqvist, J. Phys.: Condens. Matter, 2012, 24, 424213.
153. G. Darling and S. Holloway, J. Chem. Phys., 1990, 93, 9145–9156.
154. J. Juaristi, M. Alducin, R. D. Muiño, H. F. Busnengo and A. Salin, Phys. Rev. Lett., 2008, 100, 116102.
155. A. Mondal, M. Wijzenbroek, M. Bonfanti, C. Daz and G.-J. Kroes, J. Phys. Chem. A, 2013, 117, 8770–8781.
156. M. Bonfanti, M. F. Somers, C. Daz, H. F. Busnengo and G.-J. Kroes, Z. Phys. Chem., 2013, 227, 1397–1420.
157. S. Mukherjee, F. Libisch, N. Large, O. Neumann, L. V. Brown, J. Cheng, J. B. Lassiter, E. A. Carter, P. Nordlander and N. J. Halas, Nano Lett., 2012, 13, 240–247.
158. E. Hasselbrink, Surf. Sci., 2009, 603, 1564–1570.
159. B. Schindler, D. Diesing and E. Hasselbrink, J. Chem. Phys., 2011, 134, 034705.
160. B. Schindler, D. Diesing and E. Hasselbrink, J. Phys. Chem. C, 2013, 117, 6337–6345.
161. J. C. Polanyi, Science, 1987, 236, 680–690.
162. J. C. Polanyi, Angew. Chem., Int. Ed. Engl., 1987, 26, 952–971.
163. G.-J. Kroes, C. Daz, E. Pijper, R. A. Olsen and D. J. Auerbach, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 20881–20886.
164. J. Sun, A. Ruzsinszky and J. P. Perdew, Phys. Rev. Lett., 2015, 115, 036402.
165. M. Gostein, H. Parhikhteh and G. Sitz, Phys. Rev. Lett., 1995, 75, 342.
166. E. Watts and G. O. Sitz, J. Chem. Phys., 2001, 114, 4171–4179.
167. G.-J. Kroes, G. Wiesenekker, E. Baerends and R. Mowrey, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 53, 10397.
168. A. Salin, J. Chem. Phys., 2006, 124, 104704.
169. R. Olsen, G. Kroes and E. Baerends, J. Chem. Phys., 1999, 111, 11155–11163.
170. R. Olsen, P. Philipsen and E. Baerends, J. Chem. Phys., 2003, 119, 4522–4528.
171. G. Mills, G. Schenter, D. E. Makarov and H. Jónsson, Chem. Phys. Lett., 1997, 278, 91–96.
172. G. Mills, H. Jónsson and G. K. Schenter, Surf. Sci., 1995, 324, 305–337.
173. G. Mills and H. Jónsson, Phys. Rev. Lett., 1994, 72, 1124.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0cp05173j

---