Chemically accurate simulation of dissociative chemisorption of D 2 on Pt(1 1 1)

Using semi-empirical density functional theory and the quasi-classical trajectory (QCT) method, a specific reaction parameter (SRP) density functional is developed for the dissociation of dihydrogen on Pt(1 1 1). The validity of the QCT method was established by showing that QCT calculations on reaction of D2 with Pt(1 1 1) closely reproduce quantum dynamics results for reaction of D2 in its rovibrational ground state. With the SRP functional, QCT calculations reproduce experimental data on D2 sticking to Pt(1 1 1) at normal and off-normal incidence with chemical accuracy. The dissociation of dihydrogen on Pt(1 1 1) is non-activated, exhibiting a minimum barrier height of −8 meV.


Introduction
The availability of accurate barriers for reactions of molecules on metal surfaces is of central importance to chemistry. Catalysis is used to make more than 80% of the chemicals produced worldwide [1], and the accurate calculation of the rate of a heterogeneously catalyzed process requires accurate barriers for the elementary surface reactions involved [2]. This is especially true for the rate controlling steps [3,4], which often are dissociative chemisorption reactions.
Chemistry would thus benefit enormously from the availability of implementations of first principles methods that would enable the chemically accurate (i.e., to within 1 kcal/mol) calculation of barriers for reactions of molecules with metal surfaces. However, presently such implementations do not yet exist [5]. Also, density functional theory (DFT) using functionals at the gradient approximation (GA) or meta-GA level, which can be used to map out potential energy surfaces (PESs) for molecules interacting with metals, is not yet capable of predicting reaction barriers for gas phase reactions with chemical accuracy [6]. This accuracy problem of DFT is reflected in the limited accuracy with which absolute rates of heterogeneously catalyzed processes over model catalysts can now be computed with empirically optimized density func-tionals (e.g., 2 orders of magnitude for ammonia production over Ru catalysts [7]).
Currently, the most viable route to chemically accurate barriers for molecules with metal surfaces [5] uses implementations [8,9] of specific reaction parameter DFT (SRP-DFT [10]). In this semiempirical version of DFT, usually a single adjustable parameter in the density functional is fitted to reproduce an experiment that is particularly sensitive to the reaction barrier height for the specific system considered. Next, the quality of the functional is tested by checking that the candidate SRP density functional for the system also reproduces other experiments on the same system, which differ from the experiment the functional was fitted to in a nontrivial way [8,9]. Using SRP-DFT we have recently started with an effort to develop a database of chemically accurate barriers for molecules reacting with metals, which can be used to benchmark implementations of first principles methods with a claim to chemical accuracy. This database now contains data for H 2 + Cu(1 1 1) [8], H 2 + Cu(1 0 0) [11], and CH 4 + Ni(1 1 1) [9].
Here, we fit an SRP density functional for H 2 + Pt(1 1 1) to dissociative chemisorption probabilities for D 2 + Pt(1 1 1) obtained from molecular beam measurements performed at normal incidence by Luntz et al. [15]. The quality of the functional is confirmed by showing that the functional also allows reaction probabilities to be reproduced with chemical accuracy for experiments performed at off-normal incidence [15]. This is a non-trivial result, as the reaction probability for D 2 + Pt(1 1 1) does not obey normal energy scaling [15], i.e., it also depends on the component of the incidence energy parallel to the surface. This dependence arises from a particular type of correlation between the height of the barriers and their distance to the surface [23], the lowest barrier being furthest from the surface [25]. In view of the successes previously achieved for systems exhibiting a van der Waals well affecting the reactivity [9,30], we adopt a SRP density functional in which the correlation functional [31] allows at least a qualitatively accurate description of the attractive part of the van der Waals interaction. The PBEa exchange functional [32] was adopted, which allows one not only to interpolate between the well-known RPBE [33] and PBE [34] functionals, but also between PBE and a functional approximating the Wu-Cohen (WC) functional [35], which turned out to be important for the present case.
This paper is set up as follows. In Section 2.1 we describe the dynamical model we used, and in Section 2.2 how the PES for H 2 + Pt(1 1 1) was obtained. Section 2.3 describes the dynamics methods employed, and Section 2.4 gives computational details. Section 3.1 describes the PES obtained with the SRP density functional. Section 3.2 considers the accuracy of the quasi-classical trajectory (QCT) method [36] with the PES employed, and the accuracy that might be achieved by performing dynamics calculations only for the rovibrational ground state of D 2 , rather than performing a complete molecular beam simulation. In Section 3.3 we discuss how a candidate SRP density functional was derived for H 2 + Pt(1 1 1) through comparison to normal incidence data. In Section 3.4 we confirm the quality of the SRP functional through comparison of calculated sticking probabilities with experiments performed for off-normal incidence. Section 4 presents our conclusions and a brief outlook.

