*Corresponding author: E-mail<EMAIL_ADDRESS>
$\dagger$ Both authors contributed equally to this work.

# The static parallel distribution algorithms for hybrid density-functional
calculations in HONPAS package

Xinming Qin2† Honghui Shang*1† Lei Xu1 Wei Hu2 Jinlong Yang2 Shigang Li1
Yunquan Zhang1 11affiliationmark:  State Key Laboratory of Computer
Architecture, Institute of Computing Technology, Chinese Academy of Sciences,
Beijing, China
22affiliationmark:  Hefei National Laboratory for Physical Sciences at
Microscale, Department of Chemical Physics, and Synergetic Innovation Center
of Quantum Information and Quantum Physics, University of Science and
Technology of China, Hefei, Anhui 230026, China

###### Abstract

Hybrid density-functional calculation is one of the most commonly adopted
electronic structure theory used in computational chemistry and materials
science because of its balance between accuracy and computational cost.
Recently, we have developed a novel scheme called NAO2GTO to achieve linear
scaling (Order-N) calculations for hybrid density-functionalsShang et al.
(2011). In our scheme, the most time-consuming step is the calculation of the
electron repulsion integrals (ERIs) part. So how to create an even
distribution of these ERIs in parallel implementation is an issue of
particular importance. Here, we present two static scalable distributed
algorithms for the ERIs computation. Firstly, the ERIs are distributed over
ERIs shell pairs. Secondly, the ERIs is distributed over ERIs shell quartets.
In both algorithms, the calculation of ERIs is independent of each other, so
the communication time is minimized. We show our speedup results to
demonstrate the performance of these static parallel distributed algorithms in
the Hefei Order-N packages for ab initio simulations (HONPAS)Qin et al.
(2014).

##  1 INTRODUCTION

The electronic structure calculations based on density functional theory (DFT)
Parr and Yang (1989); Hohenberg and Kohn (1964); Kohn and Sham (1965) are the
work-horse of computational chemistry and materials science. However, widely
used semi-local density functionals could underestimate the band gaps because
of its inclusion of the unphysical self-interaction Mori-Sánchez et al.
(2008). A possible solution is to add the nonlocal Hartree-Fock type exchange
(HFX) into semi-local density-functionals to construct hybrid functionalsBecke
(1993); Stephens et al. (1994); Janesko et al. (2009); Paier et al. (2009);
Monkhorst (1979); Delhalle and Calais (1987); Gell-Mann and Brueckner (1957);
Heyd et al. (2003, 2006); Krukau et al. (2006); Frisch et al. (2009); Paier et
al. (2006). However, the drawback of hybrid density-functionals is that it is
significantly more expensive than conventional DFT. The most time-consuming
part in hybrid density-functional calculations becomes construction of HFX
matrix, even with the appearance of fast linear scaling algorithms that
overcome the bottlenecks encountered in conventional methodsSchwegler and
Challacombe (1996); Burant et al. (1996); Schwegler et al. (1997); Ochsenfeld
et al. (1998); Polly et al. (2004); Tymczak and Challacombe (2005); Sodt and
Head-Gordon (2008); Guidon et al. (2010); Merlot et al. (2014). As a result,
hybrid density-functional calculations must make efficient use of parallel
computing resources in order to reduce the execution time of HFX matrix
construction.

The implementation of hybrid density-functionals for solid state physics
calculations are mostly based on plane waves (PW)Gonze et al. (2002, 2016);
Paier et al. (2006) or linear combination of atomic orbitals (LCAO)Krukau et
al. (2006); Frisch et al. (2009); Dovesi et al. (2006) method. The atomic
orbitals basis set is efficient for real-space formalisms, which have
attracted considerable interest for DFT calculations because of their
favorable scaling with respect to the number of atoms and their potential for
massively parallel implementations for large-scale calculations Delley (1990);
Soler et al. (2002); Blum et al. (2009); Havu et al. (2009); Ren et al.
(2012); Enkovaara et al. (2010); Mohr et al. (2014); Frisch et al. (2009);
Shang et al. (2011). Unlike plane wave method, when constructing HFX matrix
within LCAO method, we must first calculate the ERIs via the atomic orbitals.
There are currently two types of atomic orbits that are most commonly used.
The first is gaussian type orbital(GTO), as adopted in GaussianFrisch et al.
(2009) and CRYSTALDovesi et al. (2006), its advantage is to calculate ERIs
analytically. The second is numerical atomic orbital (NAO), which is adopted
in SIESTASoler et al. (2002), DMOLDelley (1990), OPENMXOzaki (2003), et al..
The advantage of NAO is its strict locality, which naturally leads to lower
order scaling of computational time versus system size. In order to take
advantages of both types of atomic orbitals, we have proposed a new scheme
called NAO2GTOShang et al. (2011), in which GTO can be used for analytical
computation of ERIs in a straightforward and efficient way, while NAO can be
employed to set the strict cutoff for atomic orbitals. After employing several
ERI screening techniques, the construction of HFX matrix can be very efficient
and scale linearlyShang et al. (2011); Qin et al. (2014).

Parallelization of HFX matrix construction faces two major problems of load
imbalance and high communication cost. The load imbalance arises from the
irregularity of the independent tasks available in the computation, which is
due to the “screening” procedure and different types of shell quartets
distributed among processes. The high communication cost is from
interprocessor communication of the density and/or HFX matrices, which is
associated with the data access pattern. It is well known that NWChemValiev et
al. (2010) and CP2K/QuickstepVandeVondele et al. (2005) are the most
outstanding softwares in the field of high performance parallel quantum
chemical computing, and both of them use GTOs to construct HFX matrix. In
NWChem, the parallelization of HFX matrix construction is based on a static
partitioning of work followed by a work stealing phaseLiu et al. (2014); Chow
et al. (2015). The tasks are statically partitioned throughout the set of
shell (or atom) quartets, and then the work stealing phase acts to polish the
load balance. As a result, this parallel implementation gives very good
parallel scalability of Hartree–Fock calculationsLiu et al. (2014). In
CP2K/Quickstep, the HFX parallelization strategy is to replicate the global
density and HFX matrix on each MPI process in order to reduce communication. A
load balance optimization based on simulated annealing and a binning procedure
to coarse grain the load balancing problem have been developedGuidon et al.
(2008). However, this approach may limit both system size and ultimately
parallel scalability.

As the ERIs calculation is the most computationally demanding step in the
NAO2GTO scheme, the development of the new parallel algorithms is of
particular importance. Previously, for codes using localized atomic orbitals,
the parallelization of ERIs are mainly implemented to treat finite, isolated
systems Schmidt et al. (1993); Alexeev et al. (2002); Liu et al. (2014); Chow
et al. (2015), but only a few literature reports exist for the treatment of
periodic boundary conditions with such basis sets Bush et al. (2011); Guidon
et al. (2008), in which the Order-N screening for the ERIs calculations has
not been considered. The purpose of this work is to present the static
parallel distribution algorithms for the NAO2GTO schemeShang et al. (2011)
with Order-N performance in HONPAS codeQin et al. (2014). In our approaches,
the calculations of ERIs are not replicated, but are distributed over CPU
cores, as a result, both the memory and the CPU requirements of the ERIs
calculation are paralleled. The efficiency and scalability of these algorithms
are demonstrated by benchmark timings in periodic solid system with 64 silicon
atoms in the unit cell.

