Adsorption and vibrational spectroscopy of CO on the surface of MgO from periodic local coupled-cluster theory

Abstract

The adsorption of CO on the surface of MgO has long been a model problem in surface chemistry. Here, we report periodic Gaussian-based calculations for this problem using second-order perturbation theory (MP2) and coupled-cluster theory with single and double excitations (CCSD) and perturbative triple excitations [CCSD(T)], with the latter two performed using a recently developed extension of the local natural orbital approximation to problems with periodic boundary conditions. The low cost of periodic local correlation calculations allows us to calculate the full CCSD(T) binding curve of CO approaching the surface of MgO (and thus the adsorption energy) and the two-dimensional potential energy surface (PES) as a function of the distance from the surface and the CO stretching coordinate. From the PES, we obtain the fundamental vibrational frequency of CO on MgO, whose shift from the gas phase value is a common experimental probe of surface adsorption. We find that CCSD(T) correctly predicts a positive frequency shift upon adsorption of +14.7 cm⁻¹, in excellent agreement with the experimental shift of +14.3 cm⁻¹. We use our CCSD(T) results to assess the accuracy of MP2, CCSD, and several density functional theory (DFT) approximations, including exchange correlation functionals and dispersion corrections. We find that MP2 and CCSD yield reasonable binding energies and frequency shifts, whereas many DFT calculations overestimate the magnitude of the adsorption energy by 5–15 kJ mol⁻¹ and predict a negative frequency shift of about −20 cm⁻¹, which we attribute to self-interaction-induced delocalization errors that are mildly ameliorated with hybrid functionals. Our findings highlight the accuracy and computational efficiency of the periodic local correlation for the simulation of surface chemistry with accurate wavefunction methods.

---

I. Introduction

The adsorption of a molecule on to the surface of a substrate is the first step of many interfacial physicochemical processes such as heterogeneous catalysis^1,2 and gas storage and separation.^3–5 Accurate determination of the adsorption energy, and more generally the potential energy surface (PES), is thus vital for understanding the structure and dynamics of the adsorption process, which in turn facilitates the interpretation of relevant interfacial experiments (e.g., molecular beam,^1,6,7 temperature-programmed desorption,^8–11 and surface vibrational spectroscopies^12,13). Nonetheless, this has been proven a challenging computational task over the past two decades.^14,15 The major difficulty stems from the inherently weak yet correlated nature of surface–adsorbate interactions, which necessitates achieving sub-kJ mol⁻¹ accuracy for reliable comparison with experiments^14–17 and development of accurate force fields for molecular dynamics.^18–21 The demanded high accuracy is typically out of reach by density functional theory²² (DFT) whose performance can depend sensitively on the choice of approximate exchange–correlation functionals (especially because of self-interaction errors) and dispersion correction methods.^1,23–26

In principle, correlated wavefunction theories that include the correct many-body physics hold promise for computational surface science. In practice, however, their applicability is often hindered by the high computational cost required to reach convergence with respect to the one-particle basis set size, the supercell and slab size, and the level of electron correlation treatment. In this work, we focus on CO adsorption on the MgO (001) surface, which is an extensively studied model system for high levels of theory.^16,17,27,28 Recently, Shi and co-workers¹⁷ applied coupled-cluster theory with singles, doubles, and perturbative triples²⁹ [CCSD(T)], commonly known as the “gold standard” of main-group quantum chemistry, within an embedded-cluster framework³⁰ to produce their best estimate of the adsorption energy, −19.2 ± 1.0 kJ mol⁻¹, which was in good agreement with that calculated by canonical periodic CCSD(T) and diffusion Monte Carlo, as well as the value estimated from temperature-programmed desorption experiments¹¹ (about −20.0 kJ mol⁻¹).

In a preprint version of the present manuscript,³¹ we applied our recently developed periodic local CCSD(T) method^32,33 to calculate an adsorption energy in good agreement with the above values. Due to the use of local correlation theory and recent developments in periodic integral evaluation^34–36 and correlation-consistent Gaussian basis sets,³⁷ these calculations were relatively affordable: they were generated from scratch in a few days, requiring about 18k CPU hours. Because of this low cost, we now extend our work to calculate a smooth, two-dimensional PES for CO on MgO, from which we calculate the vibrational frequency shift of CO stretching upon adsorption to be +14.7 cm⁻¹, which is in nearly quantitative agreement with the infrared experimental result³⁸ of +14.3 cm⁻¹. Using the CCSD(T) PES, we assess the accuracy of lower-level wavefunction theories and various DFT protocols that combine commonly used exchange–correlation functionals and dispersion corrections with respect to the calculated binding curve (Section III B) and vibrational frequency shift of CO stretching upon adsorption (Section III C). We find that none of the tested DFT protocols can achieve high accuracy in both the binding energy and the vibrational frequency. By contrast, MP2 reproduces the entire CCSD(T) PES within about 1 kJ mol⁻¹ and the vibrational frequency shift within about 1 cm⁻¹, while CCSD performs slightly worse with an error of about 4 kJ mol⁻¹ and 3 cm⁻¹, respectively. We conclude the work in Section IV with a discussion on potential future developments.

