Applying the Coupled-Cluster Ansatz to Solids and Surfaces in the Thermodynamic Limit

Modern electronic structure theories can predict and simulate a wealth of phenomena in surface science and solid-state physics. In order to allow for a direct comparison with experiment, such ab initio predictions have to be made in the thermodynamic limit, substantially increasing the computational cost of many-electron wave-function theories. Here, we present a method that achieves thermodynamic limit results for solids and surfaces using the"gold standard"coupled cluster ansatz of quantum chemistry with unprecedented efficiency. We study the energy difference between carbon diamond and graphite crystals, adsorption energies of water on h-BN, as well as the cohesive energy of the Ne solid, demonstrating the increased efficiency and accuracy of coupled cluster theory for solids and surfaces.


INTRODUCTION
Modern ab-initio methods to solve the electronic Schrödinger equation for real solids and molecules such as density functional theory or wave function based methods are becoming increasingly accurate and efficient [1][2][3][4]. However, in contrast to molecular systems, properties of solids and surfaces need to be calculated in the thermodynamic limit. The convergence towards the thermodynamic limit with respect to the number of particles is very slow, often exceeding the computational resources of even modern super computers. This is particularly the case for many-electron wave function based theories that allow for a systematic improvability upon the description of the electronic correlation energy. Nonetheless these methods are becoming increasingly popular in theoretical physics as well as chemistry to treat electronic correlation in periodic condensed matter systems with high accuracy [3][4][5][6][7][8][9][10][11][12][13].
Electronic correlation is for the most part a shortranged phenomenon. The proper description of the wave function shape at short interelectronic distances allows for capturing the largest fraction of the correlation energy in solids [2,14]. Significant progress has been achieved for many-electron wave function based theories by exploiting the locality of electronic correlation in large molecules and solids. The development of so-called local correlation methods and embedding theories has improved their computational efficiency considerably [15][16][17][18][19][20][21][22][23][24]. However, theories that approximate long-range correlation effects * andreas.grueneis@tuwien.ac.at such as van der Waals interactions must carefully be checked for convergence with respect to the employed cutoff parameters to allow for accurate and predictive ab-initio studies of real materials. This is of particular importance in condensed matter systems where the accumulation of weak van der Waals interactions can become a non-negligible contribution to the property of interest as for example in the case of the energy difference between carbon diamond and graphite or the adsorption of a water molecule on an h-BN sheet. Pairwise additive interatomic van der Waals interactions cause an 1/N convergence of the electronic correlation energy per unit cell in insulating three dimensional systems, where N is the number of explicitly correlated atoms. However, in general the exact form of these scaling laws depends on the dimensionality and electronic response properties of the system, e.g. the adsorption energy of molecules on two dimensional insulating surfaces exhibits an 1/N 2 convergence. Moreover we note that collective phenomena such as plasmons in metallic systems can also modify the observed scaling laws [25]. For these reasons robust and reliable approximations to long-range correlation effects are non-trivial.
Many-body methods such as coupled cluster or configuration interaction theory can describe both long-and short-ranged electronic correlation effects with high accuracy. However, the scaling of the computational complexity of these theories with respect to system size is either of a high-order polynomial or even exponential form. Therefore it is difficult to treat long-range correlation effects in a computationally efficient manner using these theories. This has lead to the development of various techniques that partition the correlation problem arXiv:2004.06584v1 [cond-mat.mtrl-sci] 14 Apr 2020 according to a predefined criterion such as the distance between electron pairs or fragment size. Local theories employ correlation energy expressions that depend on localized electron pairs, making it possible to treat longdistance pairs using computationally more efficient yet less accurate theories. Embedding theories typically aim at combining the computational efficiency of mean field theories for the long-range with the high accuracy of wave function based methods applied to small fragments only. In this work we introduce an efficient method that seamlessly integrates long-range correlation effects for solids without any predefined criteria such as cutoff distance or fragment size. Our approach is inspired by structure factor interpolation techniques as performed in the field of Quantum Monte Carlo theory [26]. However, in coupled cluster theory the structure factor, being the functional derivative of the total energy with respect to the Coulomb kernel, is not directly available. Instead we seek to interpolate the partial functional derivative of the coupled cluster correlation energy expression with respect to the Coulomb kernel. The interpolation scheme is chosen such that it is directly transferable to systems with arbitrary dimensions including solids and surfaces. Due to the adverse scaling of the computational complexity in coupled cluster theories, the proposed method allows for reducing the computational cost by several orders of magnitude without compromising accuracy compared to previous studies [4].

