Accelerating Optical Absorption Spectra and Exciton Energy Computation for Nanosystems via Interpolative Separable Density Fitting
Wei Hu Meiyue Shao Andrea Cepellotti Felipe H. da Jornada Lin Lin Kyle Thicke Chao Yang Steven G. Louie
physics.comp-ph








Density fitting for GW and BSE calculations









Wei Hu et al.




^1Computational Research Division,
Lawrence Berkeley National Laboratory,

Berkeley, California 94720, United States



^2Department of Physics,
University of California, Berkeley

Berkeley, California 94720, United States



^3Department of Mathematics,
University of California, Berkeley

Berkeley, California 94720, United States



^4Department of Mathematics,
Duke University,

Durham, NC 27708, United States












Accelerating Optical Absorption Spectra and Exciton Energy Computation for Nanosystems via Interpolative Separable Density Fitting
    Wei Hu^1
Meiyue Shao^1 Andrea Cepellotti^2 Felipe
H. da Jornada^2 Lin Lin^3, 1 Kyle Thicke^4 Chao Yang^1 Steven G. Louie^2
    
==================================================================================================================================






We present an efficient way to solve the Bethe–Salpeter equation
(BSE), a model for the computation of absorption spectra in molecules and solids
that includes electron–hole excitations. Standard approaches
to construct and diagonalize the Bethe–Salpeter Hamiltonian require
at least (N_e^5) operations, where N_e is proportional to the number
of electrons in the system, limiting its application to small systems.
Our approach is based on the interpolative separable density fitting (ISDF) technique to
construct low rank approximations to the bare and screened exchange
operators associated with the BSE Hamiltonian. This approach
reduces the complexity of the Hamiltonian construction to
(N_e^3) with a much smaller pre-constant. Here, we implement the
ISDF method for the BSE calculations within the Tamm–Dancoff
approximation (TDA) in the BerkeleyGW software package. We show that
ISDF-based BSE calculations in molecules and solids reproduce
accurate exciton energies and optical absorption spectra with
significantly reduced computational cost.









§ INTRODUCTION
 


Many-Body Perturbation Theory <cit.> is a
powerful tool to describe one-particle and two-particle excitations,
and to obtain excitation energies and absorption spectra in
molecules and solids <cit.>. Within many-body
perturbation theory, Hedin's GW
approximation <cit.> has been successfully used
to compute quasi-particle (one-particle) excitation energies and the
Bethe–Salpeter equation (BSE) <cit.> describes
the interaction of a electron–hole pair (two-particle excitation)
produced during the spectral absorption in molecules and
solids <cit.>. Good agreement between theory and
experiment can only be achieved taking into account such
electron–hole interaction by solving the BSE. The BSE is an
eigenvalue problem. In the context of optical absorption, the
eigenvalues of the Bethe–Salpeter Hamiltonian (BSH) are related to
quasi-particle exciton energies and the corresponding eigenfunctions
yield the exciton wavefunctions.

The BSH has a special block structure to be shown in the next
section. It consists of bare and screened exchange kernels that
depend on the products of single-particle orbitals obtained from a
mean-field calculation. The evaluation of these kernels requires at
least (N_e^5) operations in a conventional approach, which is
very costly for large complex systems that contain hundreds or even
thousands of atoms. There are several methods, which have been
developed recently to generate a reduced basis set, to reduce such
high computational cost for the BSE
calculations <cit.>.





In this paper, we present a more efficient way to construct the BSH.
Such a construction allows the BSE to be solved efficiently by an
iterative diagonalization scheme. Our approach is based on the
recently developed interpolative separable density fitting (ISDF)
decomposition <cit.>. This ISDF decomposition has
been applied to accelerate a number of applications in computational
chemistry and materials science, including two-electron integral
computation <cit.>, correlation energy in the
random phase approximation <cit.>, density functional
perturbation theory <cit.>, and hybrid density
functional calculations <cit.>. Such a
decomposition is used to approximate the matrix consisting of
products of single-particle orbital pairs as the product of a matrix
consisting of a small number of numerical auxiliary basis vectors
and an expansion coefficient matrix <cit.>.
This low rank approximation effectively allows us to construct low
rank approximations to the bare and screened exchange kernels. The
construction of ISDF compressed BSE Hamiltonian matrix only requires
(N_e^3) operations if the rank of the numerical auxiliary
basis can be kept at (N_e) and if we keep the
bare and screened exchange kernel in a low rank factored form.
This is considerably more efficient than the (N_e^5)
complexity required in a conventional approach.





By keeping these kernels in the decomposed form, the matrix–vector
multiplications required in iterative diagonalization procedures for
computing the desired eigenvalues and eigenvectors of  can be
performed efficiently. We may also use these efficient matrix–vector
multiplications in a structure preserving Lanczos
algorithm <cit.> to obtain an approximate
absorption spectrum, without explicitly diagonalizing the
approximate .

