Transferability of the Specific Reaction Parameter Density Functional for H2 + Pt(111) to H2 + Pt(211).

The accurate description of heterogeneously catalyzed reactions may require chemically accurate evaluation of barriers for reactions of molecules at the edges of metal nanoparticles. It was recently shown that a semiempirical density functional describing the interaction of a molecule dissociating on a flat metal surface (CHD3 + Pt(111)) is transferable to the same molecule reacting on a stepped surface of the same metal (Pt(211)). However, validation of the method for additional systems is desirable. To address the question whether the specific reaction parameter (SRP) functional that describes H2 + Pt(111) with chemical accuracy is also capable of accurately describing H2 + Pt(211), we have performed molecular beam simulations with the quasi-classical trajectory (QCT) method, using the SRP functional developed for H2 + Pt(111). Our calculations used the Born-Oppenheimer static surface model. The accuracy of the QCT method was assessed by comparison with quantum dynamics results for reaction of the ro-vibrational ground state of H2. The theoretical results for sticking of H2 and D2 on Pt(211) are in quite good agreement with the experiment, but uncertainties remain because of a lack of accuracy of the QCT simulations at low incidence energies and possible inaccuracies in the reported experimental incidence energies at high energies. We also investigated the nonadiabatic effect of electron-hole pair excitation on the reactivity using the molecular dynamics with the electron friction (MDEF) method, employing the local density friction approximation (LDFA). Only small effects of electron-hole pair excitation on sticking are found.

Introduction

The heterogeneous catalysis community is highly interested in stepped surfaces because structure-sensitive-catalyzed reactions often occur at edges of nanoparticles. These edges contain low-coordinated surface atoms, which resemble the atoms present at step edges of stepped surfaces. Consequently, a number of experiments have addressed dissociative chemisorption reactions of molecules on stepped surfaces, such as NO at steps on defective Ru(0001),[1] H2 on stepped Pt surfaces,[2−8] N2 at steps on defective Ru(0001),[9,10] and methane on Pt surfaces,[11,12] to name but a few examples. A much lower number of theoretical dynamics studies have addressed dissociative chemisorption on stepped surfaces, and these studies have looked at H2 + Pt(211),[13−17] H2 + Cu(211),[18,19] H2 dissociation on defective Pd(111),[20] and at CHD3 + Pt(211).[12,21−24]

In view of the importance of dissociative chemisorption reactions on stepped surfaces to heterogeneous catalysis, it would obviously be useful to have a predictive procedure in place for accurately evaluating the interaction between a molecule and a stepped surface. Recent experimental work suggests that such a procedure may be based on experiments and dynamics calculations based on the semiempirical density functional theory (DFT) for the electronic structure, for the same molecule interacting with a low-index, flat surface of the same metal.[12] As has now been established for several systems, dynamics calculations based on electronic structure calculations with the specific reaction parameter approach to DFT (SRP–DFT) are able to reproduce sticking measurements on such systems with chemical accuracy.[12,25−28] Very recently, it has been shown that the SRP density functional (SRP-DF) for CHD3 interacting with the flat Pt(111) system is transferable to the same molecule interacting with the stepped Pt(211) system[12] (transferability of the SRP DF from H2 + Cu(111)[12] to H2 + Cu(100),[26] that is, among systems in which the same molecule interacts with different flat, low-index surfaces, had been established earlier[26]). However, this finding just concerned only one specific system, and it is important to check whether this finding also holds for other systems. The main goal of this work is to investigate whether the SRP-DF recently determined for H2 + Pt(111)[28] is also capable of yielding chemically accurate results for H2 + Pt(211).

The system of interest to our study (H2 + Pt(211)) has first been studied theoretically. Olsen et al.[13] computed a six-dimensional (6D) potential energy surface (PES) for the system with DFT, using the GGA functional due to Becke[29] and Perdew[30] (BP), and interpolating the DFT results with the corrugation reducing procedure (CRP).[31] They next performed classical trajectory studies on this PES within the Born–Oppenheimer and static surface (BOSS) approximations. On the basis of these calculations, they were able to show that a trapping mechanism contributes a component to the sticking probability which is high at low incidence energy (Ei) and decreases monotonically with Ei.[13] In this mechanism, H2 gets trapped at an unreactive site, that is, at the bottom of the step, and then diffuses to an atom at the top of the step edge, where it subsequently reacts.

Next, McCormack et al. also analyzed the other contributing mechanisms to the sticking of H2 on Pt(211).[14] Their classical trajectory calculations using the same PES as used before showed two additional mechanisms. A mechanism in which H2 reacts directly at the step is nonactivated and contributes equally at all Ei. In an additional mechanism, H2 reacts on the terrace. In this mechanism, the reaction is activated, yielding a contribution to the sticking that rises monotonically with increasing Ei. By scaling the contributions from the different mechanisms according to the different lengths of the (111) terraces in the Pt(211) and Pt(533) surfaces (both exhibiting (111) terraces and (100) steps), they[14] were able to obtain good agreement with previous experiments on H2 + Pt(533).[6]

In two subsequent studies using the same PES, Luppi et al.[15] investigated rotational effects with classical trajectory calculations, whereas Olsen et al.[17] made a comparison between quantum dynamics (QD) and classical dynamics results for the reaction of (ν = 0, j = 0) H2. According to the classical trajectory studies of Luppi et al., the trapping-mediated contribution to the reaction, which leads to a high sticking probability at low Ei, but which contribution then quickly decreases with Ei, should be present for low rotational states (j = 0 and 1), but should disappear for states with intermediate j. The reason they provided is that energy transfer to rotation should cause trapping for j = 0 and j = 1, whereas energy transfer from rotation should instead hinder trapping. Olsen et al. found that quasi-classical trajectory (QCT) calculations were in good agreement with QD results for high Ei (in excess of 0.1 eV). However, the QCT study overestimated the trapping-mediated contribution to the reaction at low Ei, which was attributed to one mechanism operative for trapping in the classical calculations (excitation of the rotation) not being allowed in QD, as the trapping well should not support rotationally excited bound states for their PES.[17]

H2 + Pt(211) has also been studied experimentally by Groot et al.[7,8,32] Their molecular beam sticking probabilities[7] were in reasonable agreement with the QD results for (ν = 0, j = 0) H2 of Olsen et al.,[17] although the QD results based on the BP functional overestimated the sticking at high Ei. Likewise, there were discrepancies at low Ei, with the computed trapping-mediated contribution to the sticking being too low compared to the experimental result. In two subsequent papers, Groot et al. showed that the sticking on surfaces with longer (111) terraces and (100) steps (Pt(533) and Pt(755)) can successfully be modeled based on the contributing mechanisms to sticking at the step and at the terrace on Pt(211),[8,32] much like McCormack et al. had done before for Pt(533).[14] They also used their results to analyze the contributions of facets and edges of Pt nanoparticles to H2 dissociation proceeding on these nanoparticles.[8]

