Reaction barriers on non-conducting surfaces beyond periodic local MP2: Diffusion of hydrogen on α-Al₂O₃(0001) as a test case

Surfaces are natural and in many cases very efficient catalysts for various types of reactions. As a result, the field of surface chemistry enjoys an intense focus from both experimental and theoretical perspectives. Due to the complexity of the processes, however, and extended nature of involved systems, the accurate theoretical description of reaction energetics is very challenging. This is especially critical for accurate calculation of reaction barriers, which, if estimated poorly, may completely undermine the subsequent kinetic modeling or reaction rate estimates. Recent attempts to achieve chemical accuracy (<1 kcal/mol) also for surface reactions are summarized in Ref. 1.

Standard computational protocols for surface reactions employ periodic super-cell models, treated with plane-wave basis sets and density functional theory (DFT), within the generalized gradient approximation (GGA). However, GGA functionals are notoriously known to substantially underestimate reaction barriers, thus in many cases providing not more than a rough estimate of the picture at hand. In a recent publication,² some of us demonstrated that a hybrid-DFT treatment, or Møller–Plesset second-order perturbation theory (MP2), can deliver a more realistic barrier for hydrogen diffusion on aluminum oxide surfaces. Nevertheless, the scattering of the results observed among these models was still noticeable, so a precise theoretical prediction was yet not possible. The MP2 value for the barrier, however, was between the results of different hybrid DFT functionals, allowing for a speculation that the actual barrier is of that ballpark. However, MP2, although generally known for being accurate for oxide and ionic systems, cannot be considered a safe benchmark. In order to reach certainty about the value of the barrier (at least within reasonable margins), one has to step much higher in the hierarchy of the quantum chemical models and reach at least the coupled-cluster singles, doubles, and perturbative triples [CCSD(T)] level of theory, combined with a basis set of suitable quality.

The CCSD(T) method itself, the “gold standard” of quantum chemistry, is computationally very expensive and, when combined with good basis sets, its applicability is restricted to small molecular systems. For somewhat larger systems, usually, a hierarchical approach is employed, where the basis set limit is reached for computationally less demanding models, while the CCSD(T) correction is then computed using a more modest (although still of decent quality) basis set.³ For even larger systems, the standard (canonical) CCSD(T) quickly becomes prohibitively expensive, even when using modest basis sets (due to the unfavorable N⁷ scaling with the system size) and, thus, can only be applied with additional approximations. One of the most successful approximate CCSD(T) models applicable to extended molecules is the local CCSD(T) model.^4–12 Periodic systems add a further level of complication for high-level quantum chemical models, both of technical and fundamental nature. Until now, several effective schemes and protocols have been developed that do allow for a CCSD(T)-level treatment. First, there are some fully periodic implementations,^13,14 yet these are still very computationally demanding, especially for large unit cells, as in our case. Another possible strategy is to employ representative finite clusters, either exclusively^15,16 or for correcting a low-level periodic treatment.^17–28 However, the applicability of this approach depends on the existence of clusters that can adequately mimic the parent solid. For the system under study here—aluminum oxide—construction of such clusters is highly problematic.

In this work, we employ a different approach: the embedded fragment CCSD(T) correction, which is a powerful method for accurate description of local chemical processes in solids, for example, reactions of adsorbed species. It offers a possibility to define a fragment seamlessly via localized orbitals rather than explicit atoms, so no bond-cutting is required to create an atomic cluster. One possible realization of such an approach is to embed a wave function treatment in a DFT environment.^29–42 An alternative to the DFT embedding is the so-called Hartree–Fock (HF) embedding, where the fragment becomes implicitly embedded in the frozen Coulomb and exchange field of the true periodic HF solution. The HF description of the environment is usually not as accurate as a DFT one, especially so for metals. On the other hand, HF embedding allows us to treat the whole system in a unified seamless quantum chemical framework, operating with frozen and active orbital spaces. Conceptually, it is similar to the standard frozen-core approximation. The only difference is that here the frozen/active spaces are determined not according to the orbital energies in the canonical representation but rather on the basis of spatial proximity to the feature of interest in the direct space representation.

In recent years, periodic HF embedding has drawn substantial interest and several groups have developed various flavors of this approach.^43–48 In this work, we employ our recent interface to the Molpro code,^49,50 which passes the one-electron Hamiltonian matrix and two-electron integrals in the basis of the fragment orbitals to the molecular code.

