---
title: "Haque et al. develop highly parallel real-space RPA method in SPARC, Georgia Tech/LLNL; 8× reduction in computational prefactor."
sdDatePublished: "2026-09-16T04:09:00Z"
source: "https://arxiv.org/pdf/2609.16408"
topics:
  - name: "scientific research"
    identifier: "medtop:20000735"
  - name: "scientific innovation"
    identifier: "medtop:20000736"
  - name: "scientific publication"
    identifier: "medtop:20000740"
  - name: "technology and engineering"
    identifier: "medtop:20000756"
  - name: "information technology and computer science"
    identifier: "medtop:20000763"
locations:
  - "Atlanta"
  - "United States"
---


Haque et al. develop highly parallel real-space RPA method in SPARC, Georgia Tech/LLNL; 8× reduction in computational prefactor.

Highly Parallel Real-Space Random Phase Approximation Using Lanczos Quadrature and Interpolation

Highly Parallel Real-Space Random Phase
Approximation Using Lanczos Quadrature and
Interpolation
Abir Haque*1, Edmond Chow1, Shikhar Shah2, Andrew J. Medford3, John E.
Pask4, and Phanish Suryanarayana1,3
1College of Computing, Georgia Institute of Technology, Atlanta, GA, USA
2Department of Mathematics, Emory University, Atlanta, GA, USA
3College of Engineering, Georgia Institute of Technology, Atlanta, GA, USA
4Physics Division, Lawrence Livermore National Laboratory, Livermore, CA,
USA
*Email: abirhaque@gatech.edu
Abstract
We present a highly parallelizable, matrix-free, real-space method for computing the
random phase approximation (RPA) correlation energy within Kohn-Sham density func-
tional theory. In particular, we avoid the explicit construction and eigen-decomposition
of the response function matrix and use Lanczos quadrature to evaluate the trace of a
matrix function on a real-space grid. We also show it is possible to exploit the spatial
smoothness of the RPA correlation energy density in real-space to reduce computational
cost. Specifically, we compute the energy density on a coarse grid, then reconstruct an
approximation to the full, fine grid via interpolation. We implement this formulation
within the SPARC electronic structure package and demonstrate its convergence, ac-
curacy, agreement with planewave results, and scaling. Interpolation can enable an 8×
reduction in the computational prefactor. Given the embarrassingly parallel nature of
the proposed method, we achieve near-ideal speedups and near-cubic scaling, thus al-
lowing us to compute the RPA correlation energy for a silicon system with 256 valence
electrons at chemical accuracy in less than 30 minutes on 4,096 CPU cores.
1
Introduction
Electronic structure calculations, particularly those based on Kohn-Sham density functional
theory (DFT) [1, 2], are widely used in chemistry and materials science because they offer an
excellent balance between computational cost and accuracy. Nevertheless, the computational
cost of solving the generalized Kohn-Sham equations generally grows at least cubically with
1
arXiv:2609.16408v1 [physics.comp-ph] 14 Sep 2026

