Coupled cluster theory was initially proposed by Fritz Coester and Hermann Kümmel in the field of nuclear physics (Coester, 1958; Coester and Kümmel, 1960). In the 1960s Jiri Cizek and Josef Paldus introduced the method for electron correlation (Cizek, 1966; Cizek and Paldus, 1971) and since then it has become a widely used electronic structure theory method for quantum chemical calculations on systems that do not exhibit strong static correlation (Bartlett and Musiał, 2007). For a more detailed introduction to coupled cluster theory we refer the reader to Helgaker et al. (2000); Bartlett and Musiał (2007); Crawford and Schaefer (2007); Shavitt and Bartlett (2009). Coupled cluster theory uses an exponential ansatz of cluster operators for the many-electron wavefunction

where Φ^HF is a single Slater determinant constructed from the Hartree-Fock (HF) one-electron orbitals that best approximates the ground state energy of a many-electron system. In passing we stress, however, that CC theory can in principle employ any single reference determinant, though this can affect the accuracy. The cluster operator in the exponent is defined by T^=∑m=1nT^m, where n corresponds to the order of the coupled cluster approximation. T^m is an m-fold excitation operator that generates m-fold excited Slater determinants |Φi1..ima1..am〉 multiplied by corresponding amplitudes (ti1..ima1..am) when applied to the HF ground state (|Φ^HF〉):

We will return to the discussion on the procedure for determining the amplitudes later. The mean-field wavefunction obtained from HF theory serves as a single reference determinant for the CC wavefunction and electronic correlation effects are accounted for explicitly using excitation operators and corresponding amplitudes. The indices i and a refer to occupied and unoccupied orbitals, respectively. We stress that depending on the choice of orbitals and their symmetry, the wavefunction amplitudes will reflect a different degree of sparsity. In the case of Bloch orbitals, as commonly used in periodic systems, the amplitudes vanish unless the sum of the momenta of the occupied orbitals is equal to the sum of the momenta of the unoccupied orbitals modulo a reciprocal lattice vector. We note that this property holds only for momentum-conserving Hamiltonians.

In coupled cluster singles and doubles theory the cluster operator is approximated using T^≈T^1+T^2. Due to the exponential ansatz of the CCSD wavefunction, the coefficients of all i-fold excited Slater determinants with i ≥ 2 are approximated using anti-symmetrized products of single and double excitation amplitudes. For two-electron systems, 2-fold excited Slater determinants can be generated at most and therefore CCSD theory becomes exact. Higher orders of coupled cluster theory become exact for systems with the corresponding number of electrons. The cluster amplitudes tia and tijab are obtained by inserting the wavefunction ansatz in the Schrödinger equation, projecting on the left by 〈Φia|e−T^ and 〈Φijab|e−T^, and solving the amplitude equations 〈Φia|e-T^HeT^|ΦHF〉=0 and 〈Φijab|e-T^HeT^|ΦHF〉=0, respectively. For practical computer implementations the latter equations are not useful but need to be recast in expressions that depend explicitly on Hamiltonian matrix elements and amplitudes, which can be achieved using further algebraic transformations including the Baker-Campbell-Hausdorff formula and Wick's theorem, Slater rules or diagrammatic techniques (Helgaker et al., 2000; Bartlett and Musiał, 2007; Crawford and Schaefer, 2007;Shavitt and Bartlett, 2009).

For most applications a good trade-off between computational cost and high accuracy is achieved by truncating the cluster operator at doubles (CCSD) and accounting for the effect of triples in a perturbative manner. This approach is referred to as CCSD(T) theory (Raghavachari et al., 1989) and achieves chemical accuracy for the prediction of reaction energies and barrier heights for a wide range of chemical reactions (Helgaker et al., 2000; Bartlett and Musiał, 2007). We note that CCSD and CCSD(T) theory include all diagrams for the correlation energy that occurs in third- and fourth-order perturbation theory, respectively. However, as a consequence of the single-reference approximation in CC theory, the treatment of strong correlation problems is extremely limited. Examples for strong correlation problems include molecular dissociation problems where the HF approximation to the wavefunction fails dramatically due to the multideterminant nature of the true wavefunction. Conventional single-reference CC theory are accurate only for wavefunctions that are dominated by single-reference determinants. We note in passing that multireference coupled cluster theories aim at the treatment of strong correlation problems (Köhn et al., 2013; Evangelista, 2018).

