Range-separated stochastic resolution of identity: Formulation and application to second-order Green’s function theory

Many-body perturbation theory (MBPT) based on Green’s function (GF) approaches [e.g., the Møller–Plesset (MP) perturbation theory,¹ the second-order Green’s function (GF2) approach,² and the GW³ approximation] have proven very useful in predicting the ground state properties beyond the limitations of density functional theory (DFT) and the Hartree–Fock (HF) method as well as in predicting quasi-particle and neutral excitation. In these methods, correlations are treated systemically by expanding the self-energy (which contains the information of correlations) in the Coulomb^2,4 or screened Coulomb^5–7 interactions. MBPT has been applied to a variety of molecular and bulk systems in predicting, e.g., correlation energies, ionization potentials and electron affinities,^8–22 and excited states.^7,15,23–26 Excluding several recent applications,^27–33 MBPT has been limited to relatively small systems due to the steep computational scaling with the system size.

A particularly interesting implementation of MBPT, relevant to the applications reported below, is based on a second-order approximation to the electron self-energy,^2,4,34 which has received increasing attention in recent years.^8–10,35 In contrast to the GW approximation,¹ dynamical exchange correlations are included explicitly in the GF2 self-energy to the second order in Coulomb interactions, providing accurate ground state energies^36,37 and quasi-particle energies.^9,10,38,39 Although the results of recent studies are extremely promising, the GF2 approach suffers from a high computational cost [O(N⁵)], limiting its application to relatively small system sizes.

To overcome this limitation, two stochastic formulations were recently introduced to reduce the computational scaling. Neuhauser et al.⁴⁰ developed a stochastic decomposition of the imaginary time GF to reduce the overall scaling of GF2 to O(N³). Takeshita et al.⁴¹ and Dou et al.⁴² proposed an approach that builds on the stochastic resolution of identity (SRI) for the electron repulsion integrals (ERIs)⁴³ to describe both ground and quasi-particle excited states. Similar to the deterministic resolution of identity (RI),^44–48 the SRI decouples the four-index ERIs. While the number of auxiliary bases increases with the system size for the RI, the number of stochastic orbitals in the SRI is independent of the system size, resulting in an overall O(N³) scaling. However, the SRI technique comes at the cost of introducing a statistical error in the energy and nuclear forces,^43,49–52 which can be controlled by increasing the number of stochastic realizations, N_s. While the overall scaling of the stochastic formulations of GF2 is similar to DFT and HF, achieving chemical accuracy requires a large number of stochastic realizations, resulting in increasingly longer computational time, even for small systems.^41,42

In this work, we develop a range-separated stochastic resolution of identity (RS-SRI) approach to decouple the four-index ERIs, where the short-range ERIs (larger values) are treated deterministically using the resolution of identity (RI)^44–48 and the remaining terms are treated using the stochastic resolution of identity (SRI).⁴³ The RS-SRI approach allows for a significant reduction in the statistical error without the need to increase the number of stochastic realizations while maintaining the overall O(N³) scaling. We apply the RS-SRI technique to GF2 theory and demonstrate its ability to reduce the overall computational scaling from O(N⁵) to O(N³) as well as increase the sampling efficiency by nearly two orders of magnitude as compared to the SRI technique.

Consider a generic many-body electronic Hamiltonian in the second-quantization representation,

where $âi†$ and â_i are the Fermionic creation and annihilation operators, respectively, for an electron in orbital χ_i(r). In the applications below, χ_i(r) is chosen to be an atomic orbital, but we do not use the locality of the basis to reduce the scaling or do we introduce a cutoff to compute the ERIs [see Eq. (2)] or the overlap matrix [see Eq. (3)]. Therefore, the formalism and the resulting scaling reported below are general for any choice of basis. The creation and annihilation operators obey the following anti-commutation relationship:

where S_ij = ∫χ_i(r)χ_j(r)dr is the matrix element of the overlap matrix S. In Eq. (1), h_ij is the matrix element of the one-body Hamiltonian and $v$_ijkl is the four-index ERI (⁠$v$_ijkl ≡ (ij|kl)),

Describing correlations within a many-body perturbation technique beyond the mean-field approximation relies on the contraction of $v$_ijkl (or powers of $v$_ijkl), a task that becomes computationally intractable with increasing levels of accuracy. A common approach to reduce the computational complexity is based on the resolution of identity (RI), where the four-index ERIs in Eq. (3) are approximated by products of three-index ERIs and two-index ERIs,^44,46–48

Here, χ_A(r) and χ_B(r) are the auxiliary orbitals, and (ij|A) and V_AB are the three-index and two-index ERIs respectively,

For convenience, we define a new set of three-index ERIs $KijQ$⁠,

such that the four-index ERI can be expressed in terms of three-index ERIs only,

The advantage of the above decomposition is that the resolution of identity reduces the number of two-body ERIs from O(N⁴) to O(N²N_aux), where N is the size of the atomic basis and N_aux is the size of the auxiliary basis. However, since N_aux increases nearly linearly with the size of the atomic basis N and since the calculation of $KijQ$ scales as O(N⁴), the approach does not always reduce the computational scaling of the correlation energy for, e.g., MP2 and GF2.^41–43 

Recently, we have introduced a stochastic version of the resolution of identity, which provides a framework to reduce the scaling for contraction within many-body perturbation techniques at the account of introducing a controlled statistical error in the calculated observables (e.g., the forces on the nuclei and the energy per electron). The balance between accuracy and efficiency is controlled by the number of stochastic realizations (N_s) according to the central limit theorem. The stochastic RI approach utilizes the same set of three-index ERIs (ij|A) while circumventing the need to directly compute $KijQ$ by introducing a set of N_s stochastic orbitals, {θ^ξ}, ξ = 1, 2, …, N_s. The stochastic orbitals are defined as arrays of length N_aux with randomly selected elements 1 or −1, i.e., $θAξ=±1$⁠. Defining

the expression for $v$_ijkl can be reduced to

where $⋯θ$ implies a statistical average over the stochastic orbitals, {θ}. The overall computational scaling of the $Rijξ$ matrices is O(N_sN³), but N_s is found to be independent of the system size for different applications.^{25,29,41–43,49,53} The SRI technique has been successfully used to reduce the scaling of the correlation energy within MP2 and GF2 theories, from O(N⁵) to O(N³).

The above approach has been implemented for simple molecules and for hydrogen chains of different lengths in order to assess its accuracy for large systems.^41–43 To converge the results to chemical accuracy required a rather large number of stochastic orbitals (N_s ≈ 1000), which limits the application of the SRI technique to relatively small systems (due to the large “prefactor”), with N → 1000, still exceedingly larger than the deterministic approach.^41–43 In order to reduce the number of stochastic orbitals and to allow for a smaller statistical error, we first sort the ERIs (ij|A) according to their magnitude and keep only those that are larger than a threshold,

Here, ${|(ij|A)|}jmax$ is the maximal value of |(ij|A)| for each j and ϵ′ is a predefined parameter. The superscript L (or S) denotes large (or small) values. By setting the cutoff threshold to depend on $ϵ′N$⁠, the number of nonzero elements in (ij|A)^L for each j scales as O(N) [rather than O(N²) if no threshold is used or O(1) if a fixed threshold ϵ′ is used]. This implies that the total non-vanishing elements in (ij|A)^L scale as O(N²). We then define $[KijQ]L$ as

and keep only the terms that are larger than a predefined threshold, namely, we set $[KijQ]L=0$ for values below the threshold according to

The calculation of $[KijQ]L$ using the above procedure scales as O(N³). We proceed by defining