The outline of this paper is as follows: In Section 2, we begin with a
description of the theory of hybrid functionals. In Section 3, we describe the
detail implementation of our parallel distribution . In Section 4, we present
the benchmark results and the efficiency of our scheme.

##  2 Fundamental Theoretical Framework

Before addressing the parallel algorithms, we recall the basic equations used
in this work. A spin-unpolarized notation is used throughout the text for the
sake of simplicity, but a formal generalization to the collinear spin case is
straightforward. In Kohn-Sham DFT, the total-energy functional is given as

$E_{\textrm{\scriptsize KS}}=T_{\textrm{\scriptsize
s}}[n]+E_{\textrm{\scriptsize ext}}[n]+E_{\textrm{\scriptsize
H}}[n]+E_{\textrm{\scriptsize xc}}[n]+E_{\textrm{\scriptsize nuc-nuc}}\;.$ (1)

Here, $n(\mathbf{r})$ is the electron density, $T_{\textrm{\scriptsize s}}$ is
the kinetic energy of non-interacting electrons, while $E_{\textrm{\scriptsize
ext}}$ is external energy stemming from the electron-nuclear attraction,
$E_{\textrm{\scriptsize H}}$ is the Hartree energy, $E_{\textrm{\scriptsize
xc}}$ is the exchange-correlation energy, and $E_{\textrm{\scriptsize nuc-
nuc}}$ is the nucleus-nucleus repulsion energy.

The ground state electron density $n_{0}(\mathbf{r})$ (and the associated
ground state total energy) is obtained by variationally minimizing Eq. (1)
under the constraint that the number of electrons $N_{e}$ is conserved. This
yields the chemical potential $\mu=\delta E_{KS}/\delta n$ of the electrons
and the Kohn-Sham single particle equations

$\hat{h}_{\textrm{\scriptsize KS}}\psi_{i}=\left[\hat{t}_{\textrm{\scriptsize
s}}+v_{\textrm{\scriptsize ext}}(r)+v_{\textrm{\scriptsize
H}}+v_{\textrm{\scriptsize xc}}\right]\psi_{i}=\epsilon_{p}\psi_{i}$ (2)

for the Kohn-Sham Hamiltonian $\hat{h}_{\textrm{\scriptsize KS}}$. In Eq. (2),
$\hat{t}_{\textrm{\scriptsize s}}$ denotes the kinetic energy operator,
$v_{\textrm{\scriptsize ext}}$ the external potential, $v_{H}$ the Hartree
potential, and $v_{xc}$ the exchange-correlation potential. Solving Eq. (2)
yields the Kohn-Sham single particle states $\psi_{p}$ and their eigenenergies
$\epsilon_{p}$. The single particle states determine the electron density via

$n(\mathbf{r})=\sum_{i}f_{i}|\psi_{i}|^{2},$ (3)

in which $f_{i}$ denotes the Fermi-Dirac distribution function.

To solve Eq. (2) in numerical implementations, the Kohn-Sham states are
expanded in a finite basis set. For periodic systems, the crystalline orbital
$\psi_{i}(\mathbf{k,r})$ normalized in all space is a linear combination of
Bloch functions $\phi_{\mu}(\mathbf{k,r})$, defined in terms of atomic
orbitals $\chi_{\mu}^{\mathbf{R}}(\mathbf{r})$.

$\psi_{i}(\mathbf{k,r})=\sum_{\mu}c_{\mu,i}(\mathbf{k})\phi_{\mu}(\mathbf{k,r})$
(4)
$\phi_{\mu}(\mathbf{k,r})=\dfrac{1}{\sqrt{N}}\sum_{\mathbf{R}}\chi_{\mu}^{\mathbf{R}}(\mathbf{r})\mathbf{e}^{i\mathbf{k}\cdotp(\mathbf{R+r_{\mu}})}$
(5)

where the Greek letter $\mu$ is the index of atomic orbitals, $i$ is the
suffix for different bands, $\mathbf{R}$ is the origin of a unit cell, N is
the number of unit cells in the system.
$\chi_{\mu}^{\mathbf{R}}(\mathbf{r})=\chi_{\mu}(\mathbf{r-R-r_{\mu}})$ is the
$\mu$-th atomic orbital, whose center is displaced from the origin of the unit
cell at $\mathbf{R}$ by $\mathbf{r}_{\mu}$. $c_{\mu,i}(\mathbf{k})$ is the
wave function coefficient, which is obtained by solving the following
equation,

$H(\mathbf{k})c(\mathbf{k})=E(\mathbf{k})S(\mathbf{k})c(\mathbf{k})$ (6)
$[H(\mathbf{k})]_{\mu\nu}=\sum_{\mathbf{R}}H_{\mu\nu}^{\mathbf{R}}\mathbf{e}^{i\mathbf{k}\cdotp(\mathbf{R+r_{\nu}-r_{\mu}})}$
(7)
$H^{\mathbf{R}}_{\mu\nu}=<\chi_{\mu}^{\mathbf{0}}|\hat{H}|\chi_{\nu}^{\mathbf{R}}>$
(8)
$[S(\mathbf{k})]_{\mu\nu}=\sum_{\mathbf{R}}S_{\mu\nu}^{\mathbf{R}}\mathbf{e}^{i\mathbf{k}\cdotp(\mathbf{R+r_{\nu}-r_{\mu}})}$
(9)
$S^{\mathbf{R}}_{\mu\nu}=<\chi_{\mu}^{\mathbf{0}}|\chi_{\nu}^{\mathbf{R}}>$
(10)

In Eq. (8), $H^{\mathbf{R}}_{\mu\nu}$ is a matrix element of the one-electron
Hamiltonian operator $\hat{H}$ between the atomic orbital $\chi_{\mu}$ located
in the central unit cell $0$ and $\chi_{\nu}$ located in the unit cell
$\mathbf{R}$.

It should be noted that, the exchange-correlation potential $v^{xc}$ is a
local and periodic in semi-local DFT, while in Hartree-Fock and hybrid
functionals, the Hartree-Fock exchange (HFX) potential matrix element is
defined as:

$[V^{X}]_{\mu\lambda}^{\mathbf{G}}=-\frac{1}{2}\sum_{\nu\sigma}\sum_{\mathbf{N,H}}P_{\nu\sigma}^{\mathbf{H-N}}\mathbf{[(\chi_{\mu}^{0}\chi_{\nu}^{N}|\chi_{\lambda}^{G}\chi_{\sigma}^{H})]}$
(11)

where $\mathbf{G}$, $\mathbf{N}$, and $\mathbf{H}$ represent different unit
cells. The density matrix element $P_{\nu\sigma}^{\mathbf{N}}$ is computed by
an integration over the Brillouin zone (BZ)

