Quantum Monte Carlo calculations on dissociative chemisorption of H₂ + Al(110): Minimum barrier heights and their comparison to DFT values

In order to correctly model molecular interactions on metal surfaces, accurate barrier heights E_b for the surface reactions are required. They are important to both the study of heterogeneously catalyzed processes¹ and the elementary surface reactions of which these processes consist.^2,3 By tuning the energy of the rate-controlling state, it may be possible to lower or increase the rate of the associated elementary reaction.⁴ In heterogeneous catalysis, the rate-controlling state⁵ is frequently a transition state of a dissociative chemisorption reaction. This is the case, for example, in steam reforming⁶ and ammonia production.^7,8 Since, worldwide, the greater part of chemical production involves heterogeneous catalysis,⁹ benchmark data for dissociative chemisorption reactions are thus of great interest.

To provide theoretical data, density functional theory (DFT) is often employed as the main computational method. In DFT, one often uses the generalized gradient approximation (GGA) and, increasingly, the meta-GGA level of theory to study the interactions of molecules with metal surfaces. The drawback of the corresponding semi-local functionals is that barrier heights computed for molecule–metal surface systems may vary strongly depending on the specific functional used,² allowing for a large variation in the predicted reactions rates of heterogeneously catalyzed reactions.

To overcome these shortcomings, a semi-empirical DFT approach called the specific reaction parameter approach to DFT (SRP-DFT) has been implemented and demonstrated to provide chemically accurate barrier heights E_b (accuracy better than 1 kcal/mol) for molecule–metal surface reactions.¹⁰ The SRP-DFT approach has predictive power to the extent that an SRP functional developed by reproducing one experiment on a particular system is required to describe at least one other experiment on the same system with quantitative accuracy and will usually also describe additional experiments on the same system with such accuracy.¹⁰ In addition, SRP density functionals developed for one specific system have been shown to be transferable to a chemically similar system in a number of instances.^11,12 A drawback of the SRP-DFT method is that the experimental data required to fit a SRP functional may not be available for all systems of interest and, where such data are available, details concerning, e.g., the velocity and internal state distributions of the reacting molecules may be lacking,¹³ and/or different experiments may yield differing results.^13,14 It is therefore also important to develop and/or to test first principles methods that are capable of high accuracy calculations on molecules interacting with metals. Even if those calculations may only be achievable for one or for a few geometries, the results may serve as benchmarks.

One such approach developed with the aim of achieving higher accuracy for molecule–metal surface interactions is embedded correlated wave function (ECW) theory. Although several different embedding schemes exist,^15–19 the general idea of ECW theory is to divide a system into a local region of interest, i.e., the cluster, which is treated with correlated wave function theory, and its environment, which is treated with DFT.¹⁷ ECW theory in its various flavors has been applied already to a few interfacial systems^18,20–23 and has been shown to be capable of providing an accurate description of systems in which charge transfer, multiconfigurational character, and excited states play an important role.²⁴ Recently, for example, a potential energy surface (PES) based on ECW calculations was used successfully²⁵ in quasi-classical trajectory calculations to semi-quantitatively reproduce sticking probabilities for a system for which DFT at the GGA level notoriously fails, i.e., O₂ + Al(111). Although this result²⁵ is promising, some drawbacks of this theory are that increasing the number of atoms in a cluster makes calculations rapidly become impractical due to the computational cost and that the error associated with incomplete convergence of the energy with the size of the embedded cluster may be several kcal/mol.²² 

Within the family of DFT itself, an interesting method that potentially allows for higher accuracy for molecule–metal interactions is the random phase approximation (RPA) or, to be more precise, adiabatic-connection fluctuation–dissipation density functional theory within the RPA.^26,27 The RPA may be considered as a rung 5 functional (the GGA and meta-GGA functionals being rung 2 and rung 3 functionals, respectively).²⁸ The RPA uses the exact exchange (or the Hartree–Fock) energy as computed from the occupied Kohn–Sham orbitals.²⁹ The RPA correlation energy is calculated separately, using both occupied and virtual orbitals. The RPA has already been employed to study the adsorption of molecules on metal surfaces,^30–32 and it has even been argued that RPA energies can serve as benchmarks for calculations on adsorption of molecules on metal surfaces.³³ Several successes have been achieved with the RPA for molecule–surface interactions, for example, in solving the CO adsorption puzzle.³⁴ However, given that the performance of the method for a database of 10 chemisorption energies with a mean unsigned error of 4.8 kcal/mol fell halfway between results obtained with the RPBE and BEEF-vdW results,³³ we argue that presently evidence is lacking that the method is accurate enough to provide benchmark results. Improvements may be achievable by using RPA including single excitations³⁵ or second order screened exchange.³⁶ Both methods, however, tend to worsen the description of barrier heights,^37,38 and more modern solutions (e.g., Ref. 39) may be needed to achieve an improvement.

Another interesting development is the implementation of coupled-cluster theory for solids,^40,41 for which certain tests exhibit good agreement with experimental findings. Although these results show promise for future development,⁴² the method currently has not yet been rigorously used to benchmark the interactions of molecules with metal surfaces.

Other promising methods come from the family of stochastic electronic structure theories called quantum Monte Carlo (QMC). In principle, the QMC method is exact.^43–45 Additionally, many QMC methods show a favorable scaling with system size (especially if localized orbitals are used⁴⁶), and large-scale parallelization is feasible.⁴⁵ Of this family of methods, Diffusion Monte Carlo (DMC) has already been applied to a few molecule–metal surface systems. In 2008, Pozzo and Alfè⁴⁷ utilized DMC to study H₂ dissociation on the Mg(0001) surface. The study itself provided useful insights into certain issues, such as the role of system size, but the accuracy of the DMC results could not be verified due to the lack of experimental data for the specific system studied. In more recent work, some of us⁴⁸ used DMC to study the benchmark H₂ + Cu(111) system. After systematic corrections were applied, the DMC result agreed within 1.6 ± 1.0 kcal/mol with a chemically accurate semi-empirical value of the reaction barrier height. QMC has also been used to study N₂ + Cu(111),⁴⁹ CO and H₂O + Cu(100),⁵⁰ H₂ + Pt(111),⁵¹ and the site-preference of CO on Rh(111), Ir(111), Pt(111), and Cu(111).⁵² 