The goal of the present paper is to test whether the SRP-DF for H2 + Pt(111) is transferable to H2 + Pt(211). For this reason, we will put emphasis on comparison of sticking probabilities computed with a PES obtained with the SRP-DF for H2 + Pt(111) with the experimental results of ref (8), taking the experimental conditions (velocity distributions of the beams, nozzle temperatures Tn used) into account as fully as possible. Our calculations are done within the BOSS model and mainly use the QCT method for the dynamics. We will not reanalyze the mechanisms contributing to the reaction, simply noting that the dependence of the computed sticking probabilities on Ei is in accordance with conclusions arrived at earlier by Olsen et al.[13] and McCormack et al.[14] We find that overall, the computed sticking probability is in good agreement with the experiment for both H2 and D2 + Pt(211), suggesting that the transferability may well hold. However, at present, this conclusion is not yet certain because of uncertainties in the parameters needed to describe the molecular beams used in the experiments. Our results suggest that once more precisely defined experimental results become available, the comparison with the experiment should be revisited on the basis of QD calculations.

This paper is setup as follows. Section describes the dynamical model, and Sections and 2.3 describe the construction of the PES and the PES interpolation method. The dynamics methods used here are explained in Sections and 2.5. Section describes how we calculate the observables. Section provides computational details. In Section , the results of the calculations are shown and discussed. Section describes the computed PES. In Section , we compare the QCT results with the QD results. The isotope effect of the QCT results for reaction of (ν = 0, j = 0) H2 and D2 is shown and discussed in Section . Section provides theoretical results on molecular beam sticking probabilities and comparison with the experimental data. In Section , the effect of electron–hole pair excitation on the reactivity is discussed and the molecular dynamics with electron friction (MDEF) results are compared with the MD results for sticking. Conclusions are provided in Section .

Theoretical Methodology

Dynamical Model

The dynamics simulations presented in the following approach the true reaction dynamics of the system by assuming the reaction to take place on an ideal rigid Pt(211) surface at zero coverage. During the entire dynamics, the surface atoms are fixed at their initial equilibrium positions as obtained from DFT calculations. The dynamical degrees of freedom (DOF) treated here are the six DOF of H2. These are the center-of-mass (COM) position given by Cartesian coordinates X, Y, Z relative to a surface atom, the interatomic H–H distance r, and the angular orientation of the molecule defined with respect to the macroscopic surface plane. As usual, X and Y are the lateral components of the COM position and Z is the molecule–surface distance. The orientation of the molecule is specified by the polar angle θ ∈ [0, π] and the azimuthal angle ϕ ∈ [0, 2π]. The corresponding coordinate system is visualized in Figure .

Figure 1

Coordinate systems for H2 on Pt(211). (a) Top view of the (1 × 1) unit cell also showing the dissociated reference geometry of H2 used to converge the computational setup with respect to the adsorption energy Eads. First and second layer Pt atoms are in silver and dark gray, respectively. H atoms are blue colored. (b) Side view of the slab model. The Z-axis (molecule–surface distance) in the standard coordinate system drawn in black is aligned with the normal to the macroscopic surface. X and Y are the lateral components of the COM position of H2 indicated by a red dot. Furthermore, r is the interatomic H–H distance (not shown), and the angular orientation is specified by the polar angle θ ∈ [0, π] and the azimuthal angle ϕ ∈ [0, 2π] (not shown). The angular orientation of H2 in the internal coordinate system is defined with respect to the normal of the (111) terrace, as shown in red. The two coordinate systems include an angle χ of 20°. The corresponding angular coordinates are {θ′, ϕ′}. The surface lattice constants are L = 6.955 Å and L = 2.839 Å.

Electronic Structure Calculations

In this work, electronic structure calculations are carried out using periodic DFT as implemented in the Vienna ab initio simulation package (VASP).[33−36] Specifically, we employ an exchange-correlation functional of the formwhich contains PBEα exchange[37] and the vdW-DF2-functional of Lundquist and Langreth and co-workers.[38] The latter accounts for long-range van der Waals (vdW) interactions. The α-value was set to 0.57 according to our previous work[28] where we have determined this value to be suitable in order to bring computed and measured[39] sticking probabilities for D2 on Pt(111) in quantitative agreement. At first sight, the strategy of fitting a DFT functional to an experiment performed on a particular system might lead to a functional that is too specific to be accurate also for other systems, even though they might appear very similar chemically. However, recent theoretical work on the dissociation of molecular hydrogen on different facets of Cu[25,26] and methane dissociation on nickel[27] and platinum[12] surfaces has shown that so-optimized functionals may indeed be transferable among different but chemically similar systems. This suggests that the SRP functional designed for the D2 + Pt(111) system might be of similar accuracy for the D2(H2) + Pt(211) system.

The DFT calculations on the D2 + Pt(211) system presented here are based on a Pt(211) slab model with four layers using a (1 × 2) supercell. As often done for hydrogen + metal systems, we here assume effects resulting from surface atom motion on the dissociation dynamics to be negligible at the relevant experimental conditions to which we will compare our simulations. Consequently, we content ourselves with a representation of the interaction potential for a frozen Pt(211) surface. The surface atom positions of the three uppermost layers are initially optimized by relaxing the Pt slab but then kept frozen for all subsequent calculations on the system. We took care that the mirror axis was not affected by the geometry optimization of the slab. The resulting slab model obeys the symmetry of the p1m1 plane group.[18] This is helpful in reducing the computational burden associated with the construction of the 6D PES, as we will show below. Similar to ref (18), the vacuum gap separating periodic slab images is about 16.2 Å. We use a Γ-centered 7 × 7 × 1 k-point mesh generated according to the Monkhorst grid scheme.[40] The energy cutoff, EPAW, used in the projector augmented wave (PAW) method was set to 450 eV. We employ Fermi smearing with a width of 0.1 eV. The optimal number of k-points and surface layers and the optimal EPAW value were determined by convergence calculations, as summarized in Table . There we list the adsorption energy Eads computed as difference between the minimum energy of H2 at its equilibrium distance req ≈ 0.74 Å in the gas phase (here about 6 Å away from the surface and parallel to the surface) and the dissociatively adsorbed configuration of H2 on Pt(211), as depicted in Figure . Eads-values are listed in Table for different slab thicknesses, k-point meshes, and cutoff energies. The lattice constants of the rectangular (1 × 1) surface unit cell are L = 6.955 Å along the X-axis and L = 2.839 Å along the Y-axis, corresponding to a bulk lattice constant D of 4.016 Å. The latter value compares reasonably well with the experimental value (D = 3.916 Å[41]).

Table 1

Adsorption Energies Eads in eV for H2 on Pt(211) Computed Using Different k-Point Meshes, Cutoff Energies EPAW, and Number of Layers in the Slaba

		four-layer slab				five-layer slab
	EPAW [eV]	350	400	450	500	350	400	450	500
k-points	5 × 5 × 1	0.951	0.940	0.934	0.931	0.951	0.939	0.934	0.931
	6 × 6 × 1	0.952	0.941	0.935	0.932	0.951	0.940	0.934	0.931
	7 × 7 × 1	0.962	0.952	0.945*	0.943	0.962	0.951	0.945	0.942
	8 × 8 × 1	0.963	0.953	0.947	0.944	0.953	0.952	0.946	0.943

The Eads-value obtained with a converged computational setup is marked by an asterisk. The reference geometry of dissociated H2 used to determine Eads is shown in Figure .