To establish a connection between coupled cluster methods and the random-phase approximation, we now turn to coupled cluster doubles (CCD) theory. CCD theory approximates the cluster operator using T^≈T^2. The cluster amplitudes tijab are obtained by solving the quadratic amplitude equations 〈Φijab|e-T^2HeT^2|ΦHF〉=0 that in a canonical spin-orbital basis read

In the above equation we use Einstein summation convention. The indices i, j, k, and l label occupy orbital indices, whereas a, b, c, and d label virtual orbital indices. ϵ corresponds to the HF one-electron energies and the anti-symmetrized electron repulsion integrals (ERIs) are defined as 〈ij||ab〉=〈ij|r12-1|ab〉-〈ij|r12-1|ba〉, where the ERIs are defined by

In the above expression the spin-orbitals χ depend on the space-spin coordinate x = (r, σ) and the spatial coordinates are integrated over all space. In section 3 we will return to the discussion of the evaluation of these integrals using different frameworks of the independent particle basis sets. Equation (3) is solved for the amplitudes in an iterative manner by updating the amplitudes in every iteration using the right-hand side of Equation (3). To accelerate convergence standard techniques such as direct inversion of the iterative subspace (DIIS) can be employed (Pulay, 1980). Once the amplitudes are obtained, the CCD correlation energy can be calculated by

As a consequence of the large number of virtual orbitals, N_v, compared to the number of occupied orbitals, N_o, needed to obtain converged correlation energies even when extrapolation techniques are employed, the computation and storage of ERIs as well as their contraction with amplitudes constitute the main source of computational cost and memory in canonical CC calculations. In terms of computational cost and memory the so-called particle-particle ladder term 〈cd||ab〉tijcd exhibits the most unfavorable scaling with respect to the number of virtual orbitals O(Nv4No2). In passing we note that a number of techniques have been suggested to reduce the computational cost of the 〈cd||ab〉tijcd term in CC calculations (Hummel et al., 2017; Dutta et al., 2018) and that local schemes can treat this term very efficiently through construction, since it couples a single electron-pair to itself only. The distribution of memory and computational load in parallel computer implementations as well as the possible on-the-fly calculation of the required ERIs make it possible to study systems containing several hundreds of virtual orbitals routinely. Furthermore, modern tensor framework libraries greatly simplify the development of compact and also fully parallel coupled cluster theory computer codes (Solomonik et al., 2014).

Equation (3) arranges the terms of the amplitude equations in the same manner as Shepherd et al. (2014), who investigated different approximations to coupled cluster doubles theory for the uniform electron gas. The various terms are labeled in agreement with the corresponding diagrams such that the top line of Equation (3) represents the driver and the ring terms, the second line ladder terms, the third line crossed-ring terms, and the bottom-line mosaics (Scuseria et al., 2013). The quadratic amplitude equations of coupled cluster doubles theory couple all electron pairs of a many-electron system using perturbation theory diagrams of a certain type to infinite order. This is illustrated by the iterative solution for the amplitudes in Equation (3): in the first iteration it follows that tijab=〈ij||ab〉ϵi+ϵj-ϵa-ϵb. In subsequent iterations the approximate first-order amplitudes are coupled via the corresponding ladder, ring, crossed-ring, and mosaic terms to each other. In this manner a summation of all possible couplings between different diagrams to the infinite order is performed. The included diagrams balance and account for important physical effects of many-electron systems. All terms are accompanied by so-called exchange terms that are needed to correct for exclusion principle violating (EPV) contributions and account for the fermionic character of the many-electron wavefunction. EPV contributions cause self-interaction errors that are not necessarily restricted to approximate DFT methods. Furthermore, the magnitude of ring and ladder contributions to the electronic correlation energy is system dependent. For homogeneous systems the dimensionality and electronic density plays a crucial role. In the high-density limit the ring approximation becomes exact for three-dimensional systems, whereas ladder terms become more important for lower-densities and lower dimensions (Freeman, 1977, 1978, 1983; Bishop and Lührmann, 1978). The inclusion of ring diagrams is important to describe collective polarization effects and obtain effectively screened inter-electronic interactions that yield convergent correlation energies in metals.