Dynamical model
The calculations use the so-called Born-Oppenheimer static surface (BOSS) model [8]. As discussed in for instance Ref. [29], this model allows accurate calculations on reactive scattering of H 2 from metal surfaces. With the model, the calculation of reaction probabilities is split in two parts: First, the PES is calculated (Section 2.2), and next the PES is used in dynamics calculations (Section 2.3). In the PES and the dynamics calculations, only the six molecular degrees of freedom of the H 2 molecule are taken into account. The coordinates to describe the motion of the molecule are shown in Fig. 1a.

Calculation of the PES
The ground state PES was calculated using DFT. The exchangecorrelation functional used to compute the PES may be written as i.e., we use the PBEa exchange functional [32], with a being the adjustable parameter, and the van der Waals DF2 functional of Langreth and Lundqvist and coworkers [31]. With the choice a = 1 the PBEa functional corresponds [32] to the PBE functional [34], while for a ! 1 the PBEa functional corresponds [32] to the RPBE functional [33]. For a = 0.52 a functional is obtained that closely resembles [32] the WC functional that performs well in solid state calculations [35]. The use of PBEa in semi-empirical applications would seem to be especially advantageous if interpolation is required between PBE and a less repulsive exchange functional; if the goal is to interpolate between PBE and RPBE exchange we suggest using a weighted average of these two [9,37], as using a ! 1 in PBEa to obtain the RPBE limit is a bit awkward for this purpose.
To obtain a global expression for the PES, the accurate corrugation reducing procedure (CRP) [38] was used to interpolate points calculated on a grid with DFT. The procedure used is exactly the same as used in Ref. [29]. The p3m1 plane group symmetry [39] associated with the Pt(1 1 1) surface was used.

Dynamics calculations of reaction probabilities
Reaction probabilities were calculated for the (v = 0, j = 0) state of D 2 with the time-dependent wave packet (TDWP) method [40] in an implementation for dihydrogen scattering from surfaces with hexagonal symmetry that is fully described in Ref. [25]. Dissociation probabilities of D 2 colliding with Pt(1 1 1) for comparison with molecular beam experiments on the same system [15] were calculated with the QCT method [36] in an implementation described in Ref. [29]. Earlier calculations predicted that even for the lighter H 2 molecule the QCT method yields dissociative chemisorption probabilities for hydrogen dissociation on Pt(1 1 1) that are in excellent agreement with quantum dynamics results [24]. For the best comparison with experiments, the calculations include Monte-Carlo averaging over the velocity distributions of the hydrogen beams, and Boltzmann averaging over the rovibrational states of hydrogen, as fully described in Ref. [29]. An important assumption made in our calculations is that the molecular beams used in the experiments of Luntz et al. [15] are quite similar to hydrogen beams produced in experiments of Juurlink and co-workers [41], and we used the beam parameters presented in table 3 of Ref.
[30] to simulate D 2 beams in our work on the basis of this assumption.

Computational details
The DFT calculations were performed with the VASP (version 5.2.12) programme [42][43][44]. Standard projector augmented wave (PAW) potentials [45] were used. First, the bulk fcc lattice constant was determined in the same manner as used previously for H 2 + Au (1 1 1) [46], using a 20 Â 20 Â 20 C-centered grid of k-points. With the optimized SRP density functional (using a = 0.57, see Section 3) a lattice constant of 4.015 Å was obtained, in reasonable agreement with the experimental value of 3.91 Å. Next, a relaxed 5-layer slab was obtained, again in the same manner as used before for H 2 + Au (1 1 1) [46], using a 20 Â 20 Â 1 C-centered grid of k-points. After having obtained the relaxed slab, single point calculations were carried out on H 2 + Pt(1 1 1), using a 9 Â 9 Â 1 C-centered grid of k-points, and a plane wave cut-off of 400 eV, in a super cell approach in which 13 Å of vacuum length was used for the spacing between the Pt(1 1 1) slabs. The grid of the points for which the H 2 + Pt(1 1 1) calculations were done, and other details of the calculations, were taken the same as in Ref. [29]. The CRP PES was extrapolated to the gas-phase potential of H 2 in the same way as used in Ref. [29].
In the QCT calculation of dissociative chemisorption probabilities for comparison with molecular beam experiments, 10,000 trajectories were run for each (v,j) state with the vibrational quantum number v 6 3 and the rotational quantum number j 6 20. For each j, uniform sampling was performed of the magnetic rotational quantum number m j . The centre-of-mass of H 2 was originally placed at Z = 9 Å, with the velocity directed towards the surface and sampled from appropriate velocity distributions for D 2 beams (see Table 3 of Ref. [30]). The molecule is considered dissociated once r > 2.25 Å, and considered scattered once Z > 9 Å. Other computational details of the QCT calculations are the same as in Ref. [29]. The surface lattice constant (i.e., the nearest neighbor Pt-Pt distance) used in the QCT calculations (and in the TDWP calculations) was taken as the computed Pt lattice constant divided by ffiffiffi 2 p (i.e., as 2.84 Å).
In the TDWP calculations on (v = 0, j = 0) D 2 + Pt(1 1 1), two separate wave packet calculations were performed to cover the collision energy range E i = 0.05-0.55 eV. This procedure avoids problems that may arise from the interaction of the low energy components of the wave packet with optical potentials if only one broad wave packet is used to cover a very large translational energy range. The input parameters we used in the TDWP calculations are listed in Table S1. Convergence tests carried out suggest that, with the use of these parameters, the reaction probabilities computed for (v = 0, j = 0) D 2 are converged to within better than 2% of their values (i.e., relative errors 62%).