Here, we use DMC to study the dissociative chemisorption of H₂ on Al(110). Unlike Cu, the electronic structure of Al does not have electrons in d-orbitals, making the calculations less computationally expensive. As a result of the lower cost, it is possible to study barriers associated with multiple geometries of the system and to investigate how the variation of barrier height with geometry compares with the variation obtained with several widely used density functionals. Unfortunately, directly benchmarking such results against experimental values is not possible as barrier heights are not directly accessible from experiment.² For H₂ + Al(110), however, molecular beam sticking experiments have been performed,⁵³ ultimately allowing for a validation of the results via dynamics calculations. Sticking in this system has been studied early on with theory^54,55 on the basis of the local density approximation and the PW91 density functional,⁵⁶ respectively, in both cases applying a simple dynamical model (the hole model⁵⁷). Both theoretical studies^54,55 resulted in overestimation of the measured sticking probabilities. While presently it should not yet be possible to develop a PES on the basis of DMC calculations only, it may be possible to fit a SRP density functional to the DMC results presented here. This SRP functional may be used in dynamics calculations with an appropriate dynamical model to derive sticking probabilities for comparison with the experimental results to validate the DMC results presented here. As this requires additional dynamical calculations, this will be part of a future publication.

The goals of the present paper are thus to present DMC results for the system H₂ + Al(110) as well as to compare these with the DFT results obtained with selected functionals for the same system. We focus on the energies of the two lowest lying transition states as determined with DFT and on four more reaction barrier heights computed in reduced dimensionality. Together, these energies describe how the barrier height varies across the surface and with molecular orientation. As noted above, a direct comparison with sticking experiments on H₂ + Al(110) is outside the scope of the present paper as this requires a dynamical calculation. However, an additional goal of our paper is to pave the way for such a comparison by providing DMC benchmark results that can be used to fit an SRP density functional, which, in turn, can be used to compute a PES for the system under investigation. As noted in a recent perspective paper,² requirements for such a PES are that it accurately describes not only the minimum barrier height but also the variation of the barrier height with impact site and orientation.⁵⁷ 

This paper is organized as follows: Section II briefly describes the QMC methods used. Section III details the setup for the geometries and the parameters in the DFT and QMC calculations. Section IV discusses the DMC results and some of the associated errors and compares with the DFT results. Section V concludes the paper and also provides an outlook.

In this section, we briefly describe the QMC methods used in this paper, namely, the variational Monte Carlo (VMC) method and the diffusion Monte Carlo (DMC) method.

The DMC method uses the imaginary time Schrodinger equation to evolve a set of electronic configurations from Ψ_T toward the ground state. A typical DMC calculation consists of two parts; the first is an equilibration phase in order to represent the ground-state wave function, and the second is a further propagation to accumulate statistics for properties such as the ground-state electronic energy. Due to the fermion sign problem,⁴⁴ the DMC method makes use of the fixed-node approximation⁶⁰ in which the nodes of the wave function are fixed to those of Ψ_T. The DMC method then produces the lowest-energy state possible. This value may, however, not be exact as it contains fixed-node errors. For reviews on this method, see Refs. 44 and 61.

As discussed in the Introduction, the goal of the DMC calculations is threefold: (a) To provide accurate values for the two energetically lowest lying barriers and for four additional barriers in reduced dimensionality for the dissociative chemisorption reaction of H₂ on Al(110); (b) to thereby provide benchmark data to check DFT functionals against; and (c) to pave the way for the development of an SRP density functional based on first principles (DMC) calculations, which can be used to validate the DMC results by comparing experimental sticking probabilities with the values computed with dynamics calculations using a PES based on SRP-DFT data. These three goals have mostly determined the choices that we have made in (i) how we choose the points in the PES to include in the investigation, and (ii) how the DMC and the DFT calculations are otherwise set up. While we present many details of the computational setup (for example, the exact setup of the DFT calculations and the form of the Jastrow function) in the supplementary material, we want to discuss here the main points influencing the accuracy, reliability, and significance or meaning of the values that we compute.

Most of the choices we describe in this section follow directly from the requirement that the DMC data resulting from this study should allow fitting an SRP density functional or choosing an appropriate functional as basis for a dynamics study that will allow a comparison of computed reaction probabilities with experiment. In activated dissociative chemisorption, the sticking probability is governed by the entrance channel and barrier regions of the PES.⁶² The most important points to address with DMC therefore concern reaction barrier geometries and their energies relative to that of the molecule at asymptotic distances from the surface, and these are what we focus on.

The H₂ on Al(110) system is characterized by two similarly low lying transition states: one above the long bridge site (denoted TS1 in the following) and one above the short bridge site (denoted TS2). Since the energetic difference between these two transition states is very small (around 1 kcal/mol–2 kcal/mol, depending on the functional) compared to the total barrier height (around 20 kcal/mol), it is clearly of interest to include both transition states in our investigation: both of them will likely be dynamically important in the dissociative chemisorption reaction. Details on how these transition states were determined can be found in Sec. II a of the supplementary material. Here, we only want to stress that in optimizing TS1 and TS2 (as well as TS3–TS6 later on), the metal surface was kept fixed at the geometry of the relaxed metal surface in vacuum. The reason for this choice is simply that, in highly activated dissociation studied in UHV, at high collision energies of interest, the surface atoms do not have time to respond to the motion of the incoming molecule.⁶³ The TS geometry of interest is therefore the one in which the metal surface corresponds to the relaxed metal surface in vacuum.

To ultimately allow an accurate description of the variation of the barrier height across the surface and with H₂ orientation, four more points on the PES (denoted TS3–TS6 in the following) were added to TS1 and TS2. TS3–TS6 correspond to transition states in restricted two-dimensional (2D) cuts through the PES at high symmetry impact points. The H₂ orientations relative to the surface chosen in this study are shown in Fig. 1 and listed in Table I. Details on how the exact TS geometries were obtained can be found in Sec. II b of the supplementary material. TS3–TS6 have strongly differing total energies that allow us to test the performance of DFT functionals for barrier heights over a wide range of energies.