We have implemented the ISDF based BSH construction in the
BerkeleyGW software package <cit.>. We
demonstrate that this approach can produce accurate exciton energies
and optical absorption spectra for molecules and solids, and reduce
the computational cost associated with the construction of the BSE
Hamiltonian significantly compared to the
algorithms used in the existing version of the BerkeleyGW software
package.





§ BETHE–SALPETER EQUATION





The Bethe–Salpeter equation is an eigenvalue problem of the form

    H_BSE X = X E,

where X is an exciton wavefunction, E is the corresponding exciton
energy and H_BSE is the Bethe–Salpeter Hamiltonian that
has the following block structure

    H_BSE =
      [  D + 2V_A - W_A      2V_B - W_B;     -2V_B + W_B - D -2V_A + W_A;                 ],








where D(i_vi_c,j_vj_c) = (ϵ_i_c -
ϵ_i_v)δ_i_vj_cδ_i_cj_c is an
(N_vN_c)×(N_vN_c) diagonal matrix with ϵ_i_v,
i_v=1,2,…,N_v being the quasi-particle energies associated
with valence bands and ϵ_i_c,
i_c=N_v+1,N_v+2,…,N_v+N_c being the quasi-particle energies
associated with conduction bands. These quasi-particle energies are
typically obtained from the so-called GW
calculation <cit.>. The V_A and V_B matrices
represent the (bare) exchange of electron–hole pairs, and the
W_A and W_B matrices are often referred to as the direct
terms that describe screened exchange of electron–hole pairs. These
matrices are defined as follows:

    V_A(i_vi_c,j_vj_c)    =∫_i_c(r)ψ_i_v(r)
    V(r,r')_j_v(r')ψ_j_c(r') dr dr',
    
    
    V_B(i_vi_c,j_vj_c)    =∫_i_c(r)ψ_i_v(r)
    V(r,r')_j_c(r')ψ_j_v(r') dr dr',
    
    
    W_A(i_vi_c,j_vj_c)
       =∫_i_c(r)ψ_j_c(r)W(r,r')
    _j_v(r')ψ_i_v(r')  dr dr', 
    
     W_B(i_vi_c,j_vj_c)
       =∫_i_c(r)ψ_j_v(r)W(r,r')
    _j_c(r')ψ_i_v(r')  dr dr',









where ψ_i_v and ψ_i_c are usually taken to be the
valence and conduction single-particle orbitals obtained from a
Kohn–Sham density functional theory (KSDFT) calculation
respectively, and V(r,r') and
W(r,r') are the bare and screened Coulomb operators.
Both V_A and W_A are Hermitian, whereas V_B and W_B are
complex symmetric in general. In the so-called Tamm–Dancoff
approximation (TDA) <cit.>, both V_B and W_B are
neglected and set to zeros in Equation (<ref>). As
a result,  is Hermitian in this case, and we only need to be
concerned with the upper left block of H_BSE.




Let M_cc(r) = {ψ_i_c_j_c},
M_vc(r) = {ψ_i_c_i_v}, and
M_vv(r) = {ψ_i_v_j_v} be matrices that
consist of products of discretized orbital pairs in real space, and
M̂_cc(G), M̂_vc(G),
M̂_vv(G) be the reciprocal space representation of
these matrices. Equations (<ref>) can then be succinctly
written as












    V_A    = M̂_vc^*V̂M̂_vc, 
    
    W_A    = (M̂_cc^*ŴM̂_vv),

where V̂ and Ŵ are reciprocal space representations
of the bare and screened exchange operators V and W,
respectively, and the reshape function is used to map the
(i_cj_c,i_vj_v)th element on the right-hand side of
(<ref>) to the (i_ci_v,j_cj_v)th element of
W_A.  A similar set of equations can be derived for V_B
and W_B. However, in this paper, we will mainly focus on the TDA
model.

The reason that the right-hand sides
of (<ref>) are computed in the reciprocal space is
that V̂ is diagonal and an energy cutoff is often adopted to limit
the number of the Fourier components used to represent ψ_i. As
a result, the leading dimension of M̂_cc, M̂_vc and
M̂_cc, which we denote by N_g, is often much smaller than
that of M_cc, M_vc and M_vv, which we denote by N_r.


















In additional to performing (N_e^2) Fast Fourier transforms
(FFTs) to obtain M̂_cc, M̂_vc and M̂_vv
from M_cc, M_vc and M_vv, respectively, we need to
perform at least

    (N_gN_c^2N_v^2 + N_g^2N_cN_v)

floating-point operations to obtain V_A and W_A using
matrix–matrix multiplication operations.

Note that, to achieve high accuracy with a large basis set, such as
the plane wave basis set, N_g is typically much larger than N_c
or N_v. The number of occupied bands is either N_e or N_e/2
depending on how spin is counted. In many existing calculations,
the actual number of conducting band N_c included in the
calculation is often a small multiple of N_v, whereas N_g
is often 100–10000 × N_e (N_r ∼ 10 × N_g). Therefore, the
second term in Equation (<ref>), which accounts for the cost
of multiplying Ŵ_A with Ẑ in
Equation (<ref>) can be much larger than other
parts under the Tamm–Dancoff approximation (TDA) in the BSE
calculations.





§ INTERPOLATIVE SEPARABLE DENSITY FITTING (ISDF) DECOMPOSITION


In order to reduce the computational complexity, we seek to minimize
the number of integrals we need to perform in
Equation (<ref>). This is possible if we can rewrite the
matrix M_ij, where the labels i and j are indices of
either valence or conducting orbitals,
as the product of a matrix Θ_ij that contains a set of
N_ij^t linearly independent auxiliary basis vectors with
N_ij^t ≈ tN_e ≪(N_e^2) (t is a small constant
referred as a rank truncation
parameter) <cit.> and an expansion
coefficient matrix C_ij. For large problems, the number of
columns of M_ij, which is either (N_vN_c), (N_v^2), or
(N_c^2), is typically larger than the number of grid points
N_r on which ψ_n(r) is sampled, i.e., the number of
rows in M_ij. As a result, N_ij^t should be much smaller than the
number of columns of M_ij. Even when a cutoff is used to limit
the size of N_c or N_v so that the number of columns in M_ij
is much less than N_g, we can still approximate M_ij by
Θ_ij C_ij with a Θ_ij that has a smaller rank
N_ij^t ∼ t√(N_iN_j).







To simplify our discussion, let us drop the subscript of M,
Θ and C for the moment, and describe the basic idea of
ISDF. The optimal low rank approximation of M can be obtained from
a singular value decomposition. However, the complexity of this
decomposition is at least (N_r^2 N_e^2) or (N_e^4). An
alternative decomposition, which is close to optimal, but has a much
favorable complexity has recently been developed. This type of
decomposition is called interpolative separable density fitting
(ISDF) <cit.>, which we will now describe.

In ISDF, instead of computing Θ and C simultaneously,
we fix the coefficient matrix C first, and
determine the auxiliary basis matrix Θ by solving a linear
least squares problem

    min M - Θ C _F^2,

where each column of M is given by ψ_i()_j()
sampled on a dense real space grids {_i}_i=1^N_r, and
Θ = [ζ_1, ζ_2, …, ζ_N^t] contains the
auxiliary basis vectors to be determined, ·_F denotes the
Frobenius norm.

We choose C to be a matrix that consists of
ψ_i()_j() evaluated at a subset of N^t
carefully chosen real space grid points, with N^t ≪ N_r and N^t
≪ N_e^2, i.e., each column of C indexed by (i,j) is given by

    [ ψ_i(1)_j(1), ⋯,
        ψ_i(k)_j(k), ⋯,
        ψ_i(N^t)_j(N^t)]^𝖳.


If the minimum of the objective function in Equation (<ref>) is
zero, the product of Θ and a column of C can be viewed as an
interpolation of corresponding function
{ψ_i(r)_j(r)} in M. However, in
general, we cannot expect the minimum of (<ref>) be
zero. The least squares minimizer is given by

    Θ = MC^* (CC^*)^-1.


It may appear that the matrix–matrix multiplications MC^* and
CC^* take (N_e^4) operations because the size of M is
N_r ×(N_e^2) and the size of C is N^t ×(N_e^2). However, both multiplications can be carried out with
fewer operations due to the separable structure of M and
C <cit.>. As a result, the computational
complexity for computing the interpolation vectors is
(N_e^3).




Intuitively, the least squares problem in
Equation (<ref>) is easier to solve when the rows of
C, which can be selected from the rows of M, are maximally
linearly independent. This task can be achieved by performing a QR
factorization of M^𝖳 with column pivoting
(QRCP) <cit.>. In QRCP, we choose a
permutation Π such that the factorization










    M^𝖳Π = QR

yields a unitary matrix Q and an upper triangular matrix R with
decreasing matrix elements along the diagonal of R. The magnitude
of each diagonal element R indicates how important the
corresponding column of the permuted M^𝖳 is, and whether the
corresponding grid point should be chosen as an interpolation point.
The QRCP decomposition can be terminated when the (N^t+1)-st
diagonal element of R becomes less than a predetermined threshold.
The leading N^t columns of the permuted M^𝖳 are considered to be
maximally linearly independent numerically. The corresponding grid
points are chosen as the interpolation points. The indices for the
chosen interpolation points N^t can be obtained from indices
of the nonzero entries of the first N^t columns of the permutation
matrix Π.



Roughly speaking, the QRCP moves matrix columns
of M^𝖳 with large norms to the left, and pushes matrix columns
of M^𝖳 with small norms to the right. Note that the square of the
vector 2-norm of the column of M^T associated with  is just

    ∑_i,j=1^N|ψ_i() _j()|^2 =
      (∑_i=1^N|ψ_i()|^2)
      (∑_j=1^N|ψ_j()|^2).

In the case when {ψ_i} and {ψ_j} both belong to the
set of occupied orbitals, the norm of each column of M^𝖳 is
simply the electron density. Hence the interpolation points chosen
by QRCP tend to locate in areas where the electron density is
relatively large. Once a column is selected, all other columns are
immediately orthogonalized with respect to the chosen column. Hence
nearly linearly dependent matrix columns will not be selected
repeatedly. As a result, the interpolation points chosen by QRCP are
well separated spatially. Notice that the standard QRCP procedure has
a high computational cost of (N_e^2N_r^2) ∼(N_e^4). But
it can be combined with the randomized sampling
method <cit.> so that its cost is reduced to
(N_rN_e^2) ∼(N_e^3).




§ LOW RANK REPRESENTATIONS OF BARE AND SCREENED EXCHANGE OPERATORS VIA ISDF








Applying the ISDF decomposition to M_cc, M_vc and M_vv
yields





























    M_cc   ≈Θ_cc C_cc,  
    
    M_vc   ≈Θ_vc C_vc, 
    
    M_vv   ≈Θ_vv C_vv.

It follows from Equations (<ref>), (<ref>)
and (<ref>) that the exchange and direct terms of the
BSE Hamiltonian can be written as

    V_A     = C_vc^*V_A C_vc, 
    
    W_A     = (C_cc^*W_A C_vv),

where V_A=Θ̂_vc^*V̂Θ̂_vc and
W_A = Θ̂_cc^*ŴΘ̂_vv are the projected
exchange and direct terms under the auxiliary
basis Θ̂_vc, Θ̂_cc and
Θ̂_vv. Here, Θ̂_vc, Θ̂_cc
and Θ̂_vv are reciprocal space representations
of Θ_vc, Θ_cc and Θ_vv, respectively, that
can be obtained via FFTs,

Note that the dimension of the matrix C_cc^*W_A
C_cc on the right-hand side of Equation (<ref>) is
N_c^2 × N_v^2.  It needs to be reshaped into a matrix of
dimension N_v N_c × N_vN_c according to the mapping
W_A(i_cj_c,i_vj_v) → W_A(i_vi_c,j_vj_c) before it can
combined with V_A matrix to construct the BSH.

Once the ISDF approximations for M_vc, M_cc and
M_vv are available, the cost for constructing a
low rank approximation to the exchange and direct terms reduced to
that of computing the projected exchange and direct kernels
Θ̂_vc^*V̂Θ̂_vc and
Θ̂_cc^*ŴΘ̂_vv,
respectively.  If the ranks of Θ_vc, Θ_cc
and Θ_vv are N_vc^t, N_cc^t and N_vv^t,
respectively, then the computational complexity for computing the
compressed exchange and direct kernels is
(N_vc^t N_vc^t N_g + N_cc^tN_vv^t N_g + N_vv^t N_g^2), which
is significantly lower than the complexity of the conventional
approach given in (<ref>).  When N_vc^t ∼ t√(N_vN_c), N_cc^t ∼ t√(N_cN_c) and
N_vv^t ∼ t√(N_vN_v) are on the order of N_e, the complexity of constructing
the compressed kernels is (N_e^3).





























































































§ ITERATIVE DIAGONALIZATION OF THE BSE HAMILTONIAN


In the conventional approach, exciton energies and wavefunctions
can be computed by using the recently developed BSEPACK
library <cit.>
to diagonalize the BSE Hamiltonian H_BSE. When TDA
is adopted, we may also just use a standard diagonalization
procedure implemented in the ScaLAPACK <cit.>
library.

When ISDF is used to construct low rank approximations to
the bare and screened exchange operators V_A and W_A, we
should keep both matrices in the factored form given
by Equation (<ref>). This is because that multiplying the
matrices on the right-hand sides of Equation (<ref>) would
require at least (N_e^5) operations, which is higher than the cost
for using the ISDF procedure to construct low rank
approximations to the BSH. Instead of using BSEPACK or ScaLAPACK
which has a complexity of (N_e^6), we propose to use
iterative methods to diagonalize the approximate BSH
constructed via the ISDF decomposition.



When TDA is used, several iterative methods such as the
Lanczos <cit.> and LOBPCG <cit.> algorithms can be used to compute a
few desired eigenvalues of the H_BSE. In each step of these
algorithms, we need to multiply H_BSE with a vector x of
size N_vN_c. When V_A is kept in the factored form given by
(<ref>), V_Ax can be evaluated as three matrix vector
multiplications performed in sequence, i.e.,

    V_A x ← C_vc^* [V_A (C_vc x)].

The complexity of these calculations is (N_vN_c N_vc^t).
If N_vc^t is on the order of N_e, then each V_A x can be
carried out in (N_e^3) operations.

Because C_cc^* W_A C_vv cannot be multiplied
with a vector x of size N_vN_c before it is reshaped,
a different multiplication scheme needs to be used. It follows
from the separable nature of C_vv and C_cc that this
multiplication can be succinctly described as

    W_A x = reshape[ Ψ_c^*(W⊙
    (Ψ_c X Ψ_v^*))Ψ_v
    ],

where X is a N_c × N_v matrix reshaped from the vector x,
Ψ_c is a N_cc^t × N_c matrix containing
ψ_i_c(r̂_k) as its elements, Ψ_v is a N_vv^t × N_v
matrix containing ψ_i_v(r̂_k) as its elements, and
⊙ denotes componentwise multiplication (Hadamard product).
The reshape function is used to turn the N_c× N_v matrix–matrix
product back into a size N_vN_c vector. If N_vv^t and N_cc^t
are on the order of N_e, then all matrix–matrix multiplications
in Equation (<ref>) can be carried out in (N_e^3) operations.
This make the complexity of each step of the iterative method
(N_e^3). If the number of iterative steps required to reach
convergence is not excessively large, then the ISDF enabled
iterative diagonalization can be carried out in (N_e^3) operations.













































§ ESTIMATING OPTICAL ABSORPTION SPECTRA WITHOUT DIAGONALIZATION

If we have all eigenpairs of H_BSE, we can easily obtain
the optical absorption spectrum, which is the imaginary part of the
dielectric function defined as


    ε_2(ω) = [
    8π^2e^2/Ω
    d_r^H((ω-iη)I-H_BSE)^-1d_l
    ],


where Ω is the volume of the primitive cell, e is the
elementary charge, d_r and d_l are the right and left optical
transition vectors, and η is a broadening factor used to
account for the lifetime of excitation.  However, it can become
prohibitively expensive to use an iterative diagonalization method
to compute all eigenpairs of H_BSE.

Fortunately, it is possible to use a structure preserving iterative
method to estimate the optical absorption spectrum without explicitly
computing all eigenpairs of H_BSE.
In Ref. <cit.>, we developed
a structure preserving Lanczos algorithm for estimating
the optical spectrum. The algorithm has been implemented in the
BSEPACK <cit.> library. When TDA is used,
the standard Lanczos algorithm can be used to estimate
the absorption spectrum.  When our objective is
to obtain the basic shape of the absorption spectrum to identify
where the major peaks are, it is not necessary to compute
all eigenpairs of H_BSE. As a result, the accuracy
required to construct approximate bare and screened exchange operators
used in BSE can possibly be lowered, thereby allowing us to
use a more aggressive truncation threshold in ISDF to further reduce
the cost of H_BSE construction. We will demonstrate
this possibility in the next section.



§ NUMERICAL RESULTS
 

In this section, we demonstrate the accuracy and efficiency of the ISDF
method when it is used to compute exciton energies and optical
absorption spectrum in the BSE framework.  We implemented the ISDF based BSH
construction in the BerkeleyGW software package <cit.>.
BerkeleyGW is a massively parallel computational package that uses
a many-body perturbation theory and Green's function formalism
to study quasi-particle excitation energies and optical absorption of
nanosystems. We use the ab initio software package
Quantum ESPRESSO (QE) <cit.> to compute the ground-state
quantities required in the GW and BSE calculations. In our Quantum ESPRESSO
density functional theory (DFT) based electronic structure calculations, we
use Hartwigsen–Goedecker–Hutter (HGH) norm-conserving
pseudopotentials <cit.> and the
LDA <cit.> exchange–correlation functional.
All the calculations were carried out on a single core at the Cori
[https://www.nersc.gov/systems/cori/]
systems at the National Energy Research Scientific Computing Center
(NERSC).

We performed calculations for three systems. They consist of a bulk
silicon Si_8 system and two molecular systems: carbon monoxide
(CO) and benzene (C_6H_6) as plotted in
Fig. <ref>. All systems are closed shell systems, and
the number of occupied bands is N_v = N_e/2, where N_e is the
valence electrons in the system.






 §.§ Accuracy


We first measure the accuracy of the ISDF method by comparing
the computed eigenvalues of the BSH matrices constructed with and without the ISDF
decomposition.

In our test, we set the plane wave energy cutoff required in the QE
calculations to E_cut = 10 Ha, which is relatively low.
However, this is sufficient for assessing the effect of ISDF. Such a
choice of E_cut results in N_r = 35937 and N_g = 2301
for the Si_8 system, N_r = 19683 and N_g = 1237 for the CO
molecule (N_v = 5), N_r = 91125 and N_g = 6235 for the benzene
molecule. We only include the lowest N_c conducting bands in the
BSE calculation. The number of active conduction bands (N_c) and
valence bands (N_v), the number of reciprocal grids and the
dimensions of the corresponding BSE Hamiltonian H_BSE for
these three systems are listed in Table <ref>.



In Fig. <ref>, we plot the singular values of the
matrices M_vc(r) = {ψ_i_c(r)_i_v(r)},
M_cc(r) = {ψ_i_c(r)_j_c(r)}
and M_vv(r) =
{ψ_i_v(r)_j_v(r)} associated with the
CO molecule. We observe that the singular values of these matrices
decay rapidly. For example, the leading 500 (out of 3600) singular
values of M_cc(r) decreases rapidly towards zero. All
other singular values are below 10^-4. Therefore, the
numerical rank N_cc^t of M_cc is roughly 500 (t = 8.3), or roughly 15% of the
number of columns in M_cc. Consequently, we expect that the rank
of Θ_cc produced in ISDF decomposition can be set to 15%
of N_c^2 without sacrificing the accuracy of the computed
eigenvalues.



This prediction is confirmed in Fig. <ref>, where
we plot the absolute difference between the lowest the exciton
energy of Si_8 computed with and without using ISDF to construct
H_BSE. To be specific, the error in the desired
eigenvalue is computed as ΔE = E_ISDF -
E_BGW, where E_ISDF is computed from the
H_BSE constructed with ISDF approximation, and
E_BGW is computed from a standard H_BSE
constructed without using ISDF. We first vary one of the ratios
N^t_cc/N_cc, N^t_vc/N_vc and N^t_vv/N_vv while
holding the others at a constant of 1. We observe that the error in
the lowest exciton energy (positive eigenvalue) is around 10^-3
Ha, when either N^t_cc/N_cc or N^t_vc/N_vc is set to
0.1 while the other ratios are held at 1. However, reducing
N^t_vv/N_vv to 0.1 introduces a significant amount of error
(0.1 Ha) in the lowest exciton energy. This is likely due to the
fact that N_v=16 is too small. We then hold N^t_vv/N_vv at
0.5 and let both N^t_cc/N_cc and N^t_vc/N_vc vary. The
variation of ΔE with respect to these ratios is also
plotted as in Fig. <ref>. We observe that the error in
the lowest exciton energy is still around 10^-3 Ha even when
both N^t_cc/N_cc and N^t_vc/N_vc are set to 0.1.





















We then check the absolute error ΔE (Ha) of all the exciton
energies computed with the ISDF method by comparing them with the ones obtained from a conventional BSE calculation implemented in BerkeleyGW
for the CO and benzene molecules. As we can see from  Fig. <ref>, the errors associated with these eigenvalues are all below 0.002 Ha when N^t_cc/N_cc is 0.1.






 §.§ Efficiency


At the moment, we have only implemented a sequential version of the
ISDF method within the BerkeleyGW software package. Therefore, our
efficiency test is limited in the size of the problem as well as the
number of conducting bands (N_c) we can include in the bare and
screened exchange operators. As a result, our performance
measurement does not fully reflect the computational complexity
analysis presented in the previous sections. In particular, taking
benzene as an example, N_g = 6235 is much larger than N_v = 15
and N_c = 60, therefore the computational cost of N_g^2N_v^2
∼(N_e^4) term is much higher than the N_gN_v^2N_c^2 ∼(N_e^5) term in the conventional BSE calculations.



Nonetheless, in this section, we will demonstrate the benefit of
using ISDF to reduce the cost for constructing the BSE Hamiltonian
H_BSE. In Table <ref>, we focus on the
benzene example and report the wall-clock time required to construct
the ISDF approximations of the M_vc,
M_cc, and M_vv matrices at different rank truncation
levels. Without using ISDF, it takes 746.0 seconds to
construct the reciprocal space representations of M_vc,
M_cc, and M_vv in BerkeleyGW. Most of the timing is spent
in the large number of FFTs applied to M_vc,
M_cc, and M_vv to obtain the reciprocal space representation of these matrices. We can clearly see that if only N^t_cc/N_cc is changed from 0.5 (t = 30.0) to 0.1 (t = 6.0), the wall-clock time used to construct a low rank approximation to M_cc can be reduced from 578.9 to 34.3 seconds. Furthermore, the total cost of computing M_vc, M_cc and M_vv can be further
reduced to around 1/19th that in a conventional approach
(39.3 vs. 746.0 seconds) if N^t_vc/N_vc,
N^t_vv/N_vv and N^t_cc/N_cc are all set to
0.1.















Note that because ISDF decomposition is carried out on a real space
grid, its measured cost reflect the cost for performing
QRCP in real space. Even though QRCP with random sampling has
(N_e^3) complexity, it has a relatively large pre-constant
compared to the size of the problem. Hence the measured
cost of QRCP is relatively high in this case. Recently,
Dong et al. proposed a new
approach to find the interpolation points based on the centroidal
Voronoi tessellation (CVT) method <cit.>, which
offers a much less expensive alternative to the QRCP procedure when
the ISDF method is used in hybrid functional calculations. We will
explore this CVT method in the BSE-ISDF calculations in the future.

In Table <ref>, we report the wall-clock time required
to construct the projected bare and screened exchange matrices
V_A and W_A that appear in
Equation (<ref>) once the ISDF approximations of
M_vc, M_vv, and M_cc become available.
Without ISDF, it takes
1.574 + 4.198 = 5.772 seconds to construct both W_A and V_A.
When N_cc^t/N_cc is set to 0.1 only,
the construction cost for W_A, which is the dominant cost, is reduced by
a factor of 2.8. Furthermore, if N^t_vc/N_vc, N^t_vv/N_vv
and N^t_cc/N_cc are all set to 0.1. We reduce the
cost for constructing V_A and W_A by
a factor of 63.0 and 10.1, respectively.
Note that the original implementation of the W_A and V_A in BerkeleyGW
is much slower because the elements of W_A and V_A are integrated one by
one. For benzene, it take takes 103,154 seconds (28.65 hours) to construct the BSE
Hamiltonian H_BSE in the original BerkeleyGW code.















































 §.§ Optical absorption spectra


One important application of BSE is to compute the optical
absorption spectrum, which is determined by optical dielectric
function in Equation (<ref>).
Fig. <ref> plots the optical absorption spectra
for both CO and benzene obtained from approximate H_BSE
constructed with the ISDF method and the H_BSE
constructed in a conventional approach implemented in BerkeleyGW.
When the rank truncation ratio N_cc^t/N_cc
is set to be only 0.10 (t = 6.0), the absorption
spectrum obtained from the ISDF approximate H_BSE
is nearly indistinguishable from that produced from the conventional
approach.  When N_cc^t/N_cc is set to
0.05 (t = 3.0), the absorption spectrum obtained from ISDF approximate H_BSE
still preserves the main features (peaks) of the absorption spectrum
obtained in a conventional approach even though some of the peaks are slightly shifted, and the height of some peaks are slightly off.






§ CONCLUSION AND OUTLOOK
 


In summary, we have demonstrated that the interpolative separable density
fitting (ISDF) technique can be used to efficiently and accurately
construct the Bethe–Salpeter Hamiltonian matrix. The ISDF method
allows us to reduce the complexity of the Hamiltonian construction
from (N_e^5) to (N_e^3) with a much smaller pre-constant.
We show that the ISDF based BSE calculations in molecules and solids
can efficiently produce accurate exciton energies and optical
absorption spectrum in molecules and solids.

In the future, we plan to replace the costly QRCP procedure 
with the centroidal Voronoi tessellation (CVT)
method <cit.> for selecting the interpolation points in
the ISDF method. The CVT method is expected to significantly reduce the
computational cost for selecting interpolating point in the ISDF procedure for the BSE calculations.

The performance results reported here are based on a sequential
implementation of the ISDF method. In the near future, we
will implement a parallel version suitable for large-scale
distributed memory parallel computers. Such an implementation
will allow us to tackle much larger problems for which
the favorable scaling of the ISDF approach is much more
pronounced.
















§ ACKNOWLEDGMENTS


This work is partly supported by 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. DE-AC02-05CH11231, as part of the Computational
Materials Sciences Program, and by the
Center for Applied Mathematics for Energy Research Applications
(CAMERA) (L. L. and C. Y.). The authors thank the National Energy
Research Scientific Computing (NERSC) center for making
computational resources available.

10

ScaLAPACK_2011
T. Auckenthaler, V. Blum, H. J. Bungartz, T. Huckle, R. Johanni, L. Krämer,
  B. Lang, H. Lederer, and P. R. Willems.
Parallel solution of partial symmetric eigenvalue problems from
  electronic structure calculations.
Parallel Comput., 37(12):783–794, 2011.

JCP_334_221_2017
P. Benner, S. Dolgov, V. Khoromskaia, and B. N. Khoromskij.
Fast iterative solution of the Bethe–Salpeter eigenvalue problem
  using low-rank and QTT tensor approximation.
J. Comput. Phys., 334:221–239, 2017.

JCTC_11_5197_2015
J. Brabec, L. Lin, M. Shao, N. Govind, Y. Saad, C. Yang, and E. G. Ng.
Efficient algorithms for estimating the absorption spectrum within
  linear response TDDFT.
J. Chem. Theory Comput., 11(11):5197–5208, 2015.

SIAM_13_727_1992_QRCP
T. F. Chan and P. C. Hansen.
Some applications of the rank revealing QR factorization.
SIAM J. Sci. Statist. Comput., 13:727–741, 1992.

CPC_183_1269_2012_BerkeleyGW
J. Deslippe, G. Samsonidze, D. A. Strubbe, M. Jain, M. L. Cohen, and S. G.
  Louie.
BerkeleyGW: A massively parallel computer package for the
  calculation of the quasiparticle and optical properties of materials and
  nanostructures.
Comput. Phys. Commun., 183(6):1269–1289, 2012.

arXiv_2017_CVT
K. Dong, W. Hu, and L. Lin.
Interpolative separable density fitting through centroidal voronoi
  tessellation with applications to hybrid functional electronic structure
  calculations.
arXiv:1711.01531, 2017.

JPCM_21_395502_2009_QE
P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni,
  D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso,
  S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann,
  C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri,
  R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia,
  S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M.
  Wentzcovitch.
QUANTUM ESPRESSO: A modular and open-source software project for
  quantum simulations of materials.
J. Phys.: Condens. Matter, 21(39):395502, 2009.

PRB_54_1703_1996_LDA
S. Goedecker, M. Teter, and J. Hutter.
Separable dual-space Gaussian pseudopotentials.
Phys. Rev. B, 54:1703, 1996.

PRB_58_3641_1998_HGH
C. Hartwigsen, S. Goedecker, and J. Hutter.
Relativistic separable dual-space gaussian pseudopotentials from H
  to Rn.
Phys. Rev. B, 58:3641, 1998.

PR_139_A796_1965_GW
L. Hedin.
New method for calculating the one-particle Green's function with
  application to the electron–gas problem.
Phys. Rev., 139:A796, 1965.

JCTC_13_5420_2017_ISDF
W. Hu, L. Lin, and C. Yang.
Interpolative separable density fitting decomposition for
  accelerating hybrid density functional calculations with applications to
  defects in silicon.
J. Chem. Theory Comput., 13(11):5420–5431, 2017.

MolPhys_114_1148_2016
P. B. V. Khoromskaia and B. N. Khoromskij.
A reduced basis approach for calculation of the BetheCSalpeter
  excitation energies by using low-rank tensor factorisations.
Mol. Phys., 114:1148–1161, 2016.

SIAMJSC_23_517_2001_LOBPCG
A. V. Knyazev.
Toward the optimal preconditioned eigensolver: Locally optimal block
  preconditioned conjugate gradient method.
SIAM J. Sci. Comput., 23(2):517–541, 2001.

JRNBS_45_255_1950_Lanczos
C. Lanczos.
An iteration method for the solution of the eigenvalue problem of
  linear differential and integral operators.
J. Res. Nat. Bur. Standards, 45:255–282, 1950.

LinXuYing2017
L. Lin, Z. Xu, and L. Ying.
Adaptively compressed polarizability operator for accelerating large
  scale ab initio phonon calculations.
Multiscale Model. Simul., 15:29–55, 2017.

PRB_92_075422_2015
M. P. Ljungberg, P. Koval, F. Ferrari, D. Foerster, and
  D. Sánchez-Portal.
Cubic-scaling iterative solution of the Bethe–Salpeter equation
  for finite systems.
Phys. Rev. B, 92:075422, 2015.

LuThicke2017
J. Lu and K. Thicke.
Cubic scaling algorithms for RPA correlation using interpolative
  separable density fitting.
J. Comput. Phys., 351:187–202, 2017.

JCP_302_329_2015
J. Lu and L. Ying.
Compression of the electron repulsion integral tensor in tensor
  hypercontraction format with cubic scaling cost.
J. Comput. Phys., 302:329–335, 2015.

PRB_95_075415_2017
M. Marsili, E. Mosconi, F. D. Angelis, and P. Umari.
Large-scale GW-BSE calculations with N^3 scaling: Excitonic
  effects in dye-sensitized solar cells.
Phys. Rev. B, 95:075415, 2017.

PR_46_618_1934_DMPT
C. Møler and M. S. Plesset.
Note on an approximation treatment for many-electron systems.
Phys. Rev., 46:618, 1934.

RMP_74_601_2002
G. Onida, L. Reining, and A. Rubio.
Electronic excitations: Density-functional versus many-body
  Green's-function approaches.
Rev. Mod. Phys., 74:601, 2002.

JCP_133_164109_2010
D. Rocca, D. Lu, and G. Galli.
Ab initio calculations of optical absorption spectra: Solution of the
  BetheCSalpeter equation within density matrix perturbation theory.
J. Chem. Phys., 133:164109, 2010.

PRB_62_4927_2000
M. Rohlfing and S. G. Louie.
Electron–hole excitations and optical spectra from first principles.
Phys. Rev. B, 62:4927, 2000.

PR_84_1232_1951_BSE
E. E. Salpeter and H. A. Bethe.
A relativistic equation for bound-state problems.
Phys. Rev., 84:1232, 1951.

arXiv_1611.02348_2016
M. Shao, F. H. da Jornada, L. Lin, C. Yang, J. Deslippe, and S. G. Louie.
A structure preserving Lanczos algorithm for computing the optical
  absorption spectrum.
SIAM J. Matrix. Anal. Appl., to appear.

LinearAlgebraAppl_488_148_2016
M. Shao, F. H. da Jornada, C. Yang, J. Deslippe, and S. G. Louie.
Structure preserving parallel algorithms for solving the
  BetheCSalpeter eigenvalue problem.
Linear Algebra Appl., 488:148–167, 2016.

BSEPACK_UserGuide_2016
M. Shao and C. Yang.
BSEPACK user's guide, 2016.
https://sites.google.com/a/lbl.gov/bsepack/.




