Periodic Density Matrix Embedding for CO Adsorption on the MgO(001) Surface.

The adsorption of simple gas molecules to metal oxide surfaces is a primary step in many heterogeneous catalysis applications. Quantum chemical modeling of these reactions is a challenge in terms of both cost and accuracy, and quantum-embedding methods are promising, especially for localized chemical phenomena. In this work, we employ density matrix embedding theory (DMET) for periodic systems to calculate the adsorption energy of CO to the MgO(001) surface. Using coupled-cluster theory with single and double excitations and second-order Møller-Plesset perturbation theory as quantum chemical solvers, we perform calculations with embedding clusters up to 266 electrons in 306 orbitals, with the largest embedding models agreeing to within 1.2 kcal/mol of the non-embedding references. Moreover, we present a memory-efficient procedure of storing and manipulating electron repulsion integrals in the embedding space within the framework of periodic DMET.

Magnesium oxide (MgO) surface plays an important role in several heterogeneous catalytic reactions, such as the partial oxidation of methane,[1] the Guerbet reaction at a low pressure,[2] the synthesis of 2-amino-2-chromenes using benign reactants,[3] the conversion of ethane to ethylene,[4] and the formation of carbonates from carbon monoxide in the presence of oxygen.[5] Modeling surface adsorption of simple molecules, for example, carbon monoxide (CO), to metal oxide surfaces, like MgO, is an important step for theorists toward the understanding of heterogeneous catalysis, but it is challenging.[6−9] The CO molecule binds to the MgO surface, preferentially with a C–Mg interaction.[10] The CO/MgO adsorption energy is relatively small, and a range of values has been obtained by different experimental techniques and theoretical methods. An adsorption energy of 3.23 kcal/mol was obtained from thermodesorption experiments by Wichtendahl et al.,[11] whereas temperature-programmed desorption (TPD) experiments performed by Dohnálek et al.[12] accounted for an adsorption energy of 4.84 kcal/mol. An experimental study by Xu et al. reported an interaction energy of 3.0 kcal/mol.[13] For a more extensive description of the rich experimental history of the MgO/CO adsorption, readers are referred to the review by Spoto et al.[14]

Computationally, the challenge posed by this system is the weak interaction between CO and the surface, mainly arising from van der Waals (vDW) forces. Many local and semi-local Kohn–Sham density functionals[15,16] are unable to account for vDW interactions in such cases.[17−20] The accurate estimation of the adsorption energy, therefore, requires an extensive testing of density functional theory (DFT) functionals and the incorporation of dispersion corrections.[21−25] On the other hand, size-consistent correlated wave function (CWF)-based ab initio methods can model vDW interactions,[17] and, in the past few years, their application to periodic systems has gained momentum.[26−33] An attractive feature of CWF methods is their systematic improvability; however, their steep computational cost scaling with system size poses an obstacle.[17,34] This becomes apparent in applications where one cannot exploit translational symmetry as a result of the presence of irregularities in the crystal, like point defects or surface adsorbates.

Among the most recent wave function theoretical studies, Staemmler computed an adsorption energy of 2.86 kcal/mol using the method of local increments,[17] whereas, using their combined MP2-CCSD(T) approach embedded in a potential of point charges, Boese et al.[8] calculated an adsorption energy of 5.0 kcal/mol, but also pointed out the wide range of numbers that can be obtained using different electronic structure theories. Valero et al.[6] showed that the Minnesota functionals M06-2X and M06-HF provide adsorption energies of around 6 kcal/mol.

These systems are usually modeled by either cutting a cluster from an extended system (cluster modeling) or assuming periodic boundary conditions (PBCs). Defining a cluster involves choosing an appropriate cluster size and saturating the free valencies using hydrogen atoms, which can create spurious electronic states at the boundary. Previously used cluster models[8,9,35−39] were surrounded with point charges or periodic potentials to replicate the environment. On the other hand, modeling surface adsorption with PBCs using CWF becomes prohibitively costly as a result of the apparent need of large supercells (often hundreds of atoms). To overcome the cost and maintain the accuracy of the parent method, the models can be subjected to fragmentation/embedding approaches.[8,31,32,40,41] Quantum embedding methods use a high-level quantum chemistry solver to represent a small region of interest (here referred to as the fragment/impurity), whereas the rest of the system (generally referred to as “environment”) is represented using a mean-field method, such as Kohn–Sham density functional theory (KS-DFT)[15,16] or Hartree–Fock (HF).[34,42] Modeling the adsorption of CO to a MgO surface is therefore ideal to investigate the performance of wave function-in-wave function quantum embedding approaches.

