#  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

Shizhe Jiao  Hefei National Research Center for Physical Sciences at the
Microscale, and Anhui Center for Applied Mathematics, University of Science
and Technology of China, Hefei, Anhui 230026, China  Jielan Li  Hefei National
Research Center for Physical Sciences at the Microscale, and Anhui Center for
Applied Mathematics, University of Science and Technology of China, Hefei,
Anhui 230026, China<EMAIL_ADDRESS>Xinming Qin  Hefei National
Research Center for Physical Sciences at the Microscale, and Anhui Center for
Applied Mathematics, University of Science and Technology of China, Hefei,
Anhui 230026, China  Lingyun Wan  Hefei National Research Center for Physical
Sciences at the Microscale, and Anhui Center for Applied Mathematics,
University of Science and Technology of China, Hefei, Anhui 230026, China  Wei
Hu  Hefei National Research Center for Physical Sciences at the Microscale,
and Anhui Center for Applied Mathematics, University of Science and Technology
of China, Hefei, Anhui 230026, China<EMAIL_ADDRESS>Jinlong Yang  Key
Laboratory of Precision and Intelligent Chemistry, and Department of Chemical
Physics, University of Science and Technology of China, Hefei, Anhui 230026,
China

###### Abstract

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. 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 apply the K-means clustering into
hybrid AIMD simulations for complex-valued Kohn-Sham orbitals and use an
improved weight function defined as the sum of the square modulus of complex-
valued Kohn-Sham orbitals in K-means clustering. Numerical results demonstrate
that this improved weight function in K-means clustering algorithm 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 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.

##  1 Introduction

The transposed Khatri-Rao product 1, 2 (also known as face-splitting product)
$Z=\\{z_{ij}:=\phi_{i}(\mathbf{r})\psi_{j}^{\ast}(\mathbf{r})\\}_{1\leq i\leq
N_{\phi},1\leq j\leq N_{\psi}}\in\mathbb{C}^{N_{r}\times(N_{\phi}N_{\psi})}$
of Kohn-Sham orbitals $\phi_{i}(\mathbf{r})$ and $\psi_{j}(\mathbf{r})$ in
real space $\\{\mathbf{r_{i}}\\}_{i=1}^{N_{r}}$ is inevitable for the multi-
center integrals of advanced electronic structure calculations in density
functional theory (DFT), 3, 4 especially for the Hartree-Fock (HF) 5, 6 and
post-HF electronic structure theory, such as time-dependent density functional
theory (TDDFT), 7, 8 GW approximation 9, 10, 11, 12 plus Bethe-Salpeter
equation (BSE), 13 Second-order Møller-Plesset perturbation theory (MP2), 14,
15, 16 and random phase approximation (RPA). 17 In order to reduce such high
computational cost and memory usage of such multi-center integrals in the
Kohn-Sham DFT calculations, several low rank approximation algorithms have
been proposed, such as the Cholesky decomposition, 18, 19 resolution of
identity (RI), 20, 21, 22, 23 tensor hypercontraction (THC) 24, 25, 26 and
pseudospectral decomposition. 27, 28 However, it is difficult to apply these
low rank approximation algorithms to the cases of atomic forces and
vibrational frequencies especially for periodic systems within plane-wave
basis sets, which can be used for a wide range of applications such as
geometry optimization and ab initio molecular dynamics (AIMD) simulation.

Recently, Lu et al. proposed a new tensor hypercontraction (THC) algorithm by
using the randomized QR factorization with column pivoting (QRCP) procedure
29, namely interpolative separable density fitting (ISDF), 30, 31 which can
achieve an effectively low rank approximation of the transposed Khatri-Rao
product of the Kohn-Sham orbitals ($\phi_{i}$ and $\psi_{j}$) and compress
their redundant information with cubic-scaling computational cost of
$O(N_{r}N_{\phi}N_{\psi})$. The transposed Khatri-Rao product of the Kohn-Sham
orbitals can be expressed by
$z_{ij}=\phi_{i}(\mathbf{r})\psi_{j}^{\ast}(\mathbf{r})\approx\sum_{\mu=1}^{N_{\mu}}\zeta_{\mu}(\mathbf{r})\phi_{i}(\mathbf{r_{\mu}})\psi_{j}^{\ast}(\mathbf{r_{\mu}})$,
where $\\{\mathbf{r}_{\mu}\\}_{\mu=1}^{N_{\mu}}$ are a set of interpolation
points from grid points $\\{\mathbf{r}_{i}\\}_{i=1}^{N_{r}}$ in real space,
$N_{\mu}$ is proportional to $\sqrt{N_{\phi}N_{\psi}}$
($N_{\mu}=t\sqrt{N_{\phi}N_{\psi}}$, $t$ is the rank truncation constant), and
$\zeta_{\mu}(\mathbf{r})$ is the auxiliary basis functions (ABFs). The ISDF
decomposition has already been applied successfully in several types of multi-
center integrals in the Kohn-Sham DFT calculations within Gaussian-type
orbitals (GTOs), numerical atomic orbitals (NAOs), and plane-wave (PW) basis
sets, such as hybrid DFT calculations, 31, 32, 33, 34 RPA correlation, 30
quantum Monte Carlo (QMC) simulations, 35 TDDFT, 36 MP2, 37 GW 38, 39 and BSE
40 calculations, for molecular and periodic systems.

The ISDF decomposition can be divided into two key steps, 31 including
selecting the interpolation points (IPs) and computing interpolation vectors
(IVs). The IVs can be computed easily by a least-squares fitting procedure
when the IPs are selected. The IPs mean the selection of a set of nonuniform
grid points $\\{\mathbf{r}_{\mu}\\}_{\mu=1}^{N_{\mu}}$ where the values of the
orbital pairs evaluated are almost consistent with that evaluated in all grid
points $\\{\mathbf{r}_{i}\\}_{i=1}^{N_{r}}$. Two approaches have been proposed
to select the IPs. The standard approach is the randomized QRCP as mentioned
previously, 29 which is accurate but expensive in the ISDF decomposition
process. 31 Another approach is the centroidal Voronoi tessellation (CVT)
algorithm proposed by Dong et al., 32 which only requires the information from
the electron density in the DFT calculations. The CVT method can be performed
easily by K-means clustering algorithm, a classical unsupervised machine
learning algorithm. Such K-means clustering algorithm aims at dividing $N_{r}$
data points into $N_{\mu}$ clusters. Each cluster includes all points whose
distances from the centroid of the cluster are smaller than that of other
clusters. Since K-means clustering only converges to a local optimal solution,
the accuracy of K-means clustering strongly depends on the selection of
initial centroids and the definition of distance and centroids. 41, 42, 32, 34
For real-valued orbitals, recent numerical results 31, 32, 33, 34 have
demonstrated that K-means clustering algorithms with relatively simple weight
definitions can yield reasonably high accuracy at a much lower cost.

However, the Kohn-Sham orbitals are usually represented as complex-valued
Bloch functions with k-point sampling 43, 44, 45 and noncollinear spin density
functional theory 46 for periodic systems and time-dependent wavefunctions 47,
48, 12, 49, 50 for real-time time-dependent density functional theory (RT-
TDDFT). For example, the transposed Khatri-Rao product of complex-valued Kohn-
Sham orbitals
$Z=\\{z_{ij}:=\sum_{k,l=1}^{N_{g}}\phi_{i,k}\psi_{j,l}^{\ast}e^{i(G_{k}-G_{l})\cdot\mathbf{r}}\\}$,
where $N_{g}$ is the number of plane-wave basis sets. $\phi_{i,k}$ and
$\psi_{j,l}^{\ast}$ are expansion coefficients within plane-wave basis sets
for $\phi_{i}$ and $\psi_{j}^{\ast}$, respectively. In this case, conventional
K-means clustering algorithms get into trouble 45, 46, 50 because the
centroids and distance cannot be well-defined for complex-valued points. An
alternative solution is to separate the complex data into real and imaginary
parts, then perform twice K-means clustering calculations and merge the real
and imaginary centroids together. This solution does not work for complex-
valued Kohn-Sham orbitals, since it is difficult to merge the real and
imaginary parts for their one-to-one correlation. To the best of our
knowledge, there are few works on complex-valued K-means clustering algorithms
yet. 51 In particular, Zhang et al. used complex-valued K-means to learn
filters in conventional neural networks for sleep stage classification. 52
They defined complex-valued centroids and inner product instead of Euclidean
distance. However, the complex number was converted from the real number in
this work and the real and imaginary parts are dependent on each other. The
approach is too expensive to be suitable for complex-valued Kohn-Sham orbitals
in the DFT calculations to deal with the transposed Khatri-Rao product
directly.

In the ISDF method, we introduce a weighted K-means clustering algorithm.
Because the real and imaginary parts of complex-valued Kohn-Sham orbitals are
discrete on the same grids in real space, we can perform K-means clustering on
the grids in real space and define different centroids from that in
conventional K-means clustering algorithm, which is defined as the weighted
average of all points belonging to the cluster. We convert the information of
complex-valued Kohn-Sham orbitals into the weight function. Therefore, it is
important to choose an appropriate weight function to compute the weighted
average for complex-valued Kohn-Sham orbitals.

In this work, we present an improved weight function desirable for complex-
valued Kohn-Sham orbitals in the DFT calculations. We apply successfully the
improved K-means clustering algorithm into complex-valued hybrid functional
calculations with plane-wave basis sets and achieve the acceleration of large-
scale hybrid density functional calculations containing thousands of atoms. In
particular, the ab initio molecular dynamics for complex-valued hybrid DFT
calculations in molecular and solid systems can be performed, such as liquid
water molecules 53, 53, 54, 55, 56, 57, 58, 32 and aluminium-silicon alloy,
59, 60 which is important but expensive within plane basis sets. 61, 62, 63,
64, 65, 66 We demonstrate that the energy potential calculated using this
improved weight function shows smoother, smaller energy drift and longer time
steps compared to the previous weight function in the K-means algorithm.
Therefore, this improved K-means clustering algorithm can accurately and
efficiently accelerate large-scale and long-time AIMD with complex-valued
hybrid DFT calculations.

This work is organized as follows. Section 2 gives a brief description of the
theoretical methodology, including the ISDF method, the complex-valued K-means
clustering algorithm in ISDF, the combination of hybrid DFT calculations with
ISDF method as well as their parallel implementation. Section 3 validates the
numerical accuracy and computational efficiency of the complex-valued K-means
clustering algorithm for ISDF decomposition to accelerate the hybrid DFT
calculations. A summary and outlook is given in Section 4.

##  2 Methodology

###  2.1 Interpolative separable density fitting

The ISDF decomposition is a new THC algorithm proposed by Lu and Ying firstly
29 and then promoted by Hu et al. for large-scale hybrid DFT calculations 31.
This algorithm can achieve low rank approximation of transposed Khatri-Rao
product (also known as face-splitting product)
$Z=\\{z_{ij}:=\phi_{i}(\mathbf{r})\psi_{j}^{\ast}(\mathbf{r})\\}_{1\leq i\leq
N_{\phi},1\leq j\leq N_{\psi}}\in\mathbb{C}^{N_{r}\times(N_{\phi}N_{\psi})}$
of Kohn-Sham orbitals $\phi_{i}$ and $\psi_{j}$.