II. Computational details

The system being studied is shown in Fig. 1A. A CO molecule is adsorbed with its C end on a five-fold coordinated Mg site of the pristine MgO (001) surface. The physisorption results in only slight perturbations to the surface and negligible change of the CO bond length. We model the surface using a two-layer MgO slab, which is sufficient to converge the calculated adsorption energy to accuracy within 0.3 kJ mol⁻¹ (Table S5†). The equilibrium geometries for CO, MgO, and MgO + CO are obtained using DFT with the Perdew–Burke–Ernzerhof (PBE) functional⁴¹ and the D3 dispersion correction⁴² using Quantum Espresso.^43,44 Atoms in the bottom layer of the MgO slab are fixed during geometry optimization.

Fig. 1 (A) Atomic structures of a CO molecule physically adsorbed on a MgO (001) surface. The Mg–C distance is highlighted. (B) Convergence of the HF interaction energy with surface and basis set sizes. (C) Convergence of the correlation part of the MP2 interaction energy with surface and basis set sizes. The TDL results estimated by extrapolating the surface size with A = 2 × 2 and 3 × 3 and with A = 3 × 3 and 4 × 4 are shown as dashed and solid lines, respectively. The CBS results from (DZ,TZ) and (TZ,QZ) extrapolation are also shown. (D) Isosurface plot of a representative localized occupied orbital on the CO molecule (obtained using the generalized Pipek–Mezey method^39,40) and the density of the corresponding unoccupied LNOs with tightening truncation threshold η (generated for a 3 × 3 surface using the TZ basis set). The number of unoccupied LNOs are listed at the bottom, which is only a fraction of the total number of unoccupied orbitals, 1083. (E) Convergence of the correlation part of the frozen-core LNO-CCSD/CCSD(T) interaction energy with the LNO truncation threshold η for a 2 × 2 surface with (DZ,TZ)-extrapolated CBS limit. Both the regular LNO-CCSD/CCSD(T) results (hollow circles) and those corrected using eqn (7) (filled circles) are shown. The corrected LNO-CCSD/CCSD(T) interaction energies are seen to converge much faster to the canonical CCSD/CCSD(T) results, as indicated by the green horizontal line (the shaded area indicates an error of ±0.5 kJ mol⁻¹). (F) Convergence of the correlation part of the final LNO-CCSD/CCSD(T) interaction energy with the LNO truncation threshold η. The (DZ,TZ) extrapolated CBS limit is employed. Results from two different schemes for the TDL extrapolation, A = 2 × 2 and 3 × 3 (filled circles) and A = 3 × 3 and 4 × 4 (hollow circles), agree well at convergence. The converged LNO-CCSD/CCSD(T) values are indicated by the solid horizontal line, with the shaded area indicating an error of ±0.5 kJ mol⁻¹.

We calculate the adsorption energy of CO on MgO as follows

E_ads = E_int + Δ_geom (1)

where E_int is the interaction energy between CO and MgO calculated with their respective geometries fixed to be those in the MgO + CO composite system, and Δ_geom is the geometry relaxation energy. Following ref. 17, we use our DFT-determined Δ_geom (1.1 kJ mol⁻¹ from PBE+D3, which is in good agreement with previous work^16,17,27) in eqn (1) and evaluate E_int using correlated wavefunction methods.

All wavefunction calculations are performed with the PySCF code.^45,46 The GTH-cc-pVXZ basis sets,³⁷ augmented by diffuse functions for the CO molecule and nearby surface atoms, are employed with the GTH pseudopotential optimized for HF.^47–49 All interaction energies are corrected for basis set superposition error. The two-electron integrals are treated by the range-separated density fitting algorithm,^34,35 and the integrable divergence of the HF exchange is treated using a Madelung constant correction.^50–52 The O(N⁶)/O(N⁷) cost scaling of canonical CCSD/CCSD(T) is avoided by a periodic extension of the local natural orbital (LNO) approximation, which we have recently applied to study the dissociation of water on the surface of Al₂O₃ and TiO₂.³² Our periodic implementation, detailed in ref. 32, closely follows the molecular LNO-CCSD/CCSD(T) originally developed by Kállay and co-workers.^53,54 We calculate the correlation energy of a supercell containing N electrons as a sum of contributions from all N localized occupied orbitals i,