Representation of the PES

In order to construct a continuous electronic ground state PES for molecular hydrogen interacting with a rigid Pt(211) system, we adopt the CRP[31] which allows for a fast and accurate interpolation of DFT data points. The 6D PES accounts only for the six DOF of molecular hydrogen, as shown in Figure . Details about the CRP algorithm and its implementation in our in-house computer code are presented elsewhere.[18] In the following, only a few principles of the CRP will be explained and a few details will be presented concerning the structure of the DFT data set. The interpolation of realistic globally defined PESs can become considerably error-prone when small geometrical alterations lead to strong changes of the system’s potential energy. Using the CRP, this problem can be avoided by first reducing large differences within the original DFT data points, VDFT. The resulting reduced data set, IDFTis better suited for an interpolation which will yield the smooth function I(Q⃗) used to compute the final PES according toHere, Q⃗ = (X, Y, Z, r, θ, ϕ)T is a discrete coordinate vector, labeled with the multidimensional index i, in the 6D space Q⃗ = (X, Y, Z, r, θ, ϕ)T. For the reference function, Vref(Q⃗), we are here using the sum of the two H + Pt(211) interaction potentials which are also obtained via the CRP. They describe most of the repulsive features of the PES and are therefore particularly suitable for reducing the corrugation of the PES in the CRP as explained in refs.[18,31]

In order to keep the number of DFT points to be computed as low as possible, we perform DFT calculations for specific angular orientations of H2 labeled by {θ′, ϕ′} in the following. They are defined in a modified coordinate system which is aligned with the vector normal to the (111) terrace and not with the vector normal to the macroscopic surface as is the case for the angular coordinates {θ, ϕ}. The corresponding transformations between the two coordinate systems were previously presented in ref (42) and the Supporting Information of ref (18). In Tables and 3, we list details about the DFT grid representation of the PES for the H(D) + Pt(211) system as well as for the H2(D2) + Pt(211) system (see also Figure.). The latter is required to provide the reference PES Vref(Q⃗) in eqs and 3. Note that with the coordinate system chosen for the DFT calculations, for H + Pt(111), a low minimum value of Z is needed to map out the interaction of H with Pt(211) at the bottom of the step (see Table ). In the CRP, this is required in order to remove the repulsive interaction in the H2 + Pt(211) PES over the whole interpolation range before interpolation is carried out. Because of the (100) step, the surface roughness is increased and small molecule–surface distances need to be taken into account (here, Zmin = −2.2 Å). Considering that we also describe molecular configurations in which H2 stands perpendicular to the surface and that we represent large interatomic distances (rmax = 2.5 Å), atomic repulsions must then also be represented for small atom–surface distances, down to Z = −3.45 Å (Figure.).

Table 2

Specification of the DFT Grid Used To Represent the Atomic Reference H + Pt(211) Interaction Potentiala

quantity	value	unit	remark
grid range along X on IW	[0, LX]	Å
grid range along Y on IW	[0, LY/2]	Å
grid range along Z	[−3.65, 7.05]	Å
NX number of grid points in X on IW	18		equidistant
NY number of grid points in Y on IW	3		equidistant
NZ number of grid points in Z	109		equidistant
ΔX grid spacing along X	LX/18	Å
ΔY grid spacing along Y	LY/4	Å
ΔZ grid spacing along Z	0.1	Å
representation of Vtop reference potential
grid range along Z	[0, 7.05]	Å
NZtop number of grid points in Z	576		non-equidistant

The grid along Y is defined for the irreducible wedge (IW) which makes up only the half of the Pt(211) (1 × 1) unit cell (see Figure ).

Table 3

Specification of the DFT Grid Used To Represent the H2(D2) + Pt(211) Interaction Potentiala

quantity	value	unit	remark
range of X	[0, LX]	Å
range of Y	[0, LY/2]	Å
range of Z	[−2.2, 6.6]	Å
range of r	[0.4, 2.5]	Å
range of θ′	[0, π/2]	rad
range of ϕ′	[−π/4, π/2]	rad
NX number of grid points along X	9		equidistant
NY number of grid points along Y	3		equidistant
NZ number of grid points along Z	53		equidistant
Nr number of grid points along r	22		equidistant
Nθ′ number of grid points along θ′	2		equidistant
Nϕ′ number of grid points along ϕ′	3–4(*)		equidistant
ΔX grid spacing of X	LX/9	Å
ΔY grid spacing of Y	LY/4	Å
ΔZ grid spacing of Z	0.15	Å
Δr grid spacing of r	0.1	Å
Δθ′ grid spacing of θ′	π/2	rad
Δϕ′ grid spacing of ϕ′	π/4	rad

The grid along Y is specified for the IW which equals here the lower half of the Pt(211) (1 × 1) unit cell (see Figure ). Because of symmetry, the ϕ′ dependence of the PES along the top and the brg line can be represented with three points (here at ϕ′ = 0, 45 and 90°). Because of the absence of a mirror axis associated with the t2b line, we needed an additional point (here at ϕ′ = 315°) to sample the PES along ϕ′.

Figure 2

Top view of a (1 × 1) unit cell of Pt(211). Indicated is the IW by a blue plane, and the blue dots represent the positions of H and of the center of mass of H2, at which DFT energy points were calculated in order to construct the 3D/6D PES. A few selected sites are labeled with top, brg (bridge), and t2b (top to bridge) and are further distinguished by numbers. Red dots indicate periodic images at the edge of the IW.

We apply the following interpolation order to generate a smooth function IDFT(Q⃗). First, we interpolate along the interatomic H–H distance r and the molecule–surface distance Z using a two-dimensional spline interpolation. Second, we interpolate along the polar angle θ′ using a trigonometric interpolation. Finally, we interpolate along the lateral positions X, Y and the azimuthal angle ϕ′ using a symmetry-adapted three-dimensional Fourier interpolation. The resulting PES is smooth, fast to evaluate and provides analytical forces.

MD Simulations

In this work, the dissociation dynamics of molecular hydrogen on Pt(211) is modeled using the QCT method,[43] that is, with MD simulations. The quantum mechanical ro-vibrational energy of incident H2/D2 is sampled by a Monte-Carlo procedure outlined in ref (44), and the occupation of the associated ro-vibrational levels is determined by the molecular beam parameters, as discussed below. We distinguish between standard MD simulations and MDEF.[45] The latter method allows one to study nonadiabatic effects on the dissociation dynamics because of the creation of electron–hole pairs in the surface region. For a N-dimensional system, the general equation to be solved in the following is the Langevin equation[46] which readsHere, m is the mass associated with a generalized coordinate q, and η is an element of the friction tensor which yields a dissipative term due to the coupling of the nuclear DOF of molecular hydrogen with the electronic DOF of the Pt(211) surface. Finally, R(T) is a white noise random force resulting from the electronic bath at temperature T ≔ Tel = TS, which here corresponds to the surface temperature TS. At T = 0 K, the random force disappears and only the frictional force remains in the dissipative part of eq . In the absence of electronic friction (η = 0), the Langevin equation obeys Newton’s equation of motion and the evolution of the system depends then only on the gradient of the PES. The methodology used to solve eq is described in refs.[44,47]