In this work, we use the density matrix embedding theory (DMET) algorithm to calculate the adsorption energy of a CO molecule to the MgO(001) surface. DMET, a wave function-in-wave function embedding technique,[43] was originally proposed as a promising alternative to dynamical mean-field theory (DMFT)[44] to treat strongly correlated fermions in the one-dimensional Hubbard model. Several theoretical developments and targeted applications have followed since.[45−54] DMET uses the Schmidt decomposition[55] of a mean-field wave function to model the environment of a given impurity space using an effective bath. Pham et al.[56] and Cui et al.[57] extended the DMET algorithm to periodic systems. Here, we use the coupled-cluster theory with single and double excitations (CCSD) and second-order Møller–Plesset perturbation theory (MP2) as high-level solvers within the DMET formalism and compare them to non-embedding Γ-point CCSD and MP2 references.

The DMET calculations are performed with the periodic DMET (pDMET) code,[58,59] which uses the electron integrals and quantum chemical solvers from the PySCF package.[60,61] Similar to the workflow in ref (62), we first perform a HF calculation to obtain the mean-field wave function. Next, we define the impurity region using a set of localized orbitals in real space. Here, we use the maximally localized Wannier functions (MLWFs),[63,64] implemented in the wannier90[65] code via the pyWannier90 interface.[66] Because adsorbates (i.e., perturbations to the pristine crystal) are introduced, the Brillouin zone is sampled at the Γ point and a subset of the MLWFs (which we label as Nimp) at the chemical region of interest, for example, those around the adsorbate, are chosen to define the impurity. The bath is defined using the Schmidt decomposition,[46] where the environment block (Denv) of the one-body reduced density matrix (1-RDM) is diagonalized as follows:where λ is a diagonal matrix of eigenvalues λ (i = 0, 1, ..., Nenv, where Nenv is the number of the environment orbitals). The columns of the unitary matrix U corresponding to λ other than zero or two define the entangled bath orbitals; the remainder orbitals are treated as a frozen core in the embedding calculation. The mean-field wave function after the Schmidt decomposition thus has the following form:where , , and are single determinants in the Fock spaces of the Nimp impurity orbitals, Nbath ≤ Nimp bath orbitals, and Ncore frozen core orbitals, respectively, and i = 0, 1, ..., Nimp. The high-level wave function, , diagonalizes the embedding Hamiltonian, where is the partial trace of the Hamiltonian over the determinant; its operator terms involve only Nemb = Nimp + Nbath ≤ 2Nimp embedding orbitals. For calculations of energy differences, it is important to choose the same number of Nemb embedding orbitals for the different geometries. As discussed later, the computational cost of the high-level method is thus reduced by not requiring to have Ncore orbitals in the effective Hamiltonian.

We use a density fitting (DF) approach based on the Cholesky decomposition,[67−70] where four-center electron repulsion integrals (ERIs) in the embedding space can be reconstructed in terms of the three-center ERIs aswhere P and Q represent auxiliary basis functions, M = (P |Q) is the Coulomb metric, and and are the Cholesky vectorswhere the expansion coefficients are obtained by solving a linear equationThe auxiliary basis set contains even-tempered basis (ETB) functions generated with a progression factor β = 2.0 for the auxiliary expansion of the polarized double-ζ basis set (GTH-DZVP) and polarized triple-ζ basis set (GTH-TZVP) bases and will be represented using P and Q. The prefix GTH is used because these basis sets are consistent with the Goedecker–Teter–Hutter pseudopotentials[71,72] that have been used for all of the calculations. Our approach is agnostic to whether we use pseudopotentials or an all-electron basis sets because surface adsorption primarily depends upon valence electrons and orbitals.

