Surface science using coupled cluster theory via local Wannier functions and in-RPA-embedding: The case of water on graphitic carbon nitride

First-principles modeling plays a pivotal role in the understanding of surface phenomena at the atomistic level. Adsorption processes, for instance, can only be understood with a quantum mechanical theory based on the many-electron Schrödinger equation. Density functional theory (DFT) calculations in the Kohn–Sham framework of approximate exchange and correlation (xc) energy functionals are the most popular approach for many theoretical surface science studies due to their good trade-off between accuracy and computational cost. However, the presence of uncontrolled approximations in the employed xc functionals often makes it difficult to achieve the desired level of accuracy. Alternative methods, such as many-electron perturbation theories, achieve in many systems a more controllable accuracy and could complement DFT calculations favorably as long as the number of particles does not become computationally intractable. However, despite the fact that molecular adsorption can often be understood as a relatively local process, it is challenging to take full advantage of locality and realize a computational model that requires a quantum mechanical description of only a small number of particles in the vicinity of the adsorption site without compromising accuracy. In this work, we apply a recently presented quantum mechanical embedding approach based on a plane-wave basis set implementation to a prototypical molecular adsorption problem. We obtain highly accurate adsorption energies on the level of quantum chemical many-electron theories up to coupled cluster with single, double, and perturbative triple particle–hole excitation operators [CCSD(T)].

Here, we investigate the case of the adsorption of a single water molecule on graphitic carbon nitride (gC₃N₄). Graphitic carbon nitride is a promising compound for catalysis and energy storage.^1–4 It exhibits a high thermal stability and is easily processable from abundant materials. In addition, its bandgap is sufficiently large to overcome the endothermic character of the water-splitting reaction pathway,⁵ which makes it interesting for applications in photocatalysis. This has also motivated theoretical studies of molecular adsorption processes.^6,7 An initial step in such a reaction pathway is the adsorption of a single water molecule on the surface. Theoretical studies can provide estimates of the adsorption energy but have so far been limited to calculations on the level of DFT. This can partly be attributed to the fact that the two dimensional surface is corrugated and needs to be modeled using relatively big supercells, making it difficult to apply computationally more expensive but accurate methods, which are typically limited to unit cells containing only a few tens of atoms. Therefore, it can be assumed that embedding approaches, as described in the present work, help to remove limitations imposed by the simulation cell size.

The theoretical study of the adsorption process of a single water molecule on a two dimensional material unveils strengths and weaknesses of many widely used electronic structure theories. The importance of weak van der Waals (vdW) interactions in such systems cannot be overestimated. Unfortunately, the ambiguous choice of the xc functional in DFT calculations and the availability of numerous approximations that aim at accounting for the effect of vdW interactions can make it difficult to achieve a faithful ab initio prediction of adsorption energies. Benchmark results can provide helpful insight into which functional is the most suitable for a particular system. Although high-level methods such as CCSD(T) or quantum Monte Carlo (QMC) theories have successfully produced accurate benchmark results for periodic systems, the studied surface models often have been limited to a few tens of atoms only.^8–10 Here, we advance our recently published embedding strategy,¹¹ which enables a systematic path to approach the accuracy of CCSD(T) for regional phenomena in large systems of more than thousand electrons or periodic supercells with edge lengths of more than 20 Å.