where $Rijξ$ is defined above in Eq. (9) and the computational scaling for both terms, $[Rijξ]L$ and $[Rijξ]S$⁠, is O(N³). Using these definitions, the four-index tensor $v$_ijkl can be rewritten as

Equation (16) is referred to as the range-separated stochastic resolution of identity (RS-SRI). The RS-SRI reduces to the SRI for ϵ = 1 and to the deterministic RI for ϵ = 0. This suggests that ϵ can be used as a control parameter balancing the computational efficiency and statistical errors. For optimal choices of ϵ, the contribution of $∑QNaux[KijQ]L[KklQ]L$ in Eq. (16) must be larger than the other terms.

We now apply the above formalism to the second-order Matsubara Green’s function (GF2) theory.^40–42 The main entity in the GF2 theory is the Matsubara single-particle, finite temperature, Green’s function given by (we set ℏ = 1 unless otherwise stated)

where â_i and $âj†$ are defined above in Sec. II, T_c is a time ordering operator, and τ is an imaginary time point along τ ∈ (0, −β). In the above, we have used the Heisenberg picture for the operators: $âi(τ)=e(Ĥ−μN^)τâie−(Ĥ−μN^)τ$⁠, where $N^=∑ijSijâi†âj$ is the number operator and Ĥ is the many-body Hamiltonian defined in Eq. (1). The average is taken with respect to the grand canonical partition function $⟨⋯⟩=Z−1Tr(⋯)e−β(Ĥ−μN^)$⁠, where $Z=Tre−β(Ĥ−μN^)$ is the normalization factor, β = 1/k_BT is the inverse temperature, and μ is the chemical potential.

The Matsubara GF obeys the following Dyson equation:

where F is the Fock matrix given by

and Σ is the self-energy. In the second-order Born approximation, the self-energy (in the closed shell case) is given by

The above form scales as O(N⁵) using the appropriate contraction.

The Matsubara Green’s function for the Fermionic systems obeys the following anti-symmetric relationship: G(τ) = −G(τ + β). The anti-symmetry feature allows for a Fourier representation of G(τ) in imaginary frequency,

Here, $iωn=i(2n+1)πβ$ are the Matsubara frequencies, and the inverse Fourier transform is defined by

The Dyson equation [cf. Eq. (18)] can then be solved in the frequency domain,

where $Σ̃(iωn)$ is the Fourier transform of the self-energy [Eq. (20)] and $G̃0(iωn)$ is the non-interacting GF,

Since the self-energy $Σ̃(iωn)$ depends on $G̃(iωn)$ itself, Eq. (23) as well as Eq. (20) must be solved self-consistently. This is done by first performing a Hartree–Fock calculation to obtain the overlap matrix S, the Fock matrix F, and the chemical potential μ. The Fock matrix can then be used for constructing the non-interacting GF [cf. Eq. (24)] that serves as our initial guess of $G̃(iωn)=G̃0(iωn)$⁠. The next step involves the calculation of the self-energy, which is preformed in the imaginary time domain [Eq. (20)]. The self-energy is then used to update the GF in Eq. (23), and the latter is used to update the Fock matrix in Eq. (19). It is often necessary to conserve the number of particles N_e = −2∑_ijG_ij(τ = β⁻)S_ij. This can be achieved by tuning the chemical potential μ.

The computational bottleneck in GF2 is the calculation is the self-energy, which scales formally as O(N⁵). Using the RS-SRI representation for $v$_ijkl given by Eq. (16), the self-energy can be written as

In Sec. IV, we apply the RS-SRI to a series of hydrogen chain molecules and compare the results to deterministic RI as well as to SRI. We find, in practice, that the RS-SRI scales even better than the upper theoretical limit of O(N³) and at the same time reduces the statistical error by about an order of magnitude as shown below.

