Towards efficient and accurate \emph{ab initio} solutions to periodic systems via transcorrelation and coupled cluster theory

We propose a streamlined combination scheme of the transcorrelation (TC) and coupled cluster (CC) theory, which not only increases the convergence rate with respect to the basis set, but also extends the applicability of the lowest order CC approximations to strongly correlated regimes in the three dimensional uniform electron gas (3D UEG). With the correct physical insights built into the correlator used in TC, highly accurate ground state energies with errors $\leq 0.001 $ a.u./electron relative to the state-of-the-art quantum Monte Carlo results can be obtained across a wide range of densities. The greatly improved efficiency and accuracy of our methods hold great promise for strongly correlated solids where many other methods fail.


I. INTRODUCTION
The coupled cluster (CC) methodologies [1][2][3] at the level of singles and doubles (CCSD) and perturbative triples (CCSD(T)) [4] have become the de facto standard of single-reference ab initio quantum chemistry, and can be applied to systems consisting of hundreds of electrons [5][6][7][8]. In the past few years, these methods have also shown promise in applications to the solid state [9][10][11][12][13][14], although significant challenges remain before they can be routinely applied, as for example density functional theories are. On the one hand, because of quite steep computational scaling (O(N 6 ) and O(N 7 ) for CCSD and CCSD(T) respectively), it is desirable to keep the methods at the lowest possible CC level, namely CCSD, whilst maintaining accuracy. The more accurate CCSD(T), as a perturbative correction to CCSD, additionally fails for metals [15]. It is also desirable that the CC methods can be extended to more strongly correlated systems, where the single reference nature of these approximations breaks down. There have been various attempts to develop modified CCSD methods with a higher accuracy for weakly [16][17][18][19] and strongly [20][21][22][23][24] correlated systems. The distinguishable cluster singles and doubles (DCSD) [25,26] is one such method, which has shown promise in improving CCSD in weakly and strongly correlated molecular systems [27][28][29].
In a separate development, there has been renewed interest in so-called transcorrelated (TC) methods [30][31][32][33][34][35][36][37][38][39][40][41][42][43][44], based on Jastrow factorisation of the electronic wavefunction, which result in effective similarity transformed (ST) Hamiltonians [40,43]. Although TC methods were originally proposed as a way to accelerate basis set conver-gence in electronic wavefunctions, it has become apparent that such similarity transformations can also be extremely helpful in the context of strongly correlated systems. For example, in the repulsive 2D Fermi-Hubbard model, it was found that, with a suitable Gutzwiller correlator, extremely compact forms of ground state right eigenvectors of the ST Hubbard Hamiltonian could be obtained [42], dominated by the Hartree-Fock (HF) determinant. Since single-reference coupled cluster methods work best when the wavefunction is dominated by the HF determinant, and furthermore, since the coupled cluster method can itself be cast in terms of a similarity transformation of the Hamiltonian, it is natural to ask if the two concepts -coupled cluster and transcorrelation -can be usefully combined into a single framework, whereby the compactification generated by the TC method is exploited by the CC method, to extend its range of applicability into more strongly correlated systems. The purpose of this paper is to report such an investigation, applied to the uniform electron gas (UEG), over a broad range of densities. The 3D UEG model assumes that the background is evenly and positively charged, and that the electrons interact with each other via the Coulomb interaction. As simple as it is, the UEG possesses an intricate phase diagram [45,46], which can only be accurately described by theories that perform consistently well over a broad range of densities. Historically, the UEG model has also played an important role in the development of many useful approximations. For example, several successful local and gradient-corrected density functionals [47][48][49] are based on the UEG; the random phase approximation (RPA) [50,51] was developed in a pursuit of understanding metals using the UEG as a model. In recent years, the UEG has attracted studies from various highly accurate ab initio methods and spurred the development of several new methods [15,35,40,[52][53][54][55][56][57].
When applying CC to the UEG, we work in a plane arXiv:2103.03176v1 [cond-mat.str-el] 4 Mar 2021 wave basis; momentum conservation then excludes all single excitations from the CC ansatz, greatly simplifying the resulting amplitude equations. As a result, the TC Hamiltonian can be treated with relative ease, allowing us to investigate whether the CC method can be beneficially applied to the TC Hamiltonian. We will investigate the CCD and DCD approximations, in the context of the TC Hamiltonian, and show that with a suitable form of the correlator, the basis set convergence can be greatly accelerated (as expected), but in addition highly accurate energies can be obtained across a broad range of densities 0.5 ≤ r s ≤ 50, covering both the weakly and strongly correlated regimes. This gives us confidence that the method, once suitably generalised to real systems (which will need to include the singles contribution), will allow a highly accurate yet efficient methodology for the solid state.
In the rest part of this paper, we review the UEG model, CC and TC theories in Sec. II A, II B and II C, respectively; in Sec. II D we discuss the important approximations made to the TC Hamiltonian; we demonstrate our scheme for choosing the optimal parameters in the correlator in Sec. III; we showcase and discuss our TC-CCD/DCD results in comparison with benchmark data in Sec. IV; and finally we conclude the paper in Sec. V with some outlooks for future directions.

A. Three Dimensional Uniform Electron Gas
The 3D UEG is the simplest model for realistic periodic solids, of which the Hamiltonian in real space readŝ where the const. includes the interactions between electrons and the homogenous positive background charge, and the interactions between the electrons and their own periodic images, which is termed as the Madelung constant and will disappear as the size of the simulation cell goes to infinity. Atomic units are used to simplify the equations. When plane wave basis functions and a simple cubic simulation cell of volume Ω = L 3 are used, we can reformulate the Hamiltonian in a second-quantized form, where for simplicity we use p, q, r, s . . . indices as a compact form for the general momentum (plane wave basis function) indices k p , k q , k r , k s . . . and hereon we use the two terms plane wave basis function and orbital equivalently. We stress that due to momentum conservation, i.e. k ≡ k r −k p = k q −k s , there are only three free indices among pqrs, and the interactions with the homogenous positive background charge are cancelled by the divergent Coulomb potential at k = 0, which is defined as V rs pq ≡ V (k) = 4π Ωk 2 . We also ignore the Madelung contribution in the Hamiltonian which can be added posteriorly to the ground state energy. The electron density of the system can be described by the Wigner-Seitz radius r s = 3 4πN 1/3 L, where N is the number of electrons.

B. Coupled/Distinguishable Cluster Doubles
In the CC ansatz, we let the many-electron ground state wavefunction to be where Φ 0 is the HF wavefunction andT is a cluster operator. In the case of the UEG, we work in a plane wave basis and Φ 0 is given by the Fermi sphere. We will investigate the CCD approximation and its distinguishable variant (DCD) [25,58], which is based on a modification of the CCD amplitude equations by neglecting intercluster exchange diagrams and ensuring the particle-hole symmetry and exactness for two electrons. Alternatively, DCD can be derived from screened Coulomb considerations [59]. We start with the canonical CCD and later highlight the differences between CCD and DCD. In CCD, the full cluster operator is approximated by the doubles excitations only, where T ij ab are the doubles amplitudes. Following convention, we use i, j, k . . . and a, b, c . . . to represent occupied and unoccupied orbitals in the reference determinant, respectively. Again, the momentum conservation ensures that only 3 indices of the amplitude tensor are free, saving a great deal in storing them in the computer memory.
The T 2 amplitudes are obtained by solving the projective doubles amplitude equations, where Φ ab ij are doubly substituted determinants. To be specific, a functional form of the residual, which unifies CCD and DCD for closed shell systems, can be written as where we define the permutation operatorP(ia; jb)T ij ab ≡ T ij ab + T ji ba and the following intermediates, We note that in this case the Fock matrix f p q is diagonal, with the diagonal elements being the orbital energies p . A straightforward way to update the T 2 amplitudes at iteration n + 1 will be Of course, more advanced iterative schemes can be used, e.g. DIIS [60,61], to accelerate convergence rate. Using the converged T 2 amplitudes, the correlation energy is expressed as and the total energy is where E HF = Φ 0 |Ĥ |Φ 0 .

C. Transcorrelation
In the transcorrelation framework the many-electron wavefunction is written as is a correlator consisting of pair correlations u(r i , r j ), whose form will be discussed later. Φ should satisfy the similarity-transformed eigenvalue equation It is worth pointing out that at this stage, no approximations have been made, and the spectra E ofĤ tc are the same as of the original Hamiltonian. As shown in Ref. [40], the second quantized form of theĤ tc in a plane wave basis iŝ The TC Hamiltonian has additional 2-body and 3-body interactions. Due to one of the additional 2-body interactions, the TC Hamiltonian is non-hermitian. This fact can pose some difficulties for variational methods, but not so for projection methods such as full configuration interaction Monte Carlo (FCIQMC) [9,62] and CC.

D. Approximations to the 3-body operator
The additional 3-body operator when treated without approximations will increase the computational scaling of CCD or DCD from N 6 to N 7 . To seek a good balance between the computational cost and the accuracy, we include only up to effective 2-body operators arising from normal-ordering the 3-body operator. In this approximation, only the normal-ordered 3-body operator is excluded. We can justify this approximation by analogy to the HF approximation, which constructs a mean-field solution by including only the single and double contractions from the Coulomb operator. In cases where the mean-field approximation is reasonably good, the contribution of the missing normal-ordered Coulomb operator is small, compared to the single and double contractions. In contrast to the HF approximation, the parameters in the correlator in general allow a tuning of the strength of the missing normal-ordered 3-body operator, which we will discuss in the next section.
In general, we can write our approximated Hamilto-nian asĤ where E T refers to the triply-contracted 3-body operator contribution,ω p is the doubly-contracted 3-body integral andω rs pq is the singly-contracted 3-body integral. The curly brackets indicate that the operators are normalordered with respect to the HF vacuum (Fermi sphere). We emphasize that in Eq. (22)Ẽ HF and˜ p are calculated now with the modified 2-body integralsṼ rs pq = w rs pq + V rs pq . For clarity, we outline the procedures for our TC-CCD/DCD framework.
1. Evaluating ω rs pq and V rs pq and combining them intõ V rs pq ← w rs pq + V rs pq ; 4. Evaluatingω p , and defining p ←˜ p +ω p ; 5. Evaluating the singly-contracted 3-body integral w rs pq and redefining V rs pq ←Ṽ rs pq +w rs pq ; 6. Solving the usual CCD/DCD amplitude equations using p and V rs pq for T 2 and obtaining E c ; 7. The total energy is E =Ẽ HF + E T + E c .
For details on the mathematical expressions for the contractions of the 3-body operator, we refer to the Supplementary Information.

III. CHOICE AND OPTIMIZATION OF THE CORRELATOR
Past experience with the TC method has shown that the form of the correlatorτ is of extreme importance in the TC method, otherwise the benefit of transcorrelation is lost -Φ can be simpler than Ψ only if the correlator captures the correct physics of the pair correlations in the system. An inappropriate correlator can actually lead to a harder problem than the original Schrödinger equation. In our previous study of the exact TC method in the UEG, we proposed a form of the correlator (shown below) which was found to work successfully in accelerating convergence to the basis set limit, without changing the correlation that could be captured with the basis set by a FCI level Φ function. In the present study, since we will be approximating the Φ with the CC ansatz, we additionally require the correlatorτ to capture some of the correlation inside the Hilbert space.
Here we propose a physically motivated correlator that mimics the behavior of the correlation hole between two unlike-spin electrons as r s varies in 3D UEG. The correlation hole can be examined by the pair-correlation function g(r 12 ) in real space, as studied in Ref. [63], which shows that the correlation hole between two unlike-spin electrons grows deeper and wider as the Wigner-Seitz radius r s increases or as the electron density decreases. Fig. 1 provides a sketch of the Jastrow factor with our proposed correlator u(r 12 ) as the exponent, which captures the desired behavior. We point out that the func- tional form of this correlator, which reads in real and reciprocal space respectively as where si(x) = − ∞ x sin(t) t dt, was first reported in Ref. [40] to satisfy the cusp condition between two electrons with opposite spins at short inter-electron distance and its influence is reduced to nonexistence as the complete basis set (CBS) limit is reached. This was done by choosing k c to be the same as the plane wave cutoff momentum, k F , which defines how many plane waves are included as basis functions. In contrast, to mimic the behavior of the pair-correlation function, as a first attempt in the present study we choose the parameter in this correlator such that the first nonzero root of Eq. (23) is fixed to be at r s , irrespective of the basis set. This is achieved by setting where R 1 ≈ 2.322502989. This choice can be rationalized by the physical picture that at lower densities, electrons prefer to stay further away from each other. Furthermore, this correlator, regardless of the choice of k c , retains the cusp condition for two electrons with unlikespins at r = 0 [40] and should increase the convergence rate of the computed energies with respect to the employed basis set towards the CBS limit.
To further justify the choice of this correlator, we show that for UEG with 14 electrons at r s = 5, where traditional CCD exhibits a large error, the most compact expansion of the wavefunction in Slater determinant space is reached at this value of k c . In Fig. 2, we show the weights of the HF determinant extracted from TC-FCIQMC simulations using different k c values, without making approximations to the 3-body operator. We note that due to the discrete momentum mesh as a result of using a finite simulation cell, the possible choices are k c = 2π  However, we find that this intuitive choice of k c is not always the optimal, especially at extremely low density regimes. It is reasonable to expect that the optimal k c for those systems should deviate slightly from R1 rs . So we scan a range of k c values around it to locate the one that minimizes the norm of the closed-shell amplitudes for double excitations of electrons with opposite spins, T ↓↑ 2 , in the TC-CCD/DCD calculations with a small basis set, see Fig. 3 [65].
Ideally, two separate correlators should be used for electrons with parallel and anti-parallel spins, and their parameters should be optimized simultaneously using the norm of the full amplitudes in a similar manner. For the present study, we argue that the correlations between two parallel-spin electrons are dominated by the exchange effects, which are already captured by the anti-symmetry in the Slater determinants. Therefore, we focus on capturing the correct physics between electrons with opposite spins in the correlator, i.e. the changing depth and width of the correlation hole as a function of r s [63], and minimizing the corresponding amplitudes in the CC ansätze. Indeed, we found in practice the minima in T ↓↑ 2 as a function of k c are more pronounced, and thus easier to spot than those in the norm of the full amplitudes, T 2 . We stress that this compact form of wavefunction at the optimal choice of k c should greatly benefit approximate methods like CCD and DCD, whose accuracy relies on the assumption that the true ground state wavefunction is compact around the reference determinant, which is normally chosen to be the HF determinant. , as a function of kc, calculated by TC-DCD method. The 14-electron and 54-electron systems use a basis set including 57 and 257 plane waves, respectively. All possible contractions from the 3-body interactions are included, excluding the normal-ordered 3-body interactions. The solid horizontal color lines represent the T ↓↑ well under control in that it scales approximately as u 2 (k)k 2 ∼ 1 k 6 . We note in passing that if we choose correlators that do not truncate at small k, such as the Yukawa-Coulomb correlator in Ref. [66] or the Gaskell correlator in Ref. [67], the iterative solution of the amplitude equations becomes too unstable to converge at low densities. We attribute this instability in these cases to the large missing normal-ordered 3-body interactions, similar to the instability in a HF self-consistent solution when the missing normal-ordered Coulomb interactions are large.
The FCIQMC calculations are carried out using the NECI program [68]. The CCD and DCD along with the TC integrals are implemented in a Python program using the automatic tensor contraction engine CTF [69] and the NumPy package [70]. We first examine the basis set convergence behavior of TC-CCD/DCD compared to the canonical ones. In Fig. 4 we present the total energy per electron relative to the extrapolated value for each method, retrieved as a function of the inverse of the employed number of plane waves, 1/M . As mentioned before, our correlator satisfies the cusp condition at the coalescence point of two electrons with opposite spins. So the accelerated convergence behavior in the TC methods compared to the canonical ones is not surprising. The acceleration is the most obvious at high density regimes, since at low densities the required number of basis functions to reach convergence in both the TC and non-TC methods is relatively small. These observations are consistent with those of the TC-FCIQMC reported in Ref. [40].
The optimal k c values, (TC-)CCD/DCD energies at CBS and the benchmark data are listed in Table I and II for the 14-and 54-electron 3D UEG, respectively. In Fig. 5 we present the errors of total energies per electron calculated by TC-CCD, TC-DCD, CCD and DCD relative to the most accurate FCIQMC [40,53,55] and backflow DMC (BF-DMC) [54,71] results on the 14-and 54-electron 3D UEG. The finite basis set errors in our methods have been carefully eliminated by extrapolation to the CBS limit.
In general the accuracies of the TC methods are greatly improved compared to their canonical counterparts, especially in regions (r s = 5 − 50) where the latter exhibit the largest errors. More importantly, the improved accuracies are retained when going from the 14-to the 54electron system. We highlight that the TC-DCD achieves an accuracy of ≤ 0.001 a.u./electron across a wide range of densities, i.e. r s = 0.5 − 50 for the 14-electron and r s = 0.5 − 20 for the 54-electron 3D UEG, with an exception at r s = 10 for the latter where it drops slightly out of the 0.001 a.u./electron accuracy. We argue that with the next possible smaller value of k c = √ 6, which yields a marginally lower T ↓↑ 2 , instead of the current choice of k c = 2 √ 2, the 0.001 a.u./electron accuracy at r s = 10 can be regained. The discrete grid of the k-mesh makes it hard to pick the optimal k c in Fig. 3(b). However, as the system gets larger and the k-mesh gets finer, the T ↓↑ 2 as a function of k c will also be smoother, and the choice of the optimal k c will become more definite. We use colorful shaded areas in Fig. 5 to reflect the uncertainties due to the possible choices of k c which yield similar T ↓↑ 2 values in Fig. 3. At high densities, i.e. r s = 0.5 − 2, the canonical DCD is already very accurate, and the main benefit from TC there is in accelerating the basis set convergence. Overall, DCD exhibits smaller errors than CCD, which agrees with earlier comparative studies between DCSD and CCSD [25][26][27].

V. CONCLUSIONS
We demonstrated that the correlator Eq. (23), used with transcorrelated coupled cluster theory, drastically improves the accuracy of approximate methods, i.e. CCD and DCD, for 3D UEG across a wide range of densities. The basis set convergence rate is also improved thanks to the fact that the correlator satisfies the cusp condition at the coalescence point of two unlike-spin electrons. We have explored the mechanism behind the improved accuracy of the TC-CCD and TC-DCD methods, which is related to a compactification of the many-electron wavefunction in Slater determinant space when the dominant pair correlations between electrons with unlike spins are directly included in the correlator. The optimization of the parameter in the correlator is seamlessly incorporated within the TC-CCD/DCD framework, without requiring an external algorithm. We notice that a range-separation  scheme of CCD can also achieve similar accuracy in 3D UEG, but without improving the basis set convergence rate [54]. Comparatively speaking, our methods are systematically improvable, in that a more flexible form of the correlator can be designed by a combination of a series of functions [43] or by a general function approximator, e.g. an artificial neural network, to include further correlation effects such as nuclear-electron correlations and correlations between parallel pairs of electrons. Other systematic ways of optimizing the correlator in combination with VMC [43] can also be explored. When going to real periodic solids, TC-CCSD and TC-DCSD will be needed; extra efforts are also required to compute the additional integrals besides the Coulomb integrals, where the most computationally demanding part is the singly-contracted 3-body integrals which scales like O(N 2 o N 4 v ), where N o and N v are the number of occupied and unoccupied orbitals, respectively. Fortunately, the computation of the extra integrals scales no worse than the CCSD or DCSD algorithm and it can be compensated by the accelerated convergence rate towards CBS limit in the TC framework. These perspectives will be important in extending the encouraging performance of the current TC-CCD and TC-DCD methods from the UEG to real periodic solids with moderate to strong correlation.
Competing Interests The authors declare that they have no competing financial interests.  [40] are used and for the rest BF-DMC data [54] are used for benchmark. (b) Results for 54 electrons. For rs = 0.5 − 1, TC-FCIQMC data [40] are used and for the rest BF-DMC data [71] are used for benchmark. The grey shaded areas stand for the ±0.001 a.u./electron accuracy relative to the reference data. The colorful shaded areas reflect the uncertainties in the TC-CCD and TC-DCD energies due to slightly different choices of the kc values.