Regional embedding techniques for first-principles methods, in general, form an active research field.^12–21 This strategy was likewise adapted for many-electron correlation methods applied to molecules as well as periodic systems.^11,22–36 However, regional phenomena in extended (periodic) environments exhibit long-range correlation effects, such as lattice relaxations or electron dispersion, which often exceed the size of a local fragment. As we show in this paper, the bonding of the water molecule to gC₃N₄ cannot exclusively be explained by H–N hydrogen bonds, but also long-range van der Waals interactions play a considerable role. Additionally, ionic relaxations alter the electron density far away from the adsorption site. Hence, restricting electron correlation to local fragments only introduces sizable errors. A proper treatment of electron correlation between the local fragment and the environment is therefore indispensable in surface science and other regional phenomena, such as defect or polaron formation in solids.^11,35 To achieve this, our plane-wave based scheme seamlessly embeds a high-level into a low-level correlation method. A convergence of the adsorption energy with respect to the local fragment size corresponds to a systematic approach to the high-level accuracy for the full system. Rapid convergence can be achieved if the low-level and high-level dispersions are quantitatively similar, as is observed for embedding coupled cluster theory (CC) into the direct random phase approximation (RPA).¹¹ The choice of the RPA as a low-level correlation method is particularly favorable since it can be evaluated efficiently with cubic scaling algorithms using a plane-wave basis.³⁷ Furthermore, we note that many implementations of the RPA in the framework of different basis sets and boundary conditions are available and have been successfully applied to a wide range of materials and properties.^38–43 

The paper is organized as follows: In Sec. II, we outline our scheme to embed a high-level into a low-level correlation method in detail, including the algorithm to construct localized Wannier orbitals as well as technical details to the many-electron correlation methods. The performance and results of our embedding scheme applied to the adsorption energy of H₂O on gC₃N₄ are presented and discussed in Secs. III A and III B, respectively.

We model the adsorption by means of a 3 × 3 single layer of graphitic carbon nitride (gC₃N₄) in a supercell with periodic boundary conditions. We consider the free energy (equivalently, the ground state energy) difference

where three equally shaped supercells are employed, containing both the adsorbed water on the layer (H₂O@gC₃N₄), the layer alone (gC₃N₄), and the molecule alone (H₂O). The supercell for H₂O@gC₃N₄ contains 584 valence electrons from 54 carbon atoms, 72 nitrogen atoms, and the water molecule using a vacuum of 20 Å, as depicted in Fig. 1. This corresponds to 292 spin-restricted occupied Bloch bands. The equilibrium positions of the atoms are obtained from an ionic relaxation using the DFT-PBE functional for all three unit cells and can be found in the supplementary material. The computationally most expensive steps of these free energy calculations arise from the many-electron correlation calculations, which are treated in an efficient embedding framework. We explain this embedding procedure in detail in Sec. II A.

All calculations start from a periodic spin-restricted Hartree–Fock (HF) reference of the entire supercell at the Γ-point. We use the plane-wave based Vienna Ab initio Simulation Package (VASP)⁴⁴ to calculate the HF ground state, natural orbitals,^45,46 as well as the many-electron correlation at the second-order Møller–Plesset perturbation theory (MP2)^47,48 and RPA³⁷ level of theory. The basis set is defined by the single cutoff parameter $ENCUT=500eV$ and leads to about 183 000 plane-wave basis vectors, which correspond to more than 600 basis vectors per occupied band. VASP makes use of the projector augmented wave (PAW) method⁴⁹ in combination with frozen core potentials; see Table I for details. The coupled cluster with singles and doubles (CCSD) and CCSD(T) calculations are performed with the coupled cluster for solids (cc4s) code.^10,50 All levels of many-electron correlation calculations are based on the exact same set of two-electron integrals, avoiding implementation-specific numerical issues. The calculations are performed on the Vienna Scientific Cluster equipped with 48 Intel Platinum 8174 processors and 384 GB per node.

Our embedding scheme allows for the controllable regional embedding of a high-level into a low-level correlation method for calculations in large supercells. To this end, we partition the electron density of the entire supercell into two fragments, a local region of interest (F) and the rest of the cell (R). In the considered adsorption process, the fragment F consists of the water molecule (in the case of the H₂O@gC₃N₄ supercell) and nearby parts of the densities of the gC₃N₄ layer. The construction of the fragments is achieved by selecting localized Wannier orbitals. In the following, we restate the individual technical steps of our embedding scheme as presented in our previous work.¹¹ The main difference consists in a new implementation of an algorithm for the construction of localized Wannier functions of various types for periodic systems in the plane-wave basis, which we discuss in detail in Sec. II B.