In this section, we assess the performance of the RS-SRI-GF2 approach and compare the results to deterministic GF2 and SRI-GF2 for hydrogen dimer chains $HNH$ of length N_H and on a series of water clusters. We begin with the hydrogen chain that has been the canonical model system to assess different electronic structure methods. The distance between strongly bonded hydrogen atoms was set to 0.74 Å, and the distance between weakly bonded hydrogen atoms was set to 1.26 Å. For each hydrogen, we used the STO-3G basis and the CC-pVDZ-RI fitting basis for the resolution of identity in evaluating the self-energy as well as the CC-pVDZ-JKFIT fitting basis in evaluating the Fock matrix in Eq. (19). The inverse temperature used for the calculation of the GFs was set to β = 50 inverse Hartree, sufficient to converge the results due to the large quasi-particle gap. We used the approach developed in Ref. 40 to perform the discrete Fourier transform with 20 000 Matsubara frequencies and 300 imaginary-time points.

We first test how ϵ′ and ϵ [in Eqs. (11) and (13)] affect the final results. In Fig. 1, we plot the correlation energy per electron, defined as⁴⁰ 

for H₁₀ as a function of ϵ for two values of ϵ′. The statistical errors were estimated from the standard deviation σ over N_run = 10 independent RS-SRI-GF2 runs. The RS-SRI-GF2 results agree with the deterministic GF2 results within the statistical error for all values of ϵ and ϵ′. In general, we find that the statistical errors decrease with the decreasing values of ϵ and ϵ′, consistent with the deterministic limit defined by ϵ → 0 and ϵ′ → 0, respectively. However, for the larger value of ϵ′, we find an optimal value of ϵ, suggesting that the approach to the deterministic limit may not be monotonic. Clearly, one needs to be careful when choosing the optimal thresholds for range separation. Below, for hydrogen molecule chains, we will set ϵ′ = 0.02 and ϵ = 0.1 as our thresholds for RS-SRI calculations unless otherwise noted.

In Fig. 2, we plot the correlation energy per electron for a series of hydrogen dimer chains. We compare the results obtained using the RS-SRI-GF2 with SRI-GF2 and for small systems, with deterministic calculations. We find, as expected, that the correlation energy per electron is roughly independent of the length of the chain. Furthermore, both RS-SRI-GF2 with SRI-GF2 agree with the deterministic results within their statistical error. However, the statistical error for the same number of stochastic orbitals (N_s) is significantly smaller (by nearly an order of magnitude) for RS-SRI-GF2 compared to SRI-GF2 for the entire range of system sizes. The error bar was estimated as the standard deviation σ, where we have used N_run = 10 runs to estimate the statistical fluctuations.

In Fig. 3, we plot the correlation energy per electron as a function of the inverse of the number of stochastic orbitals (⁠$1Ns$⁠) for H₈₀, H₂₀₀, and H₅₀₀. We find that the statistical fluctuations decrease as $1Ns$⁠, indicated by the decrease in the magnitude of the error bars. For N_s = 800, we compare the RS-SRI-GF2 with the SRI-GF2 (red symbol, Fig. 3) for H₅₀₀. Clearly, the statistical noise is much larger (by about a factor of 10) compared to the RS-SRI-GF2 result (green symbols). We also find that the statistical fluctuations in the correlation energy per electron are independent of the system size, consistent with N_s being independent on N, similar to other implementation of stochastic orbital techniques.^{25,29,41–43,49,53}

A closer examination of the results shown in Fig. 3 suggests that for H₅₀₀, the correlation energy per electron depends linearly on 1/N_s. This dependence has been discussed and analyzed in Ref. 40 and was attributed to the existence of a bias resulting from the nonlinear dependence of the correlation energy on the stochastic process (see the  Appendix for a complete derivation of the bias). In comparison to a previous work,⁴⁰ we find that the bias is rather small, well within the statistical errors, and thus, its existence is questionable.

