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

Specific reaction parameter density functionals (SRP-DFs) that can describe dissociative chemisorption molecular beam experiments of hydrogen (H 2 ) on cold transition metal surfaces with chemical accuracy have so far been shown to be only transferable among diﬀerent facets of the same metal, but not among diﬀerent metals. We design new SRP-DFs that include non-local vdW-DF2 correlation for the H 2 + Cu(111) system, and evaluate their transferability to the highly activated H 2 + Ag(111) and H 2 + Au(111) systems and the non-activated H 2 + Pt(111) system. We design our functionals for the H 2 + 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 2 + Cu(111) with chemical accuracy can also describe such experiments for H 2 + 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 2 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 2 + 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 2 on Cu(111) is still well described quasi-classically.


Introduction
Hydrogen (H 2 ) 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][2][3][4][5] Industrially H 2 dissociation is an important step in the production of methanol from CO 2 over a Cu/ZnO/Al 2 O 3 catalyst. [6][7][8] Additionally H 2 dissociation is an important process in the production of syngas and ammonia. 4 Increasing the predictive power of theoretical models of heterogeneous catalysis potentially has an important financial impact on the chemical industry. 9 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][11][12][13] i.e. DFs that can describe reaction barrier heights to within 1 kcal mol À1 . 5 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 2 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 14 requires experimental data as reference data. 1,[10][11][12][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 4 dissociation on Ni (111), 13 which could also describe the dissociation of CH 4 on Pt (111) with chemical accuracy, 16 where Ni and Pt belong to the same group. So far, for SRP-DFs fitted to reproduce molecular beam adsorption experiments for H 2 interacting with transition metals transferability was shown among systems in which H 2 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 4 + Ni (111) system to the CH 4 + Pt (111) system suggests that nonlocal 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 2 + 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 2 + Cu(111) system, since theoretically 10,12, and experimentally [45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][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 2 dissociation on cold metals. 23,[27][28][29]62 Additionally we evaluate the performance of a SRP-DF that was fitted to the H 2 + Pt(111) system 11 for the H 2 + 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 2 + Cu(111) system we will consider such experiments for the H 2 + Pt(111) system 63,64 and the H 2 + Ag (111) system. 54 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 2 + Cu(111) system's large body of experimental work [45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][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 2 + 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) 68 and QD results compared to the good agreement that was obtained for these two methods for the H 2 + Cu(211) system 36 using the SRP48 DF, 10,62 which does not employ non-local correlation.
For the H 2 + Cu (111) system it is known that the effect of surface motion cannot readily be ignored for specific observables at high surface temperature 62 (T s ). Analysing associative desorption and dissociative chemisorption experiments as linked through detailed balance, 69 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. 70 for the H 2 + 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 0 (n, J)) parameters. 47 Even though a complementary molecular beam dissociative chemisorption experiment is not available for the associative desorption experiment of Shuai et al. 70 on H 2 + 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][72][73][74][75] and that this has also been done for H 2 and D 2 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 39 to unconverged DFT calculations 40 using the PW91 DF, 76 and the statistical accuracy of the results of the later work 75 was limited by the number of ab initio molecular dynamics (AIMD) trajectories that could be run. However, an interesting aspect of the later work 75 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 2 + 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. 24 For the H 2 + Ag(111) system the only molecular beam dissociative chemisorption experiments we are aware of are those of Hodgson and coworkers. 54 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 79 using quantum dynamics calculations on a permutation invariant polynomial (PIP) neural network potential 80 and in work done in our group. 12,19 Earlier experimental work on the H 2 + Ag (111) system suggested that H 2 prefers to physisorb on silver surfaces, [81][82][83] and that the dissociation of H 2 on silver is endothermic, 84 exhibiting a relatively low barrier 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 2 (D 2 ) + metal systems. Initial-state resolved reaction probabilities are presented in Section 3.3, and Section 3.4 presents parameters, E 1/2 (n, 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 2 relative to the surface (parallel or end-on) are presented in Section 3.5 for the D 2 + Cu (111) system. Rotationally and vibrationally inelastic scattering of H 2 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 (n, J) parameters with values extracted from associative desorption experiments on H 2 + Cu (111) and Au (111). We also compare initial-state resolved reaction probabilities computed with theory with experimental data for H 2 + 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 2 + Cu (111). The inelastic scattering of H 2 from Cu(111) is discussed in Section 4.5 The quality of the QCT results is evaluated by comparison with QD results for H 2 + Cu (111) in Section 4.6. Section 4.7 contains an overview of the comparison theory-experiment for all the H 2 -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.

Coordinate system
All calculations in this work are carried out using the BOSS model. 10 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 2 (D 2 ) 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 2 bond length r, the polar orientation angle with respect to the surface normal y and the azimuthal angle f. 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.

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 96 or vdW-DF2 non-local correlation 97 (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. 10 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 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 a is the SRP mixing parameter, r is the three dimensional electron density and rr is the gradient of the electron density. E 1 X (r,rr) and E 2 X (r,rr) 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)).

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). 120 The method we use is analogous to that used by Wijzenbroek et al., 41 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 120 PESs with respect to the underlying DFT calculations. Details are presented in the ESI. †