Coupled cluster theory is based on a systematically improvable many-electron wavefunction ansatz and derived without using uncontrolled approximations. However, various approximations to the amplitude equations have been investigated and have proven to sometimes be more accurate and even more stable than their parent method. The distinguishable cluster approximation (DCSD) (Kats and Manby, 2013; Kats et al., 2015) demonstrates that disregarding certain diagrams in the amplitude equations and reweighting others such that the resultant method remains exact for two-electron systems yields results that clearly outperform CCSD in terms of accuracy. Furthermore, it is possible to show that by a change of representation DCSD exhibits a computational complexity that scales with respect to system size only as O(N5) instead of O(N6) as in the case of CCSD (Mardirossian et al., 2018). This illustrates that the key to efficient and accurate many-electron perturbation theories lies in carefully chosen and well-balanced truncations of the many body perturbation series. In this regard the random-phase approximation is a showcase for an efficient and reasonably accurate perturbation theory.

The random phase approximation (RPA) to the correlation energy dates back to the 1950s. It was first introduced by Macke to predict convergent correlation energies (Macke, 1950) in the uniform electron gas and was also developed by Pines and Bohm (1952) for the collective description electron interactions. In the field of ab-initio computational materials science the exact-exchange plus correlation in the random-phase approximation has attracted renewed and widespread interest in the last two decades (Paier et al., 2012). This is due to the fact that computationally increasingly efficient implementations have become available making this method capable of describing all interatomic bonding situations reasonably well: ionic, covalent, metallic, and even van der Waals bonding. The computational complexity can even be lowered to O(N3) in real space formulations (Kaltak et al., 2014). Thus, the complexity of an RPA calculation does not exceed that of a canonical hybrid density functional theory calculation, the pre-factor is however considerably larger. The RPA correlation energy can be derived from many-electron Green's function theory, or using the adiabatic-connection fluctuation-dissipation theorem (ACFDT), or from coupled cluster theory.

As shown in Scuseria et al. (2008), it is possible to transform the RPA equations, that are usually expressed in a general eigenvalue problem, to a quadratic Riccati equation that reads

We stress that in the above equation ϵ corresponds to the DFT one-electron energies. Once the amplitudes are obtained, the RPA correlation energy can be calculated by

Although the above formulation does not allow for an efficient computer implementation of the RPA, it illustrates that the RPA and CCSD are closely related. In the rings-only approximation, the second, third and fourth lines of Equation (3) are disregarded. Furthermore, the random-phase approximation includes the direct rings only. This implies that instead of using the (double bar) anti-symmetrized integrals, only 〈ij|ab〉 integrals are employed in the RPA amplitude and energy equations, making it necessary to employ a different pre-factor in the correlation energy expression to stay consistent with the many-body perturbation theory. Consequently, RPA cannot be viewed as a wavefunction theory, though it can be obtained from the CC amplitude equations as explained above. The close relationship between these approaches has motivated a number of post-RPA corrections, of which we will discuss only a small selection.

In the RPA, the magnitude of total correlation energies is significantly overestimated and binding energies are systematically underestimated compared to experiment even for weakly interacting systems. The origin of these shortcomings can in some cases be understood by comparing the RPA amplitude and energy expressions in Equations (6, 7) to the coupled cluster doubles amplitude and energy expressions in Equations (3, 5), respectively.

