Interpolative Separable Density Fitting Decomposition for Accelerating Hartree-Fock Exchange Calculations within Numerical Atomic Orbitals

Use of Multigrids to Reduce the Cost of Performing Interpolative Separable Density Fitting

In this article, we present an interpolative separable density fitting (ISDF)-based algorithm to calculate the exact exchange in periodic mean field calculations. In the past, decomposing the two-electron integrals into the tensor hypercontraction (THC) form using ISDF was the most expensive step of the entire mean field calculation. Here, we show that by using a multigrid-ISDF algorithm, both the memory and the CPU cost of this step can be reduced. The CPU cost is brought down from cubic scaling to quadratic scaling with a low computational prefactor which reduces the cost by almost 2 orders of magnitude. Thus, in the new algorithm, the cost of performing ISDF is largely negligible compared to other steps. Along with the CPU cost, the memory cost of storing the factorized two-electron integrals is also reduced by a factor of up to 35. With the current algorithm, we can perform Hartree-Fock calculations on a diamond supercell containing more than 17,000 basis functions and more than 1500 electrons on a single node with no disk usage. For this calculation, the cost of constructing the exchange matrix is only a factor of 4 slower than the cost of diagonalizing the Fock matrix. Augmenting our approach with linear scaling algorithms can further speed up the calculations.

Sparse-Stochastic Fragmented Exchange for Large-Scale Hybrid Time-Dependent Density Functional Theory Calculations

We extend our recently developed sparse-stochastic fragmented exchange formalism for ground-state near-gap hybrid DFT to calculate absorption spectra within linear-response time-dependent generalized Kohn-Sham DFT (LR-GKS-TDDFT) for systems consisting of thousands of valence electrons within a grid-based/plane-wave representation. A mixed deterministic/fragmented-stochastic compression of the exchange kernel, here using long-range explicit exchange functionals, provides an efficient method for accurate optical spectra. Both real-time propagation as well as frequency-resolved Casida-equation-type approaches for spectra are presented, and the method is applied to large molecular dyes.

Machine Learning K-Means Clustering in Interpolative Separable Density Fitting Algorithm: Advancing Accurate and Efficient Cubic-Scaling Density Functional Perturbation Theory Calculations within Plane Waves

Density functional perturbation theory (DFPT) is a crucial tool for accurately describing lattice dynamics. The adaptively compressed polarizability (ACP) method reduces the computational complexity of DFPT calculations from O(N4) to O(N3) by combining the interpolative separable density fitting (ISDF) algorithm. However, the conventional QR factorization with column pivoting (QRCP) algorithm, used for selecting the interpolation points in ISDF, not only incurs a high cubic-scaling computational cost, O(N3), but also leads to suboptimal convergence. This convergence issue is particularly pronounced when considering the complex interplay between the external potential and atomic displacement in ACP-based DFPT calculations. Here, we present a machine learning K-means clustering algorithm to select the interpolation points in ISDF, which offers a more efficient quadratic-scaling O(N2) alternative to the computationally intensive cubic-scaling O(N3) QRCP algorithm. We implement this efficient K-means-based ISDF algorithm to accelerate plane-wave DFPT calculations in KSSOLV, which is a MATLAB toolbox for performing Kohn-Sham density functional theory calculations within plane waves. We demonstrate that this K-means algorithm not only offers comparable accuracy to QRCP in ISDF but also achieves better convergence for ACP-based DFPT calculations. In particular, K-means can remarkably reduce the computational cost of selecting the interpolation points by nearly 2 orders of magnitude compared to QRCP in ISDF.

Complex-Valued K-Means Clustering of Interpolative Separable Density Fitting Algorithm for Large-Scale Hybrid Functional Enabled Ab Initio Molecular Dynamics Simulations within Plane Waves