Once the density of the fragment F and the rest R are defined by two disjoint sets of localized Wannier functions, the correlation energy calculation using method M can be restricted to orbitals in F or R (see the six steps in the following), denoted as $EFM$ or $ERM$⁠, respectively. We can then decompose the correlation energy of the entire cell as

This decomposition is exact since $EF↔RM$ is implicitly defined by the above-mentioned equation, while all other quantities are explicitly defined. The contribution $EF↔RM$ can be interpreted as the inter-fragment correlation. Embedding the high-level method M₁ into the low-level method M₂ (in short M₁:M₂) can then intuitively be defined by

Provided the low-level method M₂ can be applied to all electrons in the entire supercell, the presented embedding scheme accounts for all many-electron correlation effects but at mixed levels of theory.

The actual restriction of a correlation method to a local fragment of the density is by no means trivial nor unique. The following six steps explain our approach, which involves an efficient compression of the unoccupied space and allows for a systematic convergence of the correlation energy with respect to the size of the fragment F.

• Hartree–Fock ground state

  Self-consistent calculation of the N_occ occupied HF Bloch orbitals.

• (ii)

  Unoccupied space

  Diagonalize the HF Hamiltonian in the full plane-wave basis to obtain all N_vir unoccupied Bloch orbitals.

• (iii)

  First compression of unoccupied space

  Calculate the natural orbitals at the level of the RPA⁴⁶ for the entire supercell. Rediagonalize the Fock matrix using the first $Nvirpre=60⋅Nocc$ natural orbitals (NOs) in order to precompress the unoccupied space by a factor of ten as compared to the full plane-wave basis.

• (iv)

  Build fragments

  Transform all occupied Bloch orbitals to N_occ Wannier orbitals. Here, we use the intrinsic bonding orbitals (IBOs), which are a special case of the Pipek–Mezey (PM) orbitals; see Sec. II B. To each atom, we assign a set of orbitals based on the partial charge of the Wannier orbitals. The fragment is then defined by selecting atoms in the gC₃N₄ layer, based on the radial distance from the adsorbed water molecule, which is included in each fragment. Finally, we diagonalize the HF Hamiltonian in the basis of the $Noccloc$ selected Wannier orbitals (re-canonicalization). This simple selection scheme is sufficient as long as the local process—as in our case—is not accompanied by major changes in geometry. In other cases, adapted selection schemes, such as the even-handed selection,⁵¹ could be considered.

• (v)

  Second compression of unoccupied space

  Calculate the natural orbitals of the fragment only, based on an approximated MP2 scheme⁴⁵ using $Nvirpre$ pre-compressed orbitals to further compress the unoccupied space to $Nvirloc=Y⋅Noccloc$ orbitals. For MP2 and RPA, we choose Y = 60, and for CCSD and CCSD(T), we choose Y = 20 augmented by a basis set correction procedure.

• (vi)

  Local correlation energy

  Decompose the electron repulsion integrals (ERI) into auxiliary three index quantities using N_F auxiliary field variables, i.e., auxiliary basis functions^11,50 Calculate the correlation energy (MP2, RPA, CCSD, CCSD(T), …) for $Noccloc$ local occupied and $Nvirloc$ unoccupied orbitals and employing N_F auxiliary field variables.

Please note that steps (i) and (ii) are common preparatory steps for regular correlation energy calculations. The embedding procedure mainly consists in steps (i) + (ii) + (iv) + (vi). The compression of the unoccupied space in steps (iii) and (v) primarily serves to increase the computational efficiency. Further technical details to the correlation methods are outlined in Sec. II C.

Here, we present the theory and implementation of the algorithm to construct localized Wannier functions of various types for large supercells in the Vienna Ab Initio Simulation Package (VASP). The algorithm essentially follows the work of Abrudan et al.,⁵² which has already been used for orbital localization of molecules and solids in Refs. 54–55. Localized Wannier orbitals $Wi$ are obtained by a unitary transformation of the (delocalized) spin-orbitals $χj$⁠,