$\phi_{i}(\mathbf{r})\psi_{j}^{\ast}(\mathbf{r})\approx\sum_{\mu=1}^{N_{\mu}}\zeta_{\mu}(\mathbf{r})C_{\mu}^{ij}$
(1)

where the transposed Khatri-Rao product can be approximately decomposed into
auxiliary basis functions (ABFs) $\zeta_{\mu}(\mathbf{r})$ and expansion
coefficients $C_{\mu}^{ij}$. The number of auxiliary basis functions
$N_{\mu}=t\sqrt{N_{\phi}N_{\psi}}$ can be regarded as the numerical rank of
the decomposition, where $t$ is a small parameter to achieve the compromise
between the numerical accuracy and computational efficiency. The key to ISDF
decomposition is to solve the expansion coefficients. The tensor
hypercontraction (THC) 67, 68, 26, 69 algorithms have made a success on it.
The ISDF algorithm provides a new tensor hypercontraction to obtain the
coefficients

$C_{\mu}^{ij}=\phi_{i}(\mathbf{r_{\mu}})\psi_{j}^{\ast}(\mathbf{r_{\mu}})$ (2)

where $\\{r_{\mu}\\}_{\mu=1}^{N_{\mu}}$ are a set of interpolation points from
grid points $\\{r_{i}\\}_{i=1}^{N_{r}}$ in real space. Therefore, the
transposed Khatri-Rao product can be expressed in the following form

$\phi_{i}(\mathbf{r})\psi_{j}^{\ast}(\mathbf{r})\approx\sum_{\mu=1}^{N_{\mu}}\zeta_{\mu}(\mathbf{r})\phi_{i}(\mathbf{r_{\mu}})\psi_{j}^{\ast}(\mathbf{r_{\mu}})$
(3)

Thus the ISDF method can be divided into two steps. 31 The first one is to get
the expansion coefficients $C_{\mu}^{ij}$, namely to compute the IPs, which
can be achieved by using the QRCP procedure or the weighted K-means clustering
algorithm. From the matrix point of view, IPs procedure is to select $N_{\mu}$
rows from $Z$ for fitting the entire matrix $Z$. Because Z matrix is not full
rank, the QRCP procedure can achieve the low rank decomposition of $Z$ as
follows,

$Z^{T}\Pi=Q\begin{bmatrix}R_{11}&R_{12}\\\
0&0\end{bmatrix}_{N_{\phi}N_{\psi}\times N_{r}}$ (4)

where $Z^{T}$ is the transpose of matrix $Z$, $\Pi$ is the permutation matrix
whose first $N_{\mu}$ columns give the IPs
$\\{\mathbf{r}_{\mu}\\}_{\mu=1}^{N_{\mu}}$. $Q$ and $R_{11}$ denote orthogonal
matrix and non-singular upper triangular matrix, respectively. The absolute
values of the diagonal entries for matrix $R_{11}$ follow a non-increasing
order. As the standard method, the computational cost of the QRCP procedure
scales as $O(N_{\mu}^{2}N_{r})$, while the cost of K-means clustering
algorithm is $O(N_{r}N_{\mu})$, as Table 1 shows.

The second step is to get the ABFs $\zeta_{\mu}(\mathbf{r})$, namely to
compute the IVs, which can be obtained by the least-squares procedure. After
we obtain the expansion coefficients $C_{\mu}^{ij}$, eq 1 can be written in
matrix form

$Z\approx\Theta C$ (5)

where $Z$ is $N_{r}\times N_{e}^{2}$ matrix ($N_{\phi}\approx N_{\psi}\sim
O(N_{e})$, $N_{e}$ is the number of electrons), which comes from
$\phi_{i}(\mathbf{r})\psi_{j}^{\ast}(\mathbf{r})$ sampled on a set of dense
real space grids $\\{\mathbf{r}_{i}\\}_{i=1}^{N_{r}}$.
$\Theta=[\zeta_{1},\zeta_{2},...,\zeta_{N_{\mu}}]$, $N_{r}\times N_{\mu}$
matrix, namely IVs.
$C=[\phi_{i}(\mathbf{r}_{1})\psi_{j}^{\ast}(\mathbf{r}_{1}),...,\phi_{i}(\mathbf{r}_{\mu})\psi_{j}^{\ast}(\mathbf{r}_{\mu}),...,\phi_{i}(\mathbf{r}_{N_{\mu}})\psi_{j}^{\ast}(\mathbf{r}_{N_{\mu}})]^{T}$.
Thus the IVs $\Theta$ can be given by

$\Theta=ZC^{T}(CC^{T})^{-1}$ (6)

where $ZC^{T}$ and $CC^{T}$ both need $O(N_{e}^{4})$ floating point
operations. Nevertheless, the separable structure of $Z$ and $C$ can reduce
the operations dramatically. 31 As is well known,

$\sum_{i,j}\phi_{i}\psi_{j}=(\sum_{i}\phi_{i})(\sum_{j}\psi_{j})$ (7)

Thus the $\mu$-th row, $\nu$-th column element of $ZC^{T}$ can be written as

$P^{\phi}(\mathbf{r}_{\mu},\mathbf{r}_{\nu})P^{\psi}(\mathbf{r}_{\mu},\mathbf{r}_{\nu})$
(8)

where $P^{\phi}(\mathbf{r}_{\mu},\mathbf{r}_{\nu})$ and
$P^{\psi}(\mathbf{r}_{\mu},\mathbf{r}_{\nu})$ are defined as

$\begin{split}P^{\phi}(\mathbf{r}_{\mu},\mathbf{r}_{\nu})&=\sum_{i}^{N_{\phi}}\phi_{i}(\mathbf{r}_{\mu})\phi_{i}^{\ast}(\mathbf{r}_{\nu})\\\
P^{\psi}(\mathbf{r}_{\mu},\mathbf{r}_{\nu})&=\sum_{i}^{N_{\psi}}\psi_{i}(\mathbf{r}_{\mu})\psi_{i}^{\ast}(\mathbf{r}_{\nu})\end{split}$
(9)

Here to compute $P^{\phi}$ and $P^{\psi}$ takes $O(N_{e}^{3})$ floating point
operations and the multiplication of $P^{\phi}$ and $P^{\psi}$ only needs
$O(N_{e}^{2})$ floating point operations. This conclusion also applies to
$CC^{T}$. Therefore, we can reduce the computational complexity of IVs from
$O(N_{e}^{4})$ to $O(N_{e}^{3})$.

Table 1: Computational cost and memory usage of IPs and IVs in ISDF decomposition. Notice that $N_{r}\approx 1,000\times N_{e}$, and $N_{\phi}\approx N_{\psi}\approx N_{\mu}\sim O(N_{e})$ in the plane-wave basis sets. Step | Algorithm | Computation | Memory
---|---|---|---
IPs | QRCP | $O(N_{\mu}^{2}N_{r})$ | $O(N_{r}N_{\mu})$
K-means | $O(N_{r}N_{\mu})$ | $O(N_{r}N_{\mu})$
IVs | Least-squares | $O(N_{r}N_{\mu}N_{e})$ | $O(N_{r}N_{\phi}N_{\psi})$

###  2.2 Complex-valued K-means clustering in ISDF

As an unsupervised machine learning algorithm, the K-means clustering
algorithm has been demonstrated to be much cheaper than the QRCP procedure.
32, 34 Ideally, conventional K-means clustering algorithm is for seeking the
solution to the following optimization problem.

$argmin\sum_{\mu=1}^{N_{\mu}}\sum_{\mathbf{r_{k}}\in
C_{\mu}}||Z(\mathbf{r}_{k})-Z(\mathbf{r}_{\mu})||^{2}$ (10)

Here we divide the $N_{r}$ data points into $N_{\mu}$ clusters
$\\{C_{\mu}\\}_{\mu=1}^{N_{\mu}}$. Each cluster can be denoted by its
centroid, namely the IPs we need.

${C_{\mu}}=\\{\ Z(\mathbf{r}_{i})\ |\
dist(Z(\mathbf{r}_{i}),Z(\mathbf{r}_{\mu}))\leq
dist(Z(\mathbf{r}_{i}),Z(\mathbf{r}_{m}))\ for\ all\ m\neq\mu\ \\}$ (11)

where the $dist(Z(\mathbf{r}_{i}),Z(\mathbf{r}_{\mu}))$ means distance between
data points $Z(\mathbf{r}_{i})$ and centroid $Z(\mathbf{r}_{\mu})$. The
optimization problem can only be solved by iterative calculations and K-means
clustering converges to a local minimum. The accuracy of K-means clustering
strongly depends on the selection of initial centroids and definition of
distance and centroids. 41, 42, 32, 34 The conventional K-means clustering is
mainly used for real-valued data points. As mentioned above, due to the
dependence of the real and imaginary parts, there is a dilemma when we apply
the conventional K-means clustering algorithm to the complex-valued Kohn-Sham
orbitals. In addition, it is quite expensive to perform direct K-means
clustering on the transposed Khatri-Rao product.

In the ISDF method, we effectively avoid these problems. We perform K-means
clustering on the grids in real space where the real and imaginary parts of
complex-valued Kohn-Sham orbitals map the same grids instead of the transposed
Khatri-Rao product. APPENDIX A verifies the feasibility of this strategy. In
ISDF method, the optimization problem is reduced to

$argmin\sum_{\mu=1}^{N_{\mu}}\sum_{\mathbf{r_{k}}\in
C_{\mu}}w(\mathbf{r}_{k})||\mathbf{r}_{k}-\mathbf{r}_{\mu}||^{2}$ (12)

where the $w(\mathbf{r}_{k})$ is the weight function. The distances between
the grid points are calculated by the Euclidean metric. The centroids
$\mathbf{r}_{\mu}$ can be defined by the weighted average of all points that
belong to the corresponding cluster as follows

$\mathbf{r}_{\mu}=\frac{\sum_{\mathbf{r}_{j}\in
C_{\mu}}\mathbf{r}_{j}w(\mathbf{r}_{j})}{\sum_{\mathbf{r}_{j}\in
C_{\mu}}w(\mathbf{r}_{j})}$ (13)

The information of complex-valued Kohn-Sham orbitals is converted into the
weight function. Therefore, it is very important to select and define the
empirical weight function. Different empirical weight functions have been
proposed for different data types in practice. 32, 34

For hybrid functional electronic structure calculations within numerical
atomic orbitals (NAOs), we have proposed the norm of the row of orbital pairs
as the weight function, 34 namely

$w(\mathbf{r})=\sum_{i,j=1}^{N_{b}}|\varphi_{i}(\mathbf{r})||\varphi_{j}(\mathbf{r})|$
(14)

where $\\{\varphi_{i}\\}_{i=1}^{N_{b}}$ denote the real-valued NAOs and
$N_{b}$ denotes the number of NAOs. It should be noticed that such two sets of
real-valued orbitals involved in the transposed Khatri-Rao product are the
same.