The lack of exchange-like terms in the RPA leads to the inclusion of EPV contributions (such as 〈ii|ab〉tiiab) causing self-correlation errors; for example, RPA yields non-zero electron correlation energies for one-electron systems. This has motivated the introduction of various post-RPA corrections, such as second-order screened exchange (SOSEX), AC-SOSEX, or approximate exchange kernel methods (AXK). The inclusion of SOSEX to the random-phase approximation dates back to Monkhorst and Freeman in the 70s (Monkhorst and Oddershede, 1973; Freeman, 1977). Freeman showed in Freeman (1977) that absolute correlation energies of the uniform electron gas become significantly improved when compared to more accurate methods if SOSEX correlation energy contributions are included. Furthermore, the underestimation of cohesive energies calculated in the RPA originates from the fact that EPV contributions are larger for spin-polarized atoms than for non-spin-polarized solids. Therefore, the RPA + SOSEX approximation yields more accurate cohesive energies compared to experiment than the RPA. We note that similar findings have been obtained using related corrections to the RPA, such as AXK (Chen et al., 2018).

Another important difference between RPA and coupled cluster theory is the choice of the reference determinant and orbital energies. In the framework of ACFDT-RPA, the KS-DFT reference determinant is well-justified. However, from the perspective of the CC theory the KS-DFT orbitals violate the Brillouin's theorem and the correlation energy expression would therefore have to include terms that depend on the off-diagonal Fock matrix elements. These considerations have motivated the renormalized singles excitations (rSE) and related methods (Paier et al., 2012; Ren et al., 2013; Klimeš et al., 2015). The inclusion of these contributions improves upon the description of weakly-bound molecules and solids compared to the RPA.

In this section we discuss size extensivity and methods that accelerate the convergence to the thermodynamic limit (TDL). Both concepts are highly relevant for the application of quantum chemical methods to solids. In contrast to molecular systems, properties of solids or surfaces need to be calculated in the TDL in order to allow for a direct comparison to experiment and truly model an infinite periodic system. The TDL can be approached in different ways using; for example, (i) sampling of the Brillouin zone with increasingly dense k-point meshes and in periodic boundary conditions, (ii) studying increasingly large supercells in periodic boundary conditions, or employing (iii) increasingly large clusters with open boundary conditions and/or embedding methods. Once the TDL is approached with respect to the number of k-points or the number of atoms in the cluster, intensive properties, such as the correlation energy per atom are converged to a constant value.

An important advantage of truncated coupled cluster theories compared to truncated configuration interaction methods is their size extensivity. Size extensivity is a concept of particular importance in quantum chemistry, which judges if the calculated quantities have the correct asymptotic size dependence or not. For extensive quantities, like the (correlation) energy, a given size extensive method should yield the asymptotic N¹ dependence with N being the number of unit cells or wave vector sampling points in the Brillouin zone (Hirata, 2011). Obviously, the methods with incorrect asymptotic N^α dependence of α < 1, like the truncated configuration interaction methods, lead to the total energy per unit cell equal to that of the HF mean-field approximation, which is clearly useless for condensed-matter systems. The size extensivity of coupled cluster theories can also be understood via either the diagrammatic criteria (Bartlett, 1981) or the super-molecule criterion (Szabo and Ostlund, 1996). Moreover, it was argued that approximate post-HF correlation methods cannot capture the variational and size-extensive properties simultaneously (Hirata and Grabowski, 2014).

The TDL is approached as the number of particles becomes infinite in the simulation (super-)cell while the density is kept constant. Once the TDL is approached, correlation energies per atom need to be converged to a constant for periodic systems, corresponding to α = 1. Finite size errors are defined as the difference between the TDL and the finite simulation cell results. However, converging calculated properties with respect to the system size can be very slow, requiring substantial computational resources due to the steep scaling of the computational complexity of coupled cluster methods with respect to system size. We stress that many properties, such as the binding energy of molecules on surfaces converge even slower than their counterparts calculated on the level of mean-field theories, such as DFT. Correlated wavefunction based methods capture long-range electronic correlation effects; for example, dispersion interactions, explicitly. Although the respective contribution to the electronic correlation energy can be small, the accumulation of such interactions can become a non-negligible contribution to the property of interest. Various different strategies have been developed to correct for finite size errors that are defined as the difference between the TDL and the finite simulation cell results. These strategies often involve extrapolation methods or range-separation techniques. Local theories that employ correlation energy expressions depending on localized electron pairs, can approximate correlation energy contributions of long-distant pairs using computationally more efficient yet less accurate theories. Alternatively, local theories can account for electron pairs that are disregarded based on a distance criterion by using an R⁻⁶-type extrapolation (Usvyat et al., 2012). Canonical implementations of periodic post-HF methods employ scaling laws for extrapolations to the TDL that are based on an analog rationale (Del Ben et al., 2013; Booth et al., 2013; McClain et al., 2017). 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 (Kwee et al., 2008).