system size, limiting the range of systems that can be studied. This limitation is particularly
severe for nonlocal exchange-correlation approximations, whose evaluation requires informa-
tion beyond the electron density and its derivatives, thereby increasing both algorithmic
complexity and computational expense relative to local/semilocal functionals. Among such
approaches, the random phase approximation (RPA) is among the most accurate approxi-
mations to the exchange-correlation functional in condensed matter systems [3]. It includes
exact exchange and captures correlation through the spatially nonlocal, frequency-dependent
density response within the adiabatic-connection fluctuation-dissipation framework, rather
than relying on a local/semilocal approximation to the correlation hole. Consequently, RPA
can provide benchmark results due to its ability to accurately describe van der Waals inter-
actions, small-gap systems, and metals [4, 5].
Conventional RPA formulations, however, remain impractical for even moderate-sized
calculations because of their quartic scaling and substantial prefactor [6–8]. This cost arises
primarily from the explicit dependence on both occupied and unoccupied states: construc-
tion of the independent-particle density response function matrix requires all unoccupied
states, while evaluation of the correlation energy additionally entails numerical integration
over frequency. A variety of methods and software implementations have been developed to
perform RPA calculations [9–15], but each imposes modeling restrictions or computational
bottlenecks that limit its applicability. In particular, most of these approaches are formu-
lated within a planewave basis, which inherently restricts the boundary conditions to be
periodic [16], and which require additional approximations for isolated systems or partially
periodic systems. Among specific implementations, BerkeleyGW [17] and VASP [11] require
the precalculation of both occupied and unoccupied states. WEST reduces this requirement
to occupied orbitals alone, but still relies on global Fourier transforms, which have limited
scalability on distributed architectures because of the extensive inter-process communica-
tion required for data redistribution [18–21]. ABINIT instead computes the RPA correlation
energy through the conventional approach of eigendecomposing an explicitly constructed re-
sponse function matrix [6], thereby inheriting the same quartic-scaling limitations described
above. Overall, developing sub-quartic-scaling RPA formulations with a low prefactor, min-
imal global communication, and flexible boundary-conditions remains critical to simulating
large, complex chemical systems with chemical accuracy.
Unlike planewave-based approaches, real-space methods are not restricted to periodic
boundary conditions and can instead naturally accommodate Dirichlet, Neumann, and mixed
boundary conditions, allowing isolated, semi-infinite, and finite systems to be described
directly without artificial periodic replication or vacuum padding [22–30]. Such flexibility
is especially valuable for surfaces, nanostructures, and charged or polar systems, where
periodicity can introduce spurious electrostatic interactions between periodic images. Real-
space discretizations are also well suited to high-performance computing, as their inherent
locality supports efficient large-scale parallelism, in contrast to the dense operators and
global fast Fourier transform (FFT) communication of planewave methods.
Within the
context of RPA, Shah et al. [31] proposed a cubic-scaling real-space formulation that avoids
explicit construction of the response function matrix, employing a short-term recurrence
block Krylov subspace solver combined with subspace iteration and polynomial filtering.
Similarly, Zhang et al. [32] formulated a cubic-scaling real-space method that uses low-rank
approximations to avoid explicit construction of the response function matrix. Although
2

both methods achieve cubic scaling and support complex boundary conditions, they are
not massively parallelizable, as their eigensolvers rely on orthogonalization and subspace
projection steps that require frequent global communication. Perhaps more importantly,
both methods truncate the spectrum of the density response function for computational
efficiency, an approximation that can be inaccurate for systems whose spectra exhibit long
tails.
We present a highly parallelizable, matrix-free calculation method for the RPA correlation
energy, where Lanczos quadrature is utilized to directly evaluate the trace of the matrix
function on a real-space grid. Specifically, Lanczos quadrature computes individual diagonal
entries of the matrix function, where each entry represents the correlation energy density at
a given real-space grid point. This results in the method being embarrassingly parallel up to
the total number of real-space grid points. To significantly reduce the overall prefactor, we
exploit the spatial smoothness of these correlation energy densities. Within this framework,
Lanczos quadrature is applied only to a coarse set of grid points, and we subsequently recover
approximations for the full fine grid via interpolation. We also show that error cancellation
in the global sum of the interpolated correlation energy densities helps achieve chemical
accuracy in cases where point-wise errors may be high.
Furthermore, we show that the
proposed calculation method can achieve up to a 100× speedup compared to ABINIT, a
well-established planewave code, on comparable systems.
The remainder of this paper is organized as follows. First, we present the theoretical
foundations of computing the RPA correlation energy in real-space via Lanczos quadrature,
with and without interpolation of diagonal entries. Then, we present numerical experiments
that discuss the accuracy, scalability, and limitations of our method, which we implement
within SPARC, a real-space finite-difference electronic structure software package [23, 24].
Finally, we present our conclusions and future directions for our proposed framework for
RPA. Note that atomic units are used throughout unless otherwise specified.
2
Background
The RPA correlation energy is written as [33]:
Ec = 1
2π
Z ∞
0
Tr[log(I −νχ0(iω)) + νχ0(iω)] dω,
(1)
where ν = 4π(∇2)−1 is the Coulomb operator and χ0(iω) is the non-interacting Kohn-Sham
density response function at an imaginary frequency iω. Neglecting spin and Brillouin zone
integration, the response function matrix is explicitly written as [34, 35]:
χ0(r, r′, iω) = 2
Nd
X
j,k=1
(fj −fk)ψ∗
j(r)ψk(r)ψ∗
k(r′)ψj(r′)
εj −εk −iω
,
(2)
where r, r′ are real-space coordinates, Nd is the total number of finite difference points in
a given computational domain, ∗denotes complex conjugation, and fi are the occupations
with corresponding eigenpairs ψi, εi of the Hamiltonian H, satisfying:
Hψi = εiψi.
(3)
3