For real-valued hybrid DFT calculations within plane-wave basis sets,32 we
have defined the weight function as

$w(\mathbf{r})=\sum_{i,j=1}^{N_{e}}|\phi_{i}(\mathbf{r})|^{2}|\psi_{j}(\mathbf{r})|^{2}=(\sum_{i=1}^{N_{e}}|\phi_{i}|^{2})(\sum_{j=1}^{N}|\psi_{j}|^{2})$
(15)

The weight function is the product of the square modulus (abbreviated as PSM)
of real-valued Kohn-Sham orbitals $\\{\phi_{i}\\}_{i=1}^{N_{e}}$ and
$\\{\psi_{j}\\}_{j=1}^{N}$, where $N$ is the number of Kohn-Sham orbitals. It
should be noticed that such two sets of real-valued Kohn-Sham orbitals
involved in the transposed Khatri-Rao product are different because a two-
level self-consistent field (SCF) iteration procedure 70 is used in hybrid DFT
calculations within plane-wave basis sets.

However, the PSM weight function defined in Eq.15 is prone to more zero
elements, which makes the sampling points more localized in the real system as
shown in FIG. 11 (e, f), and the selected points are concentrated around the
ball and stick model. Because the rows of Z selected at the IPs are linearly
independent, the IPs should be delocalized as much as possible. Therefore, in
the case of complex-valued Kohn-Sham orbitals in the DFT calculations with
plane-wave basis sets, we use an improved weighted function for the transposed
Khatri-Rao product of two different sets of complex-valued Kohn-Sham orbitals
defined as

$w(\mathbf{r})=(\sum_{i=1}^{N_{e}}|\phi_{i}|^{\alpha})+(\sum_{j=1}^{N}|\psi_{j}|^{\alpha})$
(16)

when the weight function $\alpha=$ 2.0 is the sum of the square modulus
(abbreviated as SSM) of Kohn-Sham orbitals. It should be noticed that the SSM
weight function forms the second term of wavefunction and the PSM weight
function forms the fourth term of the wavefunction, because of the value of
weight function is less than 1, such improved weight function is an
alternative weight function that could show better results, superior to the
PSM weight function defined in Eq.15. The $\alpha=$ 1.0, 3.0, and 4.0
represent that the weight function is the sum of the modulus (SM), cubic
modulus (SCM) and quartic modulus (SQM) of Kohn-Sham orbitals, respectively.
Eq.16 is also valid for real-valued Kohn-Sham orbitals because real numbers
are subsets of complex numbers. FIG. 1 demonstrates the IPs selected by the
K-means with SSM and PSM as well as QRCP procedures for BH3NH3 with different
B-N distances. It is obvious that the K-means with SSM can yield IPs which are
more dispersed and delocalized than the K-means with PSM and QRCP, which
demonstrates that SSM is more suitable to simulate electron density as a
weight function. In addition, we verify the feasibility of SSM as the weight
function by demonstrating that the interpolation points using K-means with SSM
approximately minimize the residual for the ISDF decomposition (See APPENDIX
A).

Figure 1: Comparison of the IPs selected by the K-means with SSM and PSM as
well as QRCP procedures for BH3NH3 with different B-N distances of
$d_{\textrm{B-N}}$ = 1.6 Åand 2.8 Å, including (a, b) the electron density
(yellow isosurfaces), (c, d) the IPs (yellow pentagrams) by the K-means with
SSM, (e, f) the IPs (yellow pentagrams) by the K-means with PSM and (g, h) the
IPs (yellow pentagrams) by the QRCP procedure. Algorithm 1 K-means Clustering
Algorithm to Compute Interpolation Points in ISDF

Input: Grid points $\\{\mathbf{r}_{i}\\}_{i=1}^{N_{r}}$, Weight function
$w(\mathbf{r})$.

Output: Interpolation points $\\{\mathbf{r}_{\mu}\\}_{\mu=1}^{N_{\mu}}$

1: Initialize centroids $\\{\mathbf{r}_{\mu}^{\\{0\\}}\\}$, set $t\leftarrow
0$.

2: while convergence not reached do

3: Classification step: Assign $N_{r}$ points
$\\{\mathbf{r}_{i}\\}_{i=1}^{N_{r}}$ to the cluster $C_{\mu}^{(t)}$

4: Compute new centroids:
$\mathbf{r}_{\mu}^{(t+1)}\leftarrow{\sum_{\mathbf{r}_{j}\in
C_{\mu}^{(t)}}\mathbf{r}_{j}w(\mathbf{r}_{j})}/{\sum_{\mathbf{r}_{j}\in
C_{\mu}^{(t)}}w(\mathbf{r}_{j})}$

5: set $t\leftarrow t+1$

6: end while

7: Update
$\\{\mathbf{r}_{i}\\}_{i=1}^{N_{\mu}}\leftarrow\\{\mathbf{r}_{\mu}^{\\{t\\}}\\}$
.

###  2.3 Low rank approximation of complex-valued hybrid DFT via ISDF

The key spirit of DFT is to solve the Kohn-Sham equations expressed as

$H\psi_{j}=(-\frac{1}{2}\Delta_{\mathbf{r}}+V_{\text{ion}}+V_{\text{H}}[\rho]+V_{\text{XC}}[\\{\phi_{i}\\}])\psi_{j}=\epsilon_{j}\psi_{j}$
(17)

where $H$ is the Kohn-Sham Hamiltonian, $\psi_{j}$ is the $j$-th Kohn-Sham
orbital, $\epsilon_{j}$ is the corresponding orbital energy and
$V_{\text{ion}}$ is the ionic potential. In real space, the Hartree potential
$V_{\text{H}}$ is defined as

$V_{\textrm{H}}[\rho](\mathbf{r})=\int{\dfrac{\rho(\mathbf{r^{\prime}})}{|\mathbf{r}-\mathbf{r^{\prime}}|}}d\mathbf{r^{\prime}}$
(18)

the electron density is given as

$\rho(\mathbf{r})=\sum_{i=1}^{N_{e}}|\psi_{i}(\mathbf{r})|^{2}$ (19)

It should be noticed that the accuracy of the Kohn-Sham DFT calculations
strongly depends on the exchange-correlation potential
$V_{XC}[\\{\phi_{i}\\}]$, which is defined as

$V_{XC}[\\{\phi_{i}\\}]=V_{X}[\\{\phi_{i}\\}]+V_{C}[\rho]$ (20)

where $\\{\phi_{i}\\}$ denote the occupied orbitals. $V_{X}[\\{\phi_{i}\\}]$
and $V_{C}[\rho]$ represent the exchange and correlation potentials,
respectively. For complex-valued hybrid DFT, the Hartree-Fock exchange
operator is

$V_{X}[\\{\phi_{i}\\}](\mathbf{r},\mathbf{r^{\prime}})=-\sum_{i=1}^{N_{e}}{\frac{\phi_{i}^{\ast}(\mathbf{r})\phi_{i}(\mathbf{r^{\prime}})}{|\mathbf{r}-\mathbf{r^{\prime}}|}}$
(21)

When the Hartree-Fock exchange operator is applied to the orbitals

$(V_{X}[\\{\phi_{i}\\}]\psi_{j})(\mathbf{r})=-\sum_{i=1}^{N_{e}}\phi_{i}(\mathbf{r})\int{\frac{\phi_{i}^{\ast}(\mathbf{r^{\prime}})\psi_{j}(\mathbf{r^{\prime}})}{|\mathbf{r}-\mathbf{r^{\prime}}|}{d}\mathbf{r^{\prime}}}$
(22)

For large basis sets for discretizing the Kohn-Sham equations, such as the
plane-wave basis sets, it is more efficient to use an iterative
diagonalization procedure to compute the eq. 17. In practice, several DFT
packages, such as Quantum Espresso 71 and PWDFT 70, separate the self-
consistent field (SCF) iteration of all occupied orbitals into inner SCF
iteration and outer SCF iteration, called two-level SCF produce. In the inner
SCF iteration, the exchange operator $V_{X}$ defined by occupied orbitals
$\\{\phi_{i}\\}$ in eq. 21 is fixed, so that the Hamiltonian operator only
relies on the electron density $\rho$, which has to be updated constantly. In
the outer SCF iteration, the output orbitals from the inner SCF iteration will
be used for updating the exchange operator until it converges. In each inner
SCF iteration, we must solve the two-center integrals of eq. 22. The practical
numerical solution is to solve $O(N_{e}^{2})$ Poisson-like equations, which is
the most expensive part for hybrid DFT calculations within plane-wave basis
sets.

Under the ISDF decomposition of transposed Khatri-Rao product of complex-
valued Kohn-Sham orbitals
$Z=\\{z_{ij}:=\phi_{i}^{\ast}(\mathbf{r})\psi_{j}(\mathbf{r})\\}\in\mathbb{C}^{N_{r}\times
N_{e}^{2}}$, we can substitute Eq.(1) into Eq.(22)

$\begin{split}(V_{X}[\\{\phi_{i}\\}]\psi_{j})(\mathbf{r})&\approx-\sum_{i=1}^{N_{e}}\phi_{i}(\mathbf{r})\int{\frac{\sum_{\mu=1}^{N_{\mu}}\zeta_{\mu}(\mathbf{r^{\prime}})\phi_{i}^{\ast}(\mathbf{r_{\mu}})\psi_{j}(\mathbf{r_{\mu}})}{|\mathbf{r}-\mathbf{r^{\prime}}|}{d}\mathbf{r^{\prime}}}\\\
&=-\sum_{\mu=1}^{N_{\mu}}({\int{\frac{\zeta_{\mu}(\mathbf{r^{\prime}})}{|\mathbf{r}-\mathbf{r^{\prime}}|}}{d}\mathbf{r^{\prime}}{{\sum_{i=1}^{N_{e}}\phi_{i}(\mathbf{r})\phi_{i}^{\ast}(\mathbf{r_{\mu}})\psi_{j}(\mathbf{r_{\mu}})}}})\\\
&=-\sum_{\mu=1}^{N_{\mu}}({V_{\mu}(\mathbf{r}){{\sum_{i=1}^{N_{e}}\phi_{i}(\mathbf{r})\phi_{i}^{\ast}(\mathbf{r_{\mu}})\psi_{j}(\mathbf{r_{\mu}})}}})\end{split}$
(23)

where the projected Hartree-Fock exchange integral under the ISDF
decomposition is defined as

$\begin{split}V_{\mu}(\mathbf{r})=\int{\frac{\zeta_{\mu}(\mathbf{r^{\prime}})}{|\mathbf{r}-\mathbf{r^{\prime}}|}}{d}\mathbf{r^{\prime}}\end{split}$
(24)

As a consequence, the number of Poisson-like equations to be solved is reduced
to $N_{\mu}\sim O(N_{e})$ from $O(N_{e}^{2})$.

