Quantum Monte Carlo Calculations on a Benchmark Molecule-Metal Surface Reaction: H2 + Cu(111).

Accurate modeling of heterogeneous catalysis requires the availability of highly accurate potential energy surfaces. Within density functional theory, these can-unfortunately-depend heavily on the exchange-correlation functional. High-level ab initio calculations, on the other hand, are challenging due to the system size and the metallic character of the metal slab. Here, we present a quantum Monte Carlo (QMC) study for the benchmark system H2 + Cu(111), focusing on the dissociative chemisorption barrier height. These computationally extremely challenging ab initio calculations agree to within 1.6 ± 1.0 kcal/mol with a chemically accurate semiempirical value. Remaining errors, such as time-step errors and locality errors, are analyzed in detail in order to assess the reliability of the results. The benchmark studies presented here are at the cutting edge of what is computationally feasible at the present time. Illustrating not only the achievable accuracy but also the challenges arising within QMC in such a calculation, our study presents a clear picture of where we stand at the moment and which approaches might allow for even more accurate results in the future.

Introduction

Accurate calculations of reaction barrier heights, Eb, for elementary reactions of molecules with metal surfaces are crucial in allowing for predictive modeling of elementary surface reactions[1−3] and heterogeneous catalysis.[4−7] In spite of the importance of heterogeneous catalysis, a first-principles method capable of predicting Eb with chemical accuracy (errors ≤1 kcal/mol) has not yet been demonstrated. Currently, most calculations targeting such systems use density functional theory (DFT) at the generalized gradient and meta-generalized gradient approximation levels (GGA and meta-GGA, respectively) to compute electronic energies. As a result of the inaccuracy of present day GGA and meta-GGA functionals,[5] reaction probabilities computed for elementary surface reactions with different functionals may show major discrepancies,[3,8] resulting in order(s) of magnitude differences in the reaction rates.[4]

At present, the only DFT approach that has been demonstrated to provide chemically accurate values of Eb for reactions of molecules with metal surfaces, is a novel implementation[9] of the specific reaction parameter approach to DFT (SRP-DFT).[3] This approach has yielded accurate Eb values for the dissociative chemisorption of H2 on Cu(111),[3] Cu(100),[10] and Pt(111)[11] and of CHD3 on Ni(111).[12] Nevertheless, the current state of affairs is unsatisfactory for at least two reasons. First, SRP-DFT is semiempirical: an adjustable parameter in the density functional is fitted such that supersonic molecular beam experiments on the system of interest are reproduced.[3] Second, DFT energies cannot be compared to molecular beam experiments directly. Instead, intermediate dynamics simulations are necessary. This makes the fitting procedure indirect and can introduce uncertainties due to (the simplified description of) phonon and electron–hole pair excitations in the surface and the (lack of) quantum-classical correspondences. Although calorimetric measurements and temperature-programmed desorption allow nowadays a good determination of reaction energies[13] (final state minus initial state), this experimental data can still only provide very limited information on the potential landscape and would likely be insufficient to devise an SRP-DFT functional that accurately describes intermediate barrier heights. It is thus desirable to find an ab initio method that provides chemically accurate values for (selected) points on the potential landscape, including the reaction barrier height Eb.

Modern embedding theories, in which a (small) cluster is described at a high level of accuracy and the environment is taken into account via an embedding potential, constitute such an alternative.[14−16] These methods have provided valuable insight into several interesting problems in surface science, but, with these methods, it is still hard to converge the energy with respect to cluster size. Furthermore, due to the (limited) basis set in the correlated wave function calculation of the cluster, they can suffer from significant basis set errors. Also, their accuracy has not yet been benchmarked for molecule–metal surface reactions by rigorous comparison of molecular beam experiments with dynamics results on the basis of potential surfaces obtained with this method. A second alternative is presented by modern stochastic electronic structure theories, whose emergence has opened the possibility to tackle larger and larger systems and thus also solids and surfaces directly.[17,18] Several different types of these so-called quantum Monte Carlo (QMC) methods have been applied to problems in solid-state physics in recent years, including AFQMC[19] and i-FCIQMC.[20] Nevertheless, one of the most established QMC methods for solids remains diffusion Monte Carlo (DMC). The favorable scaling of DMC with system size,[21,22] the fact that the computational effort involved in a DMC calculation does not scale with the size of the basis set used to express the trial wave function if it is recast in a set of B-splines,[23] the possibility to reduce single-particle finite-size effects by twist averaging[24] and the quality of results previously achieved, make DMC a highly promising candidate for providing accurate ab initio interaction energies for molecule–metal surface interactions. It is our goal here to establish its accuracy for H2 dissociating on Cu(111).

DMC has already been applied to molecule–metal surface reactions in the past: Pozzo and Alfè studied the dissociation of H2 on Mg(0001) in 2008.[25] While this study gives much interesting insight, for example into finite-size effects, there are currently no experimental data available that allow benchmarking the accuracy of these calculations. Additionally, many industrially relevant catalysts are based on transition metals rather than on alkali or earth-alkali metals. It is therefore our goal to tackle the computationally much more challenging problem of H2 dissociating on Cu(111) for which accurate semiempirical benchmarks are also available.

Very preliminary DMC calculations for H2 dissociating on Cu(111) have been published on arXiv by one of us.[26] However, this earlier work showed severe shortcomings in the methodology and the pseudopotentials used and the good agreement with benchmark data may be due to fortuitous error cancellation.

It is the goal of this paper to present a benchmark result for the reaction barrier height of H2 dissociatively chemisorbing on Cu(111) through state-of-the-art DMC calculations. We will address in detail the setup of the geometric structure, the use of pseudopotentials and the errors resulting thereof, residual corrections to the QMC results and finite-size effects. Finally, we will discuss what accuracy we can reach, how reliable the QMC results really are and which approaches could be used in the future to improve further on our results. The paper is organized as follows: In section , we give a brief introduction to diffusion Monte Carlo. In section , the computational setup and schemes for corrections of systematic errors are explained. Thereafter, we present our results and the discussion in section , and the conclusion and outlook in section .