The expected Brillouin zone integration⁵⁶ is neglected here, since we currently aim for large supercells using a Γ-point-only sampling. The unitary matrix u_ij is determined by the optimization (maximization/minimization) of a cost functional $L$⁠, which ultimately defines the localized Wannier orbitals. For instance, the PM localization⁵⁷ is defined by maximizing the atomic partial charges,

where p ≥ 2 and $PA=∑μ∈Aμμ$ is a projector onto atom centered functions $μ$ at the atom A. In the case of spin-polarization, a separate cost functional is considered for each spin. While our implementation allows for the optimization of any cost functional $L$⁠, which is a polynomial in u_ij (including also other optimizations beside orbital localization), we choose the PM localization cost functional and its variants as a specific example. In particular, we consider the IBOs, which are obtained by choosing intrinsic atomic orbitals (IAOs) $μIAO$ for the projector P_A in the construction of the PM cost functional.⁵⁸ 

The algorithm to maximize/minimize the considered cost functional is essentially an iterative unitary matrix constrained optimization (UMCO) procedure. Starting from a stochastic unitary matrix U = (u_ij), each iteration step propagates the unitary matrix according to

The Riemannian derivative G = ΓU^† − UΓ^† is anti-Hermitian, G^† = −G, and is calculated from the Euclidean derivative $Γ=∂L/∂U*$⁠. A conjugate-gradient technique based on the Polak–Ribière–Polyak update factor is applied in order to speed-up the convergence. We refer the reader to Ref. 52 for details. The iterations are considered as converged, once the norm obeys ⟨G, G⟩ = tr(GG^†)/2 < δ where we choose δ = 10⁻¹⁰. For the case of the PM functional, we find

The optimal step size μ obeys $∂∂μL[exp(±μG)U]=0$ and is estimated by a polynomial line search approximation as described in Ref. 52. A positive sign +μ is used for the maximization (e.g., for the PM functional), while −μ corresponds to a minimization of $L$⁠.

The plane-wave basis and the specifics of the PAW method have to be respected only in the preparatory calculation of the overlaps (scalar product) between ⟨χ_i|μ⟩ to compute the action of the atomic projectors P_A. The calculation of these overlaps is parallelized and has previously been implemented for an interface between VASP and wannier90 in Ref. 59. In contrast, the iterations of Eq. (6) are not parallelized, instead run individually on each rank but with separate stochastic initialization. This redundancy allows for a simple consistency check of the reached extremum of the cost functional in order to reveal “local minimum” issues, if present. The computational complexity of the iterations is dominated by matrix multiplications; hence, follow a cubic scaling in the number of orbitals to localize. Timings are presented in Sec. III. The cubic scaling of the evaluation of $ΓijPM$ from Eq. (7) is not obvious and can be achieved using the following equivalent formulation:

The atomic overlap matrix $(XA)iμ=⟨χi|μ⟩$ with μ centered on the atom A only and its unitary rotated variant $YA=XA†U$ were used. For each atom, matrix multiplications have to be performed, resulting in a complexity of $Natoms⋅Nocc2(Nμ/Natoms)=Nocc2Nμ$⁠, where N_μ/N_atoms is the number of atom centered functions per atom.

As atom centered functions, we choose IAOs, |μ⟩ → |μ^IAO⟩, which can be obtained from atomic functions $fμ$ by a simple projection procedure,

where we follow the formulation of Janowski.⁶⁰ The atomic functions $fμ$ may be HF or DFT solutions of an isolated atom, but also hydrogen-type or Slater-type functions are a valid choice. The two projectors are defined by

While $O$ is simply the projection onto the occupied space, $Õ$ projects onto the space spanned by

