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.

M agnesium 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 Dohnaĺek 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 functionin-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 N imp ) 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 (D env ) of the one-body reduced density matrix (1-RDM) is diagonalized as follows: where λ is a diagonal matrix of eigenvalues λ i (i = 0, 1, ..., N env , where N env is the number of the environment orbitals). The columns of the unitary matrix U corresponding to λ i 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: (2) where |f i , |b i , and |core are single determinants in the Fock spaces of the N imp impurity orbitals, N bath ≤ N imp bath orbitals, and N core frozen core orbitals, respectively, and i = 0, 1, ..., N imp . The high-level wave function, | emb , diagonalizes the embedding Hamiltonian, H emb where H emb is the partial trace of the Hamiltonian over the |core determinant; its operator terms involve only N emb = N imp + N bath ≤ 2N imp embedding orbitals. For calculations of energy differences, it is important to choose the same number of N emb 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 N core 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 as where P and Q represent auxiliary basis functions, M PQ = (P | Q) is the Coulomb metric, and B ij P and B kl P are the Cholesky where the expansion coefficients C ij R are obtained by solving a linear equation The 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 4. 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 O N ( ) In other words, the right-hand side of eq 4 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 O N N ( ) aux imp 2 , 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, E eq (C−Mg bond distance of 2.479 Å), 8 and at a separated geometry, E sep (C−Mg bond distance of 6 Å), as indicated in Figure 1a. 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 1b. 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 + MgO 4 ), and finally (d) the 4 next to nearest Mg atoms in addition to choice (c) (denoted as CO + Mg 5 O 4 ). The orbitals localized at the highlighted atoms in the embedding clusters (Figure 1b) 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.
In Figure 2, we report the relative energy E rel 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 E rel 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 + MgO 4 impurity subspace.
For the MP2 reference method, E rel reaches an asymptotic value at a C−Mg distance of 6 Å and differs from E rel at 5 Å by only 0.3 kcal/mol, thereby suggesting that 6 Å is a reasonable choice for a separated geometry. Using the CO + MgO 4 fragment, the DMET E rel 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 E rel values deviate significantly from the MP2 references. Morover, the RHF E rel 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 3. The numbers are reported in Tables S2 and S3 of    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 + Mg 5 O 4 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  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 MgO 4 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 selfconsistent field (CASSCF) 74−76 and n-electron valence state 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.
The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter 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.
■ ASSOCIATED CONTENT
Mean absolute deviations (MADs), adsorption energies, relative energies for C−Mg distances, memory savings for DF integrals, sample input and output files, and total energies (PDF) Transparent Peer Review report available (PDF)