Each local contribution E^LNO-CC_i,corr is evaluated independently in a truncated set of occupied and unoccupied LNOs that are optimized for local orbital i.^53,54 The accuracy of the LNO approximation can be systematically improved by adjusting a single parameter: the threshold η used to truncate the LNOs by their occupation numbers (following ref. 54, we use a threshold for unoccupied orbitals that is fixed to be 10 times smaller than that of the occupied orbitals; all reported thresholds henceforth are for unoccupied orbitals). A representative localized occupied orbital is shown in Fig. 1D, together with the density of the corresponding unoccupied LNOs, whose spatial extension increases as the truncation threshold tightens. As η → 0, the LNO-CCSD/CCSD(T) results converge to the result of canonical CCSD/CCSD(T) calculations. To expedite this convergence, we apply a standard MP2 correction evaluated at the same η.⁵³ This is equivalent to correcting full MP2 by a short-range CC correction evaluated within the LNO subspace, i.e.,

E^LNO-CC+ΔPT2_corr(η) = E^MP2_corr + ΔE^LNO-CC_corr(η) (3)

whereEqn (3) is also used to account for other effects such as core-electron correlation. More computational details can be found in the ESI.†

III. CO on MgO (001) surface

A. Equilibrium adsorption energy

We first validate our periodic LNO-CC methodology by calculating the equilibrium adsorption energy E_ads for a single CO molecule adsorbed on the MgO (001) surface. Reference values of E_ads can be derived from the desorption activation energy measured by temperature-programmed desorption experiments.^9–11 This was done earlier by Boese and Sauer¹⁶ who obtained a reference value of −20.6 ± 2.4 kJ mol⁻¹ based on the experimental work by Dohnálek et al.,¹¹ and more recently by Shi and co-workers¹⁷ who derived a modified value of −19.2 ± 1.0 kJ mol⁻¹ by averaging two experimental results.^9,11 A correct theoretical prediction of E_ads, however, necessitates the use of correlated wavefunction methods. Boese and Sauer¹⁶ employed cluster-based MP2-in-DFT embedding and found an E_ads of −20.9 ± 0.7 kJ mol⁻¹, and Shi and co-workers¹⁷ calculated E_ads using embedded-cluster LNO-CCSD(T) to be −19.2 ± 1.0 kJ mol⁻¹. Clearly, both values are in good agreement with experimental estimates.

We revisit this problem to establish a protocol to achieve sub-kJ mol⁻¹ accuracy for the CCSD(T) interaction energy. We use a fixed d(Mg–C) = 2.460 Å for this purpose, which also facilitates comparison with previous works.^16,17,28,55 Using eqn (3), the CC interaction energy can be calculated as

where E^MP2_int is the MP2 interaction energy and ΔE^LNO-CC(FC)_int,corr is the LNO-CC correction (4); we evaluate the latter with the [2s²2p⁶] semicore electrons of Mg being frozen. Both terms in eqn (5) need to be converged to the thermodynamic limit (TDL) with respect to the surface size A and to the complete basis set (CBS) limit with the basis set size X, while ΔE^LNO-CC(FC)_int,corr also needs to be converged to the full CC limit with the LNO truncation threshold η.

The MP2 interaction energy can be further decomposed as a sum of the HF part (E^HF_int) and the MP2 correlation part (E^MP2_int,corr). Fig. 1B shows the convergence of E^HF_int(A,X). We see that using a 4 × 4 surface and a QZ basis set essentially attains both the TDL and the CBS limit, giving E^HF_int(A_∞,X_∞) ≈ 1.8 ± 0.1 kJ mol⁻¹, where half the difference between TZ and QZ results is taken as an estimate of the remaining basis set incompleteness error. Fig. 1C shows the convergence of E^MP2_int,corr(A,X), which is slower than that of E^HF_int(A,X) in both parameters. However, reliable extrapolations to both limits can be performed based on the asymptotic behaviors

E^MP2_int,corr(A,X) ≈ E^MP2_int,corr(A_∞,X) + c₁A⁻¹ (6a)

E^MP2_int,corr(A,X) ≈ E^MP2_int,corr(A,X_∞) + c₂X⁻³ (6b)

A good estimate of E^MP2_int,corr(A_∞,X_∞) is obtained by extrapolating the surface size with A = 3 × 3 and 4 × 4 and the basis set size with TZ (X = 3) and QZ (X = 4). This gives E^MP2_int,corr(A_∞,X_∞) ≈ −21.7 ± 0.3 kJ mol⁻¹, where the uncertainty arises from the remaining basis set incompleteness error and the pseudopotential error (Table S2†). Notably, repeating the calculations without correlating Mg's semicore electrons gives E^MP2(FC)_int,corr(A_∞,X_∞) ≈ −19.5 kJ mol⁻¹, revealing a sizable contribution of about −2.2 kJ mol⁻¹ from the semicore electrons. Summing the HF and MP2 correlation contributions, we arrive at our final estimate of E^MP2_int, which is −19.9 ± 0.3 kJ mol⁻¹.

