An optimized twist angle to find the twist-averaged correlation energy applied to the uniform electron gas

TL;DR

This paper proposes a single, optimized twist angle for coupled cluster doubles calculations that approximates twist averaging, significantly reducing computational cost while maintaining accuracy in the thermodynamic limit for the uniform electron gas.

Contribution

It introduces the concept of connectivity to identify a special twist angle that mimics full twist averaging in CCD calculations, improving efficiency.

Findings

The special twist angle yields results comparable to full twist averaging.

The method reduces computational cost from multiple CCD calculations to one CCD plus several MP2 calculations.

It enhances the accuracy of finite electron number calculations in representing the thermodynamic limit.

Abstract

We explore an alternative to twist averaging in order to obtain more cost-effective and accurate extrapolations to the thermodynamic limit (TDL) for coupled cluster doubles (CCD) calculations. We seek a single twist angle to perform calculations at, instead of integrating over many random points or a grid. We introduce the concept of connectivity, a quantity derived from the non-zero four-index integrals in an MP2 calculation. This allows us to find a special twist angle that provides appropriate connectivity in the energy equation, and which yields results comparable to full twist averaging. This special twist angle effectively makes the finite electron number CCD calculation represent the TDL more accurately, reducing the cost of twist-averaged CCD over $N_\mathrm{s}$ twist angles from $N_s$ CCD calculations to $N_s$ MP2 calculations plus one CCD calculation.

Equations

$$E_{\mathrm{HF}}({\mathbf{k}}_s)=\sum _i{T}_i({\mathbf{k}}_s)-\frac{1}{2}\sum _{ij}{v}_{ijji}({\mathbf{k}}_s)$$

$$⟨E_{\mathrm{HF}}{⟩}_{{\mathbf{k}}_{\mathbf{s}}}=\frac{1}{N_{\mathrm{s}}}\sum _{{\mathbf{k}}_s}^{N_{\mathrm{s}}}\sum _i{T}_i({\mathbf{k}}_s)-\frac{1}{N_{\mathrm{s}}}\sum _{{\mathbf{k}}_s}^{N_{\mathrm{s}}}\frac{1}{2}\sum _{ij}{v}_{ijji}({\mathbf{k}}_s)$$

$$⟨E_{\mathrm{HF}}{⟩}_{{\mathbf{k}}_{\mathbf{s}}}=\sum _i\left[\frac{1}{N_{\mathrm{s}}}\sum _{{\mathbf{k}}_s}^{N_{\mathrm{s}}}{T}_i({\mathbf{k}}_s)\right]-\frac{1}{2}\sum _{ij}\left[\frac{1}{N_{\mathrm{s}}}\sum _{{\mathbf{k}}_s}^{N_{\mathrm{s}}}{v}_{ijji}({\mathbf{k}}_s)\right].$$

$$⟨E_{\mathrm{HF}}{⟩}_{{\mathbf{k}}_{\mathbf{s}}}=\sum _i⟨T_i{⟩}_{{\mathbf{k}}_{\mathbf{s}}}-\frac{1}{2}\sum _{ij}⟨v_{ijji}{⟩}_{{\mathbf{k}}_{\mathbf{s}}}$$

$$⟨E_{\mathrm{corr}}{⟩}_{{\mathbf{k}}_{\mathbf{s}}}=\frac{1}{N_{\mathrm{s}}}\sum _{{\mathbf{k}}_s}^{N_{\mathrm{s}}}\frac{1}{4}\sum _{ijab}{\bar{t}}_{ijab}({\mathbf{k}}_s){\bar{v}}_{ijab}({\mathbf{k}}_s),$$

$${\bar{t}}_{ijab}({\mathbf{k}}_s){\bar{v}}_{ijab}({\mathbf{k}}_s)=\frac{|{\bar{v}}_{ijab}({\mathbf{k}}_s){|}^2}{{ϵ}_i({\mathbf{k}}_s)+{ϵ}_j({\mathbf{k}}_s)-{ϵ}_a({\mathbf{k}}_s)-{ϵ}_b({\mathbf{k}}_s)}$$