Diffusion Monte Carlo in a Nutshell

Diffusion Monte Carlo (DMC) is a projector method that takes advantage of the imaginary time Schrödinger equation to project the electronic ground-state from an initial trial wave function. Excellent reviews on this method have been written elsewhere.[17,18] We will therefore only summarize the method very briefly.

In DMC, the imaginary time Schrödinger equation is recast into a drift-diffusion and branching equation of a set of particle configurations (“walkers”) in imaginary time via a stochastic implementation. After the average walker population has been equilibrated to represent the ground-state wave function, the walkers can be further propagated in imaginary time to accumulate statistics and to determine properties such as the electronic ground-state energy.

In principle, DMC is an exact method. For electronic structures, to avoid an exponential scaling of its computational cost with system size, however, the fixed-node approximation has to be applied, in which the wave function nodes are fixed to the nodes of a trial wave function. The trial wave function typically has the form of a Slater–Jastrow function, i.e., a Slater wave function (from CAS, DFT, ...) multiplied by a Jastrow factor that can introduce additional correlation into the wave function. The Jastrow factor is usually parametrized in terms of electron–electron correlations, electron–nucleus correlations, and electron–electron–nucleus three-body correlations and is optimized in a preceding variational Monte Carlo (VMC) calculation.

The propagation of walkers in DMC is not constrained in real space by any basis set limitations. Errors resulting from the use of finite-basis sets and basis-set superposition errors can only enter indirectly via the trial wave function due to the fixed-node constraint and the locality approximation.[27,28] Such basis-set related errors can be kept negligibly small by using highly converged plane wave calculations for the trial wave function generation without significant additional cost.[23] Furthermore, when using the fixed-node approximation, DMC scales as + , where c is a small constant and N is the number of particles, if a constant statistical error bar in the total energy is sought and localized basis functions are used to express the trial wave function.[21,22] Both of these properties make DMC a particularly interesting method for studying metallic systems.

Computational Setup

The Geometry

Since our goal is to predict the true barrier height Eb of the dissociation reaction of H2 on Cu(111) as accurately as possible, it is crucial to describe the true barrier geometry as exactly as possible. In DFT calculations, the transition-state geometries are typically obtained by optimizing the geometry to correspond to a stationary point in the potential energy landscape. Although geometry optimizations and the determination of minimum energy pathways are nowadays possible within variational Monte Carlo,[29−31] optimizing geometries for calculations of this size is currently too costly at the QMC level. We therefore need to rely on experimental or on DFT geometries.

The Copper Slab

DFT calculations with GGA functionals that are suitable for calculating adsorption energies of simple diatomics on transition metals generally overestimate the lattice constant by a few percent.[32] This is related to a fundamental problem that exists for DFT with GGA functionals: within the GGA, no functional has been found so far that accurately describes both the adsorption of a molecule at a metal surface and metallic lattice constants.[33,34]

Since QMC calculations are likely to be sensitive to an expansion of the surface,[35−37] they should not be based on structures obtained from DFT employing GGA functionals that are (reasonably) good for adsorption energies: If an expanded GGA-DFT lattice geometry were used, QMC would be expected to underestimate the H2 + Cu(111) reaction barrier height.[35,36] Therefore, for both the transition-state and the asymptotic geometries of the H2 + Cu(111) reaction, the experimental geometry of the Cu(111) surface is adopted, with an experimental room temperature bulk lattice constant[38]a3D = 3.61 Å. For the (111) surface, this yields a surface lattice constant of a2D = a3D/√2 = 2.553 Å. The bulk interlayer distance can be computed as d = a3D/√3 = 2.084 Å. At the surface, the interlayer distance is decreased: Medium Energy Ion Scattering (MEIS) experiments[39] show that the room temperature distances d1,2 and d2,3 between the first and second layers and between second and third layers are reduced by 1.0 ± 0.4% and 0.2 ± 0.4%, respectively. We therefore take d1,2 = 2.063 Å and d2,3 = 2.080 Å, as shown in Figure .

Figure 1

Geometry used for the barrier height calculation of H2 dissociating on Cu(111). Left: asymptotic geometry; right: transition-state geometry. Note that the DMC calculations are performed using 2D periodicity (i.e., the super cell is not repeated in the direction normal to the surface).

Based on convergence tests, the Cu(111) surface is modeled with a slab consisting of four layers. Residual systematic errors resulting from this simplification are estimated based on DFT calculations and included in a correction scheme as described below.

As vacuum distance dv (i.e., as distance between the upper layer of copper atoms and the lowest layer of copper atoms in the next periodic image) we use dv = 13.0 Å, which allows for converged DFT results. Errors that may result from this limited value are also discussed below.

The Transition-State Geometry

Having defined the copper surface geometry, the coordinates of H2 relative to the Cu(111) surface still need to be defined for the transition state and the asymptotic state. Although specific observables from reactive scattering experiments are sensitive to certain aspects of the reaction barrier geometry,[3,40] reaction barrier geometries can usually not be determined unambiguously from experiments and have to be determined via electronic structure calculations.[41] We base the reaction barrier geometry of H2 relative to the surface on previous SRP-DFT calculations.[3]

