On the duality of ring and ladder diagrams and its importance for many-electron perturbation theories

TL;DR

This paper reveals the duality of ring and ladder diagrams in many-electron perturbation theories, showing how their understanding can improve computational efficiency and accuracy in electronic structure calculations.

Contribution

It introduces a diagrammatic decomposition highlighting the duality of ring and ladder diagrams and demonstrates an effective approximation for ladder diagrams using second-order perturbation theory.

Findings

Ladder diagrams can be approximated by second-order perturbation theory.

The approach reduces coupled cluster calculation costs by two orders of magnitude.

The approximation maintains accuracy in electronic correlation energy calculations.

Abstract

We present a diagrammatic decomposition of the transition pair correlation function for the uniform electron gas. We demonstrate explicitly that ring and ladder diagrams are dual counterparts that capture significant long- and short-ranged interelectronic correlation effects, respectively. Our findings help to guide the further development of approximate many-electron theories and reveal that the contribution of the ladder diagrams to the electronic correlation energy can be approximated in an effective manner using second-order perturbation theory. We employ the latter approximation to reduce the computational cost of coupled cluster theory calculations for insulators and semiconductors by two orders of magnitude without compromising accuracy.

Tables

Table 1. Table 1: CCSD and CCSD-PPL valence electron correlation energies for the Ne atom in m E h subscript 𝐸 ℎ E_{h} . The CBS limit value is -315.7 m E h subscript 𝐸 ℎ E_{h}   Shiozaki  et al. ( 2008 ) . DZ, TZ and QZ refer to the employed atom-centered Gaussian basis sets aug-cc-pVXZ (X=D,T and Q). The MP2 CBS limit energy was computed using a [Q5] extrapolation.

	DZ	TZ	QZ
CCSD	-210.15	-274.09	-297.76
CCSD-PPL	-308.08	-314.95	-317.27

Equations

$$E_{\mathrm{c}}=\int \text{d}{\mathbf{r}}_{12}g({\mathbf{r}}_{12})v({\mathbf{r}}_{12}).$$

$$g({\mathbf{r}}_{12})=\frac{1}{(2\pi {)}^3}\sum _{abij}{\mathrm{e}}^{\mathrm{i}\left({\mathbf{k}}_i-{\mathbf{k}}_a\right)\cdot {\mathbf{r}}_{12}}(2t_{ij}^{ab}-t_{ji}^{ab})$$

$$\begin{array}[]{r c rc cc cc cc cc cc cc ccc}&&&\overbrace{\rule{25.83325pt}{0.0pt}}^{\text{driver}}&&\lx@intercol\hfil\overbrace{\rule{103.33301pt}{0.0pt}}^{\text{ring}}\hfil\lx@intercol&&\overbrace{\rule{34.44434pt}{0.0pt}}^{\text{ppl}}&&\lx@intercol\hfil\overbrace{\rule{146.38843pt}{0.0pt}}^{\text{rest}}\hfil\lx@intercol\\ \raisebox{2.15277pt}{\includegraphics{T}}&=&&\raisebox{2.15277pt}{\includegraphics{V}}&+&\raisebox{2.15277pt}{\includegraphics{VT}}&+&\raisebox{2.15277pt}{\includegraphics{TV}}&+&\raisebox{2.15277pt}{\includegraphics{TVpp}}&+&\raisebox{2.15277pt}{\includegraphics{TVT}}&+&\raisebox{2.15277pt}{\includegraphics{TVph}}&+&\ldots\\[25.83325pt] t^{ab}_{ij}&=&(&v^{ab}_{ij}&+&v^{ak}_{ic}t^{cb}_{kj}&+&v^{kb}_{cj}t^{ac}_{ik}&+&v^{ab}_{cd}t^{cd}_{ij}&+&v^{kl}_{cd}t^{ac}_{ik}t^{db}_{lj}&-&v^{kb}_{id}t^{ad}_{kj}&+&\ldots&)&/&\Delta^{ij}_{ab}\end{array}$$