We stress that the problem of slow thermodynamic limit convergence and concomitantly large finite size errors is a common feature of Quantum Monte Carlo (QMC) and many-electron quantum chemical methods due to their explicit treatment of electronic correlation effects. Recently structure factor interpolation and twist averaging (TA) techniques have been introduced in Gruber et al. (2018) to correct efficiently for finite size errors in CC calculations of solids. These techniques are inspired by related methods used in QMC calculations (Chiesa et al., 2006; Holzmann et al., 2016). These finite size (FS) corrections allows for thermodynamic limit results of solids and surfaces to be achieved using quantum chemical wavefunction theories in a very efficient manner (Liao and Grüneis, 2016; Gruber et al., 2018), reducing the computational cost significantly. As shown in Figure 1, it is possible to achieve a much more rapid TDL convergence for simple solids, such as carbon graphite and diamond by virtue of these corrections. Furthermore, the fact that the converged energies per atom remain constant as the k-mesh density is increased, reflects that CC theories are size extensive.

In this section, we briefly review the calculation of two-electron repulsion integrals in the framework of the projector augmented wave (PAW) method using a plane wave basis set. These integrals are the most important quantities in addition to the cluster amplitudes for coupled cluster theory calculations. In the PAW method, the all-electron orbitals (|ψ_n〉) are obtained from the pseudo orbitals (|ψ˜n〉) using a linear transformation (Blöchl, 1994),

The index n, labeling the orbitals ψ, is understood to be shorthand for the band index and the Bloch wave vector k_n, while the index i is a shorthand for the atomic site R_i, the angular momentum quantum numbers l_i and m_i, and an additional index ϵ_i denoting the linearization energy. The wave vector is conventionally chosen to lie within the first Brillouin zone. The pseudo orbitals are the variational quantities of the PAW method and are expanded in reciprocal space using plane waves, 〈r|ψ˜n〉=1Ω∑GCGnei(kn+G)r. In this framework, it is possible to evaluate two-electron repulsion integrals approximately using the following expression:

where ṽ is the diagonal Coulomb kernel in reciprocal space 4πG2. We note the reciprocal Coulomb kernel exhibits a singularity a G = 0. We employ a correction for the Coulomb kernel at the singularity that is obtained using a scheme introduced by Gygi and Baldereschi (1986). In the present PAW implementation, the Fourier transformed codensities ρia(G) are approximated using Equation (2.87) of Harl (2008) as originally implemented by Kresse et al. for the calculation of correlation energies within the random phase approximation (Harl and Kresse, 2008). The codensities and the diagonal Coulomb kernel make it possible to calculate the ERIs in a computationally efficient on-the-fly manner. We note that the memory footprint of the transformed codensities scales cubic with respect to the system size. Furthermore, it is also possible to further reduce the dimension of the auxiliary plane-wave index G significantly for system where the large number of degrees of freedom provided by the plane-waves are not needed to describe the employed codensities. A significant down-folding of the plane wave basis set size is possible without compromising accuracy for calculations of surfaces as well as atoms or molecules in a box. A convenient and computational efficient method to achieve this downfolding is based on singular value decomposition and outlined in Hummel et al. (2017). We note that our more recent implementation of CC theory that calculates the ERIs according to Equation (2) (cc4s) employs the cyclops tensor framework (CTF, Solomonik et al., 2014) and will be released in the near future.

The implementation of quantum-chemistry methods in the numeric atom-centered orbital (NAO) framework utilizes the resolution-of-identity (RI) technique (also known as “density fitting”) to calculate the two-electron repulsion integrals which scale as N⁴ in memory with respect to system size N. The key idea is to decompose the four-rank repulsion integrals in terms of three-rank tensors with an N³ scaling of memory and is therefore better suited to be pre-stored:

where the index μ runs over an auxiliary basis {P_μ(r)}, and mprμ is a decomposed three-rank tensor in the molecular orbital basis. In RI-V approximation, mprμ is determined by directly minimizing the errors in the repulsion integral of atomic orbitals using Equations (47, 55) of Ren et al. (2012). A hybird-RI algorithm has been designed to balance the computing cost and communication in the CCSD(T) implementation for molecules in the Fritz Haber Institute ab initio molecular simulation (FHI-aims) package. Together with a domain-based distributed-memory strategy, it allows for an effective utilization of the quickly increasing memory bandwidth of today's supercomputers to avoid the on-the-fly disk storage and minimize interconnect communication, particularly for the tensor contraction in the evaluation of the particle-particle terms. As a result, an excellent strong scaling can be achieved up to over 10,000 cores (Shen et al., 2018). The parallel efficiency is competitive with the CC implementations in state-of-the-art high-performance computing computational chemistry packages.

Furthermore, a local variant of RI-V, namely RI-LVL, has been developed in the numeric atom-centered orbital framework. Unlike the standard RI-V approximation in Equation (10), RI-LVL expands the products of basis functions only in the subset of those auxiliary basis functions {P(IJ)} which are located at the same atoms of I and J as the basis functions of χ_i and χ_j:

where the three-center (rp|λ) and two-center (μ|ν) integrals are defined as

and

respectively. L^RP = (μ|ν)⁻¹ with μ, ν∈P(RP). This choice further reduces the memory scaling to N² without compromising the accuracy (Ihrig et al., 2015). For molecules, it is a useful alternative to the standard RI-V, which, however, is the only option for solids, since the memory consumption of RI-V quickly becomes unaffordable with the increase of k-grid numbers in periodic boundary conditions (Levchenko et al., 2015). For MP2 and RPA, we demonstrated that the NAO-based periodic implementation can provide a smooth and consistent convergence toward the complete k-grid limit (Zhang et al., 2018). The periodic CCSD in FHI-aims has been implemented, but not yet fully optimized.

We now turn to a brief overview of ab-initio computational materials science studies that have been performed using periodic coupled cluster theory. This overview is by no means complete but reflects the general purpose and scope of the CC implementation based on the PAW formalism and its relation to other approaches, such as the method of increments. A selection of the investigated systems is depicted in Figure 2. PAW-based periodic CC calculations have initially been limited to insulating solids with two-atomic unit cells, such as the LiH crystal or the uniform electron gas simulation cell (Grüneis et al., 2011; Booth et al., 2013; Shepherd and Grüneis, 2013). However, recent methodological improvements make the study of surfaces and supercells containing approximately thirty atoms possible. In particular the development of finite size corrections (Gruber et al., 2018), improved compact approximations to the virtual orbital manifold (Grüneis et al., 2011; Booth et al., 2016) and auxiliary basis sets (Hummel et al., 2017) have substantially expanded the scope of periodic CC calculations. The advancement of explicitly correlated techniques in combination with a plane wave basis set holds the promise to reduce the computational cost further in the near future (Grüneis et al., 2013, 2017).

