Chemically Accurate Adsorption Energies: CO and H2O on the MgO

Article Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

pubs.acs.org/JCTC

Chemically Accurate Adsorption Energies: CO and H2O on the MgO(001) Surface Maristella Alessio, Denis Usvyat,* and Joachim Sauer* Institut für Chemie, Humboldt-Universität zu Berlin, Unter den Linden 6, 10099 Berlin, Germany

J. Chem. Theory Comput. Downloaded from pubs.acs.org by IOWA STATE UNIV on 01/05/19. For personal use only.

S Supporting Information *

ABSTRACT: Hybrid MP2:DFT-D structure optimizations are performed at BSSE-free CBS-extrapolated potential energy surfaces for molecule−oxide surface interactions (BSSE, basis set superposition error; CBS, complete basis set limit). Subsequently single point MP2 calculations are performed to estimate the effects of increasing the basis set size in the CBS extrapolation and increasing the cluster model size. The resulting estimates of the periodic MP2 limit agree within 1 kJ/mol with Local MP2 calculations using periodic boundary conditions. Single point CCSD(T) calculations are performed to determine ΔCC = CCSD(T) − MP2 energy differences. The final hybrid MP2:DFT-D+ΔCC estimate for CO on the MgO(001) surface at low coverage, −21.2 ± 0.5 kJ/mol, is in close agreement with the reference energy derived from temperature-programmed desorption experiments, −20.6 ± 2.4 kJ/mol. For H2O on MgO(001), at limiting zero coverage, we predict an adsorption energy of −53.7 ± 4.2 kJ/mol which falls in the range of values, −55.8 ± 12.2 kJ/mol, derived from a high coverage low energy electron diffraction experiments and estimated lateral interactions. (T)) “single point” calculations are carried out for smaller cluster models and basis sets. CCSD(T) correlation effects are then estimated as ΔCC = CCSD(T) − MP2 differences. With our hybrid MP2:DFT-D energies, we are aiming at approximating the MP2 limit for the periodic structure, and with the hybrid MP2:(DFT-D)+ΔCC energies, we are aiming at the CCSD(T) limit for the periodic structure. The problem is that neither the exact CCSD(T) nor the exact MP2 results are known for periodic systems with cell sizes required to describe molecule−surface interactions. Within our hybrid MP2:DFT‑D approach, the periodic MP2 limit can only be estimated by performing single point calculations on a series of cluster models of increasing size.2,4 Here, we perform single point local MP2 (LMP2) calculations with periodic boundary conditions using the CRYSCOR code.11 We use the same computational settings as for the hybrid MP2:DFT-D calculations, which allows a judgment on how accurately the latter approximates periodic MP2 results, and we find that the differences are well within chemical accuracy limits (±4.2 kJ/ mol). In the absence of CCSD(T) results, comparison with experiments is the way to judge the performance of our hybrid MP2:(DFT-D)+ΔCC method. Making this comparison, first, we need to be aware of the accuracy limits of CCSD(T) itself, for which ±2 kJ/mol may be a fair judgment.12,13 Moreover, in the area of molecule−surface interactions, this comparison