Potential energy surface
Two-dimensional cuts (so-called elbow plots) through the PES used in the dynamics calculations on H 2 + Pt(1 1 1) are shown in Fig. 2, in all cases for H 2 oriented parallel to the surface. With the optimized SRP density functional (using a = 0.57, see Section 3.3), the dissociation is non-activated in the sense that the transition state has an energy that is 8 meV below the gas phase minimum energy of H 2 (the early barrier for dissociation above the top site, see also Table 1, which lists the geometries and barrier heights corresponding to the results shown in Fig. 2). With the functional used, the barrier height (E b ) shows a larger energetic corrugation (i.e., a greater variation with impact site) than previously obtained with the Becke-Perdew functional (Ref. [25] and references therein). This is what should be expected for a functional accurately describing the experiments on dissociative chemisorption of D 2 on Pt(1 1 1) [15], as the previously computed sticking probability vs. incident translational energy curve was too steep [18,25]. Note that previous experience with H 2 -metal systems suggests that the use of Lundqvist-Langreth van der Waals correlation functionals, as employed here, yield PESs with larger energetic corrugation than ordinary GGA correlation functionals [29,30].   Fig. 3 shows a plot of the potential at r = r e , after averaging over X, Y, h, and /, with r e being the minimum H-H distance of gasphase H 2 . This averaged potential curve shows a van der Waals minimum well depth of 72 meV, in excellent agreement with the range of values found in experiments (i.e., 55 meV [14,18], and 76 meV [47]). Getting the van der Waals attractive interaction right may be important to obtaining a correct value for the energy of the ''early" transition state (which occurs at Z = 2.2 Å, see Table 1) and is probably also important to the calculation of probabilities for diffractive scattering, for which detailed experimental results are available [18].
3.2. Quantum vs. quasi-classical dynamics, and the importance of simulating the molecular beam Fig. 4a shows a comparison of reaction probabilities computed for D 2 in its initial (v = 0, j = 0) state for specific incidence energies with quantum dynamics and with quasi-classical dynamics. The calculations used the optimized SRP density functional (i.e., with a = 0.57, see Section 3.3). Even in the absence of averaging over initial rovibrational states and over the distribution of energies, as would be appropriate for comparisons with molecular beam experiments, the quantum and QCT results are in excellent agreement with one another. In the following, we will therefore use the QCT method to compute sticking probabilities for comparison to the molecular beam experiments of Luntz et al. [15] (see Fig. 5). Fig. 4b shows a comparison of reaction probabilities computed with the QCT method for D 2 in its initial (v = 0, j = 0) state for specific incidence energies with QCT results obtained with full averaging over the rovibrational state populations and velocity distributions that are typical for molecular beam experiments using pure D 2 beams [30,41]. The comparison of Fig. 4b suggests that it should not really be necessary to take the effect of the velocity distribution and the rovibrational state distribution into account, in broad agreement with an earlier theoretical study of H 2 + Pt(1 1 1) [27]. This is in sharp contrast with findings for the highly activated H 2 + Cu(1 1 1) reaction [8,48]; for this system, taking into account the velocity distribution is necessary for accurate results, because the reactivity may come entirely from incidence energies above the average incidence energy of the beam, and above the high reaction threshold. Even though taking into account the beam conditions should be much less important for D 2 + Pt (1 1 1), in the following we will always represent computational results with full averaging over the incidence energy and rovibrational state population of the D 2 beams, to obtain the best possible comparison with the molecular beam experiments of Luntz et al. [15].