In the previous implementations of periodic DMET,[56,57] the quantum impurity solvers used four-center two-electron integrals obtained by contracting the three-center electron repulsion integrals (ERIs) in the embedding space as in eq . This eliminated the full-basis 4-index ERI array but still required the storage of the 4-index ERIs in the embedding basis, whose memory cost scales with the size of the impurity as .[56,57] On the other hand, in the current implementation, we use the DF for the MP2 and CCSD high-level solvers, in which the programmable equations for the energy are implemented in terms of the Cholesky vectors themselves and the (ij|kl) integrals in the embedding basis are not required. In other words, the right-hand side of eq is not evaluated but is algebraically substituted into the energy expressions in the high-level solver implementations. The formal memory cost scaling of this approach (with respect to the size of the impurity) is , which, in practice, is much more favorable for our applications, as shown later.

We compute the adsorption energy, ΔE, as the difference between the energy at the equilibrium geometry, Eeq (C–Mg bond distance of 2.479 Å),[8] and at a separated geometry, Esep (C–Mg bond distance of 6 Å), as indicated in Figure a. The 4 × 4 × 2 slab model of MgO has a vacuum of approximately 16 Å in the vertical direction (collinear to CO) above the MgO surface to avoid the interaction between neighboring images. We consider four choices of impurity clusters, as shown in Figure b. For these four choices, the orbitals are localized on (a) only the CO molecule, (b) the CO molecule and the nearest Mg atom on the substrate (denoted as CO + Mg), (c) the CO molecule, the nearest Mg atom, and the 4 nearest O atoms on the substrate (denoted as CO + MgO4), and finally (d) the 4 next to nearest Mg atoms in addition to choice (c) (denoted as CO + Mg5O4). The orbitals localized at the highlighted atoms in the embedding clusters (Figure b) have been considered as the fragment. We do not correct for basis set superposition error (BSSE) in our calculations, because this would require the use of ghost basis functions. The Schmidt decomposition is unable to produce bath orbitals entangled to the unoccupied ghost orbitals; therefore, these correction calculations would be systematically deficient in bath orbitals compared to the calculation being corrected. A proper way to account for the most entangled orbitals from the environment is desired especially for physical/chemical phenomena, where BSSE is non-negligible and is currently an area upon which we are working.

Figure 1

(a) CO at a distance of 2.479 Å (left) representing the geometry at equilibrium (referred as eq) and 6 Å (referred as sep) from the MgO surface (right) representing the geometry when there is no interaction between the substrate and adsorbate. Magnesium (Mg) atoms are shown in red; oxygen (O) atoms are shown in blue; and carbon (C) atoms are shown in gray. (b) Atoms highlighted in yellow form the impurity clusters used for DMET calculations.

In Figure , we report the relative energy Erel of the CO + MgO model as a function of the distance between C (in CO) and Mg (in MgO) from 2 to 6 Å. We take as a reference value the total energy at the C–Mg distance of 2.479 Å, and Erel at all other geometries are reported relative to this reference. The results are obtained using periodic MP2 calculations, restricted HF (RHF), and DMET-MP2. The DMET-MP2 calculations are performed using the smaller CO + Mg impurity subspace and the larger CO + MgO4 impurity subspace.

Figure 2

Relative energies Erel (kcal/mol) obtained using non-embedding MP2 (blue diamond), DMET-MP2 with the CO + Mg impurity cluster (red circles), DMET-MP2 with the CO + MgO4 impurity cluster (dark blue circles) and RHF (purple crosses). The abscissa represents the Mg–C distances in angstroms. All Erel values are reported as differences with respect to the value at the C–Mg distance of 2.479 Å. All calculations are performed using the DZVP basis set.

For the MP2 reference method, Erel reaches an asymptotic value at a C–Mg distance of 6 Å and differs from Erel at 5 Å by only 0.3 kcal/mol, thereby suggesting that 6 Å is a reasonable choice for a separated geometry. Using the CO + MgO4 fragment, the DMET Erel values at each geometry are within 2 kcal/mol of the MP2 references. Using the CO + Mg fragment, the DMET value at the C–Mg bond distance of 2 Å has a large disagreement (ca. 5 kcal/mol) with the non-embedding reference, suggesting the importance of using a larger fragment space. The RHF Erel values deviate significantly from the MP2 references. Morover, the RHF Erel values at 3, 4, and 5 Å are negative, thereby indicating the presence of a minimum C–Mg bond length significantly away from the literature value of 2.479 Å.[8] This is consistent with cluster HF calculations by Nygren et al.[73] DMET on the other hand reproduces the binding energy (to within 1.5 kcal/mol) that is predicted by the reference.