$$v_{ijab}=\frac{4\pi }{L^3}\frac{1}{({\mathbf{k}}_i-{\mathbf{k}}_a{)}^2}{\delta }_{{\mathbf{k}}_i-{\mathbf{k}}_a,{\mathbf{k}}_b-{\mathbf{k}}_j}{\delta }_{{\sigma }_i{\sigma }_a}{\delta }_{{\sigma }_j{\sigma }_b}.$$

$$\underset{{\mathbf{k}}_{\mathbf{s}}}{\min }\sum _x\frac{1}{x^2}{\left(h_x-⟨h_x{⟩}_{{\mathbf{k}}_{\mathbf{s}}}\right)}^2$$

Full text

I Introduction

In recent years, the use of wavefunction-based post-Kohn–Sham or post-Hartree–Fock methods to solve problems in materials science has proliferated. Müller and Paulus (2012) This is in part driven by an interest in obtaining precise energies (accurate to within 1mHa) for complex systems using hierarchies of methods found in quantum chemistry such as coupled cluster theory. While growing in popularity, wavefunction methods have yet to see widespread adoption, in large part due to their significant computational cost scaling with system size. This is especially of note in coupled cluster theory using a plane wave basis, and as a result, some authors are seeking methods to control finite size errors in order to run calculations using smaller system sizes. Gruber et al. (2018)

Finite size errors arise when attempts are made to simulate an infinite system Hamiltonian with a periodic supercell containing a necessarily finite particle number.Fraser et al. (1996); Drummond et al. (2008) The finite size of a supercell places a limitation on the minimum momenta in Fourier sums (e.g., with a cubic box of length $L$, the smallest momentum transfer is $2\pi/L$). These limitations ultimately lead to errors in the correlation energy; Gruber et al. (2018); Ruggeri, Ríos, and Alavi (2018) this has been attributed to long range van der Waals forces. Gruber et al. (2018); Gruber and Grüneis (2018)

Since these finite size errors are large and slowly converging with increasing supercell size, which has been analyzed in detail for coupled cluster theory, McClain et al. (2017) there has been significant interest in developing wavefunction methods with reduced computational cost to circumvent finite size error and allow the treatment of larger supercells. These include embedding methods,Sun and Chan (2016) such as density matrix embedding,Knizia and Chan (2012, 2013); Bulik, Chen, and Scuseria (2014); Bulik, Scuseria, and Dukelsky (2014); Ricke et al. (2017); Zheng et al. (2017); Pham, Bernales, and Gagliardi (2018) wavefunction-in-DFT embedding,Henderson (2006); Tuma and Sauer (2006); Gomes, Jacob, and Visscher (2008); Sharifzadeh, Huang, and Carter (2009); Huang, Pavone, and Carter (2011); Libisch, Huang, and Carter (2014); Manby et al. (2012); Goodpaster et al. (2014); Chulhai and Goodpaster (2018) electrostatic embedding,Hirata et al. (2005); Dahlke and Truhlar (2007); Hirata (2008); Leverentz and Truhlar (2009); Bygrave, Allan, and Manby (2012) QM/MM-inspired schemes,Shoemaker and Gordon (1999); Sherwood et al. (2003); Herschend, Baudin, and Hermansson (2004); Beran and Nanda (2010); Chung et al. (2015) and others.Eskridge, Krakauer, and Zhang (2018); Lan, Kananenka, and Zgid (2015); Rusakov et al. (2019); Voloshina, Gaston, and Paulus (2007); Masur et al. (2016) Local correlation methodsCollins and Bettens (2015); Usvyat, Maschio, and SchÃŒtz (2018) such as fragment-based schemes,Gordon et al. (2012); Li, Li, and Jiang (2007); Li, Ni, and Li (2016); Rolik and Kállay (2011); Li and Li (2004); Kobayashi, Imamura, and Nakai (2007); Kristensen et al. (2011); Ghosh et al. (2010); Kitaura et al. (1999); Fedorov and Kitaura (2007); Netzloff and Collins (2007); Ziółkowski et al. (2010) incremental methods,Stoll (1992); Paulus (2006); Friedrich, Hanrath, and Dolg (2007); Stoll and Doll (2012); Friedrich and Hänchen (2013); Voloshina and Paulus (2014); Kállay (2015); Fertitta et al. (2018) and heirarchical methods,Deev and Collins (2005); Manby, Alfè, and Gillan (2006); Nolan et al. (2009); Collins (2011) break the system into smaller subsystems, then extrapolate or stitch together the energies. Some methods take advantage of range separationToulouse et al. (2009); Bruneval (2012); Shepherd, Henderson, and Scuseria (2014a) or other distance-based schemesSpencer and Alavi (2008); Maurer et al. (2013); Kats and Manby (2013); Kats (2016); Ayala, Scuseria, and Savin (1999) to reduce computational cost.