The position-dependent friction coefficients in eq are computed using the local-density friction approximation (LDFA) with the use of the independent atom approximation (IAA).[48] As a consequence, only the diagonal elements of the friction tensor η remain and off-diagonal elements vanish. In the LDFA model, η is a function of the electron density ρ(x, y, z) embedding the ion with position (x, y, z). In accordance with previous results,[49] we assume that the embedding density corresponds to a good approximation to the unperturbed electron density of the bare Pt(211) surface which is here obtained from a single DFT calculation. To compute the friction coefficient for the H(D) atom, we adopt the relation[44]where the parameters are a = 0.70881ℏ/a0, b = 0.554188, and c = 0.68314a0–1 and were previously fitted in ref (44) to ab initio data.[50] The Wigner–Seitz radius rs = (3/(4πρ))1/3 depends on the density ρ(x, y, z) embedding the hydron at position (x, y, z). It is convenient to solve eq in Cartesian coordinates and to use proper coordinate transformations to compute the potential and forces as functions of the six molecular coordinates presented in Figure .

Following previous studies on the reactive scattering of diatomic molecules from metal surfaces,[44,51] the effect of electron–hole pair excitation on the reaction of H2(D2) on Pt(211) can also be studied by scaling the LDFA–IAA friction coefficients. Here, we consider scaling factors of 1 (η = ηLDFA) and 2 (η = 2 × ηLDFA). Here, we investigate what happens if the friction coefficients are multiplied by a factor of 2 because the LDFA–IAA friction model is approximate, ignoring the possible effects of the electronic structure of the molecule. Friction coefficients computed with the orbital-dependent friction model tend to come out larger.[52−54] In the former case, we have performed calculations for TS = Tel = 0 and 300 K, whereas in the latter case, we only performed calculations at TS = Tel = 0 K, that is, in the absence of random forces.

QD Simulations

6D QD calculations are performed with the time-dependent wave packet method[55,56] using our in-house wave packet propagation code by solving the time-dependent Schrödinger equationHere, Ψ(Q⃗; t) is the corresponding nuclear wave function of molecular hydrogen at time t. The Hamilton operator used in eq accounts for the motion in the six molecular DOFs of H2 and readswhere ∇⃗ is the nabla operator, Ĵ(θ, ϕ) is the angular momentum operator for the hydrogen molecule, M is the molecular mass, and μ is the reduced mass of H2(D2). The initial nuclear wave function is represented as a product of a wave function describing initial translational motion and a ro-vibrational eigenfunction Φν,(r, θ, ϕ) of gaseous H2(D2) characterized by the vibrational quantum number ν, the angular momentum quantum number j, and the angular momentum projection quantum number m. Therefore, the initial wave function readswhere k⃗0 = (k0, k0, k0)T is the initial wave vector. The wave function describing initial translational motion is given byHere, the initial wave packet β(k0) is characterized by a half-width parameter σ according towith k̅ being the average momentum and Z0 is the position of the center of the initial wave packet.

The equations of motion were solved using the split-operator method.[57] The motion in X, Y, Z, and r was represented using Fourier grids. Quadratic optical potentials[58] were used to absorb the wave function at the edges of the grid in r and Z. A nondirect product finite basis representation was used to describe the rotational motion of H2.[59,60] To compute reaction probabilities, first S-matrix elements were computed for diffractive and ro-vibrationally elastic and inelastic scattering, using the scattering matrix formalism of Balint-Kurti et al.[61] These were used to compute probabilities for diffractive and ro-vibrationally elastic and inelastic scattering. The sum of these probabilities yields the reflection probability and subtracting from 1 then yields the reaction probability.

Computation of Observables

Using the quasi-classical method, we aim to model the sticking of H2(D2) on Pt(211) at conditions present in experiments we compare with by taking into account the different translational and ro-vibrational energy distributions characterizing the different molecular beams. At a nozzle temperature Tn, the probability Pbeam of finding molecular hydrogen in a specific ro-vibrational state ν, j with a velocity v + dv in the beam iswhere the flux-weighted velocity distributionis normalized by a normalization constant A and characterized by a width parameter α and the stream velocity v0. The ro-vibrational state distribution is given by

The weight w(j) accounts for the different nuclear spin configurations of ortho- and para hydrogen molecules. For H2, w(j) = 1/4 (3/4) for even (odd) j-values and for D2, w(j) = 2/3 (1/3) for even (odd) values of j. The function f(ν, j, Tn) is defined as

The appearance of a factor of 0.8 in the rotational energy distribution reflects that rotational and nozzle temperatures assume the relation Trot = 0.8Tn because of rotational cooling upon expansion of the gas in the nozzle.[62] The experimental beam parameters for the H2/D2 + Pt(211) systems are listed in Table .

Table 4

Parameters Used for the Molecular Beam Simulations of H2 and D2 on Pt(211)a

	⟨Ei⟩ [eV]	v0 [m/s]	α [m/s]	Tn [K]	n in ⟨Ei⟩ = nkBTn	⟨Ecorr⟩ [eV]
H2	0.004	626.5	55.9	293
	0.009	943.5	127.8	293
	0.013	1085.1	111.6	293
	0.014	1145.2	118.7	293
	0.025	1531.4	96.6	293
	0.035	1747.5	293.9	293
	0.043	2031.2	80.6	293
	0.132	3392.1	578.0	500	3.07	0.116
	0.181	3959.8	690.8	700	3.00	0.163
	0.169	4009.0	185.2	1300
	0.233	4442.8	862.5	900	3.00	0.210
	0.282	4838.8	1022.9	1100	2.98	0.255
	0.338	5223.2	1215.6	1300	3.02	0.302
	0.413	5617.0	1535.8	1500	3.20	0.348
	0.454	5790.7	1711.3	1700	3.10	0.395
D2	0.008	626.8	49.8	293
	0.027	1103.2	134.9	293
	0.040	1379.2	75.5	293
	0.054	1555.6	218.1	600
	0.076	1860.0	237.5	800
	0.110	2239.9	267.3	1100
	0.130	2430.7	311.9	500	3.03	0.116
	0.140	2531.1	294.5	1100
	0.234	3191.0	548.1	900	3.02	0.209
	0.276	3628.1	160.3	1300
	0.346	3814.8	772.0	1300	3.09	0.303
	0.457	4304.1	989.4	1700	3.12	0.395

The parameters were obtained from fitting the experimental TOF data to eq 6 in the Supporting Information. For pure H2 and D2 beams, we also provide n in ⟨Ei⟩ = nkBTn and the corrected average incidence energy ⟨Ecorr⟩ = 2.7kBTn.

The quasi-classical initial conditions are prepared using a Monte Carlo procedure described in ref (44) and sample directly the probability distribution Pbeam. The resulting probability Pi for dissociative adsorption, scattering, and non-dissociative trapping of an ensemble of molecules is determined by the ratiowhere N stands for the number of adsorbed, dissociated, or trapped trajectories (Nads, Ndiss, Ntrap, respectively) and N is the total number of trajectories computed for a specific energy point ⟨Ei⟩, where ⟨Ei⟩ denotes the average translational incidence energy of the molecule.

Computational Details