In Fig. 4, we plot the computational wall time of the different GF2 approaches (deterministic GF2, RS-SRI-GF2, and SRI-GF2) as a function of the length of the hydrogen atom chain, N_H. All calculations are performed on a single node with the 32-core Intel-Xeon processor E5-2698 v3 (“Haswell”) at 2.3 GHz. The deterministic GF2 scales as O(N^5.1), the SRI-GF2 scales as O(N^3.1), and the current approach, for the same level of accuracy as in the SRI-GF2, scales as O(N^2.2), slightly better than the theoretical limit of O(N³). Note that the RS-SRI-GF2 approach has a much smaller total wall time compared to the other approaches, across the entire system range studied. As additional checks, the inset of Fig. 4 shows the scaling of computing $[KijQ]L$ as well as the scaling of the deterministic portion of the self-energy (terms that only involve $[KijQ]L$ but not $RijL$ or $RijS$⁠). The former scales as O(N³), and the latter is found to scale as O(N²).

Finally, we demonstrate that the range-separated stochastic resolution of identity method is also useful for three dimensional systems. In Fig. 5, we plot the correction function per electron for a series of water clusters. We set ϵ′ = 0.005 and ϵ = 0.025 in the RS-SRI-GF2 calculations. Note that, again, the SRI-GF2 as well as the RS-SRI-GF2 results agree with the deterministic GF2 results within the statistical errors, which are significantly smaller for the RS-SRI-GF2 approach. The statistical error is estimated by the standard deviation over N_run = 10 runs. The computational scaling for the water cluster is found numerically to be O(N^2.4) (not shown), slightly better than the theoretical limit of (N³). This implies that the range-separated stochastic resolution of identity approach is also useful to reduce the computational scaling as well as statistical error for 3D structures.

We have developed a range-separated stochastic resolution of identity approach to decouple the four-index electron repulsion integrals and implemented the approach within the second-order Green’s function formalism, GF2. The RS-SRI technique can be viewed as a hybridization of the RI and SRI techniques, leveraging from both the accuracy of the RI and the reduced computational complexity of the SRI approaches. Results calculated for hydrogen dimer chains of varying lengths and for water clusters show an improved scaling of O(N^2.2) with the size of the basis, N. In comparison to our previous fully stochastic approach, the RS-SRI-GF2 approach significantly reduces the statistical error, resulting in computational wall times that are nearly two orders of magnitude shorter compared to the SRI-GF2. While we focused in this work on the specific implementation of the RS-SRI, the approach lends itself to higher-order approximations to the self-energy and for going beyond ground state properties. Future work should assess the performance of this RS-SRI technique for a wider range of geometries as well as its applicability to the calculation of excited state properties.

W.D. and M.C. contributed equally to this work.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

D.N. and E.R. are grateful for support from the Center for Computational Study of Excited State Phenomena in Energy Materials (C2SEPEM) at the Lawrence Berkeley National Laboratory, which is funded by the U.S. Department of Energy, Office of Science, Basic energy Sciences, Materials Sciences and Engineering Division under Contract No. DEAC02-05CH11231 as part of the Computational Materials Sciences Program. R.B. is grateful for support from the U.S.–Israel Binational Science Foundation (Grant No. BSF-2020602). The resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231, are greatly acknowledged.

In this appendix, we provide a short analysis of the origin of the bias and its interplay with the statistical error. Consider the general case where g is the observable of interest given by

where G is some smooth function and ⟨x⟩ is the expected value of a random variable x, which also has the variance σ² = ⟨(x − ⟨x⟩)²⟩,

Using N_s independent samples, x_n, n = 1, …, N_s, we define $yNs$ as

which has the following properties:

We now Taylor expand $G(yNs)$ near ⟨x⟩,

According to the central limit theorem, the expectation value of g = G(⟨x⟩) is then given by

The second term on the right-hand side of the above equation is what we refer to as bias, and the last term is the statistical error. Note that bias depends linearly on 1/N_s, whereas the statistical error depends on $1/Ns$⁠.