In addition to work on developing or modifying electronic structure methods, much work on reducing the cost of wavefunction methods has been focused on modifying basis sets in order to accelerate convergence and decrease computation time. Local orbital methods have been popular,Pisani et al. (2005); Ayala, Kudin, and Scuseria (2001); Usvyat, Maschio, and SchÃŒtz (2015); Werner, Manby, and Knowles (2003); Flocke and Bartlett (2004); Werner et al. (2015); Rolik et al. (2013); Förner et al. (1985); SchÃŒtz and Werner (2000); Neese, Wennmohs, and Hansen (2009); Sun et al. (2017); Booth et al. (2016); Blum et al. (2009); Subotnik and Head-Gordon (2005) often based on the local ansatz of Pulay and SaeboSaebo and Pulay (1993) or Stollhoff and Fulde.Stollhoff and Fulde (1977) Other common methods include progressive downsampling,Shimazaki and Hirata (2009); Hirata and Shimazaki (2009); Ohnishi and Hirata (2010) downfolding,Purwanto, Zhang, and Krakauer (2013) use of explicitly-correlated basis setsAdler and Werner (2009); Shiozaki and Hirata (2010); Grüneis et al. (2013); Usvyat (2013); Grüneis (2015) or natural orbitals,Grüneis et al. (2011) and tensor manipulations.Hohenstein, Parrish, and Martínez (2012); Benedikt, Böhm, and Auer (2013); Hummel, Tsatsoulis, and GrÃŒneis (2017); Peng and Kowalski (2017); Motta et al. (2018) Discussion of the details and relative merits of these methods is beyond the scope of this paper; for a review, we direct the interested reader to Refs. Huang and Carter, 2008; Müller and Paulus, 2012; Beran, 2016; GrÃŒneis, 2018.

However, there has been some work on developing corrections for finite size errors.Fraser et al. (1996); Kent et al. (1999); Kwee, Zhang, and Krakauer (2008); Drummond et al. (2008); Holzmann et al. (2016); Liao and Grüneis (2016) Many-body methods can sometimes be integrated to the thermodynamic limit (TDL),Gell-Mann and Brueckner (1957); Nozières and Pines (1958); Onsager, Mittag, and Stephen (1966); Bishop and Lührmann (1982, 1978); Bishop (1991); Ziesche (2007) allowing for the derivation of analytic finite-size correction expressions.Chiesa et al. (2006) Several studies from the last year have particular relevance to our work here. Grüeneis *et al.*Gruber et al. (2018); Gruber and Grüneis (2018) employed a grid integration within periodic coupled cluster for ab initio Hamiltonians with applications to various solids. In another study, Alavi *et al.*Ruggeri, Ríos, and Alavi (2018) devised a novel extrapolation relationship that links different electron gas calculations through the density parameter. Both of these papers use a technique known as twist averaging to try to remove finite size error.

Twist averaging is a method that attempts to control finite size errors by first offsetting the $k$-point grid by a small amount, ${\bf k}_{s}$, and then averaging over all possible offsets.Lin, Zong, and Ceperley (2001)

We refer to ${\bf k}_{s}$ here as a twist angle. One of the main purposes of twist averaging is to provide for a smoother extrapolation to the thermodynamic limit by reducing severe energy fluctuations as the particle number varies.

When performed with a fixed particle number and box length, this process is referred to as twist averaging in the canonical ensemble, which is what we study here. When employed in stochastic methods, such as variational Monte Carlo,Lin, Zong, and Ceperley (2001) diffusion Monte CarloDrummond et al. (2008) or full configuration interaction quantum Monte Carlo,Ruggeri, Ríos, and Alavi (2018); Shepherd (2013) the grid can be stochastically sampled at the same time as the main stochastic algorithm, and both stochastic error and error in twist-averaging related to approximate integration can be removed at the same time. As a result, the scaling with the number of twist angles sampled is extremely modest. Unfortunately, the same cost savings cannot be realized for deterministic methods. In this case, in order to achieve a reasonable estimate for the average, one must use a large number of individual energy calculations. This results in the cost scaling linearly with the number of twist angles used, although the lessening of finite size effects with rising electron number would alleviate this scaling to some extent.McClain et al. (2017)