We use the embedded fragment to evaluate the CCSD(T) correction to the periodic local MP2 (LMP2) treatment, i.e., to correct the MP2 method error. All other significant errors of the periodic LMP2 treatment are analyzed and corrected in the periodic regime. The goal of this paper is to provide a CCSD(T) benchmark for the hydrogen diffusion barrier on α-Al₂O₃(0001) and formulate a protocol for high precision calculations of quantities relevant to the description of the kinetics of catalytic reactions on non-conducting surfaces.

The initial and final states of the considered reaction represent two stable geometries of dissociative water adsorption on α-Al₂O₃(0001) surface, which according to the nomenclature of Refs. 2 and 51 are marked as “1-4” and “1-2,” respectively (see Fig. 1). The 1-2 structure is a global minimum for water adsorption on this surface with 1-4, according to LMP2 or B3LYP-D3, being about 0.4 eV higher in energy.² The reaction from 1-4 to 1-2 represents the diffusion of the hydrogen atom from one surface oxygen atom to another.

The barrier for this reaction was calculated in Ref. 2 using DFT and LMP2 methods. Even when discarding the Hartree–Fock (HF) outlier, the resulting free energy values (at 300 K) for the barrier scatter considerably: from 0.27 eV (PW91) to 0.58 eV (B3LYP-D3). The MP2 barrier was found between these values, closer to the B3LYP result. The main difference in the results was originating from the electronic contribution to the barrier, as the zero point energy (ZPE) and thermal contributions, computed with B3LYP and Perdew–Burke–Ernzerhof (PBE), are comparable (−0.11 eV vs −0.15 eV, respectively). Consequently, in this work, we focus solely on the electronic contribution to the barrier,

The α-Al₂O₃(0001) surface in our calculations is represented by a nine-layer slab. Since we employ atomic-orbital (AO) basis sets, we do not replicate the slab in the $z⃗$-direction, commonly done in plane-wave calculations. The 1-4 and transition state (TS) structures of Eq. (1) were taken from Ref. 2, where they were optimized with B3LYP within a 63-atom super-cell. The explicit structure specifications can be found in the supplementary material.

We note that in this work, we do not aim at investigating the possible errors due to the slab model of the aluminum oxide surface or inaccuracy of DFT structure optimizations. First, it is very difficult to go beyond a slab-model and DFT for geometry optimization, ZPEs, and thermal effects. At the same time, DFT geometries and vibrational frequencies are usually not inaccurate. Hence, the main focus is commonly put on the electronic part of the barrier energy evaluated from two single point calculations on DFT optimized structures. If it is captured with high accuracy, the overall free energy of the barrier becomes a good reference for experiment.

Second, a high precision electronic energy difference alone is a valuable benchmark for DFT. As DFT remains the only practical option for large scale calculations, benchmarking the performance of the employed functional or functionals on key quantities becomes absolutely instrumental.

In this section, we will describe in detail approximations and possible sources of error, namely, basis set choice, thresholds and truncations inherent to the approximations employed, level of theory, and fragment choice.

In this work, we calculate the barrier [Eq. (1)] using a hierarchy of quantum chemical models. We start from the periodic AO-based HF method of the Crystal code.⁵² Unfortunately, rich standard molecular basis sets are usually inapplicable for periodic calculations without additional adjustment, as even slightly diffuse orbitals cause quasi-linear dependencies among basis functions and numerical instabilities in the SCF procedure. In this work, for most of the calculations, we use a VTZ-quality basis set, optimized for periodic systems, which was already employed in Ref. 2 (basis set AO3). We refer to this basis set as VTZ. The explicit specification of this basis set as well as all other basis sets used here is given in the supplementary material.

Yet, it may become possible to employ rich molecular basis sets (at least some of them) in periodic calculations, if the corresponding orbitals are placed not everywhere but only on a few selected atoms. This opens the possibility to pick up an atomic fragment spatially close to the reaction site and use a molecular basis set hierarchy for the atoms of this fragment only. In this way, we employ the cc-pVTZ and cc-pVQZ basis sets of the Dunning basis set family⁵³ on fragments defined below. The fragments only involve oxygen and hydrogen atoms, as cc-pVXZ basis sets for aluminum are very diffuse and not applicable, even when put on just a few atoms. Due to technical reasons (the AOs of higher angular momentum than f are not yet implemented in Cryscor), the g-AOs were omitted from the cc-pVQZ basis of oxygen. All atoms outside the chosen fragment still feature the VTZ basis.