with the overlap matrix S_μν = ⟨f_μ|f_ν⟩. The subsequent orthogonalization is denoted as “orth.” The orbitals decorated with a tilde, $χ̃i$⁠, can be understood as the occupied orbitals as described by the certainly incomplete “minimal basis” spanned by ${fμ}$⁠. Note that the occupied orbitals $χi$ originate from a preceding mean-field calculation using the plane-wave basis, which we consider as the complete basis. The action of $O−Õ$ on $fμ$ in Eq. (9) constructs the necessary augmentation such that the IAOs $μIAO$ form an exact basis of the occupied space as described by the complete basis. This is true as long as no single occupied orbital is orthogonal to the atomic functions, ∑_ν⟨f_ν|χ_i⟩ ≠ 0, ∀i. Additionally, if the density spillage of the orbitals with tilde is not too large, $μIAO$ exhibits almost the same spatial shape as the atomic function $fμ$⁠. In this case, the IAOs can be considered as a minimal but exact atom centered and local basis for the occupied space. Furthermore, the atomic partial charge in the PM cost functional of Eq. (5) is stable against different choices of atomic functions $fμ$ if the IAOs are used for the projectors P_A, as shown in Ref. 58. The resulting Wannier orbitals are thus the IBOs. The exponential decay of the gradient norm ⟨G, G⟩ using the presented UMCO algorithm applied to water on gC₃N₄ is shown in Fig. 2.

For the post-HF calculations of the localized Bloch orbitals (local Wannier orbitals) in the fragment, we employ MP2, RPA, CCSD, and CCSD(T) theory. The localized orbitals are re-canonicalized using VASP, making it possible to employ canonical formulations of all many-electron theories.

The virtual orbitals are computed using the following two-step procedure, as already outlined in Ref. 11. We first compute the RPA natural orbitals (RPANOs) for the full supercell prior to the localization procedure of the occupied orbitals.⁴⁶ We truncate the RPANOs to a subset constructed from 60 natural orbitals per occupied spatial orbital and re-canonicalize the Fock matrix in this subspace. We refer to the corresponding total number of unoccupied orbitals at this stage by $Nvirpre$⁠. In the second compression step, we compute approximate MP2 natural orbitals for a set of localized occupied orbitals. To this end, we calculate the virtual–virtual (⁠$Nvirpre×Nvirpre$⁠) block of an approximate MP2 reduced density matrix defined by Eq. (2) in Ref. 45. In this step, the sum over the occupied orbitals for the expression of the reduced density matrix is restricted to a semi-canonicalized set of occupied orbitals that are localized on the fragment. As a consequence, the resulting approximate MP2NOs with the largest occupation numbers are also localized on the fragment. We refer to the corresponding number of employed NOs by $Nvirloc$⁠.

Once an optimized set of localized occupied, and virtual orbitals has been obtained for each fragment, the required information for the post-HF calculations is passed to the cc4s code, which performs MP2, CCSD, and CCSD(T) calculations. In addition to the HF orbital energies, the decomposed ERI are required for the MP2 and CC calculations. We note that certain integrals are computed on-the-fly in order to minimize the memory footprint of the CCSD calculations. For the decomposed ERI, we employ a procedure that was previously described in Ref. 50 and involves auxiliary three index quantities using N_F optimized auxiliary basis functions. All Coulomb integrals, $Vsrpq$⁠, needed by coupled cluster theories can be computed using the following expression:

where p, q, r, and s refer to (occupied or virtual) orbital indices. F refers to auxiliary basis functions that have been obtained using a singular value decomposition outlined in Ref. 50. For the present work, due to the local character of both occupied and virtual orbitals, significant reductions of the auxiliary basis set size are possible without compromising the precision of computed correlation energies. We have carefully checked that the computed correlation energies are converged to within meV with respect to the size of the optimized auxiliary basis set. We stress that as a result of the optimization procedure described above, all dimensions that contribute to the scaling of the computational complexity of the electron correlation methods exhibit a linear scaling with fragment size: $NF∝Nvirloc∝Noccloc$⁠, despite using a plane-wave basis set.