Here, we seek to remedy the linear scaling of twist averaging for deterministic methods by devising a way to provide an energy that is as accurate as twist-averaging, but with single-calculation cost.

In principle, it is possible to find a single twist angle which exactly reproduces the total twist-averaged energy by recognizing that it is an integral of the energy over the twist angles for a system. This was the same logic used in analysis by Baldereschi to find a special $k$-pointBaldereschi (1973) and has been used by others in the QMC community to find a special twist angle.Dagrada et al. (2016); Rajagopal et al. (1994, 1995) We are motivated similarly and wish to find a single twist angle that yields an energy approximately equal to the full twist-averaged energy for CCD and related wavefunction methods.

We take advantage of the similarity between the MP2 and CCD correlation energy expressions, using the much cheaper MP2 method to find a single twist angle that produces a system with the most similar number of allowed excitations to the twist-averaged system. We refer to this set of allowed excitations as the ‘connectivity’. We then use this twist angle to calculate the CCD energy, which is in good agreement with the fully twist-averaged CCD energy. Finally, we compare our energies to those obtained using one twist angle at the Baldereschi point.Baldereschi (1973) We do not seek to completely remedy the whole of the finite size error, instead noting that other authors have come up with corrections or extrapolations that can be used after twist-averaging is applied.Chiesa et al. (2006); Drummond and Needs (2009); Gruber and Grüneis (2018)

II Twist averaging & Connectivity

Both continuum/real-space and basis-set twist averaging have been used effectively in quantum Monte Carlo calculations;Lin, Zong, and Ceperley (2001); Drummond et al. (2008); Ruggeri, Ríos, and Alavi (2018); Shepherd (2013) however, twist averaging remains relatively rare in coupled cluster calculations. In Fig. 1, the total $\Gamma$-point CCD energy ($N=38$ to $N=922$) and twist-averaged CCD energy ($N=38$ to $N=294$) are plotted alongside the extrapolation to the TDL for the uniform electron gas ($0.609(3)$ Ha/electron, where the error in the last digit is in parentheses). The CCD calculation is performed in a finite basis that is analogous to a minimal basis.Shepherd and Grüneis (2013) The $\Gamma$-point energy is highly non-monotonic; it does not fit well with the extrapolation. The twist-averaged data shows a much better fit with the extrapolation, resulting in a better estimate of the TDL. The drawback of twist averaging, however, is that it costs $N_{\mathrm{s}}\,\mathcal{O}\mathrm{[CCD]}$ for $N_{\mathrm{s}}$ twist angles (here, 100). The twist-averaged energy becomes too costly to calculate with CCD for system sizes above 294 electrons.

Figure 1 is a clear statement of the problem we wish to resolve here. Twist averaging resolves some finite size errors that are present at an individual particle number $N$, and allows for improved extrapolation to the thermodynamic limit. That said, the scaling with the number of twist angles is cost-prohibitive. We aim to develop an approximation to twist averaging that gives comparable accuracy at a fraction of the cost. We begin by analyzing how the Hartree-Fock energy and the MP2 correlation energy are modified by twist averaging. This analysis then allows us to build an algorithm that produces CCD twist-averaged accuracy/results for only MP2 cost.

II.1 Hartree-Fock and single-particle eigenvalues

A finite-sized electron gas at the $\Gamma$-point is only closed-shell at certain so-called magic numbers, which are determined by the symmetry of the lattice (for example $N=2$, 14, 38, and 54). One of the reasons that the $\Gamma$-point calculations are so noisy (Fig. 1) is that there are degeneracies in the HF eigenvalues, which can be seen in Fig. 2 and has long been recognized. Drummond and Needs (2009) This can be partially remedied by modifying the Hartree–Fock eigenvalues. The starting-point for this is writing the HF energy as follows:

[TABLE]