Next, for the accurate post-HF description of dispersion, diffuse high-angular momentum orbitals are needed in the basis.⁵⁴ We include them via the dual basis set technique,⁵⁵ i.e., only for the post-HF treatment. We denote such basis sets as (A)VTZ or (aug-)cc-pVXZ:⁵⁶ the latter correspond to the cases when the Dunning basis sets are used on the fragment atoms.

Within the dual basis set scheme, it is also possible to perturbatively estimate the contribution of the appended orbitals to the HF energy through the so-called first-order singles,^55,57

Here, N_k is the number of k-points in the Brillouin zone sampling, ϵ_a(k) and ϵ_i(k) are the virtual and occupied orbital energies, respectively, and f_ai(k) are the occupied–virtual elements of the Fock matrix, which are not all zero in the dual basis set scheme. This approach is also employed to add the tight orbitals from the cc-pwCVTZ basis set⁵⁸ for Al atoms for evaluating the correlated core contribution to the barrier (see Sec. III B 1).

For the periodic post-HF treatment, we use a local MP2 model. The closed-shell local MP2 energy may be written as

Here, i and j are the localized occupied orbitals, a and b are the localized virtual orbitals, and (ia|jb) are two-electron integrals in the chemical notation,

and $Tabij$ are the MP2 amplitudes obtained by solving the MP2 equations.^59,60 First, we note that in the local (direct space) representation, the orbitals are real, so there is no need to explicitly specify the complex conjugation in (4). Second, in contrast to the canonical basis, the local MP2 equations are not trivially invertible. However, these equations are linear and an iterative LMP2 equation solver usually converges quickly.

In (3), the summation ranges are denoted as $P$⁠, which represents the pair-list for the occupied orbitals, and [ij] for the virtual orbitals specific for each occupied orbital pair ij. In the direct space representation for the correlation energy per cell, one of the indices (e.g., “i”) has to be restricted to the reference cell. An additional restriction of other indices implies an approximation to the MP2 method. In this work, we will explore the effect of the involved approximations that truncate the pair list and virtual space.

We start with the frozen core approximation, very common in correlated calculations, which restricts the pair-list $P$ to the valence occupied orbitals only. It assumes cancellation of the core–valence and core–core correlation contribution in energy differences related to chemical processes. Although this assumption is generally reasonable, the influence of the core correlation associated with the high core orbitals of the metal atoms (in our case 2s2p-core of Al) may not be entirely negligible. We investigate this by comparing the periodic LMP2 barriers with and without the 2s2p-core of Al atoms correlated,

where “FC” stands for the “frozen core,” and “FC_1s” denotes that only 1s-core of oxygen and aluminum atoms is kept frozen. These calculations employ the (C)VTZ and (AC)VTZ basis sets that contain the special tight orbitals from the cc-pwCVTZ basis of the Al atoms.⁵⁸ 

Next, we focus on the so-called local approximations. The first approximation of this kind is the pair-approximation that further restricts the pair-list. In molecular local correlation methods, the pair-list is truncated for efficiency reasons (in particular, to reach linear scaling), as the pair-energy contributions decay with the inverse sixth power of the distance between the centers of the local orbitals.^9,61 In a periodic system, it is in principle not possible to include all Wannier function (WF) pairs in the pair-list $P$⁠, even formally, as their number is infinite. In order to investigate the influence of the pair approximation, we progressively expand the pair-list until convergence of $ΔEcorrLMP2$ is achieved.

The second essential local approximation is the so-called domain approximation, which restricts the virtual space for each pair ij to the pair domains [ij]. A pair domain is thus a pair-specific set of virtual orbitals.

In this work, we use two types of virtual orbitals and, therefore, two types of pair domains. The virtuals of the first type are the projected atomic orbitals (PAOs).^57,62 These orbitals are AOs, which are projected out from the occupied manifold and so are atom-centered. A PAO pair domain [ij] consists of all PAOs, centered on the atoms, that are spatially close to the occupied orbitals i or j of a pair. It is convenient to define the minimal PAO domains, which include PAOs from only one atom (for lone-pair WFs) or two atoms (for bond WFs) per occupied orbital. The domains can then be systematically expanded by adding the PAOs from nearest-neighbor atoms (iext = 1), second-nearest-neighbor atoms (iext = 2), etc. In aluminum oxide, which is to a large extent ionic, the WFs are centered on oxygens, so the corresponding minimal domains consist of the PAOs of just single oxygen atoms. Extension to the nearest-neighbors (iext = 1) adds a shell of the aluminum atoms surrounding the initial oxygen, and a second nearest-neighbor (iext = 2) is the next shell of the oxygen atoms.