K-means clustering, as a classic unsupervised machine learning algorithm, is the key step to select the interpolation sampling points in interpolative separable density fitting (ISDF) decomposition for hybrid functional electronic structure calculations. Real-valued K-means clustering for accelerating the ISDF decomposition has been demonstrated for large-scale hybrid functional enabled ab initio molecular dynamics (hybrid AIMD) simulations within plane-wave basis sets where the Kohn-Sham orbitals are real-valued. However, it is unclear whether such K-means clustering works for complex-valued Kohn-Sham orbitals. Here, we propose an improved weight function defined as the sum of the square modulus of complex-valued Kohn-Sham orbitals in K-means clustering for hybrid AIMD simulations. Numerical results demonstrate that the K-means algorithm with a new weight function yields smoother and more delocalized interpolation sampling points, resulting in smoother energy potential, smaller energy drift, and longer time steps for hybrid AIMD simulations compared to the previous weight function used in the real-valued K-means algorithm. In particular, we find that this improved algorithm can obtain more accurate oxygen-oxygen radial distribution functions in liquid water molecules and a more accurate power spectrum in crystal silicon dioxide compared to the previous K-means algorithm. Finally, we describe a massively parallel implementation of this ISDF decomposition to accelerate large-scale complex-valued hybrid AIMD simulations containing thousands of atoms (2,744 atoms), which can scale up to 5,504 CPU cores on modern supercomputers.

Benchmarking the accuracy of the separable resolution of the identity approach for correlated methods in the numeric atom-centered orbitals framework

Four-center two-electron Coulomb integrals routinely appear in electronic structure algorithms. The resolution-of-the-identity (RI) is a popular technique to reduce the computational cost for the numerical evaluation of these integrals in localized basis-sets codes. Recently, Duchemin and Blase proposed a separable RI scheme [J. Chem. Phys. 150, 174120 (2019)], which preserves the accuracy of the standard global RI method with the Coulomb metric and permits the formulation of cubic-scaling random phase approximation (RPA) and GW approaches. Here, we present the implementation of a separable RI scheme within an all-electron numeric atom-centered orbital framework. We present comprehensive benchmark results using the Thiel and the GW100 test set. Our benchmarks include atomization energies from Hartree-Fock, second-order Møller-Plesset (MP2), coupled-cluster singles and doubles, RPA, and renormalized second-order perturbation theory, as well as quasiparticle energies from GW. We found that the separable RI approach reproduces RI-free HF calculations within 9 meV and MP2 calculations within 1 meV. We have confirmed that the separable RI error is independent of the system size by including disordered carbon clusters up to 116 atoms in our benchmarks.

Low-Scaling Algorithm for the Random Phase Approximation Using Tensor Hypercontraction with k-point Sampling

We present a low-scaling algorithm for the random phase approximation (RPA) with k-point sampling in the framework of tensor hypercontraction (THC) for electron repulsion integrals (ERIs). The THC factorization is obtained via a revised interpolative separable density fitting (ISDF) procedure with a momentum-dependent auxiliary basis for generic single-particle Bloch orbitals. Our formulation does not require preoptimized interpolating points or auxiliary bases, and the accuracy is systematically controlled by the number of interpolating points. The resulting RPA algorithm scales linearly with the number of k-points and cubically with the system size without any assumption on sparsity or locality of orbitals. The errors of ERIs and RPA energy show rapid convergence with respect to the size of the THC auxiliary basis, suggesting a promising and robust direction to construct efficient algorithms of higher order many-body perturbation theories for large-scale systems.

Even Faster Exact Exchange for Solids via Tensor Hypercontraction

Hybrid density functional theory (DFT) remains intractable for large periodic systems due to the demanding computational cost of exact exchange. We apply the tensor hypercontraction (THC) (or interpolative separable density fitting) approximation to periodic hybrid DFT calculations with Gaussian-type orbitals using the Gaussian plane wave approach. This is done to lower the computational scaling with respect to the number of basis functions (N) and k-points (Nk) at a fixed system size. Additionally, we propose an algorithm to fit only occupied orbital products via THC (i.e., a set of points, NISDF) to further reduce computation time and memory usage. This algorithm has linear scaling cost with k-points, no explicit dependence of NISDF on basis set size, and overall cubic scaling with unit cell size. Significant speedups and reduced memory usage may be obtained for moderately sized k-point meshes, with additional gains for large k-point meshes. Adequate accuracy can be obtained using THC-oo-K for self-consistent calculations. We perform illustrative hybrid density function theory calculations on the benzene crystal in the basis set and thermodynamic limits to highlight the utility of this algorithm.

Efficient implementation of analytical gradients for periodic hybrid functional calculations within fitted numerical atomic orbitals from NAO2GTO