Figure 2: Three different types of data partition for the matrix used in the
ISDF formulation for hybrid density functional calculations: (a) 2D block
cyclic partition ($I_{R}\times I_{C}$ MPI processor grid), (b) 1D column
cyclic partition (1 × $P_{n}$ MPI processor grid), and (c) 1D row cyclic
partition ($P_{n}$ × 1 MPI processor grid). $P_{n}$ is total computational
cores and $I_{R}\times I_{C}=P_{n}$.

###  2.4 Parallel implementation

We implement this complex-valued K-means clustering algorithm for large-scale
hybrid DFT calculations within plane-wave basis sets in PWDFT, 70 which is an
open source plane-wave based electronic structure calculations software. We
also realize a parallel implementation of such low-scaling hybrid DFT
calculations in PWDFT as shown in FIG. 3.

Figure 3: Flowchart of the ISDF formulation in PWDFT. Red, blue and orange
boxes denote 1D column block cyclic partition, 2D block cyclic partition and
1D row block cyclic partition, respectively.

The discretized Kohn-Sham orbitals can be denoted by
$\Phi=[\phi_{1},\phi_{2},...,\phi_{N_{e}}]\in\mathbb{C}^{N_{r}\times N_{e}}$
and $\Psi=[\psi_{1},\psi_{2},...,\psi_{N}]\in\mathbb{C}^{N_{r}\times N}$. Thus
the parallel is implemented easily with the aid of ScaLAPACK library when
$P_{n}$ processors are used. There are three different data partition types
including 1D column cyclic partition, 1D row cyclic partition and 2D block
cyclic partition in our program, as shown in FIG. 2. It is easy for parallel
implementation when we apply the Hamiltonian operator to the orbitals and use
a sequential fast Fourier transformation (FFT) library by the 1D column block
cyclic partition. The 2D block partition is suitable for QRCP method and
matrix inversion, while the 1D row block partition should be adopted for
matrix multiplication and K-means method. The transform among the different
data partition types is achieved by the pdgemr2d subroutine in the ScaLAPACK
library. FIG. 3 demonstrates the flowchart of the ISDF formulation in PWDFT.
Firstly, the orbitals $\Phi$ and $\Psi$ can be stored using the 1D column
cyclic partition as the input. The interpolation points can be computed by the
K-means clustering algorithm. The weight functions are computed from $\Phi$
and $\Psi$. The $N_{r}$ grid points are equally distributed to each core to
compute the distances between grid points and centroids for parallel
implementation of the K-means clustering part. Then the quasi density matrices
$P^{\phi}(\mathbf{r}_{\mu},\mathbf{r}_{\nu})\in\mathbb{C}^{N_{\mu}\times
N_{\mu}}$ and
$P^{\psi}(\mathbf{r}_{\mu},\mathbf{r}_{\nu})\in\mathbb{C}^{N_{\mu}\times
N_{\mu}}$ in Eq. 9 are transformed into the 2D block cyclic partition to
construct the IVs. The matrix $ZC^{T}$ in Eq. 6 can be calculated in parallel
fashion and the matrix $CC^{T}$ is exactly subsampling rows of $ZC^{T}$. After
we obtain the IVs
$\Theta=[\zeta_{1},\zeta_{2},...,\zeta_{N_{\mu}}]\in\mathbb{C}^{N_{r}\times
N_{\mu}}$ by linear equation solver in ScaLAPACK, the data partition type
should be converted from 2D block cyclic partition to 1D column cyclic
partition for computing the
$V=[V_{1},V_{2},...,V_{N_{\mu}}]\in\mathbb{C}^{N_{r}\times N_{\mu}}$.

Table 2: Accuracy of complex-valued ISDF (SSM, $\alpha=2$) based HSE06 calculations by using K-means clustering algorithm to compute the IPs with respect to the rank parameter $t$ for liquid water molecules (H2O)64, semiconducting solid Si216 and metallic aluminium-silicon alloy Al176Si24, including the VBM $E_{\textrm{VBM}}$ (eV), CBM $E_{\textrm{CBM}}$ (eV), the energy gap $E_{\textrm{g}}$ (eV), the absolute errors of Hartree-Fock exchange energy $\Delta E_{\textrm{HF}}$ (Ha/atom), total energy $\Delta E$ (Ha/atom) and atomic forces $\Delta F$ (Ha/Bohr). The ACE-enabled HSE06 calculations are used for the reference. $t$ | $E_{\textrm{VBM}}$ | $E_{\textrm{CBM}}$ | $E_{\textrm{g}}$ | $\Delta E_{\textrm{HF}}$ | $\Delta E$ | $\Delta F$
---|---|---|---|---|---|---
Liquid water molecules (H2O)64 ($N_{\textrm{band}}=255$)
4.0 | -3.8136 | 2.4196 | 6.2333 | $2.22\times 10^{-04}$ | $3.12\times 10^{-04}$ | $1.30\times 10^{-02}$
6.0 | -3.8261 | 2.2910 | 6.1170 | $2.63\times 10^{-05}$ | $4.39\times 10^{-05}$ | $5.59\times 10^{-04}$
8.0 | -3.8291 | 2.2740 | 6.1031 | $3.46\times 10^{-06}$ | $1.49\times 10^{-05}$ | $3.22\times 10^{-04}$
10.0 | -3.8298 | 2.2713 | 6.1011 | $1.26\times 10^{-07}$ | $1.86\times 10^{-05}$ | $6.54\times 10^{-05}$
12.0 | -3.8301 | 2.2708 | 6.1010 | $5.68\times 10^{-07}$ | $2.81\times 10^{-05}$ | $5.01\times 10^{-05}$
16.0 | -3.8302 | 2.2706 | 6.1008 | $6.65\times 10^{-07}$ | $2.87\times 10^{-05}$ | $5.68\times 10^{-05}$
20.0 | -3.8300 | 2.2706 | 6.1006 | $3.65\times 10^{-09}$ | $4.17\times 10^{-07}$ | $2.11\times 10^{-05}$
Ref | -3.8299 | 2.2705 | 6.1004 | $0.00\times 10^{00}$ | $0.00\times 10^{00}$ | $0.00\times 10^{00}$
Semiconducting bulk silicon solid Si216 ($N_{\textrm{band}}=432$)
4.0 | 6.7419 | 8.3454 | 1.6035 | $2.53\times 10^{-03}$ | $2.96\times 10^{-03}$ | $2.73\times 10^{-03}$
6.0 | 6.6641 | 8.1544 | 1.4903 | $3.66\times 10^{-04}$ | $4.64\times 10^{-04}$ | $6.47\times 10^{-04}$
8.0 | 6.6513 | 8.1040 | 1.4527 | $7.45\times 10^{-05}$ | $9.97\times 10^{-05}$ | $1.82\times 10^{-04}$
10.0 | 6.6480 | 8.0959 | 1.4479 | $1.87\times 10^{-05}$ | $2.83\times 10^{-05}$ | $1.99\times 10^{-04}$
12.0 | 6.6472 | 8.0942 | 1.4470 | $4.73\times 10^{-06}$ | $9.66\times 10^{-06}$ | $4.56\times 10^{-05}$
16.0 | 6.6469 | 8.0932 | 1.4463 | $1.77\times 10^{-07}$ | $1.51\times 10^{-06}$ | $1.60\times 10^{-05}$
20.0 | 6.6468 | 8.0930 | 1.4462 | $4.27\times 10^{-07}$ | $3.66\times 10^{-07}$ | $5.35\times 10^{-06}$
Ref | 6.6468 | 8.0930 | 1.4462 | $0.00\times 10^{00}$ | $0.00\times 10^{00}$ | $0.00\times 10^{00}$
Metallic aluminium-silicon alloy Al176Si24 ($N_{\textrm{band}}=312$)
4.0 | 7.9370 | 8.0369 | 0.0999 | $3.93\times 10^{-03}$ | $4.20\times 10^{-03}$ | $4.09\times 10^{-03}$
6.0 | 7.8118 | 7.9168 | 0.1050 | $6.91\times 10^{-04}$ | $7.29\times 10^{-04}$ | $1.22\times 10^{-03}$
8.0 | 7.7760 | 7.8759 | 0.0999 | $8.50\times 10^{-05}$ | $8.82\times 10^{-05}$ | $3.20\times 10^{-04}$
10.0 | 7.7714 | 7.8703 | 0.0989 | $1.76\times 10^{-05}$ | $1.88\times 10^{-05}$ | $1.02\times 10^{-04}$
12.0 | 7.7708 | 7.8695 | 0.0987 | $6.13\times 10^{-06}$ | $7.04\times 10^{-06}$ | $4.82\times 10^{-05}$
16.0 | 7.7706 | 7.8693 | 0.0986 | $1.09\times 10^{-06}$ | $1.77\times 10^{-06}$ | $2.96\times 10^{-05}$
20.0 | 7.7705 | 7.8691 | 0.0986 | $1.73\times 10^{-07}$ | $6.25\times 10^{-07}$ | $1.18\times 10^{-05}$
Ref | 7.7705 | 7.8691 | 0.0986 | $0.00\times 10^{00}$ | $0.00\times 10^{00}$ | $0.00\times 10^{00}$

##  3 Results and discussion

In this section, we demonstrate the numerical accuracy and computational
efficiency of the K-means clustering algorithm of complex-valued Kohn-Sham
orbitals for ISDF to accelerate hybrid density functional calculations. All
calculations are implemented in the PWDFT 70 software package and Message
Passing Interface (MPI) is used for handling data communication. We use the
Hartwigsen-Goedecker-Hutter (HGH) norm-conserving pseudopotentials72 and the
HSE06 73 functional to describe the electronic structures of molecules and
solids.

Figure 4: Atomic structures of (a) insulator liquid water molecules (H2O)64,
(b) semiconducting bulk silicon solid Si216 and (c) metallic aluminium-silicon
alloy Al176Si24. The red, white, yellow and purple circles denote O, H, Si and
Al atoms, respectively.

We first benchmark the numerical accuracy of ISDF-K-means with weight function
SSM by comparing standard HSE06 results. Then we compare different numerical
accuracy of ISDF-K-means with SSM, ISDF-K-means with PSM and ISDF-QRCP.
Furthermore, we compare the AIMD results from K-means with different weight
functions. Finally, we demonstrate the computational scaling of ISDF-K-means
as well as ISDF-QRCP and the parallel scalability of ISDF-enabled hybrid
density functional calculations on modern heterogeneous supercomputers.

###  3.1 Numerical accuracy

####  3.1.1 Total energy and atomic forces

Firstly, we test the accuracy of the ISDF method for complex-valued Kohn-Sham
orbitals taking liquid water molecules (H2O)64, semiconducting bulk silicon
Si216 ($E_{gap}$ = 1.45 eV) and metallic disordered aluminium-silicon alloy
Al176Si24 ($E_{gap}$ $\textless$ 0.1 eV) as examples, whose crystal structures
are demonstrated in FIG. 4. The cutoff energies for (H2O)64, Si216 and
Al176Si24 are 60.0, 20.0 and 20.0 Ha, respectively. For (H2O)64 system, the
DFT-D2 method is used to account for the van der Waals (VdW) interaction. 74
We obtain all reference results utilizing adaptively compressed exchange (ACE)
algorithm, 75, 70 which reduces the cost of applying Hartree-Fock exchange
operator into Kohn-Sham orbitals without loss of accuracy.