Within the DFT approach, the lowest barrier to dissociation corresponds to the bridge-to-hollow dissociation geometry, with H2 parallel to the surface, its molecular center of mass positioned above a bridge site and the dissociating H-atoms moving to the nearby hcp and fcc hollow sites (see also Figure , top views). This geometry is also physically plausible as it represents the shortest route of the H-atoms to the energetically most favorable 3-fold hollow sites.[42−44] In the SRP-DFT barrier geometry,[3] the H–H distance at the barrier is given by rb = 1.032 Å and the molecule–surface distance is given by Zb = 1.164 Å as shown in Figure . We believe the SRP-DFT estimate for rb to be accurate: The value of rb is intimately connected to the dependence of the effective barrier height E0(ν) on the vibrational quantum number ν. This dependence can be extracted from associative desorption experiments[45,46] and is well reproduced by theoretical calculations based on potential energies from the SRP-DFT approach for both D2 and H2 on Cu(111),[40,47] thus establishing the reliability of rb. To estimate the influence of possible inaccuracies in the values of Zb (and rb), we performed DFT calculations with the SRP48 functional,[48] which was designed to reproduce the SRP-DFT minimum barrier height. Simultaneously varying r and Z according to rb – δ < r < rb + δ, and Zb – δ < Z < Zb + δ, with δ = 0.05 Å, we obtained barrier heights which differed from the SRP-DFT value by less than 0.47 kcal/mol, clearly demonstrating the robustness of our calculations toward possible inaccuracies in rb and Zb. The limited vacuum distance will have no effect on the DMC calculations of the transition-state geometry, since the DMC calculations are performed using 2D periodicity only (i.e., the slab is not periodically repeated in the direction orthogonal to the surface).

The Asymptotic Geometry

For the asymptotic geometry, the value of the H–H distance is taken equal to the experimental value[49] in H2: ra = 0.741 Å (see Figure a). Due to the limited vacuum distance dv, the H2–surface distance is limited to Za = dv/2 = 6.5 Å. Larger vacuum distances are not necessary to converge the DFT results and would lead to wave function files too large to deal with for our purpose. This limited value of Za may not be sufficiently large to allow the true interaction energy between the molecule and the surface to become negligible (as intended in the asymptotic geometry). Since SRP-DFT and PBE do not account for van der Waals interaction, this limited value will not introduce an error in these DFT cases. In DMC, however, van der Waals interactions are included and the residual interaction may introduce a systematic error. Fortunately, as discussed below, the associated error is both small and fairly well known, so that it can be corrected for (see section ).

An alternative route would have been to split the determination of the electronic energy corresponding to the asymptotic geometry into two calculations: one for the molecule and one for the slab. It has been shown, however, that error cancellation in DMC will be better if the system size is not changed within one calculation.[50] We therefore prefer the former approach that leads to readily assessable errors.

The Coverage

Instead of modeling an isolated H2 molecule on a surface, we use a finite H2 coverage of 1/4 of a monolayer (i.e., in DFT, we use a 2×2 repetition of the primitive Cu(111) cell covered by one H2 molecule). This is necessary for computational reasons: In DFT, it will reduce the number of electrons and of plane waves, thus also limiting the size of the wave function file to fit into the memory of standard computing nodes (note that we need a fairly high plane wave cutoff for our pseudopotentials). In QMC, there is no saving in terms of number of electrons compared to a coverage of 1/16 (one molecule on a 4×4 primitive cell), since we need to use k-point unfolding to a 4×4 unit cell to avoid extensive finite-size errors in the QMC calculations. (The larger supercell obtained via k-point unfolding will reduce both many-body finite-size effects and single-particle finite-size effects.) The higher coverage nevertheless allows for a considerable saving in computational cost: For each QMC calculation in the 4×4 supercell, there are four molecules interacting with the metal surface. To obtain the same statistical error bar per molecule, the total error bar ΔE can thus be 4 times larger, reducing the computational cost, which scales as C ∝ 1/ΔE2, by a factor of roughly 16. Errors incurred due to the high coverage are corrected for via DFT calculations (see section ).

The Pseudopotentials

For highly accurate results, the choice of pseudopotential (PPs) is very important. While the use of Ar-core PPs for Cu is often sufficient in DFT calculations, this will in general not be true for high-level quantum chemistry methods such as DMC. Furthermore, since DMC can suffer from locality errors—especially when heavy atoms are present[51−53]—special care has to be taken in the choice of PPs applied and in the way the Jastrow function is optimized. In principle, it would be desirable to perform the calculations using Ne-core PPs for copper throughout the slab. For a 4×4 unit cell with four layers, however, this would involve the description of more than a thousand electrons. To avoid this, we use Ne-core PPs in the first layer and Ar-core PPs in the second to fourth layer. This approach and the specific choice of PPs is scrutinized in the following.

A traditional choice for a small-core PP for Cu that can be used in QMC would be the Mg-core PP by Trail and Needs.[54,55] For this PP, however, the CuH binding energy has been shown in previous work[51] to deviate from the experimental result by more than 6 kcal/mol. More recently, Burkatzki, Filippi, and Dolg have developed Ne-core PPs for 3d transition metals for QMC calculations with Gaussian basis-set based trial wave functions.[56] Their PPs, however, cannot easily be used in plane wave codes since they would require too high energy cutoffs. In 2016, Krogel et al. have developed a new database of comparably soft Ne-core PPs for 3d metals in QMC.[57] Their Cu PP suffers, however, from comparably large locality errors if the local channel is left at l = 1, as suggested for their PP (a problem that can be alleviated if changing the local channel in the QMC calculations). Furthermore, these PPs were developed based on LDA, which would be inconsistent with our GGA based calculations. Using the Opium code and a similar approach as that described by Krogel et al.,[57] we have therefore developed a new Ne-core Rappe–Rabe–Kaxiras–Joannopoulos (RRKJ) PP.[58] This PP has angular momentum channels s, p, andd, and uses lloc = 0 as local channel, which showed the smallest errors when transforming to the fully nonlocal (Kleinman-Bylander) representation. The cutoffs are set to 0.8 au for the 3s, 3p, and 3d orbitals and 2 au for the 4s and 4p orbitals. In DFT, this PP gives excellent bulk parameters for Cu and shows excellent agreement for the barrier height of the dissociation of H2 on Cu(111) with all electron results (see Tables and 2). For the DMC calculations, the local channel was set to lloc = 2, to avoid unnecessary errors arising from the angular integration. Applying this PP, we find good agreement with the experimental binding energy of CuH and acceptable locality errors (see Figure ). (Note that the case of CuH dissociation can be viewed as “worst case” scenario in terms of locality errors since the wave function at the Cu core is expected to change much more drastically in this reaction than from asymptotic geometry to transition-state geometry in the present case.)