At first sight, it may seem that, ideally, the full geometry of the TSs in the PES should be determined self-consistently, i.e., by optimizing the full geometry using the method in question (DMC or DFT with a specific functional), in all calculations. While there is something to be said for such a choice, this is (i) not possible in DMC at the present point in time as there is no periodic DMC code available yet that can compute forces at the DMC level using pseudopotentials and a plane-wave (or B-spline) basis and (ii) not desirable even in the DFT calculations for reasons discussed in the following. Let us consider the short-bridge transition state (TS2) as an example. Depending on the exchange–correlation functional used in the transition state search with DFT, this transition state will either lie parallel to the surface (as in PBE-vdW-DF^64–67) or in a slightly tilted geometry (as in PBE). Additionally, the H–H binding distance as well as the distance of the molecule to the surface will slightly change. Comparing barrier heights for self-consistently optimized transition states would allow us to compare the actual (self-consistently optimized) DFT barrier heights, but it will not allow us to compare the energies obtained with different functionals and/or with QMC for one specific molecular geometry in the 6D PES of the molecule. The latter is, however, much more desirable when trying to choose the best functional for a subsequent molecular dynamics calculation through a comparison with DMC. In the present paper, we have therefore chosen a different route, sketched in Fig. 2: we use fully self-consistent PBE-DFT transition state searches exclusively to determine geometries of the molecule relative to the surface that are likely to be dynamically relevant for dissociative chemisorption. The PBE functional was chosen because it often represents a compromise between accuracy for the metal lattice and accuracy for reaction barriers.⁶⁸ For simplicity, TS2 was restricted to a geometry in which H₂ is parallel to the surface, as suggested by the PBE-vdW-DF calculations. The difference between the PBE energies of this restricted transition state TS2 and the actual PBE short-bridge transition state was less than 0.1 kcal/mol, so this restriction should not have a large influence on the barrier height. TS3–TS6 were restricted to specific impact sites and H₂ orientations as given by Table I and do not necessarily correspond to true transition states; they can be considered as barrier geometries in a reduced 2D (r, z) space.

Having determined how to choose the molecular geometry, we now move to the geometry of the surface. It is known that different functionals predict different lattice constants and interlayer spacings. In the case of aluminum, for example, PBE predicts a lattice constant of a₀ = 4.04 Å, RPBE-vdW-DF predicts a₀ = 4.08 Å, while the extrapolated 0 K and anharmonic zero-point vibration corrected lattice constant extracted from experiment is a₀ = 4.02 Å⁶⁹ (see Sec. I c of the supplementary material). In the DFT calculations comparing to QMC, it is preferable to use the self-consistently optimized DFT lattice parameters and interlayer distances for the slab rather than experimental values for several reasons. First, lattice strain is known to have a strong influence on the barrier height for small molecules reacting on metal surfaces.^70–72 Transferring the exact same geometry (including the same slab geometry) to a calculation with a different functional may thus result in errors due to lattice strain. Second, while the static surface approximation is often used in describing sticking of H₂ on metal surfaces, the small weight of Al suggests that this is not sufficient when describing the dissociation of H₂ on Al(110) and that surface temperature effects (leading to a thermally corrugated surface) may be important. These temperature effects can, however, only be captured in dynamics calculations if one can define a 0 K slab in its relaxed geometry. To overcome both issues, in the DFT calculations, we use the self-consistently determined PBE molecular coordinates for each TS and, for each functional, a slab geometry that is optimized self-consistently for the metal surface interacting with vacuum (i.e., using the functional in question with no molecule present). For the DMC calculations, we cannot determine the surface structure self-consistently due to the absence of forces in the relevant codes. However, an obvious choice follows from the notion that DMC should, in principle, be exact. In this case, one can use the experimental geometry of the surface in the DMC calculations. The parameters that we used to describe the slab in the VMC and DMC calculations (and obviously in the DFT calculation preparing the initial wave function for the QMC calculations) are shown in Table II. These parameters are described in Sec. I c of the supplementary material.

The molecular coordinates (see Table I) used in the QMC and DFT calculations are thus extracted from self-consistent PBE-DFT calculations by determining the H–H bond distance r, the molecule-surface distance z, the impact site (top, hollow, …), and the orientation of the molecule relative to the surface described by the angles θ and φ (see Fig. 3). This molecular geometry is then used in all QMC and DFT calculations, each of which uses their own slab geometry. For the asymptotic geometry, we use a self-consistently optimized geometry for the H₂ molecule for all DFT calculations used in comparing with QMC and the experimental H₂ bond distance for the DMC calculation.

The above has consequences for what we call the TS geometries and TS energies: only the PBE energies used to compare with QMC are energies of (restricted) TSs in the true sense, i.e., these PBE TS geometries and energies have been obtained for geometries that were self-consistently computed with PBE-DFT for both the molecule and the surface. All other TS geometries and energies computed here are so in only an approximate sense, but they do represent an equivalent geometry in the 6D PES of the molecule. In fact, the molecular geometries used in the QMC calculations together with the DFT surface geometries may be thought of as defining full geometries that serve as “anchor points,” which we can use in subsequent SRP-DFT calculations to “nail” the DFT PES to the “QMC PES,” by demanding that the DFT energies are close to the QMC values for these points.

In the present study, we consider a surface coverage of $14$⁠, i.e., one H₂ molecule per 2 × 2 Al(110) surface unit cell. As shown in Sec. III c of the supplementary material, for this cell size (i.e., 2 × 2), we can expect the barrier heights to be converged to less than 1 kcal/mol. In comparing DMC with DFT data, this somewhat incomplete convergence in coverage should not have a big influence as we compare the results obtained at the same coverage.

