Interpolative Separable Density Fitting through Centroidal Voronoi Tessellation with Applications to Hybrid Functional Electronic Structure Calculations

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.

Machine Learning K-Means Clustering of Interpolative Separable Density Fitting Algorithm for Accurate and Efficient Cubic-Scaling Exact Exchange Plus Random Phase Approximation within Plane Waves

The exact-exchange plus random-phase approximation (EXX+RPA) method has emerged as a crucial tool for precisely characterizing electronic structures in molecular and solid systems. We present an accurate and efficient implementation of EXX+RPA calculations that scale cubically and are conducted within plane waves. Our approach incorporates the interpolative separable density fitting (ISDF) algorithm, effectively mitigating the computational challenges associated with the plane wave basis set. To overcome the constraints of the conventional ISDF algorithm, characterized by the exceptionally high prefactor in QR factorization for interpolation point selection, we introduce an enhanced machine learning K-means method. This method incorporates a novel empirical weight function called "SSM+" for more precise interpolation point selection, capturing physical information more accurately across diverse systems. Our machine learning approach offers a quasiquadratic scaling alternative, effectively replacing the computationally demanding cubic-scaling QRCP algorithm in plane-wave-based EXX+RPA calculations. Furthermore, we enhance the method's capabilities by optimizing GPU acceleration using MATLAB's integrated GPU toolkit. In particular, our approach reduces the computational scaling of χ0 from 3.80 to 2.13 and the overall computational scaling of EXX from 2.74 to 2.10. We achieve a remarkable GPU acceleration speedup of up to 35×. Regarding CPU computation time, the standard quartic-scaling method requires 22 h to compute Si128, while QRCP completes the calculation in only around 1 h, achieving a speedup up to 20×. However, the utilization of the K-means algorithm reduces the time to 800 s, a substantial improvement of 100× compared to the standard algorithm. By employing the K-means algorithm, the computational time for interpolative point calculation using QRCP decreases from 1 h to 1 min, resulting in a 55× speed increase. With this improved algorithm, we successfully computed the dissociation curve of H2 and the equilibrium polyynic geometry of C18 molecules.

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.

Robust Tensor Hypercontraction of the Particle-Particle Ladder Term in Equation-of-Motion Coupled Cluster Theory

One method of representing a high-rank tensor as a (hyper-)product of lower-rank tensors is the tensor hypercontraction (THC) method of Hohenstein et al. This strategy has been found to be useful for reducing the polynomial scaling of coupled-cluster methods by representation of a four-dimensional tensor of electron-repulsion integrals in terms of five two-dimensional matrices. Pierce et al. have already shown that the application of a robust form of THC to the particle-particle ladder (PPL) term reduces the cost of this term in couple-cluster singles and doubles (CCSD) from O ( N 6 ) to O ( N 5 ) with negligible errors in energy with respect to the density-fitted variant. In this work, we have implemented the least-squares variant of THC (LS-THC) which does not require a nonlinear tensor factorization, including the robust form (R-LS-THC), for the calculation of the excitation and electron attachment energies using equation-of-motion coupled cluster methods EOMEE-CCSD and EOMEA-CCSD, respectively. We have benchmarked the effect of the R-LS-THC-PPL approximation on excitation energies using the comprehensive QUEST database and the accuracy of electron attachment energies using the NAB22 database. We find that errors on the order of 1 meV are achievable with a reduction in total calculation time of approximately 5 ×.

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.

High-Throughput Condensed-Phase Hybrid Density Functional Theory for Large-Scale Finite-Gap Systems: The SeA Approach