Table 1

DFT Tests for the Small-Core Cu PP: Bulk Parametersa

	lattice constant (Å)	bulk modulus (GPa)	cohesive energy (kcal/mol)
all-electron PBE[59]	3.629	143	—
all-electron PBE[32]	3.632	—	—
experiment[32,38]	3.596	144	80
small-core Cu PP, PBE	3.628	144	82

The experimental value for the lattice constant is corrected for zero-point anharmonic expansion. For the PP tests, the PBE functional[60] was used. Lattice constant and bulk modulus are obtained from a fit to a Birch–Murnaghan equation of state. Due to the influence of the approximate xc-functional on the lattice constant, lattice constants obtained in this work should be compared with all-electron PBE results.

Table 2

DFT Tests for the Small-Core Cu PP: Barrier Height for the Dissociation Energy of H2 on Cu(111) Calculated Using the PBE Exchange Correlation Functional[60]a

	barrier height (kcal/mol)
PAW12hpv	10.8
all-electron	10.3
small-core Cu PP	10.4

PAW12hpv = projector augmented wave (PAW) potentials[61] released by VASP[62] in 2012, with a hard H PAW and a Cu PAW with the semi-core p-electrons treated as valence electrons.

Figure 2

DMC tests on the CuH binding energy using the small-core Cu PP, dependent on the number of parameters used in the parametrization of the Jastrow function: blue, using T-moves;[28,63] green, using the locality approximation;[27] red line, experimental value.[64]

As large-core Cu PP, we use an in-house Troullier–Martins Ar-core PP created using the fhi98PP code.[65] This PP has comparably small cutoffs in the pseudoization radius especially for the s-channel (1.67, 2.29, 2.08, and 2.08 au for the s, p, d, and f channels, respectively; lloc = 3) and has previously been shown to exhibit somewhat smaller locality errors in QMC than the standard PP from the Fritz-Haber Institute using the Troullier–Martins scheme.[66] Since we only use this PP from the second layer onward, the accurate reproduction of the CuH binding energy is considerably less important than for the Ne-core PP used in the first layer. The H-atom is also represented by an in-house Troullier–Martins PP[67] (s channel only; cutoff radius = 0.6 au).

Having chosen a set of PPs that can be expected to give reliable results, we also ascertain the reliability of the mixing of small- and large-core Cu PPs by performing several tests. First, to ensure that the mixing of PPs does not lead to excessive charge transfer between Cu-atoms described by different PPs, we performed DFT calculations of the charge distribution in Cu2. Using a Bader charge analysis we found only minimal charge transfer (see Table ). Additionally, we computed the DFT barrier height for the dissociation energy of H2 on Cu(111) using our Ne-core PP only and using the mixed-core PP approach and found negligible changes of 0.02 kcal/mol. Due to the positive results of these test calculations and since the largest changes in density are confined to the uppermost Cu layer (see Figure ), we expect our mixed-core PP approach to yield negligible errors.

Table 3

Bader Charge Analysis of the Charges in Cu2, When Describing One Cu-atom by an Ar-Core PP and One by a Ne-Core PP

	charge (e)
	large-core Cu	small-core Cu
nominal	11	19
obtained	10.96	19.04

Figure 3

Change in electronic density between the asymptotic geometry and the transition-state geometry based on DFT calculations using the small-core Cu PP. The inset on the left shows the position of the 2D cuts shown on the right: dark blue atoms, first layer; medium blue, second layer; light blue, third layer.

Setup of the DFT Calculations

DFT calculations are used to provide the Slater wave functions subsequently used in QMC. These calculations are performed with the pwscf code from the quantum espresso suite[68] (version 5.1 with minor modifications that enable the correct CASINO-type output for wave functions using different pseudopotentials for the same atom type and that account for an error in the original implementation of the plane wave to CASINO converter present in version 5.1 of the quantum espresso suite). All calculations use the PBE exchange-correlation functional[60] and a high plane wave cutoff of 350 Ry to ensure convergence for the very hard PPs used in this calculation. To facilitate convergence, we use a Marzari–Vanderbilt smearing with a smearing value of 0.0074 Ry. For k-point converged DFT results, a 16×16×1 Γ-centered k-point grid is used in the 2×2 supercell. To obtain the Slater part of the DMC trial wave function at each twist (see QMC setup in section ), separate, self-consistent calculations are performed in DFT: For the QMC calculations performed on the 2×2 supercell, we use a self-consistent calculation at the k-point corresponding to the twist in question. For the QMC calculations on the 4×4 supercell, we employ a 2×2×1 k-point grid in the DFT calculations that is shifted by the k-value of the twist and then use k-point unfolding to the 4×4 supercell.

Setup of the QMC Calculations

The subsequent DMC calculations are performed with the CASINO[69] software package (version 2016-04-28 (beta)) with minor modifications to correctly include two different pseudopotentials for the same atom type.