THEORY
The electronic correlation energy can be calculated in a plane wave basis set using the following expression [27] In the above equation G corresponds to a plane wave vector that is defined as G = g + ∆k, where g is a reciprocal lattice vector and ∆k is the difference between any two Bloch wave vectors that are conventionally chosen to sample the first Brillouin zone. v(G) is the Coulomb kernel in reciprocal space that diverges at G = 0, making it numerically necessary to disregard this contribution to the sum as indicated by the apostrophe. Thus, S(G) is the partial functional derivative of the correlation energy with respect to v(G) and we will return to its explicit definition later as well as in the supplementary information [27].
The thermodynamic limit is approached as N → ∞, where N is the number of particles in the simulation cell while the density is kept constant. Finite size errors are defined as the difference between the thermodynamic limit and the finite simulation cell results. For electronic correlation energies obtained using many-electron perturbation theories these errors typically decay as 1/N as a consequence of long-range interatomic van der Waals forces. In the thermodynamic limit G of Eq. (1) is replaced by G . Therefore finite size errors in the correlation energy of periodic systems originate from two sources [28]: (i) quadrature errors in the summation over G, and (ii) the slow convergence of S(G) with respect to the employed supercell size or k-point mesh. In the following we will discuss how to reduce both errors substantially.
We first seek to discuss finite size errors originating from the quadrature in the summation over G. These contributions can be partitioned into the G = 0 volume element contribution and the remaining terms. We note that as a result of the Coulomb divergence, the integrable contribution of S(0)v(0) to the correlation energy is usually neglected in computer implementations of Eq. (1) [11,27]. However, this is the dominant contribution to the finite size error of the correlation energy of insulators. A Taylor expansion of S(G) around G = 0 shows that S(G) exhibits a quadratic behavior close to zero, explaining the 1/N decay of the finite size error for three dimensional insulators [27]. An estimate of S(0)v(0) can be obtained by spherically averaging S(G) and interpolating around G = 0. Subsequently, the interpolated function is multiplied with the analytic Coulomb kernel and integrated over a sphere around G = 0, yielding an estimate of S(0)v(0) [27]. However, this approach is not well defined for anisotropic systems because it requires a spherical cutoff parameter. In this work we propose to interpolate S(G) using a tricubic interpolation without spherical averaging. Once obtained the interpolation of S(G) and the analytic expression for the Coulomb kernel allows for integrating over G on a very fine grid, simulating the thermodynamic limit integration. This approach accounts for the S(0)v(0) contribution to the correlation energy and reduces quadrature errors originating from too coarse a Brillouin zone sampling. We will refer to coupled cluster correlation energies obtained using this interpolation strategy as CC-FS.
To illustrate the importance of the interpolation method we consider the following example. Figure 1 shows slices and an isosurface of the interpolated S(G) for carbon graphite in the ABC stacking. Black dots indicate sampling points of S(G) obtained using coupled cluster singles and doubles theory and a 4 × 4 × 4 kpoint mesh. This figure illustrates that S(G) is very anisotropic around G = 0. Furthermore we show that even a 4 × 4 × 4 k-point mesh sampling indicated by the black dots corresponds to a relatively coarse grid, causing non-negligible quadrature errors. We will return to the discussion of the results for the correlation energy later.
We now turn to the discussion of finite size errors in semiconductors and metals. We stress that small gap systems suffer from a relatively slow convergence of S(G) with respect to the studied system size. This behavior can be understood by considering the definition of S(G) Partial derivative of correlation energy with respect to Coulomb kernel for graphite. Slices and isosurface of S(G) for carbon graphite in the ABC stacking. Darker colors indicate more negative values. White corresponds to zero. The blue isosurface is a hexagonally shaped torus and reflects the anisotropic shape for S(G). Black dots represent the sampling points using a 4 × 4 × 4 k-point mesh.
in second-order Møller-Plesset perturbation (MP2) theory where i correspond to one-electron energies usually obtained from Hartree-Fock theory. The indices i, j and a, b label occupied and virtual orbitals respectively and are understood to be a shorthand for the Bloch wave vector k i and a band index n i . Due to momentum conservation k b can be calculated from the other Bloch wave vectors in the above equation. Γ ab ij (G) is defined in the supplementary information. The summation over Bloch vectors in Eq. (2) introduces quadrature errors that cause the slow convergence of S(G) towards the thermodynamic limit. In the case semiconductors or metals these errors can become significant because 1/( i + j − a − b ) varies strongly depending on k i , k j , k a and k b . In particular materials with a Dirac cone at the Fermi surface such as graphene exhibit a large variation of the denominator between zero and several eV depending on k. As a result Eq. (2) needs to be calculated using a finer k-point mesh to reduce quadrature errors. We will show in this work that the above quadrature errors can be substantially reduced by calculating and averaging S(G) for a set of shifted k-point meshes. Note that the vectors G are not affected by shifting the employed k-mesh because G depends only on the difference between any two Bloch wave vectors ∆k. We replace S(G) in Eq. (1) with an average obtained for N t different k-meshes shifted from Γ by t i such thatS The shifts t i are chosen such that they sample the first Brillouin zone uniformly. Coupled cluster theory calculations for different shifts can be performed independently from each other and the computational complexity scales only linearly with respect to N t . Coupled cluster theory energies that have been obtained using this twistaveraging technique will be referred to as CC-TA or CC-TA-FS if the interpolation method has been employed as well.
We note that quantum Monte Carlo (QMC) methods employ finite size corrections that share similarities with the methods outlined above [28][29][30]. However, QMC methods such as diffusion Monte Carlo are real-space theories that provide estimates of total energies rather than partitioning the energy into a Hartree-Fock and an electronic correlation contribution. An advantage of the partitioning ansatz is that Hartree-Fock energy contributions can be converged to the thermodynamic limit independently from the correlation energy at little extra computational cost. Consequently, finite size corrections are only required for the comparatively smaller correlation energy contributions. In passing we note that auxiliary field quantum Monte Carlo theory employs finite size corrections that are based on parametrized density functionals obtained from finite uniform electron gas simulation cells [31]. Such corrections work for solids but have not yet been applied to surfaces or molecular crystals where they are expected to be less accurate.