Figure 5: Accuracy of complex-valued ISDF based hybrid functional calculations
(HSE06) obtained by using the K-means (SMS and PMS) and QRCP procedures to
select the interpolation points, with varying rank parameter t from 4 to 20
for Si216, Al176Si24 and (H2O)64, including the error of ((a), (b), (c))
Hartree-Fock exchange energy $\Delta E_{\textrm{HF}}$ (Ha/atom), ((d), (e),
(f)) total energy $\Delta E$ (Ha/atom) and ((g), (h), (i)) atomic forces
$\Delta F$ (Ha/Bohr).

Table 2 demonstrates the valence band maximum (VBM) energy level, the
conduction band minimum (CBM) energy level, the energy gap $E_{\textrm{g}}$,
the absolute errors of Hartree-Fock exchange energy $\Delta E_{\textrm{HF}}$,
total energy $\Delta E$ (Ha/atom) as well as atomic forces $\Delta F$
(Ha/Bohr) by ISDF-based HSE06 calculations using K-means with weight function
SMS with respect to the rank parameter $t$. Here we define

$\displaystyle\Delta
E_{\textrm{HF}}=(E_{\textrm{HF}}^{ISDF}-E_{\textrm{HF}}^{Ref})/N_{A}$ (25)
$\displaystyle\Delta E=(E^{ISDF}-E^{Ref})/N_{A}$ $\displaystyle\Delta
F=\max\limits_{I}|(F_{I}^{ISDF}-F_{I}^{Ref})|$

where $N_{A}$ is the total number of atoms, and $I$ is the atom index.

Figure 6: Variation of eigenvalue error in hybrid HSE06 DFT calculations using
K-means with SSM and PSM as well as QRCP method with respect to the rank
parameter $t$ ($t$ = 6 and $t$ = 12) for (a) Si216, (b) Al176Si24 and (c)
(H2O)64.

We remark that above physical quantities are approaching the reference value
gradually as the rank parameter $t$ is increasing for three studied systems.
Besides, for all tested systems, the error of energy gap is less than 0.01 eV
and the error of total energy per atom is less than $10^{-4}$ Ha/atom when $t$
is set to 8.0, which suggests that a small rank parameter can yield enough
accurate results. The accuracy of standard method QRCP as well as K-means
method with SSM and PSM are compared in FIG. 5. When the same rank parameter
$t$ is set, the K-means method with SSM yields higher accuracy than K-means
with PSM for Hartree-Fock exchange energy, total energy and atomic forces in
most cases. In order to further compare the accuracy of K-means with different
weight functions and QRCP, we demonstrate variation of eigenvalue error with
different rank parameters $t$ = 6.0 and 12.0 for Si216, Al176Si24 and (H2O)64,
as shown in FIG. 6. The error of the $i$-th eigenvalue is calculated by
$\Delta\epsilon_{i}=\epsilon_{i}^{ISDF}-\epsilon_{i}^{Ref}$. Similarly, it is
obvious that the K-means with SSM shows a lower error of eigenvalues for all
tested systems.

####  3.1.2 AIMD

In order to further verify the accuracy of the ISDF method with the improved
K-means algorithm, we perform AIMD optimized by hybrid DFT calculations for
bulk silicon system Si64, aluminium-silicon alloy Al176Si24 and silicon
dioxide SiO2 under the NVE ensemble and liquid water molecules (H2O)64 under
the NVT ensemble. For NVT ensemble, a single level Nose-Hoover thermostat 76,
77 is used at 295 K of temperature and the mass of the Nose-Hoover thermostat
is 85,000 a.u.

Figure 7: Energy potential of hybrid HSE06 DFT AIMD simulations by using the
ISDF-K-means with different weight functions and ISDF-QRCP methods as well as
ACE-enabled procedure as the reference on the bulk silicon Si64.

FIG. 7 demonstrates that ISDF-K-means with different weight functions and
ISDF-QRCP method results in the total potential energy along the MD trajectory
with a time step of 1.0 fs on bulk silicon system Si64. We remark that the
absolute error from K-means method with SSM is smaller than that from K-means
method with PSM. Moreover, the K-means method with SSM yields a much smoother
potential energy change compared with that yielded by the K-means method with
PSM. The possible reason is that the interpolation points selected by weight
function SSM represent the distribution of electron density better than that
by weight function PSM. Therefore, weight function SSM will perform better in
some cases where a smooth potential energy surface is required, such as
geometry optimization. In addition, we also compare the potential energy
change along the MD trajectory of K-means method with SM ($\alpha=1$), SSM
($\alpha=2$), SCM ($\alpha=3$) as well as SQM ($\alpha=4$). We remark that the
weight function SM exhibits a similar performance to the weight function SSM.
Note that the K-means method with SQM yields almost the same potential energy
change as that from the K-means method with PSM. It is because both weight
function SQM and PSM are essentially quartic modulus of wavefunctions in
hybrid functional calculations where $\Psi$ is the same as $\Phi$ when the
unoccupied orbitals are not computed.

Figure 8: Relatively energy drift of hybrid HSE06 AIMD simulations by using
the ISDF-K-means with different weight functions and ISDF-QRCP methods as well
as ACE-enabled procedure as the reference on the bulk silicon Si64.

Energy drift represents the change of total energy compared to the initial
total energy along the MD trajectory, which can be defined by
$E_{drift}(t)=(E_{tot}(t)-E_{tot}(0))/E_{tot}(0)$. FIG. 8 shows the controlled
energy drift by the ISDF-K-means with different weight functions and ISDF-QRCP
on bulk silicon system Si64. We remark that the K-means with SSM exhibits less
loss of accuracy than the K-means with PSM and QRCP methods. Besides,
divergence occurs for the energy drift of K-means with PSM while that of
K-means with SSM maintains stable accuracy. On the other hand, the performance
of the K-means method with SM is similar to that of the K-means method with
SSM, while the performance of the K-means method with SQM is similar to that
of the K-means method with PSM, which is consistent with the results from the
energy potential change. In addition, for AIMD with a time step of 1.0 fs on
aluminium-silicon alloy Al176Si24, as FIG.S1 and FIG.S2 show, the K-means
method with SSM also yields slightly smoother potential energy and less loss
of accuracy of energy drift than the K-means method with PSM. Similarly, for
liquid water molecules, the K-means with SSM exhibits higher accuracy than the
K-means with PSM on the whole, as shown in FIG.S5. From FIG.S4, the potential
energy yielded by the K-means with SSM is more consistent with the reference
results compared with that by the K-means with PSM for water molecules, which
suggests a stronger stability of SSM. In order to further check the stability
of the K-means with SSM method, we perform a long MD simulation for 21.5 ps
with a time step of 12.0 fs on Si64, as FIG. 9 shows. K-means with SSM
exhibits less loss of accuracy and better stability than K-means with PSM.

Figure 9: Relatively energy drift of hybrid HSE06 AIMD simulations with a time
step of 12.0 fs for Si64.

We also sample the oxygen-oxygen radial distribution functions (RDF) by hybrid
HSE06 AIMD simulations for 0.37 ps with a time step of 1.0 fs on liquid water
system (H2O)64 under NVT ensemble at 295 K using K-means with SSM and K-means
with PSM, as FIG. 10 shows. We remark that the root mean square error (RMSE)
of K-means with SSM (0.022) is much smaller than that of K-means with PSM
(0.042), which demonstrates better accuracy for SSM compared with PSM.
Therefore, this improved K-means clustering algorithm in the ISDF
decomposition can accurately and efficiently realize low-scaling, large-scale
and long-time AIMD with complex-valued hybrid DFT calculations.

Figure 10: Oxygen-oxygen radial distribution functions gOO(r) of the liquid
water system (H2O)64 at 295 K obtained from hybrid HSE06 + DFT-D2 AIMD
simulations for 0.37 ps with a time step 1.0 fs using the ISDF-K-means with
SSM (red solid line) and ISDF-K-means with PSM (blue solid line) methods and
the ACE algorithm (black dashed line as the reference). The root mean square
error (RMSE) of SSM and PSM with respect to the reference is 0.022 and 0.042,
respectively.

In addition to studying the oxygen-oxygen RDF of liquid water system (H2O)64,
we also calculate the power spectrum 78 of the crystal silicon dioxide SiO2 by
using HSE06 AIMD simulations under NVE ensemble at 300 K using K-means with
SSM and K-means with PSM. As shown in FIG. 11, the black vertical line shows
three different vibration frequencies in experiment 79, which is Si-O
asymmetrical bending vibration with 464 cm -1, Si-O symmetrical bending
vibration with 695 cm-1 and Si-O asymmetrical stretching vibration with
1080-1175 cm-1, respectively, we could see it is closer to the asymmetrical
bending vibration frequency and asymmetrical bending vibration frequency when
using K-means with SSM. Furthermore, we could find a small peak with 695 cm-1
using K-means with SSM, however, we couldn’t see this in K-means with PSM.
Therefore, the improved K-means clustering with SSM in ISDF algorithm could
simulate more accurate power spectral vibrational frequencies in hybrid AIMD
simulations.

Figure 11: Power spectrum of the crystal silicon dioxide system SiO2 at 300K
obtained from hybrid HSE06 AIMD simulations for 2ps with a time step 1.0 fs
using the ISDF-K-means with SSM (red solid line) and ISDF-K-means with PSM
(blue solid line) methods and the ACE algorithm (black dashed line as the
reference). The experimental results are shown as four black vertical lines.

###  3.2 Computational efficiency

####  3.2.1 Computational scaling

To compare the computational efficiency of K-means and QRCP method, we
demonstrate the computational time of IPs selected by QRCP and K-means method
and IVs for complex-valued wavefunctions in FIG.S5. The tested systems are
bulk silicon systems including from 8 to 512 atoms. The cutoff energy is set
to 10 Ha. From the fitting curves, we remark that selecting IPs by K-means is
much faster than that by QRCP. For periodic systems where $N_{r}$ is
proportional to the number of atoms, the K-means algorithm scales as
$O(N^{2.0})$ while the QRCP method scales as $O(N^{3.0})$, which is consistent
with the conclusion in numerical atomic orbitals. 34 The IVs procedure by the
least-square method scales as $O(N^{2.4})$. For ISDF-QRCP method, the IPs
procedure is the most time-consuming part due to the expensive QRCP. However,
the IVs procedure becomes the dominant part when the K-means method replaces
the QRCP.

####  3.2.2 Parallel scalability

Figure 12: (a) The change of wallclock time in ISDF, inner SCF and outer SCF
iteration with respect to the number of cores for the Si1000 system. (b) The
change of wallclock time in ISDF, inner SCF and outer SCF iteration with
respect to system size.