High-throughput electronic structure calculations (often performed using density functional theory (DFT)) play a central role in screening existing and novel materials, sampling potential energy surfaces, and generating data for machine learning applications. By including a fraction of exact exchange (EXX), hybrid functionals reduce the self-interaction error in semilocal DFT and furnish a more accurate description of the underlying electronic structure, albeit at a computational cost that often prohibits such high-throughput applications. To address this challenge, we have constructed a robust, accurate, and computationally efficient framework for high-throughput condensed-phase hybrid DFT and implemented this approach in the PWSCF module of Quantum ESPRESSO (QE). The resulting SeA approach (SeA = SCDM + exx + ACE) combines and seamlessly integrates: (i) the selected columns of the density matrix method (SCDM, a robust noniterative orbital localization scheme that sidesteps system-dependent optimization protocols), (ii) a recently extended version of exx (a black-box linear-scaling EXX algorithm that exploits sparsity between localized orbitals in real space when evaluating the action of the standard/full-rank V^xx operator), and (iii) adaptively compressed exchange (ACE, a low-rank V^xx approximation). In doing so, SeA harnesses three levels of computational savings: pair selection and domain truncation from SCDM + exx (which only considers spatially overlapping orbitals on orbital-pair-specific and system-size-independent domains) and low-rank V^xx approximation from ACE (which reduces the number of calls to SCDM + exx during the self-consistent field (SCF) procedure). Across a diverse set of 200 nonequilibrium (H₂O)₆₄ configurations (with densities spanning 0.4-1.7 g/cm³), SeA provides a 1-2 order-of-magnitude speedup in the overall time-to-solution, i.e., ≈8-26× compared to the convolution-based PWSCF(ACE) implementation in QE and ≈78-247× compared to the conventional PWSCF(Full) approach, and yields energies, ionic forces, and other properties with high fidelity. As a proof-of-principle high-throughput application, we trained a deep neural network (DNN) potential for ambient liquid water at the hybrid DFT level using SeA via an actively learned data set with ≈8,700 (H₂O)₆₄ configurations. Using an out-of-sample set of (H₂O)₅₁₂ configurations (at nonambient conditions), we confirmed the accuracy of this SeA-trained potential and showcased the capabilities of SeA by computing the ground-truth ionic forces in this challenging system containing >1,500 atoms.

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.

Fast Exchange with Gaussian Basis Set Using Robust Pseudospectral Method

In this article, we present an algorithm to efficiently evaluate the exchange matrix in periodic systems when a Gaussian basis set with pseudopotentials is used. The usual algorithm for evaluating exchange matrix scales cubically with the system size because one has to perform O(N²) fast Fourier transform (FFT). Here, we introduce an algorithm that retains the cubic scaling but reduces the prefactor significantly by eliminating the need to do FFTs during each exchange build. This is accomplished by representing the products of Gaussian basis function using a linear combination of an auxiliary basis the number of which scales linearly with the size of the system. We store the potential due to these auxiliary functions in memory, which allows us to obtain the exchange matrix without the need to do FFT, albeit at the cost of additional memory requirement. Although the basic idea of using auxiliary functions is not new, our algorithm is cheaper due to a combination of three ingredients: (a) we use a robust pseudospectral method that allows us to use a relatively small number of auxiliary basis to obtain high accuracy; (b) we use occ-RI exchange, which eliminates the need to construct the full exchange matrix; and (c) we use the (interpolative separable density fitting) ISDF algorithm to construct these auxiliary basis sets that are used in the robust pseudospectral method. The resulting algorithm is accurate, and we note that the error in the final energy decreases exponentially rapidly with the number of auxiliary functions.

Accelerating Linear-Response Time-Dependent Hybrid Density Functional Theory with Low-Rank Decomposition Techniques in the Plane-Wave Basis

We present an efficient low-rank implementation of linear-response time-dependent density functional theory for hybrid functionals (hybrid-LR-TDDFT) within the plane-wave pseudopotential framework. The adaptively compressed exchange (ACE) operator and the natural transition orbitals (NTOs) are introduced to build the low-rank representation of the nonlocal exchange operator in the hybrid-LR-TDDFT Hamiltonian. Numerical tests demonstrate that the ACE approximation significantly reduces the computational cost of applying the nonlocal exchange operator without loss of accuracy, and the NTO approximation can further accelerate the hybrid-LR-TDDFT calculations by introducing an NTO cutoff parameter. This new method enables us to effectively study the excitonic properties of two-dimensional MoS₂ consisting of 216 atoms and ∼1900 electrons with range-separated hybrid functionals on a single graphics processing unit.

Quantum algorithms for electronic structures: basis sets and boundary conditions