Transitioning to CC, we first show in Fig. 1E the η-convergence of E^{LNO-CC+ΔPT2(FC)}_int,corr(A,X,η), extrapolated to the CBS limit based on DZ and TZ results for a 2 × 2 surface. The system is small enough (about 650 orbitals in the TZ basis set) to allow performing canonical CC calculations, whose results are also shown in Fig. 1E for comparison. We see that the error of LNO-CCSD(T) (orange hollow circles) drops fast to 0.5 kJ mol⁻¹ with a modest threshold of 10⁻⁶, but achieving a similar accuracy in LNO-CCSD (blue hollow circles) requires a tighter threshold of 10⁻⁷. While LNO-CC calculations with a tight threshold of 10⁻⁷ are feasible for larger systems, a more efficient alternative is to correct the LNO-CC results by the difference from an inexpensive canonical CC calculation using DZ basis set and 2 × 2 surface, which we denoted as

Δ^CC(FC)(η) = E^CC(FC)_int,corr(2 × 2,DZ) − E^LNO-CC(FC)_int,corr(2 × 2,DZ,η) (7)

The corrected LNO-CC results are shown as filled circles in Fig. 1E and seen to converge much faster than the uncorrected ones. In Fig. S2,† we confirmed that the same correction (7) is also transferrable for correcting larger surfaces.

To obtain our final estimate of the CC interaction energy using eqn (5), we extrapolate ΔE^LNO-CC(FC)_int,corr(A,X,η) to the TDL and the CBS limit using A = 2 × 2 and 3 × 3 and the DZ and TZ basis sets, and add this number to the converged E^MP2_int determined above. The resulting η-dependent LNO-CC interaction energy is then converged with respect to η (Fig. 1F). Our best estimate of the CCSD and CCSD(T) interaction energy, obtained with η = 10⁻⁷, is −16.8 ± 0.5 kJ mol⁻¹ and −21.0 ± 0.5 kJ mol⁻¹, where the error bar combines that of E^MP2_int (0.3 kJ mol⁻¹), of Δ^CC(FC) (0.2 kJ mol⁻¹), and of the LNO truncation error which is estimated to be 0.3 kJ mol⁻¹ using the difference between the two LNO-CCSD calculations with the tightest thresholds shown in Fig. 1F.

Our final numbers for the HF, MP2, CCSD, and CCSD(T) E_ads are summarized in Table 1, where they are compared to both the experimental reference and other computational results from literature. Compared to ref. 17, our periodic Gaussian-based approach shows sub-kJ mol⁻¹ agreement with the embedded-cluster approach for both the MP2 and CCSD(T) adsorption energy, and a slightly worse agreement of 2 kJ mol⁻¹ with the periodic plane-wave CCSD and CCSD(T), likely due to the larger error bar associated with the latter. Given the use of slightly different geometries and different strategies to attain convergence along the necessary theoretical axes, the level of agreement achieved between our work and ref. 17 is remarkable. Our final CCSD(T) adsorption energy, −20.0 ± 0.5 kJ mol⁻¹, also agrees well with corrected^16,17 temperature-programmed desorption experimental results.^9,11 By contrast, insufficient convergence along one or more of the theoretical axes can lead to significant errors in the calculated interaction energy,^28,55 which has been discussed thoroughly in previous work.¹⁷

Table 1 CO adsorption energy E_ads on the MgO (001) surface with a Mg–C distance of 2.460 Å obtained through various wavefunction methods compared to recent results in literature. The error bars of our numbers are explained in the main text, while those of the results in ref. 17 have the uncertainty in Δ_relax removed

Method	Comput. details^a	E ads (kJ mol^−1)	Reference
a Gaussian-type orbital (GTO), plane wave (PW), frozen natural orbital (FNO), diffusion Monte Carlo (DMC). b Based on TPD experiments in ref. 11. c Based on TPD experiments in ref. 9 and 11.
HF	Periodic/GTO	+2.9 ± 0.1	This work
MP2	Periodic/GTO	−18.8 ± 0.3
CCSD	Periodic/GTO/LNO	−15.7 ± 0.5
CCSD(T)	Periodic/GTO/LNO	−20.0 ± 0.5
MP2	Cluster/GTO	−18.5 ± 0.5	Ref. 17
CCSD(T)	Cluster/GTO/LNO	−19.2 ± 0.6
CCSD	Periodic/PW/FNO	−14.0 ± 2.1
CCSD(T)	Periodic/PW/FNO	−18.6 ± 2.1
DMC	Periodic	−18.1 ± 2.3
Experiments		−20.6 ± 2.4^b	Ref. 16
		−19.2 ± 1.0^c	Ref. 17

---

B. Full binding curve