We test the parallel scalability of the complex-valued ISDF method with
K-means clustering algorithm for large-scale hybrid density functional (HSE06)
calculations, as shown in FIG. 12. FIG. 12(a) shows the change of the
wallclock time of ISDF part, inner SCF and outer SCF in one outer SCF
iteration with respect to the numbers of cores for the bulk silicon system
containing 1,000 atoms, which illustrates strong parallel scalability of our
algorithm. We remark that the wallclock time of inner SCF, ISDF and outer SCF
exhibits excellent scalability up to 2,000 CPU cores. As for the weak parallel
scalability, FIG. 12(b) demonstrates the change of wallclock time with respect
to the numbers of atoms for the bulk silicon system including from 64 to 2,744
atoms. The ISDF method scales well with respect to the system size (up to
2,744 atoms) on 5,504 CPU cores. The hybrid density functional calculations
for complex-valued Kohn-Sham orbitals require more usage of computation and
memory, which can be improved by performing OpenMP parallel implementation in
our future work. 80 Therefore, this improved K-means clustering algorithm can
accurately and efficiently accelerate large-scale and long-time ab initio
molecular dynamics with complex-valued hybrid DFT calculations.

##  4 Conclusion and outlook

In conclusion, we present an improved K-means clustering algorithm for
complex-valued wavefunctions to select interpolation sampling points in the
ISDF decomposition. By applying the new K-means clustering algorithm with SSM
into hybrid density functional calculations, we demonstrate that the improved
K-means clustering algorithm yields more accurate and smoother interpolation
sampling points compared to the K-means clustering with PSM for complex-valued
Kohn-Sham orbitals. In particular, K-means with SSM exhibits less loss of
accuracy and better stability for AIMD RDF and power spectrum simulations with
hybrid density functional. Moreover, we implement the parallel ISDF
decomposition for large-scale hybrid functional calculations. We show that the
ISDF can scale up to 5,504 CPU cores for a system containing 2,744 atoms. The
complex-valued wavefunction is indispensable for multi-k points sampling DFT
and RT-TDDFT. The K-means clustering algorithm is more suitable than QRCP
method for the dynamic simulation process because of its cheap cost.
Therefore, we will apply the complex-valued K-means clustering algorithm to
excited-state RT-TDDFT with hybrid functionals in our future work.

## Acknowledgments

This work is partly supported by the Strategic Priority Research Program of
the Chinese Academy of Sciences (XDB0450101), the Innovation Program for
Quantum Science and Technology (2021ZD0303306), the National Natural Science
Foundation of China (22288201, 22173093, 21688102), by the Anhui Provincial
Key Research and Development Program (2022a05020052), the National Key
Research and Development Program of China (2016YFA0200604, 2021YFB0300600),
and the CAS Project for Young Scientists in Basic Research (YSBR-005). The
authors thank the Hefei Advanced Computing Center, the Supercomputing Center
of Chinese Academy of Sciences, the Supercomputing Center of USTC, the
National Supercomputing Center in Wuxi, and Tianjin, Shanghai, and Guangzhou
Supercomputing Centers for the computational resources.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Verification of the feasibility for K-means with SSM

In order to verify the feasibility of SSM as the weight function, here we
demonstrate that the interpolation points using K-means with SSM approximately
minimize the residual for the ISDF decomposition. For simplicity, suppose
$N=N_{\phi}=N_{\psi}$, transposed Khatri-Rao product $Z$ is
$Z(\mathbf{r})=[\phi_{i}(\mathbf{r})\psi_{j}^{\ast}(\mathbf{r})]_{i,j=1}^{N}$.
We cluster $N_{r}$ matrix rows of Z into subsets
$\\{C_{\mu}\\}_{\mu=1}^{N_{\mu}}$ and select $N_{\mu}$ matrix row
$Z(\mathbf{r}_{\mu})$ for representing each $C_{\mu}$. Thus the error of ISDF
can be approximated as

$R=\sum_{\mu=1}^{N_{\mu}}\sum_{\mathbf{r_{k}}\in
C_{\mu}}||Z(\mathbf{r_{k}})-Proj_{span\\{Z(\mathbf{r}_{\mu})\\}}Z(\mathbf{r_{k}})||^{2}$
(26)

where the projection is the $L^{2}$ inner product

$Proj_{span\\{Z(\mathbf{r}_{\mu})\\}}Z(\mathbf{r_{k}})=\frac{Z(\mathbf{r_{k}})\cdot
Z^{\ast}(\mathbf{r}_{\mu})}{Z(\mathbf{r}_{\mu})\cdot
Z^{\ast}(\mathbf{r}_{\mu})}Z(\mathbf{r}_{\mu})$ (27)

we define electron density
$\rho(\mathbf{r}_{\mu})=\sum_{i=1}^{N}|\phi_{i}(\mathbf{r}_{\mu})|^{2}=\sum_{j=1}^{N}|\psi_{j}(\mathbf{r}_{\mu})|^{2}$,
$\Phi(\mathbf{r})=[\phi_{i}(\mathbf{r})]_{i=1}^{N}$,
$\Psi(\mathbf{r})=[\psi_{j}(\mathbf{r})]_{j=1}^{N}$.

$\begin{split}Z(\mathbf{r}_{\mu})\cdot
Z^{\ast}(\mathbf{r}_{\mu})&=(\Phi(\mathbf{r}_{\mu})\cdot\Psi^{\ast}(\mathbf{r}_{\mu}))(\Phi^{\ast}(\mathbf{r}_{\mu})\cdot\Psi(\mathbf{r}_{\mu}))\\\
&=\sum_{i,j=1}^{N}|\phi_{i}(\mathbf{r}_{\mu})|^{2}|\psi_{j}(\mathbf{r}_{\mu})|^{2}\\\
&=(\sum_{i=1}^{N}|\phi_{i}(\mathbf{r}_{\mu})|^{2})(\sum_{j=1}^{N}|\psi_{j}(\mathbf{r}_{\mu})|^{2})\\\
&=\rho^{2}(\mathbf{r}_{\mu})\end{split}$ (28)

we have

$\begin{split}R&=\sum_{\mu=1}^{N_{\mu}}\sum_{\mathbf{r}\in
C_{\mu}}||Z(\mathbf{r_{k}})-\frac{Z(\mathbf{r_{k}})\cdot
Z^{\ast}(\mathbf{r}_{\mu})}{Z(\mathbf{r}_{\mu})\cdot
Z^{\ast}(\mathbf{r}_{\mu})}Z(\mathbf{r}_{\mu})||^{2}\\\
&=\sum_{\mu=1}^{N_{\mu}}\sum_{\mathbf{r}\in C_{\mu}}Z(\mathbf{r_{k}})\cdot
Z^{\ast}(\mathbf{r_{k}})[1-\frac{(Z^{\ast}(\mathbf{r_{k}})\cdot
Z(\mathbf{r}_{\mu}))(Z(\mathbf{r_{k}})\cdot
Z^{\ast}(\mathbf{r}_{\mu}))}{(Z(\mathbf{r_{k}})\cdot
Z^{\ast}(\mathbf{r_{k}}))(Z(\mathbf{r}_{\mu})\cdot
Z^{\ast}(\mathbf{r}_{\mu}))}]\\\ &=\sum_{\mu=1}^{N_{\mu}}\sum_{\mathbf{r}\in
C_{\mu}}\rho^{2}(\mathbf{r}_{k})[1-\frac{(\Phi(\mathbf{r}_{k})\cdot\Phi^{\ast}(\mathbf{r}_{\mu}))^{2}(\Psi(\mathbf{r}_{k})\cdot\Psi^{\ast}(\mathbf{r}_{\mu}))^{2}}{\rho^{2}(\mathbf{r}_{k})\rho^{2}(\mathbf{r}_{\mu})}]\\\
&=\sum_{\mu=1}^{N_{\mu}}\sum_{\mathbf{r}\in
C_{\mu}}\rho^{2}(\mathbf{r}_{k})[1-cos^{2}(\theta_{1}(\mathbf{r}_{k},\mathbf{r}_{\mu}))cos^{2}(\theta_{2}(\mathbf{r}_{k},\mathbf{r}_{\mu}))]\\\
&=\sum_{\mu=1}^{N_{\mu}}\sum_{\mathbf{r}\in
C_{\mu}}\rho^{2}(\mathbf{r}_{k})[sin^{2}(\theta_{1}(\mathbf{r}_{k},\mathbf{r}_{\mu}))+sin^{2}(\theta_{2}(\mathbf{r}_{k},\mathbf{r}_{\mu}))\quad...\\\
&\quad\quad\quad\quad-
sin^{2}(\theta_{1}(\mathbf{r}_{k},\mathbf{r}_{\mu}))sin^{2}(\theta_{2}(\mathbf{r}_{k},\mathbf{r}_{\mu}))]\\\
&\leq\sum_{\mu=1}^{N_{\mu}}\sum_{\mathbf{r}\in
C_{\mu}}\rho^{2}(\mathbf{r}_{k})[sin^{2}(\theta_{1}(\mathbf{r}_{k},\mathbf{r}_{\mu}))+sin^{2}(\theta_{2}(\mathbf{r}_{k},\mathbf{r}_{\mu}))]\end{split}$
(29)

where $\theta_{1}(\mathbf{r}_{k},\mathbf{r}_{\mu})$ and
$\theta_{2}(\mathbf{r}_{k},\mathbf{r}_{\mu})$ are the angles between the
vectors $\Phi(\mathbf{r}_{k})$ and $\Phi^{\ast}(\mathbf{r}_{\mu})$ as well as
$\Psi(\mathbf{r}_{k})$ and $\Psi^{\ast}(\mathbf{r}_{\mu})$, respectively.
Because

$\begin{split}&\rho(\mathbf{r}_{k})[sin^{2}(\theta_{1}(\mathbf{r}_{k},\mathbf{r}_{\mu}))+sin^{2}(\theta_{2}(\mathbf{r}_{k},\mathbf{r}_{\mu}))]\\\
&=\Phi(\mathbf{r}_{k})\cdot\Phi^{\ast}(\mathbf{r}_{k})sin^{2}(\theta_{1}(\mathbf{r}_{k},\mathbf{r}_{\mu}))+\Psi(\mathbf{r}_{k})\cdot\Psi^{\ast}(\mathbf{r}_{k})sin^{2}(\theta_{2}(\mathbf{r}_{k},\mathbf{r}_{\mu}))\\\
&\leq||\Phi(\mathbf{r}_{k})-\Phi^{\ast}(\mathbf{r}_{\mu})||^{2}+||\Psi(\mathbf{r}_{k})-\Psi^{\ast}(\mathbf{r}_{\mu})||^{2}\end{split}$
(30)

we can obtain