To achieve basis set converged electronic correlation energy differences on all levels of theory, we employ the following procedures. The MP2 correlation energy calculations are performed for each fragment using 60 natural orbitals per occupied orbital. In addition, an automatic basis set extrapolation to the complete basis set (CBS) limit is performed. The extrapolation procedure was previously described in Ref. 47. Here, the extrapolation is performed for each electron pair correlation energy contribution to yield $ϵijMP2$⁠. Throughout this work, we refer to CBS limit estimates when referring to the MP2 energy.

The CCSD correlation energies are computed using 20 natural orbitals per occupied orbital for each fragment. To correct for the remaining basis set incompleteness error (BSIE) of CCSD correlation energies, we employ a correction that was previously introduced in Ref. 61. This correction includes contributions from two terms: (i) a BSIE correction on the level of MP2 theory and (ii) an approximation to the BSIE of the particle–particle-ladder (PPL) term contribution to the CCSD energy. We have concluded in Refs. 62–63 that these are indeed the dominant BSIE contributions on the level of CCSD theory for a wide range of properties. In passing, we note that Refs. 64 and 65 have shown related findings. The employed PPL correction is electron pair specific and depends on $ϵijMP2$⁠. As described in Ref. 61, this approach (also referred to as FPb+Δps-ppl) yields CCSD reaction energies that deviate from the CBS limit by about 11 meV only, which is closer to the CBS limit than F12 calculations employing an aug-cc-pVTZ basis set. For the sake of brevity, we will simply refer to the CCSD CBS limit energy estimates obtained in this manner as CCSD energies.

The perturbative triples (T) correlation energy contribution is computed using 20 natural orbitals per occupied orbital, which is considered a well converged CBS limit estimate. This is partly confirmed by benchmark results published in Ref. 61.

We now turn to the results obtained using the embedding approach and their convergence with respect to the number of occupied orbitals in the studied fragments. Figure 3 shows the first seven fragments built from selected IBOs as used for the embedding scheme. The localized occupied orbitals included in the employed fragments are selected according to their distance from the water molecule. The full set of localized IBOs exhibits an average orbital spread (standard deviation to the orbital center) of 0.937 and 0.939 Å for the layer with and without the adsorbed water molecule, respectively. The IBOs with the largest spread read 1.324 and 1.323 Å, respectively. Apart from the relaxed atomic positions and formed H–N hydrogen bonds, we obtain qualitatively similar sets of IBOs in the layer for both simulation cells.

We first discuss the convergence of direct MP2 (dMP2, i.e., MP2 without exchange-like term) and RPA adsorption energies with respect to the fragment size. Figure 4 depicts water adsorption energies retrieved as a function of the number of occupied orbitals in the fragment. The horizontal dashed lines show the dMP2 and RPA reference adsorption energies obtained from calculations of the full periodic cell. We note that the adsorption energies obtained for the fragments include electron correlation contributions of the bare fragments only, i.e., only $EFM1$ from Eq. (3) without embedding into a low-level theory M₂. As apparent from Fig. 4, no convergence with respect to the fragment size can be achieved for dMP2 and RPA in the present system for the considered fragments. This lack of convergence for the non-embedded fragments was already observed in Ref. 11 for the case of the ionically relaxed surface of TiO₂ adsorbing a water molecule. In contrast, the embedding into a low-level method, M₁:M₂, leads in all cases to a clearly smoother convergence. The MP2:dMP2 embedding rapidly approaches its limit at −0.54 eV, as shown in Fig. 5. The progression to fragments with up to 120 correlated bands reveals a remaining numerical noise of about 0.02 eV. We therefore estimate the error of our embedding results to this value.

