A shortcut to the thermodynamic limit for quantum many-body calculations of metals

Computationally efficient and accurate quantum mechanical approximations to solve the many-electron Schrödinger equation are crucial for computational materials science. Methods such as coupled cluster theory show potential for widespread adoption if computational cost bottlenecks can be removed. For example, extremely dense k-point grids are required to model long-range electronic correlation effects, particularly for metals. Although these grids can be made more effective by averaging calculations over an offset (or twist angle), the resultant cost in time for coupled cluster theory is prohibitive. We show here that a single special twist angle can be found using the transition structure factor, which provides the same benefit as twist averaging with one or two orders of magnitude reduction in computational time. We demonstrate that this not only works for metal systems but also is applicable to a broader range of materials, including insulators and semiconductors.

Many-electron perturbation theories are at the heart of theoretical quantum chemistry and play a central role in understanding and predicting chemical processes in the gas phase, condensed matter systems, and on surfaces.Although finite-order perturbation theories were proposed initially [1], it has become evident that renormalized perturbation theories-resummations of certain diagrams in the many-electron perturbation series expansion-often exhibit a better trade-off between accuracy and computational cost, and are numerically more stable [2][3][4][5][6].The resulting theories have substantially contributed to the success of modern ab initio methods.In particular, these include the coupled cluster method and its different approximations, such as the widelyused coupled cluster singles and doubles plus perturbative triples theory (CCSD(T)) [7].Yet despite the great advances of modern coupled cluster theory approximations and their applications to molecules, surfaces, and solids it remains difficult to treat metallic systems at the same level of theory, with the exception of applications to simple model Hamiltonians such as the uniform electron gas [8][9][10][11][12][13][14][15][16][17][18].This constitutes a significant obstacle in modern computational materials science, where accurate benchmark results obtained with high-level theories for all types of materials are desperately needed to complement computationally more-efficient but less-accurate density functional theory calculations.The important role of metals in surface chemistry, catalysis, and many technological applications highlights the necessity to expand the scope of accurate quantum chemical perturbation theories further.Therefore we are confident that the * james-shepherd@uiowa.edu techniques outlined in this work, which make it possible to study metals using coupled cluster theory routinely, present a significant advance especially in combination with other methodological improvements such as embedding theories [19][20][21].
Despite recent successes of methods development [22][23][24][25][26][27][28][29][30], metals have still remained a formidable challenge to quantum chemistry.The quintessential property of a metal-that it has a zero gap-has led to the view of these calculations being intractable.This view comes from experiences in molecular systems with small HOMO-LUMO gaps, which are prone to strong correlation effects that are inaccessible by single reference methods like coupled cluster theories.Further, small gaps can lead to numerical instability in both Hartree-Fock and resummation methods just in the regime where one approaches the thermodynamic limit [24,[31][32][33].However, for the uniform electron gas, which is the simplest and most popular metallic model Hamiltonian, each of these barriers can be overcome.[9,[13][14][15][16][17][18] Through scrupulous accounting for error and its cancellation, we have found that it is possible to converge each component of the total energy to reach benchmark accuracy.On the journey to a universal application of single-reference coupled cluster theory to all solids, the final hurdle is to transfer these successes to real metals.
Here, we present a simple but efficient scheme to determine an effective Hamiltonian that directly enables the application of widely-used coupled cluster methods to study real metals.The proposed method determines a crystal momentum vector used in the periodic boundary conditions of the ab initio calculations, such that the electronic transition structure factor in the thermodynamic limit is best approximated by a single finite pe-  Diagram to show algorithmic details and computational efficiency of the method proposed here, structure factor twist averaging (sfTA), in comparison to the conventional method, twist averaging.For the twist averaging method in (a), 100 CCSD calculations are run at random twist angles over a set grid size for a chosen solid.From these, the TA-CCSD energies are calculated and corrections are applied to obtain the final CCSD correlation energy.For the sfTA method in (b), a single twist angle is selected based on how close to the average structure factor of the system it is.This twist angle is then used to run a single calculation on the chosen solid and grid size before the same corrections are applied to obtain the final correlation energy.In both methods, a basis set correction (∆ECBS) is calculated from the smaller grid size and applied to all grids.To correct finite size error either a correction is calculated (∆EFS), or an extrapolation to the TDL can be performed.In (c), the stochastic error in the TA-CCSD energy for a 32 atom supercell of Na is shown for an increasing number of twist angles.This error is not converged (indicating undersampling) for fewer than 40 twist angles.Here, the TA-CCSD energy at 100 twist angles was used as the reference value to calculate these errors.When our sfTA scheme is used, the systematic error is of the order 1 meV/el.By contrast, the stochastic error in twist averaging only reached this value when over 100 twist angles have been used, noting here that the reference value used to obtain these error is the TA-CCSD energy at 100 twist angles.This comparison is evidence for a 100-fold increase in computational efficiency in the CCSD calculation.
riodic system.Our method ensures that degeneracies in the ground state of the unperturbed Hartree-Fock Hamiltonian are automatically lifted, which is crucial for the numerical stability of these self-consistent field calculations of metals.The resulting Hartree-Fock ground state is used as a reference for the coupled cluster theory calculation in a numerically-stable and efficient manner.We apply the outlined method to prototypical metallic lithium and silicon phases.

RESULTS
The quantum mechanical many-electron Hamiltonian of periodic solids is invariant under any translational symmetry transformation that respects the periodicity of the attractive nuclear potential.Crystal momentum vectors k that lie in the first Brillouin zone serve as a label for translational symmetry transformations.The full Hamiltonian can be written as  A comparison between the full twist averaged transition structure factor and the sfTA structure factor as computed from the special twist angle for a 32 atom Na system using CCSD.As can be seen in the graph, the sfTA structure factor appears to lie directly on top of the twist averaged structure factor (MAE = 0.08( 13)) This demonstrates the ability of the method to reproduce the TA-CCSD results for solids.
where the one-and two-electron Hamiltonians for a particular set of k-vectors are given in second quantization by respectively.The indices p, q, r and s refer to occupied and unoccupied one-electron Bloch states characterized by crystal momentum vectors k.The two-electron integrals pq|rs are non-zero only if the corresponding kvectors conserve total crystal momentum.In practice, the thermodynamic limit (TDL), computed by integrating over all k, is approached by sampling the first Brillouin zone using a regular k-mesh with an increasing number of k-points (lim Here, we seek to define an effective Hamiltonian that can also be applied to metals and, in combination with recently proposed finite size corrections [24,34], yields an even more rapid and numerically-stable convergence of correlation energies to the TDL.To this end, we approximate the Brillouin zone integrations in Eq. ( 1) such that dkF (k) Here, F(k) refers to the terms inside the curly brackets in Eq. ( 1).The vector k s corresponds to an offset of the employed k-mesh with respect to the origin, yielding a coupled cluster correlation energy dispersion.Current approaches take an average over different offsets, which increases the computational cost of a calculation by a significant factor.Here, we find that an effective Hamiltonian can be formed and the simulation run at one specific offset only, which we call the effective twist angle.
Our approach is to select one effective twist angle, or k s , at which to perform the calculation by considering the twist angle that best represents the electronic transition structure factor within an first-order Møller-Plesset (MP1) wavefunction formalism.This is represented schematically in Fig. 1.The Fourier components of the transition structure factor can be used to express the projected correlation energy obtained from a wavefunction |Ψ such that Here, the momentum G corresponds to a plane wave vector that is defined as G = g + ∆k, where g is a reciprocal lattice vector and ∆k is the difference between any two crystal momentum vectors that are conventionally chosen to sample the first Brillouin zone.v(G) are the Fourier coefficients of Coulomb kernel with the familiar form of 4π G 2 for excitations allowed by momentum conservation and the prime on the sum implies that the G=0 contribution is treated in an approximate fashion.[27] Our approach is to determine an effective twist angle that most closely approximates the twist-averaged structure factor (S(G) = ks 1 2 .This approach is based on the idea that the single twist angle then represents the average structure factor, not just for MP2, but for other methods as well.The single twist angle is then used to calculate the correlation energy using CCSD.An example of this is illustrated in Fig. 2 and the two structure factors are visually indistinguishable from one another.This small error demonstrates that one twist angle is able to reproduce the transition structure factor averaged across k s .As this method is a way to produce the twist averaged energy based on the best representation of the structure factor, we refer to it as structure-factor-based twist averaging (sfTA).Our protocol to determine the total energy and energy differences between material phases is: (a) run sfTA CCSD calculations at increasingly dense k-point meshes keeping the number of orbitals per k-point constant; (b) for a small k-point mesh (e.g. 2 × 2 × 2) converge to the complete basis set limit using natural orbitals and apply the following correction: ∆E CBS (222) = E CBS (222)-E corr (M, 222); either (c) apply finite size corrections to each mesh using an interpolation scheme [24]; or (d) extrapolate the correlation energy to the TDL using a power-law of 1/N k ; (e) to find the total energy, the correlation energy from steps (a)-(d) is added to Hartree-Fock calculations which have been extrapolated to the TDL.As HF is significantly less expensive than CCSD, larger k-point meshes (up to 6 × 6 × 6) are used for this step.
Figure 3 shows the difference between the correlation energies for the face-centered cubic (fcc) and body-Table I. Equilibrium lattice properties for the two silicon phases, including the transition volume (Vt), the Bulk modulus (B0), pressure derivative of the Bulk modulus (B 0 ), and the volume at equilibrium (V0).The energy difference was calculated as the difference between the minimum of the curves for β-Sn and diamond structures.These fits were used to find the transition pressure (Pt).Our HF and CCSD are shown in comparison to previous quantum Monte Carlo studies [35][36][37][38].We have corrected all of the energy differences except HF and DMC+EMP-pp to include a core-polarization correction of 30 meV/atom taken from Ref. 37; these are also applied to the transition pressures.For DMC+EMP-pp, we use the data including the core-polarization corrections referenced in their paper.Also, the DMC+EMP-pp contains a correction to the energy-volume curve using the empirical (EMP) pseudopotential.Finally, for all methods, the fully corrected transition pressures contain the corrections for zero point energy and finite-temperature vibrational effects and core-polarization corrections.We include these vibrational correction terms in CCSD following the numbers from Ref. 37. Experimental numbers were given by Ref. 39    centered cubic (bcc) phases of lithium.Calculations were performed on a variety of k-point meshes for extrapolation to the TDL.The corrected energies were then extrapolated to the TDL using a power law of 1/N k , where N k is the total number of k-points.For both phases, a correction for the basis set incompleteness error was calculated based on the complete basis set (CBS) limit extrapolation for a 2 × 2 × 2 k-point mesh and a range of basis sets as described in the protocol.
After extrapolation to the TDL shown in Fig. 3 was performed, the Hartree-Fock energy obtained from extrapolating TA-HF energies to the TDL (using k-point meshes up to 6 × 6 × 6) was added.Overall, this resulted in a total energy difference of −0.06 (5) eV/atom between fcc and bcc.This number indicates that the electronic contribution at 0K does not change the relative stability of the bcc versus fcc phases of lithium within the remaining extrapolation error.Overall, the smoothness of this convergence to the TDL seen in Fig. 3 when compared with the Γ-point calculation shows a successful application of sfTA.
We now apply this method to the volume-energy curves of silicon in a diamond and β-tin lattice.These allotropes of silicon are semiconducting and metallic, constituting an excellent test case for the ability of the proposed sfTA technique to treat small-gap systems on the same footing.Furthermore these phases have been studied extensively also by other electronic structure theories including DFT and the random phase approximation (RPA) [42].
We performed calculations using sfTA-CCSD on a 3×3×3 k-point mesh.Finite size corrections were applied using a mesh-interpolation scheme and a complete basis   7) Å3 /atom.The transition between the two phases can clearly be seen to take place between the volumes around 15 to 18 Å/atom.The inset on the graph shows the correlation for both phases for the same range of volumes.The calculations shown were initially run using a 3 × 3 × 3 k-point mesh.Then, finite size and basis set corrections were applied.The β-Sn energies have a 30 meV/atom correction added from the core polarization potential (CPP).These represent the best quality estimates we have for the complete basis set/thermodynamic limit energies for CCSD that we can compute.
set correction was derived from a 2 × 2 × 2 k-point mesh run at the equilibrium volume.[24] The sfTA-CCSD correlation energies were then added to TDL extrapolated TA-HF energies (using k-point meshes up to 5 × 5 × 5) to obtain total energies.The Birch-Murnaghan equation of state was fit to the diamond and β-Sn phases of silicon, shown in Fig. 4, allowing us to obtain equilibrium lattice properties, including equilibrium volumes and bulk moduli.We find that the sfTA-CCSD energy curves are smooth both for the total and correlation energies.The fitting parameters are shown in Table I compared to several highaccuracy quantum Monte Carlo methods from previous studies [35][36][37][38].Experimental values are also included [39][40][41].
We note that the DMC findings compiled in Table I span a relatively large range for different equilibrium properties including equilibrium volumes, bulk moduli, and its derivative.In particular the equilibrium volumes of silicon diamond range from 19.75 Å3 /atom to 20.11 Å3 /atom.The experimental equilibrium volume was measured to be 20.0Å3 /atom.The difference in the DMC estimates can partly be explained by different approximations used to correct finite size errors and the dependence on the employed pseudopotentials.Similar trends as seen for the equilibrium volumes can also be observed for the bulk moduli.The most striking changes between the different DMC calculations can be observed for the energy differences between both phases, which reduces from 505 meV/atom in Ref. 37 to 329 meV/atom in Ref. 38.Due to the strong dependence of the transition pressure between the diamond to β-Sn allotropes on the energy difference, the predictions change from 17.8 GPa to 13.16 GPa, which is gradually improving compared to the experimental estimates that lie between 11.2 and 12.6 GPa.Despite the evolutionary character of the DMC findings summarized in Table I, we consider the latest DMC calculations from Ref. 38 to be the literature benchmark that can be used to assess other electronic structure calculations directly without the necessity to account for other effects that need to be considered when comparing to experiment, such as finite temperature and vibrational entropy contributions.This is corroborated by an independent study using AFQMC which found numbers in good agreement with the spread of DMC results [36].
We now turn to the discussion of results obtained using the quantum chemical methods.In this regard, HF forms the starting point and will also serve as reference for the subsequent CCSD calculations.The HF results listed in Table I demonstrate that HF theory strongly overestimates the diamond equilibrium volume and bulk modulus compared to experiment and DMC, which can mainly be attributed to the neglect of electronic correlation effect.Likewise, the energy difference between the diamond and β-Sn phase is significantly overestimated, yielding a transition pressure that is almost five times larger than the experimental findings.
CCSD theory is able to substantially improve upon the HF description by accounting for correlation effects once we have accounted for finite size as well as basis set corrections.While the CCSD equilibrium volume and bulk modulus for the diamond phase are in good agreement with experiment and latest DMC estimates, we find that the predicted transition pressure are still overestimated by several GPa.Our calculated CCSD energy difference between the two phases was found to be 0.609 eV/atom, with a transition pressure of 20.94 GPa as seen in Table I.This energy difference was found to be larger than the largest DMC energy difference (0.505 eV/atom) by 0.104 eV/atom.This is promising because CCSD is the simplest correction to HF in the coupled cluster hierarchy of methods and these calculations took only 4 days (per volume point) on 16 cores.
The remaining discrepancy between CCSD and DMC is not surprising.Although CCSD forms an important step towards chemical accuracy, it is known from quantum chemical calculations of molecules that the inclusion of perturbative triple particle-hole excitation operators, as performed by CCSD(T), is needed for chemically accurate reaction energies [43].However, we have previously shown that the method is divergent due to its perturbative component [16].One alternative possibil-ity to improve upon CCSD theory has recently arisen in the literature with the distinguishable cluster theory (DCSD) [44].When we run this at the mid-point volume we find that DCSD lowers the energy gap between the phases by 68 meV/atom to 541 meV/atom.If this were uniform across the volume curve, this is equivalent to lowering the transition pressure to 18.6 GPa.
Taken together, these numbers are in reasonable agreement with the DMC energies.The improvement from CCSD to DCSD also demonstrates the way in which quantum chemical methods can be improved.In closing, we do agree with the authors of the DMC study in Ref. 37, who conclude that the numbers from electronic structure calculations were only in agreement with experiment after accounting for vibrational and anharmonic effects.In practice, this means that the structural transition pressure should be a little higher than experiment if the electronic structure is being treated appropriately.

DISCUSSION
In conclusion, we developed a new method to determine an effective Hamiltonian for treating metallic systems using quantum chemical methods.This method, which we call sfTA, finds the single twist angle which best reproduces the twist-averaged transition structure factor.Then, a single calculation can be used to compute an energy that is very close to the twist-averaged CCSD energy with significant cost reduction.The selected twist angle can be found efficiently regardless of the band gap of the system.Our proof-of-principle results come from studying the phases of lithium and silicon.Here, we used the effective Hamiltonian in conjunction with several other recently-developed methods to converge the calculations with respect to basis set size and particle number.The fcc and bcc lithium phases were found to be degenerate with each other, indicating relative stabilities between the phases in the thermodynamic limit.For silicon, we found results that were in reasonable agreement with literature benchmarks from quantum Monte Carlo.We found that the distinguishable cluster approximation reduced the energy between the two phases resulting in agreement with DMC energies.Overall, these results demonstrate that our method is capable of obtaining good CCSD results for a range of metallic applications.
Coupled cluster theory benefits from such an approach more than quantum Monte Carlo for two significant reasons.The first is that basis set incompleteness error from a twist-averaged (or balanced) description of correlation commutes more reliably with the electron number.This provides the crucial benefit of reducing the basis set size required to treat metals and materials with small band gaps.The second is that TDL corrections depend upon the same relationship, and the increased commutiv-ity also helps these corrections become more consistent.Taken together, the small change in the Hamiltonian's symmetry affects the calculation dramatically: calculations that previously were failure-prone are now routinely feasible.

METHODS
Codes and methods:-All calculations were performed using VASP [45,46] and the projector augmentedwave (PAW) method in a plane wave basis set [47].The corrections for the finite size effects for lithium and silicon were performed using cc4s interfaced with the VASP code [24].Twist averaging was used to obtain the Hartree-Fock data.All other calculations were run using the sfTA method.Canonical Hartree-Fock (HF) orbitals were used for all MP2 calculations.For coupled cluster, we used approximate natural orbitals, estimated from MP2 natural orbitals [48], for the lithium basis set convergence calculations and all of the silicon phase calculations.The rest of the calculations used canonical HF orbitals.
Calculation details:-We used a plane wave energy cutoff (ENCUT) of 400 eV throughout.For the initial/preliminary sodium calculations (and the carbon and silicon calculations in the SI), we used an auxilliary plane wave basis set cutoff (ENCUTGW) of 150 eV.The basis set was truncated to 48 orbitals per k-point (using the NBANDS input); these calculations were all 1×1×1 supercells.For lithium, these variables were used to perform calculations at k-point meshes of 2 × 2 × 2, 2 × 2 × 3, 2 × 3 × 3, and 3 × 3 × 3. Silicon followed likewise with the original calculations at 3 × 3 × 3.
In order to calculate the complete basis set energy for each lithium and silicon phase, NBANDS was changed to give 16, 24, 32, 40 and 48 orbitals per k-point.For these calculations, ENCUTGW was set to 400 eV in order to ensure the initial basis set was large enough that the number of orbitals designated by the NBANDS input was the limiting factor for truncating the basis set in all calculations.The lithium bcc calculations used a different NBANDS of 24 orbitals per k-point in the HF calculation preceding the approximate natural orbital calculation (when varying the NBANDS used in the coupled cluster part of the calculation) which is the same number of k-points per atom as fcc.
To obtain converged energies for the CCSD and MP2 calculations, an energy difference of 1×10 −4 eV was used for the Li, Na and C calculations.The Si phases required a smaller energy difference of 1 × 10 −6 eV to obtain fully converged energies.
Structure factor:-We introduced the structure factor in the main manuscript in Eq. ( 4), where it was defined in relation to the coupled cluster wavefunction.In describing the structure factor below, we follow the notation of Ref. 27 for mathematical consistency.The CC wavefunction has amplitudes which can be written T ab ij = t ab ij + t a i t b j .Here, the i and j indices refer to occupied orbitals while a and b indices refer to unoccupied orbitals.Then t a i and t ab ij are singles and doubles amplitudes, respectively.Then the transition structure factor can be constructed as: The coefficients C a i (G) arise from the Fourier transform of the co-densities of two orbitals i and a and the electron repulsion integrals can also be written in terms of these intermediates [27].
Lattice constants:-The equilibrium lattice constants for β-tin silicon and diamond silicon (which are 4.9 Å and 5.761 Å respectively) were scaled by factors in the range 0.85 to 1.1 to produce a range of volumes.Lithium structure information, including equilibrium lattice constants, were obtained from the NOMAD Encyclopedia data base [49].
Numerical analysis:-For the extrapolations to the TDL and the Birch-Murnaghan fits the numpy and scipy libraries were used with Python 3.7.3.Error bars on the TDL extrapolations reflect standard error in the fit parameters.

Figure 1 .
Figure1.Diagram to show algorithmic details and computational efficiency of the method proposed here, structure factor twist averaging (sfTA), in comparison to the conventional method, twist averaging.For the twist averaging method in (a), 100 CCSD calculations are run at random twist angles over a set grid size for a chosen solid.From these, the TA-CCSD energies are calculated and corrections are applied to obtain the final CCSD correlation energy.For the sfTA method in (b), a single twist angle is selected based on how close to the average structure factor of the system it is.This twist angle is then used to run a single calculation on the chosen solid and grid size before the same corrections are applied to obtain the final correlation energy.In both methods, a basis set correction (∆ECBS) is calculated from the smaller grid size and applied to all grids.To correct finite size error either a correction is calculated (∆EFS), or an extrapolation to the TDL can be performed.In (c), the stochastic error in the TA-CCSD energy for a 32 atom supercell of Na is shown for an increasing number of twist angles.This error is not converged (indicating undersampling) for fewer than 40 twist angles.Here, the TA-CCSD energy at 100 twist angles was used as the reference value to calculate these errors.When our sfTA scheme is used, the systematic error is of the order 1 meV/el.By contrast, the stochastic error in twist averaging only reached this value when over 100 twist angles have been used, noting here that the reference value used to obtain these error is the TA-CCSD energy at 100 twist angles.This comparison is evidence for a 100-fold increase in computational efficiency in the CCSD calculation.

Figure 2 .
Figure2.A comparison between the full twist averaged transition structure factor and the sfTA structure factor as computed from the special twist angle for a 32 atom Na system using CCSD.As can be seen in the graph, the sfTA structure factor appears to lie directly on top of the twist averaged structure factor (MAE = 0.08(13)) This demonstrates the ability of the method to reproduce the TA-CCSD results for solids.

Figure 3 .
Figure 3.Comparison between Γ-point and sfTA correlation energy for the two phases of lithium.The comparison, shown as a difference between the face-centered cubic (fcc) and body-centered cubic (bcc) phases, is graphed against kpoints ranging from 2 × 2 × 2 to 3 × 3 × 3.All k-pionts have a basis set correction applied.As can be seen from this comparison, the sfTA results show a much smoother convergence of the difference with a difference of −0.07(2) eV/atom at the TDL.The numbers in parentheses are the errors which come from the fitting algorithm.

Figure 4 .
Figure 4. Total energies for the diamond and β-tin phases of Si across a range of volumes.We can clearly see the minima for both phases, with diamond reaching a minimum at 20.0(4)Å3 /atom and β-tin reaching a minimum at 15.97(7) Å3 /atom.The transition between the two phases can clearly be seen to take place between the volumes around 15 to 18 Å/atom.The inset on the graph shows the correlation for both phases for the same range of volumes.The calculations shown were initially run using a 3 × 3 × 3 k-point mesh.Then, finite size and basis set corrections were applied.The β-Sn energies have a 30 meV/atom correction added from the core polarization potential (CPP).These represent the best quality estimates we have for the complete basis set/thermodynamic limit energies for CCSD that we can compute.