$\begin{split}R&\leq\sum_{\mu=1}^{N_{\mu}}\sum_{\mathbf{r}\in
C_{\mu}}\rho(\mathbf{r}_{k})[||\Phi(\mathbf{r}_{k})-\Phi^{\ast}(\mathbf{r}_{\mu})||^{2}+||\Psi(\mathbf{r}_{k})-\Psi^{\ast}(\mathbf{r}_{\mu})||^{2}]\\\
&\approx\sum_{\mu=1}^{N_{\mu}}\sum_{\mathbf{r}\in
C_{\mu}}\rho(\mathbf{r}_{k})(||\nabla_{r}\Phi(\mathbf{r}_{\mu})||^{2}+||\nabla_{r}\Psi(\mathbf{r}_{\mu})||^{2})||\mathbf{r}_{k}-\mathbf{r}_{\mu}||^{2}\\\
&=\frac{1}{2}\sum_{\mu=1}^{N_{\mu}}\sum_{\mathbf{r}\in
C_{\mu}}(\sum_{i=1}^{N}|\phi_{i}(\mathbf{r}_{k})|^{2}+\sum_{j=1}^{N}|\psi_{j}(\mathbf{r}_{k})|^{2})\quad...\\\
&\quad\quad\quad\quad(||\nabla_{r}\Phi(\mathbf{r}_{\mu})||^{2}+||\nabla_{r}\Psi(\mathbf{r}_{\mu})||^{2})||\mathbf{r}_{k}-\mathbf{r}_{\mu}||^{2}\end{split}$
(31)

Thus the minimization criterion of weighted K-means with SSM can be derived
when the spatial inhomogeneity of the gradient $\Phi(\mathbf{r})$ and
$\Psi(\mathbf{r})$ is neglected.