Having observed that the embedding into a low-level method substantially improves the convergence, we now seek to better understand the source of the slow convergence of the non-embedded correlation energies in Fig. 4. Due to relaxations of the atoms caused by the water adsorption, the electron density is considerably altered in the entire supercell, as can be observed at the right bottom in Fig. 3. Note that the dipole of the water molecule alone would only alter the density in a much smaller region of the layer. The density change in the entire cell, however, diminishes the similarity of the fragments with and without the adsorbed water molecule, albeit the fragments were built using the same selection of atomic sites in the layer (step iv) of the corresponding supercells. The diminished similarity manifests itself in the fragment’s surfaces. As apparent from the lower right plot in Fig. 3, this cannot be eased using larger and larger fragments. Instead, the surface effect is only pushed along as the fragment size increases but is abruptly corrected as soon as the fragments reach the size of the supercell, i.e., when the surfaces melt into their periodic images. Conversely, the convergence of the non-embedded fragments is restored if the atomic sites are not relaxed, i.e., if the positions of the atoms in the layer are equivalent in the supercell with and without the adsorbed water. We show the analog to Fig. 4 for the non-relaxed case in the supplementary material. There, we also present the individual contributions $EFM$⁠, $ERM$⁠, and $EF↔RM$ from Eq. (2) for M = RPA. From this, we conclude that the change in atomic positions in our studied material prohibits a local fragment-only approach to capture the electron correlation effects of the considered adsorption process. Instead, the M1:M2 embedding is essential in order to achieve convergence with respect to the fragment size at all.

In the following, we discuss the performance of our embedding scheme for the coupled cluster methods, as shown in Fig. 6. As mentioned above, we believe that embedding a high-level into a low-level correlation method allows us to remedy two central issues of local fragments in periodic environments: (i) it accounts for the missing non-local long-range correlation and (ii) it cancels spurious long-range relaxation effects in the atomic positions. The latter is particularly important if lattice relaxations reduce the similarity between the considered supercells. Slightly displaced ions result in a change in the electron density, leading to a variation of the localized Wannier orbitals and accordingly to a different character of the fragment surface. Nevertheless, even in situations without relaxation of the atomic positions, non-local long-range dispersion interactions can prohibit a fast and smooth convergence of the correlation energy with respect to the fragment size. In Ref. 11, we demonstrated that a rapid convergence can be observed if the low-level and high-level long-range dispersions are quantitatively similar. In the language of diagrammatic perturbation theory, we expect this if the low-level method is a diagrammatic subset of the high-level theory, where the former can easily be restricted to capturing only long-range terms accurately. For the MP2:dMP2 and CCSD:RPA embedding, this is confirmed by our presented results. In the first case, the low-level method accounts for the full long-range correlation, while the second case benefits from the fact that the direct ring diagrams capture the important part of the screened dispersion correlation of CCSD. However, we have to assume equally that the long-range behavior of the perturbative triples in CCSD(T) cannot be well approximated by neither dMP2 nor RPA. The rapid convergence of CCSD(T):RPA and CCSD(T):dMP2 is therefore primarily due to the rapid decay of the triples themselves. Presumably, the seemingly faster convergence of CCSD(T):dMP2 can be attributed to the fortuitous similarity of the final energies of dMP2 and CCSD(T). Again, we believe that CCSD(T):RPA is generally more favorable since the underlying CCSD:RPA is guaranteed to converge quickly while the perturbative triples are hardly corrected.