RESULTS
We now turn to the discussion of the results obtained using the methods outlined above. The present computations were performed using the VASP code [32,33] and the projector augmented wave method [34]. The coupled cluster theory calculations were partly performed using the newly developed cc4s code [35] interfaced with VASP and employing the automated tensor contraction engine CTF [36]. More technical details are outlined in the supplementary information.
As a first application we investigate the carbon diamond and graphite crystals. Before discussing the thermodynamic limit convergence we seek to address the convergence of the calculated correlation energy differences with respect to the employed orbital basis. We employ MP2 natural orbitals that are obtained using a procedure outlined in Ref. [37]. Figure 2 shows the convergence of coupled cluster singles and doubles (CCSD) and perturbative triples (T) correlation energy differences with respect to the number of bands using a 2 × 2 × 2 k-point mesh, respectively. We find that calculations using 16 orbitals per carbon atom yield an energy difference that agrees to within 4 meV/atom compared to results obtained using 40 natural orbitals per atom. The (T) correction to CCSD converges even faster with respect to the number of orbitals and is fortuitously close to zero in the case of the 2 × 2 × 2 k-point mesh. We stress that natural orbitals allow for a significantly more systematic truncatability and improved basis set incompleteness error cancellation between different systems compared to virtual Hartree-Fock or density functional theory orbitals. Achieving the same level of accuracy requires several hundred virtual Hartree-Fock orbitals per atom. Convergence with respect to other computational parameters has been checked and is discussed in the supplementary information.
We now turn to the discussion of finite size errors in total correlation energies. The top and middle panel in Figure 3 show CCSD correlation energies retrieved as a function of the number of k-points of graphite and diamond, respectively. We note that twist averaging (TA) is necessary for CCSD correlation energies to achieve a smooth 1/N convergence to the thermodynamic limit in particular for graphite. Accounting for quadrature errors by means of the tricubic interpolation method (CCSD-TA-FS) yields rapidly convergent correlation energies for both carbon diamond and graphite. CCSD-TA-FS correlation energies obtained using a 2 × 2 × 2 k-mesh only deviate from extrapolated CCSD thermodynamic limit energies by approximately 60 meV/atom. We note that correlation energies obtained using the same k-mesh and CCSD-TA theory exhibit finite size errors on the scale of 200-300 meV/atom. The CCSD(T) correlation energies are obtained using twist averaging for the (T) contribution and adding the correction to the CCSD-TA-FS result obtained using a 4 × 4 × 4 k-point mesh. This allows for investigating the finite size errors of the (T) correction independently from finite size errors of CCSD theory. We find that (T) converges rapidly with respect to the employed k-mesh size, reflecting its short-rangedness.
The bottom panel in Figure 3 shows the differences of the total energies of both carbon allotropes including zero point corrections retrieved as a function of the employed k-point mesh. Results obtained using CCSD theory without finite size corrections are depicted by the blue line and oscillate strongly with increasing kpoint mesh density. Using too coarse k-point meshes predicts graphite to be more stable than diamond, whereas denser k-point meshes predict diamond to be the more stable allotrope. Employing the averaging over different shifts yields CCSD-TA results that converge significantly smoother with increasing k-point mesh density as shown by the green line. Furthermore performing the newly proposed tricubic interpolation and integration in addition to the twist averaging referred to as CCSD-TA-FS yields rapidly convergent energy differences shown by the red line. We note that CCSD-TA-FS using a 2 × 2 × 2 k-point mesh is as close to the thermodynamic limit as CCSD-TA using a 4 × 4 × 4 k-point mesh. Since the computational complexity scales at least as O(N 4 k ) with respect to the number of k-points this corresponds to a reduction in the computational cost by three orders of magnitude. From these calculations we conclude I. Cohesive energies of solid Neon obtained using MP2, CCSD and CCSD(T) theory. The summarized results have been extrapolated to the complete basis set limit using pseudized aug-cc-pV(D,T)Z basis sets and corrected for basis set superposition errors using counterpoise corrections. All units in meV/atom. that CCSD theory predicts diamond to be more stable than graphite by 22 meV/atom including zero point energies. We have also performed perturbative triples calculations and added the corresponding correlation energy correction to our CCSD findings. CCSD(T) theory predicts graphite to be slightly less stable than diamond by 7 meV/atom including zero point corrections. We note that our CCSD(T) results agree with experimental findings to within the observed precision of CCSD(T) for similar applications which is generally better than 1 kcal/mol (43 meV/atom). The difference in the experimental Gibbs free energy of carbon diamond and graphite at room temperature has been reported to be 25 meV/atom [38], predicting graphite to be more stable than diamond.
Having demonstrated the ability of the proposed method to correct for finite size errors on the scale of about 100 meV/atom, we now seek to study the thermodynamic limit of the cohesive energy of the weakly bound Neon solid. In this case we need to correct for finite size errors on the scale of a few meV/atom. The dominant contribution to the attractive long range interatomic interaction of Neon atoms originates from van der Waals forces. Table I summarizes MP2, CCSD and CCSD(T) cohesive energies obtained with and without the proposed finite size correction. The correction yields MP2 and CCSD cohesive energies using 3×3×3 k-meshes that deviate from the thermodynamic limit results by approximately 1-2 meV/atom, whereas the uncorrected estimates deviate by 2-4 meV/atom. Although the finite size errors are small on an absolute scale, we stress that the corresponding relative finite size errors of the cohesive energy are non-negligible. Our best estimates for the MP2, CCSD and CCSD(T) cohesive energies using finite size corrections and a 4 × 4 × 4 k-mesh agree with results obtained using the incremental method to within 3 meV/atom [6,10]. Furthermore CCSD(T) predicts a cohesive energy of 30 meV/atom which is in good agreement with experimental estimates of 27 meV/atom corrected for zero point fluctuations [6].
As a final demonstration of the applicability of the proposed method to reach the thermodynamic limit we study the adsorption energy of a single water molecule on an h-BN sheet. The same system has recently been studied using diffusion Monte Carlo (DMC), the randomphase approximation (RPA) and dispersion functionals [39,40], as well as molecular MP2 [41] and periodic coupled-cluster theory [42], demonstrating the need for reliable methods that can account for long-range van der Waals interactions, also to provide benchmark data. Furthermore, the recent work of Al-Hamdani et al. [39] illustrates the importance of long-range correlation effects that account for approximately 25% of the reference adsorption energy computed in a (4 × 4) unit cell of h-BN. Figure 4 shows calculated adsorption energies at the level of RPA plus second-order screened exchange, MP2, CCSD and CCSD(T) theories retrieved as a function of the number of atoms in the h-BN sheet. Convergence with respect to other computational parameters has been checked and is discussed in the supplementary information. Using MP2 theory it is possible to study very large systems [39] and we find that the MP2 adsorption energy converges slowly to a thermodynamic limit value of 119 meV. We note that finite size errors for adsorption energies on two dimensional insulators are expected to decay as 1/N 2 , which is the predicted scaling from pairwise additive van der Waals interactions [39]. Applying the proposed finite size correction to MP2 theory for the (4 × 4) unit cell h-BN sheet with 32 atoms yields an adsorption energy of 113 meV in close agreement with the thermodynamic limit result. We observe a similar speedup in convergence using RPA+SOSEX-FS theory, illustrating the transferability of the proposed method. The water adsorption on the 32-atom h-BN sheet can also be studied using the more sophisticated CCSD theory [42]. CCSD with and without finite size corrections yields an adsorption energy of 83 meV and 68 meV, respectively. The finite size correction of MP2 and CCSD theory agree to within a few meV. However, we note that CCSD theory underbinds the water molecule. We estimate the (T) contribution using the 18-atom cell only and find that CCSD(T) theory yields an adsorption energy of 102 meV and 87 meV with and without finite size correction, respectively. The DMC adsorption energy was reported to be 84 meV without finite size corrections and agrees well with our CCSD(T) results using the same 32-atom h-BN sheet disregarding finite size corrections.