$$\begin{gathered} g^{\mathrm{CCSD}}({\mathbf{r}}_{12})= \\ {g}^{\mathrm{driver}}({\mathbf{r}}_{12})+g^{\mathrm{ring}}({\mathbf{r}}_{12})+g^{\mathrm{ppl}}({\mathbf{r}}_{12})+g^{\mathrm{rest}}({\mathbf{r}}_{12}). \end{gathered}$$

$$S^{\mathrm{ring}}(\mathbf{q})$$

$$S^{\mathrm{ppl}}(\mathbf{q})$$

$$\begin{gathered} E_{\mathrm{c}}^{\mathrm{CCSD}-\mathrm{PPL}}(N) \\ =E_{\mathrm{c}}^{\mathrm{CCSD}}(N)+\Delta E^{\mathrm{driver}}(N)+\Delta E^{\mathrm{ppl}}(N), \end{gathered}$$

$$\Delta E^{\mathrm{driver}}(N)$$

$$\Delta E^{\mathrm{ppl}}(N)$$

Full text

particle-particle ladder diagram; ring diagrams; coupled cluster; random phase approximation

Introduction. – The ongoing advancement of modern electronic structure theories over the last few decades led to a proliferation of computational studies and research on atoms, molecules and condensed matter systems based on first principles. Electronic structure methods allow to simulate a wide range of physically as well as chemically relevant material properties including thermodynamic stability, electric conductivity, magnetic ordering and optical properties. Yet many of the most successful methods such as approximate density functional theory or many-electron perturbation theories rely substantially on fortuitous error cancellation with far-reaching consequences. Attempting to improve upon individual terms in the many-body perturbation expansion of the electronic exchange and correlation energy by employing higher levels of theory often has a contrary effect and deteriorates the achieved level of accuracy, making it hard to approach the exact result in a systematic manner. A profound physical understanding of the individual perturbation theory contributions to the ground- and excited state energies of many-electron systems and their balance is, however, needed to design novel, more efficient and concomitantly more accurate theories. Here, we present a new method that makes it possible to decompose the many-body contributions to the correlation energies of ab-initio systems as well as the uniform electron gas in a manner that enables a deeper understanding of widely-used approximations. Furthermore we demonstrate that the obtained findings can be used to improve the efficiency of many-electron theory calculations.

In many-body physics and quantum chemistry, approximate methods are often tailored to become exact in limiting cases. Examples for such electronic structure theories include: (i) the random-phase approximation (RPA), which captures the leading order contributions to the correlation energy of the uniform electron gas in the high density limit Macke (1950); Pines and Bohm (1952); Gell-Mann and Brueckner (1957), (ii) the ladder theory which works reliable in the low-density limit Freeman (1983); Yasuhara (1972, 1974), and (iii) coupled cluster singles and doubles theory (CCSD) that combines the latter two approaches and is exact for two-electron systems Coester and Kümmel (1960); Čížek (1969). Yet all of the methods mentioned above exhibit shortcomings for real materials, respectively: systematic overcorrelation, divergence of the correlation energy for metals, and a poor trade-off between computational cost and accuracy. In a remarkable manner, however, errors like the overcorrelation of the RPA, cancel out when considering energy differences, and sometimes even mimic more sophisticated electronic correlation effects Henderson and Scuseria (2010); Harl and Kresse (2009); Dobson et al. (2006). Furthermore disregarding exchange-like interactions between electron pairs in coupled cluster theory systematically improves upon the achieved level of accuracy Kats and Manby (2013). These examples illustrate that it is imperative to fully understand the effect of individual contributions as a function of the electronic density and interelectronic distance in order to guide the further development of more accurate and efficient approximations to the electronic correlation energy in strongly as well as weakly correlated systems.

Theory. – In this work we introduce a decomposition of the electronic transition pair correlation function and the corresponding correlation energy contributions. For the sake of brevity we will restrict the following discussion to the uniform electron gas (UEG) and turn to ab-initio systems afterwards. The correlation energy can be written as an integral of the Coulomb potential $v(\mathbf{r}_{12})=1/|\mathbf{r}_{12}|$ and a function $g(\mathbf{r}_{12})$ of the inter-electronic radius $\mathbf{r}_{12}$

[TABLE]

The function $g(\mathbf{r}_{12})$ is called transition pair correlation function (PCF). It can be explicitly given by

[TABLE]