We first review studies that aim at assessing the precision of periodic CC calculations. An important task of method development includes demonstrating convergence of calculated properties, such as the ground state energy with respect to the employed computational parameters; for example, basis set and k-mesh, and achieve agreement with results obtained using entirely different computational approaches. In this regard the calculation of the cohesive energy of the LiH solid on the level of HF, MP2, and CC theory has become common practice for benchmarking various implementations. The first quantum chemical calculations of the LiH crystal were performed by Egil Hylleraas in 1930 (Hylleraas, 1930). In 2009, Nolan et al. have presented MP2 and coupled cluster ground state energy calculations of the LiH crystal converged with respect to basis set and system size through a combination of periodic and finite-cluster electronic structure calculations (Nolan et al., 2009). These findings have served as a reliable reference for other approaches including (local) periodic and incremental methods used to calculate MP2 or CC cohesive energies (Marsman et al., 2009; Grüneis et al., 2011; Stoll and Doll, 2012; Del Ben et al., 2013; Usvyat, 2013). Table 1 summarizes the obtained cohesive energies, indicating the achieved level of precision of the various calculations is better than chemical accuracy. Furthermore, the CCSD(T) estimates agree with the experimental value (4.98 eV/LiH) to within a few ten meV. In a similar context, the adsorption energy of a single water molecule on the LiH surface has been studied and results have been compared among different theories and implementations to test their precision and accuracy (Tsatsoulis et al., 2017). The corresponding adsorption energies are summarized in Table 2 and demonstrate an agreement of a few meV between different implementations. Such benchmark calculations are useful despite their lack of reference to experiment. However, to further assess the scope of an implementation it is necessary to investigate different bonding situations. Weakly bound molecular and rare-gas crystals have been studied intensively using the incremental method with high accuracy (Rościszewski et al., 1999; Schwerdtfeger et al., 2010). These systems are composed of many weakly interacting fragments, making them ideally suited for an expansion of the correlation energy into few-body incremental contributions. In contrast to the incremental method, a particular challenge for fully periodic and canonical approaches originates from the long-rangedness of correlation energy contributions to the binding energy, requiring dense k-mesh calculations or reliable finite size corrections (Gruber et al., 2018). However, by using recently proposed finite size corrections, it is possible to achieve well-converged CC cohesive energies of the neon solid in good agreement with results obtained using the method of increments despite employing relatively coarse k-meshes (Schwerdtfeger et al., 2010; Gruber et al., 2018). The results obtained for the cohesive energy of the neon solid are summarized in Table 3. The periodic CCSD(T) estimate for the cohesive energy agrees to within 3 meV/atom with experimental value of 27 meV/atom, which has been corrected for zero-point fluctuations (Rościszewski et al., 1999).

To benchmark the accuracy of periodic ab-initio methods in a systematic manner it is common practice to calculate cohesive energies, lattice constants and bulk moduli for a range of solids exhibiting different chemical bonding including van-der Waals, ionic, covalent and metallic systems. In molecular quantum chemistry similar calculations are routinely performed for a range of test sets that reflect a range of different chemical bonding situations as discussed in the following section. In solids, the systematic comparison to experimental findings, which have been corrected for beyond Born-Oppenheimer approximation effects, allows to identify and understand systematic errors for different levels of theory; for example, MP2 theory overestimates correlation energies for systems with small band gaps, which results in an overestimation of corresponding cohesive energies (Grüneis et al., 2010). However, the computational cost of coupled cluster calculations using the PAW-based periodic implementation limits the number of benchmark systems so far. The range of simple solids for which cohesive energies have been calculated in the complete basis set limit on the level of CC theory include LiH (rock-salt), Ne (fcc), C (diamond), BN (zinc-blende), and AlP (zinc-blende) (Booth et al., 2013; Gruber et al., 2018). Furthermore, phase transitions and energy differences of different LiH, C, and BN allotropes have also been investigated (Grüneis, 2015a; Gruber and Grüneis, 2018; Gruber et al., 2018). The study of pressure-temperature phase diagrams allows to investigate the relative level of accuracy of electronic structure theories for different chemical bonds that are present in the investigated phases. As such the predicted equilibrium phase boundaries often change significantly with respect to the employed electronic structure theory. In Gruber and Grüneis (2018) we have employed periodic CCSD(T) theory to calculate relative enthalpies and pressure-temperature phase diagrams for C and BN allotropes. In contrast to currently available DFT methods, quantum chemical wavefunction theories allow for achieving reliable and systematically improvable benchmark results for the relative stability of C and BN allotropes.