The advantages of quantum computers are believed to significantly change the research paradigm of chemical and materials sciences, where computational characterization and theoretical design play an increasingly important role. It is especially desirable to solve the electronic structure problem, a central problem in chemistry and materials science, efficiently and accurately with well-designed quantum algorithms. Various quantum electronic-structure algorithms have been proposed in the literature. In this article, we briefly review recent progress in this direction with a special emphasis on the basis sets and boundary conditions. Compared to classical electronic structure calculations, there are new considerations in choosing a basis set in quantum algorithms. For example, the effect of the basis set on the circuit complexity is very important in quantum algorithm design. Electronic structure calculations should be performed with an appropriate boundary condition. Simply using a wave function ansatz designed for molecular systems in a material system with a periodic boundary condition may lead to significant errors. Artificial boundary conditions can be used to partition a large system into smaller fragments to save quantum resources. The basis sets and boundary conditions are expected to play a crucial role in electronic structure calculations on future quantum computers, especially for realistic systems.

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.

Enabling Large-Scale Condensed-Phase Hybrid Density Functional Theory-Based Ab Initio Molecular Dynamics II: Extensions to the Isobaric-Isoenthalpic and Isobaric-Isothermal Ensembles

In the previous paper of this series [Ko, H.-Y. et al. J. Chem. Theory Comput. 2020, 16, 3757-3785], we presented a theoretical and algorithmic framework based on a localized representation of the occupied space that exploits the inherent sparsity in the real-space evaluation of the exact exchange (EXX) interaction in finite-gap systems. This was accompanied by a detailed description of exx, a massively parallel hybrid message-passing interface MPI/OpenMP implementation of this approach in Quantum ESPRESSO (QE) that enables linear scaling hybrid density functional theory (DFT)-based ab initio molecular dynamics (AIMD) in the microcanonical/canonical (NVE/NVT) ensembles of condensed-phase systems containing 500-1000 atoms (in fixed orthorhombic cells) with a wall time cost comparable to semi-local DFT. In this work, we extend the current capabilities of exx to enable hybrid DFT-based AIMD simulations of large-scale condensed-phase systems with general and fluctuating cells in the isobaric-isoenthalpic/isobaric-isothermal (NpH/NpT) ensembles. The theoretical extensions to this approach include an analytical derivation of the EXX contribution to the stress tensor for systems in general simulation cells with a computational complexity that scales linearly with system size. The corresponding algorithmic extensions to exx include optimized routines that (i) handle both static and fluctuating simulation cells with non-orthogonal lattice symmetries, (ii) solve Poisson's equation in general/non-orthogonal cells via an automated selection of the auxiliary grid directions in the Natan-Kronik representation of the discrete Laplacian operator, and (iii) evaluate the EXX contribution to the stress tensor. Using this approach, we perform a case study on a variety of condensed-phase systems (including liquid water, a benzene molecular crystal polymorph, and semi-conducting crystalline silicon) and demonstrate that the EXX contributions to the energy and stress tensor simultaneously converge with an appropriate choice of exx parameters. This is followed by a critical assessment of the computational performance of the extended exx module across several different high-performance computing architectures via case studies on (i) the computational complexity due to lattice symmetry during NpT simulations of three different ice polymorphs (i.e., ice Ih, II, and III) and (ii) the strong/weak parallel scaling during large-scale NpT simulations of liquid water. We demonstrate that the robust and highly scalable implementation of this approach in the extended exx module is capable of evaluating the EXX contribution to the stress tensor with negligible cost (<1%) as well as all other EXX-related quantities needed during NpT simulations of liquid water (with a very tight 150 Ry planewave cutoff) in ≈5.2 s ((H₂O)₁₂₈) and ≈6.8 s ((H₂O)₂₅₆) per AIMD step. As such, the extended exx module presented in this work brings us another step closer to routinely performing hybrid DFT-based AIMD simulations of sufficient duration for large-scale condensed-phase systems across a wide range of thermodynamic conditions.

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.

Cubic-Scaling All-Electron GW Calculations with a Separable Density-Fitting Space-Time Approach

We present an implementation of the GW space-time approach that allows cubic-scaling all-electron calculations with standard Gaussian basis sets without exploiting any localization or sparsity considerations. The independent-electron susceptibility is constructed in a time representation over a nonuniform distribution of real-space locations {r_k} optimized within a separable resolution-of-the-identity framework to reproduce standard Coulomb-fitting calculations with meV accuracy. The compactness of the obtained {r_k} distribution leads to a crossover with the standard Coulomb-fitting scheme for system sizes below a few hundred electrons. The needed analytic continuation follows a recent approach that requires the continuation of the screened Coulomb potential rather than the much more structured self-energy. The present scheme is benchmarked over large molecular sets, and scaling properties are demonstrated on a family of defected hexagonal boron-nitride flakes containing up to 6000 electrons.

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.

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