$P_{\nu\sigma}^{\mathbf{N}}=\sum_{j}\int_{BZ}c_{\nu,j}^{*}(\mathbf{k})c_{\sigma,j}(\mathbf{k})\theta(\epsilon_{F}-\epsilon_{j}(\mathbf{k}))\mathbf{e}^{i\mathbf{k}\cdotp\mathbf{N}}d\mathbf{k}$
(12)

where $\theta$ is the step function, $\epsilon_{F}$ is the fermi energy and
$\epsilon_{j}(\mathbf{k})$ is the $j$-th eigenvalue at point $\mathbf{k}$.

In order to calculate the following ERI in Eq. (11)

$\mathbf{(\chi_{\mu}^{0}\chi_{\nu}^{N}|\chi_{\lambda}^{G}\chi_{\sigma}^{H})=\int\int\frac{\chi_{\mu}^{0}(r)\chi_{\nu}^{N}(r)\chi_{\lambda}^{G}(r^{\prime})\chi_{\sigma}^{H}(r^{\prime})}{|r-r^{\prime}|}drdr^{\prime}}$
(13)

we use NAO2GTO scheme described in the following section.

Following the flowchart in Fig.1, the NAO2GTO scheme is to firstly fit the NAO
with GTOs, and then to calculate ERIs analytically with fitted GTOs. A NAO is
a product of a numerical radial function and a spherical harmonic

$\phi_{Ilmn}(\mathbf{r})=\varphi_{Iln}(r)\mathbf{Y}_{lm}(\hat{r})$ (14)

The radial part of the numerical atomic orbital $\varphi_{Iln}(r)$ is
calculated by the following equation:

$(-\frac{1}{2}\mathbf{\frac{1}{r}}\frac{\mathrm{d^{2}}}{\mathrm{d}r^{2}}r+\frac{l(l+1)}{2r^{2}}+V(r)+V_{cut})\varphi_{Iln}(r)=\epsilon_{l}\varphi_{Iln}(r)$
(15)

where $V(r)$ denotes the electrostatic potential for orbital
$\varphi_{Iln}(r)$, and $V_{cut}$ ensures a smooth decay of each radial
function which is strictly zero outside a confining radius $r_{cut}$.

##  3 METHODS

##  4 NAO2GTO scheme

Figure 1: The flowchart of the NAO2GTO scheme in HONPAS.

In our NAO2GTO approach, the radial part of the NAO $\varphi_{Iln}(r)$ is
fitted by the summation of several GTOs, denoted as $\chi(r)$

$\chi(r)\equiv\sum_{m}{D_{m}r^{l}\exp(-\alpha_{m}r^{2})}$ (16)

Parameters $\alpha_{m}$ and $D_{m}$ are determined by minimizing the residual
sum of squares of the difference

$\sum_{i}{(\chi(r_{i})/r_{i}^{l}-\varphi_{Iln}(r_{i})/r_{i}^{l})^{2}}$ (17)

In practice of the solid system calculation, too diffused basis set may cause
convergence problemGuidon et al. (2008), as a result the exponents smaller
than 0.10 are usually not needed, and we made a constraint, i.e.
($\alpha>0.1$) during our minimal search. First we use constrained genetic
algorithmGoldberg (1989); Conn et al. (1991) to do a global search for initial
guess and then do a constrained local minimal search using trust-region-
reflective algorithm, which is a subspace trust-region method and is based on
the interior-reflective Newton method described in Coleman and Li (1994) and
Coleman and Li (1996). Each iteration involves the approximate solution of a
large linear system using the method of preconditioned conjugate gradients. We
make N ($N>500$) global searches to make sure to have a global minimal.

Algorithm 1 The algorithm of NAO2GTO fitting scheme.

for $iter=1$ to $N$ do

constrained genetic algorithm get initial $\alpha^{iter}_{m}$ and
$D^{iter}_{m}$

constrained local minimal search to get $\alpha^{iter}_{m}$ and $D^{iter}_{m}$

$err=\sum_{i}{[\sum_{m}{D_{m}\exp(-\alpha_{m}r_{i}^{2})}-\varphi_{Iln}(r_{i})/r_{i}^{l}]^{2}}$

if $iter=1$ .or. $err<best\\_err$ then

$best\\_err=err$

$\alpha_{m}=\alpha^{iter}_{m}$ and $D_{m}=D^{iter}_{m}$

end if

end for

Using the above fitting parameters, we build all ERIs that needed for the HFX.
In our implementation, the full permutation symmetry of the ERIs has been
considered for solid systems:

$(\mu^{\mathbf{0}}\nu^{\mathbf{H}}|\lambda^{\mathbf{G}}\sigma^{\mathbf{N}})=(\mu^{\mathbf{0}}\nu^{\mathbf{H}}|\sigma^{\mathbf{N}}\lambda^{\mathbf{G}})=$
$(\nu^{\mathbf{0}}\mu^{\mathbf{-H}}|\lambda^{\mathbf{G-H}}\sigma^{\mathbf{N-H}})=(\nu^{\mathbf{0}}\mu^{\mathbf{-H}}|\sigma^{\mathbf{N-H}}\lambda^{\mathbf{G-H}})=$
$(\lambda^{\mathbf{0}}\sigma^{\mathbf{N-G}}|\mu^{\mathbf{-G}}\nu^{\mathbf{H-G}})=(\lambda^{\mathbf{0}}\sigma^{\mathbf{N-G}}|\nu^{\mathbf{H-G}}\mu^{\mathbf{-G}})=$
$(\sigma^{\mathbf{0}}\lambda^{\mathbf{G-N}}|\mu^{\mathbf{-N}}\nu^{\mathbf{H-N}})=(\sigma^{\mathbf{0}}\lambda^{\mathbf{G-N}}|\nu^{\mathbf{H-N}}\mu^{\mathbf{-N}})$
(18)

In this way, we save a factor of 8 in the number of integrals that have to be
calculated.

When calculating the ERIs with GTOs, the atomic orbitals are grouped into
shells according to the angular momentum. The list need to be distributed in
parallel is in fact the shell quartet. For a shell with a angular momentum of
$l$, the number of atomic orbital basis functions is $2l+1$, so in a shell
quartet integral $(IJ|KL)$, we calculate in total
$(2l_{I}+1)(2l_{J}+1)(2l_{K}+1)(2l_{L}+1)$ atomic basis orbital ERIs together.
As a result, the computational expense is strongly dependent on the angular
momenta of the shell quartet. It is a challenge to distribute these shell ERIs
not only in number but also considering the time-weight.

In our NAO2GTO scheme, two shell pair lists (list-IJ and list-KL) are firstly
preselected according to Schwarz screeningHäser and Ahlrichs (1989),

$|(\mu\nu|\lambda\sigma)|\leqslant\sqrt{(\mu\nu|\mu\nu)(\lambda\sigma|\lambda\sigma)}$
(19)

and only the shell list indexes with $(IJ|IJ)>\tau$ or $(KL|KL)>\tau$ (here
$\tau$ is the drop tolerance) are stored. As shown in Eq. (11), the first
index I runs only within the unit cell, while the indexes (J,K,L) run over the
whole supercell, so the list-IJ is smaller than the list-KL. Then in the ERIs
calculations, the loops run over these two shell lists.