We now employ the established protocol to calculate the full binding curve for CO adsorption on the MgO (001) surface. An accurate binding curve is important for understanding the adsorption/desorption dynamics,⁵⁶ which can be directly related to molecular beam experiments.¹ Here, we focus on calculating the variation of E_int as a function of d(Mg–C) (with atoms in the surface kept frozen). The results from the wavefunction methods discussed in Section III A are shown in Fig. 2A. We see that our protocol generates smooth energy curves for all methods including CCSD(T). This is remarkable given the small energy scale involved here (about 2 kJ mol⁻¹ near equilibrium; see the inset of Fig. 2A) and highlights the extremely high fidelity of our periodic LNO-CCSD(T) methodology. Performance-wise, CCSD and MP2 underestimate the CCSD(T) binding energy by about 4 and 1 kJ mol⁻¹ for a wide range of d(Mg–C) near the equilibrium geometry, which is consistent with the trend discussed in Section III A for the equilibrium E_ads. All three methods smoothly approach the correct dissociation limit. Interestingly, the HF binding curve does show a shallow well of about −3 kJ mol⁻¹ at a large d(Mg–C) of nearly 3 Å.

Fig. 2 The full binding curve for CO adsorption on the MgO (001) surface calculated using selected (A) wavefunction and (B and C) DFT methods. The energy data are denoted by the markers. The solid curves are obtained by fitting the data to the Morse potential. The inset in panel (A) shows a zoomed-in view for d(Mg–C) ranging from 2.25 to 2.75 Å.

Armed with CCSD(T) results, we can now assess the performance of different DFT methods. Fig. 2B presents the binding curves obtained using the PBE functional with five popular dispersion correction models: the exchange-hole dipole moment (XDM) model,⁵⁷ the Tkatchenko–Scheffler (TS) model,⁵⁸ the many-body dispersion (MBD) model,⁵⁹ and two different generations of Grimme's DFT-D model (D2 (ref. 60) and D3 (ref. 42)). We see that all five DFT protocols significantly overestimate the equilibrium binding energy by about 6 kJ mol⁻¹ (D2) to 16 kJ mol⁻¹ (TS) and predict the equilibrium Mg–C separation to be too small by 5–10 pm. To investigate the effect of different exchange–correlation functionals, in Fig. 2C we present binding curves calculated using PBE and its two extensions, revPBE⁶¹ and PBE0,⁶² all corrected by the D3 dispersion correction. It is evident that both revPBE and PBE0 ameliorate PBE's overbinding behavior, but the improvement is way too small. As seen in Fig. 2C, we find that revPBE+D4 (ref. 63) gives the best performance near the equilibrium geometry, in agreement with ref. 17, but the D4 dispersion energy does not decay properly to zero in the dissociation limit.

C. Frequency shift of CO stretching

We now employ the established protocol to calculate the vibrational frequency of CO on the MgO (001) surface. Experimentally, this has been measured by infrared spectroscopy³⁸ to be ν_CO(ads) = 2157.5 cm⁻¹. Compared to the gas-phase frequency of ν_CO(g) = 2143.2 cm⁻¹, this corresponds to a positive vibrational frequency shift Δν_CO = ν_CO(ads) − ν_CO(g) = 14.3 cm⁻¹. Such frequency shifts are commonly used for probing both the structural and electronic information of a surface.^10,64 Accurately resolving such a small frequency shift serves as a stringent test for electronic structure theory.

A minimal theoretical description necessitates calculating a two-dimensional potential energy surface for the CO bond length d(C–O) and the surface–CO distance, chosen here to be D(Mg–CO), the distance between the CO center of mass and the surface Mg atom (the surface is kept frozen at its equilibrium geometry following previous work⁶⁵). We calculate this two-dimensional PES using various electronic structure methods by performing a series of single-point calculations over a total of 130 different geometries and fitting the resulting energy data using a sixth-order polynomial. The fitting error is found to be less than 0.5 kJ mol⁻¹ for all electronic structure methods discussed below. We then extract ν_CO(ads) in the harmonic approximation from the second derivatives of the fitted PES. We calculate ν_CO(g) in a similar manner from a one-dimensional PES of d(C–O) for a free CO molecule. See the ESI for more details.†

The PES for the surface adsorbed CO obtained using LNO-CCSD(T) is shown in Fig. 3A. The contour of the PES displays a characteristic elliptical shape for the coupling of a strong vibrational mode to a weak one. From the PES of the adsorbed CO and that of a free CO (not shown), we extract ν^CCSD(T)_CO(ads) = 2173.9 cm⁻¹ and ν^CCSD(T)_CO(g) = 2159.2 cm⁻¹, both of which are about 16 cm⁻¹ higher than the experimentally observed frequencies. As a result, the frequency shift predicted by CCSD(T), Δν^CCSD(T)_CO = 14.7 cm⁻¹, is in quantitative agreement with experiment, with an error less than 0.5 cm⁻¹. From a practical point of view, the harmonic CCSD(T) overestimation of vibrational frequencies of small molecules is well known and usually corrected by applying an empirical scaling factor.^66,67 The nearly perfect agreement between the computational and experimental Δν_CO thus indicates a highly transferrable scaling factor for CO in gas and condensed phases. Indeed, applying the scaling factor 0.9926 determined by our gas-phase calculation and the experiment, we obtain a scaled vibrational frequency for adsorbed CO, ν^{CCSD(T)-scaled}_CO(ads) = 2157.8 cm⁻¹, which agrees nearly perfectly with the experimental frequency.