Concerning the number of layers, we consider a metal slab with 10 layers. This is, considering other studies of molecules on metal surfaces, an astonishingly high value. However, we deemed this necessary as large oscillations were observed in the DFT barrier height when increasing the number of layers for H₂ on Al(110). These oscillations in the barrier height as a function of the number of layers,⁷³ which may be due to quantum size effects⁷⁴ (analogous to standing waves that may or may not form in a box of a certain size or analogous to Friedel oscillations that arise from surface inhomogeneities⁷⁵), slowly died out to about ±1 kcal/mol when using 10 layers or more (see Sec. III c of the supplementary material). Obviously, these oscillations have a direct impact on the barrier heights that we want to extract for TS1 and TS2. For TS3–TS6, we are mainly interested in the energetic differences between DFT and DMC calculations. As we found the oscillations in barrier height to be only weakly influenced by changes in the lattice constant, the choice of functional, and the specific adsorption site studied (see Sec. III c of the supplementary material), incomplete convergence in the number of layers should have a small effect on the comparison as all calculations are performed using the same number of layers. However, since we are interested in absolute numbers for TS1 and TS2, 10 layers were deemed to be a good choice, as the 10-layer result in DFT seems to coincide well with the results for very large amounts of layers. We expect possible errors in the barrier heights of TS1 and TS2 resulting from the restricted number of layers to be considerably less than 1 kcal/mol (see Sec. III c of the supplementary material).

As the metal slabs used in our calculations contain many Al layers, the number of layers used dominates the computational cost. To reduce this cost in the (computationally expensive) DMC calculations, we used a trick. As discussed in more detail in Sec. II c of the supplementary material, in our DMC calculations, it is computationally advantageous to use a setup in which the molecules are not only placed on the top of the slab but also added at the bottom of the slab. The molecular geometry originally extracted for on top adsorption is thereby replicated to the bottom of the slabs used in our DMC calculations as well as in the DFT calculations that are used in the comparison with DMC.

A well-motivated choice of pseudopotentials is essential in order to achieve high accuracy and to allow for a good comparison of DFT and DMC calculations. We use pseudopotentials in both the plane-wave DFT and the DMC calculations. The PAW Al_GW pseudopotential with three valence electrons treated explicitly, which is provided in the VASP software package, yields very good DFT results for the lattice constant and the bulk modulus of aluminum when compared to all electron calculations (see Sec. IV of the supplementary material). However, it cannot be used in the DMC calculations because it is fully non-local. Hence, several different semi-local pseudopotentials (all with a Ne core for Al and one electron treated explicitly in H) were tested for their accuracy, including the Al and H pseudopotentials developed by Burkatzki et al.,⁷⁷ the Al and H pseudopotentials developed by Trail and Needs,^78,79 and the Al and H Trouiller–Martins-type fhi98PP pseudopotentials provided on the Abinit website.^80,81 The performance of the various pseudopotentials was tested by comparing the DFT-PBE bulk lattice constant, the bulk modulus, and the TS1 and TS2 barrier heights to all electron PBE calculations and by comparing the AlH binding energy as computed in DMC to all-electron results using coupled-cluster singles, doubles and perturbative triples [CCSD(T)] with the AVQZ basis set. The details of these tests are given in Sec. IV b of the supplementary material. Our tests suggest that the Trail–Needs pseudopotentials give a very good performance for the problem at hand (lattice constant correct to within 0.01 Å, bulk modulus correct to within 1%, DFT barrier heights for TS1 and TS2 as well as Al–H binding energies accurate to within 0.2 kcal/mol). Similar observations hold true for the PAW GW pseudopotentials distributed with VASP and discussed in this paper (see Sec. IV of the supplementary material for further specification) but not for some of the other pseudopotential combinations. All DMC calculations and the DFT calculations for the trial wavefunction generation for VMC were hence performed using the Trail–Needs pseudopotentials. All other DFT calculations were performed using the PAW GW pseudopotentials distributed with VASP as the use of Trail–Needs pseudopotentials is highly inefficient in plane-wave DFT calculations due to the high plane-wave cutoff required.

Details on the setup of the DFT calculations used in optimizing the slab geometry, in determining the various transition states, and in performing the final DFT calculations used in the comparison with DMC data are given in Secs. I a, I b, and V of the supplementary material. The choices made are motivated based on convergence tests (see Secs. III a and III b of the supplementary material).

Details on the setup of the QMC calculations are given in Sec. VI of the supplementary material. Briefly summarized, QMC calculations were performed using the CASINO^43,82 package v.2.13.709 (beta), and the trial wave functions were generated using Slater determinant wave functions from PBE-DFT calculations performed with the QUANTUM ESPRESSO⁸³ package. VMC was used to optimize a Jastrow factor (containing up to three-body terms) for the subsequent production runs of the DMC calculations. The optimization was performed by minimizing the energy expectation value.⁸⁴ The time step as well as the number of walkers used in the DMC calculation was determined from convergence tests as detailed in Sec. IV. Pseudopotentials were treated using the T-move scheme.⁸⁵ 

Single-particle finite-size effects were minimized via twist averaging. Twist averaging is a procedure in which several different k-points are taken into account, achieving a result somewhat similar to k-point averaging in DFT. We used 32 symmetry unique twists in the 2 × 2 supercell and 8 unique twists in the 4 × 4 supercell resulting from a regular 8 × 8 (4 × 4 for the 4 × 4 supercell) k-point grid. Details on the twist angles used can be found in Sec. VI a of the supplementary material. Many-body finite-size effects were minimized by extrapolating the results obtained for a supercell containing 2 × 2 surface atoms and 4 × 4 surface atoms to infinite system size. To keep a coverage of $14$⁠, we thereby consider two H₂ molecules in the 2 × 2 case (one on top and one on the bottom of the slab) and eight H₂ molecules for the 4 × 4 slab (4 on top and 4 on bottom). Mathematical details on the twist averaging procedure as well as on the finite-size extrapolation can be found in Sec. VII of the supplementary material. The results presented below suggest that this procedure allows us to extract reliable energy differences (i.e., molecule–surface interaction energies) for the system at hand.

Since we are interested in comparatively high accuracy (<1 kcal/mol), we present briefly our convergence tests for the number of walkers and the time step used in the DMC imaginary-time propagation. These convergence tests were performed for the TS1 barrier height at the Γ-point only, using a 2 × 2 supercell. Due to the large number of twists used in the calculations and the associated requirements in equilibrating the walkers and in obtaining reliable error bars, we would like to limit the number of walkers to a comparatively small number. The results shown in Fig. 4 suggest that using a time step of 0.06 a.u. and 2500 walkers leads to a negligible bias in the energy, and these settings were used in all subsequent calculations.