Then before the calculation of every ERI, we use Schwarz inequality Eq. (19)
again to estimate a rigorous upper bound, that only the ERIs with non-
negligible contributions are calculated, we note this screening method as
Schwarz screening. Using the exponential decay of the charge distributions,
the Schwarz screening reduces the total number of ERIs to be computed from
$O(N^{4})$ to $O(N^{2})$. The Schwarz screening tolerance is set to 10-5 in
the following calculation. In addition, the NAO screening is also adopted as
the NAO is strictly truncatedShang et al. (2011). The NAO screening is safe in
the calculation of the short-range ERI because in this case the Hartree-Fock
exchange (HFX) Hamiltonian matrix is sparse due to the screened Coulomb
potentialIzmaylov et al. (2006). As a result, we store this HFX Hamiltonian
with a sparse matrix data structure.

In practice, it also should be noted that as the angular part of the NAOs is
spherical harmonic while the GTOs are Cartesian Gaussians, we need to make a
transformation between Cartesian and Spherical harmonic functions. The
difference between these two harmonic functions is the number of atomic
orbitals including in the shells whose angular momentum are larger than 1. For
example, a d-shell has 5 Spherical orbitals, but have 6 Cartesian orbitals. A
detailed transformation method can be find in Ref.Schlegel and Frisch (1995).

###  4.1 Parallel schemes

The ERIs have four indexed that can be paralleled. One possible parallel
scheme to make distribution of just one shell index, however, as the index
number in one shell is too small to make distribution over CPU cores, such
shell-one distribution may cause serious load imbalance.

The other parallel scheme is to make distribution for shell-pair (list-KL). It
is a straightforward way to parallelize the ERIs as in practice we loop over
two pair lists. However, although the shell-pair can be distributed evenly
before ERIs calculations, ERI-screening is needed during the ERIs
calculations, which also causes load imbalance. The practical implementation
of the described formalism closely follows the flowchart shown in Algorithm 2.
After the building list-KL is completed, we distributed it into CPU cores with
list-KL-local at every cores. Then in the following ERIs screening and
calculation, only loops over list-KL-local is needed for every core. The
advantage of this scheme is that it naturally bypasses the the communication
process, and every CPU core only go over and compute its assigned list-KL-
local. However, although the list-KL is distributed evenly over processors,
the ERI screening is located after the parallelization, which causes different
number of shell ERIs that need to be calculated in every processor. Such
different number of shell ERIs makes load imbalance.

In order to achieve load balance, distribution of individual shell-
quartet$(IJ|KL)$ after the ERI screening process a possible choice. Even if
the computational time is nonuniformity in the case of the different shell
type, this distribution can also yield an even time over CPU cores because of
its smallest distribution chunks. The practical implementation of the shell-
quartet algorithm is shown in Algorithm 3. Every CPU core will go over the
global pair lists( list-IJ and list-KL), and make the ERI screening to
determine which ERIs are needed to be computed. Then a global counter is set
to count the number of computed ERIs, this counter is distributed over CPU
cores to make sure the number of calculated ERIs is evenly distributed. The
disadvantage of this algorithm is that every processor need to make the whole
ERI screening, while in the above parallel-pair algorithm, only the ERI
screening in its local lists is needed. Such globally calculated ERIs
screening decreases the parallel efficiency.

Algorithm 2 Flowchart of the parallel-pair algorithms for ERIs. Here the
shell-pair-list-KL-local means to distribute the shell-list-KL over CPU cores
at the beginning. The description of Schwarz screening and NAO screening are
in the text.

get shell-pair-list-KL-local

for list-IJ in shell-pair-list-IJ do

for list-KL in shell-pair-list-KL-local do

if Schwarz screening.and. NAO screening then

compute shell ERI $(IJ|KL)$

end if

end for

end for

Algorithm 3 Flowchart of the parallel-quartet algorithms for ERIs. Here N
refers to the total number of the CPU cores, current-core refers to the index
of the current processor. The description of Schwarz screening and NAO
screening are in the text.

for list-IJ in shell-pair-list-IJ do

for list-KL in shell-pair-list-KL do

if Schwarz screening.and. NAO screening then

i ++

if i mod N eq current-core then

compute shell ERI $(IJ|KL)$

end if

end if

end for

end for

##  5 RESULTS AND DISCUSSION

In order to demonstrate the performance of the above two static parallel
schemes, we use silicon bulk contained 64 atoms in the unit cell as a test
case as shown in Fig. 2. Norm-conserving pseudopotentials generated with the
Troullier-MartinsTroullier and Martins (1991) scheme, in fully separable form
developed by Kleiman and ByladerKleinman and Bylander (1982), are used to
represent interaction between core ion and valence electrons. The screened
hybrid functional HSE06Krukau et al. (2006) was used in the following
calculations. Both single-zeta (SZ) contained s and p shells and double-zeta
plus polarization (DZP) basis set contained s,p and d shells are considered.
All calculations were carried out on Tianhe-2 supercomputer located in
National Supercomputer Center in Guangzhou, China, the configuration of the
machine is shown in Table 1. The Intel Math Kernel Library(version 10.0.3.020)
is used in the calculations.

Table 1: The parameters for each node of Tianhe-2. Component | Value
---|---
CPU | Intel(R) Xeon(R) CPU E5-2692
Freq. (GHz) | 2.2
Cores | 12

For the parallel-pair and parallel-quartet algorithms, which are fully
parallelized and involve no communication, load imbalance is one of factors
that may affect the parallel efficiency. To examine the load balance, the
timing at every cores are shown in Fig. 3-Fig. 6. It is clearly shown that the
parallel-pair algorithm (red line) is load imbalance, for SZ basis set, the
time difference between cores is around 10% in 12 cores (Fig. 3) case and
around 80% in 192 cores (Fig. 4). For DZP basis set, d shells have been
considered, which caused more serious load imbalance, the time difference
between cores is even around 22% in 12 cores (Fig.5) and around 100% in 192
cores (Fig.6). On the other hand, the load balance in parallel-quartet
algorithm is quite well, the time difference between cores is within 1%.
However, in this algorithm, as the ERIs screening part is made for the whole
ERIs by all the CPU cores. which is a constant time even with the increasing
CPU cores, there are replicate calculations for the ERIs screening which
decrease the parallel efficiency. As shown in Fig.4), for small basis set at
large number of CPU cores, the average CPU time of the parallel-quartet is
around twice as the parallel-pair algorithm.

| basis set | shells | NAOs
---|---|---|---
Si64 | SZ | 128 | 256
Si64 | DZP | 320 | 832