Fig. 3 (A) CCSD(T) potential energy surface for the CO stretching mode coupled with its adsorption onto the MgO (001) surface, calculated using LNO-CCSD(T) with the protocol established in Section III A. d(C–O) is the CO bond length and D(Mg–CO) is the distance between the CO center of mass and the surface Mg atom. The number labelling each contour line is the energy relative to the energy minimum (denoted by a green star) in kJ mol⁻¹. (B) The CO stretching frequency for a free CO molecule (red) and a CO molecule adsorbed on the MgO surface (blue) predicted by different methods. The scaled CCSD(T) frequencies, using a scaling factor based on our gas-phase calculation and experiment, as discussed in the main text, are also shown. The experimental values are shown as dashed horizontal lines of corresponding color. For each computational method, a gray arrow indicates the direction of the predicted frequency shift upon adsorption. (C) Shift of the CO stretching frequency upon adsorption calculated from the frequency data shown in panel (B). Solid blue and hatched red bars indicate blue and red shift, respectively. The experimental reference is shown as the green shaded area.

The CCSD(T) frequency and frequency shift discussed above are displayed in Fig. 3B and C. In the same plots, we also show the absolute frequency and frequency shift predicted by several different electronic structure methods. MP2 slightly underestimates ν_CO(ads) by about 20 cm⁻¹, while CCSD overestimates it by about 90 cm⁻¹. However, both methods also predict a similar trend for ν_CO(g), making their frequency shift, 15.8 cm⁻¹ from MP2 and 17.8 cm⁻¹ from CCSD, in reasonable agreement with the experiment.

In contrast, DFT with semilocal functionals, such as PBE and revPBE, predicts a large negative Δν_CO ranging from −15 to −20 cm⁻¹, regardless of the dispersion correction being used (Fig. 3C). This qualitative error is caused by the imbalanced error in the adsorbed and gas-phase calculations: while both ν_CO(ads) and ν_CO(g) are underestimated, the former is underestimated to a much larger extent (Fig. 3B), likely due to a larger delocalization error in condensed-phase systems. Improvement is seen by using hybrid functionals such as PBE0 and B3LYP, which suffer less from the delocalization error due to partial inclusion of the HF exchange energy.⁶⁸ However, the predicted Δν_CO, 3.7 cm⁻¹ from PBE0+D3(0) and 0.3 cm⁻¹ from B3LYP + D3(0), are still too small compared to the experiment. Despite having qualitatively correct shapes, the PESs obtained using these DFT methods exhibit clear quantitative deviations from the CCSD(T) PES (Fig. S5†). This highlights the challenge for achieving spectroscopic accuracy for vibrations of adsorbed molecules with commonly used DFT methods.

IV. Conclusion

To summarize, we computationally studied the adsorption of a single CO molecule on the MgO (001) surface using the “gold-standard” method of quantum chemistry, CCSD(T), made possible by our recently developed periodic LNO-CCSD(T) method. After validating our method using the equilibrium adsorption energy, we leverage the low computational cost and high fidelity of LNO-CCSD(T) to calculate a smooth, two-dimensional PES for CO on MgO, from which we extract the vibrational frequencies and frequency shift of CO stretching upon adsorption, finding quantitative agreement with infrared spectroscopy. Using the CCSD(T) PES, we assess the accuracy of lower-level wavefunction theories and various DFT protocols that combine commonly used exchange–correlation functionals and dispersion corrections. We find that none of the tested DFT protocols can achieve high accuracy in both the binding energy and the vibrational frequency. By contrast, MP2 reproduces the entire CCSD(T) PES quite accurately, while CCSD performs slightly worse but still qualitatively correctly.

The ability to generate a smooth PES with high accuracy is a longstanding challenge in computational surface science for benchmarking and improving lower-level theories like DFT^{1,14,56,69,70} and more recently, for training high-quality machine learning potentials that enable accurate ab initio molecular dynamics study of surface reactions.^18,20,21 For insulating and semiconducting surfaces with weak electron correlation, the smooth and accurate PESs attainable through our periodic local CCSD(T) method are thus very promising to fulfill both goals. But for chemistry on the surface of metals or strongly correlated solids, for which CCSD(T) is inapplicable,^71–76 new methods are sorely needed.

Conflicts of interest

The authors declare no competing conflicts of interest.

Acknowledgements

This work was supported by the National Science Foundation under Grant No. CHE-1848369 and by the Columbia Center for Computational Electrochemistry. We acknowledge computing resources from Columbia University's Shared Research Computing Facility project, which is supported by NIH Research Facility Improvement Grant 1G20RR030893-01, and associated funds from the New York State Empire State Development, Division of Science, Technology and Innovation (NYSTAR) Contract C090171, both awarded April 15, 2010.

References