where $T_{i}$ is the kinetic energy of orbital $i$ and $v_{ijji}$ is the exchange integral between electrons in orbitals $i$ and $j$. Here, we have included the explicit form of the dependence on the twist angle, ${\bf k}_{s}$.

The twist-averaged energy is found by summing Eq. (1) over all possible ${\bf k}_{s}$:

[TABLE]

where $N_{s}$ indicates the number of twist angles used. Swapping the sums yields:

[TABLE]

Therefore, twist averaging the HF energy is numerically identical to twist averaging the individual matrix elements:

[TABLE]

Overall, then, we can use twist-averaged HF eigenvalues in place of twist-averaging the HF energy, obtaining a more reasonable density of states Fig. 2. We will use this in our subsequent scheme.

II.2 Beyond Hartree–Fock

The above approach does not generalize to correlated theories because they have more complex energy expressions. For example, averaging the second-order Møller-Plesset theory (MP2) correlation energy over all possible twist angles can be written:

[TABLE]

where $i$ and $j$ refer to occupied orbitals and $a$ and $b$ refer to unoccupied orbitals. The symbols $\bar{v}$ and $\bar{t}$ refer to the antisymmetrized electron-repulsion integral and amplitude respectively. For MP2:

[TABLE]

Even though MP2 diverges in the thermodynamic limit, the energy expression (Eq. (5)) has a similar structure to coupled cluster theory, the random phase approximation, and even full configuration interaction quantum Monte Carlo. As such, we can make generalized observations using the MP2 energy expression, and then use these observations to derive a scheme to find an optimal $k_{s}$ twist angle that works for all of these methods.

II.3 The connectivity approach

The MP2 correlation energy can vary substantially as the twist angle is changed. For example, in the $N=14$ electron system with a basis set of $M=38$ orbitals, the MP2 energy can vary between $-0.0171$ Ha/electron to $-0.0001$ Ha/electron. This arises, in particular, because the number of low-momentum excitations (minimum $|{\bf k}_{i}-{\bf k}_{a}|$) will vary significantly. Since the contribution of each excitation to the MP2 sum is $|{\bf k}_{i}-{\bf k}_{a}|^{-4}$, there is a rapid decay of an excitation’s contribution to the correlation energy beyond the minimum vector.

This effect arises because, when the twist angle is changed, different orbitals now fall into the occupied ($ij$) space, and different orbitals fall into the virtual ($ab$) space. This changes the value of the sum over both occupied and virtual orbitals, since many individual terms in the sum are now substantively different. We illustrate this using a diagram in the Supplementary Information.

By contrast, the integrals themselves do not change; to show this, the integral can be written:

[TABLE]

The Kronecker deltas, $\delta$, ensure that momentum and spin symmetry (denoted $\sigma$) are conserved. On changing ${\bf k}_{p}\rightarrow{\bf k}_{p}+{\bf k}_{s}$ for all ${\bf k}$’s, the difference in the denominator here does not change, since $({\bf k}_{i}+{\bf k}_{s}-{\bf k}_{a}-{\bf k}_{s})^{2}=({\bf k}_{i}-{\bf k}_{a})^{2}$. In general, our calculations were set up using details which can be found in our prior work e.g. Ref. Shepherd et al., 2012.

At this stage, we conjecture that if one of the mechanisms by which twist averaging is affecting the MP2 energy (and other correlation energies) is to smooth out the inconsistent contributions between different momenta, then it might be possible for us to find a ‘special twist angle’ where the number of low-momentum states for that single twist angle is a good match to the average number of momentum states across all twist angles. Further, we will show this special twist angle is transferable to other, more sophisticated methods such as coupled cluster doubles theory.

To find this special twist angle, we proceed as follows:

For a given twist angle ${\bf k}_{s}$, loop over the same $ijab$ as the MP2 sum $\sum_{ijab}$. For each $ijab$ set:

1. (a)

Determine the momentum transfer $x=|{\bf n}_{i}-{\bf n}_{a}|^{2}$ where ${\bf n}_{a}$ is the integer equivalent of the quantum number: ${\bf k}_{a}=\frac{2\pi}{L}{\bf n}_{a}$. 2. (b)

Increment a histogram element $h_{x}$ by one. 2. 2.