The time integration of eq is done in Cartesian coordinates using a time step of Δt = 2.0ℏ/Eh (≈0.0484 fs) with the stochastic Ermak–Buckholz propagator,[63] which also works accurately in the non-dissipative case. Further technical details are given in refs.[44,47] The maximal allowed propagation time for each trajectory is tf = 10 ps. In the non-dissipative case, our QCT setup usually leads to an energy conservation error of smaller than 1 meV. All trajectories start at a molecule–surface distance of 7 Å and initially sample the ensemble properties of the experimental molecular beam, that is, we model the ro-vibrational state distribution according to the nozzle temperature as well as the translational energy distribution of the incidence beam. The parameters characterizing the molecular beam are given in Table , and details about their experimental determination are given in the Supporting Information. The initial conditions used in the quasi-classical simulations are determined using the Monte-Carlo algorithm explained in ref (44).

We compute N = 10 000 trajectories per energy point and count trajectories as dissociatively adsorbed if they assume an interatomic H–H distance larger than 2.5 Å during the dynamics. Scattered trajectories are characterized by a sign change in the Z-component of the total momentum vector and have to pass a molecule–surface distance of Zsc = 7.1 Å. We call a trajectory trapped if the total propagation time of 10 ps is reached and neither dissociation nor scattering has occurred.

The dissociative chemisorption of H2(ν = 0, j = 0) on Pt(211) is quantum mechanically investigated over a translational energy range of Ei ∈ [0.05, 0.75] eV using two different wave packet propagations. The analysis line used to evaluate the scattered fraction of the wave packet was put at ZstartCAP = 6.6 Å. This is a suitable value because the PES is r-dependent only for all values Z ≥ 6.6 Å, so it allows representing the wave function on a smaller grid using N points in Z for all channels but the channel representing the initial state (called the specular state and represented on a larger grid called the specular grid, using N points). These parameters, and other parameters discussed below, are presented in Table .

Table 5

Characterization of the Two Different Wave Packet (WP) Calculations for (ν = 0, j = 0)H2 Incident Normally on Pt(211) for Translational Energies of Ei ∈ [0.05, 0.75] eVa

property	WP1	WP2	unit
WP grid parameters
range of X	[0, LX]	[0, LX]	a0
NX grid points in X	36	36
range of Y	[0, LY]	[0, LY]	a0
NY grid points in Y	12	12
range of Z	[−2.0, 19.45]	[−2.0, 17.10]	a0
NZ	144	192
ΔZ	0.15	0.10	a0
NZspec	210	220
range of r	[0.80, 9.05]	[0.80, 7.85]	a0
Nr	56	48
Δr	0.15	0.15	a0
jmax = mjmax	22	32
CAPs
ZCAP range	[12.55, 19.45]	[12.50, 16.90]	a0
ZCAP optimum	0.05	0.08	eV
specular grid
ZspecCAP start	22.75	16.10	a0
ZspecCAP end	29.35	19.90	a0
ZspecCAP optimum	0.05	0.08	eV
rCAP range	[4.10,9.05]	[4.55, 7.85]	a0
rCAP optimum	0.05	0.20	eV
propagation
Δt	2.00	2.00	ℏ/Eh
tf	3870.21	1741.60	fs
initial wave packet
energy range, Ei	[0.05, 0.25]	[0.20, 0.75]	eV
center of WP, Z0	16.45	14.30	a0

Specified are the grid parameters for the wave function and the PES, and parameters defining the complex absorbing potential in r and Z, the center position Z0 of the initial wave packet, and the corresponding translational energy range Ei covered.

The grids in Z start at Z = Zstart and share the same grid spacing. The grid in r is described in a similar way by the parameters rstart, Nr, and Δr. The numbers of grid points used in X and Y (N and N) are also provided, as are the maximum value of j and m used in the basis set (jmax and m). The optical potentials used [also called complex absorbing potentials (CAPs)] are characterized by the value of the coordinate at which they start and end, and the value of the kinetic energy for which they should show optimal absorption;[58] these values were taken differently for the regular and the specular grid in Z. The time step Δt used in the split operator propagation and the total propagation time tf are also provided. The initial wave packet is centered on Z0 and is constructed in such a way that 95% of the norm of the initial wave function is associated with kinetic energies in motion toward the surface between Emin and Emax, as also provided in Table .

Results and Discussion

Static DFT Calculations

Before we come to the dynamics calculations, we here first present general features of the interaction potential of atomic and molecular hydrogen and a Pt(211) surface. In Figure , we plot the minimum potential energy values for atomic H, assuming the optimal atom–surface distance Z over the full (1 × 1) unit cell. The resulting H–Pt(211) PES resembles the PES earlier developed by Olsen et al.[42] on the basis of DFT energy point calculations using a B88P86 functional.[29,30] For example, the most stable adsorption site for a single hydrogen atom on Pt(211) is located near the brg1 position at the step edge (see also Figure ). Additional minima are found close to the top2 and the top3 sites. In agreement with Olsen et al.,[42] we also obtain the largest diffusion barrier to be ≈0.60 eV above the global minimum in the vicinity of the brg2 site. The specific position of the global minimum for H adsorption suggests the minimum barrier for H2 dissociation to be on top of the step edge at the top1 site because the top1-to-brg1 path represents a short route for H atoms to assume their most favorable geometry on the surface. In general, the abstraction of atomic hydrogen from Pt(211) requires large amounts of energies as is known also for other H + transition-metal systems.[64,65] The value of Eads ≈ 3.7 eV computed here is, however, overestimated by ∼0.7 to 1.0 eV because we did not perform spin-polarized DFT calculations, which are not relevant to the comparison with the work of Olsen et al.,[42] to the reaction paths for H2 dissociation and to the dynamics of H2 dissociation.

Figure 3

Minimum potential energy for H on Pt(211) for geometry-optimized atom–surface distances Zopt on a (1 × 1) supercell. The energies are given relative to the most stable configuration of H on Pt(211) which is here near to the brg1 position (see Figure ). Because our DFT calculations do not include spin-polarization, the corresponding highest adsorption energy of 3.74 eV for a single H atom should to our experience be overestimated by ∼0.7 eV. The contour line spacing is 0.03 eV.

In Figure , we present different two-dimensional (2D) PES cuts along the H–H and the molecule–surface distances (r, Z) for H2 approaching Pt(211) with orientations parallel to the surface at different impact sites and azimuthal orientations as shown in the insets of the figure. As can be seen from Figure a, the dissociation of H2 proceeds indeed nonactivated directly over the top1 site, that is, over a Pt atom at the step edge. Following the color code of the figure, H2 can spontaneously dissociate after passing an early, but shallow barrier of E† = −83 meV (barrier is below the classical gas-phase minimum) in the entrance channel. The two H atoms are then accommodated exothermally on the surface. This result is in agreement with the previous work of McCormack et al.[14] where a nonactivated route to dissociation was revealed for impacts near the top1 site and with H atoms dissociating to brg1 sites. This result also matches up with the above analysis of the topology of H on Pt(211) PES that suggested the lowest barrier to be close to the top1 site. Furthermore, the associated barrierless path enables the contribution of a direct nonactivated mechanism for reaction at all incidence energies, as found experimentally[8] as well as theoretically.[14] Interestingly, already small changes of the molecular geometry lead to significant changes of the topology of the PES. For example, moving H2 from the step edge to the bottom of the step while retaining its orientation, as shown in Figure c, yields a 2D-PES that has a large activation barrier of E† = 556 meV and dissociation appears to be endothermic. Aligning now the molecular axis with the X-axis of the surface unit cell, as shown in Figure b, reduces somewhat the barrier but the PES becomes strongly repulsive for very large values of r (r > 2 Å). This suggests that not only the dissociation of H2 on Pt(211) may be accompanied by a strong angular reorientation dynamics but also that associative desorption may set in after the molecule has experienced large interatomic stretches.