To obtain a trial wave function for DMC, the Slater wave function for the Γ-point transition-state geometry is multiplied by a Jastrow function. This function is chosen to contain electron–electron terms u, electron–nucleus terms χ, and electron–electron–nucleus three-body terms f. These terms are parametrized by polynomials of degree N multiplied by a smooth cutoff-function. For the electron–electron term, we use a polynomial of degree N = 5 and one function each for same-spin and opposite-spin electron pairs. In the electron–nucleus terms, we use Nχ = 5, again with spin-dependent functions. In the three-body terms, the expansion in electron–electron distances Nee and electron–nucleus distances Nen is of degree 2 with no spin dependence. These parameters, as well as the cutoff radii, are initially optimized in a VMC calculation at the transition-state geometry (at the Γ-point) by minimizing the variance of the local energies[70] for samples of 22 000 configurations. The cutoff radii are then fixed, and the remaining preoptimized parameters are optimized for the transition-state geometry and the asymptotic geometry separately (Γ-point only), this time with respect to energy.[71] The sample size for these optimization runs was 50 000 configurations, and we used four converged iteration cycles with changing configurations to allow the configurations to adapt to the optimized wave function. Optimizing with respect to energy has proven beneficial in minimizing locality errors in the DMC calculations.[51] Furthermore, using the same preoptimized parameters for both the transition-state geometry and the asymptotic geometry will ensure the best possible error cancellation of locality errors when taking energy differences. The influence of residual locality errors will be discussed later.

To minimize single-particle finite-size effect in the DMC calculations, we use twist averaging:[24] Since we use periodic boundary conditions to emulate the macroscopic properties of the slab, the Hamiltonian is symmetric with respect to the translation of any electron by a multiple of the supercell’s in-plane lattice vectors a. In effective single-particle theories, this leads to Bloch’s theorem, and the k-point dependence of the single-particle energies is taken into account by integrating over a k-point grid. In many-particle theories such as DMC, this symmetry translates into a many-body generalization of Bloch’s theorem. The wave function is periodic with respect to a translation of an electron along a:andwhere the twist K can be assumed to lie within the supercell’s Brillouin zone and u is periodic under a supercell lattice translation.[72] In analogy to k-point averaging in DFT, in DMC, one can average over results that are solutions to trial wave functions with different twist-angles in order to approach the thermodynamic limit faster (i.e., with smaller supercell sizes). We generate these trial wave functions by multiplying the Slater wave functions corresponding to each twist by the optimized Jastrow function. For the DMC calculations in the 2×2 supercell, we use 12 randomly chosen twists to average over and to extract the single-particle finite-size corrections. For the DMC calculations on the 4×4 supercell, we use a slightly different approach and use twists corresponding to a 4×4×1 Γ-centered k-point grid. This grid results in seven symmetry-independent twists with weights ranging between 1 and 4. The weights are taken into account in the twist averaging procedure (see eq , below), while the residual single-particle finite-size effects are computed in the same way as for the 2×2 supercell (see below). All DMC calculations use a time step of 0.005 au and a configuration size of more than 6000 walkers. Convergence tests and the possible influence of residual finite-time-step errors and single-particle finite-size effects will be discussed later. The PPs are treated using the T-move scheme[28,63] that has been proven to be beneficial in terms of the size of locality errors.[51]

Correction of Systematic Errors

The final QMC energy barrier EbDMC is given bywhere ΔE̅DMC is the twist-averaged energy difference between transition-state and asymptotic geometry. Dsp-fs and Dmb-fs are corrections accounting for single-particle finite-size effects and finite-size effects arising from long-range correlations, often referred to as many-body finite-size effects, respectively.

Adding a correction for systematic errors, Dsys, gives a corrected DMC barrier height,All quantities used above are detailed in the following.

The Twist-Averaged Energy

The twist-averaged energy difference is given bywhere M is the number of symmetry-independent twists (i.e., wave vectors[73,74]), w is the weight of the ith twist (w = 1 for the stochastically chosen k-points of the 2×2 supercell and w is the number of symmetry-equivalent k-points for the 4×4 k-point grid used for the 4×4 supercell). The DMC results for the transition-state geometry and the asymptotic geometry at twist i are given by EDMC,ts and EDMC,asy, respectively.

Single-Particle Finite-Size Corrections

Single-particle finite-size effects are very efficiently reduced by averaging over several twists,[24] as described in eq . For metallic systems, however, the number of twists needed to reach convergence is larger than what can be reasonably afforded computationally for the system under consideration: the more twists used, the higher the relative influence of the equilibration phase in the computational cost. Additionally, a minimum number of DMC steps is necessary at each twist to allow an accurate determination of error bars. Therefore, instead of using more and more twists, we correct for the remaining difference via the relationwhere ΔE̅DFT is the k-point converged DFT barrier height and ΔE̅twistsDFT is the DFT barrier height obtained in the same way as ΔE̅DMC in eq , simply replacing QMC results by DFT results. The scaling factor α results from a linear regression model of the relation of (EDMC,ts – EDMC,asy) with respect to (EDFT,ts – EDFT,asy) .

Many-Body Finite-Size Effects

Explicitly correlated methods such as DMC suffer from a second type of finite-size error, Dmb-fs. This type of error results from the contribution of long-range interactions to the kinetic energy and the interaction energy. Such an error is not present in DFT, where these contributions are implicitly taken care of by the exchange-correlation functional. Several correction schemes exist for these errors (see, e.g., ref (75) for an overview), but many of them are not strictly speaking applicable to 2D periodic systems.[75,76] A simpler and maybe more robust way to correct for these errors is by increasing stepwise the size of the supercell and extrapolating to infinite system size. For 2D periodic systems, the dependence of Dmb-fs on the system size N is suggested to scale as[75]The residual finite-size correction can thus be obtained fromwhere c and EbDMC follow from linear regression. We use the results of ΔE̅sp-fsDMC for the 2×2 supercell and the 4×4 supercell to extrapolate to infinite system size using the above scaling law. Since we are only using two system sizes, the linear regression suggested in eq thus breaks down towhere (N) is proportional to the number of particles in the system.

Systematic Effects

During the description of the computational setup, we have mentioned three possible sources of systematic errors:The total systematic error in our calculations is given by the sum of all these contributions:

