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

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