The NAO2GTO scheme provides an efficient way to evaluate the electron repulsion integrals (ERIs) over numerical atomic orbitals (NAOs) with auxiliary Gaussian-type orbitals (GTOs). However, the NAO2GTO fitting will significantly impact the accuracy and convergence of hybrid functional calculations. To address this issue, here we propose to use the fitted orbitals as a new numerical basis to properly handle the mismatch between NAOs and fitted GTOs. We present an efficient and linear-scaling implementation of analytical gradients of Hartree-Fock exchange (HFX) energy for periodic HSE06 calculations with fitted NAOs in the HONPAS package. In our implementation, the ERIs and their derivatives for HFX matrix and forces are evaluated analytically with the auxiliary GTOs, while other terms are calculated using numerically discretized GTOs. Several integral screening techniques are employed to reduce the number of required ERI derivatives. We benchmark the accuracy and efficiency of our implementation and demonstrate that our results of lattice constants, bulk moduli, and band gaps of several typical semiconductors are in good agreement with the experimental values. We also show that the calculation of HFX forces based on a master-worker dynamic parallel scheme has a very high efficiency and scales linearly with respect to system size. Finally, we study the geometry optimization and polaron formation due to an excess electron in rutile TiO2 by means of HSE06 calculations to further validate the applicability of our implementation.

Interpolative Separable Density Fitting for Accelerating Two-Electron Integrals: A Theoretical Perspective

Low-rank approximations have long been considered an efficient way to accelerate electronic structure calculations associated with the evaluation of electron repulsion integrals (ERIs). As an accurate and efficient algorithm for compressing the ERI tensor, the interpolative separable density fitting (ISDF) decomposition has recently attracted great attention in this context. In this perspective, we introduce the ISDF decomposition from the theoretical aspects and technique details. The ISDF decomposition can construct a fully separable low-rank approximation (tensor hypercontraction factorization) of ERIs in real space with a cubic cost, offering great flexibility for accelerating high-scaling electronic structure calculations. We review the typical applications of ISDF in hybrid functionals, time-dependent density functional theory, and GW approximation. Finally, we discuss the promising directions for future development of ISDF.

Analytic gradients for local density fitting Hartree-Fock and Kohn-Sham methods

We present analytic gradients for local density fitting Hartree-Fock (HF) and hybrid Kohn-Sham (KS) density functional methods. Due to the non-variational nature of the local fitting algorithm, the method of Lagrange multipliers is used to avoid the solution of the coupled perturbed HF and KS equations. We propose efficient algorithms for the solution of the arising Z-vector equations and the gradient calculation that preserve the third-order scaling and low memory requirement of the original local fitting algorithm. In order to demonstrate the speed and accuracy of our implementation, gradient calculations and geometry optimizations are presented for various molecular systems. Our results show that significant speedups can be achieved compared to conventional density fitting calculations without sacrificing accuracy.

Reproducibility of Hybrid Density Functional Calculations for Equation-of-State Properties and Band Gaps

Hybrid density functional (HDF) approximations usually deliver higher accuracy than local and semilocal approximations to the exchange-correlation functional, but this comes with drastically increased computational cost. Practical implementations of HDFs inevitably involve numerical approximations─even more so than their local and semilocal counterparts due to the additional numerical complexity arising from treating the exact-exchange component. This raises the question regarding the reproducibility of the HDF results yielded by different implementations. In this work, we benchmark the numerical precision of four independent implementations of the popular Heyd-Scuseria-Ernzerhof (HSE) range-separated HDF on describing key materials' properties, including both properties derived from equations of state (EOS) and band gaps of 20 crystalline solids. We find that the energy band gaps obtained by the four codes agree with each other rather satisfactorily. However, for lattice constants and bulk moduli, the deviations between the results computed by different codes are of the same order of magnitude as the deviations between the computational and experimental results. On the one hand, this means that the HSE functional is rather accurate for describing the cohesive properties of simple insulating solids. On the other hand, this also suggests that the numerical precision achieved with current major HSE implementations is not sufficiently high to unambiguously assess the physical accuracy of HDFs. It is found that the pseudopotential treatment of the core electrons is a major factor that contributes to this uncertainty.

Low-Rank Approximations Accelerated Plane-Wave Hybrid Functional Calculations with k-Point Sampling