Figure 4

2D potential cuts through the 6D PES for dissociative adsorption of H2 on Pt(211) along r and Z at the nine different sites on the (1 × 1) unit cell. In all cases, H2 approaches parallel to the macroscopic surface (θ = 90°). Top views of the molecular configurations are shown as insets. Contour levels are given in the energy range of [−1, 2] eV with a spacing of 0.2 eV. The zero value of the PES is set equal to the gas-phase minimum energy. Negative (positive)-valued contour lines are plotted in blue (black), and the zero-valued contour line is shown in red. Green circles indicate the position of the reaction barrier, and barrier heights E† are also shown.

The different impact sites and initial orientations of the molecule do not only affect how large the barrier toward bond cleavage is and the length of the path toward a favorable adsorption state. They also influence the way in which vibrational and translational energy plays in favor of reaction. Throughout the nine plots presented in Figure , one recognizes the typical elbow form of the PES along the r, Z coordinates. On the one hand, the curvature of the minimum energy paths in the elbows controls the vibration–translation (V–T) coupling,[66] which may facilitate dissociation in quasi-classical simulations artificially due to the zero-point energy conversion effect: the higher the curvature, the more coupling. On the other hand, the Polanyi rules[67] relate the efficiency of translational and vibrational excitation of the incident molecule for reaction to the position of the barrier. In late-barrier system resembling the product state reaction is promoted vibrationally, whereas in early-barrier systems, reaction is more enhanced by translational excitation. For the H2 + Pt(211) system, vibrationally nonadiabatic V–T processes as well as the Polanyi rules are expected to come into play during the reaction dynamics. For example, we find relatively early barriers for impact situations shown in Figure b–d, suggesting a preference of translational excitation for reaction. Impact sites associated with a late barrier are shown in Figure f,h,i. In impacts on these sites, reaction is more likely to be promoted by initial vibrational excitation.

Reaction barrier energies and associated geometries for the nine incidence situations outlined in Figure are specified in Table . Although the barriers to dissociation could be decreased somewhat when optimized with respect to θ for cases in which H2 does not dissociate parallel to the step (Figure b,d,e,g,h), Figure and Table nevertheless provide a good view of the H2–Pt(211) interaction. We find the latest (r† = 1.62 Å) and highest barrier (E† = 692 meV) for molecules incident at the brg4 site (see also Figure h). This indicates a considerable range of activation energies (∼700 meV) for the dissociation process. The Z†-values reported in Table range from 0.51 Å at the top2 site (bottom of the step) to 2.79 Å at the top1 site (top of the step edge). This resembles to some extent the overall shape of the Pt(211) surface because step-top and step-bottom Pt atoms are displaced by ΔZ = 1.27 Å.

Table 6

Barrier Heights and Geometries for H2 on Pt(211) for the Geometries Shown in Figure a

Configuration	r† [Å]	Z† [Å]	E† [eV]
top1 (ϕ = 90°), Figure 4a	0.75	2.79	–0.083
top2 (ϕ = 0°), Figure 4b	0.90	0.59	0.396
top2 (ϕ = 90°), Figure 4c	0.88	0.51	0.556
top3 (ϕ = 90°), Figure 4i	1.00	0.99	0.118
brg1 (ϕ = 0°), Figure 4d	0.80	1.75	0.186
brg3 (ϕ = 0°), Figure 4e	0.94	0.73	0.639
brg4 (ϕ = 30°), Figure 4h	1.62	0.75	0.692
brg5 (ϕ = 120°), Figure 4g	0.89	1.37	0.318
t2b1 (ϕ = 90°), Figure 4f	1.34	1.53	0.035

Energies are given relative to the gas-phase minimum energy of H2.

The vdW-DF2 functional employed here yields not only rather large activation energies for the direct dissociation process but also considerable physisorption wells of ∼72 meV located comparably far away from the surface. The presence of such wells may additionally contribute to the trapping dynamics of small molecules or may even increase the chance of redirecting the molecule toward non-dissociative pathways. Baerends and co-workers[13,14] previously reported on the importance of trapping as a mechanism for indirect dissociation of H2 on Pt(211). They used a PES that was constructed on the basis of standard GGA–DFT calculations, and the authors found only a shallow physisorption well for impacts at the bottom-step. When using the DF2-functional in the description of the dynamics of molecular hydrogen on Pt(211), as done in this work, the trapping mechanism may become more substantial, which may affect the computation of sticking probabilities for slow molecules.

Comparison of QCT and QD

Figure shows the comparison between the QCT and QD results for H2 (ν = 0, j = 0). As already discussed in the introduction, the shape of the reaction probability curve in both the QD and the QCT dynamics arises from the presence of a trapping mechanism, which yields a contribution to the reactivity that decreases with incidence energy, and an activated mechanism, the contribution of which increases with incidence energy. As a result, the reaction of the H2 molecule on Pt(211) exhibits a nonmonotonic behavior as a function of the collision energy. The reaction probability curve shows very high dissociation probabilities at very low collision energies. The minimum value of the reaction probability is at an intermediate value of the collision energy, and the slope of the reaction probability curve becomes positive at higher collision energies.

Figure 5

Initial-state resolved reaction probability for H2 (ν = 0, j = 0) dissociation on Pt(211) calculated with QD in comparison with the QCT results.

As noted by McCormack et al.,[14] with a GGA PES, nonactivated indirect dissociation may occur when a molecule hits the lower edge of the step on a nonreactive site, which showed the presence of a shallow chemisorption well on their PES. A difference with our PES is that physisorption can occur anywhere at the surface because of the presence of vdW wells for the PES computed with the vdW-DF2 correlation functional.

The QCT calculations reproduce the QD results at the higher incidence energies reasonably well. At low and intermediate energies, in the QD results, the trapping mechanism manifests itself by the occurrence of peaks in the reaction probabilities, with the peak energies corresponding to the energies of the associated metastable quantum resonance (trapped) states. The comparison suggests that at low and at intermediate energies (up to 0.2 eV), the QCT results tend to overestimate the reactivity a bit. This could be due to two reasons. First, the increase of the reaction probability with decreasing energy at the lower incidence energy is understood to occur as a result of trapping of molecules entering the potential well, in which energy from the motion perpendicular to the surface is transferred into rotation and translational motion parallel to the surface.[17] In the QD calculations, trapping should only be due to energy transfer to the motion parallel to the surface.[17] However, classically, it is also allowed that energy is transferred from the motion toward the surface to the rotational DOFs.[17] Second, the QCT calculations may suffer from an artificial effect called zero-point-energy leakage, that is, in QCT calculations, the quantization of vibrational energy may be lost and the original vibrational zero point energy may be transferred to other DOF.