Figure 2: The silicon bulk contained 64 atoms in the unit cell that used as benchmark system in this work. In the upper table, the number of shells as well as the number of the NAO basis functions for different basis sets are listed. Figure 3: The load balance for bulk silicon supercell with 64 atoms using SZ basis set at 12 CPU cores. Figure 4: The load balance for bulk silicon supercell with 64 atoms using SZ basis set at 192 CPU cores. Figure 5: The load balance for bulk silicon supercell with 64 atoms using DZP basis set at 12 CPU cores. Figure 6: The load balance for bulk silicon supercell with 64 atoms using DZP basis set at 192 CPU cores. Table 2: The CPU time (in seconds) for the calculation of the ERIs of Si 64 system with SZ basis set using different parallel schemes.  cores | parallel-pair | parallel-quartet
---|---|---
12 | 110.9 | 108.2
24 | 55.6 | 57.0
96 | 15.5 | 18.2
192 | 8.6 | 11.6
Table 3: The CPU time (in seconds) for the calculation of the ERIs of Si 64 system with DZP basis set using different parallel schemes.  cores | parallel-pair | parallel-quartet
---|---|---
12 | 1645.1 | 1572.9
24 | 904.0 | 806.1
96 | 251.3 | 225.9
192 | 129.6 | 128.0
Figure 7:  (Color online) Parallel Speedups and efficiency for ERIs
calculation formation using different parallel schemes. Speedups were obtained
on Tianhe-2 for bulk silicon supercell with 64 atoms using SZ basis set. The
speedup is referenced to a run on 12 CPUs. Figure 8:  (Color online) Parallel
Speedups and efficiency for ERIs calculation using different parallel schemes.
Speedups were obtained on Tianhe-2 for bulk silicon supercell with 64 atoms
using DZP basis set The speedup is referenced to a run on 12 CPUs.

Such global ERI screening in parallel-quartet algorithm also contributes
significantly to lowering the parallel speedup and efficiency for SZ basis
set, as shown in Table2 and Fig.7. The parallel efficiency is only 58% at 192
CPU core for parallel-quartet while holds around 80% for parallel-pair
algorithm.

For DZP basis set, the load imbalance caused by d shells become another factor
of lowering the parallel efficiency, so in this case, as shown in Table2 and
Fig.8, the parallel speedup and efficiency are similar for both parallel-pair
and parallel-quartet algorithms, that around 80% at 192 CPU cores.

##  6 CONCLUSIONS

In summary, we have shown our two static parallel algorithms for the ERIs
calculations in NAO2GTO method. We have also analyzed the performance of these
two parallel algorithms for their load balance and parallel efficiency. On the
basis of our results, the static distribution of ERI shell pairs, produces
load imbalance that causes the efficiency to decrease, limiting the number of
CPU cores that can be utilized. On the other hand, the static distribution of
ERI shell quartet can yield very high load balance, however, because the need
of the global ERI screening calculation, the parallel efficiency has been
dramatically reduced for small basis set. We have also tried another static
method that firstly create a need-to-calculate ERIs list by considering all
the screening methods as well as the eight-fold permutational symmetry and
secondly distribute the ERIs in the need-to-calculate list over a number of
processes. However, we find the time to build the need-to-calculate ERIs list
is even larger than the global ERI screening calculation. On the next step, we
need to distribute the ERI screening calculation while keep the load balance
of the ERI calculation, and a dynamic distribution could enables load balance
with little loss of efficiency.

This work was supported by the National Key Research $\&$ Development Program
of China (2017YFB0202302), the Special Fund for Strategic Pilot Technology of
Chinese Academy of Sciences (XDC01040000), the National Natural Science
Foundation of China (61502450), the National Natural Science Foundation of
China (21803066), Research Start-Up Grants (KY2340000094) from University of
Science and Technology of China, and Chinese Academy of Sciences Pioneer
Hundred Talents Program. The authors thank the Tianhe-2 Supercomputer Center
for computational resources.