We now discuss the converged adsorption energies obtained from each method considered. The results of each method are compiled in Table II. The HF result of −206 meV provides a reference to measure the contribution of the electron correlation of the water adsorption since the HF captures the full electrostatics but lacks electron correlation by definition. The PBE result indicates significant correlation effects when comparing to the uncorrelated HF result. As corrected functionals, we employed PBE-D3 (the method of Grimme et al.⁶⁶ with Becke–Jonson damping), PBE-TS (the Tkatchenko–Scheffler method⁶⁷), optB86B (the method of Langreth and Lundqvist in the optimized form as presented in Ref. 68), and PBE0-MBD (the many-body dispersion energy method as presented in Ref. 69 using the scaling parameter β = 0.85). These corrections for long-range correlation suggest an even stronger impact such that the total electron correlation contribution to the adsorption energy can be assumed to be of a similar, if not larger, order of magnitude as the electrostatic energy. This suggestion is supported by the many-electron correlation methods dMP2 and RPA, which could be applied to the entire supercells. However, the large difference between dMP2 and RPA of about 100 meV prohibits a consistent picture. In passing, we note that we perform RPA calculations using a HF reference, which is in contrast to the more prevalent approach of performing RPA calculations on top of a Kohn–Sham reference. The discrepancy between dMP2 and RPA comes as no surprise since the accuracy of both methods is diminished in extended systems by the following systematic problems. To mention only the most eminent issues, the MP2 (as well as dMP2) method suffers from severe overestimations of the long-range correlation⁷⁰ due to the missing electron screening, which, on in contrast, is considerably more accurately captured by the RPA, as observed for many materials.^{38,40–43,72–75} Yet the RPA results may be deteriorated from the self-correlation problem of short-range correlation and from too deep a correlation hole.^77–78 The estimated CCSD(T) result finally confirms the well-known overcorrelation and undercorrelation trends of MP2 and RPA, respectively. With about 0.05 eV, the triples contribute 16% to the correlation contribution and 10% to the full adsorption energy. Having computed a reliable CCSD(T) adsorption energy, we can conclude that PBE-D3 and PBE0-MBD achieve an excellent trade-off between accuracy and cost for the present system. We note, however, that this is in contrast to a related study of water adsorption on a single layer of hexagonal boron nitride, where PBE-D3 and PBE0-MBD overestimated the adsorption energy compared to more accurate CCSD(T) and QMC findings.⁹ 

We also report the timings of the calculations in the steps (i)–(vi) from Sec. II A. All timings correspond to the largest system, i.e., H₂O@gC₃N₄. The HF iterations (step i) converged in about 210 core hours. The calculation of the unoccupied bands (step ii) takes about 3300 core hours. Natural orbitals at the RPA level (step iii) cost about 21 000 core h. The localization to construct the IBOs (step iv) finished in 0.2 core h after 87 iteration steps. Timings for the calculation of the second compression (natural orbitals at the level of MP2, step v) as well as for step (vi), MP2 and CCSD(T), can be found in Fig. 7 for different fragment sizes.

Apparently, the construction of the localized Wannier functions is the least expensive step.

We have calculated the adsorption energy of water on graphitic carbon nitride employing coupled cluster theory including singles, doubles, and perturbative triples, CCSD(T). Results at the level of MP2, RPA, and CCSD were also presented. This was made possible by an in-RPA-embedding that leverages the numerical similarity between CC and RPA for the electronic long-range correlation. Recently developed techniques were employed to reach the complete basis set limit for all correlation methods.

Furthermore, we presented the implementation of an algorithm that enables various orbital localizations for periodic systems, provided that a proper cost functional can be defined. The construction of localized Wannier functions using the Pipek–Mezey cost functional (in the variant for intrinsic bonding orbitals) was demonstrated and used to define the local fragments for the embedding procedure.

Our findings show that fragment based techniques for black box correlation calculations of regional phenomena invariably need to correct for effects far outside the region of interest. These effects may be due to long-range dispersion correlation and distorted local fragments caused by long-range density alterations, which, in turn, are triggered by far reaching ionic relaxations of the regional process. We demonstrated that such effects are reliably captured by the RPA.

See the supplementary material for the employed atomic structures (POSCAR files) and adsorption energies without preceding ionic relaxations.

The authors acknowledge support and funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. 715594). Andreas Grüneis and Andreas Irmler acknowledge support from the European Union’s Horizon 2020 research and innovation program under Grant Agreement No. 951786 (The NOMAD CoE). T.S. extends his special thanks to Henrique Miranda and Manual Engel for the enhanced parallelization efficiency to calculate the scalar products (overlaps) between the Bloch orbitals and local trial functions for the new UMCO module. We thank Livia Getzner for providing the relaxed atomic structure used in the present study. The computational results presented have been achieved, in part, using the Vienna Scientific Cluster (VSC). The authors acknowledge TU Wien Bibliothek for financial support through its Open Access Funding Programme.

The authors declare no conflicts of interest.

The data that support the findings of this study are available within the article and its supplementary material.