In the following, we present the raw DMC data taking TS1 as an example to exemplify the procedure of twist averaging and finite-size extrapolation. Together with the trends observed for the other TS (data provided in Sec. VIII of the supplementary material), this will allow us to demonstrate why we are confident that the finite-size errors are under control, i.e., they are within a fraction of a kcal/mol.

We start with a discussion of the twist averaging procedure. Particularly, for metallic systems, twist averaging is highly important, very much in the same way as converging the k-point grid in DFT calculations is. Using TS1, we demonstrate why and to what extent we are confident that our twist averaging procedure leads to sufficiently converged results. Figure 5 shows the raw DMC data together with the corresponding DFT data for each twist considered in the 2 × 2 supercell. Although the scatter in the data for different twists is very large (more than 10 kcal/mol), the twist-averaged DFT result is very close to the k-point converged DFT result (see Tables III and S25–S29 in the supplementary material). The largest difference was found for TS4 with a 1.2 kcal/mol shift (see Table S27 of the supplementary material). This shows that the amount of twists considered is reasonable. Additionally, a clear linear trend between the DMC and the DFT data at different twists is visible. Although the differences between the individual DMC data and the fitted line can be up to a few kcal/mol (and are thus clearly larger than the statistical error for each individual twist), this allows us to fit the DMC vs DFT energy via a linear regression. The resulting slope of the fit curve is m = 1.2(1) for TS1 (see Table III). The uncertainty in m is computed from the covariance matrix of the linear fit. Comparing this value with the results for TS2–TS6, the slope m in TS1 actually shows an untypically large deviation from m = 1, where m = 1 may be viewed as “ideal” behavior.

For the calculations using the (2 × 2) cell, resulting single-particle finite-size corrections $Dsp−fs$⁠, which should correct for the difference between twist-averaged and k-point converged result, [see Eqs. (S1)–(S3) in the supplementary material] range from −1.2 kcal/mol (for TS4) to 0.4 kcal/mol (for TS1). Keeping in mind the targeted accuracy of well below 1 kcal/mol, the size of these corrections may be considered large, especially when considering the non-perfect linear fit. Fortunately, errors that remain in the energies obtained for the 2 × 2 supercell after correction for single-particle finite-size are not projected 1:1 to the final finite-size corrected result. Instead, an offset of 1 kcal/mol in the 2 × 2 cell will lead to comparatively small offsets of around 0.2 kcal/mol in the fully final finite-size corrected results due to the scaling relation in Eq. (S4). Possible uncertainties in the 2 × 2 supercell are thus not at all worrisome, and we will therefore focus on the results from the 4 × 4 cell in the following.

The results for the 4 × 4 supercell are shown in Fig. 5(b), again taking TS1 as an example. The spread in the DMC (and DFT) values for different twists is strongly reduced compared to the 2 × 2 supercell. For TS1, the spread is around 4 kcal/mol; for TS2–TS6, the spread is often even smaller, around 2 kcal/mol (see Figs. S6–S10 in the supplementary material). While the spread in values is strongly decreased, this is not true for the deviation of the DFT twist-averaged and DFT k-point converged energy difference, $ΔEk−pointconv.DFT−ΔĒtwistsDFT$⁠, which still lies in the range of 0.3 kcal/mol–0.5 kcal/mol for TS1–TS6. This may be expected since the number of twists considered is smaller in the 4 × 4 supercell. What is surprising, though, is the fact that this difference is often negative in the 2 × 2 cell, while we observe it to be positive in the 4 × 4 cell (Tables III and S25–S29 in the supplementary material). Since the twists used in the 2 × 2 cell correspond to those used in the 4 × 4 cell after k-point unfolding, one might have expected this shift to be similar in the 2 × 2 and the 4 × 4 cell. The reason for this offset has to lie in the fact that the twist averaging is performed canonically (i.e., constant number of electrons at each twist), and this is, in our view, a strong sign that supercells larger than 2 × 2 should be considered (going to cells even larger than 4 × 4 is unfortunately not possible with present day computational power).

The discussion above clearly shows that if sub kcal/mol accuracy is sought, single-particle finite-size effects should be corrected for, even with the 4 × 4 supercell. In the 4 × 4 cells, the DMC data typically follow the linear trend within one or two standard deviations σ of the raw DMC data. Unfortunately, the relatively large error bars (∼0.4 kcal/mol) in the results of each twist, together with the small range of values available, make estimating m difficult. Indeed, for the 4 × 4 supercell, we often observe m values much smaller than 1, all associated with comparatively large error bars [e.g., m = 0.3(4) for TS4]. While these error bars are considered in the error propagation, we found it interesting to investigate the limiting case of setting m = 1. Doing so typically only lead to small shifts in the final results for $EbDMC$⁠. The maximum deviations were observed for TS3 and TS4, where differences between $EbDMC$ when using m = 1 to perform the single-particle finite-size correction differed from those obtained using the fitted value for m by 0.30(30) kcal/mol and 0.36(30) kcal/mol, respectively. Although these shifts are considerably smaller than 2σ, one may want to keep them in mind as these shifts seem to be systematically positive (only exception: TS5). Nevertheless, even when considering these shifts as possible systematic errors, the resulting twist averaging procedure should be accurate on a sub kcal/mol scale and the resulting barriers rather under than overestimated. We also note that the values of m for the lowest two transition states, i.e., TS1 and TS2, are close to one. The remaining analysis is performed with m taken from the fitting procedure. Values obtained for m = 1 are given in Sec. IX of the supplementary material.

The last correction to discuss is the many-body finite-size correction. As shown in Table IV, the many-body finite-size correction resulting from the extrapolation to infinite supercell size, $Dmb−fs=EbDMC$ − $ΔEsp−fsDMC(4×4)$⁠, is small for TS1 [see Eqs. (S4)–(S6) in the supplementary material for mathematical details]. This holds true for all other TS as well, for which corrections of no more than 0.2 kcal/mol were observed. Interestingly, the extrapolations to infinite cell size show a positive slope for all TSs except TS5 (bottom graph of Fig. S9 in supplementary material), which has a negative slope. This slope would actually be positive if the analysis were done using values for $ΔĒsp−fsDMC$ obtained for m = 1. In any case, the small absolute values of $Dmb−fs$ are reassuring as small errors in the scaling relation given in Eq. (S4) in the supplementary material should hence not matter. This suggests that the use of 4 × 4 supercells (together with extrapolation techniques) is sufficient in order to obtain accurate data.