Henceforth, we will consider these operators as discretized within the real-space finite differ-
ence approach, wherein the above expressions remain unchanged. To avoid direct calculations
of χ0(iω), a matrix-free approach involving matrix-vector operations between χ0(iω) and a
vector u can be utilized. In particular, one may compute r = χ0(iω)u by utilizing the den-
sity functional perturbation theory (DFPT) framework [18, 36, 37]. This involves solving
the Sternheimer equation for Ns occupied states:
(H −εnI −iωI)∆ψn = −ψn ⊙u,
n = 1, . . . , Ns,
(4)
which contains a complex symmetric linear system with ⊙as the Hadamard product. Once
∆ψn for each occupied state is obtained, r is calculated by accumulating all corresponding
perturbations in electron density across all states.
3
Trace Approximations
3.1
Lanczos Quadrature
Our primary contribution lies in utilizing interpolation and Lanczos quadrature to approxi-
mate the trace of the function of a matrix. To achieve this, we first compute a subset of the
diagonal entries of the matrix function. While computing all elements of the diagonal of a
matrix via Lanczos quadrature, given an efficient matrix-vector product, is a well-established
technique within linear algebra and computational chemistry [31, 38–41], we show how it
can be used with interpolation. Alternative approaches to approximating the trace, such as
stochastic Lanczos quadrature [39], are available, but we do not utilize them in this work
due to their reliance on random sampling, which introduces statistical variance. Addition-
ally, low-rank approximations to the trace truncate the spectrum, which may be inaccurate
for spectra with long tails. In this section, we discuss the approximation techniques uti-
lized within the proposed method. Generally, one can directly approximate eT
j f(A)ej in the
expression
Tr[f(A)] =
n
X
j=1
eT
j f(A)ej
(5)
via Lanczos quadrature, where ej is the jth standard basis vector. Each term eT
j f(A)ej is
approximated via m steps of the Lanczos algorithm for symmetric A and starting vector ej,
where A is typically used in a matrix-free manner via an efficient matrix-vector multiplica-
tion. The Lanczos algorithm then returns a matrix Vm with orthonormal columns and an
m × m tridiagonal matrix Tm satisfying
AVm = VmTm.
(6)
This yields the approximation:
eT
j f(A)ej ≈[f(Tm)]1,1
(7)
where f(Tm) is quickly computed via an eigenvalue decomposition since m is small.
4

We are interested in the spectrum of νχ0(iω), a non-symmetric matrix. One may apply
a similarity transform on νχ0(iω) to obtain ν1/2χ0(iω)ν1/2, a symmetric matrix that shares
the same spectrum as νχ0(iω). This allows us to utilize the Lanczos method for quadrature,
which has significantly reduced computational complexity compared to the Arnoldi method
for no