## References

  * Alexeev et al. (2002) Alexeev Y, Kendall RA and Gordon MS (2002) The distributed data SCF.  _Computer Physics Communications_ 143(1): 69–82.  10.1016/s0010-4655(01)00439-8.  URL https://doi.org/10.1016%2Fs0010-4655%2801%2900439-8.
  * Becke (1993) Becke AD (1993) Density-functional thermochemistry. III. the role of exact exchange.  _The Journal of Chemical Physics_ 98(7): 5648–5652.  10.1063/1.464913.  URL https://doi.org/10.1063%2F1.464913.
  * Blum et al. (2009) Blum V, Gehrke R, Hanke F, Havu P, Havu V, Ren X, Reuter K and Scheffler M (2009) Ab initio molecular simulations with numeric atom-centered orbitals.  _Comput. Phys. Commun._ 180(11): 2175–2196.  10.1016/j.cpc.2009.06.022.  URL http://linkinghub.elsevier.com/retrieve/pii/S0010465509002033.
  * Burant et al. (1996) Burant JC, Scuseria GE and Frisch MJ (1996) A linear scaling method for hartree–fock exchange calculations of large molecules.  _The Journal of Chemical Physics_ 105(19): 8969–8972.  10.1063/1.472627.  URL https://doi.org/10.1063%2F1.472627.
  * Bush et al. (2011) Bush IJ, Tomic S, Searle BG, Mallia G, Bailey CL, Montanari B, Bernasconi L, Carr JM and Harrison NM (2011) Parallel implementation of the ab initio CRYSTAL program: electronic structure calculations for periodic systems.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 467(2131): 2112–2126.  10.1098/rspa.2010.0563.  URL https://doi.org/10.1098%2Frspa.2010.0563.
  * Chow et al. (2015) Chow E, Liu X, Smelyanskiy M and Hammond JR (2015) Parallel scalability of hartree–fock calculations.  _The Journal of Chemical Physics_ 142(10): 104103.  10.1063/1.4913961.  URL https://doi.org/10.1063%2F1.4913961.
  * Coleman and Li (1994) Coleman TF and Li Y (1994) On the convergence of interior-reflective newton methods for nonlinear minimization subject to bounds.  _Mathematical Programming_ 67(1-3): 189–224.  10.1007/bf01582221.  URL https://doi.org/10.1007%2Fbf01582221.
  * Coleman and Li (1996) Coleman TF and Li Y (1996) An interior trust region approach for nonlinear minimization subject to bounds.  _SIAM Journal on Optimization_ 6(2): 418–445.  10.1137/0806023.  URL https://doi.org/10.1137%2F0806023.
  * Conn et al. (1991) Conn AR, Gould NIM and Toint P (1991) A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds.  _SIAM Journal on Numerical Analysis_ 28(2): 545–572.  10.1137/0728030.  URL https://doi.org/10.1137%2F0728030.
  * Delhalle and Calais (1987) Delhalle J and Calais JL (1987) Direct-space analysis of the hartree-fock energy bands and density of states for metallic extended systems.  _Physical Review B_ 35(18): 9460–9466.  10.1103/physrevb.35.9460.  URL https://doi.org/10.1103%2Fphysrevb.35.9460.
  * Delley (1990) Delley B (1990) An all-electron numerical method for solving the local density functional for polyatomic molecules.  _The Journal of Chemical Physics_ 92(1): 508–517.  10.1063/1.458452.  URL https://doi.org/10.1063%2F1.458452.
  * Dovesi et al. (2006) Dovesi R, Saunders V, Roetti C, Orlando R, Zicovich-Wilsona CM, Pascale F, Civalleri B, Doll K, Harriso NM, Bush I, D’Arco P and Llunell M (2006) Crystal06 user’s manual.  University of Torino, Torino, 2006.
  * Enkovaara et al. (2010) Enkovaara J, Rostgaard C, Mortensen JJ, Chen J, Dułak M, Ferrighi L, Gavnholt J, Glinsvad C, Haikola V, Hansen HA, Kristoffersen HH, Kuisma M, Larsen AH, Lehtovaara L, Ljungberg M, Lopez-Acevedo O, Moses PG, Ojanen J, Olsen T, Petzold V, Romero NA, Stausholm-Mø ller J, Strange M, Tritsaris GA, Vanin M, Walter M, Hammer B, Häkkinen H, Madsen GKH, Nieminen RM, Nø rskov JK, Puska M, Rantala TT, Schiø tz J, Thygesen KS and Jacobsen KW (2010) Electronic structure calculations with GPAW: a real-space implementation of the projector augmented-wave method.  _J. Phys. Condens. Matter_ 22(25): 253202.  URL http://stacks.iop.org/0953-8984/22/i=25/a=253202.
  * Frisch et al. (2009) Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, Scalmani G, Barone V, Mennucci B, Petersson GA, Nakatsuji H, Caricato M, Li X, Hratchian HP, Izmaylov AF, Bloino J, Zheng G, Sonnenberg JL, Hada M, Ehara M, Toyota K, Fukuda R, Hasegawa J, Ishida M, Nakajima T, Honda Y, Kitao O, Nakai H, Vreven T, Montgomery JA, Peralta JE, Ogliaro F, Bearpark M, Heyd JJ, Brothers E, Kudin KN, Staroverov VN, Kobayashi R, Normand J, Raghavachari K, Rendell A, Burant JC, Iyengar SS, Tomasi J, Cossi M, Rega N, Millam JM, Klene M, Knox JE, Cross JB, Bakken V, Adamo C, Jaramillo J, Gomperts R, Stratmann RE, Yazyev O, Austin AJ, Cammi R, Pomelli C, Ochterski JW, Martin RL, Morokuma K, Zakrzewski VG, Voth GA, Salvador P, Dannenberg JJ, Dapprich S, Daniels AD, Farkas, Foresman JB, Ortiz JV, Cioslowski J and Fox DJ (2009) Gaussian 09, Revision B.01.
  * Gell-Mann and Brueckner (1957) Gell-Mann M and Brueckner KA (1957) Correlation energy of an electron gas at high density.  _Physical Review_ 106(2): 364–368.  10.1103/physrev.106.372.  URL https://doi.org/10.1103%2Fphysrev.106.372.
  * Goldberg (1989) Goldberg DE (1989) _Genetic Algorithms in Search, Optimization and Machine Learning_.  1st edition. Boston, MA, USA: Addison-Wesley Longman Publishing Co., Inc.  ISBN 0201157675.
  * Gonze et al. (2002) Gonze X, Beuken JM, Caracas R, Detraux F, Fuchs M, Rignanese GM, Sindic L, Verstraete M, Zerah G, Jollet F, Torrent M, Roy A, Mikami M, Ghosez P, Raty JY and Allan DC (2002) First-principles computation of material properties: The ABINIT software project.  _Comput. Mater. Sci._ 25(3): 478–492.  10.1016/S0927-0256(02)00325-7.  URL https://doi.org/10.1016%2Fs0927-0256%2802%2900325-7.
  * Gonze et al. (2016) Gonze X, Jollet F, Abreu Araujo F, Adams D, Amadon B, Applencourt T, Audouze C, Beuken JM, Bieder J, Bokhanchuk A, Bousquet E, Bruneval F, Caliste D, Côté M, Dahm F, Da Pieve F, Delaveau M, Di Gennaro M, Dorado B, Espejo C, Geneste G, Genovese L, Gerossier A, Giantomassi M, Gillet Y, Hamann D, He L, Jomard G, Laflamme Janssen J, Le Roux S, Levitt A, Lherbier A, Liu F, Lukačević I, Martin A, Martins C, Oliveira M, Poncé S, Pouillon Y, Rangel T, Rignanese GM, Romero A, Rousseau B, Rubel O, Shukri A, Stankovski M, Torrent M, Van Setten M, Van Troeye B, Verstraete M, Waroquiers D, Wiktor J, Xu B, Zhou A and Zwanziger J (2016) Recent developments in the ABINIT software package.  _Comput. Phys. Commun._ 10.1016/j.cpc.2016.04.003.  URL http://linkinghub.elsevier.com/retrieve/pii/S0010465516300923.
  * Guidon et al. (2010) Guidon M, Hutter J and VandeVondele J (2010) Auxiliary density matrix methods for hartree-fock exchange calculations.  _Journal of Chemical Theory and Computation_ 6(8): 2348–2364.  10.1021/ct1002225.  URL https://doi.org/10.1021%2Fct1002225.
  * Guidon et al. (2008) Guidon M, Schiffmann F, Hutter J and VandeVondele J (2008) Ab initio molecular dynamics using hybrid density functionals.  _The Journal of Chemical Physics_ 128(21): 214104.  10.1063/1.2931945.  URL https://doi.org/10.1063%2F1.2931945.
  * Havu et al. (2009) Havu V, Blum V, Havu P and Scheffler M (2009) Efficient integration for all-electron electronic structure calculation using numeric basis functions.  _J. Comput. Phys._ 228(22): 8367–8379.  10.1016/j.jcp.2009.08.008.  URL http://linkinghub.elsevier.com/retrieve/pii/S0021999109004458.
  * Heyd et al. (2003) Heyd J, Scuseria GE and Ernzerhof M (2003) Hybrid functionals based on a screened coulomb potential.  _The Journal of Chemical Physics_ 118(18): 8207–8215.  10.1063/1.1564060.  URL https://doi.org/10.1063%2F1.1564060.
  * Heyd et al. (2006) Heyd J, Scuseria GE and Ernzerhof M (2006) Erratum: “hybrid functionals based on a screened coulomb potential” [j. chem. phys. 118, 8207 (2003)].  _The Journal of Chemical Physics_ 124(21): 219906.  10.1063/1.2204597.  URL https://doi.org/10.1063%2F1.2204597.
  * Hohenberg and Kohn (1964) Hohenberg P and Kohn W (1964) Inhomogeneous electron gas.  _Physical Review_ 136(3B): B864–B871.  10.1103/physrev.136.b864.  URL https://doi.org/10.1103%2Fphysrev.136.b864.
  * Häser and Ahlrichs (1989) Häser M and Ahlrichs R (1989) Improvements on the direct SCF method.  _Journal of Computational Chemistry_ 10(1): 104–111.  10.1002/jcc.540100111.  URL https://doi.org/10.1002%2Fjcc.540100111.
  * Izmaylov et al. (2006) Izmaylov AF, Scuseria GE and Frisch MJ (2006) Efficient evaluation of short-range hartree-fock exchange in large molecules and periodic systems.  _The Journal of Chemical Physics_ 125(10): 104103.  10.1063/1.2347713.  URL https://doi.org/10.1063/1.2347713.
  * Janesko et al. (2009) Janesko BG, Henderson TM and Scuseria GE (2009) Screened hybrid density functionals for solid-state chemistry and physics.  _Phys. Chem. Chem. Phys._ 11(3): 443–454.  10.1039/b812838c.  URL https://doi.org/10.1039%2Fb812838c.
  * Kleinman and Bylander (1982) Kleinman L and Bylander DM (1982) Efficacious form for model pseudopotentials.  _Physical Review Letters_ 48(20): 1425–1428.  10.1103/physrevlett.48.1425.  URL https://doi.org/10.1103%2Fphysrevlett.48.1425.
  * Kohn and Sham (1965) Kohn W and Sham LJ (1965) Self-consistent equations including exchange and correlation effects.  _Physical Review_ 140(4A): A1133–A1138.  10.1103/physrev.140.a1133.  URL https://doi.org/10.1103%2Fphysrev.140.a1133.
  * Krukau et al. (2006) Krukau AV, Vydrov OA, Izmaylov AF and Scuseria GE (2006) Influence of the exchange screening parameter on the performance of screened hybrid functionals.  _The Journal of Chemical Physics_ 125(22): 224106.  10.1063/1.2404663.  URL https://doi.org/10.1063%2F1.2404663.
  * Liu et al. (2014) Liu X, Patel A and Chow E (2014) A new scalable parallel algorithm for fock matrix construction.  In: _2014 IEEE 28th International Parallel and Distributed Processing Symposium_. IEEE.  10.1109/ipdps.2014.97.  URL https://doi.org/10.1109%2Fipdps.2014.97.
  * Merlot et al. (2014) Merlot P, Izsák R, Borgoo A, Kjærgaard T, Helgaker T and Reine S (2014) Charge-constrained auxiliary-density-matrix methods for the hartree-fock exchange contribution.  _The Journal of Chemical Physics_ 141(9): 094104.  10.1063/1.4894267.  URL https://doi.org/10.1063%2F1.4894267.
  * Mohr et al. (2014) Mohr S, Ratcliff LE, Boulanger P, Genovese L, Caliste D, Deutsch T and Goedecker S (2014) Daubechies wavelets for linear scaling density functional theory.  _J. Chem. Phys._ 140(20): 204110.  10.1063/1.4871876.  URL http://www.ncbi.nlm.nih.gov/pubmed/24880269.
  * Monkhorst (1979) Monkhorst HJ (1979) Hartree-Fock density of states for extended systems.  _Physical Review B_ 20(4): 1504–1513.  10.1103/physrevb.20.1504.  URL https://doi.org/10.1103%2Fphysrevb.20.1504.
  * Mori-Sánchez et al. (2008) Mori-Sánchez P, Cohen AJ and Yang W (2008) Localization and delocalization errors in density functional theory and implications for band-gap prediction.  _Physical Review Letters_ 100(14).  10.1103/physrevlett.100.146401.  URL https://doi.org/10.1103%2Fphysrevlett.100.146401.
  * Ochsenfeld et al. (1998) Ochsenfeld C, White CA and Head-Gordon M (1998) Linear and sublinear scaling formation of Hartree–Fock-type exchange matrices.  _The Journal of Chemical Physics_ 109(5): 1663–1669.  10.1063/1.476741.  URL https://doi.org/10.1063%2F1.476741.
  * Ozaki (2003) Ozaki T (2003) Variationally optimized atomic orbitals for large-scale electronic structures.  _Physical Review B_ 67(15).  10.1103/physrevb.67.155108.  URL https://doi.org/10.1103%2Fphysrevb.67.155108.
  * Paier et al. (2009) Paier J, Diaconu CV, Scuseria GE, Guidon M, VandeVondele J and Hutter J (2009) Accurate hartree-fock energy of extended systems using large gaussian basis sets.  _Physical Review B_ 80(17).  10.1103/physrevb.80.174114.  URL https://doi.org/10.1103%2Fphysrevb.80.174114.
  * Paier et al. (2006) Paier J, Marsman M, Hummer K, Kresse G, Gerber IC and Ángyán JG (2006) Screened hybrid density functionals applied to solids.  _The Journal of Chemical Physics_ 124(15): 154709.  10.1063/1.2187006.  URL https://doi.org/10.1063%2F1.2187006.
  * Parr and Yang (1989) Parr RG and Yang W (1989) _Density Functional Theory of Atoms and Molecules_.  New York: Oxford University Press.
  * Polly et al. (2004) Polly R, Werner HJ, Manby FR and Knowles PJ (2004) Fast hartree–fock theory using local density fitting approximations.  _Molecular Physics_ 102(21-22): 2311–2321.  10.1080/0026897042000274801.  URL https://doi.org/10.1080%2F0026897042000274801.
  * Qin et al. (2014) Qin X, Shang H, Xiang H, Li Z and Yang J (2014) HONPAS: A linear scaling open-source solution for large system simulations.  _International Journal of Quantum Chemistry_ 115(10): 647–655.  10.1002/qua.24837.  URL https://doi.org/10.1002%2Fqua.24837.
  * Ren et al. (2012) Ren X, Rinke P, Blum V, Wieferink J, Tkatchenko A, Sanfilippo A, Reuter K and Scheffler M (2012) Resolution-of-identity approach to Hartree–Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions.  _New J. Phys._ 14(5): 053020.  10.1088/1367-2630/14/5/053020.  URL http://stacks.iop.org/1367-2630/14/i=5/a=053020?key=crossref.351b343783c2c1df1596219a941a74eb.
  * Schlegel and Frisch (1995) Schlegel HB and Frisch MJ (1995) Transformation between cartesian and pure spherical harmonic gaussians.  _International Journal of Quantum Chemistry_ 54(2): 83–87.  10.1002/qua.560540202.  URL https://doi.org/10.1002%2Fqua.560540202.
  * Schmidt et al. (1993) Schmidt MW, Baldridge KK, Boatz JA, Elbert ST, Gordon MS, Jensen JH, Koseki S, Matsunaga N, Nguyen KA, Su S, Windus TL, Dupuis M and Montgomery JA (1993) General atomic and molecular electronic structure system.  _Journal of Computational Chemistry_ 14(11): 1347–1363.  10.1002/jcc.540141112.  URL https://doi.org/10.1002%2Fjcc.540141112.
  * Schwegler and Challacombe (1996) Schwegler E and Challacombe M (1996) Linear scaling computation of the hartree–fock exchange matrix.  _The Journal of Chemical Physics_ 105(7): 2726–2734.  10.1063/1.472135.  URL https://doi.org/10.1063%2F1.472135.
  * Schwegler et al. (1997) Schwegler E, Challacombe M and Head-Gordon M (1997) Linear scaling computation of the fock matrix. II. rigorous bounds on exchange integrals and incremental fock build.  _The Journal of Chemical Physics_ 106(23): 9708–9717.  10.1063/1.473833.  URL https://doi.org/10.1063%2F1.473833.
  * Shang et al. (2011) Shang H, Li Z and Yang J (2011) Implementation of screened hybrid density functional for periodic systems with numerical atomic orbitals: Basis function fitting and integral screening.  _The Journal of Chemical Physics_ 135(3): 034110.  10.1063/1.3610379.  URL https://doi.org/10.1063%2F1.3610379.
  * Sodt and Head-Gordon (2008) Sodt A and Head-Gordon M (2008) Hartree-Fock exchange computed using the atomic resolution of the identity approximation.  _The Journal of Chemical Physics_ 128(10): 104106.  10.1063/1.2828533.  URL https://doi.org/10.1063%2F1.2828533.
  * Soler et al. (2002) Soler JM, Artacho E, Gale JD, García A, Junquera J, Ordejón P and Sánchez-Portal D (2002) The siesta method for ab initio order-$n$ materials simulation.  _Journal of Physics: Condensed Matter_ 14(11): 2745–2779.  10.1088/0953-8984/14/11/302.  URL https://doi.org/10.1088%2F0953-8984%2F14%2F11%2F302.
  * Stephens et al. (1994) Stephens PJ, Devlin FJ, Chabalowski CF and Frisch MJ (1994) Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields.  _The Journal of Physical Chemistry_ 98(45): 11623–11627.  10.1021/j100096a001.  URL https://doi.org/10.1021%2Fj100096a001.
  * Troullier and Martins (1991) Troullier N and Martins JL (1991) Efficient pseudopotentials for plane-wave calculations.  _Physical Review B_ 43(3): 1993–2006.  10.1103/physrevb.43.1993.  URL https://doi.org/10.1103%2Fphysrevb.43.1993.
  * Tymczak and Challacombe (2005) Tymczak CJ and Challacombe M (2005) Linear scaling computation of the fock matrix. VII. periodic density functional theory at the $\gamma$ point.  _The Journal of Chemical Physics_ 122(13): 134102.  10.1063/1.1853374.  URL https://doi.org/10.1063%2F1.1853374.
  * Valiev et al. (2010) Valiev M, Bylaska E, Govind N, Kowalski K, Straatsma T, Dam HV, Wang D, Nieplocha J, Apra E, Windus T and de Jong W (2010) NWChem: A comprehensive and scalable open-source solution for large scale molecular simulations.  _Computer Physics Communications_ 181(9): 1477–1489.  10.1016/j.cpc.2010.04.018.  URL https://doi.org/10.1016%2Fj.cpc.2010.04.018.
  * VandeVondele et al. (2005) VandeVondele J, Krack M, Mohamed F, Parrinello M, Chassaing T and Hutter J (2005) Quickstep: Fast and accurate density functional calculations using a mixed gaussian and plane waves approach.  _Computer Physics Communications_ 167(2): 103–128.  10.1016/j.cpc.2004.12.014.  URL https://doi.org/10.1016%2Fj.cpc.2004.12.014.