Quasi-classical dynamics
The QCT calculations presented in this work are carried out on six dimensional PESs and assume quasi-classical initial conditions. 68 This means that we take into account the quantum mechanical energies of the impinging H 2 and D 2 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. 122 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.

Quantum dynamics
Six dimensional quantum dynamics (QD) calculations are performed by solving the time-dependent Schrödinger equation, i h dCðQ; tÞ dt ¼ĤCðQ; tÞ; using the time-dependent wave packet (TDWP) method 123,124 with our in-house computer package. 66,67 Here, C( -Q,t) denotes the nuclear wave function of H 2 at time t with -Q being the Fig. 1 The COM coordinate system used for the description of the H 2 (D 2 ) 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 2 (D 2 ) 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 (y = 901,f = 01) 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. position vector. Furthermore we employ the following Hamiltonian in order to take into account the six degrees of freedom of H 2 : Here, M and m are the mass and reduced mass of H 2 ,Ĵ 2 (y,f) is the angular momentum operator and V( -Q) is the six dimensional PES. The scattered wave packet is analysed using the scattering matrix formalism, 125 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 P nJm J !n 0 J 0 m 0 J nm ðEÞ 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 0 is computed using S 0 ðhEiÞ ¼ X n;J ð P v; n; J; T n ð Þ P deg ðE; n; JÞdE: Here, hEi 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 n and the angular momentum quantum number J and to have a velocity between v and v + dv is denoted by: P(v,n, J,T n )dv = P flux (v;T n )dv Â P int (n, J,T n ).
The flux-weighted velocity distribution P flux is a function of T n and is determined by the width parameter a and the stream velocity v 0 according to 44 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 (n, J,T n ) = (2J + 1) Â e (À(En,0ÀE0,0)/kBTvib) Â e (À(En,JÀEn,0)/kBTrot) .
Here, k B is the Boltzmann constant and E n,J is the energy of the quantum state characterized by n 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 , 51 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 2 , g N is 1/4 (3/4) for even (odd) values of J, and for D 2 , 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, n, 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, n, J), as: where P r (E, n, 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, n, 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. 36 2.6.2 Comparing to experimental E 0 (n, J) parameters. Experimentally, for H 2 -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. 69 Experiments on H 2 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, n, J), may be related to the degeneracy averaged initial-state resolved reaction probability, using: P des ðE; n; JÞ / Ee À E k B T S P deg ðE; n; JÞ: The extracted reaction probabilities are usually fitted to a sigmoid function, e.g. the function involving the error function: Here, the A n,J values are the saturation values of the extracted degeneracy averaged reaction probabilities, and the effective barrier height (E 0 (n, J)) is the incidence energy at which P deg (E, n, J) first becomes equal to 1 2 A n;J . Using eqn (11), E 0 (n, J) also is the inflexion point of the reaction probability curve if the saturation value A n,J corresponds to the absolute saturation value. W n,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 n,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 2 + Cu(111) 47,51 and D 2 + 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 2 + 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 2 and D 2 + Au (111), in which A n,J was set to one for (n = 0,J = 6) H 2 and for (n = 0,J = 0) D 2 , and the A n,J values for different (n, J) states of H 2 and D 2 , respectively, reflected the values of the reaction probabilities relative to these reference states. 70 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 (n, 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,n, 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 (n, J). The E 1/2 (n, 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 2 , D 2 and HD desorbing from Cu surfaces, 47 and for H 2 and D 2 desorbing from Au (111). 70 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 (2) 0 (E,n, J), is computed from initial-state resolved reaction probabilities as follows: 126 The rotational quadrupole alignment parameter is a measure of the dependence of the reaction on the alignment of H 2 relative to the surface. 2.6.4 Rovibrational state populations of H 2 and D 2 desorbing from Au (111). State distributions of desorbing molecules are calculated in the following manner: 70 Here T S is the surface temperature, and E max (n, J) is the maximum kinetic energy sensitivity of the experiment, 70 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,n, 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(n,J) . The error function fits derived in ref. 70 are only valid below E max(n,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(n, J) populations are normalized to the total n = 0 population according to: The ratio n = 1 : n = 0 can then calculated as: The upper limits to J used in eqn (14) and (15) are discussed below.

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. 131 The standard VASP projector augmented wave (PAW) potentials 132 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 135 and PBE correlation 115 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 G-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 G-centered k-point grid and a vacuum distance of 16 Å.   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][57][58][137][138][139][140][141] 3.1.2 H 2 + metal surface PESs. Barrier heights and geometries for H 2 + Cu(111) for high symmetry geometries are shown in Table 3. The energetic corrugation x, 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.

Results
Barrier heights and positions for H 2 + Ag(111), Au (111) and Pt (111) are shown in Tables S5, S6 and S7 (ESI †) respectively. With respect to the H 2 + Pt(111) system the most striking result is that only the PBEa57-DF2 11 and MS-PBEl 12 DFs exhibit a double barrier structure for the top-to-bridge (t2b) geometrie whereas the other DFs tested do not. The PBEa57-DF2 11 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 120 PES for the B86SRP68-DF2 DF for H 2 + Cu(111) using B4900 randomly sampled geometries. Based on all the randomly sampled points taken together our CRP 120 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 2 with the surface lower then 4 eV the rms error reduces to 8 meV (B0.2 kcal mol À1 ). Our CRP 120 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 2 in a parallel orientation (f = 01,y = 901) above a top atom of Cu (111) for different DFs, comparing with experimental results. 59 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 96 and vdW-DF2 97 DFs. Agreement with experiment is best for the B86SRP68-DF2 DF.
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 2 + Ag (111), H 2 + Au (111), and H 2 + Pt(111) we find depths that are in good agreement with experimental work 61,90,99,110 for the B86SRP68-DF2 DF.

Molecular beam sticking probabilities
Molecular beam sticking probabilities computed with five DFs for H 2 and D 2 reacting on Cu (111) are shown in Fig. 4, comparing to experimental results of Auerbach and coworkers 45,51 and Rendulic and coworkers. 46 Fig. S5 (ESI †) shows the comparison with an additional experiment for the five DFs discussed here, namely for the pure D 2 molecular beams of Auerbach and coworkers. 45 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 2 + Cu(111) shown in Fig. 4 with chemical accuracy. Fig. 5 shows comparisons to two additional sets of For the B86SRP68-DF2 DF QD results are also shown for the experiments concerning H 2 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 60 have actually reported their molecular beam parameters. Fig. 8a shows that the B86SRP68-DF2 DF describes the experiment 64 with overall chemical accuracy, albeit that the energy shifts are larger than 4.2 kJ mol À1 at the lowest energies. However, Fig. S7a (ESI †) shows that the experiments of Cao et al. 60 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.
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 †).

Initial-state resolved reaction probabilities
Initial-state resolved reaction probabilities will be presented for H 2 reacting on both Cu (111) and Ag (111). Fig. 9 shows fully initial-state resolved reaction probabilities for H 2 reacting on Cu (111), as obtained with the QD and QCT methods. At the Table 3 Barrier heights for H 2 reacting on Cu(111) for high symmetry geometries. For the bridge site f = 901 and y = 901, for the t2b and hcp geometries f = 01 and y = 901. 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, x, is also shown in eV   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 o 3 rovibrational states. Fig. 10 shows degeneracy averaged initial-state resolved reaction probabilities for H 2 and D 2 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 2 and when using the MS-PBEl mGGA. 12 Note, however, that the P deg (E, n, J) extracted from experiments were not normalized, but simply assumed to saturate at unity, making it hard to perform an accurate comparison with experiment.
3.4 E 1/2 (n, J) parameters E 1/2 (n, J) parameters calculated for H 2 (D 2 ) + Cu(111) using Methods A1 and B1 are shown in Fig. 11 and 12, respectively, also comparing with experiment. 47 Using Method B1, the DFs that include non-local correlation qualitatively reproduce the rotational hindering observed experimentally for (n = 0), i.e. the increase of E 1/2 (n, 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 (n, J) parameters, which describes the dependence of the E 1/2 (n, 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 †). E 1/2 (n, J) parameters calculated for H 2 (D 2 ) + Au (111) using Method B2 are shown in Fig. 13, also comparing to experiment. 70 Accompanying MAD and mean signed deviations (MSD) values of the computed E 1/2 (n, J) parameters from the experimental values are presented in Table S10 (ESI †) for both H 2 (D 2 ) + Cu (111) and H 2 (D 2 ) + Au(111).

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 2 desorbing from Cu (111). 49 Note that a positive A (2) 0 (n, J) indicates a preference for parallel reaction, a negative value a preference for perpendicular reaction, and zero means the reaction proceeds independent of D 2 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. (111) Vibrationally inelastic scattering probabilities, P(n = 0,J-n = 1,J = 3) for H 2 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 (n = 0,J = 5) rovibrational state, except that the SRP48 DF yields smaller vibrational excitation probabilities for E 4 0.8 eV (panel c). For the (n = 0,J = 1) and (n = 0,J = 3) initial rovibrational states the differences are larger (panels a and b). Fig. 16 shows the ratio of rotationally elastic and inelastic scattering probabilities P(n = 1,J = 0n = 1,J = 2)/P(n = 1, J = 0n = 1,J = 0) computed with two DFs and comparing with    (111). Experiment is shown in black. 144 QCT results are shown for the following DFs: SRP48 (purple), 19 MS-PBEl (red), 12 B86SRP68-DF2 (blue), and PBEa57-DF2 (green). The dotted line represents the experimental curve shifted by 4.2 kJ mol À1 , denoting the limit to chemical accuracy. Experimental results are shown in red, 64 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. experiment. 54 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. 54 3.7 Rovibrational state populations of H 2 and D 2 desorbing from Au(111) Fig. 17 shows rovibrational state populations of H 2 and D 2 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 (n, J) state (see eqn (13)) and g N (2J + 1) being the statistical weight for rotational level J. 70 For D 2 , g N = 2 for even J and 1 for odd J; for H 2 , g N = 1 for even J and 3 for odd J. In such a plot a Boltzmann distribution will appear as a straight line. 70 Here we integrate eqn (13) up to E max (n, J).

Inelastic scattering of H 2 from Cu
Fig. S10 (ESI †) shows the rovibrational state populations of H 2 desorbing from Au(111) as reported in by Shuai et al. 70 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. 70 as outlined in a private communication. 146 Table 4 shows the ratio of the fluxes of molecules desorbing in the first excited and in the vibrational ground state for both H 2 and D 2 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 (n, J).    (111). Red circles represent the SRP48 values, 62 green circles the MS-B86bl values, 12 blue circles represent the B86SRP68-DF2 values with the solid blue circles corresponding to QD calculations, magenta circles represent the PBEa57-DF2 values, and solid black circles represent experimental results. 47 Our aim has been to develop new SRP-DFs that include nonlocal correlation for the H 2 + 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 2 and D 2 on Cu (111). [45][46][47][48][49][50][51][52][53][54][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 2 (D 2 ) + 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 2 (D 2 ) + 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 2 (D 2 ) + 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 2 dissociation on cold metals. 23,[27][28][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    (111). Experimental results are shown in black. 49 Theoretical results using the QCT method are shown for the B86SRP68-DF2 (red), SRP48 10 (green), PBEa57-DF2 11 (blue), and the MS-B86bl 12 (magenta) DFs.
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][149][150] Recently, this has been done for H 2 + Cu(111) 75 using ab initio molecular dynamics with electronic friction (AIM-DEF) calculations employing the PBE 115 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.

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. 151 Table S4 (ESI †) shows the calculated lattice constants for all DFs investigated in this work. Highly accurate results are only achieved using PBEsol 117 and our previously developed mGGA MS-B86bl. 12 All our candidate SRP-DFs yield improved results over the vdW-DF1, 96 vdW-DF2, 97      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. 12 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.

Static PES properties
The experimental van der Waals well depth that has been measured for H 2 + Cu (111) 59 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 96 and vdW-DF2, 97 yield wells that are too deep, especially for vdW-DF1, 96 as also found in earlier work 152 (Fig. 3b, see also Table S8, ESI †). The two previous SRP-DFs for this system produce a very tiny (SRP48 10 ) or no (MS-B86bl 12 ) van der Waals well at all (Fig. 3b). Of the four DFs considered in Fig. 3b, the vdW-DF2 97 DF produces the best results. Other exploratory calculations carried out by us showed that using vdW-DF1 96 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. 97 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 2 on Cu(111), as can be seen from the MAD values in Fig. 4. The PBEa57-DF2 DF has originally been fitted to reproduce experiments for H 2 (D 2 ) + Pt(111), 11 and the optPBE-DF1 114 DF has previously been shown to be chemically accurate H 2 (D 2 ) + Cu (111). 41 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. 114 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 97 (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 PBEa57-DF2 to SRPsol63-DF2 the calculated van der Waals well more closely resembles experiment. 59 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 2 + Cu(111), 59,61 H 2 + Ag(111), 90 H 2 + Au (111) 61 and H 2 + 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. 61 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, 104 and it has been remarked that the results obtained later 59 are in fact also consistent with the earlier measurements. 104 Additionally we suspect that the reported van der Waals wells for H 2 + Ag (111) 90 and H 2 + Au(111) 61 might possibly be too close to the surface 104 for the same reason. This drawback of using the RMSA technique 98 might be alleviated by redoing the potential inversion on the basis of the original data on Feshbach resonances with more advanced theoretical models, 153 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 98 Poelsema et al. 110 presented a combined thermal energy atom scattering/thermal desorption spectroscopy (TEAS/TDS) study of the H 2 + Pt(111) system, obtaining van der Waals well geometries that were subsequently accurately reproduced by theory. 11 For the H 2 + 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 2 + Cu(111), H 2 + Ag(111) and H 2 + Au (111). For the weakly activated H 2 + 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 (PBEa57-DF2) yields a value in the range of the experimental values.  Fig. 4 compares results of three sets of molecular beam sticking probabilities, S 0 , with results computed with the SRP48, 10 B86SRP68-DF2, SRPsol63-DF2, PBEa57-DF2 11 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 0 data point to the interpolated experimental data point with the same S 0 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 2 (D 2 ) + 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. 45 and Rettner et al. 53 considered here employed a low surface temperature of 120 K, the experiment of Berger et al. 46 was reportedly done on a 'cold' surface, and from the literature it is known that the BOSS model works well for activated H 2 dissociation on cold metals. 23,[26][27][28][29] Another advantage of fitting a SRP-DF to these sets of molecular beam experiments is that they cover both H 2 and D 2 for very different experimental conditions with respect to the nozzle temperature, the average collision energy, and the width of the velocity distributions. 10 The B86SRP68-DF2 DF exhibits a MAD of 1.4 kJ mol À1 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 2 + Cu(111) system. From this point onward we will mainly focus on the results obtained with the newly selected B86SRP68-DF2 DF and the PBEa57-DF2 11 DF, and on how the performance of these two DFs compares to that of the original SRP48 10 DF and of our previously developed mGGA DFs. 12 In the ESI † one additional comparison to experiment is made in Fig. S5, concerning pure D 2 molecular beams. 45 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 0 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. (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 2 + 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). 36 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 À1 to 1.6 kJ mol À1 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 51 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 2 + Cu(211) when using the SRP48 DF. 36 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 2 + Cu(211). 36 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 2 + Cu(111) is well described quasi-classically.

Molecular beam sticking of H 2 + Cu
Note that for a nozzle temperature of 2000 K we take into account all rovibrational states that have a Boltzmann weight 40.001. Highly excited rovibrational states, either with high n, 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 40.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. (111). Fig. 8 shows a comparison of calculations on D 2 + Pt (111) to the molecular beam experiments of Luntz et al. 64 The PBEa57-DF2 DF was originally fitted 11 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 À1 ), 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).

Molecular beam sticking of D 2 + Pt
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 PBEa57-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. 60 are not as well described by the B86SRP68-DF2 DF, as can be seen in  Table S9, ESI †). Earlier work from our group has shown that the experimental results of Luntz et al. 64 and Cao et al. 60 are in good agreement with each other for the lower incidence energies but somewhat diverge for the higher incidence energies. 65 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. 60 are most likely somewhat underestimated compared to the results of Luntz et al. 64 because the average incidence energies are somewhat underestimated in the experiment of Cao et al. 60 4.3.4 Molecular beam sticking of D 2 + Ag (111). We now make a comparison to molecular beam experiments on D 2 + Ag (111). 144 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. 19 Fig. 7 shows the computed S 0 for D 2 + Ag(111). Hodgson and coworkers 144 have reported translational energy distributions that were symmetric in the energy domain. In our view, and as discussed in previous work from our group, 19 the symmetric translational energy distributions are somewhat unphysical. Therefore we opted to use the molecular beam parameters of pure D 2 reacting on Cu (111) reported by Auerbach and coworkers, 45 which likewise describe beams that are narrow in translational energy. Here we see that the PBEa57-DF2 and B86SRP68-DF2 DF are similar in accuracy to our previously developed MS-PBEl mGGA DF, 12 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 PBEa57-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 2 + metal systems that incorporate GGA exchange. The MS-PBEl mGGA 12 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.

Associative desorption
4.4.1 Comparing to experimental E 0 (n, J) parameters. When comparing to experimental E 0 (n, J) parameters by calculating E 1/2 (n, 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 (n, 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 0 (n, 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 36 the MDEF method shifted the E 1/2 (n, 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 (n, J) parameters to slightly higher energies. However, as mentioned in our previous work, 36 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 (n, J) parameters from dissociative chemisorption calculations applying electronic friction would shift our E 1/2 (n, 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 36 only shifted the trend in E 1/2 (n, 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 2 from Cu(111) 75 using AIMD and AIMDEF calculations and the PBE DF, which find little effect of ehp excitations modeled at the local density friction approximation (LDFA) 154 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. (111). The comparison of the measured E 0 (n, J) parameters values with the E 1/2 (n, J) parameters computed using Method B1 is shown in Fig. 12. The E max (n, 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 0 (n, J) parameters with theoretical E 1/2 (n, J) parameters extracted using Method A1 is presented in Fig. 11. In this method we used A = 0.325 for H 2 and A = 0.513 for D 2 , as obtained in ref. 47. The use of Method A1 for H 2 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 2 . Nonetheless, SRP48 performs best for D 2 , 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 2 (0.325). This A value is much lower than the A value extracted for D 2 in Method A1 (0.513). This may well be a simple artifact resulting from the method followed: the A value determined for D 2 is likely to be more accurate because the sticking experiments were done for a kinetic energy up to 0.83 eV, 45 whereas the sticking experiments for H 2 only went up to about 0.5 eV. 51 This suggest that the A-value for D 2 (Fig. S1c, ESI †) is much more accurate than for H 2 (Fig. S1a, ESI †), and it is not clear why the A value for H 2 should differ much from it. Furthermore, the A value established for H 2 in Method A1 is much lower than the A-values established for H 2 in Method B1 (see also Fig. S1a and c, ESI †), while for D 2 the A values extracted with the two methods resemble each other, and the A values extracted with Method B1 for H 2 , much more (see Fig. S1b-d, ESI †).

E 1/2 (n, J) parameters Cu
Finally, we note at this stage the E 1/2 (n, J) parameters computed using Method A1 do, in general, not reproduce the subtle trend found experimentally that the E 0 (n, J) parameters first increase somewhat with J (see Fig. 11) (attributed to rotational hindering 45,47,51 ).
The E 1/2 (n, 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 (n, 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 nonlocal 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 2 , the agreement with the experimental dependence of the E 0 (n, 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 (n, 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 (n, J) parameters calculated using both Methods A1 and B1 do reproduce the clear trend ( Fig. 11 and 12) that at high J the E 0 (n, 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 39 used by Dai and Light 37 for sixdimensional QD calculations. As discussed in the ESI, † we have investigated whether the rotational hindering observed by Dai and Light 37 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 o 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 39 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 2 reacting on Cu (111). 23,62 Increased reactivity at low translational energies might lower calculated E 1/2 (n, J) parameters for both Methods A1 and B1. We have extracted E 1/2 (n, J) parameters using both Methods A1 and B1 from Nattino et al. 62 The results are shown in Fig. 18. The reported data was obtained using the SRP48 DF. 62 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. 47 To validate the obtained E 1/2 (n, J) parameters we have also extracted E 1/2 (n, J) parameters using the same methodology for the reported logistics function fits to BOSS direct dynamics data. 62 The good agreement for both Methods A1 and B1 between our SRP48 E 1/2 (n, 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 (n = 0,J = 11) suggests that a small decrease can be expected from allowing the surface to move, for the (n = 1,J = 6) point this is not so clear. There is however a dramatic decrease of the E 1/2 (n, J) parameter for (n = 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 (n, J) parameters, and this might thus be partly responsible for the observed rotational hindering. The decrease of the E 1/2 (n = 1,J = 4) value with the introduction of surface motion is less pronounced when using the A1 method.
Note however that the E 0 (n, J) parameters are extremely sensitive to the quality of the logistics function fits. The agreement between our SRP48 E 1/2 (n, 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 (n, 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. 75 4.4.1.2 E 1/2 (n, J) parameters Au (111). The comparison of the measured E 0 (n, J) values with the E 1/2 (n, J) parameters computed using Method B2 is shown in Fig. 13. The E max (n, J) parameters needed in this method (see eqn (S4), ESI †) have been obtained from a private communication, 146 and have been plotted in Fig. S2 (ESI †) against the E 0 (n, J) parameters extracted from the measurements. 70 We note that the E max (n, J) parameters typically equal E 0 ðn; JÞ þ 1 3 W for the n = 0 states and E 0 ðn; JÞ þ 2 3 W for the n = 1 states, where W = 0.31 eV for n = 0 and 0.29 eV for n = 1 70 (see Fig. S2, ESI †). This means that Method B1 will yield unreliable values for the E 1/2 (n, 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 (n, 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 2 + Au (111). The following tentative conclusions can be drawn: the E 1/2 (n, J) parameters computed with the mGGA DF tested here and with all DFs employing van der Waals correlation substantially overestimate the measured E 0 (n, J) parameters, with MAD values of approximately 0.1 eV for H 2 . The PBE DF would appear to perform best, with a MAD value of 46 meV for H 2 . This conclusions agrees with the observation that PBE reaction probabilities 142 allowed better fits of the measured time-of-flight spectra of H 2 and D 2 desorbing from Au(111) 70 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 (n, 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 2 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 (B1 kcal mol À1 ), 142 which would lead to slightly higher E 1/2 (n, J) parameters and therefore to increased disagreement with the measured values.
We are not able to resolve the contradiction posed by Shuai et al. 70 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 0 (n, 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 2 on Au(111) is coupled to the electronic degrees of freedom of the metal. [157][158][159][160] In line with Shuai et al. 70 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 2 (n = 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 2 + 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., 70 namely that the experimentally observed lower translational energy distributions compared to theoretical predictions (see Fig. 1 of Shuai et al. 70 ) is most likely due to ehp excitations in the desorption dynamics.
4.4.2 Rovibrational state populations of H 2 and D 2 desorbing from Au (111). Fig. 17 shows the rovibrational state populations of H 2 and D 2 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 (n = 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. 70 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 2 + Au (111). 161,162 Fig. S10 (ESI †) shows the state distributions of molecules desorbing from Au(111) as reported by Shuai et al. 70 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., 70 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. 70 Also when we shift the SRP48 (n = 0) and (n = 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 (n = 0) and (n = 1) rovibrational states is not affected by this discrepancy. The difference between our work and Shuai et al. 70 with respect to the n = 1 : n = 0 ratios arises because we use E max (n, J) as the upper integration limit in eqn (13). Shuai et al. 70 have used 5 eV as the upper integration limit. 146 The only reason we choose E max (n, J) as the upper integration boundary is because the reported error function fits are only reliable below E max (n, 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 n = 1: n = 0 ratios reported by Shuai et al. 70 In our view only integrating up to E max (n, J) is a more fair way of calculating the N(n, J) populations, though on the scale of Fig. 17 the difference between integrating to E max (n, 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 (n, J). Note also that we calculate the n = 1:n = 0 ratios only using the rovibrational states shown in Fig. 17, which is the same set of rovibrational states used by Shuai et al. 70 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 2 , although the qualitative agreement between theory and experiment is reasonable for both H 2 and D 2 . It can be seen in Table 4 that there is only a reasonable agreement for D 2 with respect to the n = 1 : n = 0 population ratios. The theoretical n = 1 : n = 0 population ratio for H 2 is however too low, a result similar to what was reported by Shuai et al. 70 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 n = 1 : n = 0 population ratio. 4.4.3 Initial-state resolved reaction probabilities Ag (111). In the case of molecular beam sticking results for D 2 + Ag(111) the MS-PBEl function performed similarly well as the GGAexchange 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 2 , 12 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 2 , 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. 12 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 2 seeded in H 2 and going op to a translational energy of 0.8 eV as done for H 2 + Cu (111) 45 ) 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 2 + Cu (111). In Fig. 14 we compare calculated rotational quadrupole alignment parameters for D 2 to experimental ones measured for D 2 desorbing from Cu (111). 49 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 2 and D 2 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. 25 With increasing translational energy the rotational quadrupole alignment parameter approaches zero since all molecules, irrespective of their orientation, will have enough energy to react. 49 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 2 + Cu (111).
Note that the B86SRP68-DF2, PBEa57-DF2 11 and MS-B86bl 12 DFs are in good agreement with each other for both the (n = 0,J = 11) state (Fig. 14a) and the (n = 1,J = 6) state (Fig. 14b), but that these three DFs predict slightly higher rotational quadrupole alignment parameters than the SRP48 DF. 26 Given that the SRP48 rotational quadrupole alignment parameters were decreased, but still somewhat too large when surface atom motion was introduced, 26 the present results suggest that the SRP48 DF yields the best description of this observable.

Inelastic scattering of H 2 from Cu(111)
In this section we will discuss inelastic scattering results for H 2 , obtained with QD. We start with the vibrationally inelastic scattering results for H 2 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 163 has been shown not to hold. 24 Here we discuss the inelastic scattering probability P(n = 0, Jn = 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 (n = 0,J) to (n = 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(n = 1, J = 0n = 1, J = 2)/P(n = 1, J = 0n = 1, J = 0) shown in Fig. 16 it is clear that the B86SRP68-DF2 DF performs not as good as the SRP48 DF. 10 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. 54 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. 54 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.

QD vs. QCT for H 2 + 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 2 + Cu(211), 36 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 4 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 2 reacting on Cu(211), 36 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 2 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. 47 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 2 + Cu (111). 47 The translational energy range used in our calculations overlaps with the translational energies at which the slow channel reactivity was observed. 47 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 2 + Cu(111), as we did before for H 2 + Cu(211). 36 We therefore propose, as done earlier for H 2 + Cu(211), 36 that the slow reaction channel reported by Kaufmann et al. 47 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][149][150] Galparsoro et al. 75 likewise did not yet find evidence for the slow reaction channel in their AIMD calculations.

Overall description of systems
When looking at the H 2 + Ag (111) and H 2 + Au(111) systems considered in this work together, one stark realization is that further development of chemically accurate DFs for H 2 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 164 DF for H 2 + Cu (111). 12 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 2 on transitions metals there is now only a single well studied system, namely H 2 + Cu(111) (and maybe H 2 + Pt(111) 65 ). 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 2 + Pd (111) and H 2 + Ag (111) or H 2 + 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 2 (D 2 ) + Cu (111). The H 2 + 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 (n, J) parameters obtained from reported DFMD data 62 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 0 (n, 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][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61]165,166 and theoretical work 10,12,18,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. 12 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 2 + 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 (n, 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 PBEa57-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 2 + Ag (111). For the D 2 + 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. 19 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. 12 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 2 + Ag(111) system. More and better defined experiments will allow us to refine our description of this system. 4.7.3 H 2 (D 2 ) + Au (111). With respect to the H 2 + 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 (n, 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. 70 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. 70 4.7.4 D 2 + Pt (111). Ghassemi et al. 11 have previously designed a SRP-DF for the D 2 + Pt(111) system. Although the B86SRP68-DF2 DF could describe the experiments of Luntz et al. 64 to within chemical accuracy, the description of the experiments of Cao et al. 60 was not as good. As was the case for the H 2 + 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. 65 Overall this system is best described by the PBEa57-DF2 DF that was specifically designed for this system. 11

Transferability
So far SRP-DFs fitted to reproduce molecular beam adsorption experiments for H 2 and D 2 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 2 reacting on Cu (111) can also describe the non-activated early barrier system of D 2 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 4 dissociation on Ni (111) 13 to CH 4 dissociation on Pt (111). 16 The SRP48 DF for H 2 + Cu (111) is not transferable to the H 2 + Ag(111) system 19 or to the H 2 + Pt(111) system. We have previously shown that a SRP-DF based on the mGGA that does not include non-local correlation (MS-PBEl 12 ) greatly improves the transferability from H 2 + Cu (111) to H 2 + Ag(111), 12 but Fig. 8 shows that this DF is not transferable to the weakly activated H 2 + 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 2 + Pt(111), PBEa57-DF2, 11 can describe H 2 + Cu(111) with overall chemical accuracy and that a new SRP-DF fitted to H 2 + Cu(111), B86SRP68-DF2, can describe the D 2 + Pt(111) experiments of Luntz et al. 64 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, 132 which at this accuracy level might be important. 169,170 Both the B86SRP68-DF2 and PBEa57-DF2 DFs more or less predict the same reactivity for the H 2 + 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 0 and E 1/2 (n, 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. 22 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 (n, 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 2 + Cu(111) we know from DFMD calculations 62 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 (n, 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 2 + 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][28][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 2 + 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 2 + Ag(111) system would allow for a systematic investigation into the effect of ehp excitations on reactivity for highly activated late barrier reactions.

Conclusions
We have constructed new SRP-DFs that include non-local correlation for the H 2 (D 2 ) + Cu(111) system and assessed the transferability of these DFs to the H 2 (D 2 ) + Ag (111), H 2 (D 2 ) + Au (111) and H 2 (D 2 ) + 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 2 (D 2 ) 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 2 + Cu(111), H 2 + Ag (111) and H 2 + Au (111).
SRP-DFs that include non-local correlation, namely B86SRP68-DF2 and PBEa57-DF2, are transferable from the highly activated late barrier H 2 + Cu(111) system to the weakly activated earlier barrier H 2 + 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 2 + Ag (111) and H 2 + 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 2 + 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 (n, 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 2 + 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 2 + 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 2 + Cu(111) is very well described quasiclassically 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 2 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.