Fit of the SRP density functional to molecular beam data for normal incidence
Before the SRP functional for H 2 + Pt(1 1 1) could be fit, a choice had to be made concerning which experimental dataset for normal Table 1 Barrier heights (E b ), the distance to the surface of the barrier (Z b ), and the H-H distance at the barrier (r b , in Å) are given for four different dissociation geometries defined by the impact site and the angle / (see Fig. 1), for dissociation of H 2 over Pt(1 1 1) with H 2 parallel to the surface (h = 90°). The results have been obtained with the PBEa-vdW-DF2 functional with a = 0.57. For the top site, results are given for two barrier geometries. The E b values in brackets correspond to the 6D PES computed with the Becke-Perdew functional (see Ref. [25]).   incidence theoretical results should be compared with. In the literature, two sets of molecular beam data are available for dihydrogen normally incident on Pt (1 1 1), i.e., those of Luntz et al. [15] and those of Samson et al. [16]. The work of Luntz et al. focused on the dihydrogen + Pt(1 1 1) system, looking at the effects of the angle of incidence h i , surface temperature T s , isotopic mass, and nozzle temperature T n in great detail, and producing data for D 2 + Pt(1 1 1) at T s = 300 K for a large range of incidence energies E i by also using seeding of D 2 in H 2 to achieve high E i . In contrast, Samson et al. only published data for D 2 + Pt(1 1 1) for normal incidence, for one value of T s (150 K), for the more limited range of E i available with pure D 2 beams only, in a paper focused on how alloying varying amounts of Sn into the surface affects the sticking. Furthermore, Luntz et al. explicitly stated that their ''incidence energies" (labeled E i in their work) were energy averaged over the TOF distribution of the beams they used, whereas Samson et al. simply assumed that the average incidence energy (which we will label as hE i i) is given by hE i i ¼ 2:75kT n . For these reasons, we have chosen to fit our SRP functional to the normal incidence data of Luntz et al., assuming that these would represent the most accurate dataset.
The  [8]) this should be reflected in the accuracy of the extracted SRP functional and minimum barrier height. Problems with the interpretation of results of molecular beam experiments due to lacking or incomplete specification of the velocity distributions have hampered efforts to obtain accurate SRP functionals and benchmark data before [49]. However, the problem noted here for H 2 + Pt(1 1 1) is not as severe as for H 2 + Pd(1 1 1) [49].
To obtain an SRP functional, first tests were performed combining the PBE functional for exchange [34] with the Lundqvist-Langreth functional of Dion et al. (vdW-DF1) [50]. With this functional, the van der Waals well was too deep compared with experimental results, and the computed reaction probabilities were shifted to too high energies and did not exhibit chemical accuracy (results not shown). For these reasons, we switched to the improved Lundqvist-Langreth functional of Lee et al. (vdW-DF2) [31], and to the PBEa functional [32], adjusting a by trial and error to obtain agreement with the sticking experiments of Luntz et al. [15]. By choosing a = 0.57, agreement with the experiments for normal incidence could be obtained to within chemical accuracy, by which we mean that the computed sticking probabilities are displaced along the energy axis from the interpolated experimental curve by no more than 1 kcal/mol (see Fig. 5). The resulting SRP-DFT PES shows a minimum barrier height of À8 meV (%1 kJ/mol), suggesting the reaction to be non-activated if the molecule hits the surface at the right site (the top site, see Table 1). The ''activated appearance" of the reaction probability curve comes from the molecule also hitting the surface at other impact sites and orientations for which higher barriers are encountered (see for instance Table 1 and Fig. 2), as already suggested by Luntz et al. at the time of their work [15].