and corresponds to the Fourier transform of the transition structure factor $S(\mathbf{q})$, recently studied in Refs. Liao and Grüneis (2016); Gruber et al. (2018). The indices $i$, $j$ and $a$, $b$ label occupied and virtual spatial orbitals, respectively. In the UEG, orbitals correspond to plane waves with respective wave vectors $\mathbf{k}_{i}$, $\mathbf{k}_{j}$ and $\mathbf{k}_{a}$, $\mathbf{k}_{b}$. The amplitudes $t_{ij}^{ab}$ are obtained by solving the underlying amplitude equations on the respective level of many-electron theory such as second-order perturbation theory Møller and Plesset (1934), coupled cluster singles and doubles (CCSD) theory Coester and Kümmel (1960); Čížek (1969), the random-phase approximation Macke (1950); Pines and Bohm (1952); Gell-Mann and Brueckner (1957) or ladder theory Yasuhara (1972, 1974). In coupled cluster doubles theory, the $t_{ij}^{ab}$’s are obtained by solving the recursive amplitude equations Szabo and Ostlund (1996); Shavitt and Bartlett (2009) that include all terms also present in the ring and the ladder approximation

[TABLE]

where $\Delta^{ij}_{ab}=\varepsilon_{i}+\varepsilon_{j}-\varepsilon_{a}-\varepsilon_{b}$ are the one-electron energy differences in the Hartree–Fock approximation. The amplitude equation above illustrates that coupled cluster methods perform a resummation of certain diagrams including particle-hole ring and particle-particle ladder (ppl) diagrams to infinite order coupling different diagrammatic channels. We note that single particle-hole excitation amplitudes ($t_{i}^{a}$) make no contribution to the wavefunction and correlation energy of the uniform electron gas Shepherd et al. (2014a). Given the amplitudes $t_{ij}^{ab}$ that solve the above equation we can replace $t^{ab}_{ij}$ in Eq. (2) by the right-hand-side of Eq. (3) to arrive at a decomposition of the transition pair correlation function into corresponding diagrammatic contributions

[TABLE]

Analogously, the transition structure factor $S^{\mathrm{CCSD}}(\mathbf{q})$ can be decomposed into its diagrammatic contributions. We note that a similar labeling of terms in the amplitude equations was employed in Refs. Shepherd et al. (2014a, b); Bishop and Lührmann (1978, 1982); Bishop (1991). In contrast to previous work, we perform a decomposition only in one order of perturbation of the otherwise fully coupled amplitudes.

Results. – In the following we will discuss numerical results of the decomposed pair correlation function obtained for the UEG. The results shown in Fig. 1(a) have been computed employing a 54 electron gas simulation cell at a density corresponding to $r_{\mathrm{s}}=5\,a_{0}$ identical to the system used in Ref. Grüneis et al. (2013). The amplitudes are expanded in an orbtial basis composed of 1863 plane waves. We first discuss $g^{\mathrm{driver}}(\mathbf{r}_{12})$, which captures the second-order correlation energy contribution to Eq. (1) and exhibits a minimum at $\mathbf{r}_{12}=\mathbf{0}$. Due to the employed finite one-electron basis set approximation, $g^{\mathrm{driver}}(\mathbf{r}_{12})|_{\mathbf{r}_{12}=\mathbf{0}}$ lacks a derivative discontinuity as required by the cusp condition Kimball (1973). In passing we note that $g^{\mathrm{driver}}(\mathbf{r}_{12})$ of a single electron pair can be well approximated in the complete basis set (CBS) limit at short interelectronic distances using a Slater-type correlation function, $-\gamma^{-1}{\mathrm{e}}^{-\gamma|\mathbf{r}_{12}|}$, where $\gamma$ is a parameter that increases with increasing density Grüneis et al. (2013). In decreasing order of absolute magnitude, the largest contributions to $g^{\mathrm{CCSD}}(\mathbf{0})$ originate from $g^{\mathrm{driver}}(\mathbf{r}_{12})$, $g^{\mathrm{rest}}(\mathbf{r}_{12})$, $g^{\mathrm{ppl}}(\mathbf{r}_{12})$ and $g^{\mathrm{ring}}(\mathbf{r}_{12})$. In contrast to $\mathbf{r}_{12}=\mathbf{0}$, we find that the relative contribution of $g^{\mathrm{ring}}(\mathbf{r}_{12})$ to $g^{\mathrm{CCSD}}(\mathbf{r}_{12})$ increases with increasing interelectronic distance, whereas $g^{\mathrm{ppl}}(\mathbf{r}_{12})$ decays to zero rapidly. This observation demonstrates that the employed decomposition allows for an analysis of the interelectronic correlation strength and its dependence on the distance for individual diagrammatic contributions. We also note that the importance of ring diagrams for the correlation energy at large interelectronic distances is known and of particular significance to the correct description of dispersion interactions in semiconductors using the RPA Harl and Kresse (2009); Dobson et al. (2006).