Isotope Effects in QCT Results for Reaction of (ν = 0, j = 0) H2 and D2

Comparison between the computed QCT reaction probability curves for H2 and D2 shows that the reaction probability of H2 is higher than that of D2 at the same incidence energy (see Figure ). We attribute this to a zero-point energy effect. H2 has more energy in zero point vibrational motion than D2, so there is a higher probability that a given amount of this energy is transferred to motion along the reaction coordinate. Gross and Scheffler[68] for H2 dissociation on Pd(100) showed that in classical dynamics (no initial zero-point energy), there is no isotope effect between H2 and D2 in the sticking probabilities. At first sight, one might expect that steering is less effective for D2 because of its higher mass and therefore less reaction for D2 than H2. On the other hand, D2 is slower than H2 at the same kinetic energy, so there is more time for the steering force to redirect the D2 molecule to a nonactivated path. However, they found the quantum dynamical sticking probabilities of D2 to be smaller than those of H2. They suggested that this small difference should be a quantum dynamical effect and that the larger vibrational zero point energy of H2 can more effectively be used to cross the reaction barrier.

Figure 6

Initial-state-resolved reaction probabilities for the dissociation of H2(D2) on Pt(211) surface are shown with red (black) symbols for the ground rotational and vibrational state. The results are obtained with the QCT method.

No isotopic dependence and also no surface temperature dependence for the sticking probability were reported by the experimentalists,[8] as shown in Figure where we show the sticking probability as a function of average incidence energy. (In ref (8), the sticking probabilities were shown as a function of the incidence energy corresponding to the most probable energy for a density-weighted incidence energy distribution, see the Supporting Information).

Figure 7

Experimental[8] sticking probability of H2 (red symbols) and D2 (black symbols) on Pt(211) as a function of average collision energy.

Comparison of Molecular Beam Sticking Probabilities with Experiment

Parameters used for the molecular beam sticking simulations (previously extracted from experiments as discussed in the Supporting Information) of H2 and D2 on Pt(211) are given in Table .

The sticking probabilities extracted from molecular beam simulations for H2 dissociation on Pt(211) are shown in Figure with a comparison to the experimental results. In the figure, the red circles show the theoretical results obtained from simulating the experimental beam conditions. The black circles display the experimental results reported by Groot et al.[8]Figure shows the same comparison for D2 dissociation on Pt(211). In both cases, in the lower-energy regime, the theoretical results overestimate the experimental reaction probabilities. For H2 on Pt(211), at higher energies, the theoretical results also overestimate the experimental results. However, overestimation happens only at the highest incidence energy for D2 + Pt(211). The energy shift (the distance along the energy axis between experimental data points and the interpolated theoretical curve) is [7-92] meV for H2 + Pt(211) and [3-55] meV for D2 + Pt(211). On this basis, our results for H2 + Pt(211) do not yet agree with the experiment to be within chemical accuracy (≈43 meV). To find the mean deviation of the theoretically calculated sticking probability curve from the experimental results, we also calculated the mean absolute error (MAE) and mean signed error (MSE). We obtained a MAE of 40.8 meV and a MSE of 9.8 meV for H2 and a MAE of 32.4 meV and a MSE of −0.4 meV for D2. On this basis, the errors in the theoretical data in both cases are less than 1 kcal/mol ≈43 meV.

Figure 8

Sticking probability for molecular beam of H2 on Pt(211) simulated with QCT. For comparison, experimental results reported by Groot et al. (black symbols: experimental data from ref (8)) are plotted beside the theoretical results (red symbols). The arrows and accompanying numbers show the collision energy difference between the interpolated theoretical results and experimental data.

Figure 9

Sticking probability for molecular beam of D2 on Pt(211) simulated with QCT. For comparison, experimental results reported by Groot et al. (black symbols: experimental data from ref (8)) are plotted beside the theoretical results (red symbols). The arrows and accompanying numbers show the collision energy difference between the interpolated theoretical results and experimental data.

As already stated, the comparison between experimental and theoretical results is not yet good at the lower incidence energies. Two reasons might be involved, which are related to that being an important contribution to sticking from a trapping-mediated mechanism. The first reason concerns the inability of the QCT method to describe the sticking probability accurately when trapping contributes to reaction. The QCT results overestimate the contribution of trapping because of translation-to-rotation energy transfer, which is not allowed in QD descriptions[69] (see Section ). The QD calculations of Figure suggest that for reaction of H2 (ν = 0, j = 0), the reaction probability decreases faster with energy at low incidence energies if quantum effects are included, which goes in the right direction for getting better agreement with experiment. The other effect that could be important is surface temperature, which we do not include in our calculations. The initial reaction probability was experimentally determined at the surface temperature of 300 K. However, the experimentalists did not observe any surface temperature dependence.[8] In our view, this makes it unlikely that the static surface approximation we used here is responsible for the discrepancy with experiment at low incidence energy.

Especially for H2, our QCT results overestimate the experimental sticking probability at high average energies, as computed from the beam parameters available from fitting experimental TOF spectra (see the Supporting Information). One question we addressed is whether this could be due to errors arising from fitting these parameters, which is critically difficult especially at high incidence energies associated with short flight times. Now it is rather well known that for pure H2 beams, the average translational energy should not exceed 2.7kBTn, as no vibrational cooling occurs, and only about 20% rotational cooling.[62,70,71] Comparing the average incidence energies of the pure H2 beams in Table with 2.7kBTn, we however find that in most cases, the average incidence energies exceed 3kBTn, and this also holds true for pure D2 beams (see also Table ). This suggests that the experimental average incidence energies extracted from the beam parameters were too high. By re-plotting the experimental results using average incidence energies Ecorr equal to 2.7kBTn, we can redo the comparison with the computed sticking probabilities, if we assume that the computed values do not much depend on the nozzle temperature through altered ro-vibrational state distributions. This is likely to hold true for nonactivated or weakly activated dissociation. As Figure shows, this approach tremendously improves the agreement with experiment for the higher incidence energies at which the sticking is dominated by activated dissociation and for which the QCT results should be accurate (see Section ): the agreement with experiment is now within chemical accuracy for these energies and pure H2 beam conditions. For D2, the agreement is not as good as for H2 for the lower incidence energies in the high-energy range (see Figure ), which is perhaps due to the rotational cooling being somewhat more efficient for D2 than for H2, due to the lower rotational constant of D2. This means that in Figure , the experimental data could move somewhat to the right (to higher energies), thereby improving the agreement with the experiment. Note also that in principle, the fits of the beam parameters are expected to be less error prone for H2 than for D2 because of longer flight times of D2.

Figure 10

Sticking probability for molecular beam of H2 on Pt(211) simulated with QCT. For comparison, experimental results reported by Groot et al. (black symbols: experimental data from ref (8).) are plotted beside the theoretical results (red symbols). The arrows and accompanying numbers show the collision energy difference between the interpolated theoretical results and experimental data. In plotting the experimental results, we have assumed that the average incidence energy in the experiments was equal to 2.7kBTn.

Figure 11