The low-rank approximations of the adaptively compressed exchange (ACE) operator and interpolative separable density fitting (ISDF) algorithms significantly reduce the computational cost and memory usage of hybrid functional calculations in real space, but the lack of k-point sampling hinders their implementation in reciprocal space for periodic systems with the plane-wave basis set. Here, we combine the ACE operator and ISDF decomposition into a new ACE-ISDF algorithm for periodic systems in reciprocal space with k-point sampling. On the basis of the ACE-ISDF algorithm with the improved reciprocal space ACE operator and k-point Fourier convolution, the time complexity of the hybrid functional calculation is reduced from O(Ne4Nk2) to O(Ne3Nklog(Nk)) (N_e and N_k are the number of electrons and k-points, respectively) with a much smaller prefactor and much lower memory consumption compared to the standard method for periodic systems with a plane-wave basis set.

Realizing Effective Cubic-Scaling Coulomb Hole Plus Screened Exchange Approximation in Periodic Systems via Interpolative Separable Density Fitting with a Plane-Wave Basis Set

The GW approximation is an effective way to accurately describe the single-electron excitations of molecules and the quasiparticle energies of solids. However, a perceived drawback of the GW calculations is their high computational cost and large memory usage, which limit their applications to large systems. Herein, we demonstrate an accurate and effective low-rank approximation to accelerate non-self-consistent GW (G₀W₀) calculations under the static Coulomb hole plus screened exchange (COHSEX) approximation for periodic systems. Our approach is to adopt the interpolative separable density fitting (ISDF) decomposition and Cauchy's integral to construct low-rank representations of the dielectric matrix ϵ and self-energy matrix Σ. This approach reduces the number of floating point operations from O(Ne4) to O(Ne3) and requires a much smaller memory footprint. Two methods are used to select the interpolation points in ISDF, including the standard QR factorization with column pivoting (QRCP) procedure and the machine learning K-means clustering (K-means) algorithm. We demonstrate that these two methods can yield similar accuracy for both molecules and solids at much lower computational cost. In particular, K-means clustering can significantly reduce the computational cost of selecting the interpolation points by an order of magnitude compared to QRCP, resulting in an overall speedup factor of about ten times ISDF accelerated the static COHSEX calculations compared with conventional COHSEX approximation.

Efficient Hybrid Density Functional Calculations for Large Periodic Systems Using Numerical Atomic Orbitals

We present an efficient, linear-scaling implementation for building the (screened) Hartree-Fock exchange (HFX) matrix for periodic systems within the framework of numerical atomic orbital (NAO) basis functions. Our implementation is based on the localized resolution of the identity approximation by which two-electron Coulomb repulsion integrals can be obtained by only computing two-center quantities-a feature that is highly beneficial to NAOs. By exploiting the locality of basis functions and efficient prescreening of the intermediate three- and two-index tensors, one can achieve a linear scaling of the computational cost for building the HFX matrix with respect to the system size. Our implementation is massively parallel, thanks to a MPI/OpenMP hybrid parallelization strategy for distributing the computational load and memory storage. All these factors add together to enable highly efficient hybrid functional calculations for large-scale periodic systems. In this work, we describe the key algorithms and implementation details for the HFX build as implemented in the ABACUS code package. The performance and scalability of our implementation with respect to the system size and the number of CPU cores are demonstrated for selected benchmark systems up to 4096 atoms.

Machine Learning K-Means Clustering Algorithm for Interpolative Separable Density Fitting to Accelerate Hybrid Functional Calculations with Numerical Atomic Orbitals

The interpolative separable density fitting (ISDF) is an efficient and accurate low-rank decomposition method to reduce the high computational cost and memory usage of the Hartree-Fock exchange (HFX) calculations with numerical atomic orbitals (NAOs). In this work, we present a machine learning K-means clustering algorithm to select the interpolation points in ISDF, which offers a much cheaper alternative to the expensive QR factorization with column pivoting (QRCP) procedure. We implement this K-means-based ISDF decomposition to accelerate hybrid functional calculations with NAOs in the HONPAS package. We demonstrate that this method can yield a similar accuracy for both molecules and solids at a much lower computational cost. In particular, K-means can remarkably reduce the computational cost of selecting the interpolation points by nearly two orders of magnitude compared to QRCP, resulting in a speedup of ∼10 times for ISDF-based HFX calculations.