Next, DMET calculations with the embedding clusters are compared to the periodic Γ-point CCSD and MP2 calculations (termed as the non-embedding references). The energy differences ΔE calculated using DMET-CCSD and DMET-MP2 with different basis sets are shown in Figure . The numbers are reported in Tables S2 and S3 of the Supporting Information. Four different basis set compositions have been used. They are divided as either DZVP on all of the atoms or TZVP on important atoms and DZVP on all of the others. TZVP (X) refers to TZVP applied on the X set of atoms and DZVP on all others.

Figure 3

Adsorption energies (ΔE) between the equilibrium (2.479 Å) and separated (6 Å) geometries calculated using different basis sets and impurity cluster models. TZVP (X) refers to X being treated with the TZVP basis set and the rest of the system using the DZVP basis set. The solid lines correspond to the periodic Γ-point CCSD/MP2 calculations with red/dark blue color coding, respectively.

With the two largest impurity clusters, there is a closer agreement with the non-embedding references. This suggests that a larger number of surface atoms in the DMET impurity space is necessary for better accuracy. In Figure S1 of the Supporting Information, we plot the mean absolute deviations (MADs) from the non-embedding references and report them in Table S3 of the Supporting Information.

The requirement of bigger impurity clusters implies the storage and manipulation of a higher number of electron repulsion integrals (ERIs). With our DF implementation, we observe a severe reduction in memory requirements. For the test case of the CO + Mg5O4 fragment at the equilibrium geometry with more than 200 orbitals, the 4c–2e calculation requires 200 Gb of memory on a AMD EPYC 7502 32-core processor, while the DF integral calculation requires 30 Gb of memory. This is because the previous implementation in refs (56 and 62) required storage of two-electron integrals, whereas the new implementation requires storage of decomposed intermediates, where Nimp is the number of impurity fragment orbitals and Naux is the number of auxiliary density-fitting orbitals. For large impurity fragments, and the storage saving becomes significant. All of the calculations use multithreading (i.e., shared memory) only; there is no multiprocessing. It should be noted that Cui et al.[57] also mentioned the possible memory efficiency that can be achieved by incorporating a similar algorithm. The current algorithm does not involve a formal speedup in terms of time required for a particular calculation. In this particular example, the 4c–2e calculation requires a wall time of 2 h and 21 mins, whereas the DF integral calculation requires a wall time of 1 h and 55 min.

Now, we examine the scaling of CCSD-DMET compared to the reference CCSD calculations. The computational cost of DMET is mainly dominated by the cost of the CCSD calculations within the embedding space. The most expensive term required in CCSD has a scaling of , where Nvir is the number of virtual orbitals and Nocc is the number of occupied orbitals.[42] The primary contribution to the large scaling arises from the virtual orbitals. In the current framework, most of the virtual orbitals are part of the environment (i.e., a part of Ncore orbitals), thereby significantly reducing the cost. If all Ncore unentangled orbitals are doubly occupied, the scaling of CCSD-DMET becomes ; if they are all empty, the scaling of CCSD-DMET becomes .

In summary, we have used a periodic implementation of DMET to calculate the adsorption energy of the CO molecule with the MgO(001) surface. We have investigated two widely used quantum chemical solvers, CCSD and MP2, as high-level methods. We infer that DMET-CCSD and DMET-MP2 can be used to obtain adsorption energies with high accuracy and at a significantly lower cost compared to the non-embedding references. We additionally observed that an impurity cluster including at least a MgO4 moiety on the MgO surface is required for accurate adsorption energies. Therefore, we implemented an efficient way to store and manipulate the memory-intensive ERIs within the periodic DMET algorithm. We envision that, with our recent implementation of the multireference solvers,[62] such as complete active space self-consistent field (CASSCF)[74−76] and n-electron valence state second-order perturbation theory (NEVPT2),[77−80] this approach can allow us to study bond-breaking phenomena of multireference systems on surfaces, at an affordable cost, which would be otherwise non-trivial for mean-field and single reference methods.