Create a vector ${\bf h}$, whose elements are $h_{x}$, which correspond to the number of of $v_{ijab}$ matrix elements with magnitude $\frac{1}{\pi L}\frac{1}{x}$ that are encountered during the MP2 sum. 3. 3.

Average ${\bf h}$ over all twist angles, yielding $\langle{\bf h}\rangle_{\bf{k}_{s}}$ 4. 4.

Loop over the twist angles again, and find the single ${\bf h}$ (and corresponding twist angle) that best matches $\langle{\bf h}\rangle_{\bf{k}_{s}}$ using:

[TABLE]

The weight term $1/x^{2}$ was chosen empirically to diminish the contributions of large numbers of high-momentum weights that contribute relatively little to the energy.

Looking at Eq. (5), there are two ways to proceed. We could either use this special ${\bf k}_{s}$ for all aspects of the calculation (e.g. for both the integral evaluation and the eigenvalue difference), or we could use the special ${\bf k}_{s}$ for the integral only, and twist-average the eigenvalues before performing the CCD calculation. We found that the latter was more numerically effective for $N=14$ and decided to use this approach to generate the results presented here. In general, though, for larger systems it does not make a large difference.

In practice, we implemented this algorithm within an MP2 and CCD code; we call the MP2 calculation at each twist angle and then the CCD calculation once at the end. For the remainder of this work, we will call this application of the above algorithm the “connectivity scheme,” referencing the idea that the pattern of non-zero matrix elements $v_{ijab}$ resembles a connected network.

III Results

We demonstrate the effectiveness of this algorithm for coupled cluster calculations on the uniform electron gas in Fig. 3. In general, our results as show that the connectivity scheme works for different electron numbers, basis sets, and $r_{s}$ values. Furthermore, evaluation of the connectivity scheme is approximately 100x cheaper than twist averaging.

In Fig. 3(a), we compare the connectivity scheme to full twist-averaging for CCD calculations on the uniform electron gas. Energy differences from the $\Gamma$-point energy are plotted for each electron number. Our results show that the connectivity scheme delivers comparable accuracy (mean absolute deviation = 0.3 mHa/electron) to twist averaging, with the benefit of being much faster to compute. The connectivity scheme is substantially cheaper than the twist-averaging scheme: the $N=294$ twist-averaged calculation, for example, costs $58$ hours, which is about the same time it takes to run the $N=922$ connectivity scheme calculation. A complete set of timings is provided in the Supplementary Information.

In Fig. 3(b), we compare our connectivity scheme to full twist-averaging over a range of basis set sizes ($M=36-2838$ orbitals) for 54 electrons. In Fig. 3(c), we compare the connectivity scheme to full twist-averaging over a range of $r_{s}$ values ($0.01-50.0$ a.u.) for 54 electrons. In both cases there is good agreement between the two methods for all system sizes, proving that the connectivity scheme delivers good accuracy when compared with twist averaging for a range of both basis set sizes (mean absolute deviation $<$ 0.35 mHa/electron) and $r_{s}$ values (mean absolute deviation $<$ 0.25 mHa/electron) at a decreased cost.

In Fig. 4, we show the extrapolation of our connectivity scheme CCD correlation energy to the thermodynamic limit for the $r_{s}=1.0$ uniform electron gas. We perform calculations up to $N=922$ electrons, and fit these results to the equation $E=a+bN^{-1}$, as proposed by other authors. Drummond et al. (2008) We then use this fit to extrapolate the correlation energy to the thermodynamic limit. We also performed the same extrapolation for the twist-averaged data set up to $N=294$ electrons (not shown). The extrapolations predict the TDL energy to be $-0.0340(8)$ Ha/electron for the connectivity scheme and $-0.033(4)$ Ha/electron for the twist-averaged scheme, a difference of $0.001(4)$ Ha/electron. The numbers in parentheses are errors in the final digit. These agree within error, and the connectivity scheme has an improved error due to having more data points.

Next, we demonstrate how to use this method to obtain a complete basis set and thermodynamic limit estimate for the uniform electron gas. Connectivity scheme CCD energies were collected for the $N=54$ electron system with basis sets varying from $M=922$ to $M=2838$ orbitals, and for systems with electron numbers varying between $N=162$ to $610$, with $M\approx 4N$. These data allow us to extrapolate to both the complete basis set limit and the thermodynamic limit by using the numerical approach set out in our previous work.Shepherd (2016) This yields an energy that is 0.0566(6), with the error in parentheses resulting from the extrapolations; this is in good agreement with our prior estimate with significantly less error.Shepherd (2016) For more details the reader is referred to the Supplementary Information.