The resulting DMC values for $EbDMC$ (including all finite-size corrections) are reported in Table V for TS1—TS6. Note that the error bars stated include uncertainties in the single-particle finite-size corrections stemming from uncertainties in the slope m [Figs. 5(a) and 5(b)] as well as uncertainties in the many-body finite-size errors arising from the finite-size extrapolation [Fig. 5(c), see also Sec. VII of the supplementary material for details on the error propagation and error analysis].

Interesting to note are the barrier heights TS1 and TS2. DMC predicts these values to be 25.1(2) kcal/mol and 26.7(2) kcal/mol, respectively. As discussed further in Sec. III e of the supplementary material, a crude error analysis suggests that the systematic errors in these values should be lower than 1 kcal/mol when considering the influence of (i) the restricted number of layers used, (ii) the high coverage used, (iii) the restricted vacuum distance and the restricted distance of the molecule from the surface in the asymptotic geometry, (iv) the accuracy in the TS geometries extracted using PBE, (v) the transferability of pseudopotentials, (vi) time step error, and (vii) assuming the fixed-node error in the reaction barrier height to be +1 kcal/mol on the basis of an analysis of DMC results for gas phase reactions (see below). Additionally, one may want to consider the influence of time step bias and the finite-size effects in DMC discussed above and locality errors. These errors are also estimated to lie well below 1 kcal/mol, as discussed in Secs. IV a and IV b, and in Sec. III e of the supplementary material, respectively, and the same should be true for the total statistical error in the final DMC energies.

DFT predicts TS1 and TS2 to lie very close in energy, with TS1 generally lying just below TS2. This behavior is confirmed in the DMC calculations. The energetic difference between TS1 and TS2 is 1.6(3) kcal/mol on the DMC level and is somewhat larger than that predicted by DFT when using the PBE functional (0.4 kcal/mol). This variation in energy will be an important parameter when judging the quality of various exchange–correlation functionals in the following.

In the following, we compare the DMC data obtained in Sec. IV C with the DFT results obtained for a range of commonly used exchange–correlation functionals. This will allow us to assess their performance and will enable the selection of candidate functionals for the later development of an SRP functional. To obtain the DFT data that we present here, the slabs are optimized self-consistently for each functional, as explained above. In the supplementary material, we present similar DFT results, but these were obtained using experimental slab geometries throughout. Qualitatively, the results can be considered very similar.

The main results are summarized in Fig. 6(a) and Table VI; a more detailed view for TS1 and TS2 is shown in Fig. 7.

From Fig. 6(a) as well as from Table VI, it is immediately striking that almost all functionals considered (with the RPBE-vdW-DF functional being an exception) underestimate the DMC barriers [see Fig. 6(a)]. This might seem surprising insofar as PBE and RPBE (both significantly underestimating the DMC barrier height here) often bracket the true barrier height extracted from SRP-DFT calculations fitted to molecular beam experiments on dissociative chemisorption on a metal surface.^10,86 However, examples where DFT calculations with the RPBE exchange–correlation functional apparently underestimate barrier heights are known, such as O₂ + Al(111)^22,87 and HCl + Au(111).^88,89 In addition, a similar result as obtained here has been found in an earlier study on H₂ + Mg(0001) when comparing RPBE with the DMC results: the DMC barrier (27.2 ± 0.7 kcal/mol) computed for this reaction was higher than the RPBE value by 2.5 kcal/mol.⁴⁷ A very recent study found evidence that the RPBE functional underestimates barrier heights for dissociative chemisorption on metals if (W-EA) is smaller than ∼7 eV, where W is the work function of the metal and EA is the electron affinity of the molecule.⁹⁰ As a result of a low work function, charge transfer to the molecule may occur in the transition state, and this would lead to underestimation of the TS energy when using a density functional containing semi-local exchange.⁹¹ For H₂ + Al(110), (W-EA) ≈ 6.8 eV–6.9 eV,⁹⁰ so it is not surprising that the DMC barrier height exceeds the RPBE value. In summary, the comparison of the QMC with the DFT barrier heights is in line with what one might expect.

Nevertheless, it seems worthwhile to ask the question whether errors in the DMC calculations could cause a systematic offset. Studies of gas phase reaction barrier heights suggest that there could be a systematic error in DMC results due to fixed-node errors. Since DMC is variational when using the T-move scheme, the DMC barrier height would be overestimated if the fixed-node error is larger for the transition state than for the asymptotic geometry. Zhou and Wang indeed found DMC to overestimate barrier heights on average by 0.9 kcal/mol for a database of 38 hydrogen atom transfer reactions.⁹² Similarly, Krongchon and co-workers⁹³ found DMC based on a PBE wave functions to overestimate gas phase reaction barrier heights for a dataset of 38 barrier heights for non-hydrogen-transfer reactions by about 1 kcal/mol. Analyzing their results, Wagner et al.⁹³ found the size of the fixed-node error to correlate with the difference in energy between the LUMO and the HOMO (the lowest unoccupied and the highest occupied molecular orbitals, respectively). They also found the fixed-node error to be larger in the transition states than in the product and reaction states, explaining why the barrier heights (which are energy differences of variational calculations) are overestimated. If the latter finding also holds for molecule–metal surface reactions, fixed-node error would be expected to contribute about +1 kcal/mol to the overall systematic error in the barrier heights we computed for H₂ + Al(110). These conclusions are tentative as drawing conclusions for molecule–surface reactions from gas phase reactions may be dangerous. Additionally, since the fixed-node error seems to be sensitive to bond-breaking events,^92,93 one might also expect the fixed-node error to vary somewhat for the six geometries considered as the corresponding barriers occur at considerably different H–H bond stretching distances and molecule–surface heights. It is thus clear that future research of the size of fixed-node errors in barriers for dissociative chemisorption on metal surfaces should be of high interest.