Restriction of the virtual space to the PAO pair-domains leads to a compact pair-specific virtual space, and, at the same time, a common global set of virtual orbitals remains, which facilitates calculation of integrals and other intermediates. Yet, a known issue of PAO-based domains is a possible imbalance between the virtual spaces for different points on the potential energy surface, at least when small domains are used. The reason is that the virtual space can be varied only rather coarsely, by adding or removing all the PAOs of a given atom.

As an alternative to PAOs, in periodic local MP2, one can also use the so-called orbital specific virtuals (OSVs).^63,64 OSVs for a given occupied orbital i are defined as the MP2 virtual natural orbitals for a diagonal pair ii. To truncate the virtual space by including only the essential OSVs, an energy cutoff is usually used (in Cryscor 10⁻⁵ hartree by default). An OSV pair-domain for a non-diagonal pair ij is defined as a union of the OSVs corresponding to the orbitals i and j. There is also one aspect that is specific to our OSV implementation: in order to counterbalance the bias of the OSVs toward the short-range correlations, the most diffuse PAOs of the corresponding minimal PAO domain (one shell of PAOs per angular momentum per center) are added to the set of the OSVs.

The OSVs and, thus, the OSV domains automatically adapt to the change in the structure, so the imbalance and discontinuity of the potential surface are much weaker than with PAO domains. On the downside however, with OSVs, it turns out to be difficult to explore convergence with domain size in periodic calculations, as tightening of the cutoff energy threshold requires progressive expansion of the fitting basis.⁶⁴ 

The next approximation in the periodic local MP2 method concerns the evaluation of the electron repulsion integrals (ia|jb) of Eq. (4). These integrals are used both in the LMP2 energy expression [Eq. (3)] and in the LMP2 amplitude equations.⁶⁰ In order to gain efficiency, they are not evaluated and transformed as 4-index AO integrals but rather assembled from 3- and 2-index quantities by means of the density fitting (DF) approximation. In our approach, we use local robust density fitting, which is very efficient in bulk systems.^65,66 Yet, this is an approximation, as the integrals are not exact. DF is known to be generally very accurate for energy differences. It is, however, of interest to check its accuracy also for this system.

The DF approximation in MP2 implies decomposition of occupied–virtual orbital products ρ_ia(r) = φ_i(r)φ_a(r) in an auxiliary (or fitting) basis {ϕ_P},^65,66

where [ia]_fit is the fitting domain that includes the fitting functions centered on the atoms supporting the product density ρ_ia(r). As the fitting functions, we use Poisson-type orbitals (PTOs),^65,67,68 directly converted from the MP2-optimized Gaussian-type-orbital (GTO) fitting basis sets.⁶⁵ In addition to the momentless PTOs, in order to capture the moments of the fitted product densities, the PTO fitting basis sets are appended with a few GTOs: one shell per angular momentum per center in the fitting-domain (s-type GTOs are not included, as the densities ρ_ia are chargeless).

In the calculations presented here, we used the PTO basis converted from the GTO fitting basis set optimized for MP2 with aug-cc-pVTZ orbital basis.⁶⁹ In order to check the accuracy of the DF approximation, we compared these results against those calculated with aug-cc-pV5Z-optimized fitting basis sets. The deviation in $ΔEcorrMP2$ between these two fitting basis sets was less than 10⁻⁵ eV, suggesting that the error due to the density fitting approximation is absolutely minute.

The local correlation scheme allows not only for a highly efficient correlated treatment, but it also enables partitioning of correlation part of the interaction energies into physically interpretable components. Usually, such partitioning is used in the context of inter-molecular interactions, where one can isolate the inter-molecular as well as intra-molecular components of the interaction energies. In our approach, we also employ a similar partitioning scheme to investigate how much of the correlated correction to the barrier stems from the correlation inside the fragments (see Secs. III D and III C).