Finally, in Fig. 5, we compare the CCD energies from full twist-averaging, our connectivity scheme, and performing single calculation using the Baldereschi point as a twist angle. This point, first developed for insulators, is well known for the role it played in developing efficient thermodynamic integrationsBaldereschi (1973); Chadi and Cohen (1973); Cunningham (1974); Monkhorst and Pack (1976) and was subsequently used for twist-averaging as the center-point of uniform grid twist-averaging by Drummond *et al.*Drummond and Needs (2009). At higher electron numbers ($N\geq 162$) the difference between BP and the TA energies falls below 1mHa/electron as all of the approaches converge to the same energy. At small electron numbers, however, the Baldereschi point significantly deviates from the twist-averaged energy, while the connectivity scheme is a much better approximation.

IV Discussion & concluding remarks

Our results show that a finite electron gas is best able to reproduce the twist-averaged total and correlation energies when a special $\bf{k}_{s}$-point is chosen to minimize the differences between the momentum connectivity of the finite system and a reference (here, a twist-averaged finite system). Our interpretation of the connectivity-derived special $\bf{k}_{s}$-point’s utility is that the low-momentum two-particle excitations from HF often suffer from finite size errors due to the shape of the Fermi surface in $k$-space. By finding a particularly representative $k_{s}$-point, we aim to take the ‘best case’ of a representative shape–or, at least, as best as can be managed by a truly finite system. When we examine the occupied orbitals in $k$-space at the special $\bf{k}_{s}$-point, they adopt low-symmetry patterns that tend more toward the shape of a sphere than the $\Gamma$-point distribution.

Though we have made significant progress here towards ameliorating finite size error, there are still two open questions. First, could our method be modified in order to minimize the energy difference to the thermodynamic limit rather than just to the twist-averaged energy? The second open question surrounds the extrapolation – in particular, what is the actual form of the energy as the system size tends to infinity? We could investigate this source of error by comparing with the known high-density limit of RPA, which CCD is expected to be able to capture. We leave both of these investigations for future work.

Overall, the results here should improve our ability to understand infinite-sized model systems that are necessarily represented as finite systems, such as the electron gas with varying dimensions, the Hubbard model, and the models of nuclear matter we previously studied. Shepherd (2016); Baardsen et al. (2016); Baardsen (2014) This communication is timely due to a resurgence of interest in the uniform electron gas Neufeld and Thom (2017); White and Chan (2018); Spencer et al. (2018); McClain et al. (2016); Spencer et al. (2019); Malone et al. (2016); Shepherd and Grüneis (2013); Shepherd, Henderson, and Scuseria (2014a); Grüneis et al. (2013) and of twist-averaged coupled cluster calculations. Gruber et al. (2018); Hagen et al. (2014) We expect this work can immediately be applied to improve calculations.

Our long-term goals are to use this approach to study realistic systems. Though calculations are left for future manuscripts, we expect to follow a similar approach to our prior work in this area. In particular, we start by observing the similarity between how twist-averaging works in plane wave ab initio calculations where the energy is still obtained as a sum over matrix elements $v_{ijab}$ (as in Eq. (7)) which are offset by a twist angle. Specifically, then, it should be possible to choose the twist angle in the same way as we propose here, so for a cubic system with $N$ electrons and a box length of $L$, the same twist angle as used here should work. As such, we will soon be applying this to real solids and leave this for a future study.

Supplementary Material.–** ** The reader is directed to the supplementary material for raw data tables and illustrations mentioned in the text.

Acknowledgements.–** ** JJS and TM acknowledge the University of Iowa for funding. JJS thanks the University of Iowa for an Old Gold Award. ARM was supported by the National Science Foundation Graduate Research Fellowship under Grant No. 1122374. The code used throughout this work is a locally modified version of a github repository used in previous work Shepherd, Henderson, and Scuseria (2014a, b): https://github.com/jamesjshepherd/uegccd.