Apart from the fixed-node error, possible deficiencies in the finite-size corrections are a possible source of error. However, as discussed in the analysis of finite-size corrections, we expect that applying these corrections leads to underestimated DMC barrier heights rather than overestimated values due to the uncertainty in the slope m of the fit of DMC vs DFT values for each twist. Based on Refs. 92 and 93 and on the discussion above, we thus believe that the systematic underestimation of the barrier heights by most DFT functionals when compared to DMC values is in large part not an artifact, but real, and speculate that this finding may be related to the low work function of Al surfaces. In the following, we will thus consider the DMC values as reference values.

Using the DMC values as reference values, we can compute different statistical observables [the root-mean square deviation (RMSD), the mean signed deviation (MSD), and the mean absolute deviation (MAD)] that are indicative of the discrepancies between DFT and DMC energies of TS1–TS6 for the functionals considered. The result of this analysis is summarized in Table VI. Clearly, the error measures depend strongly on the functional, as is already clear from Fig. 6, where an approximate downward trend in the barrier heights is visible for all TS considered (note that we sorted the functionals in Fig. 6 to show this trend, while trying to keep “sister”-functionals, such as PBE and PBE-vdW-DF, next to each other).

The poor performance in terms of RMSD of the strongly constrained SCAN⁹⁴ functional for H₂ + Al(110) (RMSD = 7.0 kcal/mol) is in line with its poor performance on adsorption and chemisorption of molecules on metal surfaces^95,96 and on the dissociative chemisorption of H₂ on Cu(111).⁹⁷ However, in our particular case, PBE, which is one of the workhorses in surface studies, does even worse (RMSD = 7.3 kcal/mol). Barrier heights calculated with RPBE,⁹⁸ the often used meta-GGA functional revTPSS,⁶⁸ and the recently derived meta-GGA functional MS-B86bl⁹⁷ show a clear improvement relative to PBE⁶⁴ and SCAN. However, with an RMSD of 2.5 kcal/mol, 3.9 kcal/mol, and 2.8 kcal/mol, respectively, these functionals are still far from being chemically accurate for the DMC barrier heights. In particular, revTPSS still severely underestimates the dynamically relevant lowest lying reaction barrier heights, TS1 and TS2 (see Table VI). According to the RMSD, the BEEF-vdW⁹⁹ barrier heights (RMDS = 0.9 kcal/mol) are consistently better than those obtained with the other functionals. The good performance of the BEEF-vdW functional for H₂ + Al(110) is in line with other theoretical investigations for surface processes and is probably due to the fact that BEEF-vdW has been semi-empirically fitted to both gas phase reaction barrier heights and adsorption energies of molecules on transition metal surfaces. The good performance of BEEF-vdW is, however, quite closely matched by RPBE-vdW-DF^65,98 (RMSD = 1.4 kcal/mol), albeit with the disadvantage that the energetic corrugation seems to be described worse in RPBE-vdW-DF (see below). The finding that BEEF-vdW and the RPBE-vdW-DF functional yield barrier heights in quite good agreement with the DMC values is encouraging: it suggests that an SRP-DF for H₂ + Al(110) may be found that is based on GGA exchange and a van der Waals correlation functional.

Figure 6(a) shows a distinct trend for the functionals tested for all barrier heights: functionals that strongly underestimate TS1 also tend to underestimate TS2–TS6 very strongly. This indicates that the energetic corrugation (i.e., the variation of the barrier height with impact sites—tested in TS1–TS4) and the anisotropy of the barrier height (i.e., how the barrier height varies with orientation—tested by the comparison of TS1 with TS5 and of TS2 with TS6) may be captured quite well by all functionals tested. This is further supported by Figs. 6(b) and 6(c), where we show measures indicative of the energetic corrugation and anisotropy. As shown in Fig. 6(b), all functionals seem to describe the change in barrier height from TS1 to the other TSs quite well [note the different scale on the y axis between Figs. 6(a) and 6(b)]. As a matter of fact, nearly all errors observed in the energetic corrugation and anisotropy are below 2 kcal/mol [see Fig. 6(c)]. The results also indicate that, interestingly, the offset does not scale with barrier height. This is an important observation since many physical observables, such as the width of the sticking probability vs collision energy curve, depend sensitively on the energetic corrugation and the anisotropy of the barrier height,¹⁰⁰ as suggested also by application of the hole model.⁵⁷ This effect (i.e., that barrier heights are shifted by nearly constant values when using different functionals) has been observed earlier [e.g., for H₂ + Cu(111)¹⁰⁰]. Based on Bayesian statistics, Ref. 101 suggested this to be a general feature of GGA functionals for the energetic corrugation of molecular adsorption. The essential contribution of the present work is that this effect can also be found when comparing with first principles theory, i.e., with DMC.

Interestingly, functionals that do better in terms of RMSD for the barrier heights do not necessarily perform better for the RMSD of the energetic corrugation (see, for example, Table VI for RPBE-vdW-DF) and the other way round. BEEF-vdW, however, seems to perform well with respect to both measures, i.e., the RMSD in the barrier heights and the RMSD in the energetic corrugation. The main drawback of this functional for the system at hand thus seems to be the fact that it underestimates TS1 and TS2 by 0.8 kcal/mol and 1.7 kcal/mol, respectively, which is more than would be desirable when quantitatively modeling molecular beam experiments. However, none of the other functionals tested performed any better.

Although the energetic corrugation is overall well described by all functionals considered, Fig. 7 points to another interesting fact: except for revTPSS, which breaks the trend, all functionals considered seem to lie more or less on one line when plotting the difference between the DFT and DMC barrier height for TS1 vs that for TS2. This is strongly in line with the previous observation that all functionals seem to capture the energetic corrugation. A striking observation, however, is that all functionals seem to underestimate TS2 by ∼1 kcal/mol more than they underestimate TS1. The statistical uncertainty in the DMC results for TS1 and TS2 is too small to explain such a large offset. We therefore analyze the possible reasons for this effect:

• TS2 might have a larger fixed-node error than TS1, thus increasing the barrier height of TS2 compared to that of TS1. Considering that TS2 lies further away from the surface in a geometry more reminiscent of the asymptotic geometry, this seems unlikely although further tests would be necessary to rule out this option, as suggested in Sec. V.

• Either TS1 or TS2 could be influenced more strongly by the lattice constant (which is considerably smaller in experiment than predicted by DFT). A similar (and even stronger) trend is, however, observed when comparing the DFT results obtained with the experimental slab geometry (see Sec. X of the supplementary material). Hence, this seems an unlikely explanation.

