Stochastic self-consistent second-order Green's function method for correlation energies of large electronic systems
Daniel Neuhauser Roi Baer Dominika Zgid
physics.chem-ph cond-mat.mes-hall cond-mat.mtrl-sci cond-mat.str-el physics.comp-ph


The second-order Matsubara Green's function method (GF2) is a robust
temperature dependent quantum chemistry approach, extending beyond
the random-phase approximation. However, till now the scope of GF2
applications was quite limited as they require computer resources
which rise steeply with system size. In each step of the self-consistent
GF2 calculation there are two parts: the estimation of the self-energy
from the previous step's Green's function, and updating the Green's
function from the self-energy. The first part formally scales as the
fifth power of the system size while the second has a much gentler
cubic scaling. Here, we develop a stochastic approach to GF2 (sGF2)
which reduces the fifth power scaling of the first step to merely
quadratic, leaving the overall sGF2 scaling as cubic. We apply the
method to linear hydrogen chains containing up to 1000 electrons,
showing that the approach is numerically stable, efficient and accurate.
The stochastic errors are very small, of the order of 0.1% or less
of the correlation energy for large systems, with only a moderate
computational effort. The first iteration of GF2 is an MP2 calculation
that is done in linear scaling, hence we obtain an extremely
fast stochastic MP2 (sMP2) method as a by-product. While here we consider
finite systems with large band gaps where at low temperatures effects
are negligible, the sGF2 formalism is temperature dependent and general
and can be applied to finite or periodic systems with small gaps at
finite temperatures.




§ INTRODUCTION


Second-order Green's function (GF2) is a temperature-dependent self-consistent
perturbation approach where the Green's function is iteratively renormalized.
At self-consistency the self-energy which accounts for the many-body
correlation effects is a functional of the Green's function, Σ(G).
The GF2 approximation as implemented here is described by the diagrams
in Fig. <ref> and employs Matsubara Green’s
functions that are temperature dependent and expressed on the imaginary
axis.<cit.> The implementation
we discuss, for total energies, relies on thermal Matsubara Green's
functions instead of real time Green's functions.<cit.>
This offers advantages in terms of stability and smoothness of the
self-energy. 

Upon convergence the GF2 method includes all second order skeleton
diagrams dressed with the renormalized second order Green's function
propagators, as illustrated in Fig. <ref>. Specifically,
as shown in Ref. , GF2, which at convergence
is reference independent, preserves the desirable features of M�ller-Plesset
perturbation theory (MP2) while avoiding the divergences that appear
when static correlation is important. Additionally, GF2 possesses
only a very small fractional charge and spin error,<cit.>
less than either typical hybrid density functionals or RPA with exchange,
therefore having a minimal many-body self-interaction error. In solids
GF2 describes the insulating and Mott regimes and recovers the internal
and free energy for multiple solid phases.<cit.>
Moreover, GF2 is useful for efficient Green's function embedding techniques
such as in the self-energy embedding method (SEET).<cit.>

The formal advantages of GF2 come, however, with a price tag. The
calculation of the self-energy matrix scales as O(n_τN^5),
where n_τ is the size of the imaginary time grid and N
the number of atomic orbitals (AOs). This leads to steep numerical
costs which prevent application of GF2 to systems larger than a few
dozen electrons. The application to larger systems requires therefore
a different paradigm and here we therefore develop a statistical formulation
of GF2 that calculates the self-energy matrix in linear-scaling.

The key to the present development, distinguishing it from previous
work <cit.>,
is the conversion of nested summations into stochastic averages. Our
method draws from previous work on stochastic electronic structure
methods, including stochastic- density functional theory (sDFT),<cit.>,
sDFT with long-range exact exchange,<cit.> multi-exciton
generation,<cit.> Moller-Plesset perturbation theory (sMP2),<cit.>
random-phase approximation (sRPA),<cit.> GW approximation
(sGW),<cit.> time-dependent
DFT (sTDDFT),<cit.> optimally-tuned range separated hybrid
DFT <cit.> and Bethe-Salpeter equation (sBSE).<cit.>
Among these, the closest to this work are the stochastic version of
sMP2 in real-time plane-waves,<cit.> and MO-based
MP2 with Gaussian basis sets.<cit.> The stochastic method
presented here benefits from the fact that the GF2 self-energy is
a smooth function of imaginary time and is therefore naturally amenable
to random sampling. 



§ METHOD




 §.§ Brief review of GF2


Our starting point is a basis of N real single-electron non-orthogonal
atomic-orbital (AO) states ϕ_i(𝐫), with
an N× N overlap matrix S_ij=⟨ϕ_i|ϕ_j.⟩.
Such states could be of any form, Gaussian, numerical, etc., but for
efficiency should be localized. We then use second quantization creation
a_i^† and annihilation a_i operators with respect
to the non-orthogonal basis ϕ_i(𝐫). The
non-orthogonality is manifested only in a modified commutation relation, 


    {a_i,a_j^+}=(S^-1)_ij.

The Hamiltonian for the interacting electrons has the usual form

    Ĥ=∑_ijh_ija_i^†a_j+1/2∑_ijklv_ijkla_i^†a_k^†a_la_j,