We now discuss the convergence of the decomposed pair correlation function with respect to the employed orbital basis set. $g^{\mathrm{driver}}(\mathbf{r}_{12})$ converges slowly to the complete basis set limit. This can be seen from $\Delta g^{\mathrm{driver}}(\mathbf{r}_{12})$ as depicted in Fig. 1(b), which corresponds to the difference of $g^{\mathrm{driver}}(\mathbf{r}_{12})$ calculated using 1021 and 1863 plane wave orbitals. This analysis demonstrates that the ppl term converges at a similar rate albeit with an opposite sign as indicated by $\Delta g^{\mathrm{driver}}(\mathbf{r}_{12})$ and $\Delta g^{\mathrm{ppl}}(\mathbf{r}_{12})$ in Fig. 1(b). Moreover, $\Delta g^{\mathrm{driver}}(\mathbf{r}_{12})$ and $\Delta g^{\mathrm{ppl}}(\mathbf{r}_{12})$ are largest in magnitude for short interlectronic distances, illustrating that the increasing basis set is required to capture short-ranged electronic correlation effects. Note that other than $g^{\mathrm{driver}}(\mathbf{r}_{12})$ and $g^{\mathrm{ppl}}(\mathbf{r}_{12})$, the ring diagrams converge rapidly with respect to the employed number of orbitals. Likewise we find that $g^{\mathrm{rest}}(\mathbf{r}_{12})$ converges significantly faster than $g^{\mathrm{driver}}(\mathbf{r}_{12})$ and $g^{\mathrm{ppl}}(\mathbf{r}_{12})$ and we will not perform a more detailed analysis of this term in the present work.

Based on the decomposition of the electron pair correlation function introduced above, we now discuss the corresponding correlation energy contributions and their convergence to the complete basis set limit. To this end we have estimated the CBS limit reference energies by employing 1863 virtual orbitals, which provides sufficiently well converged values for the purpose of the present discussion. We study an identical simulation cell to the one described above with a density corresponding to $r_{\mathrm{s}}=5\,a_{0}$. In addition we discuss results obtained for $r_{\mathrm{s}}=1.8\,a_{0}$. Figures 1(c) and (d) show the convergence of the basis set incompleteness errors for the decomposed correlation energy retrieved as a function of the employed orbital basis set size ($\Delta E^{\mathrm{driver}}$, $\Delta E^{\mathrm{ring}}$, $\Delta E^{\mathrm{ppl}}$ and $\Delta E^{\mathrm{rest}}$). A comparison between Fig. 1(c) and Fig. 1(d) reveals that the convergence behavior of most terms is similar for both densities. Only $\Delta E^{\mathrm{rest}}$ decays to zero in one case from above and in the other case from below. We find that the slowest decay of the basis set incompleteness error is observed for $\Delta E^{\mathrm{driver}}$ and $\Delta E^{\mathrm{ppl}}$, in agreement with the discussion of the convergence of the corresponding pair correlation functions carried out in the paragraph above. $\Delta E^{\mathrm{ring}}$ and $\Delta E^{\mathrm{rest}}$ contributions exhibit a rapid convergence.