• There could be a structural deficiency present in most exchange–correlation functionals used in the present work. Indeed, the vdW functionals typically follow the energetic corrugation for TS1 vs TS2 better than the non-vdW functionals (compare PBE vs PBE-vdW-DF and RPBE vs RPBE-vdW-DF in Fig. 6). Additionally, BEEF-vdW, which is based on DF2, which typically shows a stronger van der Waals well than functionals based on DF1, does particularly well with respect to the energetic corrugation between TS1 and TS2. This is therefore a plausible explanation.

We have presented DMC barrier height calculations for six single points corresponding to barrier geometries in the PES of H₂ on Al(110).

Four of the single-points correspond to transition states in restricted space for the dissociative chemisorption of H₂ on Al(110), while two of the single-points correspond to a good degree of accuracy to “real” transition states for the system. The corresponding barrier heights of the latter two were found to be 25.1(2) kcal/mol and 26.7(2) kcal/mol. Systematic errors of these barrier heights obtained with a simplified analysis were somewhat smaller than 1 kcal/mol and mainly arise from incomplete convergence of the energies with the number of layers and the supercell size and from the value of the fixed-node error that was estimated on the basis of results for gas phase reactions.

Using the DMC values as reference values, the performance of several DFT functionals was investigated. It was found that most functionals underestimate the reaction barrier heights, but they correctly capture the energetic corrugation. None of the six functionals for which results are presented here yield barrier heights agreeing with the DMC results to within chemical accuracy for all six barriers investigated. BEEF-vdW performed well with a maximum deviation of 1.7 kcal/mol, closely followed by RPBE-vdW-DF with a maximum deviation of 2.1 kcal/mol.

In the future, more functionals will be tested in order to derive an SRP functional. This functional should ideally show chemical accuracy for the energetically most relevant barriers, TS1 and TS2, and preferably also describe TS3–TS6 with high accuracy. Finding such a functional will allow for a PES to be generated such that dynamics calculations can be performed. This allows for sticking probabilities to be computed and compared to experiment. Such an analysis will give insight into the extent to which DMC at the present state of the art can be useful as benchmark to model the detailed dynamics of dissociative chemisorption reactions on simple metals. Specifically, if the fit of the DFT results to the DMC barrier energies is good with the functional selected, such a comparison with experiment will also allow one to assess the accuracy of DMC in describing molecule–metal surface interactions. In this case, the SRP-DFT-QMC procedure described may also provide a handle on the size of the fixed-node error as well as other systematic errors in the calculations.

From a more fundamental point of view, several open questions remain in this paper that can hopefully be addressed in the near future. First and foremost, we did not yet attempt to assess the fixed-node error in the DMC results. Reasons for this choice include computational expense and our hope to get a better handle on the actual accuracy of the DMC calculations from a comparison with experimental molecular beam data, as discussed above. Nevertheless, such an assessment would be of interest, particularly as the fixed-node error is known to often affect the description of strong interactions^103–105 such as bond-breaking and bond-making. The fixed-node error could thus influence, for example, the observed energy difference between TS1 and TS2, as discussed in the paper. A possible approach to assess this error (at least qualitatively) would be to use various trial wave functions resulting from different density functionals. Here, we suggest the use of a meta-GGA functional, such as MS-B86bl,⁹⁷ which should be capable of simultaneously describing the metal lattice as well as the molecule-metal surface interaction accurately. Possible methods to actually reduce the fixed-node error would be to apply backflow transformations¹⁰⁶ or to try to reoptimize the orbitals^84,107 at the VMC stage, which, however, would likely necessitate a switch to an atom centered basis set, and this may come with problems of its own. Second, obtaining more DMC datapoints would be of interest. For example, it would be interesting to add geometries in which H₂ is dissociatively adsorbed to the surface, i.e., in the product state, thereby studying adsorption. Another example would be to address the van der Waals well of H₂ on Al(110). The relevance thereof was discussed in Sec. IV D. Third, it would be interesting to compare the DMC results to the results of other correlated electronic structure methods, some of which are under active development, such as periodic CCSD(T), but also full CI QMC. Finally, it would be interesting to assess the accuracy of the transition state geometries used in our analysis. DMC forces are already available for molecular systems,¹⁰⁸ but so far, to our knowledge, an implementation for periodic systems with pseudo-potentials is not yet available. Access to forces would not only allow determining whether TS2 has a flat or a tilted geometry but it would also allow us to obtain more accurate values for the energies and coordinates of TS1 and TS2.

The following information can be found in the supplementary material: I. Computational details concerning the generation of the slab geometries. a. Bulk lattice constants for DFT calculations. b. Interlayer spacing for DFT calculations. c. Setup of the metal slab for the DMC calculations according to experimental data. II. Details concerning the determination of single-point geometries considered. a. Transition state searches for TS1 and TS2. b. Restricted transition state searches for TS3 to TS6. c. Notes on why we replicate our molecule to top and bottom. d. The asymptotic geometry. III. Convergence tests: plane-wave cutoff, k-points, and barrier heights of TS1 and TS2. a. plane-wave cutoff. b. k-point grid. c. coverage, #layers, vacuum distance. d. Investigation of possible influences of the transfer of molecular geometry on the barrier heights of TS1 and TS2. e. Overall expected systematic errors in TS1 and TS2. IV. Pseudopotential tests. a. Bulk lattice constant, bulk modulus, and TS1 and TS2 barrier heights. b. Binding energy of AlH. V. Setup of the DFT calculations performed to compare to the DMC calculations. VI. Setup of the QMC calculations. a. Trial wavefunction generation: The slater part. b. Trial wavefunction generation: The Jastrow function and the VMC optimization. VII. Finite-size corrections applied to the DMC data. VIII. Raw DMC data for TS2–TS6. IX. Finite-size extrapolated values, when using m = 1 in the correction of single-particle finite-size effects. X. Comparison with DFT results using the experimental geometries.

This work was financially supported through a NWO/CW TOP Grant (No. 715.007.001). Furthermore, this work was carried out on the Dutch National e-Infrastructure with the support of NWO-EW.