According to Eq. (3), the total correlation energy can be defined via the pair energies,

With this, the correlation energy and, thus, the correlation part of the barrier can be partitioned into the intra-fragment, extra-fragment, and inter-fragment-environment components,

with

Furthermore, different pair restriction thresholds can be applied to different pair classes individually. In this way, the pair classes that do not contribute substantially to the quantity of interest (in our case the barrier) can be restricted more harshly. This allows for substantial gains in efficiency without affecting the accuracy.

It is well known that the correlation energy is sensitive to the basis set quality. Therefore, it is essential to estimate the basis set limit for the correlated part of the barrier, instead of just remaining at the AVTZ basis set level (as in Ref. 2). To this end, we employ two techniques. First, we extrapolate $ΔEcorrLMP2$ using the standard inverse-cubic formula.⁷⁰ The extrapolation is based on the results with the (aug-)cc-pVTZ and (aug-)cc-pVQZ basis sets placed on the fragment atoms (see Subsection III D). The basis set extrapolation is done with both the PAO and OSV approaches. The second scheme for $ΔEcorrLMP2$ at the basis set limit is the periodic LMP2-F12 method, which, at this point, is only implemented with PAOs.⁶⁶ 

Finally, in order to approach chemical accuracy, we also need to correct for the MP2 method error itself. Here, the correlated treatment is extended to the CCSD(T) level via the embedded fragment approach. In this method, a group of localized occupied and virtual orbitals from the periodic calculation, which are spatially close to the reaction, define a fragment. The embedding of the fragment in the periodic HF solution is incorporated in the fragment’s one-electron Hamiltonian $hpqfrag$⁠, which is constructed from the periodic Fock matrix F^cryst and includes the Coulomb and exchange potential from the rest of the crystal,⁵⁰ 

where p and q (as well as r and s below) are the general orbitals of the fragment (occupied or virtual). The two-electron part of the fragment Hamiltonian elements (pq|rs) is computed in Cryscor using the localized orbitals from the actual periodic calculation.⁵⁰ Importantly, the fragment orbitals are, by construction, orthogonal to the occupied orbitals of the environment, allowing for a clean and seamless cut between the fragment and embedding spaces.

The Hamiltonian matrix elements $hpqfrag$ and (pq|rs) are transferred to the molecular code Molpro⁷¹ via a FCIDUMP interface,⁷² where any available canonical post-Hartree–Fock treatment can be applied. The method correction is then defined as

As noted above, the focal high-level method for the correction, denoted by “meth” in (14), is CCSD(T), providing δ^CCSD(T)ΔE. We will, however, also test other less expensive methods, such as CCSD and distinguishable cluster with singles and doubles (DCSD).^73–75 

As discussed above in Secs. III A and III C, we employ fragments for the basis set extrapolation of the correlation energy and for the evaluation of δ^CCSD(T)ΔE. We use three fragments of progressively increasing size, which is shown in Fig. 2. The oxygen atoms of the fragments are marked in blue.

For the basis set extrapolation, the atomic orbitals of the (aug-)cc-pVTZ or (aug-)cc-pVQZ basis set are placed on the fragment oxygen and hydrogen atoms. The rest, including all aluminum atoms, remain with the (A)VTZ basis. The standard Dunning basis sets of aluminum contain too diffuse AOs that cannot be used for periodic calculations. However, since all the valence WFs are localized on oxygen atoms, a triple-zeta basis set level for aluminum is expected to be sufficient.

For the embedded fragment treatment, the fragments are defined by the valence WFs, which are centered on the marked oxygen atoms: four per atom. The fragment’s virtual space is spanned by the PAOs from the minimal domains of these WFs (see Sec. III B 3). This includes the PAOs of the oxygen atoms of the fragment, the hydrogen atoms, and one aluminum atom bound to the adsorbed OH group. For the VTZ basis and Fragment 2, we also expanded the virtual basis to the PAOs from the next nearest-neighbor atoms, which added another four aluminum atoms (see Fig. 2). Table I summarizes the specification of the occupied and virtual space in the fragments.