Confirming the quality of the SRP density functional by comparison to molecular beam data for off-normal incidence
Strictly speaking, the functional obtained in Section 3.4 is, at this stage, only a ''candidate SRP functional": to become an SRP functional, dynamics calculations with the functional should also be able to reproduce other experiments on the same system, which differ from the experiments the functional was fit to in a nontrivial way [8,9]. For this, we chose to use the datasets obtained by Luntz et al. for off-normal incidence [15]. For H 2 + Pt(1 1 1) more recent, detailed data on molecular diffraction are also available [18], but recent work on H 2 + Ru(0 0 0 1) suggests that accurately reproducing diffraction data is fraught with difficulties [30]. This is most likely related to the need to extrapolate the experimental data to a low temperature or even static surface regime using Debye-Waller attenuation, or to simulate the effect of surface temperature in accurate quantum dynamics calculations [30].
Reaction probabilities computed for h i = 30 and 45°agree with the experimental values to within chemical accuracy (Fig. 6). Larger displacements than 1 kcal/mol of the computed reaction probabilities from the interpolated experimental sticking curve are observed for h i = 60°, but we argue that for this large an incidence angle our operational definition of chemical accuracy may not be appropriate. The slope of the measured sticking curve as a function of total incidence energy is small, so that a small error in the measured reaction probability could have a large effect of the energy displacement of the computed reaction probability to the interpolated experimental curve. In this context, we note that error bars on the measured sticking probabilities were lacking [15]. In view of the positive results for h i = 30 and 45°, we argue that our PBEa-vdW-DF2 functional is an SRP functional for H 2 + Pt (1 1 1), and that the minimum barrier data (and the barriers obtained for other impact sites shown in Table 1) can be used for benchmark purposes, i.e., they can be included in an emerging database with chemically accurate barriers for molecules interacting with transition metals [5].
Luntz et al. did not specify the incidence plane used in their experiments on off-normal incidence [15]. The computed data shown in Fig. 6 are for incidence along the h1 1 À2i direction, which corresponds to the vector bisecting the u and v vectors in Fig. 1. However, for incidence along the h1 0 À1i direction (corresponding to the direction of u in Fig. 1), the computed sticking probabilities closely reproduce the values computed for the h1 1 À2i direction, and they likewise reproduce the experimental data, for h i = 30 and 45°(not shown). For h i = 60°and incidence along the h1 0 À1i direction, the computed sticking probabilities do not quite reproduce the values computed for the h1 1 À2i direction (in agreement with earlier theoretical work on H 2 + Pt(1 1 1) [25]), but the result that for this incidence direction and large angle the computed data do not reproduce the experiments with chemical accuracy is also obtained for the h1 0 À1i direction.

Conclusions and outlook
We have obtained an SRP density functional for H 2 + Pt(1 1 1) by adjusting the a parameter in the PBEa-vdW-DF2 functional until reaction probabilities computed with the QCT method reproduced sticking probabilities measured for normally incident D 2 with chemical accuracy. In the QCT calculations, the rovibrational state populations and the velocity distributions of the incident beams were taken into account. Also, the appropriateness of the use of the QCT method for the purpose of accurately calculating reaction probabilities for D 2 + Pt(1 1 1) was established by a comparison with quantum dynamics calculations for the initial (v = 0, j = 0) state of D 2 . The quality of the SRP functional was confirmed by showing that QCT calculations using the functional also reproduced data for off-normal incidence for h i = 30 and 45°, for which the computed reaction probabilities show no dependence on the plane of incidence. The minimum barrier height obtained for the reaction is À8 meV, in agreement with the experimental observation of no, or only a small energetic threshold to reaction [15]. This value can be entered into a small [5], but growing [9] database with barriers of reactions of molecules with metal surfaces, for which chemical accuracy is claimed. Our conclusion depends on the assumption that the data of Luntz et al. are accurate and the validity of our interpretation of the average incidence energy in their experiments [15] (see Section 3.3). To confirm this, accurate new experiments on reaction of H 2 or D 2 with Pt(1 1 1) for varying incidence angles and well-defined molecular beam velocity distributions and incidence plane would be welcomed. Such experiments might also be able to confirm or falsify our prediction [25] that significantly different reaction probabilities should be obtained for different incidence planes and large incidence angles h i .
Future computational work could address the question of how the dynamical model may have to be extended to accurately reproduce the detailed molecular diffraction data available for H 2 + Pt (1 1 1) [18]. Once such a model is available it could be used in further tests of the SRP density functional for H 2 + Pt(1 1 1) by comparison to these data, and in tests of the candidate SRP density functional for H 2 + Ru(0 0 0 1), for which detailed diffraction data are also available [30]. We also suggest that the SRP functional be used to model data on the reaction of H 2 with stepped Pt surfaces [20,21], to check whether SRP functionals developed for a low index transition metal surface exhibit transferability to systems in which the same molecule interacts with a vicinal or stepped surface of that metal.   Fig. 6. Reaction probabilities computed for D 2 + Pt(1 1 1) with the SRP density functional (see text) are shown as a function of hE i i, comparing to the molecular beam results of Luntz et al. [15]. The results are for off-normal incidence at the indicated incidence angles h i of 30, 45 and 60°, along the h1 1 À2i incidence direction. The arrows and accompanying numbers show the collision energy spacing (in kJ/mol, 1 kcal/mol % 4.2 kJ/mol) between the computed sticking probabilities and the interpolated experimental sticking probability data (green circles).
(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)