where h_ij=∫ drϕ_j(r)(-1/2∇^2+v_ext(r))ϕ_i(r)
and v_ext(r) is the bare external potential
(due to the nuclei), while V̂ is the two electron-electron
(e-e) Coulomb interaction described by the 2-electron integrals 

    v_ijkl=∬ϕ_i(r)ϕ_j(r)v(|r-r^'|)ϕ_k(r^')ϕ_l(r^')drdr^',

where v(r)=1/r is the Coulomb interaction potential. 

At a finite temperature β^-1 and chemical potential μ
we employ the grand canonical density operator e^-β(Ĥ-μN̂)/Z,
where N̂=∑_ijS_ija_i^†a_i is the electron-number
operator and Z(β)=[e^-β(Ĥ-μN̂)]
is the partition function. The thermal expectation value of any operator
Â can be calculated as ⟨Â⟩ =[e^-β((Ĥ-μN̂)/Z(β)Â].
For one-body observables Â=∑_ijA_ija_i^†a_j
we write ⟨Â⟩ =∑_ijA_ijP_ij
where P_ij=⟨ a_i^†a_j⟩ is
the reduced density matrix.

The 1-particle Green's function G_jk(τ) at an imaginary
time τ is a generalization of the concept of the density matrix
and obeys an equation of motion that can be solved by perturbation
methods. Formally:

    G_jk(τ)=-⟨ Ta_j(τ)a_k^†⟩ ,

where a_j(τ)≡ e^(Ĥ-μN̂)τa_je^-(Ĥ-μN̂)τ
with -β<τ<β, and T is the time-ordering symbol: 

    Ta_j(τ)a_k^†≡θ(τ)a_j(τ)a_k^†-θ(-τ)a_k^†a_j(τ).

Note that G(τ) is a real and symmetric matrix.

Each element G_jk(τ) (and therefore the entire
matrix G(τ)) is discontinuous when going from negative
to positive times, but this discontinuity is not a problem since we
only need to treat explicitly positive times τ>0 while negative
τ's are accessible by the anti-periodic relation for G(τ)

    G(τ)=-G(τ+β),   -β<τ<0,

as directly verified by substitution in Eq. (<ref>).
Hence G(τ) can be expanded as a Fourier series involving
the Matsubara frequencies ω_n=(2n+1)π/β
:

    G(τ)     =1/β∑_n=-∞^∞G(iω_n)e^-iω_nτ

where: 

    G(iω_n)=∫_0^βG(τ)e^iω_nτdτ.


The Green's function of Eq. (<ref>) gives access to
the reduced density matrix by taking the imaginary time τ as
a negative infinitesimal (denoted as 0^-):dxn@ucla.edu

    P_kj    =2G_kj(0^-)=-2G_kj(β^-)
        =2/β∑_n=0^∞e^-iω_n0^-G_kj(iω_n).

 Hence, all thermal averages of one-electron operators are accessible
through the sum of the Matsubara coefficients. 

Perturbation theory can be used to build approximations for G(τ)
based on a non-interacting Green's function G_0(τ)
corresponding to a reference one-body Hamiltonian Ĥ_̂0̂=∑_ijF_ija_i^†a_j.
Here, F is any real symmetric “Fock” matrix such that Ĥ_̂0̂
well approximates the interacting electron Hamiltonian. The derivation
of G_0(τ) requires orthogonal combination of the
basis set, i.e., finding a matrix X that fulfills XX^T=S^-1.
Then it is straightforward to show that 


    G_0(τ)     =
        Xe^-τ(F̅-μ)[θ(-τ)/1+e^β(F̅-μ)-θ(τ)/1+e^-β(F̅-μ)]X^T

where F̅=X^TFX is the Fock matrix in the orthogonal basis
set. Note that for positive (or negative) imaginary times G_0(τ)
is a real, smooth and non-oscillatory Green's function. This is important
for us since it much easier to stochastically sample a smooth function. 



Integration of Eqs. (<ref>) yields: 

    G_0(iω_n)=((μ+iω_n)S-F)^-1.

Since we now know how to write down Green's functions for non-interacting
systems, we rewrite the unknown part of the exact Green's function
by introducing the frequency-dependent self-energy, formally defined
by:

    G(iω_n)=((μ+iω_n)S-F-Σ(iω_n))^-1,

and by construction the self-energy fulfills the Dyson equation:

    G(iω_n)=G_0(iω_n)+G_0(iω_n)Σ(iω_n)G(iω_n).

Instead of viewing these equations as a definition of the self-energy
Σ(iω_n), we can calculate this self-energy
to a given order of perturbation theory in ΔĤ=Ĥ-Ĥ_0.
Specifically, the GF2 approximation<cit.> uses
a Hartree-Fock ansatz for F, 

    F_ij    =h_ij+1/2P_kl(2v_ijkl-v_ilkj),

where an Einstein summation convention is used, summing indices that
appear in pairs (here, both k and l). The self-energy in imaginary
time Σ(τ) is then obtained by second order perturbation
theory (see Fig. (<ref>)): 

    Σ_ij(τ)     =G_kl(τ)G_mn(τ)G_pq(β-τ)v_impk(2v_jnlq-v_jlnq).

Note that Σ(τ) and Σ(iω_n) are connected
by exactly the same Matsubara relations connecting G(τ) and
G(iω_n), Eqs. <ref>-<ref>.

The self-consistent one-body Green's function governs all one-body
expectation values. Moreover, even the total two-body potential energy
is available, by differentiation of the matrix trace (denoted by [])
of the Green's function with respect to τ: ⟨V̂⟩ =-1/2lim_τ→0^-[((∂/∂τ-μ)S+h)G(τ)].
Hence, the total energy is:

    ⟨Ĥ⟩ = Tr[hP-1/2lim_τ→0^-((∂/∂τ-μ)S+h)G(τ)].

It is easy to show by plugging the definition of G(τ) to Eq.
(<ref>) that this total energy has convenient
frequency and time forms:


    ⟨Ĥ⟩    =1/2[(h+F)P]+2/β∑_n[G(iω_n)Σ^T(iω_n)]
        =1/2[(h+F)P]+2∫_0^βTr[G(β-τ)Σ(τ)]dτ.


To conclude, the combination of Eqs. (<ref>), (<ref>),
(<ref>) and (<ref>) along with the requirement
that the density matrix describes N_e electrons results in the
following self-consistent GF2 procedure:


  * Perform a standard HF calculation and obtain a starting guess for
the Fock matrix F=F_HF and the density matrix P=P_HF. Set
Σ(iω_n)=0 for the set of N_ω positive
Matsubara frequencies ω_n, n=0,1,2,…,N_ω-1.,
.

  * Given Σ(iω_n)
and F, find μ such that Tr[PS]=N_e,
where P is given in Eq. (<ref>) from G(τ=β^-)
which depends on μ through the basic definition Eq. (<ref>).

  * Calculate G(τ) (Eq. (<ref>))
and P (Eq. (<ref>)).

  * Calculate the Fock matrix F from P (Eq. (<ref>)).

  * Calculate the self-energy Σ(τ)
from Eq. (<ref>) and transform is to the Matsubara
frequency domain to yield Σ(iω_n).

  * Calculate the total energy⟨Ĥ⟩ from
Eq. (<ref>).

  * Repeat steps 2-6 until convergence of the density and the total energy.

Once converged, the GF2 correlation energy is defined as the difference
E_ corr=⟨Ĥ⟩ -E_ HF
between the converged total energy (Eq. (<ref>)) and the
initial Hartree-Fock energy, E_HF=1/2[(h+F_HF)P_HF].
Note that in the first iteration GF2 yields automatically the temperature-dependent
MP2 energy: 


    E_ MP2^corr    =∫_0^βTr[G_0(β-τ;F_HF)Σ_0(τ)]dτ,

where Σ_0(τ) is that of Eq. (<ref>)
with G_0 replacing G. This expression reduces to the familiar
MP2 energy expression at the limit β→∞ (zero temperature
limit), when evaluated in the molecular orbital basis set that diagonalizes
the matrix F_HF. 

Finally, a technical point. The representation of the Green's functions
in τ-space can be complicated when the energy range of the eigenvalues
of F is large since a function of the type e^-τ(f-μ)/(1+e^-β(f-μ))
can be spiky when f>μ and τ→β or when f<μ and
τ→0. This requires special techniques for both imaginary time
and frequency grids as discussed in Refs. .



 §.§ sGF2: Stochastic approach to GF2


Most of the computational steps in the above algorithm scale with
system size N (number of AO basis functions) as O(N_SC× N_τ× N^3)
where N_SC is the number of GF2 self-consistent iterations and
N_τ is the number of time-steps. However, the main numerical
challenge in GF2 is step <ref> (Eq. (<ref>))
which scales formally as O(N_SC× N_τ× N^5)
making GF2 highly expensive for any reasonably sized system. This
steep scaling is due to the contraction of two 4-index tensors with
three Green's function matrices. 

To reduce this high complexity, we turn to the stochastic paradigm
which represents the matrices G(τ) by an equivalent
random average over stochastically chosen vectors. Fundamentally,
this is based on resolving the identity operator. Specifically, for
each τ we generate a vector η^0 of N components randomly
set to +1 or -1. Vectors at different times τ are statistically
independent, but we omit for simplicity their τ labeling. Then,
the key, and trivial, observation is that average of the product of
different components of η^0 is the unit matrix, which we write
symbolically as 


    η_k^0η_l^0=δ_kl.

We emphasize that the equality in this equation should be interpreted
to hold in the limit of averaging over infinitely many random vectors
η^0. 

Given this separable presentation of the unit matrix, it is easy to
rewrite any matrix as an average over separable vectors. Specifically,
from η^0 we define the two vectors: 


    η=√(|G(τ)|)η^0,     η̅=sgn(G(τ))√(|G(τ)|)η^0,

and then


    G_kl(τ)=η̅_kη_l.

Here, the square-root matrix is √(|G(τ)|)=A√(|g|)A^T,
where A(τ) is the unitary matrix of eigenvectors
and g(τ) is the diagonal matrix of eigenvalues of
G(τ). 

As a side note, we have a freedom to choose other vectors; specifically,
any two vectors η̅=D̅η^0,     η=Dη^0,
will work if D̅D^T=G(τ). In principle, we can even use
the simplest choice D̅=1,     D=G(τ), corresponding
to η̅=η^0 and η=G(τ)η^0. But while this
latter choice has the advantage that G(τ) does not need to be
diagonalized, we find that it is numerically better to use Eq. (<ref>)
as it is more balanced and therefore converges faster with the number
of stochastic samples. Also note that at the first iteration, where
G(τ)=G_0(τ), there is no need to diagonalize G_0(τ)
at different times, since it is obtained directly from the eigenstates
of F̅ in Eq. (<ref>).

Going back to Eq. (<ref>), we similarly
separate the other two Green's function matrices appearing in Eq. (<ref>),
writing them as G_mn(τ)=ξ̅_mξ_n and
G_pq(β-τ)=ζ̅_pζ_q. The self-energy
in Eq. (<ref>) is then

    Σ_ij(τ)     =η̅_kξ̅_mζ̅_pv_impk(2η_lξ_nζ_qv_jnlq-η_lξ_nζ_qv_jlnq),

so that it is separable to a product of two terms

    Σ_ij(τ)=u̅_i[2u_j-w_j],

where we defined three auxiliary vectors 

    u̅_i    =η̅_kξ̅_mζ̅_pv_impk
    
    u_j    =η_lξ_nζ_qv_jnlq
    
    w_i    =η_lξ_nζ_qv_jlnq.

The self-energy in Eq. (<ref>) should be viewed as
the average, over the stochastic vectors ξ^0, η^0 and
ζ^0, of the product term (u̅_i times 2u_j-w_j).

The direct calculation of the vectors u̅, u, w by Eq. (<ref>)
is numerically expensive once M>30 . We reduce the scaling by recalling
the definition of v_jnlq in Eq. (<ref>):

    u_j    =η_lξ_nζ_q∬ϕ_j(r)ϕ_n(r)v(|r-r^'|)ϕ_l(r^')ϕ_q(r^')drdr^'
        =∬ϕ_j(r)ξ(r)v(|r-r^'|)η(r^')ζ(r^')drdr^',

where:

    η(r)=η_lϕ_l(r),

and ξ(r) and ζ(r)
are analogously defined. We can therefore write


    u_j=∫ϕ_j(r)ξ(r)v_ηζ(r)dr,

where

    v_ηζ(r)≡∫ v(|r-r^'|)η(r^')ζ(r^')dr^'

is the Coulomb potential corresponding to the random charge distribution
η(r)ζ(r).
Similar expressions apply for u̅_i and w_j. 

Equations (<ref>)-(<ref>) are performed numerically
using FFT methods on a 3D Cartesian grid with N_g grid points,
so Eq. (<ref>) is calculated with O(N_glog N_g)
operations. Since the AO basis functions ϕ_i(r)
are local in 3D space, the calculations of η(r),
ξ(r) and ζ(r)
in Eq. (<ref>) scale linearly with system size. 

Eq. (<ref>) gives an exact expression for Σ_ij(τ),
as an expected value over formally an infinite number of stochastic
orbitals η^0, ξ^0 and ζ^0. Actual calculations
use a finite number I of “stochastic iterations”, where in
each such iteration a set of stochastic vectors η^0, ξ^0
and ζ^0 (different at each τ) is generated and Σ_ij(τ)
is averaged over them. The overall scaling of this step is therefore
I× N_τ×(N_glog N_g+N^2). We note
that the typical values of N_τ and I are in the hundreds,
see the discussion of the stochastic error below. 

Finally, we note that while the stochastic vectors (η^0,ζ^0,ξ^0)
are statistically independent for each time point τ, the same
τ-dependent vectors are used at each GF2 iteration, making it
possible to converge these iterations. 






§ RESULTS




 §.§ Systems and specifics


The algorithm was tested on linear hydrogen chains, (H_M) a nearest
neighbor distance of 1�, for several sizes: M=10,100,300 and
1000. The linearity was for convenience and we emphasize that it
does not play any role in the algorithm. The smallest chain was used
to demonstrate the convergence of the approach to the basis-set deterministic
values, and the other three calculations were used to study the dependence
of the algorithm on system size. 

In all calculations, an STO-3G basis was used, so that in this case
N=M and obviously the number of electrons is also N_e=M. A
periodic spatial grid of 0.5a_0 spacing was used to represent
the wave functions, and the grids contained 10×10 points in
the direction orthogonal to chain and between 60 and 4000 points along
the chain, depending on system size. For the smallest system (H_10)
a finer, bigger grid was also used, as detailed below. 

Other, technical details: 


  * Periodic images were screened using the method of Ref. .

  * The inverse temperature was β=50E_h^-1. 

  * A Chebyshev-type imaginary-time grid with 128 time points was employed
using a spline-fit method <cit.> for the
frequency-to-time conversions of G(iω_n) and Σ(iω_n)
and for the evaluation of the two-body energy.




 §.§ Small system


In our GF2 and MP2 algorithm, we make two types of numerical discretizations.
First, we use a finite number (labeled I) of stochastic iterations
to sample the self-energy, so we must show convergence as I grows.
Second, we use grids for bypassing the need to sum over O(N^4)
two-electron integrals, hence we need to demonstrate convergence with
respect to grid quality. We therefore examine in this section a small
system, linear H_10, and make four types of GF2/MP2 correlation
energy calculations: 


  * DET: fully deterministic calculations based on the analytical 2-electron
integrals; 

  * STOC(I)-NG: stochastic calculations based on I stochastic iterations
and on the analytical two electron integrals; 

  * STOC(I)-G1 and STOC(I)-G2: stochastic calculations based on
I stochastic iterations and on a 3D grid. Here, G1 is the same
type of grid we use for the larger calculations, and includes 10×10×60
points with a spacing h=0.5a_0. G2 is somewhat denser and covers
more space, with 16×16×100 points and h=0.4a_0. 

Our strategy is to first show that STOC-NG(I) converges to the
deterministic set (DET) as I grows. Then we show that for a given
number of stochastic orbitals, I=800, both grid results are quite
close to the non-grid result, and that the somewhat better second
grid (STOC(I=800)-G2 leads to extremely close results to the non-grid
values (STOC(I=800)-NG), so that the convergence with grid is very
rapid. 

We repeat the STOC-NG calculation 10 times determining the average
correlation energy E̅ and its standard deviation σ
as a function of I. The results are shown in the left panels of
Fig. <ref> as error-bars at E̅±σ,
which shrink approximately as 1/√(I) and which include the
DET result, represented as dashed horizontal lines, showing very small
or no bias. For MP2, a bias in the stochastic calculations is not
expected since the correlation energy is calculated linearly from
the first iteration of the self-energy Σ_0 (Eq. (<ref>)).
But for GF2 such a bias may form since the the “noisy” self-energy
is used non-linearly to update the Green's function in Eq. (<ref>).
However, for this small N=10 system the stochastic MP2 and GF2
energies do not exhibit a noticeable bias. We discuss the bias in
larger systems below. 

Next, we asses the errors associated with using grid calculations
replacing the analytical 2-electron integration. In both right panels
of Fig. <ref> we show 10 blue dots, each corresponding
to a pair of stochastic energies (E_STOC-G1(800),E_STOC-NG(800))_i
, i=1,…,10, both calculated with the same random seed s_i
(of course s_i and s_j are statistically independent). We
also show 10 orange dots, each corresponding to a pair of stochastic
energies (E_STOC-G2(800),E_STOC-NG(800))_i,
also calculated with the same seed s_i as before. The use of
the same seeds for each pair of blue and orange dots allows for comparison
of the grid error (which is the horizontal distance of a point from
the diagonal) without worrying about the larger statistical error,
seen as the spread of the results along the diagonal. We see that
the grid error decreases significantly when moving from G1 to G2,
but even the error for G1 is already very small (about 0.5meV per
electron). 



 §.§ Larger systems


In the small system considered above the bias was not noticeable and
here we examine the bias in larger systems. In Fig. <ref>
we show the STOC-G1(I) correlation energies in three specific systems
composed of N=100,300 and 1000 hydrogen atoms placed on a straight
line with a nearest neighbor spacing of 1 �. 

We first study the MP2 correlation energy of each system, appearing
in the lower energy range in the figure. The starting point of the
GF2 calculation is the Hartree-Fock F_HF and P_HF matrices,
so the MP2 energy is half the correlation energy of the first self-consistent
iteration (see Eqs. (<ref>) and (<ref>)). The statistical
errors in MP2 are pure fluctuations, a random number distributed normally
with zero average and with standard deviation given by σ_0/√(I)
where σ_0 is independent of I but shrinks with chain
length: σ_0∝1/√(L), exhibiting “self averaging”. <cit.>
The stochastic MP2 errors are very small and decrease with system
size, so for N=1000 the standard deviation of the I=800 iteration
calculation is 0.07% of the total correlation energy. For perspective,
note that (deterministic) errors of larger or similar magnitude are
present in linear scaling local or divide and conquer MP2 methods
with density fitting. <cit.>

Next, we discuss the stochastic estimates of the self-consistent GF2
correlation energies. These exhibit statistical errors with two visible
components. The first is a fluctuation, similar in nature to
that of the MP2 calculation, and the second component is a bias
which decreases as I grows. In fact, we expect the bias to asymptotically
decrease inversely with I, [A bias arises whenever we plug a random variable x, having an expected
value μ and variance σ^2, into a nonlinear function
f(x). One cannot hope that f(x) will have
the expected value of f(μ) unless f is a linear
function. A simple example is f(x)=x^2, where from
the definition of variance ⟨ f(x) ⟩ = f(μ)+σ^2
. Using the Taylor expansion of f around μ, it is straightforward
to show that f(x̅), where x̅=1/I∑_i=1^Ix_i
is an average over I samples and when I is sufficiently large,
⟨ f(x̅)⟩≈ f(μ)+f”(x̅)σ^2/2I
and so the bias is proportional to the variance of x, the curvature
of f at μ and inversely proportional to the number of iterations
I. ] so we fit the numerical GF2 results to a straight line in I^-1.
Table <ref> shows the estimate of the correlation energies,
the fluctuation and the bias as a function of the number of stochastic
orbitals. The results are highly accurate, for example when I=800
is used for the largest system (N=1000), the errors in MP2 and
in GF2 are smaller than 0.1%.





Timings. The measured overall CPU time for the stochastic
self-energy calculation (performed on a XEON system) can be expressed
as

    T^Σ≈2.5× N× N_τ× I×10^-7hr,

where, as mentioned, N is the number of electrons and I the
number of stochastic orbitals in the system. The MP2 wall time calculation
is essentially equal to the self-energy time divided by the number
of cores n_CORES, since the parallelization has negligible overhead:

    T_wall^MP2≈T^Σ/n_CORES.


GF2 involves an additional step, where the Green's function is constructed
from the self-energy and this step scales cubically with system size.
Furthermore, there are N_SC self-consistent iterations. The total
time is therefore found to be:

    T_wall^GF2   ≈ N_SC(1.7× N^3× N_τ1hr+T^Σ)/n_CORES.

For the H_1000 system, with I=800 stochastic orbitals,
the MP2 calculation takes T_wall^MP2=24hr/n_CORES,
i.e. about 30min when using 48 cores. 

The GF2 calculation for this same system involves N_SC=12 iterations
and a cubic part which takes about 2 core-hours per iteration, i.e.,
the cubic part is still an order of magnitude smaller than the self-energy
sampling time for this system size. The wall time is therefore T_wall^GF2=6.5hr
with 48 cores. 

For the H_300 system we find T_wall^MP2=12min
and T_wall^GF2=2hr while for H_100
we have T_wall^MP2=3.5min and T_wall^GF2=40 min. 

Note that these timings are for a single calculation. The error estimation
uses, as mentioned, ten completely independent runs, and therefore
took 10 times longer. 

For comparison, we note that the CPU time for the deterministic calculation
in the H_100 system takes 45 min. on a single core,
which is 4 times faster than the stochastic calculation. Since the
deterministic algorithm scales steeply as O(N^5),
the crossover occurs already at H_150 and at H_1000
the deterministic calculation would take 10^4-10^5/n_CORES
wall time hours per SCF iteration, compared to 24/n_CORES
hours for the stochastic calculations. 



 §.§ Born Oppenheimer potential curves


Potential energy curves can be calculated by correlated sampling,
where at each new nuclear configuration one employs the same set of
stochastic orbitals η_0,ξ_0 and ζ_0 for
the self-energy estimation. For demonstration, the HF, MP2 and GF2
Born Oppenheimer potentials of the H_100 system are shown
in Fig. <ref> as a function of the displacement of atom
no. 25 (counting from the left). In all three methods the most stable
position of the atom is at ∼-0.1a_0, slightly displaced towards
the nearest chain end. HF theory produces an energy potential with
large variations of up to 1.5eV and large vibrational frequencies
of order of 3.4eV. The MP2 curve is much smoother and the vibrational
frequencies reduces to ∼2.4eV while the GF2 energy curve is
considerably flatter, predicting a vibrational frequency of ∼1.0
eV. 






§ SUMMARY AND CONCLUSIONS


The problem we addressed here is the reduction of the the steep O( N^5)
scaling associated with the implementation of self-consistent GF2
calculations. We developed an effective way to reduce complexity to
O( N^3) by using stochastic techniques for calculating
the self-energy. A detailed derivation was given along with a specific
algorithm. The sampling error in the overall algorithm was studied
for linear H_N systems, and the simulation showed that
the stochastic errors in the correlation energies can be controlled
to less than 0.1% for very large systems. While the studied systems
were linear, the algorithm makes no use of the linearity and applies
equally well to any geometry.

As a byproduct, since the first step in GF2 is equivalent to MP2,
we obtain a stochastic MP2 method (sMP2) performed on top of an existing
HF calculation. This approach too has a formal complexity of O(N^5)
which is reduced here to linear O( N), except
for a single overall Fock-matrix diagonalization which is often available
from the underlying HF or DFT ground-state calculation. The errors
in this well-scaling stochastic MP2 method are comparable to those
of local MP2 approaches used in quantum chemistry. 

For GF2, the method has two main stages. The first stage, as in the
MP2 case, is a linear scaling calculation of the self-energy. This
self-energy is then used in the second stage to construct the Green's
function, at an O(N^3) cost. A complication
arises in GF2 due to this second stage (but not in MP2!), where the
self-energy enters non-linearly into the expression for the Green's
function. This non-linearity gives rise to a noticeable bias which
is proportional to the system size N. To overcome this bias the
number of stochastic orbitals I used in the first step must be
increased in proportion to the system size N, and hence the self-energy
calculation in GF2 attains an O(N^2) scaling.
The overall scaling of the GF2 calculation is unaffected by this bias
problem and remains O(N^3).

The present calculations give a fully self-consistent Green's function
method for a large system with a thousand electrons described by a
full quantum chemistry Hamiltonian. Moreover, we demonstrated that
the splitting of matrices by a random average over stochastically
chosen vectors leads to small variance and that relatively few Monte
Carlo samples already yield quite accurate correlation energies. The
reason for this excellent sampling dependence is two-fold: the stochastic
sampling inherently acts only in the space of atomic orbitals while
the actual spatial integrals (Eq. (<ref>)) are evaluated
using a deterministic, numerically exact calculation; in addition,
since the Green's function matrices are smooth in imaginary time,
different random vectors can be used at each imaginary-time point
thereby enhancing the stochastic sampling efficiency.

We have shown that both sMP2 and sGF2 are suitable for calculating
potential energy curves or surfaces. Interestingly, for the H_N
systems the potential curve is much smoother and flatter than in HF
or MP2. 

As for future applications, we note that sGF2 and sMP2 methods are
automatically suitable for periodic systems, as all the deterministic
steps and the time-frequency transforms are very efficient when done
in the reciprocal (k) space. The only additional detail is that
in periodic systems one needs to choose the random vectors to be in
k-space and then convert them to real-space, as detailed in an
upcoming article.

Finally, we also note that, beyond the results presented here, it
should also be possible to achieve further reduction of the stochastic
error with an embedded fragment approach, analogous to self-energy
embedding approaches, where a deterministic self-energy is calculated
for embedded saturated fragments as introduced for stochastic DFT
applications. <cit.>

Discussions with Eran Rabani are gratefully acknowledged. D.N. was
supported by NSF grant DMR-16111382, R.B. acknowledges support by
BSF grant 2015687 and D.Z. was supported by NSF grant CHE-1453894.





@ifundefinedendmcitethebibliography
  
38








f
subitem(mcitesubitemcount)

  
  

[Baym(1962)]baym1962self
Baym, G. Self-consistent approximations in many-body systems. Phys. Rev.
  1962, 127, 1391
[Dahlen and van Leeuwen(2005)Dahlen, and van Leeuwen]Dahlen05
Dahlen, N. E.; van Leeuwen, R. Self-consistent solution of the Dyson equation
  for atoms and molecules within a conserving approximation. The Journal
  of Chemical Physics 2005, 122
[Phillips and Zgid(2014)Phillips, and Zgid]Phillips2014
Phillips, J. J.; Zgid, D. Communication: The description of strong correlation
  within self-consistent Green's function second-order perturbation theory.
  J. Chem. Phys. 2014, 140, 241101
[Phillips et al.(2015)Phillips, Kananenka, and
  Zgid]Phillips2015
Phillips, J. J.; Kananenka, A. A.; Zgid, D. Fractional charge and spin errors
  in self-consistent Green's function theory. J. Chem. Phys.
  2015, 142, 194108
[Hedin(1965)]Hedin1965
Hedin, L. New Method for Calculating the One-Particle Green's Function with
  Application to the Electron-Gas Problem. Phys. Rev. 1965,
  139, A796–A823
[Fetter and Walecka(1971)Fetter, and Walecka]Fetter1971
Fetter, A. L.; Walecka, J. D. Quantum Thoery of Many Particle Systems;
  McGraw-Hill: New York, 1971; p 299
[Onida et al.(2002)Onida, Reining, and Rubio]Onida2002
Onida, G.; Reining, L.; Rubio, A. Electronic excitations: density-functional
  versus many-body Green's-function approaches. Rev. Mod. Phys.
  2002, 74, 601–659
[Rusakov and Zgid(2016)Rusakov, and Zgid]Rusakov2016
Rusakov, A. A.; Zgid, D. Self-consistent second-order Green's function
  perturbation theory for periodic systems. J. Chem. Phys.
  2016, 144, 054106
[Welden et al.(2016)Welden, Rusakov, and Zgid]Welden16
Welden, A. R.; Rusakov, A. A.; Zgid, D. Exploring connections between
  statistical mechanics and Green's functions for realistic systems:
  Temperature dependent electronic entropy and internal energy from a
  self-consistent second-order Green's function. The Journal of Chemical
  Physics 2016, 145, 204106
[Kananenka et al.(2015)Kananenka, Gull, and
  Zgid]Kananenka2015
Kananenka, A. A.; Gull, E.; Zgid, D. Systematically improvable multiscale
  solver for correlated electron systems. Phys. Rev. B 2015,
  91, 121111
[Lan et al.(2015)Lan, Kananenka, and Zgid]Tran15b
Lan, T. N.; Kananenka, A. A.; Zgid, D. Communication: Towards ab initio
  self-energy embedding theory in quantum chemistry. The Journal of
  Chemical Physics 2015, 143
[Nguyen Lan et al.(2016)Nguyen Lan, Kananenka, and
  Zgid]Tran16
Nguyen Lan, T.; Kananenka, A. A.; Zgid, D. Rigorous Ab Initio Quantum Embedding
  for Quantum Chemistry Using Green's Function Theory: Screened Interaction,
  Nonlocal Self-Energy Relaxation, Orbital Basis, and Chemical Accuracy.
  J. Chem. Theory Comput. 2016, 12, 4856–4870
[Lan and Zgid(2017)Lan, and Zgid]Tran17
Lan, T. N.; Zgid, D. Generalized Self-Energy Embedding Theory. The
  Journal of Physical Chemistry Letters 2017, 8, 2200–2205,
  PMID: 28453934
[Zgid and Gull(2017)Zgid, and Gull]zgid_njp17
Zgid, D.; Gull, E. Finite temperature quantum embedding theories for correlated
  systems. New J. Phys. 2017, 19, 023047
[Thom and Alavi(2007)Thom, and Alavi]Thom2007a
Thom, A. J.; Alavi, A. Stochastic perturbation theory: A low-scaling approach
  to correlated electronic energies. Phys. Rev. Lett. 2007,
  99, 143001
[Kozik et al.(2010)Kozik, Van Houcke, Gull, Pollet, Prokof'ev,
  Svistunov, and Troyer]kozik2010diagrammatic
Kozik, E.; Van Houcke, K.; Gull, E.; Pollet, L.; Prokof'ev, N.; Svistunov, B.;
  Troyer, M. Diagrammatic Monte Carlo for correlated fermions. EPL
  (Europhysics Letters) 2010, 90, 10004
[Willow et al.(2013)Willow, Kim, and
  Hirata]willow2013stochastic
Willow, S. Y.; Kim, K. S.; Hirata, S. Stochastic evaluation of second-order
  Dyson self-energies. The Journal of chemical physics 2013,
  138, 164111
[Baer et al.(2013)Baer, Neuhauser, and Rabani]Baer2013
Baer, R.; Neuhauser, D.; Rabani, E. Self-Averaging Stochastic Kohn-Sham
  Density-Functional Theory. Phys. Rev. Lett. 2013, 111,
  106402
[Neuhauser et al.(2014)Neuhauser, Baer, and
  Rabani]Neuhauser2014a
Neuhauser, D.; Baer, R.; Rabani, E. Communication: Embedded fragment stochastic
  density functional theory. J. Chem. Phys. 2014, 141,
  041102
[Arnon et al.(2017)Arnon, Rabani, Neuhauser, and
  Baer]arnon2017equilibrium
Arnon, E.; Rabani, E.; Neuhauser, D.; Baer, R. Equilibrium configurations of
  large nanostructures using the embedded saturated-fragments stochastic
  density functional theory. J. Chem. Phys. 2017, 146,
  224111
[Baer and Neuhauser(2012)Baer, and Neuhauser]Baer2012
Baer, R.; Neuhauser, D. Communication: Monte Carlo calculation of the
  exchange energy. J. Chem. Phys. 2012, 137,
  051103–4
[Baer and Rabani(2012)Baer, and Rabani]Baer2012a
Baer, R.; Rabani, E. Expeditious stochastic calculation of multiexciton
  generation rates in semiconductor nanocrystals. Nano Lett.
  2012, 12, 2123
[Neuhauser et al.(2013)Neuhauser, Rabani, and
  Baer]Neuhauser2013
Neuhauser, D.; Rabani, E.; Baer, R. Expeditious Stochastic Approach for MP2
  Energies in Large Electronic Systems. J. Chem. Theory Comput.
  2013, 9, 24–27
[Ge et al.(2013)Ge, Gao, Baer, Rabani, and Neuhauser]Ge2013
Ge, Q.; Gao, Y.; Baer, R.; Rabani, E.; Neuhauser, D. A Guided Stochastic
  Energy-Domain Formulation of the Second Order Møller–Plesset Perturbation
  Theory. J. Phys. Chem. Lett. 2013, 5, 185–189
[Takeshita et al.(2017)Takeshita, de Jong, Neuhauser, Baer, and
  Rabani]takeshita2017stochastic
Takeshita, T. Y.; de Jong, W. A.; Neuhauser, D.; Baer, R.; Rabani, E. A
  Stochastic Formulation of the Resolution of Identity: Application to Second
  Order M {\o} ller-Plesset Perturbation Theory. arXiv
  preprint arXiv:1704.02044 2017, 
[Neuhauser et al.(2013)Neuhauser, Rabani, and
  Baer]Neuhauser2013a
Neuhauser, D.; Rabani, E.; Baer, R. Expeditious Stochastic Calculation of
  Random-Phase Approximation Energies for Thousands of Electrons in Three
  Dimensions. J. Phys. Chem. Lett. 2013, 4,
  1172–1176
[Neuhauser et al.(2014)Neuhauser, Gao, Arntsen, Karshenas,
  Rabani, and Baer]Neuhauser2014
Neuhauser, D.; Gao, Y.; Arntsen, C.; Karshenas, C.; Rabani, E.; Baer, R.
  Breaking the Theoretical Scaling Limit for Predicting Quasiparticle Energies:
  The Stochastic G W Approach. Phys. Rev. Lett. 2014,
  113, 076402
[Vlcek et al.(2017)Vlcek, Rabani, Neuhauser, and
  Baer]Vlcek2017
Vlcek, V.; Rabani, E.; Neuhauser, D.; Baer, R. Stochastic GW calculations for
  molecules. arXiv preprint arXiv:1612.08999 2017, 
[Vlček et al.(2017)Vlček, Baer, Rabani, and
  Neuhauser]vlcek2017self
Vlček, V.; Baer, R.; Rabani, E.; Neuhauser, D. Self-consistent band-gap
  renormalization GW. arXiv preprint arXiv:1701.02023 2017,
  
[Gao et al.(2015)Gao, Neuhauser, Baer, and Rabani]Gao2015
Gao, Y.; Neuhauser, D.; Baer, R.; Rabani, E. Sublinear scaling for
  time-dependent stochastic density functional theory. J. Chem. Phys.
  2015, 142, 034106
[Neuhauser et al.(2015)Neuhauser, Rabani, Cytter, and
  Baer]Neuhauser2015
Neuhauser, D.; Rabani, E.; Cytter, Y.; Baer, R. Stochastic Optimally-Tuned
  Ranged-Separated Hybrid Density Functional Theory. J. Phys. Chem. A
  2015, 
[Rabani et al.(2015)Rabani, Baer, and Neuhauser]Rabani2015
Rabani, E.; Baer, R.; Neuhauser, D. Time-dependent stochastic Bethe-Salpeter
  approach. Phys. Rev. B 2015, 91, 235302
[Kananenka et al.(2016)Kananenka, Phillips, and
  Zgid]Kananenka2016
Kananenka, A. A.; Phillips, J. J.; Zgid, D. Efficient temperature-dependent
  Green's functions methods for realistic systems: Compact grids for orthogonal
  polynomial transforms. J. Chem. Theory Comput. 2016,
  12, 564–571
[Kananenka et al.(2016)Kananenka, Welden, Lan, Gull, and
  Zgid]Kananenka16
Kananenka, A. A.; Welden, A. R.; Lan, T. N.; Gull, E.; Zgid, D. Efficient
  Temperature-Dependent Green's Function Methods for Realistic Systems: Using
  Cubic Spline Interpolation to Approximate Matsubara Green's Functions.
  J. Chem. Theory Comput. 2016, 12, 2250–2259
[Martyna and Tuckerman(1999)Martyna, and Tuckerman]Martyna1999
Martyna, G. J.; Tuckerman, M. E. A reciprocal space based method for treating
  long range interactions in ab initio and force-field-based calculations in
  clusters. J. Chem. Phys. 1999, 110, 2810–2821
[Werner et al.(2003)Werner, Manby, and Knowles]Werner2003
Werner, H. J.; Manby, F. R.; Knowles, P. J. Fast linear scaling second-order
  Moller-Plesset perturbation theory (MP2) using local and density fitting
  approximations. J. Chem. Phys. 2003, 118,
  8149–8160
[Baudin et al.(2016)Baudin, Ettenhuber, Reine, Kristensen, and
  Kjærgaard]Baudin2016
Baudin, P.; Ettenhuber, P.; Reine, S.; Kristensen, K.; Kjærgaard, T.
  Efficient linear-scaling second-order Møller-Plesset perturbation theory:
  The divide–expand–consolidate RI-MP2 model. The Journal of chemical
  physics 2016, 144, 054102