1. INTRODUCTION A hybrid high-level QM:low-level QM approach (QM, quantum mechanics) has been proposed to overcome the typical dilemma in computational material chemistry between the system size we need to consider and the accuracy we aim at. It combines the advantages of two different QM approaches: the more accurate and systematically improvable wave function theory applied to a finite-size cluster model of the active site with the more approximate but computationally more efficient density functional theory (DFT) for the entire system applying periodic boundary conditions.1,2 Specifically, we employ second order Møller−Plesset perturbation theory (MP2) for the calculations on the cluster models and DFT with some account of dispersion (DFT-D) for the entire periodic system.3,4 This hybrid high-level cluster:low-level periodic approach is in the spirit of Morokuma’s Integrated MO:MO (IMOMO) approach5 which later became ONIOM.6 The method we apply also uses the “subtractive scheme” for calculating energies and forces of a smaller cluster embedded in a larger system but goes beyond ONIOM as periodic boundary conditions can be applied at the low level.7,8 For structure optimizations at the hybrid MP2:DFT-D potential energy surface (PES), we employ the MonaLisa code,9 which allows to make use of energies and forces that are extrapolated to the complete basis set (CBS) limit,10 and also offers the option to apply counterpoise corrections (CPC) for the basis set superposition error (BSSE). At the hybrid MP2:DFT-D equilibrium structures, coupled-cluster including single, double, and (perturbative) triple substitutions (CCSD© XXXX American Chemical Society

Received: November 3, 2018

A

DOI: 10.1021/acs.jctc.8b01122 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation requires a careful analysis of the experiments performed as demonstrated for the CO adsorption on MgO(001) terrace sites4 for which adsorption energies with known error bars are available from TPD experiments.14 Finally, the calculated energies must be converted in enthalpy values for the conditions of the measurement using molecular statistics. For judging different computational approaches, it is more convenient to do the opposite and convert experimental enthalpies, ΔHT, in “experimentally derived” reference energies, ΔEexp, taking zero-point vibrational energies, ΔEZPV, and thermal enthalpy changes, ΔTΔH, into account.4,15 ΔEexp = ΔHT − ΔT ΔH − ΔEZPV

E HL:LL(pbc) = E LL(pbc) − E LL(C) + E HL(C) = E LL(pbc) + ΔHL (C) = E HL(C) + ΔLR (pbc,C)

(2)

where pbc indicates that calculations are performed on periodic slab models, and C stands for finite-size cluster models. The hybrid QM:QM energy is thus obtained as the low-level periodic energy, ELL(pbc), to which a high-level correction ΔHL (C) = E HL(C) − E LL(C)

(3)

is added. The latter is the difference between the high-level and low-level energies evaluated for the same cluster model. In an equivalent way, the hybrid QM:QM energy can be interpreted as the high-level energy of the cluster model, EHL(C), to which the long-range correction

(1)

The limited accuracy of the calculated vibrational energies creates uncertainties for zero-point vibrational energy and thermal energy changes that add to the uncertainty limits of the experimental values. This hybrid MP2:(DFT-D)+ΔCC method was shown to yield results that agree within chemical accuracy limits (±4.2 kJ/mol) with experimentally derived energies for a variety of experimentally well-defined molecule−surface interactions such as adsorption on flat oxide surfaces,4,10 in metal−organic frameworks (MOFs),16−19 and in zeolites,15 as well as for energy barriers in zeolite catalysis.20 These are all systems for which CCSD(T) is expected to yield chemically accurate results. For other types of surfaces such as transition metal oxides or metals, the hybrid high-level QM:low-level QM approach would still be applicable, but the high-level method needs to be different and the accuracy that can be reached would be defined by the high-level method. We are looking at two different systems: (i) We reconsider CO adsorption on MgO(001) terrace sites which has been studied before with the hybrid MP2:(DFT-D)+ΔCC method.4 The present study differs in several respects. First, instead of the QMPOT code,1,8 the MonaLisa code9 is employed, and hybrid MP2:DFT-D optimizations are performed at the BSSE-free CBS extrapolated PES. Second, the present study uses both PBE+D2 and B3LYP+D2 as low-level methods. Comparison with ref 4 shows that BSSE has a non-negligible effect on the predicted molecule−surface distance and that the hybrid MP2:DFT-D approach is stable with respect to variations of the low-level DFT-D method. Third, we have been also able to extend the basis set for CCSD(T) calculations from double-ζ to triple-ζ quality resulting in a significant change of the ΔCC difference. (ii) Having revalidated our hybrid MP2:(DFT-D)+ΔCC method, we employ it for predicting the adsorption structure and energy of H2O on the MgO(001) surface at the limiting zero coverage as a reference value. Due to attractive lateral interactions, experiments yield information only for monolayer coverage in a temperature range of 100−200 K,21,22 but need to estimate lateral interactions from studying condensed water phases.

ΔLR (pbc,C) = E LL(pbc) − E LL(C)

(4)

computed at the low-level is added. To perform structure optimizations on the hybrid QM:QM potential energy surface(PES), the forces on the nuclei are expressed by an analogous subtractive relation.7,8 While the presented implementation mechanically embeds the high-level QM cluster into the low-level QM lattice and considers the interactions between the two regions at the lowlevel only, over the years also other hybrid cluster QM:periodic QM schemes have been proposed in the literature, e.g., refs 24−28. Further, electronic embedding schemes have been proposed29,30 and developed into efficient tools29 that can be applied to large chemical systems.31 2.2. Extrapolation to the Complete Basis Set Limit and Counterpoise Correction. Calculations for molecule− surface interactions based on atom-centered basis functions, e.g., Gaussian-type orbitals (GTOs), suffer from the basis set superposition error (BSSE) 32 and from the basis set incompleteness error (BSIE).33 To improve the accuracy of the high-level QM term as part of the hybrid QM:QM method, the high-level QM energies and forces are calculated with an extrapolation scheme that uses two basis sets of consecutive cardinal numbers (X-1 and X).10 For each basis set, the energies and gradient can also be corrected for BSSE using the counterpoise correction scheme. This defines multilevel QM/ (CPC-)CBS:QM methods and allows structure optimizations on BSIE-free and, if the user decides so, also BSSE-free hybrid potential energy surfaces. For decades, it has been discussed whether BSSE-corrected binding energies are superior to uncorrected ones and if BSSE should be subtracted fully or only a fraction of it. If a series of correlation-consistent basis sets of increasing cardinal numbers is considered, both the CP-corrected and uncorrected values will converge to the same limit as shown by Dunning;33 however, convergence is faster with CP-corrected results. Hence, CBS extrapolation of both CP-corrected and uncorrected results will lead to the same limiting value, but only if performed with high cardinal quantum numbers, e.g., CBS(Q,5). Our results for the CO·Mg9O9 cluster in the Supporting Information (Tables S3 and S4) show a difference of 1.0 kJ/mol for the CBS(Q,5) extrapolations. From the average limit (−14.4 kJ/mol), the uncorrected CPC-CBS(D,T) result we use in our hybrid MP2:DFT-D scheme deviates by +0.5 kJ/mol only, whereas the uncorrected CBS(D,T) value deviates by −4.9 kJ/mol. The corresponding numbers for the H2O·Mg9O9 cluster in Tables S5 and S6 of the

2. METHODS 2.1. Subtractive Embedding Scheme. Hybrid QM:QM calculations are performed using the MonaLisa code9 which applies a “mechanical” embedding scheme5,23 under periodic boundary conditions (pbc).7 The total hybrid QM:QM energy is defined via the subtractive scheme:5,7,23 B


Article

Journal of Chemical Theory and Computation

All the contributions in eq 11 are computed considering the hybrid QM:QM equilibrium structures of S·M and S and the QM high-level equilibrium structure of the M molecule. Additionally, as seen in eq 6 for a general QM level, the hybrid QM:QM adsorption energy can be computed combining the hybrid adsorption energy, obtained at the equilibrium structure of S·M with the hybrid relaxation energy, ΔErelax,HL:LL(pbc).

Supporting Information also show a smaller deviation of the CPC-CBS(D,T) value (+1.9 kJ/mol) than of the CBS(D,T) value (−4.9 kJ/mol). Hence, for the systems studied here, we decided to optimize at the CPC-CBS(D,T) potential energy surfaces. 2.3. Adsorption Energies and Enthalpies. The adsorption energy, ΔE, is defined as ΔE = E(S·M) − E(S) − E(M)

ΔE HL:LL(pbc) = ΔE ̃HL:LL (pbc) − ΔErelax,HL:LL(pbc)

(5)

(12)

where S·M is the ad-molecule/MgO(001) system, S is the clean MgO(001) surface, and M is the molecule in the gas phase, all of them at their respective equilibrium structures. The adsorption energy can also be obtained combining single point calculations at the equilibrium structure of S·M with the relaxation energies of the MgO(001) surface and the adsorbed molecule ΔE = ΔE ̃ − ΔErelax

The final hybrid QM:QM adsorption energy, ΔEHL:LL(pbc), is not influenced by the chosen formalism, eqs 11 or 12, whereas the individual contributions ΔELL(pbc), ΔELL(C), ΔEHL(C), ΔLR, etc. can be significantly affected (Tables 2 and 6). 2.4. Computational Details. This work treats the low coverage case of CO (Θ = 1/8) and two different low coverage cases of H2O (Θ = 1/8 and 1/16) on the MgO(001) surface, see sections 3.1.1 and 3.2.1 for the description of the surface and cluster models adopted. Hybrid MP2:DFT-D calculations are performed using the MonaLisa program9 which supplies interfaces to the TURBOMOLE and VASP codes. Plane wave DFT calculations are performed with the 5.2.12 version of the VASP code.35 An energy cutoff of 600 and 400 eV is applied for the coverages Θ = 1/8 (both CO and H2O) and Θ = 1/16 (H2O only), respectively. The periodic calculations of the surface models are performed using a 2 × 2 × 1 and a 1 × 1 × 1 k-point grid each including the Γ point. The projector-augmented wave (PAW) scheme is used. For CO and H2O on MgO(001) at Θ = 1/8 coverage, standard PAW potentials for H, C, and O atoms are used, whereas the 2p shells of Mg2+ ions are considered as valence ones (Mg_pv). For H2O·MgO(001) at Θ = 1/16 coverage, standard PAW potentials are used for all the atoms. The PBE36 functional is employed. For CO·MgO(001), the B3LYP37,38 functional is used as an additional low-level method. Both PBE and B3LYP functionals are augmented by a damped atom-pair potential parametrized by Grimme39 (“D2”). For CO· MgO(001) and H2O·MgO(001) at Θ = 1/8 coverage, the dispersion contribution is computed by Ewald summation,40 and the Ne atom dispersion parameters (C6, atomic dispersion coefficient and R0, atomic vdW radii) are used for the Mg2+ ions as suggested by Tosoni and Sauer3 (considered as the default computational option). For H2O·MgO(001) at Θ = 1/ 16 coverage, the VASP standard cutoff radius of 30 Å for pair interactions is applied, and the Grimme Mg atom dispersion parameters are used for the Mg2+ ions (further referred to as “PBE+D2[Mg]”). For the Θ = 1/8 coverage of both CO and H2O, the VASP calculations of the gas phase molecules and adsorption cluster models are performed enlarging the lattice parameter from 8.48 Å (see sections 3.1.1 and 3.2.1) to 15 Å, while for the Θ = 1/16 coverage of H2O the cell volume is not altered. The MP2 calculations are carried out by using the resolution of identity (RI)-CC2 module41 available in version 6.5 of the TURBOMOLE code.42,43 The “1s frozen core” approximation is applied to all the atoms. A basis set, named “aug’-ccp(C)VXZ” in previous works,3,4,10 is used. The notation “aug’” indicates that augmentation functions (aug-cc-pVXZ)44 are employed for H and O atoms but not for the Mg2+ ions and “(C)” stands for core−valence basis sets (cc-pwCVXZ)45 adopted for the Mg2+ ions but not for other atoms. As usual, “X” is the basis set cardinal number. The CPC-CBS scheme is applied to MP2 energies and forces. The CBS limit is

(6)

The ΔẼ adsorption energy is computed as ΔE ̃ = E(S·M) − E(S//S ·M) − E(M//S·M)

(7)

where the // notation stands for “at the structure of” and means the corresponding calculations are performed at the equilibrium structure of S·M. The relaxation energy, ΔErelax, is a sum of MgO(001) surface and adsorbed molecule relaxation energies ΔErelax = ΔErelax (S) + ΔErelax (M)

(8)

Each contribution in eq 8 is defined in a generic form as follows: ΔErelax (i) = E(i) − E(i//S ·M)

(9)

with i = S, M. The partitioning of the adsorption energy into ΔẼ and ΔErelax allows one to account for BSSE by applying counterpoise corrections to eq 7. Similar to eq 5, the hybrid QM:QM adsorption energy, ΔEHL:LL, is defined as ΔE HL:LL = E HL:LL(S ·M) − E HL:LL(S) − E HL:LL(M) (10)

which requires hybrid QM:QM mechanical embedding calculations for the adsorbate−surface system S·M and the surface S, while for the molecule M a pure QM high-level calculation is performed. Each term of eq 10 is computed by applying the subtractive scheme of eq 2, which leads to an additional expression of the hybrid QM:QM adsorption energy ΔE HL:LL(pbc) = ΔE LL(pbc) − ΔE LL(C) + ΔE HL(C) = ΔE HL(C) + ΔLR (pbc,C)

(11)

Usually, the cluster for the high-level calculations can be chosen such that the long-range correction calculated at the low level is only a small fraction of the total adsorption energy. To get accurate results, the low-level cluster calculations should use the same code with the same basis set/core electron treatment as used for the low-level calculations of the periodic system, e.g., VASP with a plane wave basis set and PAW core description. Results from all-electron calculations with Gaussian basis sets and PAW/plane wave calculations will finally converge for large enough basis sets and high kinetic energy cut-offs, see, e.g., ref 34, but for the basis set size affordable for periodic systems, this is not always the case (Tables S1 and S2 in the Supporting Information). C


Article

Journal of Chemical Theory and Computation calculated by a two-basis set extrapolation scheme, with X-1 = D-zeta and X = T-zeta. At the hybrid MP2:DFT+D2 equilibrium structures of CO and H2O on MgO(001) at Θ = 1/8 coverage, for the smaller CO·Mg9O9 and H2O·Mg9O9 cluster models (see sections 3.1.1 and 3.2.1 and also Figures 1 and 3), MP2 single point energies are additionally computed with T/Q-zeta and Q/5-zeta basis set extrapolations. For both CO and H2O, additional MP2 and DFT+D2 single point calculations are performed on cluster models as large as Mg25O25. To carry out hybrid MP2:DFT +D2 structure optimizations, DFT+D2 cluster and periodic calculations use the VASP code and plane waves. For CO on MgO, DFT+D2 single point calculations going from Mg9O9 to Mg25O25 clusters are computed with the TURBOMOLE code and the “def2”-TZVP basis set46 of Weigend and Ahlrichs. This is motivated by the superior computational efficiency of the TURBOMOLE code for larger clusters and the fact that the difference between DFT results with VASP/plane waves and TURBOMOLE/Gaussian basis sets is less than 0.05 kJ/mol (Table S1 in the Supporting Information). On the contrary, for H2O, this difference is 2.1 ± 1.1 kJ/mol, which reduces just to 1.1 ± 0.6 kJ/mol using the larger “def2”-QZVP basis set46 (Table S1 in the Supporting Information). Therefore, for this system, DFT+D2 single point calculations on larger cluster models are rather performed using VASP/plane waves. This requires a further increase of the cell size to 25 Å in order to avoid artificial interactions with the cluster replicas. A comparison between plane wave and Gaussian-type orbital basis sets as a function of the basis set quality can be found in Table S2 of the Supporting Information. As the last step, CCSD(T) single point calculations are performed for both CO and H2O on MgO(001) with the same computational setting as adopted for the MP2 calculations. These calculations were feasible only for the smaller Mg9O9 cluster model cut out from the hybrid MP2:DFT+D2 equilibrium structures of CO and H2O·MgO(001) at Θ = 1/ 8 coverage and with D/T-zeta basis sets. All the DFT+D2, MP2, and CCSD(T) cluster calculations are BSSE corrected.

Figure 1. Top view of MgO(001) slab (a) and side view of the CO(1/8)·MgO(001) system (b). The CO·Mg9O9 embedded cluster is highlighted as solid spheres. Color code: magnesium, green; oxygen, red; carbon, gray.

limit, and the BSSE was corrected only a posteriori. The present work differs in two respects: the hybrid MP2:DFT+D2 optimization uses basis set extrapolation with BSSE corrections, and the MP2 high-level method is combined with both PBE+D2 and B3LYP+D2 as low-level methods, whereas previous calculations have been performed with PBE+D2 only. The DFT+D2 adsorption energies and selected bond distances are given in Table 1. The adsorption energies, ΔE, and their dispersion contributions, ΔED, are compared for PBE +D2 and B3LYP+D2 functionals. The hybrid B3LYP+D2 functional yields a smaller total adsorption energy than PBE+D2, −16.9 compared to −22.1 kJ/mol. The dispersion contribution is larger for B3LYP+D2 (−9.6 kJ/mol) than for PBE+D2 (−7.0 kJ/mol), while the pure B3LYP adsorption energy is only 50% of the PBE one. The Mg2+···CO adsorption complex distance with B3LYP+D2 is 8.8 pm longer than with PBE+D2. The C−O bond distance after adsorption gets 0.1 pm longer at the PBE+D2 level but 0.1 pm shorter with the B3LYP+D2 functional. Due to the relative weak interaction between CO and the MgO(001) surface, the relaxation energies of the surface and the CO molecule are of the order of 1 kJ/mol only (about 5% of the total interaction energy). The DFT+D2 optimized structures are refined by hybrid MP2:DFT+D2 calculations. Table 2 shows the results of three different variants: (i) MP2/CBS(D,T):PBE+D2, (ii) MP2/ CPC-CBS(D,T):PBE+D2, and (iii) MP2/CPC-CBS(D,T):B3LYP+D2. The latter two use the CPC-CBS scheme as high-level method and, thus, describe BSSE- and BSIE-free PESs. The hybrid MP2/CBS(D,T):PBE+D2 scheme yields a 7.6 pm longer Mg2+···CO bond distance than plain PBE+D2,

3. HYBRID MP2:DFT-D RESULTS 3.1. CO on MgO(001). 3.1.1. Surface and Cluster Models. The MgO(001) slab model uses a supercell of dimension 8.48 Å × 8.48 Å × 25 Å, following previous work4 (Figure 1). The two bottom layers of the four-layer slab are kept frozen during the optimization as well as the cell volume. The adopted periodic adsorption model has one CO molecule per supercell perpendicularly oriented to the surface with the C atom pointing to one of the eight pentacoordinated terrace site, [Mg2+]5c, which corresponds to a low coverage Θ = 1/8 (therefore referred to as “CO(1/8)· MgO(001)”). For this low coverage, the repulsive lateral interactions are negligible, as shown in ref 4. In all the hybrid MP2:DFT+D2 structure optimizations, a neutral two layer Mg9O9 model25 of the MgO(001) surface is adopted with one CO ad-molecule (Figure 1). We refrain from embedding in point charges because overpolarization may occur in particular with extended basis sets.4 A posteriori, to verify the effects of using a larger cluster size, the latter is increased from Mg9O9 to Mg25O25. 3.1.2. Hybrid MP2:DFT-D Results. A periodic MP2 estimate of the adsorption energy of CO on MgO(001) was obtained before by hybrid MP2:PBE+D2 calculations.4 In that study, the MP2 energies and forces were extrapolated to the CBS D


Article


Table 1. Mg2+···CO and (C-O)ads Bond Distances in pm, Adsorption Energies, ΔE, Dispersion Contributions, ΔED, and Relaxation Energies, ΔErelax, in kJ/mol for PBE+D2 and B3LYP+D2 Functionals PBE+D2 B3LYP+D2

R(Mg2+···CO)

R(C-O)ads

R(C-O)gasa

ΔE

ΔED

ΔErelaxb

240.8 249.6

114.4 113.3

114.3 113.4

−22.1 −16.9

−7.0 −9.6

−1.1 −1.0

a

Isolated CO molecule. bSee eqs 8 and 9 for definition.

Table 2. Mg2+···CO and (C-O)ads Bond Distances in pm, Hybrid MP2:DFT+D2 Adsorption Energies, ΔEHL:LL,CPC, and Their Contributions in kJ/mola Hybrid scheme

MP2/CBS(D,T):PBE+D2

MP2/CPC-CBS(D,T):PBE+D2

MP2/CPC-CBS(D,T):B3LYP+D2

R(Mg2+···C) R(C−O)adsd ΔELL(pbc) ΔELL(C) ΔEHL,CPC(C) BSSEe ΔLR(pbc,C) ΔHL(C) ΔErelax,HL:LL(pbc) ΔEHL:LL,CPC(pbc)

248.4 [240.8]b 113.2 [114.4]b −21.6 (−22.7) −21.2 (−19.0) −15.8 (−13.8) 5.5 −0.4 (−3.6) 5.4 (5.2) (−1.3) −16.2 (−17.4)

252.2 113.2 −21.3 (−22.4) −21.1 (−18.7) −16.1 (−14.0) 5.3 −0.2 (−3.7) 5.0 (4.7) (−1.4) −16.3 (−17.7)

249.4 [249.6]c 113.3 [113.3]c −14.7 (−16.1) −13.1 (−10.3) −16.8 (−14.1) 5.5 −1.6 (−5.7) −3.8 (−3.8) (−1.3) −18.5 (−19.8)

See eqs 10 and 11 for definition. In parentheses: ΔẼ values computed at the corresponding hybrid equilibrium structures of CO on MgO(001), see eq 7. bPBE+D2 values in square brackets (Table 1). cB3LYP+D2 values in square brackets (Table 1). dR(C-O)gas = 113.5 pm with MP2/ CBS(D,T). eBSSE calculated at the respective hybrid QM:QM structure from MP2/CBS(D,T) “monomer” and “ghost” calculations. a

Table 3. Contributions to Hybrid MP2/CPC-CBS(D,T):DFT+D2 Adsorption Energy (kJ/mol) for Two Cluster Sizes: CO· Mg9O9 and CO·Mg25O25a PBE+D2 ΔEHL,CPC(C) ΔELL,CPC(C)d ΔELL(pbc) ΔLR(pbc,C)d ΔHL(C)d ΔẼ HL:LL,CPC(pbc)f

B3LYP+D2

CO·Mg9O9

CO·Mg25O25

Δb

CO·Mg9O9

CO·Mg25O25

Δb

−14.0 (3.7)c −18.7 (−5.2)e

−18.9 (−0.3)c −21.5 (−6.1)e −22.4 (−6.7)e −0.9 (−0.6)e 2.6 −19.8

−4.9 (−4.0)c −2.8 (−0.9)e

−14.1 (4.3)c −10.3 (−7.4)e

−5.1 (−4.2)c −4.0 (−1.3)e

+2.8 (+0.9)e −2.1 −2.1

−5.8 (−2.2)e −3.8 −19.8

−19.1 (0.1)c −14.3 (−8.6)e −16.1 (−9.5)e −1.8 (−0.9)e −4.8 −20.9

−3.7 (−1.5)e 4.8 −17.7

+4.0 (+1.3)e −1.1 −1.1

a

All the values are BSSE corrected at the hybrid MP2/CPC-CBS(D,T):DFT+D2 equilibrium structure of CO(1/8)·MgO(001). bDifference between the two cluster models. cThe Hartree−Fock contributions to the MP2 energies are reported in parentheses. dThe low-level cluster calculations use the TURBOMOLE code and the def2-TZVP basis set. Therefore, the DFT+D2 cluster energy contributions and related terms may differ from those in Table 2 (the differences are of the order of 0.03 kJ/mol). eDispersion contributions to the DFT+D2 energies are reported in parentheses. fSee eq 7 for definition.

The hybrid MP2/CPC-CBS(D,T):B3LYP+D2 results show smaller deviations from the B3LYP+D2 results than observed for the PBE+D2 functional. The adsorption energy gets 1.6 kJ/ mol more negative, and the Mg2+···CO distance gets 0.2 pm shorter. Both B3LYP+D2 and hybrid MP2/CPC-CBS(D,T):B3LYP+D2 describe a C−O bond shortening, by 0.1 and 0.2 pm respectively. In the following, we consider energy changes beyond the hybrid results of Table 2. At the hybrid MP2/CPCCBS(D,T):PBE+D2 equilibrium structure of CO(1/8)· MgO(001), the effect of using larger basis sets for CBS limit extrapolation is evaluated on the smaller CO·Mg9O9 cluster as the difference between MP2/CPC-CBS(Q,5) and MP2/CPCCBS(D,T) adsorption energies

whereas correcting for BSSE increases the distance by another 3.8 pm, a non-negligible increment. Both hybrid calculations predict a 0.3 pm shortening of the C−O bond distance, while PBE+D2 yields a lengthening of 0.1 pm. The MP2/ CBS(D,T):PBE+D2 and MP2/CPC−CBS(D,T):PBE+D2 adsorption energies are 5.9 and 5.8 kJ/mol, respectively, less negative compared to the PBE+D2 value. Despite the basis set extrapolation, BSSE of the MP2 energies is substantial, between 5.3 and 5.5 kJ/mol. It assumes different values for each hybrid coupling variant because the equilibrium structures are different. The increase in the Mg2+··· CO distance from 248.4 pm for MP2/CBS(D,T):PBE+D2 to 252.2 pm for MP2/CPC-CBS(D,T):PBE+D2 leads to a decrease in BSSE from 5.5 to 5.3 kJ/mol, respectively, in correspondence with a smaller orbital overlap between the fragments.

CBS(Q,5) CBS(D,T) Δbasis = ΔEMP2,CPC − ΔEMP2,CPC

E

(13)


Article

Journal of Chemical Theory and Computation The MP2 results together with the Hartree−Fock and correlation contributions for each basis set are fully reported in Table S3 of Supporting Information. Additionally, Table S4 in Supporting Information shows the BSSE-corrected and uncorrected MP2 results for three different extrapolation schemes: CBS(D,T), CBS(T,Q), and CBS(Q,5). Already the CBS(D,T) scheme yields converged results for the BSSEcorrected MP2 energies with a negligible basis set extension effect, Δbasis, of just 0.05 kJ/mol. Correcting for the BSSE appears important. Extrapolating the BSSE-uncorrected MP2 results yields a CBS(Q,5) value that is 1.0 kJ/mol more binding. To examine the effect of increasing the cluster size, single point calculations are performed on a larger cluster model, CO· Mg25O25, which is cut out from the hybrid-adsorbed molecule−surface equilibrium structures. Table 3 shows the hybrid MP2/CPC-CBS(D,T):DFT+D2 adsorption energies for the CO·Mg9O9 and CO·Mg25O25 clusters and both the PBE+D2 and B3LYP+D2 functionals. For the PBE+D2 functional, the long-range correction, ΔLR, decreases in absolute value from −3.7 to −0.9 kJ/mol as the cluster size grows, since a larger share of dispersive forces is accounted for by the low-level cluster contribution (see the ΔELL(C) dispersion contribution in parentheses of Table 3). The high-level contribution, ΔEHL(C), increases with the cluster size balancing the opposite long-range correction decrease. The influence of the cluster size on the hybrid adsorption energy is about 2 kJ/mol, showing as MP2/CPCCBS(D,T) compensates the low-level description. In an analogous way and at the MP2/CPC-CBS(D,T):B3LYP+D2 equilibrium structure, the long-range correction is evaluated for the two considered different cluster sizes with the B3LYP+D2 functional. Similar to PBE+D2 results, a decrease in the longrange correction, ΔLR, with the cluster size is observed, while a better balance with the increasing of the high-level contribution shows a smaller influence on the cluster size, just 1 kJ/mol. As suggested previously, e.g., see refs 4 and 47, using two different cluster sizes allows the calculation of scaling factors for ΔLR f=

̃ (CO ·Mg O25) ΔE(MP2,pbc) = ΔEMP2 25 + f ·ΔLR (CO·Mg 25O25) + Δbasis − ΔErelax,HL:LL

where ΔEMP2 is the MP2/CPC-CBS(D,T) adsorption energy, f · ΔLR is the scaled long-range correction of CO·Mg25O25, Δbasis is the basis set extension correction, and ΔErelax,HL:LL is the hybrid relaxation energy, both computed for CO·Mg9O9. As shown in Tables S3 and S4 of the Supporting Information, the correction for the basis set extension does not play a significant role in the present study (less than 0.1 kJ/mol). The estimated periodic MP2 results including the PBE+D2 and B3LYP+D2 scaled long-range corrections, i.e., eq 15, are −19.2 and −20.0 kJ/mol, respectively (Table 4). Figure 2 compares DFT+D2, hybrid MP2:DFT+D2, and periodic MP2 estimates obtained using PBE+D2 and B3LYP

Figure 2. Adsorption energies for DFT+D2 and hybrid MP2:DFT +D2 method including the CPC-CBS scheme and periodic MP2 estimates. Color code: red, PBE+D2; blue, B3LYP+D2 functional.

+D2 as low-level methods. The hybrid MP2/CPC-CBS:DFT +D2 adsorption energies (Table 2) differ substantially from the DFT+D2 results (Table 1). The effect for PBE and B3LYP is in the opposite direction. The hybrid with MP2 makes the PBE +D2 result 5.8 kJ/mol less binding, whereas the B3LYP+D2 results becomes 1.6 kJ/mol more binding, so that the relative order is exchanged. With the scaled long-range corrections, the difference between the PBE+D2- and B3LYP+D2-based predictions gets smaller. The final uncertainty due by the choice of the low-level method decreases to 1 kJ/mol for the periodic MP2 estimates, well within chemical accuracy limits, and the predicted Mg2+···CO distances differ by 3 pm only. 3.2. H2O on MgO(001). 3.2.1. Surface and Cluster Models. A hybrid MP2:DFT+D2 scheme is applied for predicting the adsorption energy and structure of H2O on the MgO(001) surface at the limiting zero coverage, i.e., one isolated molecule in a large enough supercell. Two different low coverages corresponding to two different surface models (Figure 3) and computational settings are considered. Table 5 collects the main differences between the two systems. The first one has a periodic adsorption surface model with one of eight surface Mg2+ ions occupied, which corresponds to a coverage of Θ = 1/8 and is named “H2O(1/8)·MgO(001)”. Its computational setting resembles the one previously used for benchmarking the adsorption energies of CO,4 CH4,3,10 and C2H610 on the MgO(001) surface and also applied here for CO on MgO(001). The second periodic adsorption model is one of 16 Mg 2+ ions occupied with H2O molecules

ΔEMP2[(CO·Mg 25O25) − (CO·Mg 9O9)] ΔE DFT + D2[(CO·Mg 25O25) − (CO·Mg 9O9)]

(14)

which account for the difference in the effect of enlarging the cluster size on MP2 and DFT+D2 energies. With f = 1.8 for PBE+D2 and f = 1.3 for B3LYP+D2 (Table 4), both functionals underestimate the effects of increasing the cluster size compared to MP2 and PBE+D2 to a larger extent. The periodic MP2 adsorption energies are then estimated as Table 4. Estimated Periodic MP2 Adsorption Energies, ΔEMP2(pbc)est., in kJ/mol PBE+D2 B3LYP+D2

f

ΔLR

f · ΔLR

ΔErelax

ΔEMP2(pbc)est

1.8 1.3

−0.9 −1.8

−1.7 −2.3

−1.4 −1.3a

−19.2b −20.0b

a

(15)

a ΔErelax values computed for the CO·Mg9O9 cluster are used, see Table 2. bFor ΔE(MP2,CO·Mg25O25), see Table 3.

F


Article


site. One of the H atoms points toward one surface O2− ion, while the other H atom points away from the surface as also previously observed (Figure 4).55,56

Figure 4. Cut out from the H2O·MgO(001)periodic adsorption model. R1 and R2 are the Mg2+···O and O2−···H(OH) distances. Color code: magnesium, green; oxygen, red; hydrogen, white.

Figure 3. Top view of H2O(1/8)·MgO(001) (a) and H2O(1/16)· MgO(001) (b) systems. H2O·Mg9O9 embedded clusters are highlighted as solid spheres. Color code: magnesium, green; oxygen, red; hydrogen, white.

For H2O(1/8)·MgO(001), the MgO(001) slab model is built following the one of CO on MgO(001) (Figure 3). For the H2O(1/16)·MgO(001) system, a larger 11.83 Å × 11.83 Å × 30 Å supercell is used (Figure 3).48 The surface models for both systems are created with four MgO layers of which the two topmost are fully relaxed. The surface cell is kept frozen during structure optimizations. For both H2O(1/8)·MgO(001) and H2O(1/16)·MgO(001), hybrid optimizations use a cluster model consisting of a H2O molecule on a Mg9O9 surface model25 (Figure 3), and single point calculations are performed also on the larger H2O/ Mg25O25 cluster model. 3.2.2. Hybrid MP2:DFT-D Results. The H2O(1/8)·MgO(001) and H2O(1/16)·MgO(001) systems (see Table 5 for a comparison) are optimized with PBE+D2 and PBE+D2[Mg], respectively. Harmonic frequencies are calculated on the corresponding energy minimum structures. For the adsorption complexes, distances R1 and R2 (Figure 4), adsorption energies, ΔE, and their dispersion contributions, ΔED, are given in Table 5. R1 is the distance between the Mg2+ ion and the O atom of the H2O molecule, and R2 is the distance between the O2− ion and the H atom of the H2O molecule pointing to the surface. For H2O(1/16)·MgO(001), the distances R1 and R2 are 3.5 and 4.8 pm, respectively, larger than for H2O(1/8)·MgO(001), whereas the adsorption energy is 6 kJ/mol more binding. This

corresponding to a coverage of Θ = 1/16, in the following labeled “H2O(1/16)·MgO(001)”. It uses the same DFT computational details and surface models adopted previously by Ončaḱ et al.48 Comparing the results will tell us how stable our hybrid MP2:DFT+D2 results are with respect to reasonable changes of the surface model and the computational details. At monolayer coverage and in a temperature range of 100− 200 K, experimental21,22 and computational49−53 studies agree that H2O molecules on the MgO(001) surfaces give rise to a mixed layer of both hydroxyl groups and water molecules. At even higher coverages dissolution of the MgO surface can occur,48,54 in other words the Mg2+ ions become more inclined to leave the surface. Dissociation of H2O on the MgO surface is driven by the presence of either surface defects or neighboring adsorbed molecules, both stabilizing the dissociated products. Therefore, at the low coverages we consider, H2O adsorbs only molecularly on the regular MgO(001) surface.55 The structure optimizations show that the molecular plane is not parallel to the surface but slightly tilted (≈ 50° with respect to the normal of the MgO(001) surface). The molecule adsorbs with its oxygen atom nearly above the Mg2+

Table 5. PBE+D2 and PBE+D2[Mg] Results of H2O(1/8)·MgO(001) and H2O(1/16)·MgO(001) Systems, Respectivelya System Coverage, Θ Slab PBE+D2 (VASP): Basis set Brillouin zone sampling Dispersion (D2) Notation R1(Mg2+···O) R2(O2−···H(OH)) ΔE ΔED ΔErelaxe

H2O(1/8)·MgO

H2O(1/16)·MgO

1/8 8.48 Å × 8.48 Å × 25 Å four-layer thick

1/16 11.83 Å × 11.83 Å × 30 Å four-layer thick

PWs/600 eV and PAW (Mg_pv) 2 × 2 × 1 k-mesh and Γ point “Ne” ref 3 Ewald summation PBE+D2

PWs/400 eV and PAW Γ point “Mg” ref 39 30 Å cutoff PBE+D2[Mg]

223.4 174.9 −55.0 −10.8 (−18.8)b,c −12.5 (−9.0)f

226.9 179.7 −61.0 −20.0 (−11.3)b,d −10.1 (−7.5)f

a Mg2+···O and O2−···H(OH) distances in pm, adsorption energies, ΔE, dispersion energy contributions, ΔED, and relaxation energies, ΔErelax, in kJ/mol. bSingle point calculations in parentheses (see also Table S7 in the Supporting Information). cMg dispersion parameters for Mg2+ in parentheses. dNe dispersion parameters for Mg2+ in parentheses. eSee eqs 8 and 9 for definition. fMgO(001) surface relaxation contribution to the total relaxation energy in parentheses.

G


Article


Table 6. MP2/CPC-CBS(D,T):PBE+D2 and MP2/CPC-CBS(D,T):PBE+D2[Mg] Results of H2O(1/8)·MgO(001) and H2O(1/16)·MgO(001) System, Respectivelya H2O(1/8)·MgO Structure

b

R1(Mg2+···O) R2(O2−···H(OH)) ΔELL(pbc) ΔELL(C) ΔEHL,CPC(C) BSSEc ΔLR(pbc,C) ΔHL(C) ΔErelax,HL:LL(pbc) ΔEHL:LL,CPC(pbc)

PBE+D2

H2O(1/16)·MgO


223.4 174.9 −55.0 −51.3 −43.2 7.2 −3.7 8.1

220.7 179.7 −54.2 −50.5 −44.3 6.9 −3.7 6.2

−46.9

−48.0

(−61.3) (−57.4) (−55.6) (−3.9) (1.8) (−11.5) (−59.5)

PBE+D2[Mg]

MP2/CPC-CBS(D,T):PBE+D2[Mg]

226.9 179.7 −61.0 −58.4 −43.0 7.2 −2.6 15.4

221.3 180.5 −63.3 −56.0 −44.9 6.8 −7.3 11.1

−45.6

−52.2

(−68.2) (−62.1) (−53.1) (−6.2) (8.9) (−7.1) (−59.3)

Mg2+···O and O2−···H(OH) distances in pm, hybrid adsorption energies, ΔEHL:LL, and their contributions in kJ/mol. See eqs 10 and 11 for definition. In parentheses: ΔẼ values computed at the corresponding hybrid equilibrium structures of H2O on MgO(001), see eq 7. bReference structure at which the hybrid calculations are performed. cBSSE calculated at the respective hybrid MP2/CPC-CBS(D,T):DFT+D2 equilibrium structures from MP2/CBS(D,T) “monomer” and “ghost” calculations. a

Table 7. Contributions to Hybrid MP2/CPC-CBS(D,T):DFT+D2 adsorption energy of H2O(1/8)·MgO(001) and H2O(1/ 16)·MgO(001) (kJ/mol) for Two Cluster Sizes: H2O·Mg9O9 and H2O·Mg25O25a PBE+D2b ΔEHL,CPC(C) ΔELL(C) ΔELL(pbc) ΔLR(pbc,C) ΔHL(C) ΔẼ HL:LL,CPC(pbc)g

PBE+D2[Mg]c

H2O·Mg9O9

H2O·Mg25O25

Δ

−55.6 −57.4 (−8.1)f −61.3 (−10.0)f −3.9 (−1.9)f 1.8 −59.5

−60.1 −60.5 (−9.4)f

−4.5 −3.1 (−1.3)f

−0.8 (−0.6)f 0.4 −61.3

+3.1 (+1.3)f −1.4 −1.4

d

H2O·Mg9O9

H2O·Mg25O25

Δd

−53.1 −62.3e (−14.1)f −68.2 (−18.2)f −5.9e (−4.0)f 9.2e −59.1

−57.7 −66.7 (−16.9)f

−4.6 −4.4 (−2.7)f

−1.5 (−1.3)f 9.0 −59.2

+ 4.4 (+2.7)f −0.2 −0.2

a

All the values are BSSE-corrected. bAt the hybrid MP2/CPC-CBS(D,T):PBE+D2 equilibrium structure of H2O(1/8)·MgO(001). cAt the hybrid MP2/CPC-CBS(D,T):PBE+D2[Mg] equilibrium structure of H2O(1/16)·MgO(001). dDifference between the two cluster models. eThe low-level cluster calculations use the VASP/plane waves code and the lattice parameter is further increased to 25 Å. Therefore, the DFT+D2 cluster energy contributions and related terms may differ from those in Table 6 (differences are negligible for PBE+D2 and of the order of 0.2 kJ/mol for PBE +D2[Mg]). fDispersion contributions to the DFT+D2 energies are reported between parentheses. gSee eq 7 for definition.

the same surface (of the order of 1 pm only),10 where the interaction was much weaker and purely dispersive. DFT+D2 equilibrium structures are refined by hybrid MP2/ CPC-CBS:DFT+D2 calculations, which use D/T-zeta basis set extrapolation for the H2O·Mg9O9 cluster. Table 6 shows the results. With MP2:PBE+D2, the R1 distances are shorter than with PBE+D2, 2.7 pm for H2O(1/8)·MgO(001) and 5.6 pm for H2O(1/16)·MgO(001), while the R2 distances are 4.8 and 0.8 pm, respectively, longer. The hybrid MP2:PBE+D2 adsorption energies are 7.1 and 8.8 kJ/mol, respectively, less binding than the corresponding PBE+D2 values. As before for CO on MgO(001), we evaluate also for H2O the effects of using larger basis sets and cluster sizes on the hybrid MP2:DFT+D2 results. At the hybrid MP2/CPCCBS(D,T):PBE+D2 equilibrium structure of H2O(1/8)·MgO(001), extending the basis sets from D/T-zeta to Q/5-zeta yields a non-negligible basis set extension correction, Δbasis, of −1.2 kJ/mol on the MP2/CPC-CBS adsorption energy of H 2 O·Mg 9 O 9 (Tables S5 and S6 in the Supporting Information). Only at the MP2/CPC-CBS(T,Q) level can the adsorption energy be considered converged with Δbasis of −0.04 kJ/mol. The BSSE correction is important also when the largest Q/5-zeta basis sets are used. Indeed, the uncorrected

is ascribed to the different D2 parameters adopted, which leads to a dispersion contribution that is 9.2 kJ/mol larger for H2O(1/16)·MgO(001) than for H2O(1/8)·MgO(001). Indeed, PBE+D2 single point calculations at the PBE+D2[Mg] equilibrium structure of H2O(1/16)·MgO(001) reproduce within 1 kJ/mol the dispersion energy computed with PBE+D2 for H2O(1/8)·MgO(001) and vice versa (see also Table S7 in Supporting Information). Due to the relative strong H2O−MgO(001) surface interaction, the relaxation energies of the surface and of the H2O molecule cover almost 20% of the total adsorption energy, and over 70% of it comes from the surface relaxation. Surface reconstruction upon H2O adsorption is particularly relevant at higher coverages, when in the presence of dissociated H2O molecules surface Mg2+ ions may leave the topmost surface layer, see, e.g. refs 48 and 54. Despite these high coverage cases, the degree of rumpling of the 001 surface in the presence of H2O is not negligible even at our low coverages. The Mg2+ and O2− ions involved in the H2O− MgO(001) surface interaction relax 5−7 pm upward with respect to the clean MgO(001) surface (Figure S1 in the Supporting Information). Here, the surface flexibility is much more pronounced than for the alkane monolayer adsorption on H


Article


H2O(1/16)·MgO(001) as low-level methods. As observed already for the CO on MgO(001), applying both the high-level correction and the scaling of the long-range correction reduces the uncertainty coming from the low-level method as the difference between PBE+D2- and PBE+D2[Mg]-based results decrease from 6 (DFT+D2) to 2.5 kJ/mol (periodic MP2 estimates). Correspondingly, performing hybrid structure optimizations yields R1 and R2 distances with MP2:PBE+D2 and MP2:PBE+D2[Mg] schemes 3.5 pm closer than those observed with PBE+D2 and PBE+D2[Mg] low-level methods (Table 5). These are additional confirmations which show as well that hybrid energies and structures are converged.

MP2/CBS(Q,5) adsorption energy is 1.4 kJ/mol more binding than that of the BSSE-corrected value. Table 7 shows the results of using a larger Mg25O25 cluster model for both H2O(1/8)·MgO(001) and H2O(1/16)·MgO(001) systems at their corresponding hybrid MP2:DFT+D2 equilibrium structures. For H2O(1/8)·MgO(001), going from H2O·Mg9O9 to H2O·Mg25O25, the D2 dispersion term increases from −8.1 to −9.4 kJ/mol, correspondingly, and the PBE+D2 long-range correction, ΔLR, decreases in absolute value from −3.9 to −0.8 kJ/mol. On the contrary, the MP2 contribution increases in absolute value by 4.5 kJ/mol as the cluster size grows. The total effect of extending the cluster size on the hybrid energy is 1.4 kJ/mol only, which shows a good balance between the opposite behaviors of ΔLR and MP2 terms. A similar trend is observed for H2O(1/16)·MgO(001). Using larger cluster sizes, the decrease in the PBE+D2[Mg] long-range correction (+4.4 kJ/mol) is well counterbalanced by the opposite increase of the MP2 term (−4.6 kJ/mol), and the total effect on the final hybrid energy is negligible. For the PBE+D2 functional, the scaling factor, f, of the longrange correction (eq 14) is 1.5. In other words, PBE+D2 underestimates the effect of extending the cluster size with respect to MP2. The MP2 limit for the H2O(1/8)·MgO(001) periodic structure is estimated following eq 15

4. COMPARISON WITH PERIODIC LOCAL MP2 To investigate how close the results of the hybrid MP2:DFT +D2 scheme approach the thermodynamic MP2 limit for the CO(1/8)·MgO(001) and H2O(1/8)·MgO(001) systems, we performed single point periodic HF and periodic local MP2 (LMP2) calculations at the hybrid MP2:DFT+D2 equilibrium structures using the CRYSTAL57 and CRYSCOR11 codes. The structure of the studied systems is provided in the Supporting Information in a form of CRYSTAL input. Basis sets of augmented-triple-ζ quality are employed (in the following referred to as AVTZ). Because of the convergence problems in periodic HF calculations for rich and diffuse basis sets, to construct these basis sets the protocol of ref 58 is followed. The diffuse orbitals from the augmented part of the basis are processed only at the LMP2 stage of the calculations via the dual basis set scheme. The basis sets and technical parameters of the HF and LMP2 calculations are given in Section S6 of the Supporting Information. In order to estimate the basis set limit for the correlation energy, the periodic local F12 correction is calculated. The correction corresponds to the 3*A(FIX) variant of the LMP2F12 method.59 The LMP2 and LMP2-F12 calculations employing AVTZ basis set are carried out within the frozen core (FC) approximation; i.e., only the valence localized orbitals centered on oxygen or on O−H bonds are correlated. Next, to estimate the Mg-core contribution to binding, we calculate the corresponding correction as a difference between the core-correlated (2s, 2p-Mg part of the core) and the frozen-core LMP2 energies. In these calculations, the tight functions from the cc-pwCVTZ basis set45 are added for the Mg atoms (we note this basis as ACVTZ). For the pair-specific virtual space in the LMP2 calculations, we employ extended PAO domains. For the MgO slab, that corresponds to inclusion up to the second-nearest-neighbor atoms with respect to the oxygens, on which the corresponding localized orbitals are centered. The same extended domains are used as the RI domains in the F12 calculations. Only for the very localized Mg-core orbitals in the core-correlated calculations the domains are not extended. In the local correlation treatment of molecule−surface interactions, the localized occupied orbitals are usually well localized either on the slab or on the ad-molecule (adM). This allows one to identify the intra-slab, intra-adM, and inter-adMslab orbital pairs and thus the corresponding correlation energy contributions to the interaction energy. The intra-slab and intra-adM interaction energies represent the correlated exchange and electrostatics contributions (often repulsive). The pair energies in such contributions decay exponentially with the interorbital distance within the pair. The inter-adMslab pair energies, associated with van der Waals dispersion,

ΔE(MP2,pbc) = −60.1 − (1.5 × 0.8) − 1.2 + 11.5 = −51.0kJ/mol

(16)

The scaling factor, f, of the PBE+D2[Mg] long-range correction is instead 1.0. Adjusting the PBE+D2[Mg] longrange correction for this factor allows the estimation of the periodic MP2 adsorption energy also for H2O(1/16)·MgO(001), which yields ΔE(MP2,pbc) = −57.7 − (1.0 × 1.5) − 1.2 + 7.1 = −53.4.kJ/mol

(17)

Both the results of eqs 16 and 17 are 3.0 and 1.2 kJ/mol, respectively, more binding than the corresponding hybrid MP2:DFT+D2 values (Table 6), and they agree within 2.5 kJ/ mol. Figure 5 shows the adsorption energies for DFT+D2, MP2:DFT+D2, and periodic MP2 estimates obtained using PBE+D2 for H2O(1/8)·MgO(001) and PBE+D2[Mg] for

Figure 5. Adsorption energies for DFT+D2 and hybrid MP2:DFT +D2 method including the CPC-CBS scheme and periodic MP2 estimates. Color code: red, PBE+D2 for H2O(1/8)·MgO(001); blue, PBE+D2[Mg] for H2O(1/16)·MgO(001). I


Article


Table 8. Periodic HF and LMP2 Adsorption Energies of CO(1/8)·MgO(001) and H2O(1/8)·MgO(001) Systems in kJ/mola LMP2 HF

Intra-MgO

Intra-adM

Inter-adM-MgO

Inter-slab replicas

Total

Total HF+LMP2

−16.6 −2.3 −19.0 −2.2 −21.2

−17.1

CO(1/8)·MgO −0.5

2-layer,FC 2-layer,core-corr. ACVTZ-limit 2-layer,FC,F12 BS-limit

+5.2 +0.7 +5.9 −0.2 +5.7

−0.5 −0.5

2-layer,FC 2-layer,core-corr. ACVTZ-limit 2-layer,FC,F12 BS-limit

−37.0 −37.0 −37.0

−24.6 −3.0 −27.6 −1.8 −29.4 H2O(1/8)·MgO

+3.4 +0.0 +3.4 −0.2 +3.2 +15.0 +1.7 +16.6 −0.7 +16.0

−0.6 −0.0 −0.6 −0.6

−41.7 −5.0 −46.6 −2.9 −49.5

+10.0 +0.0 +10.0 −1.1 +8.9

−0.8 −0.0 −0.8 −0.8

−17.6 −3.3 −20.9 −4.7 −25.6

−19.4 −21.6 −54.6 −57.9 −62.5

a

LMP2 adsorption energy is partitioned into intra-slab and intra-adM contributions, inter-slab-adM energy, and very-distant-pair contributions from the slab replicas.

Table 9. Periodic LMP2 and Hybrid MP2:DFT+D2 Adsorption Energies of CO(1/8)·MgO(001) and H2O(1/8)·MgO(001) Systems in kJ/mol CO(1/8)·MgOa Periodic LMP2 ΔEMP2/ACVTZ(C) fe ΔLR(pbc,C)d f · ΔLR(pbc,C)d ΔEHF/ACVTZ(pbc) ΔEcorr/ACVTZ(pbc) ΔEMP2/ACVTZ(pbc) BS-correction (HF) BS-correction (corr.) BS-correction (MP2) ΔEcorr/CBS(pbc) ΔEMP2/CBS(pbc)

H2O(1/8)·MgOb Δ

MP2:B3LYP+D2

c

Periodic LMP2

−10.8 1.3 −5.8 −7.3

d

−0.5 −19.0 −19.4 −2.2h −2.2h −21.2 −21.6

MP2:PBE+D2

Δc

−52.0 1.5 −3.9 −5.7

−18.1f −0.6g −2.7g −3.2g

−1.3

−37.0 −20.9 −57.9

−21.3

−0.3

−4.7h −4.7h −25.6 −62.5

a

−57.7f −0.3g −4.5g −4.8g −62.5

−0.2

0.0

b

At the hybrid MP2/CPC-CBS(D,T):B3LYP+D2 equilibrium structure of CO(1/8)·MgO(001). At the hybrid MP2/CPC-CBS(D,T):PBE+D2 equilibrium structure of H2O(1/8)·MgO(001). cDeviation between periodic LMP2 results and periodic MP2 estimates by the hybrid scheme. d The ad-molecule/Mg9O9 cluster model is considered. eScaling factors of the long-range corrections are obtained with MP2/CBS(D,T) as the high-level method. fPeriodic MP2 estimate, ΔEMP2/ACVTZ(pbc)est. = ΔEMP2/ACVTZ(C) + f · ΔLR(pbc,C). gBasis set extension correction using the difference between MP2/CBS(Q,5) and MP2/ACVTZ energies of the ad-molecule/Mg9O9 cluster model. hBasis set extension correction of the correlation contribution as provided by the periodic F12 limit

decay much slower, as the sixth inverse power of the interorbital distance. Hence, a much larger cutoff distance has to be used for such pairs. In our calculations, we use conservative values for the cutoff thresholds: 14 Å for interpairs and 6 Å for the intraslab ones. The F12 correction is intrinsically short range, since the corresponding pair contributions decay exponentially. Therefore, the cutoff thresholds were reduced to 8 and 4 Å for the inter- and intrapairs, respectively. On the basis of the energies of the explicitly included and processed pairs, one can also estimate the contribution from the interpairs beyond the truncation. For that, one can on the fly compute the C6 coefficients, corresponding to the interpairs for the localized orbitals of the given type11 and interpolate the pair energies of each type outside the active region (given by the 14 Å cutoff) by means of the long-range part of the Lennard-Jones potential. Furthermore, having obtained these C6 coefficients, one can formally replicate the slab to form a semi-infinite crystal and use these slab replicas to approximate the dispersion from beyond the actual slab.60

Table 8 compiles the LMP2 adsorption energies and their individual components for CO(1/8)·MgO(001) and H2O(1/ 8)·MgO(001), respectively. The adsorption energy is defined here as the counterpoise corrected energy difference between the adM−surface system and the isolated (but accompanied by the ghost counterparts) ad-molecule and slab. For CO(1/8)·MgO(001), the HF contribution to binding is small, which is not surprising. As the dipole moment of the CO is minute, the electrostatic interaction in CO(1/8)·MgO(001) is weak and nearly fully compensated by the exchange repulsion. The decisive part of the interaction in this system comes then from the correlation and, in particular, dispersion. In H2O(1/8)·MgO(001), the electrostatic component is much more substantial, leading to a pronounced interaction already at the HF level. Moreover, the correlation part seems to be less important here. It is less than a third of the overall adsorption energy. However, the long-range part thereof, which is associated with dispersive interactions, is much larger and individually forms the biggest contribution. J


Article


Table 10. BSSE-Corrected MP2, CCSD and CCSD(T) Adsorption Energies, ΔE, in kJ/mol, for CO(1/8)·MgO(001) and H2O(1/8)·MgO(001)a System

MP2

CCSD

ΔCCSD

CCSD(T)

Δ(T)

ΔCC

CO(1/8)·MgO

DZ TZ CBS(D,T)

−10.1 −14.5 −16.1

−7.9 −12.2 −13.8

2.2 2.3 2.3

−10.9 −15.9 −17.8

−3.0 −3.7 −4.0

−0.8 −1.4 −1.6

H2O(1/8)·MgO

DZ TZ CBS(D,T)

−34.8 −41.5 −44.3

−32.3 −39.0 −41.7

2.5 2.6 2.6

−35.7 −43.7 −47.0

−3.4 −4.7 −5.3

−0.9 −2.2 −2.7

a

At the hybrid MP2/CPC-CBS(D,T):PBE+D2 equilibrium structures.

kJ/mol region or below. The largest effect came from “ITOL 2” parameter that controls the onset of the multipole approximation for the two-electron integrals for the Coulomb term in HF or DFT, which amounted to around −0.3 kJ/mol. The local approximations in periodic LMP2 did not introduce a substantial error either, at least with the chosen parameters. The error due to the pair-list truncation was eliminated by the very-distant energy extrapolation, while the extension of the domains to up to second nearest neighbor atoms (leading to pair domains with up to 24 atoms) provides virtually the converged energy with the domain size. Finally, the HF/LMP2 interaction energy with the 2-layer, 3-layer, or 4-layer slab model mutually agree within 0.2 kJ/mol suggesting that even the 2-layer slab model is adequate to represent H2O- or COMgO(001) surface interaction if the slab replication technique is used in LMP2. Estimation of the basis set incompleteness error is more problematic. For the LMP2 basis set limit, we employ the periodic F12 correction on top of the AVTZ LMP2 energy. The F12 method substantially accelerates the basis set convergence, but the remaining uncertainty can still remain in the 1 kJ/mol region. From the comparison to the hybrid scheme basis set extrapolations, the error of our F12 results seems to be of the order of 0.5 kJ/mol. In the periodic HF, to approach the basis set limit result is even more difficult. Very rich basis sets are inapplicable in periodic HF due to the convergence problems. The basis set extrapolation for the HF energy, conducted in the hybrid scheme with finite clusters, reveals that the HF/AVTZ result is close to the limit, but the basis set correction is not entirely negligible (e.g., −0.6 kJ/mol in CO·MgO; Table 9). We estimate the overall error bar of the periodic HF/LMP2 calculations to be ±1 kJ/mol.

The energy of the pairs beyond the interorbital cutoff distance of 14 Å is tiny; the overall contribution is not more than a few hundredths of a kJ/mol in both systems. At the same time, the dispersive contribution to the adsorption energy originating from beyond the two-layer slab is not absolutely negligible: −0.6 kJ/mol in CO(1/8)·MgO(001) or −0.8 kJ/mol in H2O(1/8)·MgO(001). Only beyond four layers these interactions become unimportant (about −0.1 kJ/ mol). Now, we are at the position to compare the periodic single point energies to those of the hybrid MP2:DFT+D2 scheme (see Tables S8 and S9 of the Supporting Information for their derivation). We perform MP2 single point calculations on the Mg9O9 cluster model, corresponding to the same structure as the periodic calculations and the same ACVTZ basis set (yet without the dual basis set scheme). In the hybrid method, the periodic estimate is computed as a sum of the MP2/ACVTZ adsorption energy of the cluster and the DFT+D2 long-range correction scaled by the adjustment factor f as previously obtained. For the basis set limit estimate in the hybrid scheme, we use a Q/5-zeta basis set extrapolation scheme for the cluster interaction energy (exponential for the HF energy and inverse third power for the correlation energy). Table 9 lists the results from both methodologies. Comparing the MP2 (both ACVTZ and CBS) results, we note that the agreement between two methodologies is remarkably good (within 1 kJ/mol), especially if taking into account the complete different routes to the final interaction energies. Furthermore, a comparison between the computer time of hybrid MP2:DFT+D2 and LPM2 methods can be found in Table S10 of the Supporting Information. The computational cost of a periodic HF/LMP2 calculation is comparable with that of a hybrid MP2:DFT+D2 scheme, whereas a significant difference of 3 orders of magnitude was shown with respect to a full periodic MP2 within the VASP code (see ref 15 and Supporting Information therein). To estimate the uncertainty in the periodic LMP2 results, we analyze the possible sources of errors. These can come from the (i) technical parameters of the calculations (truncation thresholds, finite k-mesh, auxiliary basis for the density fitting, etc.), (ii) the approximations that the methodology is based upon (the domain- and pair-list approximation in LMP2 or finite basis set in both HF and MP2), and (iii) inadequacy of the structure model (a few-layer slab vs a crystalline surface). Most of these errors can be quantified by studying the influence of tightening of the corresponding thresholds or releasing the approximations. The technical tolerances used in our periodic HF/LMP2 calculations have already been tightened with respect to their defaults, and their further modification affected the interaction energies only in the 0.1

5. CCSD(T) CONTRIBUTION To estimate CCSD(T) correlation effects, CCSD(T) single point calculations at the hybrid MP2:PBE+D2 structures are performed on the smaller CO/H2O·Mg9O9 cluster models only. We could afford D- and T-zeta basis sets only (aug’-ccp(C)VDZ and aug’-cc-p(C)VTZ). Table 10 shows the BSSEcorrected MP2, CCSD, and CCSD(T) adsorption energies for CO(1/8)·MgO(001) and H2O(1/8)·MgO(001), together with the energy differences, ΔCCSD = CCSD − MP2, Δ(T) = CCSD(T) − CCSD, and ΔCC = CCSD(T) − MP2. For CO, the present CCSD(T) calculations are a significant computational improvement compared to a previous ΔCC estimate,4 −0.1 ± 0.3 kJ/mol, which was obtained using smaller basis sets (D-zeta only) and smaller cluster models that are charged and, hence, require point charge embedding. The CCSD(T) correction to the MP2 value is −1.6 kJ/mol, which K


Article


Table 11. Adsorption Energy in kJ/mol, Mg2+···CO and (C−O)ads Bond Distances for (i) Hybrid QM:QM Method Including CPC-CBS Scheme, (ii) Periodic MP2 Estimate, and (iii) Final Estimate Including CCSD(T) Correction R(Mg ···CO) ΔEHL:LL,CPC(pbc) f · ΔLRb MP2(pbc) est. ΔCC Final estimate Periodic LMP2+ΔCC Exp.4,14

MP2:PBE+D2, ref 4


MP2/CPC-CBS(D,T):B3LYP+D2

Average

248.0 −17.6

−20.6 ± 2.4

−20.6 ± 2.4

249.4 −18.5 −2.3 −20.0 ± 0.2c,d −1.6 ± 0.1 −21.7 ± 0.3d,f −22.2 ± 1.1g −20.6 ± 2.4

250.8 ± 1.4a

−20.9 ± 0.7 −0.1 ± 0.3 −21.0 ± 1.0e

252.2 −16.3 −1.7 −19.2 ± 0.4c,d −1.6 ± 0.1 −20.8 ± 0.5d,f

2+

−19.6 ± 0.4a,c −1.6 ± 0.1 −21.2 ± 0.5a,f −20.6 ± 2.4

a

Average of PBE+D2 and B3LYP+D2 as low-level methods. bScaling of long-range correction, see Table 4. cIncludes the scaled long-range correction, f · ΔLR. dIncludes the uncertainty range coming from the low-level method which is estimated ±0.4 kJ/mol for PBE+D2 and ±0.2 kJ/ mol for B3LYP+D2 by the long-range correction decrease, from −0.9 to −1.7 kJ/mol and from −1.8 to −2.3 kJ/mol, respectively, upon scaling. Uncertainty coming from the basis correction is negligible (of the order of 0.05 kJ/mol). eIncludes the CCSD(T) correction (−0.1 ± 0.3 kJ/mol) of ref 4. fIncludes the CCSD(T) correction (−1.6 ± 0.1 kJ/mol) of Table 10 (MP2(pbc)est.+ΔCC). gAt the hybrid MP2/CPC-CBS(D,T):B3LYP +D2 equilibrium structure of CO(1/8)·MgO(001). It includes the CCSD(T) correction (−1.6 ± 0.1 kJ/mol) of Table 10 and the LMP2 relaxation energy (−1.0 kJ/mol).

Table 12. Zero-Point Vibrational Energy, ΔEZPV, Thermal Energy Change, ΔTΔE, Adsorption Enthalpy, ΔHT, at the Observed Temperature (203 K), and Derived Experimental Adsorption Energy of H2O on MgO(001), ΔEexp, in kJ/mol H2O(1/8)·MgO

H2O(1/16)·MgO

7.5 −0.8 5.1

9.1 −1.3 6.1

ΔEZPV Δ203ΔEa Δ203ΔE + ΔEZPV − R × 203 ΔEexpb

Average −50.2 ± 11.7

ΔH20361,62

5.6 ± 0.5 −55.8 ± 12.2

a

Thermal contribution to the internal energy at 203 K. bEquation 1.

Table 10 to the periodic MP2 estimate yields a final hybrid MP2:(DFT+D2)+ΔCC adsorption energy of −21.2 ± 0.5 kJ/ mol. The deviation from the experiment is −0.6 ± 2.9 kJ/mol, well within chemical accuracy limits. We also stress that the DFT-D contributions that survive in the final result (scaled long-range contribution) amount to just 2.0 ± 0.3 kJ/mol, only 10% of the total adsorption energy. Previous hybrid MP2:PBE+D2 calculations4 yielded a distance of 248 pm and a final MP2 estimate of −20.9 ± 0.7 kJ/mol. This includes the effect of CPC on the forces during structure optimizations (4 pm elongation of the distance; Table 2), a different CBS scheme, and a different cluster size extrapolation scheme resulting in 1.7 kJ/mol stronger binding. Since the present estimate of the CC contribution is 1.5 kJ/ mol more binding, the final hybrid MP2:(PBE+D2)+ΔCC estimates deviate only 0.2 kJ/mol from each other. This shows some compensation between different contributions in different computational hybrid protocols (both reasonable choices), but the effects are well within chemical accuracy limits. The LMP2 result of Table 8, after including the relaxation energy (−1.0 kJ/mol) and the ΔCC contribution (−1.6 ± 0.1 kJ/mol), yields a final adsorption energy of −22.2 ± 1.1 kJ/ mol, which also is well within chemical accuracy limits of the experiment (−1.6 ± 3.5 kJ/mol deviation). Among the DFT+D2 methods (Table 2 and Figure 2), PBE +D2 adsorption energies are only 1.5 kJ/mol more binding than the experimental value, but compared to the hybrid MP2:PBE+D2 calculation, the Mg2+···CO distance is 10 pm shorter. In contrast, B3LYP+D2 misses 3.7 kJ/mol on binding but yields a very good distance, which is only 1 pm shorter than the hybrid result. 6.2. H2O on MgO(001). For the adsorption of H2O on MgO(001) at low coverage, direct measurements of the H2O−

results from a positive ΔCCSD (2.3 kJ/mol) and a negative value of the triple substitution contribution Δ(T) (−4.0 kJ/ mol). Also for H2O adsorption, a positive value of ΔCCSD (2.6 kJ/mol) partly compensates the larger negative value of the triple contribution Δ(T) (−5.3 kJ/mol) resulting in a total CCSD(T) contribution of −2.7 kJ/mol. Finally, uncertainty ranges of ±0.1 and ±0.3 kJ/mol for CO and H2O adsorption, respectively, are estimated from the difference between CBS(D,T) and T-zeta CCSD(T) results.

6. COMPARISON WITH EXPERIMENT 6.1. CO on MgO(001). Experimental values for the CO· MgO(001) adsorption energy are available from temperatureprogrammed desorption (TPD) experiments for a thin MgO(001) film.14 From the observed Arrhenius desorption barriers, an experimental adsorption energy of −20.6 ± 2.4 kJ/ mol was previously obtained by extrapolating to the low coverage considered (Θ = 1/8) and properly correcting for zero-point vibrational energy (ZPVE) and thermal effects (eq 1).4 The specified uncertainty for the zero-point vibrational energy (0.4 kJ/mol)4 covers also anharmonicity effects which are just −0.15 kJ/mol.18 The calculation has been performed for CO adsorption on a metal−organic framework that features the same Mg2+(O2−)5 sites as the MgO(001) terraces.18 Table 11 shows the final hybrid results for PBE+D2 and B3LYP+D2 as low-level methods. Comparison of the longrange corrections (smaller for PBE+D2) and scaling factors f (smaller for B3LYP+D2) does not suggest preference for one or the other low-level method. Therefore, average values are considered as final hybrid estimates. This yields a Mg2+···CO distance of 251 ± 1 pm and an MP2 adsorption energy of −19.6 ± 0.4 kJ/mol. Adding the CCSD(T) contribution of L


Article


Table 13. Adsorption Energy, in kJ/mol, and Mg2+···O and O2−···H(OH) Bond Distances for (i) Hybrid QM:QM Method Including CPC-CBS Scheme, (ii) Periodic MP2 Estimate, and (iii) Final Estimate Including CCSD(T) Correction

R1(Mg2+···O) R2(O2−···H(OH)) ΔEHL:LL,CPC(pbc) f · ΔLRa MP2(pbc) est. ΔCC Final estimate Periodic LMP2+ΔCC Exp.

H2O(1/8)·MgO

H2O(1/16)·MgO


MP2/CPC-CBS(D,T):PBE+D2[Mg]

220.7 179.7 −48.0 −1.2 −51.0 ± 0.9b −2.7 ± 0.3 −53.7 ± 1.2b,c −55.4 ± 1.3d −55.8 ± 12.2e

221.3 180.5 −52.2 −1.5 −53.4 ± 0.9b −2.7 ± 0.3 −56.2 ± 1.2b,c

a

Scaled long-range correction, see eqs 16 and 17. bIncludes the basis set extension correction (Δbasis, see Tables S5 and S6 in the Supporting Information), the scaled long-range correction (f · ΔLR), and the corresponding uncertainties (±0.9 kJ/mol). cIncludes the CCSD(T) correction (−2.7 ± 0.3 kJ/mol) of Table 10 (MP2(pbc) est.+ΔCC). dAt the hybrid MP2/CPC-CBS(D,T):PBE+D2 equilibrium structure of H2O(1/8)· MgO(001). It includes the CCSD(T) correction (−2.7 ± 0.3 kJ/mol) of Table 10 and the LMP2 relaxation energy (−9.8 kJ/mol). eExperimental reference energy (see Table 12).

energy, −55.4 ± 1.3 kJ/mol, by 1.7 ± 2.5 kJ/mol only. The latter is obtained from the LMP2 result of Table 8 after including the relaxation energy (−9.8 kJ/mol) and the ΔCC correction of Table 10 (−2.7 ± 0.3 kJ/mol). Our predicted final hybrid MP2:(DFT+D2)+ΔCC estimate (−53.7 ± 1.2 kJ/mol) is expected to agree within chemical accuracy limits (±4.2 kJ/mol including the 2.0 kJ/mol CCSD(T) accuracy) with the experiment. Our final prediction for the adsorption energy of isolated H2O molecules is therefore −53.7 ± 4.2 kJ/mol, which falls into the range of values derived from experiments for low loadings (−55.8 ± 12.2 kJ/mol). Since the latter rely on assumptions about the magnitude of lateral interactions, we feel that our calculations provide a more accurate and reliable benchmark value for the energy of H2O adsorption on terrace sites of the MgO(001) surface for the limit of zero coverage. We note that the explicit PBE+D2 contribution in the final result (scaled long-range contribution) is a mere 1.2 ± 0.2 kJ/mol, only 2% of the total adsorption energy. Diffusion quantum Monte Carlo (DMC) calculations for an isolated H2O molecule on a two-layer slab model of MgO(001) yielded a binding of −46.4 kJ/mol.63 The calculations were performed with rigid structures of the MgO surface and the H2O molecule. MP2 cluster calculations with these rigid structures have also been performed, and an adsorption energy of −44.4 kJ/mol was reported as best estimate.63 We have performed full structure optimizations at the hybrid MP2/PBE+D2 level, and our best estimate from MP2 cluster calculations is −50.0 kJ/mol. This value results from CPC-CBS(D,T) calculations on H2O/Mg9O9, (Table 6) a contribution of −4.5 kJ/mol for increasing the model size to H2O/Mg25O25 (Table 7), and a contribution of −1.2 kJ/mol for extending the basis sets to CPC-CBS(Q,5) (eq 16). Most of the previous computational studies of H2O on MgO(001) were performed using DFT with some account of dispersion,49,55,56,64,65 but recently, single point calculations at PBE equilibrium structures using the random-phase approximation (RPA) have also been also reported.66 The difference of 11 kJ/mol between standard RPA (−47.5 kJ/mol exchange and correlation are computed nonself-consistently from PBE orbitals) and hybrid RPA (−58.7 kJ/mol, includes selfconsistent Fock exchange contributions) indicates a nonnegligible effect of the choice of the method. Compared to our

MgO(001) surface interaction are not available. Only an adsorption enthalpy with a known but large error bar (−50.2 ± 11.7 kJ/mol) has been derived as the difference between the heat of adsorption of H2O monolayers on the Mg(001) surface (−85.3 ± 2.1 kJ/mol) and the lateral molecule−molecule interactions in condensed H2O phases (−35.1 ± 9.6 kJ/mol), both obtained from LEED adsorption isotherms.61,62 An “experimental” adsorption energy, ΔEexp, is obtained from the estimated experimental adsorption enthalpy, ΔHT, according to eq 1. The zero-point vibrational energies, ΔEZPV, and the thermal energy changes, ΔTΔE, are computed for the PBE+D2 minimum structures with PBE+D2 and within the harmonic approximation. From the two sets of calculations, an average experimental reference energy of −55.8 ± 12.2 kJ/mol is obtained (Table 12). Table 13 shows the final hybrid MP2:(DFT+D2)+ΔCC estimates for the low coverage systems H2O(1/8)·MgO(001) and H2O(1/16)·MgO(001), −53.7 and −56.2 kJ/mol, respectively. The difference to the corresponding low-level calculations for the former is smaller (+1.3 kJ/mol, PBE+D2) than that for the latter (+4.9 kJ/mol, PBE+D2[Mg]). This provides additional evidence that using the Ne dispersion parameters for the isoelectronic Mg2+ cation3 is superior compared to using the parameters for the Mg atom also for the Mg2+ ion. In further discussions, we will focus on the hybrid MP2:(PBE+D2)+ΔCC results for the H2O(1/8)·MgO(001) model. To estimate the remaining uncertainties of our final energies with respect to the CCSD(T) limit, we consider the following effects. The CBS extrapolations of the BSSE-corrected and uncorrected values should converge to the same limiting value. The CBS(Q,5) and CPC-CBS(Q,5) results show a difference of 1.4 kJ/mol (Table S6). Hence, this range, ±0.7 kJ/mol, is taken as uncertainty of our MP2 results. The uncertainty of the long-range correction, ± 0.2 kJ/mol, we derive from the decrease in the long-range correction upon scaling, from −0.8 to −1.2 kJ/mol. Finally, an uncertainty range of ±0.3 kJ/mol is derived for the coupled cluster correction, ΔCC, from the difference between CBS(D,T) and T-zeta CCSD(T) results. The sum of these uncertainty values, ± 1.2 kJ/mol, when added to the PBE+D2-based hybrid result (−53.7 kJ/mol), yields a final CCSD(T) estimate of −53.7 ± 1.2 kJ/mol for the adsorption energy. It deviates from the final LMP2+ΔCC M


Article


Table 14. Adsorption Energies, −ΔE (kJ/mol), Obtained with Different Methods/Models Compared to Energies Derived from Experimental Enthalpies PBE+D2 PBE+D2 MP2/CPC-CBS(D,T) Best MP2 estimate ΔCC Final QM:QM estimate Experiment

Model

Structure

CH4a

CO

H2O

H2Ob

pbc Mg25O25 Mg25O25 pbc

optimized single pointc single pointc optimizedc

14.8

22.1 21.5 18.9 19.2 (1.6) 21.2 ± 0.5 20.6 ± 2.4

55.0 60.5 60.1 51.0 (2.7) 53.7 ± 1.2

61.0b 66.7b 57.7 53.4 (2.7) 56.2 ± 1.2 55.8 ± 12.2

11.4 (2.6) 14.0 ± 1.0 15.0 ± 0.6

a

Reference 10. bOriginal Mg parameters of ref 39 are used for Mg2+. cOptimization performed at MP2/CPC-CBS(D,T):PBE+D2 level.

reference value, the former result deviates by +6.2 kJ/mol, while the latter deviates by −5.0 kJ/mol.

From the point of view of application of standard software, Table 14 compares our final hybrid QM:QM estimates with the results obtained with standard approaches such as DFT-D with periodic boundary conditions (PBE+D2 in our case) and DFT-D and MP2 calculations on large cluster models (Mg25O25 in our case). We note that the latter are single point calculations that could be performed at the DFT-D structure but in our case have been performed at the hybrid MP2:PBE+D2 structure (Mg 9O 9 cluster for the MP2 calculations). We have also included results from a previous study on a monolayer of CH4 on MgO(001). For adsorption energies, PBE+D2 with periodic boundary conditions seems to perform very well for the systems studied, in particular if Ne atom dispersion parameters are used for the Mg2+ ions, but this is not true for the molecule−surface distances. In addition to the previous validation studies,3,4,10,15−20 this work provides further evidence for the hybrid MP2:(DFT +D2)+ΔCC scheme being a reliable and robust tool for structure and energy predictions for molecule−surface interactions. The hybrid scheme presented is general and scalable, as the cluster size can be enlarged and the high-level method improved when improved methods and algorithms for large systems become available, see, e.g., recent CCSD(T) developments by Neese and co-workers28,67,68 or Werner and co-workers.69

7. SUMMARY AND CONCLUSIONS A hybrid MP2:(DFT+D2)+ΔCC scheme is applied to adsorption of CO and H2O on terrace sites of the MgO(001) surface at low coverage. Hybrid MP2/CPCCBS(D,T):DFT+D2 structure optimizations are performed on BSSE-free and CBS-extrapolated potential energy surfaces using the MonaLisa code.9 Optimization at the hybrid level can have large effects on molecule−surface distances. For CO on MgO(001), with hybrid MP2/CPC-CBS:PBE+D2 optimization, the Mg2+···CO distance gets 11 pm larger than with PBE +D2. Structure relaxation on adsorption is found to be crucial for reaching chemical accuracy. For H2O on MgO(001), the relaxation energy (−11.5 kJ/mol) is 20% of the adsorption energy (−53.7 kJ/mol). Subsequent to the hybrid MP2:DFT+D2 optimization, two types of single point calculations are performed: (i) MP2CPC-CBS(Q,5) calculations to estimate the effect of increasing the basis set size in the CBS-extrapolation and (ii) MP2 and DFT calculations on larger cluster models (Mg25O25 compared to Mg9O9) providing scaling factors for the long-range correction. The resulting estimate of the periodic MP2 limit is in close agreement, within 1 kJ/mol, with the LMP2 calculations using periodic boundary conditions (Table 9). After performing CCSD(T) calculations with CBS(D,T) extrapolation, final hybrid MP2:(DFT+D2)+ΔCC results are obtained. For both CO and H2O on MgO(001), the lion’s share of the final adsorption energy comes from the MP2 calculations. For CO, out of the −21.2 kJ/mol in total, −1.6 kJ/mol are CCSD(T) contributions, and only −2.0 kJ/mol are long-range corrections evaluated at the DFT+D2 level. Similarly, for H2O, out of the −53.7 kJ/mol, −2.7 kJ/mol are CCSD(T) contributions, and as little as −1.2 kJ/mol are PBE+D2 long-range corrections. This shows that for the systems studied MP2 is a very good approximation and that cluster models can be made large enough so that effects of the low-level method are minor. Taking the estimated uncertainty due to the specific computational protocol into account, the result for CO/ MgO(001), −21.2 ± 0.5 kJ/mol, agrees well within chemical accuracy limits with the experimentally derived adsorption energy of −20.6 ± 2.4 kJ/mol. This lends further credit to our prediction of −53.7 ± 4.2 kJ/mol for the adsorption energy for H 2 O/MgO (low coverage limit) that has not been experimentally determined. It falls in the range of values, −55.8 ± 12.2 kJ/mol, derived from a high coverage LEED experiment and estimated lateral interactions.61,62

■

ASSOCIATED CONTENT

S Supporting Information *

This information is available free of charge via the Internet at The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.8b01122. Comparison between plane wave and Gaussian basis sets, basis set extension effects, Ne vs Mg parameters for Mg2+, relaxation of surface atoms, computational details of LMP2 calculations, computer time for hybrid MP2:DFT+D2 vs LMP2 calculations, structure files (VASP format), and total energies that enter the hybrid MP2:DFT+D2 scheme. (PDF)

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (D. Usvyat). *E-mail: [email protected] (J. Sauer). ORCID

Joachim Sauer: 0000-0001-6798-6212 Notes

The authors declare no competing financial interest. N


Article


■

(17) Sillar, K.; Sauer, J. Ab Initio Prediction of Adsorption Isotherms for Small Molecules in Metal−Organic Frameworks: The Effect of Lateral Interactions for Methane/CPO-27-Mg. J. Am. Chem. Soc. 2012, 134, 18354−18365. (18) Kundu, A.; Piccini, G.; Sillar, K.; Sauer, J. Ab Initio Prediction of Adsorption Isotherms for Small Molecules in Metal-Organic Frameworks. J. Am. Chem. Soc. 2016, 138, 14047−14056. (19) Sillar, K.; Kundu, A.; Sauer, J. Ab Initio Adsorption Isotherms for Molecules with Lateral Interactions: CO2 in Metal−Organic Frameworks. J. Phys. Chem. C 2017, 121, 12789−12799. (20) Piccini, G.; Alessio, M.; Sauer, J. Ab Initio Calculation of Rate Constants for molecule−surface Reactions with Chemical Accuracy. Angew. Chem., Int. Ed. 2016, 55, 5235−5237. (21) Kim, Y. D.; Stultz, J.; Goodman, D. W. Dissociation of water on MgO(100). J. Phys. Chem. B 2002, 106, 1515−1517. (22) Yu, Y.; Guo, Q.; Liu, S.; Wang, E.; Møller, P. J. Partial dissociation of water on a MgO(100) film. Phys. Rev. B: Condens. Matter Mater. Phys. 2003, 68, 115414. (23) Bakowies, D.; Thiel, W. Hybrid models for combined quantum mechanical and molecular mechanical approaches. J. Phys. Chem. 1996, 100, 10580−10594. (24) Stoll, H. Correlation energy of diamond. Phys. Rev. B: Condens. Matter Mater. Phys. 1992, 46, 6700−6704. (25) Ugliengo, P.; Damin, A. Are dispersive forces relevant for CO adsorption on the MgO(001) surface? Chem. Phys. Lett. 2002, 366, 683−690. (26) Nolan, S. J.; Gillan, M. J.; Alfè, D.; Allan, N. L.; Manby, F. R. Calculation of properties of crystalline lithium hydride using correlated wave function theory. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 165109. (27) Usvyat, D.; Sadeghian, K.; Maschio, L.; Schütz, M. Geometrical frustration of an argon monolayer adsorbed on the MgO (100) surface: An accurate periodic ab initio study. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 045412. (28) Kubas, A.; Berger, D.; Oberhofer, H.; Maganas, D.; Reuter, K.; Neese, F. Surface Adsorption Energetics Studied with “Gold Standard” Wave-Function-Based Ab Initio Methods: Small-Molecule Binding to TiO2(110). J. Phys. Chem. Lett. 2016, 7, 4207−4212. (29) Manby, F. R.; Stella, M.; Goodpaster, J. D.; Miller, T. F. A Simple, Exact Density-Functional-Theory Embedding Scheme. J. Chem. Theory Comput. 2012, 8, 2564−2568. (30) Usvyat, D.; Maschio, L.; Schütz, M. Periodic and fragment models based on the local correlation approach. WIREs Comput. Mol. Sci. 2018, 8, No. e1357. (31) Bennie, S. J.; van der Kamp, M. W.; Pennifold, R. C. R.; Stella, M.; Manby, F. R.; Mulholland, A. J. A Projector-Embedding Approach for Multiscale Coupled-Cluster Calculations Applied to Citrate Synthase. J. Chem. Theory Comput. 2016, 12, 2689−2697. (32) Boys, S. F.; Bernardi, F. Calculation of Small Molecular Interactions by Differences of Separate Total Energies - Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553−566. (33) Dunning, T. H. A road map for the calculation of molecular binding energies. J. Phys. Chem. A 2000, 104, 9062−9080. (34) Tosoni, S.; Sauer, J.; Civalleri, B.; Ugliengo, P.; Tuma, C. A Comparison between Plane Wave and Gaussian-type Orbital Basis Sets for Hydrogen Bonded Systems: Formic Acid as a Test Case. J. Chem. Phys. 2007, 127, 154102. (35) Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169−11186. (36) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (37) Lee, C. T.; Yang, W. T.; Parr, R. G. Development of the ColleSalvetti Correlation-Energy Formula into a Functional of the Electron-Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (38) Becke, A. D. Density functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652.

ACKNOWLEDGMENTS Dedicated to the memory of Keiji Morokuma, a pioneer in QM:QM methodology. This work has been supported by German Research Foundation (DFG) within a Reinhart Koselleck grant to J.S. and by the “Fonds der Chemischen Industrie”. M.A. has been a member of the International MaxPlanck Research School “ Complex Surfaces in Materials Science”.

■

REFERENCES

(1) Tuma, C.; Sauer, J. A hybrid MP2/planewave-DFT scheme for large chemical systems: Proton jumps in zeolites. Chem. Phys. Lett. 2004, 387, 388−394. (2) Tuma, C.; Sauer, J. Treating dispersion effects in extended systems by hybrid MP2: DFT calculations - protonation of isobutene in zeolite ferrierite. Phys. Chem. Chem. Phys. 2006, 8, 3955−3965. (3) Tosoni, S.; Sauer, J. Accurate quantum chemical energies for the interaction of hydrocarbons with oxide surfaces: CH4/MgO(001). Phys. Chem. Chem. Phys. 2010, 12, 14330−14340. (4) Boese, A. D.; Sauer, J. Accurate adsorption energies of small molecules on oxide surfaces: CO-MgO(001). Phys. Chem. Chem. Phys. 2013, 15, 16481−16493. (5) Humbel, S.; Sieber, S.; Morokuma, K. The IMOMO method: Integration of different levels of molecular orbital approximations for geometry optimization of large systems: Test for n-butane conformation and SN2 reaction: RCl+Cl−. J. Chem. Phys. 1996, 105, 1959−1967. (6) Chung, L. W.; Sameera, W. M. C.; Ramozzi, R.; Page, A. J.; Hatanaka, M.; Petrova, G. P.; Harris, T. V.; Li, X.; Ke, Z. F.; Liu, F. Y.; Li, H. B.; Ding, L. N.; Morokuma, K. The ONIOM Method and Its Applications. Chem. Rev. 2015, 115, 5678−5796. (7) Eichler, U.; Kölmel, C. M.; Sauer, J. Combining ab initio techniques with analytical potential functions for structure predictions of large systems: Method and application to crystalline silica polymorphs. J. Comput. Chem. 1997, 18, 463−477. (8) Sierka, M.; Sauer, J. Finding Transition Structures in Extended Systems: A Strategy Based on a Combined Quantum Mechanics − Empirical Valence Bond Approach. J. Chem. Phys. 2000, 112, 6983− 6996. (9) Bischoff, F. A.; Alessio, M.; John, M.; Rybicki, M.; Sauer, J. Multi-Level Energy Landscapes: The MonaLisa Program.; HumboldtUniversity of Berlin, 2017. https://www.chemie.hu-berlin.de/de/ forschung/quantenchemie/monalisa/ (accessed December 2018). (10) Alessio, M.; Bischoff, F. A.; Sauer, J. Chemically accurate adsorption energies for methane and ethane monolayers on the MgO(001) surface. Phys. Chem. Chem. Phys. 2018, 20, 9760−9769. (11) Pisani, C.; Maschio, L.; Casassa, S.; Halo, M.; Schütz, M.; Usvyat, D. Periodic local MP2 method for the study of electronic correlation in crystals: Theory and preliminary applications. J. Comput. Chem. 2008, 29, 2113−2124. (12) Zheng, J. J.; Zhao, Y.; Truhlar, D. G. The DBH24/08 Database and Its Use to Assess Electronic Structure Model Chemistries for Chemical Reaction Barrier Heights. J. Chem. Theory Comput. 2009, 5, 808−821. (13) Sylvetsky, N.; Peterson, K. A.; Karton, A.; Martin, J. M. L. Toward a W4-F12 approach: Can explicitly correlated and orbitalbased ab initio CCSD(T) limits be reconciled? J. Chem. Phys. 2016, 144, 214101. (14) Dohnalek, Z.; Kimmel, G. A.; Joyce, S. A.; Ayotte, P.; Smith, R. S.; Kay, B. D. Physisorption of CO on the MgO(100) surface. J. Phys. Chem. B 2001, 105, 3747−3751. (15) Piccini, G.; Alessio, M.; Sauer, J.; Zhi, Y. C.; Liu, Y.; Kolvenbach, R.; Jentys, A.; Lercher, J. A. Accurate Adsorption Thermodynamics of Small Alkanes in Zeolites. Ab initio Theory and Experiment for H-Chabazite. J. Phys. Chem. C 2015, 119, 6128−6137. (16) Sillar, K.; Hofmann, A.; Sauer, J. Ab Initio Study of Hydrogen Adsorption in MOF-5. J. Am. Chem. Soc. 2009, 131, 4143−4150. O


Article

Journal of Chemical Theory and Computation (39) Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 1787−1799. (40) Kerber, T.; Sierka, M.; Sauer, J. Application of semiempirical long-range dispersion corrections to periodic systems in density functional theory. J. Comput. Chem. 2008, 29, 2088−2097. (41) Hättig, C.; Weigend, F. CC2 excitation energy calculations on large molecules using the resolution of the identity approximation. J. Chem. Phys. 2000, 113, 5154−5161. (42) TURBOMOLE V6.5 2013, a development of University of Karlsruhe and Forschungzentrum Karlsruhe GmbH, 1989−2007, TURBOMOLE GmbH, since 2007, available from http://www. turbomole.com (accessed December 2018). (43) Hättig, C.; Hellweg, A.; Kohn, A. Distributed memory parallel implementation of energies and gradients for second-order MøllerPlesset perturbation theory with the resolution-of-the-identity approximation. Phys. Chem. Chem. Phys. 2006, 8, 1159−1169. (44) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. ElectronAffinities of the 1st-Row Atoms Revisited - Systematic Basis-Sets and Wave-Functions. J. Chem. Phys. 1992, 96, 6796−6806. (45) Prascher, B. P.; Woon, D. E.; Peterson, K. A.; Dunning, T. H.; Wilson, A. K. Gaussian basis sets for use in correlated molecular calculations. VII. Valence, core-valence, and scalar relativistic basis sets for Li, Be, Na, and Mg. Theor. Chem. Acc. 2011, 128, 69−82. (46) Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. (47) Solans-Monfort, X.; Branchadell, V.; Sodupe, M.; Sierka, M.; Sauer, J. Electron hole formation in acidic zeolite catalysts. J. Chem. Phys. 2004, 121, 6034−6041. (48) Ončaḱ , M.; Włodarczyk, R.; Sauer, J. Hydration Structures of MgO, CaO, and SrO (001) Surfaces. J. Phys. Chem. C 2016, 120, 24762−24769. (49) Giordano, L.; Goniakowski, J.; Suzanne, J. Partial dissociation of water molecules in the (3 × 2) water monolayer deposited on the MgO (100) surface. Phys. Rev. Lett. 1998, 81, 1271−1273. (50) Odelius, M. Mixed molecular and dissociative water adsorption on MgO[100]. Phys. Rev. Lett. 1999, 82, 3919−3922. (51) Delle Site, L.; Alavi, A.; Lynden-Bell, R. M. The structure and spectroscopy of monolayers of water on MgO: An ab initio study. J. Chem. Phys. 2000, 113, 3344−3350. (52) Lynden-Bell, R. M.; Delle Site, L.; Alavi, A. Structures of adsorbed water layers on MgO: an ab initio study. Surf. Sci. 2002, 496, L1−L6. (53) Włodarczyk, R.; Sierka, M.; Kwapień, K.; Sauer, J.; Carrasco, E.; Aumer, A.; Gomes, J. F.; Sterrer, M.; Freund, H. J. Structures of the ordered water monolayer on MgO(001). J. Phys. Chem. C 2011, 115, 6764−6774. (54) Ončaḱ , M.; Włodarczyk, R.; Sauer, J. Water on the MgO(001) Surface: Surface Reconstruction and Ion Solvation. J. Phys. Chem. Lett. 2015, 6, 2310−2314. (55) Carrasco, J.; Illas, F.; Lopez, N. Dynamic Ion Pairs in the Adsorption of Isolated Water Molecules on Alkaline-Earth Oxide (001) Surfaces. Phys. Rev. Lett. 2008, 100, 016101. (56) Hu, X. L.; Carrasco, J.; Klimeŝ, J.; Michaelides, A. Trends in water monomer adsorption and dissociation on flat insulating surfaces. Phys. Chem. Chem. Phys. 2011, 13, 12447−12453. (57) Dovesi, R.; Orlando, R.; Erba, A.; Zicovich-Wilson, C. M.; Civalleri, B.; Casassa, S.; Maschio, L.; Ferrabone, M.; De La Pierre, M.; D’Arco, P.; Noël, Y.; Causà, M.; Rérat, M.; Kirtman, B. CRYSTAL14: A program for the ab initio investigation of crystalline solids. Int. J. Quantum Chem. 2014, 114, 1287−1317. (58) Usvyat, D. High precision quantum-chemical treatment of adsorption: Benchmarking physisorption of molecular hydrogen on graphane. J. Chem. Phys. 2015, 143, 104704. (59) Usvyat, D. Linear-scaling explicitly correlated treatment of solids: Periodic local MP2-F12 method. J. Chem. Phys. 2013, 139, 194101.

(60) Pisani, C.; Schütz, M.; Casassa, S.; Usvyat, D.; Maschio, L.; Lorenz, M.; Erba, A. Cryscor: a program for the post-Hartree-Fock treatment of periodic systems. Phys. Chem. Chem. Phys. 2012, 14, 7615−7628. (61) Ferry, D.; Glebov, A.; Senz, V.; Suzanne, J.; Toennies, J. P.; Weiss, H. The properties of a two-dimensional water layer on MgO(001). Surf. Sci. 1997, 377−379, 634−638. (62) Ferry, D.; Picaud, S.; Hoang, P. N. M.; Girardet, C.; Giordano, L.; Demirdjian, B.; Suzanne, J. Water monolayers on MgO(100): structural investigations by LEED experiments, tensor LEED dynamical analysis and potential calculations. Surf. Sci. 1998, 409, 101−116. (63) Karalti, O.; Alfè, D.; Gillan, M. J.; Jordan, K. D. Adsorption of a water molecule on the MgO(100) surface as described by cluster and slab models. Phys. Chem. Chem. Phys. 2012, 14, 7846−7853. (64) Jung, J.; Shin, H.-J.; Kim, Y.; Kawai, M. Controlling water dissociation on an ultrathin MgO film by tuning film thickness. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 085413. (65) Kebede, G. G.; Spångberg, D.; Mitev, P. D.; Broqvist, P.; Hermansson, K. Comparing van der Waals DFT methods for water on NaCl(001) and MgO(001). J. Chem. Phys. 2017, 146, 064703. (66) Bajdich, M.; Nørskov, J. K.; Vojvodic, A. Surface energetics of alkaline-earth metal oxides: Trends in stability and adsorption of small molecules. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 91, 155401. (67) Riplinger, C.; Neese, F. An efficient and near linear scaling pair natural orbital based local coupled cluster method. J. Chem. Phys. 2013, 138, 034106. (68) Riplinger, C.; Sandhoefer, B.; Hansen, A.; Neese, F. Natural triple excitations in local coupled cluster calculations with pair natural orbitals. J. Chem. Phys. 2013, 139, 134101. (69) Ma, Q.; Werner, H.-J. Scalable Electron Correlation Methods. 5. Parallel Perturbative Triples Correction for Explicitly Correlated Local Coupled Cluster with Pair Natural Orbitals. J. Chem. Theory Comput. 2018, 14, 198−215.

P


Chemically Accurate Adsorption Energies: CO and H2O on the MgO

Recommend Documents