errors due to the limited number of layers used in the calculation, dl,

errors due to the finite coverage, dc, and

errors due to the limited distance from the H2 molecule to the surface in the asymptotic geometry, dasy.

Correction for the Finite Number of Layers

We estimate the systematic error resulting from using only four layers by testing the convergence of the DFT reaction barrier height of H2 + Cu(111) with the number of layers nl. To allow for the highest possible accuracy in these calculations, the computational setup differs from the DFT setup used to generate the Slater part of the wave functions used in DMC (see Supporting Information). For 11 ≤ nl ≤ 16, we find a slightly oscillatory behavior with energies fluctuating by about 0.1 kcal/mol. Based on the results, we estimate the finite layer correction asSee the Supporting Information for more details.

Correction for the Finite Coverage

The coverage correction is estimated from DFT calculations on the reaction barrier height of H2 on Cu(111), with coverage values ranging from 1/4 to 1/25 of a monolayer and using a large number of layers. The van der Waals interaction is taken into account via the optPBE-vdW-DF functional. This ensures that not only directly surface mediated coverage effects are taken into account but also the possible influence of substrate-mediated van der Waals interactions[77] and residual molecule–molecule interactions. From these data, we estimate the coverage correction to beSee the Supporting Information for details on the computational setup and the results.

Correction for the Limited Molecule–Surface Distance

As mentioned above, the limited molecule–surface distance of 6.5 Å may be insufficient to allow for the van der Waals interaction between the molecule and the surface to be negligible. Any correction resulting thereof can be expected to be negative (i.e., the barrier height is overestimated) since the van der Waals interaction will spuriously lower the energy at the asymptotic geometry. Estimating the size of the residual interaction using van der Waals corrections in DFT, however, is difficult. In ref (78), Lee et al. found van der Waals well depths of 0.9, 1.2, and 2 kcal/mol, depending on whether vdW-DF2, vdW-DF, or DFT-D3(PBE) was used. Extrapolating experimental resonance data, Anderson and Peterson found a well depth of 0.7 kcal/mol. Given these numbers, we can expect the van der Waals correction at 6.5 Å from the surface (i.e., at a distance that is expected to be much larger than the location of the minimum of the well) to be rather small. Due to the strong dependence of the DFT results on the exact van der Waals functional, we prefer to use experimentally motivated values for the correction. We thus use a fit of a (theoretically motivated) function for van der Waals potential to the resonance data[78] to evaluate the van der Waals interaction at 6.5 Å:

The total systematic correction, given by the sum of eqs –14, is thereforeFortunately, this value is small.

Results and Discussion

With the setup described above, the calculation involves the description of 64 Cu-atoms and thus the description of 840 correlated electrons. As discussed in detail above, the present setup has been chosen with great care, with the aim to obtain the best possible DMC description of the problem achievable at the moment. Even with a method scaling as well as DMC, however, these calculations are computationally extremely expensive, requiring a total computational time of more than 5 million core hours on Cartesius, the Dutch National supercomputer.[79] On the one hand, this seems to be a very high cost. On the other hand, the results presented here are thus also at the forefront of what is computationally achievable at the moment with a pure DMC approach. The results and the following discussion are therefore valuable in order to see where we stand at the moment, to assess the accuracy we can expect from such a calculation at the present time and to investigate the issues that are still open and how we may want to address them in the future.

DMC Results for the 2×2 Supercell

We start our discussion with the small supercell. Ultimately, the results of the 2×2 cell will “only” be used for the finite-size extrapolation. Figure shows the DMC barrier height obtained in the 2×2 supercell for different twists in comparison with the corresponding DFT barrier (shifted by the k-point converged DFT solution). Averaging over these results and correcting for residual single-particle finite-size correction Dsp-fs (see eqs and 3a), we obtain the results stated in Table with ΔE̅sp-fsDMC(2×2) = 14.8 ± 0.7 kcal/mol. The scaling factor α entering eq is thereby determined from a linear regression to the data obtained for the 2×2 supercell, as illustrated in Figure . For the 2×2 supercell, the error in ΔE̅sp-fsDMC is dominated by the statistical uncertainty in α resulting from the goodness (or badness) of the linear fit between DFT and DMC values for the individual twists.

Figure 4

DMC vs DFT barrier height for the set of randomly chosen k-points for the 2×2 supercell. The linear regression curve is given by y = αx + β, with α = 1.1 ± 0.3 kcal/mol and β = 14.4 ± 1.8 kcal/mol.

Table 4

DMC Results for the 2×2 Supercell

ΔE̅DMC(2×2)	12.4 ± 0.4 kcal/mol
α2×2	1.1 ± 0.3
ΔE̅k-point conv.DFT – ΔE̅twistsDFT(2×2)	2.2 ± 0.0 kcal/mol
α2×2(ΔE̅k-point conv.DFT – ΔE̅twistsDFT(2×2))	2.4 ± 0.6 kcal/mol
ΔEsp-fsDMC(2×2)	14.8 ± 0.7 kcal/mol

DMC Results for the 4×4 Supercell

The results for the 4×4 supercell at different twist are shown in Figure . Averaging and applying the single-particle finite-size correction Dsp-fs, we obtain ΔE̅sp-fsDMC(4×4) = 13.3(8)kcal/mol (see Table and Figure ). For the larger supercell, both the relative error of α and the DFT correction (ΔE̅k-point conv.DFT – ΔE̅twistsDFT) are smaller than in the 2×2 case. The resulting uncertainty in the single-particle finite-size error is also much smaller in the 4×4 supercell than it was in the 2×2 supercell and, for the 4×4 cell, the statistical uncertainty in ΔE̅DMC dominates the total uncertainty of the single-particle finite-size corrected barrier height ΔE̅sp-fsDMC.

Figure 5

