Johanna P. Carbone*a,
Andreas Irmlera,
Alejandro Galloa,
Tobias Schäfera,
William Z. Van Benschotenb,
James J. Shepherdb and
Andreas Grüneis*a
aInstitute for Theoretical Physics, TU Wien, Wiedner Hauptstraße 8-10/136, 1040 Vienna, Austria. E-mail: johanna.carbone@tuwien.ac.at; andreas.grueneis@tuwien.ac.at
bDepartment of Chemistry, Michigan State University, East Lansing, Michigan 48824, USA
First published on 22nd August 2024
We present an application of periodic coupled-cluster theory to the calculation of CO adsorption energies on the Pt(111) surface for different adsorption sites. The calculations employ a range of recently developed theoretical and computational methods. In particular, we use a recently introduced coupled-cluster ansatz, denoted as CCSD(cT), to compute correlation energies of the metallic Pt surface with and without adsorbed CO molecules. The convergence of Hartree–Fock adsorption energy contributions with respect to randomly shifted k-meshes is discussed. Recently introduced basis set incompleteness error corrections make it possible to achieve well-converged correlation energy contributions to the adsorption energies. We show that CCSD(cT) theory predicts the correct order of adsorption energies for the considered adsorption sites. Furthermore, we find that binding of the CO molecule to the top and fcc site is dominated by Hartree–Fock and correlation energy contributions, respectively.
The seminal work of Feibelman et al. shows that the most widely-used local and semi-local density functionals predict the wrong adsorption site preference of CO on the Pt(111) surface compared to experimental findings at low temperatures.2 This puzzle and its far reaching consequences spurred the development of many (semi-empirical) corrections to the employed exchange and correlation density functionals. However, it was shown that achieving a simultaneous and accurate description of surface energies as well as chemisorption energies is impossible using state-of-the-art (semi-)local and hybrid density functionals.3 In contrast, a many-electron perturbation theory approach based on the random-phase approximation (RPA) solves this puzzle in a satisfactory manner.3 However, despite the recent advances of RPA calculations, there still exists the need for more accurate methods that go beyond the RPA. A well-established and highly accurate electronic structure theory is diffusion Monte Carlo, which has also been applied to the CO adsorption problem.4,5 Although both DMC studies predict the correct order of stability of CO adsorption sites, a relatively large discrepancy in the adsorption energy differences is observed. Alternative accurate ab initio methods would be helpful to resolve these and other relevant discrepancies.
In the field of molecular quantum chemistry, Coupled Cluster (CC) theories have established themselves as a class of highly accurate electronic structure theories that can achieve systematically improvable results depending on the level of truncation for the underlying particle–hole excitation operators. However, until recently their application to metallic systems was quite limited and mostly restricted to the level of Coupled Cluster Singles and Doubles (CCSD) particle–hole excitation operators.6,7 As shown in molecular quantum chemistry, highly accurate reaction energies require the inclusion of triple particle–hole excitation operators.8 The most popular triples approximation used in quantum chemistry accounts for these effects in a perturbative manner and is referred to as CCSD(T) theory.9 However, CCSD(T) yields diverging correlation energies for metals.10 Recently, a modified approximation to the triples, denoted as CCSD(cT) was presented.11 This method yields convergent and highly accurate results for the uniform electron gas.11
In this work we apply CCSD(cT) theory to the study of the CO adsorption on Pt(111) and test its reliability using a 2 × 2 surface slab model with 2 layers, which is sufficient to assess the qualitative level of accuracy for a range of electronic structure theories.
Fig. 1 Employed Pt(111) surface slab model with 2 layers and the adsorbed CO molecule on the top (a) and fcc hollow site (b). |
The periodic coupled-cluster calculations reported in this work are performed using our high-performance open-source coupled cluster simulation code, Coupled Cluster For Solids (cc4s).12 The required reference wavefunction and the intermediates are obtained using the Vienna ab initio simulation package (VASP).13–15 For all calculations in VASP, a plane-wave kinetic energy cut off of Ecut = 600 eV is used. The Pt_GW, C_GW and O_GW POTCARS are employed. The smearing parameter is set to σ = 10−4 eV and a convergence criterion of ΔE = 10−6 eV is set for the self-consistent field methods. All Hartree–Fock (HF) calculations employ only integer occupation numbers and are performed using VASP.
In this work all post-HF calculations sample the first Brillouin zone using a single k-point only. Furthermore, the unoccupied HF orbitals for a given plane-wave basis have been computed by setting the number of orbitals/bands to the maximum number of plane-waves in the basis set. We note that the convergence of the CCSD correlation energy is very slow when canonical HF orbitals are employed. To accelerate the convergence to the complete basis set limit, we use approximate natural orbitals that are calculated using eqn (2) of the procedure outlined in ref. 16. After calculating all natural orbitals, a subset of them is chosen and recanonicalized. For coupled-cluster theory calculations we choose the number of unoccupied natural orbitals per occupied orbital to be in a range between 5 and 20. Additionally, for the basis set correction algorithm described below, the Moeller–Plesset Perturbation Theory (MP2) pair energies extrapolated to the complete basis set limit are needed, which are computed using the MP2 algorithm described in ref. 17. With this, VASP can provide all necessary files needed for the CCSD(cT) calculation, including the basis set error correction computed by cc4s. All calculations have been performed using about 8 compute nodes each equipped with 48 cores and 384 GB main memory.
All employed post-Hartree–Fock calculations use a finite number of particle states also referred to as virtual orbitals Nv. The truncation of the virtual orbital basis set introduces a basis set incompleteness error (BSIE). The BSIE vanishes very slowly in the limit of Nv → ∞. In order to reduce the BSIE, we use a pair-specific cusp correction for CC theory.18 This scheme is based on frozen natural orbitals (FNOs) and diagrammatically decomposed contributions to the electronic correlation energy, which dominates the BSIE. To partly account for the BSIE of the (cT) contribution to the CCSD(cT) correlation energy, we rescale the (cT) contribution using the ratio of the MP2 correlation energy from the finite basis and the extrapolated complete basis set limit estimate. This correction was investigated on the level of (T) contributions and has been denoted (T*).19
Further, we simulate the Pt surface slab model using a periodic supercell approach with a finite size, introducing a finite size error (FSE). The HF and CC ground state energy converges slowly with respect to the system or k-mesh size. This partly follows from the fact that correlated wavefunction based theories capture longer ranged electronic correlation effects such as dispersion interaction explicitly. To reduce these FSE one can employ a correction method that takes advantage of the fact that the coupled cluster correlation energy can be expressed as an integral over the electronic transition structure factor multiplied by the Coulomb kernel in reciprocal space. Since finite size errors partly originate from an incomplete sampling of this integral in reciprocal space, an interpolation technique can be used to estimate and correct the FSE. The technical details are described in ref. 20. We note that this technique is not fully justified for metallic systems, which would require a linear interpolation of the transition structure factor in the long-range limit (seen in ref. 21 and 22). Therefore we have only applied this method for the calculated adsorption energies to investigate its order of magnitude. Our computations gave finite size corrections to the adsorption energies on the scale of 0.1 eV. Although this effect is small compared to the total adsorption energies and is not expected to change the reported results of this work in a significant manner, we note that future work will employ the RPA to account for this missing contribution in a more reliable way.
Throughout this work we discuss interaction energies that are defined as
−EXY = EXCO@Y − EXCO − EXPt(111), | (1) |
Moreover, we discuss the difference in adsorption energies between top and fcc sites, which is defined as
−EXtop−fcc = EXCO@top − EXCO@fcc. | (2) |
Layers | EPBEtop−fcc | EPBEtop | EPBEfcc |
---|---|---|---|
2 | −0.12 | 1.82 | 1.94 |
4 | −0.12 | 1.61 | 1.73 |
4 | 1.68 (ref. 23) |
Having established that the two layer Pt(111) surface slab model is sufficient to determine the order of stability of the considered CO adsorption sites, we now turn to the discussion of results obtained at the level of quantum chemical many-electron theories, starting with the HF approximation. It is important to note that the reported HF and post-HF calculations employ a reciprocal representation of the Coulomb kernel that was recently presented in ref. 27. This method samples the reciprocal Coulomb kernel consistently for arbitrary shapes of reciprocal volume elements, which is especially important for surface slabs. We stress that its treatment of the Coulomb singularity in reciprocal space is needed to achieve convergent adsorption energies in the present case. Moreover, we emphasize that convergent HF energies for metallic Pt(111) surfaces with and without adsorbed CO molecules can only be obtained when sampling the first Brillouin zone using k-meshes that are randomly shifted from high-symmetry points, e.g., Γ. This behaviour is attributed to the presence of degenerate orbital energies at high-symmetry points that can not be properly accounted for in HF calculations using integer occupation numbers, leading to convergence problems of standard iterative self-consistent field solvers as implemented in VASP. However, as shown in ref. 6 and 11, we avoid these problems by averaging over HF energies computed using randomly shifted k-meshes. The standard error of the mean is employed as a measure of convergence for this Monte Carlo integration procedure over the first Brillouin zone.
Table 2 summarizes the obtained interaction energies and their difference at the level of HF theory for three different k-meshes. The first observation we report is that EHFtop−fcc converges faster than individual interaction energies with respect to the employed k-mesh size. Increasing the k-mesh from 3 × 3 × 1 to 4 × 4 × 1 changes EHFtop−fcc by 0.16 eV, which is in a similar order of magnitude as the error bars of the EHFtop−fcc 4× 4 × 1 result that originates from using about 10 random shifts in the calculations. It is noteworthy that HF theory predicts the experimentally observed order of stability for the top and fcc adsorption sites, with EHFtop−fcc ≈ 1.2 > 0 in contrast to EPBEtop−fcc = −0.12 < 0. EHFtop and EHFfcc exhibit a slightly slower convergence with k-mesh size in the beginning. However, the changes from 3 × 3 × 1 to 4 × 4 × 1 are similar to EHFtop−fcc. We note that the CO molecule is not binding to Pt(111) at the fcc site in HF theory.
k-mesh | EHFtop−fcc | EHFtop | EHFfcc |
---|---|---|---|
1 × 1 × 1 | 1.29 (0.14) | 0.69 (0.27) | −0.60 (0.28) |
3 × 3 × 1 | 1.46 (0.11) | 1.18 (0.08) | −0.27 (0.08) |
4 × 4 × 1 | 1.62 (0.12) | 1.24 (0.07) | −0.36 (0.16) |
We now discuss the basis set convergence of the correlation energy contributions to the fcc hollow and the difference between top and fcc hollow interaction energies. Although we employ recently introduced corrections to reduce the BSIE of CCSD and (cT) correlation energy contributions, it is necessary to confirm that we employ basis set sizes that allow for achieving sufficiently well converged estimates. Tables 3 and 4 summarize CCSD and (cT) correlation energy contributions for different ratios Nv/No, respectively. Here, Nv represents the number of natural orbitals, and No denotes the number of occupied orbitals. It is shown that a value of Nv/No = 10 is converged to within a few ten meV compared to slightly larger basis set sizes for ECCSD-corr.fcc and E(cT)-corr.fcc. In the case of ECCSD-corr.top−fcc and E(cT)-corr.top−fcc, this is already achieved even when compared to Nv/No = 5, illustrating that the difference in adsorption energies benefits significantly from BSIE cancellation.
Nv/No | ECCSD-corr.top−fcc | ECCSD-corr.fcc |
---|---|---|
5 | −0.74 | 1.62 |
10 | −0.76 | 1.54 |
15 | 1.55 |
Nv/No | E(cT)-corr.top−fcc | E(cT)-corr.fcc |
---|---|---|
5 | −0.44 | 0.40 |
10 | −0.43 | 0.49 |
15 | 0.52 |
Having established that using Nv/No = 10 suffices to achieve CCSD(cT) correlation energy estimates with BSIE smaller than a few ten meV, we choose to employ Nv/No = 10 and seek to address the issue of convergence with respect to the number of randomly shifted 1 × 1 × 1 k-meshes. Table 5 summarizes the obtained CCSD and (cT) correlation energy contributions to the studied interaction energies and their difference. We employ 11 random shifts for the studied systems to obtain an average contribution. The values in the parenthesis correspond to the error of the mean, illustrating that the computed interaction energies show a relatively large dependence on the chosen shift. However, using a sufficiently large number of shifts allows us to achieve a statistically meaningful estimate of EXtop−fcc, EXtop and EXfcc. It is noteworthy that the error of EXtop−fcc is smaller than the sum of the errors of EXtop and EXfcc because we use the same random shift for different structures when computing differences. The larger errors for EXtop and EXfcc indicate that these estimates depend more on the chosen random shift than EXtop−fcc.
X | EXtop−fcc | EXtop | EXfcc |
---|---|---|---|
CCSD-corr. | −0.78 (0.07) | 0.43 (0.25) | 1.21 (0.18) |
(cT)-corr. | −0.44 (0.11) | −0.02 (0.09) | 0.42 (0.09) |
CCSD(cT)-corr. | −1.21 (0.17) | 0.41 (0.33) | 1.62 (0.22) |
The correlation energy contributions on the level of CCSD(cT) summarized in Table 5 are strongly negative for Etop−fcc, indicating that correlation effects favor CO adsorption on the fcc site opposed to the HF theory preference for the top site.
Finally, we discuss the obtained CCSD(cT) estimates for the CO interaction energies on Pt(111). Table 6 and Fig. 2 summarize all individual contributions already discussed above. We briefly reiterate that binding of CO to the top site is strongly favored on the level of HF theory. The additional CCSD(cT) correlation energy contribution to the interaction energy is about 0.5 eV. In contrast to the top site, HF theory predicts that CO is not bound at the fcc hollow site. Only when adding CCSD(cT) correlation energy contributions, the corresponding interaction energy becomes attractive with a significant contribution from the (cT) approximation. We note that the (cT) contribution for the CO adsorption to the fcc hollow site is relatively large indicating that intermediate to strong correlation effects could be involved.
X | EXtop−fcc | EXtop | EXfcc |
---|---|---|---|
HF (4 × 4 × 1) | 1.62 (0.12) | 1.24 (0.07) | −0.36 (0.16) |
CCSD-corr. | −0.78 (0.07) | 0.43 (0.25) | 1.21 (0.18) |
(cT)-corr. | −0.44 (0.11) | −0.02 (0.09) | 0.42 (0.09) |
CCSD(cT) | 0.41 (0.29) | 1.65 (0.40) | 1.26 (0.38) |
Adding all contributions together gives ECCSD(cT)top−fcc ≈ 0.4 eV, in qualitative agreement with recent DMC calculations yielding ≈0.25 eV (ref. 5) and ≈0.76 eV.4 RPA calculations yield a significantly smaller difference between adsorption energies of 0.08 eV.3 We note that the main uncertainty in our CCSD(cT) estimate originates from the random k-mesh shift averaging procedure. Future work will focus on reducing this error and will hopefully help to further resolve this discrepancy.
The employed surface slab models using 2 layers are not large enough to allow for a direct comparison of all computed adsorption energies to experiment. Still using the 2 layer model, our CCSD(cT) interaction energy estimates for the top and fcc hollow site are 1.65 eV and 1.26 eV, respectively. Considering that, on the level of DFT-PBE we have estimated that the difference between 2 and 4 layers for the absolute adsorption energies are on the scale of about 200 meV, this would bring our estimates closer to experimental findings of about 1.29 eV that include zero-point vibrational effects.25 We also note that our findings for adsorption energy differences are already converged with respect to the number of layers and show a qualitative agreement with experiment and other high level theories.3–5 In particular, CCSD(cT) theory predicts the top adsorption site to be energetically more stable than the fcc hollow site by about 0.4 meV. This finding is larger than for the RPA results3 and lies between previously reported DMC findings.4,5 We emphasize that this brings periodic CC theories one step closer to play an important role in producing accurate benchmark results for technologically relevant surface chemistry problems studied in, e.g. heterogeneous catalysis.
An important observation of the present work is that the adsorption energy for the top site is dominated by electrostatic energy contributions, while the fcc site is dominated by correlation energy contributions. This in part might explain why most approximate exchange and correlation density functionals exhibit difficulties to predict the correct order of adsorption energies for different sites.2
Although the computational cost of the employed CC theories is significantly larger than that of standard DFT-PBE calculations, all required computations can be performed on a few modern multi-core compute nodes equipped with a few TB of main memory. Still, a remaining technical challenge that was identified in the present study is that standard iterative self-consistent field HF solvers converge frustratingly slow for metal surfaces, requiring up to several hundreds of steps until convergence is reached. In combination with the fact that we need to perform many HF calculations for different random k-mesh shifts, this creates a computational bottle neck that needs to be addressed in future work. Furthermore, the previously developed finite size correction method described in ref. 20, which highlighted the role of anisotropy, needs to be adapted to address the correct interpolation limit in a similar way to ref. 21 and 22 in order to account for the slower convergence of finite size error in metallic systems.11 However, we expect that the problems mentioned above can soon be resolved, paving the way for highly accurate CC theory calculations especially when combined with reduced scaling approaches such as suitable embedding techniques and local approximations.31–34
This journal is © The Royal Society of Chemistry 2024 |