## References

  * Slyusar 1998 Slyusar, V. End Products in Matrices in Radar Applications. _Radioelectron. Commun. Syst._ 1998, _41_ , 50–53
  * Khatri and Rao 1968 Khatri, C. G.; Rao, C. R. Solutions to Some Functional Equations and Their Applications to Characterization of Probability Distributions. _Sankhya_ 1968, _30_ , 167–180
  * Hohenberg and Kohn 1964 Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. _Phys. Rev._ 1964, _136_ , B864
  * Kohn and Sham 1965 Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. _Phys. Rev._ 1965, _140_ , A1133
  * Slater 1951 Slater, J. C. A Simplification of the Hartree-Fock Method. _Phys. Rev._ 1951, _81_ , 385
  * Becke 1993 Becke, A. D. A New Mixing of Hartree–Fock and Local Density-Functional Theories. _J. Chem. Phys._ 1993, _98_ , 1372–1377
  * Casida 1995 Casida, M. E. Response Theory for Molecules. _Recent Advances in Density Functional Methods:(Part I)_ 1995, _1_ , 155
  * Sternheimer 1954 Sternheimer, R. M. Electronic Polarizabilities of Ions from the Hartree-Fock Wave Functions. _Phys. Rev._ 1954, _96_ , 951
  * Hedin 1965 Hedin, L. New Method for Calculating the One-Particle Green’s Function with Application to the Electron-gas Problem. _Phys. Rev._ 1965, _139_ , A796
  * Hybertsen and Louie 1986 Hybertsen, M. S.; Louie, S. G. Electron Correlation in Semiconductors and Insulators: Band Gaps and Quasiparticle Energies. _Phys. Rev. B_ 1986, _34_ , 5390
  * Aryasetiawan and Gunnarsson 1998 Aryasetiawan, F.; Gunnarsson, O. The GW method. _Rep. Prog. Phys._ 1998, _61_ , 237
  * Onida et al. 2002 Onida, G.; Reining, L.; Rubio, A. Electronic Excitations: Density-Functional versus Many-body Green’s-Function Approaches. _Rev. Mod. Phys._ 2002, _74_ , 601
  * Salpeter and Bethe 1951 Salpeter, E. E.; Bethe, H. A. A Relativistic Equation for Bound-State Problems. _Phys. Rev._ 1951, _84_ , 1232
  * Häser 1993 Häser, M. Møller-Plesset (MP2) Perturbation Theory for Large Molecules. _Theor. Chim. Acta._ 1993, _87_ , 147–173
  * Feyereisen et al. 1993 Feyereisen, M.; Fitzgerald, G.; Komornicki, A. Use of Approximate Integrals in Ab Initio Theory. An Application in MP2 Energy Calculations. _Chem. Phys. Lett._ 1993, _208_ , 359–363
  * Bernholdt and Harrison 1996 Bernholdt, D. E.; Harrison, R. J. Large-scale Correlated Electronic Structure Calculations: the RI-MP2 Method on Parallel Computers. _Chem. Phys. Lett._ 1996, _250_ , 477–484
  * Ren et al. 2012 Ren, X.; Rinke, P.; Joas, C.; Scheffler, M. Random-Phase Approximation and Its Applications in Computational Chemistry and Materials Science. _J. Mater. Sci._ 2012, _47_ , 7447–7471
  * Beebe and Linderberg 1977 Beebe, N. H. F.; Linderberg, J. Simplifications in the Generation and Transformation of Two-electron Integrals in Molecular Calculations. _Int. J. Quantum. Chem._ 1977, _12_ , 683–705
  * Røeggen and Wisløff-Nilssen 1986 Røeggen, I.; Wisløff-Nilssen, E. On the Beebe-Linderberg Two-Electron Integral Approximation. _Chem. Phys. Lett._ 1986, _132_ , 154–160
  * Whitten 1973 Whitten, J. L. Coulombic Potential Energy Integrals and Approximations. _J. Chem. Phys._ 1973, _58_ , 4496–4501
  * Dunlap et al. 1979 Dunlap, B. I.; Connolly, J. W. D.; Sabin, J. R. On Some Approximations in Applications of X$\alpha$ Theory. _J. Chem. Phys._ 1979, _71_ , 3396–3402
  * Vahtras et al. 1993 Vahtras, O.; Almlöf, J.; Feyereisen, M. Integral Approximations for LCAO-SCF Calculations. _Chem. Phys. Lett._ 1993, _213_ , 514–518
  * Weigend 2002 Weigend, F. A Fully Direct RI-HF Algorithm: Implementation, Optimised Auxiliary Basis Sets, Demonstration of Accuracy and Efficiency. _Phys. Chem. Chem. Phys._ 2002, _4_ , 4285–4291
  * Hohenstein et al. 2012 Hohenstein, E. G.; Parrish, R. M.; Martínez, T. J. Tensor Hypercontraction Density Fitting. I. Quartic Scaling Second- and Third-Order Møller-Plesset Perturbation Theory. _J. Chem. Phys._ 2012, _137_ , 044103
  * Parrish et al. 2012 Parrish, R. M.; Hohenstein, E. G.; Martínez, T. J.; Sherrill, C. D. Tensor Hypercontraction. II. Least-squares Renormalization. _J. Chem. Phys._ 2012, _137_ , 224106
  * Hohenstein et al. 2012 Hohenstein, E. G.; Parrish, R. M.; Sherrill, C. D.; Martínez, T. J. Communication: Tensor Hypercontraction. III. Least-squares Tensor Hypercontraction for the Determination of Correlated Wavefunctions. _J. Chem. Phys._ 2012, _137_ , 221101
  * Friesner 1985 Friesner, R. A. Solution of Self-Consistent Field Electronic Structure Equations by a Pseudospectral Method. _Chem. Phys. Lett._ 1985, _116_ , 39–43
  * Friesner 1987 Friesner, R. A. Solution of the Hartree–Fock Equations for Polyatomic Molecules by a Pseudospectral Method. _J. Chem. Phys._ 1987, _86_ , 3522–3531
  * Lu and Ying 2015 Lu, J.; Ying, L. Compression of the Electron Repulsion Integral Tensor in Tensor Hypercontraction Format with Cubic Scaling Cost. _J. Comput. Phys._ 2015, _302_ , 329–335
  * Lu and Thicke 2017 Lu, J.; Thicke, K. Cubic Scaling Algorithms for RPA Correlation Using Interpolative Separable Density Fitting. _J. Comput. Phys._ 2017, _351_ , 187–202
  * Hu et al. 2017 Hu, W.; Lin, L.; Yang, C. Interpolative Separable Density Fitting Decomposition for Accelerating Hybrid Density Functional Calculations with Applications to Defects in Silicon. _J. Chem. Theory Comput._ 2017, _13_ , 5420–5431
  * Dong et al. 2018 Dong, K.; Hu, W.; Lin, L. Interpolative Separable Density Fitting Through Centroidal Voronoi Tessellation with Applications to Hybrid Functional Electronic Structure Calculations. _J. Chem. Theory Comput._ 2018, _14_ , 1311–1320
  * Qin et al. 2020 Qin, X.; Liu, J.; Hu, W.; Yang, J. Interpolative Separable Density Fitting Decomposition for Accelerating Hartree–Fock Exchange Calculations within Numerical Atomic Orbitals. _J. Phys. Chem. A_ 2020, _124_ , 5664–5674
  * Qin et al. 2020 Qin, X.; Li, J.; Hu, W.; Yang, J. Machine Learning K-Means Clustering Algorithm for Interpolative Separable Density Fitting to Accelerate Hybrid Functional Calculations with Numerical Atomic Orbitals. _J. Phys. Chem. A_ 2020, _124_ , 10066–10074
  * Malone et al. 2018 Malone, F. D.; Zhang, S.; Morales, M. A. Overcoming the Memory Bottleneck in Auxiliary Field Quantum Monte Carlo Simulations with Interpolative Separable Density Fitting. _J. Chem. Theory Comput._ 2018, _15_ , 256–264
  * Hu et al. 2020 Hu, W.; Liu, J.; Li, Y.; Ding, Z.; Yang, C.; Yang, J. Accelerating Excitation Energy Computation in Molecules and Solids within Linear-Response Time-Dependent Density Functional Theory via Interpolative Separable Density Fitting Decomposition. _J. Chem. Theory Comput._ 2020, _16_ , 964–973
  * Lee et al. 2019 Lee, J.; Lin, L.; Head-Gordon, M. Systematically Improvable Tensor Hypercontraction: Interpolative Separable Density-Fitting for Molecules Applied to Exact Exchange, Second-and Third-order Møller–Plesset Perturbation Theory. _J. Chem. Theory Comput._ 2019, _16_ , 243–263
  * Gao and Chelikowsky 2020 Gao, W.; Chelikowsky, J. R. Accelerating Time-Dependent Density Functional Theory and GW Calculations for Molecules and Nanoclusters with Symmetry Adapted Interpolative Separable Density Fitting. _J. Chem. Theory Comput._ 2020, _16_ , 2216–2223
  * Ma et al. 2021 Ma, H.; Wang, L.; Wan, L.; Li, J.; Qin, X.; Liu, J.; Hu, W.; Lin, L.; Yang, C.; Yang, J. Realizing Effective Cubic-Scaling Coulomb Hole Plus Screened Exchange Approximation in Periodic Systems via Interpolative Separable Density Fitting with a Plane-Wave Basis Set. _J. Phys. Chem. A_ 2021, _125_ , 7545–7557
  * Hu et al. 2018 Hu, W.; Shao, M.; Cepellotti, A.; Felipe, H.; Lin, L.; Thicke, K.; Yang, C.; Louie, S. G. Accelerating Optical Absorption Spectra and Exciton Energy Computation via Interpolative Separable Density Fitting. International Conference on Computational Science. 2018; pp 604–617
  * Vassilvitskii and Arthur 2006 Vassilvitskii, S.; Arthur, D. k-means++: The Advantages of Careful Seeding. Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms. 2006; pp 1027–1035
  * Singh et al. 2013 Singh, A.; Yadav, A.; Rana, A. K-means with Three Different Distance Metrics. _Int. J. Comput. Appl._ 2013, _67_ , 13–17
  * Bloch 1929 Bloch, F. Über Die Quantenmechanik Der Elektronen in Kristallgittern. _Z. Phys._ 1929, _52_ , 555–600
  * Bloch 1957 Bloch, F. Generalized Theory of Relaxation. _Phys. Rev._ 1957, _105_ , 1206
  * Wu et al. 2021 Wu, K.; Qin, X.; Hu, W.; Yang, J. Low-Rank Approximations Accelerated Plane-Wave Hybrid Functional Calculations with K-point Sampling. _J. Chem. Theory Comput._ 2021, _18_ , 206–218
  * Chen et al. 2023 Chen, S.; Wu, K.; Hu, W.; Yang, J. Low-Rank Approximations for Accelerating Plane-Wave Hybrid Functional Calculations in Unrestricted and Noncollinear Spin Density Functional Theory. _J. Chem. Phys._ 2023, _158_ , 134106
  * Runge and Gross 1984 Runge, E.; Gross, E. K. Density-Functional Theory for Time-Dependent Systems. _Phys. Rev. Lett._ 1984, _52_ , 997
  * Yabana and Bertsch 1996 Yabana, K.; Bertsch, G. Time-Dependent Local-Density Approximation in Real Time. _Phys. Rev. B_ 1996, _54_ , 4484
  * Andreani et al. 2005 Andreani, C.; Colognesi, D.; Mayers, J.; Reiter, G.; Senesi, R. Measurement of Momentum Distribution of Lightatoms and Molecules in Condensed Matter Systems using Inelastic Neutron Scattering. _Adv. Phys._ 2005, _54_ , 377–469
  * Li et al. 2023 Li, J.; Wan, L.; Jiao, S.; Hu, W.; Yang, J. Low-rank Approximations to Accelerate Hybrid Functional Enabled Real-Time Time-Dependent Density Functional Theory within Plane Waves. _Electron. Struct._ 2023, _5_ , 014008
  * Wu et al. 2021 Wu, K.-D.; Kondra, T. V.; Rana, S.; Scandolo, C. M.; Xiang, G.-Y.; Li, C.-F.; Guo, G.-C.; Streltsov, A. Operational Resource Theory of Imaginarity. _Phys. Rev. Lett._ 2021, _126_ , 090401
  * Zhang and Wu 2018 Zhang, J.; Wu, Y. Complex-Valued Unsupervised Convolutional Neural Networks for Sleep Stage Classification. _Comput. Methods Programs Biomed._ 2018, _164_ , 181–191
  * Santra et al. 2007 Santra, B.; Michaelides, A.; Scheffler, M. On the Accuracy of Density-Functional Theory Exchange-Correlation Functionals for H Bonds in Small Water Clusters: Benchmarks Approaching the Complete Basis Set Limit. _J. Chem. Phys._ 2007, _127_ , 184104
  * Santra et al. 2009 Santra, B.; Michaelides, A.; Scheffler, M. Coupled Cluster Benchmarks of Water Monomers and Dimers Extracted from Density-Functional Theory Liquid Water: The Importance of Monomer Deformations. _J. Chem. Phys._ 2009, _131_ , 124509
  * Santra et al. 2011 Santra, B.; Klimeš, J.; Alfe, D.; Tkatchenko, A.; Slater, B.; Michaelides, A.; Car, R.; Scheffler, M. Hydrogen Bonds and Van der Waals Forces in Ice at Ambient and High Pressures. _Phys. Rev. Lett._ 2011, _107_ , 185701
  * Santra et al. 2013 Santra, B.; Klimeš, J.; Tkatchenko, A.; Alfè, D.; Slater, B.; Michaelides, A.; Car, R.; Scheffler, M. On the Accuracy of Van der Waals Inclusive Density-Functional Theory Exchange-Correlation Functionals for Ice at Ambient and High Pressures. _J. Chem. Phys._ 2013, _139_ , 154702
  * DiStasio Jr et al. 2014 DiStasio Jr, R. A.; Santra, B.; Li, Z.; Wu, X.; Car, R. The Individual and Collective Effects of Exact Exchange and Dispersion Interactions on the Ab Initio Structure of Liquid Water. _J. Chem. Phys._ 2014, _141_ , 084502
  * Ko et al. 2020 Ko, H.-Y.; Jia, J.; Santra, B.; Wu, X.; Car, R.; DiStasio Jr, R. A. Enabling Large-Scale Condensed-Phase Hybrid Density Functional Theory Based Ab Initio Molecular Dynamics. 1. Theory, Algorithm, and Performance. _J. Chem. Theory Comput._ 2020, _16_ , 3757–3785
  * Khoo et al. 2011 Khoo, K.; Chan, T.-L.; Kim, M.; Chelikowsky, J. R. Ab Initio Molecular Dynamics Simulations of Molten Al1-xSix Alloys. _Phys. Rev. B_ 2011, _84_ , 214203
  * Zhang et al. 2017 Zhang, G.; Lin, L.; Hu, W.; Yang, C.; Pask, J. E. Adaptive Local Basis Set for Kohn–Sham Density Functional Theory in a Discontinuous Galerkin Framework II: Force, Vibration, and Molecular Dynamics Calculations. _J. Comput. Phys._ 2017, _335_ , 426–443
  * Chawla and Voth 1998 Chawla, S.; Voth, G. A. Exact Exchange in Ab Initio Molecular Dynamics: An Efficient Plane-Wave Based Algorithm. _J. Chem. Phys._ 1998, _108_ , 4697–4700
  * Mandal et al. 2021 Mandal, S.; Thakkur, V.; Nair, N. N. Achieving an Order of Magnitude Speedup in Hybrid-Functional-and Plane-Wave-Based Ab Initio Molecular Dynamics: Applications to Proton-Transfer Reactions in Enzymes and in Solution. _J. Chem. Theory Comput._ 2021, _17_ , 2244–2255
  * Guidon et al. 2008 Guidon, M.; Schiffmann, F.; Hutter, J.; VandeVondele, J. Ab Initio Molecular Dynamics Using Hybrid Density Functionals. _J. chem. Phys._ 2008, _128_ , 214204
  * Ko et al. 2020 Ko, H.-Y.; Jia, J.; Santra, B.; Wu, X.; Car, R.; DiStasio Jr, R. A. Enabling Large-scale condensed-Phase Hybrid Density Functional Theory Based Ab Initio Molecular Dynamics. 1. Theory, Algorithm, and Performance. _J. Chem. Theory Comput._ 2020, _16_ , 3757–3785
  * Mandal et al. 2021 Mandal, S.; Thakkur, V.; Nair, N. N. Achieving An Order of Magnitude Speedup in Hybrid-Functional-and Plane-Wave-Based Ab Initio Molecular Dynamics: Applications to Proton-Transfer Reactions in Enzymes and in Solution. _J. Chem. Theory Comput._ 2021, _17_ , 2244–2255
  * Ko et al. 2021 Ko, H.-Y.; Santra, B.; DiStasio Jr, R. A. Enabling Large-Scale Condensed-Phase Hybrid Density Functional Theory-Based Ab Initio Molecular Dynamics II: Extensions to the Isobaric–Isoenthalpic and Isobaric–Isothermal Ensembles. _J. Chem. Theory Comput._ 2021, _17_ , 7789–7813
  * Hohenstein et al. 2012 Hohenstein, E. G.; Parrish, R. M.; Martínez, T. J. Tensor Hypercontraction Density Fitting. I. Quartic Scaling Second-and Third-Order Møller-Plesset Perturbation Theory. _J. Chem. Phys._ 2012, _137_ , 044103
  * Parrish et al. 2012 Parrish, R. M.; Hohenstein, E. G.; Martínez, T. J.; Sherrill, C. D. Tensor Hypercontraction. II. Least-squares Renormalization. _J. Chem. Phys._ 2012, _137_ , 224106
  * Parrish et al. 2013 Parrish, R. M.; Hohenstein, E. G.; Martínez, T. J.; Sherrill, C. D. Discrete Variable Representation in Electronic Structure Theory: Quadrature Grids for Least-Squares Tensor Hypercontraction. _J. Chem. Phys._ 2013, _138_ , 194107
  * Hu et al. 2017 Hu, W.; Lin, L.; Banerjee, A. S.; Vecharynski, E.; Yang, C. Adaptively Compressed Exchange Operator for Large-Scale Hybrid Density Functional Calculations with Applications to the Adsorption of Water on Silicene. _J. Chem. Theory Comput._ 2017, _13_ , 1188–1198
  * Giannozzi et al. 2009 Giannozzi, P. et al. QUANTUM ESPRESSO: A Modular and Open-Source Software Project for Quantum Simulations of Materials. _J. Phys.: Condens. Matter_ 2009, _21_ , 395502
  * Hartwigsen et al. 1998 Hartwigsen, C.; Goedecker, S.; Hutter, J. Relativistic Separable Dual-space Gaussian Pseudopotentials from H to Rn. _Phys. Rev. B_ 1998, _58_ , 3641–3662
  * Heyd et al. 2006 Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Erratum: ”Hybrid Functionals Based on a Screened Coulomb Potential” [J. Chem. Phys.118, 8207 (2003)]. _J. Chem. Phys._ 2006, _118_ , 8207
  * Grimme 2006 Grimme, S. Semiempirical GGA-Type density functional constructed with a long-range dispersion correction. _J. Comput. Chem._ 2006, _27_ , 1787–1799
  * Lin 2016 Lin, L. Adaptively Compressed Exchange Operator. _J. Chem. Theory Comput._ 2016, _12_ , 2242–2249
  * Nosé 1984 Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. _J. Chem. Phys._ 1984, _81_ , 511–519
  * Hoover 1985 Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. _Phys. Rev. A: At., Mol., Opt. Phys._ 1985, _31_ , 1695
  * Thomas et al. 2013 Thomas, M.; Brehm, M.; Fligg, R.; Vöhringer, P.; Kirchner, B. Computing Vibrational Spectra from Ab Initio Molecular Dynamics. _Phys. Chem. Chem. Phys._ 2013, _15_ , 6608–6622
  * Saikia et al. 2008 Saikia, B. J.; Parthasarathy, G.; Sarmah, N. Fourier Transform Infrared Spectroscopic Estimation of Crystallinity in SiO 2 Based Rocks. _Bull. Mater. Sci._ 2008, _31_ , 775–779
  * Wan et al. 2021 Wan, L.; Liu, X.; Liu, J.; Qin, X.; Hu, W.; Yang, J. Hybrid MPI and OpenMP Parallel Implementation of Large-scale Linear-Response Time-Dependent Density Functional Theory with Plane-Wave Basis Set. _Electron. Struct._ 2021, _3_ , 024004

TOC graphic