DMC vs DFT barrier height for the 4×4 supercell. Dots represent the seven symmetry-inequivalent k-points in the 2×2×1 Γ-centered k-point grid used for twist averaging. The linear regression curve is given by y = αx + β, with α = 2.7 ± 0.3 kcal/mol and β = 12.9 ± 1.2 kcal/mol.

Table 5

DMC Results for the 4×4 Supercell

ΔE̅DMC(4×4)	10.8 ± 0.7 kcal/mol
α4×4	2.7 ± 0.3
ΔE̅k-point conv.DFT – ΔE̅twistsDFT(4×4)	0.9 ± 0.0 kcal/mol
α4×4(ΔE̅k-point conv.DFT – ΔE̅twistsDFT(4×4))	2.5 ± 0.3 kcal/mol
ΔEsp-fsDMC(4×4)	13.3 ± 0.8 kcal/mol

Figure 6

DMC results with various corrections compared to the semiempirical reference value obtained from SRP-DFT in ref (3). Gray shaded region indicates ±1 kcal/mol from the reference value. For reference, the all-electron PBE barrier height is also shown.

DMC Results Including Finite-Size and Systematic Corrections

Taken together, the data from the 2×2 supercell and the 4×4 supercell can be used to extrapolate to infinite system size ΔE̅sp-fsDMC(N→∞) = EbDMC via eq :as shown in Table and Figure . We note parenthetically that changing the exponential scaling from N–5/4, as suggested in ref (75), to N–3/2, as suggested in ref (80), changes the finite-size-extrapolated result by only 0.1 kcal/mol, clearly demonstrating the robustness of the result with respect to the exact choice of scaling factor.

Table 6

Many-Body Finite-Size Correction

cell	N	(1/N)−5/4	value (kcal/mol)
ΔE̅sp-fsDMC(2×2)	1/4	5.657	14.8 ± 0.7
E̅sp-fsDMC(4×4)	1	1.000	13.3 ± 0.8
EbDMC	∞	0	13.0 ± 1.0

Adding the corrections for systematic errors Dsys given in eq leads to the corrected DMC barrier height (see also Figure ):

Comparison with Chemically Accurate Data

We can now compare the resulting value for EbDMC and Eb,corrDMC with the available benchmark data from a potential energy surface using an SRP-DFT functional. This functional was obtained by requiring that dynamics calculations using this SRP-DFT functional in the construction of the potential energy surface reproduce measured sticking probabilities[3] for H2 on Cu(111) as well as other observables for this system. The resulting benchmark value is given by 14.5 kcal/mol (see supporting information of ref (3).). Our DMC values for EbDMC and Eb,corrDMC differ by 1.5 and 1.6 kcal/mol, respectively, from this reference value. In view of the targeted error margin of 1 kcal/mol this result seems discouraging at first sight, but it should be kept in mind that QMC inherently comes with error bars and that the benchmark value of 14.5 kcal/mol lies within 1.5 times the error bar obtained for EbDMC (i.e., the benchmark value lies within an 87% confidence interval). Furthermore, the DMC calculations are clearly in better agreement with the SRP reference value than DFT calculations using the PBE functional (PBE also being the functional used in the trial wave function generation for the DMC calculations). This positive result, however, should not stop us from asking critical questions on the expected accuracy of the DMC results presented.

The setup and possible errors resulting from it have been assessed in detail in section . We therefore now turn to DMC related errors. While DMC is formally exact, several important approximations are made. First, the nodes of the DMC wave function are fixed to the nodes obtained in the corresponding DFT calculation for a particular twist. Second, as mentioned before, the nonlocal energy contributions resulting from the use of pseudopotentials are treated approximately, giving rise to locality errors. Third, the propagation in imaginary time is done in discrete steps, leading to time-step errors and the population size is finite, leading to population errors. While care has been taken in the computational setup to keep these errors small, they may still influence the results on the order of a few kcal/mol. Assessing these errors is already challenging in small systems. In a system as large as H2 dissociating on Cu(111), where the present calculations take several million core hours in computational time, this becomes all the more complicated due to the computational restrictions we face.

To assess the fixed node and the locality errors, we would have to increase the accuracy of our trial wave function, meaning that we would have to either enhance the quality of the DFT result (by adjusting the exchange-correlation functional), or go beyond DFT–Slater–Jastrow trial wave functions (e.g., by using backflow or embedded wave function methods). Due to the use of PPs, however, backflow would most likely increase the computational cost beyond reasonable bounds.[81] Similarly, the use of multideterminant wave functions in DMC would be computationally prohibitive for the system studied. In order to get an estimate of the locality errors, we thus retreat to a different approach that will not eliminate the locality errors but that may give us a feeling of the size of the errors. To this end, we repeated our calculations for the 2×2 supercell using a simpler form of the Jastrow function without three-body terms. We thus strongly decrease the quality of our trial wave function since the presence of these terms has previously been shown to reduce the locality errors, especially for large-core Cu PPs. Specifically, for the case of CuH and the large-core PP used in this work, the three-body terms have been shown to account for about half the total locality error in the absolute energy that is incurred without three-body terms.[51] For the dissociation reaction of H2 on Cu(111), total energy changes with Jastrow parametrization of 11.7 ± 0.4 and 13.0 ± 0.4 kcal/mol are observed for the transition-state geometry and the asymptotic geometry, respectively (averaged over k-points). Most of this error can be expected to cancel out when taking energy differences. However, as already becomes clear from the numbers above, the change of Jastrow-parametrization does have an residual influence on the barrier height, increasing ΔE̅DMC by 1.2 ± 0.6 kcal/mol when going from the small to the large parametrization. This number cannot be taken as hard limit for the locality error incurred due to the comparatively large error bar and the fact that we only compare the results for different trial wave functions here but do not strictly converge the locality error to zero. Nevertheless, the 1.2 ± 0.6 kcal/mol difference between the small and the large Jastrow parametrization does indicate that residual effects of the locality approximation still play a role on the very small energy scales we are interested in.