To better understand the reason for the convergence behavior of the ring and particle-particle ladder correlation energy contributions with respect to the employed basis, we perform an analysis of a simplified system. To this end we consider two electrons in a box with a homogeneous background charge such that $\mathbf{k}_{i}=\mathbf{k}_{j}=\mathbf{0}$, $\mathbf{k}_{a}=-\mathbf{k}_{b}=\mathbf{q}$ due to momentum conservation and the Coulomb integrals are $v^{ab}_{ij}=4\pi/\mathbf{q}^{2}$. We separate the kinetic and the exchange part of one-electron energy differences $\Delta^{ij}_{ab}=-(\mathbf{q}^{2}+\delta_{q}^{2})$, where the latter is negligible compared to the kinetic term for large magnitudes of $\mathbf{q}$. For the amplitudes $t_{ij}^{ab}$ we use a first-order approximation that reads $t_{ij}^{ab}\approx-4\pi/(\mathbf{q}^{2}+\gamma^{2})^{2}$, which corresponds to the Fourier transform of the Slater-type correlation function $-\gamma^{-1}{\mathrm{e}}^{-\gamma|\mathbf{r}_{12}|}$. Replacing $t_{ij}^{ab}$ in the expressions for $S^{\mathrm{ring}}(\mathbf{q})$ and $S^{\mathrm{ppl}}(\mathbf{q})$ by the above approximation in the limit of a large simulation cell yields

[TABLE]

Note the duality relation of the terms in brackets: one is a point-wise product, the other a convolution of the Coulomb potential with the Slater-type correlation function. The latter turns into a point-wise product in real space $\gamma^{-1}{\mathrm{e}}^{-\gamma|\mathbf{r}_{12}|}/|\mathbf{r}_{12}|$, which is the Yukawa potential scaled by $\gamma^{-1}$ whose Fourier transform is $\gamma^{-1}4\pi/(\mathbf{q}^{2}+\gamma^{2})$. The important consequences of the above analysis are as follows: (i) the asymptotic behavior of the ring transition structure factor is $\lim_{|\mathbf{q}|\rightarrow\infty}S^{\rm ring}(\mathbf{q})\propto 1/\mathbf{q}^{8}$ , while (ii) the ppl transition structure factor decays as $\lim_{|\mathbf{q}|\rightarrow\infty}S^{\rm ppl}(\mathbf{q})\propto 1/\mathbf{q}^{4}$ only, which is identical to $\lim_{|\mathbf{q}|\rightarrow\infty}S^{\mathrm{driver}}(\mathbf{q})\propto 1/\mathbf{q}^{4}$. This explains qualitatively the observed correlation energy convergence in Fig. 1 for the ring and ppl term with respect to the employed plane wave basis. Moreover, the above findings indicate that the ppl contribution becomes more important with decreasing electron density because $E^{\mathrm{ppl}}\propto\gamma^{-1}$ and $\gamma\propto 1/r_{\mathrm{s}}$ Grüneis (2015). The above analysis also shows that the ppl contribution modifies the first-order coefficient in the Taylor expansion of $g(|\mathbf{r}_{12}|)$ with respect to $|\mathbf{r}_{12}|$ around $|\mathbf{r}_{12}|=0$ from above.

We now propose a finite basis set correction method for coupled cluster theory that is based on the analysis of the ppl term discussed above. Coupled cluster methods are becoming increasingly popular to perform ab-initio studies of solids and surfaces with high accuracy Gruber and Grüneis (2018); Tsatsoulis et al. (2018, 2017); McClain et al. (2017); Kubas et al. (2016); Boese and Sauer (2016); Yang et al. (2014); Booth et al. (2013), demonstrating significant potential to further expand the scope of electronic structure theory calculations in; for example, the field of surface chemistry or the study of thermodynamic stabilities of solids. However, the underlying computational cost is much larger than that of the more efficient yet less accurate approximate density functional theory calculations. The main source of the computational cost in coupled cluster calculations originates from the ppl term. Its computational complexity scales as $\mathcal{O}(N_{\mathrm{v}}^{4}N_{\mathrm{o}}^{2})$ in a canonical formulation, where $N_{\mathrm{v}}$ and $N_{\mathrm{o}}$ are the number of virtual and occupied orbitals, respectively. However, Fig. 1 and the discussion above show that the basis set incompleteness errors of the second-order correlation and the ppl term are propotional to each other albeit with an opposite sign. This observation motivates the following approximation to the CCSD correlation energy in the CBS limit that aims at accounting for the basis set incompleteness error in the driver ($E^{\mathrm{driver}}$) and ppl ($E^{\mathrm{ppl}}$) term explicitly:

[TABLE]

where

[TABLE]

We refer to the employed orbital basis set size by $N$. $E^{\mathrm{driver}}(\mathrm{CBS})$ can in practice be estimated in a computational efficient manner using basis set extrapolation or explicit correlation techniques Grüneis (2015); Shepherd et al. (2012); Schäfer et al. (2017).

We now assess the efficiency of the proposed finite basis set correction to the coupled cluster correlation energy for ab-initio systems including solids and atoms. All periodic calculations of solids have been performed using the Vienna ab-initio simulation package (VASP) Kresse and Furthmüller (1996) in the framework of the projector augmented wave method Blöchl (1994), interfaced to our coupled cluster code Hummel et al. (2017) that employs an automated tensor contraction framework (CTF) Solomonik et al. (2014). We use natural orbitals to achieve a compact approximation for the virtual orbital space Grüneis et al. (2011). Fig. 2 depicts the coupled cluster singles and doubles correlation energy convergence for a range of periodic crystals: C (diamond), BN (wurtzite), LiF (rock-salt) and Si (diamond), calculated with and without the proposed finite basis set correction. Our findings unequivocally demonstrate that the basis set convergence of the coupled cluster singles and doubles correlation energy is significantly faster including the proposed correction compared to uncorrected CCSD energies. Already 10 natural orbitals per occupied orbital agree with the CBS limit results to within chemical accuracy ($\approx$40 meV/atom) for the total energy per atom as indicated by the region between the dashed horizontal lines. 80 natural orbitals per occupied orbital serve as the approximate CBS reference in the present work. In contrast to CCSD-PPL, CCSD requires approximately 30 natural orbitals per occupied orbital to achieve a similar level of precision. We stress that a reduction in the number of virtual orbitals by a factor three, reduces the computational cost for the ppl term by two orders of magnitude.

To further assess CCSD-PPL, we have carried out calculations for the Ne atom employing a modified version of an open-source quantum chemistry code (PSI4 Parrish et al. (2017)) and atom-centered Gaussian basis sets. Table 1 shows that CCSD-PPL converges rapidly to the CBS limit. Employing a double-$\zeta$ (DZ) basis set yields CCSD-PPL correlation energies that are closer to the CBS limit than CCSD correlation energies obtained using a quadruple-$\zeta$ (QZ) basis set. Therefore the conclusions drawn from the uniform electron gas also hold in the atomic limit. We refer the reader to Ref. Irmler and Grüneis (2019) for a more extensive investigation of CCSD-PPL in atoms and molecules and a comparison to explicitly correlated methods.

Concluding Remarks. – We have presented a diagrammatic decomposition of the coupled cluster pair correlation function and performed an analysis of the ring and particle-particle ladder terms. Based on their dual structure and the observed convergence with respect to the employed basis, we have introduced an efficient finite basis set correction that allows for reducing the computational cost of coupled cluster theory calculation in atoms, molecules and solids substantially. In combination with the recently proposed finite size corrections Gruber et al. (2018) and other techniques Grüneis (2015); Hummel et al. (2017); Booth et al. (2016), this paves the way for a routine use of coupled cluster theory in electronic structure calculations of solids and surfaces. However, the significance of the present work extends far beyond the proposed basis set correction. Given exact reference results for the spin-resolved pair correlation function in ab-initio or model systems that can be obtained using; for example, full configuration interaction quantum Monte Carlo Booth et al. (2013), the errors in the diagrammatic channels of approximate many-electron theories can be analysed in detail, allowing for developing more accurate and better balanced truncations to the many body perturbation expansion of the exact electronic correlation energy. Finally, we note that the introduced decomposition scheme can also help developing more efficient embedding and local correlation theories, where the decoupling and decay of interelectronic correlation effects in different diagrammatic channels plays a crucial role and determines the accuracy of the employed approximations.

Acknowledgements. – The authors thankfully acknowledge support and funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (Grant Agreement No 715594). A. G. gratefully acknowledges many helpful discussions on F12 theory with Seiichiro Ten-no and David Tew.