We note that the actual number of the orthogonal virtual orbitals of the fragment is slightly smaller than the number of PAOs. In fact, PAOs form a non-orthogonal and even generally redundant set (as by construction there are as many PAOs as AOs). Before the in-fragment canonicalization, the fragment’s PAOs are orthogonalized, employing the same procedure as in PAO-based local correlation methods.⁷⁶ It takes the eigenvectors of the PAO overlap matrix as the orthonormal virtual set to be canonicalized. However, to get rid of linear dependencies (or quasi-linear dependencies), all the eigenvectors with the eigenvalues below a certain threshold (the default value in Cryscor is 10⁻⁴) are excluded.

We start our discussion with the HF barriers. Table II summarizes the values for ΔE_HF obtained with different basis sets. The diffuse orbitals were added only within the dual basis set approach; hence, their contribution was captured via the singles correction of Eq. (2).

As follows from Table II, the HF barrier is not very sensitive to the basis set level, at least starting from (A)VTZ. However, there is scattering of the results of about 0.02 eV, with no clear convergence pattern. For estimating ΔE_HF, we take the value of 1.13 eV, which was obtained using the (aug-)cc-pVQZ basis on the biggest fragment [and (A)VTZ on the rest]. However, due to the singles scatter, the uncertainty cannot be reduced to less than 0.02 eV. As noted above and will also be seen later, the HF barrier is much too high when compared to the correlated methods.

Next, we move on to the correlation energy contribution to the barrier, where we first focus on the periodic local MP2 treatment. In Fig. 3, we analyze the locality of the correlation effect on the barrier using LMP2(OSV)/(A)VTZ. For Fragment 1, the intra-fragment contribution is virtually zero, suggesting that the electron correlation of the electrons near the OH group does not influence the barrier. However, already for Fragment 2, ΔE_intra-frag becomes the dominant component of $ΔEcorrMP2$⁠. This suggests that with fragments 2 and 3, an accurate description of the intra-fragment correlation would also provide an accurate overall description.

As a next step, we explore the accuracy and possible sources of errors in the periodic LMP2 treatment due to the additional approximations. First, as may be seen from Table III, $ΔEcorrMP2$ converges relatively fast with the pair-list cutoff distance, suggesting that the pair approximation does not introduce any sizable error.

The domain error, however, is more difficult to estimate and eliminate, predominantly so in the PAO approach. Since the domain error and the basis set incompleteness error are of the same origin (due to restriction of the virtual space), we treat these two errors together. First, we note that, according to Table IV, the basis set extrapolated value depends very mildly on the fragment that features the actual (aug-)cc-pVTZ and (aug-)cc-pVQZ basis functions.

Conversely, the resulting energies depend much stronger on the PAO-domains. In Fig. 4, we show the basis set limit estimates for $ΔEcorrLMP2$ obtained with the basis set extrapolation and LMP2-F12 method using different PAO domains. For the minimal domains (iext = 0) and especially for the iext = 1 domains, the results depend strongly on the scheme employed to approach the limit. For iext = 2 domains, on the other hand, the extrapolated and F12 results agree as they agree also with the OSV-based calculations.

The problem occurring with smaller domain calculations seems to be related not to the dependence of the barrier on the basis set but rather to the imbalance between the virtual space of the 1-4 and transition state structures. This follows from the fact that the correction on top of the (A)VTZ value turns out to be rather small: +0.003 eV for the extrapolation and 0.016 eV for F12 (see the supplementary material). The OSV approach is much less prone to the imbalance issues, as the virtual orbitals naturally adapt to the change in the structure. We, therefore, take the OSV (aug-)cc-pV(TQ)Z-extrapolated value of −0.54 eV as the MP2 reference result. This value also does not deviate strongly from the iext = 2 PAO results: 0.006 eV from the PAO extrapolated result and 0.019 eV from the LMP2-F12 value. However, the 0.02 eV discrepancy has to be taken into account as a possible uncertainty of the result.

Next, we draw our attention to the core correlation. The frozen core approximation is standard in most post-HF applications, but when aiming at very high precision, the core correlation contribution may be non-negligible.^28,77–79 This is apparently also the case here: when the 2s2p electrons of aluminum are correlated, $δcoreΔEcorrLMP2$ adds around +0.06 eV, regardless whether PAOs, OSVs, or (AC)VTZ or (C)VTZ are used (see Table 7 of the supplementary material for the detailed results).

Finally, we correct the MP2 method error by evaluating the δ^CCSD(T) correction. Table V compiles the values for the correction computed for different fragments and basis sets. The CCSD(T)-correction to MP2 of +0.07 eV is not negligible and (if excluding Fragment 1, which is clearly too small and inadequate for this system) remains very stable across the fragments and basis sets.