An additional source of error may be the time-step error. Time-step errors can in principle easily be eliminated by going to smaller and smaller time steps. For the system under investigation, however, this approach is not feasible since the computational cost scales inversely with the time step used. The time step of τ = 0.005 au was estimated from preliminary calculations on the CuH molecule using similar Jastrow factors and the same PPs as in this work, and ultimately by requiring the acceptance probability in the DMC propagation to be well above 99%. Specifically, for the large supercell calculations, the acceptance ratio is ∼99.2%. While this is a good estimate to find a reasonable time step, it is worth checking the dependence of the results on the time step. To this end, we again used the 2×2 cell for convergence tests and found that ΔE̅DMC barely changed (by less than 0.2 kcal/mol and well below statistical error bars) when going from τ = 0.005 au to τ = 0.002 au. This excellent result, however, may be due to fortuitous error cancellation: Looking at individual k-points, changes up to 5 ± 2 kcal/mol are observed in some cases and the standard deviation is 3 kcal/mol. Additionally, while a normal probability plot[82] of the changes at different k-points (i.e., a plot of the ordered values of the changes in the DMC barrier height at different k-points versus quantiles of a normal distribution) points to a normal distribution of the changes at different k-points, the width of the distribution is about 1.5 times larger than what one may expect from the statistical uncertainty in the data, which is around 2.3 kcal/mol. Nevertheless, due to the Gaussian distribution of the errors, we may expect this error to be significantly reduced in twist-averaged calculations. Due to the large weight of some k-points in the k-point grid of the 4×4 supercell, however, we may still expect a potential influence on the order of 1 kcal/mol.

Further sources of error are finite-size-related errors. Although we reduced single-particle finite-size effects and finite-size effects stemming from long-range interactions via twist averaging and extrapolating to infinite system size, these corrections may still be insufficient. Especially for the twist averaging, we note that the relationship between DFT results at a particular twist and the corresponding DMC result is not ideally linear (especially for the 2×2 k-point grid, which enters the final result via the finite-size extrapolation). Using more k-points to reduce the DFT based correction would therefore be desirable. In the present calculations this was not possible for computational reasons: In principle, increasing the number of twists comes at no additional cost in the DMC statistics accumulation. In practice, this is not true for two reasons. First, when the number of k-points becomes larger, the relative time spent in equilibrating the walkers increases, and second, the individual twist calculations were run a minimum amount of time to allow decorrelating the error bars with sufficient reliability (error in the error smaller than 9%).

Having considered all possible sources of error in DMC calculations, we should also consider possible errors in the benchmark result of 14.5 kcal/mol. The potential energy surface produced within the SRP approach has allowed several experimental observables to be reproduced, which is certainly in favor of the benchmark value. Nevertheless, there is some uncertainty related to the fitting procedure. For example, in exploratory research investigating the influence of van der Waals interactions on reaction probabilities of H2 at metal surfaces, it has been found that certain experimental observables are equally well reproduced using the optPBE-vdW-DF functional, which, however, gives a significantly higher reaction barrier for the bridge-to-hollow dissociation of H2 on Cu(111) of Eb = 16.4 kcal/mol.[83] This functional, however, has not yet been as extensively tested on reactive and nonreactive scattering of H2 and D2 on Cu(111) as the reference functional cited above.[3,48,84] Therefore, for the time being, the SRP value of 14.5 kcal/mol should be considered the benchmark value.

Conclusion and Outlook

In conclusion, we have presented DMC calculations of the barrier-height of the bridge-to-hollow dissociation of H2 on Cu(111). The calculations are highly challenging, both in terms of the setup, and in the computational effort. We have presented a detailed discussion of the pseudopotentials used, the systematic sources of error and the choice of the geometries that should allow the most accurate possible prediction of the barrier height. Our single-determinant trial wave function based DMC barrier height lies within 1.6 kcal/mol from a semiempirical benchmark value that is fitted against experiments. This is within 1.6 times the statistical error bars of our final results. Further investigation of possible sources of errors strongly point toward a possible influence of residual locality errors and the influence of time-step errors, both of which may influence the result on the order of 1 kcal/mol, while at the same time fixed-node errors cannot be excluded. Other errors one should aim to reduce in future calculations are errors resulting from the finite number of twists.

The calculations presented here thus clearly show the future potential, but also the present weaknesses of DMC for describing molecular reactions on transition metal surfaces. In view of the computational restrictions we already faced in this study, it is a pressing question of how we can improve on these results in the future and whether QMC may indeed offer a route to chemically accurate ab initio benchmarks for molecule–metal surface reactions in the not too distant future. While the ever increasing computational power will trivially solve the computational restrictions leading to time-step errors (DMC scales linearly with the inverse time step) and finite-size errors, more computational power alone will not reduce the possible influence of the locality approximation. On the whole, the sources of error discussed (locality error, finite-size error, and time-step error) may be addressed by combining embedding methods with QMC methods. This will allow the number of layers included in the calculation to be reduced, thus limiting the influence of locality errors of (large-core) PPs used in those remote layers and at the same time the computational cost. Replacing the fourth layer only by an embedding potential would already reduce the computational cost by more than a factor of 2, thus allowing for smaller time steps and more extensive twist averaging. By combining these two state-of-the-art methods, we can thus expect to push the limit of achievable accuracy toward the chemical accuracy needed for accurate first-principle calculations and predictions of catalytic reactions and rates. A combination of embedding theory with quantum Monte Carlo methods might thereby bypass important computational restrictions (very limited cluster size, basis set superposition errors, etc.) faced in standard embedded quantum chemistry approaches. Our results, which already bring QMC within close reach of chemical accuracy and which clearly point to the current deficiencies of the approach, provide an important milestone on the way to this goal and suggest a roadmap to achieve it.