##  7 Author Biographies

Xinming Qin is a PhD student in chemical physics at USTC under the supervision
of Prof.Jinlong Yang. He received the BS degree in chemistry in July 2009 at
the same university. His main research interests are developing and applying
new algorithms for large-scale electronic structure calculations. He is the
main developer of the HONPAS.

Honghui Shang is an associate professor at the Computer Science at Institute
of Computing Technology, Chinese Academy of Sciences, China. She received her
BS degree in physics from University of Science and Technology of China in
2006, and the PhD degree in Physical Chemistry from University of Science and
Technology of China, in 2011. Between 2012 and 2018, she worked as a
postdoctoral research assistant at the Fritz Haber Institute of the Max Planck
Society, Berlin, Germany, which was responsible for the Hybrid
Inorganic/Organic Systems for Opto-Electronics (HIOS) project and for the
Novel Materials Discovery (NOMAD) project. Her main research interests are
developing the physical algorithm or numerical methods for the first-principle
calculations as well as accelerating these applications in the high
performance computers. Currently, she is the main developer of the HONPAS
(leader of of the hybrid-density-functional part) and FHI-aims (leader of the
density-functional perturbation theory part).

Lei Xu is currently working toward the BS degree in the department of physics
at the SiChuan University. His current research interests are the parallel
programming in the high performance computing domain.