CONCLUSION AND OUTLOOK
In conclusion, we have introduced an efficient and accurate thermodynamic limit correction for wave function based theory calculations of solids and surfaces that is free of adjustable parameters and easy to implement. We have demonstrated that this correction allows for reducing the computational cost by several orders of magnitude without compromising accuracy. We have studied ground state problems where the convergence to the thermodynamic limit is crucial and finite size errors span a range of 1-100 meV/atom. Despite the local character of electronic correlation we stress that a proper treatment of long-range correlation effects is of paramount importance for reliable and highly accurate many-electron theories in condensed matter systems. We have applied the proposed finite size correction in combination with the gold standard of quantum chemistry CCSD(T) theory to calculate the cohesive energy of the Ne solid, the energy difference between carbon diamond and graphite crystals as well as the adsorption of a water molecule on an h-BN sheet. In general our CCSD(T) results are in good agreement with experimental findings and DMC results. This paves the way for a routine use of highly accurate coupled cluster theories in the field of surface science and solid state physics. We believe that the ability to predict accurate benchmark results will help the entire electronic structure theory community to improve further upon computationally efficient ab − initio theories and to help interpret experimental findings more reliably. To expand the scope of the proposed techniques even further we will aim at combining them with explicit correlation and low rank factorization methods [42][43][44].
In future studies we will extend the proposed finite size corrections to the study of excited states and metallic systems. We note that excited states and spectral functions can be calculated in the framework of equation of motion coupled cluster theory for solids, yielding excited state structure factors that are expected to exhibit similar finite size errors [45]. In metallic systems the structure factor is still algebraic around G = 0. We are therefore confident that the proposed methods can also be transfered to the study of such systems and we expect the outlined twist averaging methodology will be of significant importance when approaching the thermodynamic limit. We note, however, that the perturbative triples (T) contribution requires methodological improvements when applied to metals.