Another promising area of application for quantum chemical wavefunction theories is the study of surface science problems including molecular adsorption and reactions on surfaces. CC theories are widely used for predicting benchmark results for molecular gas-phase reaction energies as well as activation barrier heights. A similar level of accuracy for molecular surface reactions can currently only be achieved by QMC calculations and it would be advantageous to fully transfer quantum chemical wavefunction based methods to this area. Currently available QC approaches are based mostly on embedding, fragment models, the incremental technique or finite-cluster calculations (Voloshina et al., 2011; Libisch et al., 2012; Usvyat et al., 2012; Boese and Sauer, 2016; Kubas et al., 2016). Recently we have performed fully periodic calculations of adsorption energies for water on h−BN and graphene sheets, confirming the expected level of accuracy (Gruber et al., 2018; Brandenburg et al., 2019). Table 4 summarizes the results of the water adsorption energies and compares them to DMC findings that agree well. Furthermore, the dissociative H₂ adsorption on the Si(001) surface has also been investigated in Tsatsoulis et al. (2018). The obtained results indicate that CC theory can predict accurate adsorption energies and study chemical reactions on surfaces in a reliable manner. However, we stress that, despite the accurate findings discussed above, the level of accuracy always depends on the electronic structure of the system and the level of many-electron wavefunction approximation. Therefore, high accuracy can only be achieved using CCSD(T) theory for systems that do not exhibit strong correlation effects.

As an alternative basis set choice instead of GTOs and plane waves, NAOs, in particular those with valence correlation consistency namely NAO-VCC-nZ, hold the promise to provide improved description of advanced correlation methods with the increase of the basis set size. NAO-VCC-nZ allows for extrapolating the results of advanced quantum-chemistry methods to the complete-basis-set (CBS) limit. NAO-VCC-nZ was generated by minimizing frozen-core RPA total energies of individual atoms from H to Ar. The consistent convergence of RPA and MP2 binding energies of small molecules is demonstrated by using NAO-VCC-nZ basis sets in the original paper (Zhang et al., 2013). The applicability of NAO-VCC-nZ for solids has been benchmarked comprehensively for MP2 and RPA calculations of cohesive energy, lattice constants, and bulk modulus for representatives of first- and second-row elements and their binaries with cubic crystal structures and various bonding characters (Zhang et al., 2018). The generalization of the periodic RPA and MP2 implementation to CCSD(T) for solids is straightforward, but strong effort should be paid in order to achieve an efficient and practical implementation of the NAO-based periodic CCSD(T) method, which is still under development.

As the first step, the CCSD(T) has been implemented in FHI-aims for molecules, which has been used together with NAO-VCC-nZ basis sets to deliver accurate CCSD(T) results in the CBS limit for manifold properties of great chemical or physical interest. For 22 bio-orient weak interactions in the S22 test sets and 10 relative energies of cysteine conformers in the CYCONF test set, the CCSD(T)/CBS results by using NAO-VCC-nZ repeat the up-to-date reference data based on GTO basis sets with the deviation of <0.1 kcal/mol on average. For the first time, the high-level theoretical reference data at the CCSD(T)/CBS level has been generated for the 34 isomerization energies in the ISO34 test set in the NAO framework. For the sake of benchmarking newly developed electronic-structure methods, the use of the CCSD(T)/CBS reference data allows for the comparison based on exactly the same molecular geometry and immune to the experimental uncertainty (Shen et al., 2018).

During the last years the scope of quantum chemical wavefunction theories and in particular coupled cluster methods have been expanded significantly in the field of complex systems and solids. Fully periodic and canonical approaches have become increasingly efficient due to methodological improvements that reduce the pre-factor of the computational cost. In this tutorial-style review we have given a brief overview of different frameworks for the implementation of periodic CC theories and their applications. Future work will focus on more systematic benchmark studies employing all recent methodological advancements including finite size corrections and explicit correlation techniques for solids. Furthermore, it is expected that local approaches, such as the use of truncated pair natural orbitals (Neese et al., 2009; Kubas et al., 2016; Ma et al., 2017) can also be transferred to solid state systems. The remaining technical and theoretical challenges are, however, significant and include the implementation of sparse tensor frameworks and the development of localization schemes that work reliable for metals as well as for insulators and semiconductors. Furthermore, the calculation of additional properties, such as gradients will eventually also be necessary in order to allow for structural relaxation. However, the high level of accuracy makes quantum chemical wavefunction theories and in particular CC theory a promising tool for future ab-initio computational materials science studies that can complement currently available DFT approaches and provide benchmark results. As discussed in this article, possible areas of application include the study of pressure-driven phase transitions, surface chemistry problems as well as optical properties of solids.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer DU declared a past co-authorship with one of the authors AG to the handling editor.