Wei Hu is currently a professor at division of theoretical and computational
sciences at Hefei National Laboratory for Physical Sciences at Microscale
(HFNL) at USTC. He received the BS degree in chemistry from USTC in 2007, and
the PhD degree in Physical Chemistry from the same university in 2013. From
2014 to 2018, he worked as a postdoctoral fellow in the Scalable Solvers Group
of the Computational Research Division at Lawrence Berkeley National
Laboratory (LBNL), Berkeley, USA. During his postdoctoral research, he
developed a new massively parallel methodology, DGDFT (Discontinuous Galerkin
Method for Density Functional Theory), for large-scale DFT calculations. His
main research interests focus on development, parallelization, and application
of advanced and highly efficient DFT methods and codes for accurate first-
principles modeling and simulations of nanomaterials.

Jinlong Yang is a full professor of chemistry and executive dean of the School
of Chemistry and Material Sciences at USTC. He obtained his PhD degree in 1991
from USTC. He was awarded Outstanding Youth Foundation of China in 2000,
selected as Changjiang Scholars Program Chair Professor in 2001 and as a
fellow of American Physical Society (APS) in 2011. He is the second awardee of
the 2005’s National Award for Natural Science (the second prize) and the
awardee (principal contributor) of the 2014’s Outstanding Science and
Technology Achievement Prize of the Chinese Academy of Sciences (CAS). His
research mainly focuses on the development of first-principles methods and
their application to clusters, nanostructures, solid materials, surfaces, and
interfaces. He is the initiator and leader of the HONPAS.

Shigang Li received his Bachelor in computer science and technology and PhD in
computer architecture from the University of Science and Technology Beijing,
China, in 2009 and 2014, respectively. He was funded by CSC for a 2-year PhD
study in University of Illinois, Urbana-Champaign. He was an assistant
professor (from June 2014 to Aug. 2018) in State Key Lab of Computer
Architecture, Institute of Computing Technology, Chinese Academy of Sciences
at the time of his achievement in this work. From Aug. 2018 to now, he is a
postdoc researcher in Department of Computer Science, ETH Zurich. His research
interests focus on the performance optimization for parallel and distributed
computing systems, including parallel algorithms, parallel programming model,
performance model, and intelligent methods for performance optimization.

Yunquan Zhang received his BS degree in computer science and engineering from
the Beijing Institute of Technology in 1995. He received a PhD degree in
Computer Software and Theory from the Chinese Academy of Sciences in 2000. He
is a full professor and PhD Advisor of State Key Lab of Computer Architecture,
ICT, CAS. He is also appointed as the Director of National Supercomputing
Center at Jinan and the General Secretary of China’s High Performance
Computing Expert Committee. He organizes and distributes China’s TOP100 List
of High Performance Computers, which traces and reports the development of the
HPC system technology and usage in China. His research interests include the
areas of high performance parallel computing, focusing on parallel programming
models, high-performance numerical algorithms, and performance modeling and
evaluation for parallel programs.