We also tested a few other methods for the method correction: CCSD, DCSD, and the spin-component-scaled-(SCS-)DCSD.⁷⁵ CCSD noticeably overestimates barrier, as δ^CCSDΔE is approximately twice as large as δ^CCSD(T)ΔE. DCSD results are much closer to CCSD(T), especially so with the SCS version of DCSD, which virtually coincides with CCSD(T). Interestingly, the dependence of the correction on the fragment is stronger for CCSD or DCSD than for CCSD(T) but is still rather weak.

Having thus examined and eliminated the main errors in the quantum chemical treatment, we are in a position to provide the benchmark value for the barrier for hydrogen diffusion from the 1-4 to 1-2 structures on the aluminum oxide surface, which amounts to 0.72 eV. In Table VI, we compile the final results and compare them to the DFT values of Ref. 2. It is clear that the HF barrier is much too high. The MP2 result is of reasonable quality, but when using the frozen core approximation, the barrier is yet substantially lower than the reference. The core correlation and CCSD(T) corrections are both noticeable and of the same sign, accumulating the error of MP2 with the frozen core.

Among the DFT functionals, the hybrids are clearly better performing than GGAs, which grossly underestimate the barrier. Interestingly, the B3LYP-D3 value is very close to our benchmark. The dispersion contribution to the barrier is not negligible but fairly small, as is revealed by the D3 component amounting to +0.04 eV. The quick convergence of the LMP2 barrier and the post-MP2 corrections with the fragment size suggests that the essential dispersion part comes from the region directly surrounding the reaction. For such a reaction, where the reactant and the product are both adsorbed on a surface, the long range dispersion contributions are indeed expected to cancel out to a large extent in the energy differences.

The goal of this contribution is two-fold. First, we report a highly accurate benchmark value of 0.72 eV for the barrier of hydrogen diffusion on α-Al₂O₃(0001) surface. The main uncertainty originates in the basis set convergence, which could not be reduced below 0.02 eV, both in HF and in MP2. However, even if these errors accumulate, it still remains within chemical accuracy. The GGA DFT barriers in comparison to this benchmark are clearly severely underestimated: by about 0.3 eV. Hybrid functionals perform better, in this case especially so B3LYP-D3. We do not have a direct experimental reference to compare to (rather, on the contrary, the value obtained in the work can serve as a reference for the experiment). However, the high barrier obtained here can be qualitatively supported by the fact that the hydrogen vibrations in positions 1-4 are observed experimentally,⁵¹ suggesting that the hydrogen atoms diffuse from 1-4 to 1-2 rather slowly.

Second, this work demonstrates a protocol for calculating reaction barriers on non-conducting surfaces. The periodic HF and local MP2 with (A)VTZ-quality basis sets are expected to be already reasonably close to the basis set limit. OSV-based LMP2 calculations are much more robust than the PAO ones as concerns the balance of the virtual space between different points on the potential surface. The basis set limit can be calculated either by extrapolation of the OSV-LMP2 energies (the corresponding Dunning basis set can be used only for the fragments surrounding the reaction) or via an LMP2-F12 calculation. Since the implementation of the OSV-based periodic local F12 method is yet underway, for the PAO-based version of F12, one should use very big domains (e.g., iext = 2); otherwise, unbalanced results may be obtained.

When metal atoms are present, the standard frozen-core approximation may not be sufficiently accurate. In our case, the correlation contributions from the upper core of Al are sizable and should be evaluated. Another important correction is the CCSD(T) correction on top of MP2. In order to evaluate this correction, an embedded fragment approach is practical and very efficient. If all the interesting steps of the reaction involve only adsorbed species, the long-range dispersion is expected to cancel to a large extent in the energy differences, reducing the essential active region to a local space around the reacting species.

See the supplementary material for the specification of the structures, basis sets, and other computational parameters, as well as a detailed compilation of the obtained results.

D.U. acknowledges the Deutsche Forschungsgemeinschaft for the support (Grant Nos. US 103/1-2 and INST 276/714-1). P.S. acknowledges the support of the Deutsche Forschungsgemeinschaft under Germany’s Excellence Strategy (Grant No. EXC 2008/1-390540038).

The authors have no conflict of interest to declare.

The data that support the findings of this study are available within the article and its supplementary material.