1. G.-J. Kroes, Phys. Chem. Chem. Phys., 2021, 23, 8962.
2. T.-N. Ye, S.-W. Park, Y. Lu, J. Li, J. Wu, M. Sasase, M. Kitano and H. Hosono, J. Am. Chem. Soc., 2021, 143, 12857.
3. K. Sumida, D. L. Rogow, J. A. Mason, T. M. McDonald, E. D. Bloch, Z. R. Herm, T.-H. Bae and J. R. Long, Chem. Rev., 2012, 112, 724.
4. O. T. Qazvini, R. Babarao and S. G. Telfer, Nat. Commun., 2021, 12, 197.
5. M. Chang, T. Yan, Y. Wei, J.-X. Wang, D. Liu and J.-F. Chen, Chem. Mater., 2022, 34, 9134.
6. A. W. Kleyn, Chem. Soc. Rev., 2003, 32, 87.
7. Z.-T. Wang, Y.-G. Wang, R. Mu, Y. Yoon, A. Dahal, G. K. Schenter, V.-A. Glezakou, R. Rousseau, I. Lyubinetsky and Z. Dohnálek, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 1801.
8. J. L. Falconer and J. A. Schwarz, Catal. Rev. Sci. Eng., 1983, 25, 141.
9. R. Wichtendahl, M. Rodriguez-Rodrigo, U. Härtel, H. Kuhlenbeck and H.-J. Freund, Phys. Status Solidi A, 1999, 173, 93.
10. M. Sterrer, T. Risse and H.-J. Freund, Appl. Catal., A, 2006, 307, 58.
11. Z. Dohnálek, G. A. Kimmel, S. A. Joyce, P. Ayotte, R. S. Smith and B. D. Kay, J. Phys. Chem. B, 2001, 105, 3747.
12. L. Zhang, C. Tian, G. A. Waychunas and Y. R. Shen, J. Am. Chem. Soc., 2008, 130, 7686.
13. F. M. Geiger, Annu. Rev. Phys. Chem., 2009, 60, 61.
14. L. Schimka, J. Harl, A. Stroppa, A. Grüneis, M. Marsman, F. Mittendorfer and G. Kresse, Nat. Mater., 2010, 9, 741.
15. J. Sauer, Acc. Chem. Res., 2019, 52, 3502.
16. A. D. Boese and J. Sauer, Phys. Chem. Chem. Phys., 2013, 15, 16481.
17. B. X. Shi, A. Zen, V. Kapil, P. R. Nagy, A. Grüneis and A. Michaelides, J. Am. Chem. Soc., 2023, 145, 25372.
18. L. Zhang, J. Han, H. Wang, R. Car and W. E, Phys. Rev. Lett., 2018, 120, 143001.
19. O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schütt, A. Tkatchenko and K.-R. Müller, Chem. Rev., 2021, 121, 10142.
20. B. Wen, M. F. C. Andrade, L.-M. Liu and A. Selloni, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2212250120.
21. Z. Zeng, F. Wodaczek, K. Liu, F. Stein, J. Hutter, J. Chen and B. Cheng, Nat. Commun., 2023, 14, 6131.
22. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864.
23. A. J. R. Hensley, K. Ghale, C. Rieg, T. Dang, E. Anderst, F. Studt, C. T. Campbell, J.-S. McEwen and Y. Xu, J. Phys. Chem. C, 2017, 121, 4937.
24. F. R. Rehak, G. Piccini, M. Alessio and J. Sauer, Phys. Chem. Chem. Phys., 2020, 22, 7577.
25. E. W. F. Smeets and G.-J. Kroes, J. Phys. Chem. C, 2021, 125, 8993.
26. G. Meng, R. Yin, X. Zhou and B. Jiang, J. Phys. Chem. C, 2021, 125, 24958.
27. M. Alessio, D. Usvyat and J. Sauer, J. Chem. Theory Comput., 2019, 15, 1329.
28. A. Mitra, M. R. Hermes, M. Cho, V. Agarawal and L. Gagliardi, J. Phys. Chem. Lett., 2022, 13, 7483.
29. K. Raghavachari, G. W. Trucks, J. A. Pople and M. Head-Gordon, Chem. Phys. Lett., 1989, 157, 479.
30. B. X. Shi, V. Kapil, A. Zen, J. Chen, A. Alavi and A. Michaelides, J. Chem. Phys., 2022, 156, 124704.
31. H.-Z. Ye and T. C. Berkelbach, arXiv, 2023, preprint, arXiv:2309.14651v2,  DOI:10.48550/arXiv.2309.14651.
32. H.-Z. Ye and T. C. Berkelbach, arXiv, 2023, preprint, arXiv:2309.14640,  DOI:10.48550/arXiv.2309.14640.
33. H.-Z. Ye and T. C. Berkelbach, arXiv, 2024, preprint, arXiv:2407.11258,  DOI:10.48550/arXiv.2407.11258.
34. H.-Z. Ye and T. C. Berkelbach, J. Chem. Phys., 2021, 154, 131104.
35. H.-Z. Ye and T. C. Berkelbach, J. Chem. Phys., 2021, 155, 124106.
36. S. J. Bintrim, T. C. Berkelbach and H.-Z. Ye, J. Chem. Theory Comput., 2022, 18, 5374.
37. H.-Z. Ye and T. C. Berkelbach, J. Chem. Theory Comput., 2022, 18, 1595.
38. G. Spoto, E. Gribov, A. Damin, G. Ricchiardi and A. Zecchina, Surf. Sci., 2003, 540, L605.
39. J. Pipek and P. G. Mezey, J. Chem. Phys., 1989, 90, 4916.
40. S. Lehtola and H. Jónsson, J. Chem. Theory Comput., 2014, 10, 642.
41. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
42. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
43. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
44. P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. Otero-de-la-Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu and S. Baroni, J. Phys.: Condens. Matter, 2017, 29, 465901.
45. Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters and G. K.-L. Chan, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1340.
46. Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov and G. K.-L. Chan, J. Chem. Phys., 2020, 153, 024109.
47. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703.
48. C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 3641.
49. J. Hutter, New optimization of GTH pseudopotentials for PBE, SCAN, PBE0 functionals. GTH pseudopotentials for Hartree-Fock. NLCC pseudopotentials for PBE, 2019, https://github.com/juerghutter/GTH.
50. J. Paier, M. Marsman, K. Hummer, G. Kresse, I. C. Gerber and J. G. Ángyán, J. Chem. Phys., 2006, 124, 154709.
51. P. Broqvist, A. Alkauskas and A. Pasquarello, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 085114.
52. R. Sundararaman and T. A. Arias, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 165122.
53. Z. Rolik and M. Kállay, J. Chem. Phys., 2011, 135, 104111.
54. Z. Rolik, L. Szegedy, I. Ladjánszki, B. Ladóczki and M. Kállay, J. Chem. Phys., 2013, 139, 094105.
55. A. Mazheika and S. V. Levchenko, J. Phys. Chem. C, 2016, 120, 26934.
56. E. W. F. Smeets, J. Voss and G.-J. Kroes, J. Phys. Chem. A, 2019, 123, 5395.
57. A. D. Becke and E. R. Johnson, J. Chem. Phys., 2005, 122, 154104.
58. A. Tkatchenko and M. Scheffler, Phys. Rev. Lett., 2009, 102, 073005.
59. A. Tkatchenko, R. A. DiStasio, R. Car and M. Scheffler, Phys. Rev. Lett., 2012, 108, 236402.
60. S. Grimme, J. Comput. Chem., 2006, 27, 1787.
61. Y. Zhang and W. Yang, Phys. Rev. Lett., 1998, 80, 890.
62. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158.
63. E. Caldeweyher, C. Bannwarth and S. Grimme, J. Chem. Phys., 2017, 147, 034112.
64. J. B. Benziger and L. Larson, J. Catal., 1982, 77, 550.
65. N. He, M. Huang and F. A. Evangelista, J. Phys. Chem. A, 2023, 127, 1975.
66. A. P. Scott and L. Radom, J. Phys. Chem., 1996, 100, 16502.
67. D. P. Tew, W. Klopper, M. Heckert and J. Gauss, J. Phys. Chem. A, 2007, 111, 11242.
68. N. Mardirossian and M. Head-Gordon, Mol. Phys., 2017, 115, 2315.
69. J. M. Boereboom, M. Wijzenbroek, M. F. Somers and G. J. Kroes, J. Chem. Phys., 2013, 139, 244707.
70. T. Tchakoua, E. W. F. Smeets, M. Somers and G.-J. Kroes, J. Phys. Chem. C, 2019, 123, 20420.
71. M. Motta, C. Genovese, F. Ma, Z.-H. Cui, R. Sawaya, G. K.-L. Chan, N. Chepiga, P. Helms, C. Jiménez-Hoyos, A. J. Millis, U. Ray, E. Ronca, H. Shi, S. Sorella, E. M. Stoudenmire, S. R. White and S. Zhang, (Simons Collaboration on the Many-Electron Problem), Phys. Rev. X, 2020, 10, 031058.
72. M. S. Church and B. M. Rubenstein, J. Chem. Phys., 2021, 154, 184103.
73. J. Lee, H. Q. Pham and D. R. Reichman, J. Chem. Theory Comput., 2022, 18, 7024.
74. T. N. Mihm, T. Schäfer, S. K. Ramadugu, L. Weiler, A. Grüneis and J. J. Shepherd, Nat. Comput. Sci., 2021, 1, 801.
75. V. A. Neufeld and T. C. Berkelbach, Phys. Rev. Lett., 2023, 131, 186402.
76. N. Masios, A. Irmler, T. Schäfer and A. Grüneis, Phys. Rev. Lett., 2023, 131, 186401.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4fd00041b

---