The high cost associated with the evaluation of Hartree-Fock exchange (HFX) makes hybrid functionals computationally challenging for large systems. In this work, we present an efficient way to accelerate HFX calculations with numerical atomic basis sets. Our approach is based on the recently proposed interpolative separable density fitting (ISDF) decomposition to construct a low-rank approximation of the HFX matrix, which avoids explicit calculations of the electron repulsion integrals (ERIs) and significantly reduces the computational cost. We implement the ISDF method for hybrid functional (PBE0) calculations in the HONPAS package. We take benzene and polycyclic aromatic hydrocarbon molecules as examples and demonstrate that hybrid functionals with ISDF yield quite promising results at a significantly reduced computational cost. Especially, the ISDF approach reduces the total cost of the evaluating HFX matrix by nearly 2 orders of magnitude compared to conventional approaches of direct evaluation of ERIs.

Enabling Large-Scale Condensed-Phase Hybrid Density Functional Theory Based Ab Initio Molecular Dynamics. 1. Theory, Algorithm, and Performance

By including a fraction of exact exchange (EXX), hybrid functionals reduce the self-interaction error in semilocal density functional theory (DFT) and thereby furnish a more accurate and reliable description of the underlying electronic structure in systems throughout biology, chemistry, physics, and materials science. However, the high computational cost associated with the evaluation of all required EXX quantities has limited the applicability of hybrid DFT in the treatment of large molecules and complex condensed-phase materials. To overcome this limitation, we describe a linear-scaling approach that utilizes a local representation of the occupied orbitals (e.g., maximally localized Wannier functions (MLWFs)) to exploit the sparsity in the real-space evaluation of the quantum mechanical exchange interaction in finite-gap systems. In this work, we present a detailed description of the theoretical and algorithmic advances required to perform MLWF-based ab initio molecular dynamics (AIMD) simulations of large-scale condensed-phase systems of interest at the hybrid DFT level. We focus our theoretical discussion on the integration of this approach into the framework of Car-Parrinello AIMD, and highlight the central role played by the MLWF-product potential (i.e., the solution of Poisson's equation for each corresponding MLWF-product density) in the evaluation of the EXX energy and wave function forces. We then provide a comprehensive description of the exx algorithm implemented in the open-source Quantum ESPRESSO program, which employs a hybrid MPI/OpenMP parallelization scheme to efficiently utilize the high-performance computing (HPC) resources available on current- and next-generation supercomputer architectures. This is followed by a critical assessment of the accuracy and parallel performance (e.g., strong and weak scaling) of this approach when AIMD simulations of liquid water are performed in the canonical (NVT) ensemble. With access to HPC resources, we demonstrate that exx enables hybrid DFT-based AIMD simulations of condensed-phase systems containing 500-1000 atoms (e.g., (H₂O)₂₅₆) with a wall time cost that is comparable to that of semilocal DFT. In doing so, exx takes us one step closer to routinely performing AIMD simulations of complex and large-scale condensed-phase systems for sufficiently long time scales at the hybrid DFT level of theory.

Accelerating Excitation Energy Computation in Molecules and Solids within Linear-Response Time-Dependent Density Functional Theory via Interpolative Separable Density Fitting Decomposition

We present an efficient way to compute the excitation energies in molecules and solids within linear-response time-dependent density functional theory (LR-TDDFT). Conventional methods to construct and diagonalize the LR-TDDFT Hamiltonian require ultrahigh computational cost, limiting its optoelectronic applications to small systems. Our new method is based on the interpolative separable density fitting (ISDF) decomposition combined with implicitly constructing and iteratively diagonalizing the LR-TDDFT Hamiltonian and only requires low computational cost to accelerate the LR-TDDFT calculations in the plane-wave basis sets under the periodic boundary condition. We show that this method accurately reproduces excitation energies in a fullerene (C₆₀) molecule and bulk silicon Si₆₄ system with significantly reduced computational cost compared to conventional direct and iterative calculations. The efficiency of this ISDF method enables us to investigate the excited-state properties of liquid water absorption on MoS₂ and phosphorene by using the LR-TDDFT calculations. Our computational results show that an aqueous environment has a weak effect on low excitation energies but a strong effect on high excitation energies of 2D semiconductors for photocatalytic water splitting.