Sticking probability for molecular beam of D2 on Pt(211) simulated with QCT. For comparison, experimental results reported by Groot et al. (black symbols: experimental data from ref (8).) are plotted beside the theoretical results (red symbols). The arrows and accompanying numbers show the collision energy difference between the interpolated theoretical results and experimental data. In plotting the experimental results, we have assumed that the average incidence energy in the experiments was equal to 2.7kBTn.

Another solution to the puzzle of why the average incidence energies calculated from the beam parameters did not correspond to 2.7kBTn for pure beams could be that the nozzle temperature was actually higher than measured. This could in principle be simulated by assuming that the nozzle temperature can be computed from the measured average incidence energy, instead of adapting the average incidence energy to the measured nozzle temperature. This was not pursued computationally, as it would only be expected to lead to a small increase of the computed sticking probability, and to somewhat larger discrepancies for H2 + Pt(211), for which the agreement with experiment was worst to start with.

Above, we have suggested that the rotational cooling in a D2 beam could be somewhat more efficient than that in the H2 beam (due to the rotational constant of D2 being lower). If this were true, this would suggest that we could have plotted the experimental data for the pure D2 beams as a function of ⟨Ei⟩ = ckBTn with c somewhat larger than 2.7 (for instance, 2.75 or 2.8) in Figure . If this would be correct, this would increase the agreement between theory and experiment in this figure, as already discussed above. However, it should also alter the conclusion regarding the absence of an isotope effect drawn originally by the experimentalists: if this assumption would be correct, the sticking probabilities measured for H2 should be somewhat higher than those for D2, at least for the results from the pure H2 and pure D2 experiments. This would bring theory and experiment in agreement also regarding the qualitative conclusion on the isotope effect.

Comparison of MD and MDEF Results for Sticking

Figures and 13 show the results of MD and MDEF calculations for H2 + Pt(211) and D2 + Pt(211). At low energies, adding electronic friction and doubling the friction coefficient increase the sticking probability for D2. Doubling the electronic friction coefficient increases the sticking probabilities of H2 only at intermediate energies. At higher incidence energies, adding electronic friction decreases the sticking probability a little bit. Adding this energy dissipation channel reduces sticking somewhat at higher incidence energies because energy in the bond stretch coordinate is nonadiabatically dissipated to electron–hole pair excitation. Also, modeling the effect of the finite electronic temperature decreases the sticking probability at lower incidence energies, but there is no dramatic effect at higher incidence energies. The effect of Tel is negligible for ⟨Ei⟩ > 0.13 eV and very small at lower incidence energy. At the lowest incidence energies, the electronic dissipative channel enhances the trapping and, therefore, the dissociation probability.[72] The dissociation process is expected to increase in the presence of a trapping mechanism because once the molecule is trapped on the surface and starts to dissipate energy, it is difficult for the trapped molecule to recover the perpendicular translational energy to escape from the surface. The effect of including electron–hole pair excitations is therefore to increase the trapping-mediated contribution to the reactivity and thereby the reactivity. However, it keeps the direct mechanism almost unchanged. Raising the electronic temperature at lower incidence energies, that is, through the presence of hot electrons, leads to collisions of the hot electrons with the molecule that can excite the molecular DOFs and provide the trapped molecule with sufficiently high energy to get desorbed from the surface to the gas phase. Taking the electronic temperature in our calculations at lower incidence energies into account diminishes the trapping effect and therefore reduces the overall reactivity.

Figure 12

Sticking probability as a function of the average incidence energy obtained from MD and MDEF calculations. Black symbols show the MD, red and purple symbols show results of MDEF calculations using friction coefficient multiplied by different factors (×1 and ×2, respectively), and green symbols show MDEF results using an electronic temperature Tel = 300 K.

Figure 13

Sticking probability as a function of the average incidence energy obtained from MD and MDEF calculations. Black symbols show the MD, red and purple symbols show results of MDEF calculations using friction coefficient multiplied by different factors (×1 and ×2, respectively), and green symbols show MDEF results using an electronic temperature Tel = 300 K.

The good agreement between the MD and MDEF results at higher incidence energies confirms that the BOSS model, which does not consider electron–hole pair excitation, may accurately describe the dissociation of H2 and D2 on Pt(211) through the direct reaction mechanism at the terrace, and therefore, at higher incidence energies.

Conclusions

To address the question whether the SRP-DF functional derived for the H2 + Pt(111) is transferable to the H2 + Pt(211) system, we have performed calculations on the dissociation of H2/D2 on the stepped Pt(211) surface. We used the VASP software package to compute the raw DFT data. The CRP interpolation method was used to accurately fit these data and construct the 6D PES based on the PBEα–vdW-DF2 functional with α set to 0.57. The potential energy for H on Pt(211) for geometry optimized atom–surface distances on a (1 × 1) supercell was discussed and was compared with the previously developed PES of Olsen et al.[42] We have also discussed features of the PES for H2 dissociation on Pt(211) and reported on minimum barrier heights and associated geometries.

We have performed calculations within the BOSS model and within the MDEF model, in order to study nonadiabatic effects on the dissociation dynamics due to the creation of electron–hole pairs in the surface. The QCT method has been used to compute the initial-state resolved reaction probability and molecular beam sticking probability. The initial-state resolved reaction probability results obtained with the QCT method were compared with the results of QD calculations. The QCT calculations reproduced the QD results at the high-energy range but not at the low-energy range. The discrepancy between the results of these two dynamics methods at the low-energy regime was discussed. We have also shown and discussed the isotope effect in the QCT results of the reaction probability of (ν = 0, j = 0) of H2 and D2.

We have computed the sticking probabilities of molecular hydrogen and deuterium on Pt(211) and compared our theoretical results with the experimental data. Our theoretical results showed that the reactivity on Pt(211) is enhanced relative to Pt(111), in agreement with experiment. The lowest barrier height for reaction was found at the upper edge of the step. Reaction on the upper edge of the step is not activated. We have simulated molecular beam sticking probabilities and compared them with the experimental data of Groot et al.[8] We have reported the energy shifts between the experimental data and the spline-interpolated theoretical data to be [7-92] meV for H2 +Pt(211) and [3-55] meV for D2 + Pt(211). Thus, chemical accuracy was not yet achieved in our theoretical results. However, it is well known that the average energy of pure H2 beams should not exceed 2.7kBTn because of the absence of vibrational cooling and the occurrence of only about 20% rotational cooling for a pure beam. Nevertheless, we found that in most cases, the average energies of the pure H2 and the pure D2 beams exceeded 3kBTn. Consequently, we have replotted the experimental results employing average energies equal to 2.7kBTn and redone the comparison with computed sticking probabilities. With this modification, the agreement between experiment and theory tremendously improved for H2. The agreement between theory and experiment for D2 was not as satisfactory as for H2 at the lower incidence energies in the high-energy range. These results suggest that the experiments should be repeated and be reported for more accurately measured beam parameters to enable a better determination of the accuracy of the theoretical results.

Finally, we have presented the comparison of MD and MDEF results for the sticking probability for both H2 and D2 and discussed the effect of adding electronic friction and doubling the friction coefficient, and the effect of electronic temperature on the sticking at low and high incidence energies.