# Exact exchange matrix of periodic Hartree-Fock theory for all-electron
simulations

Qiming Sun<EMAIL_ADDRESS>AxiomQuant Investment Management LLC,
Shanghai 200120, China

###### Abstract

This work presents an algorithm to evaluate Coulomb and exchange matrices in
Fock operator using range separation techniques at various aspects. This
algorithm is particularly favorable for the scenario of (1) all-electron
calculations or (2) computing exchange matrix for a large number of
$\mathbf{k}$-point samples. An all electron Hartree-Fock calculation with 110k
basis functions is demonstrated in this work.

##  I Introduction

In crystal electronic structure simulations, planewave basis is a natural
choice because of the periodicity in wave functionKresse and Joubert (1999);
Enkovaara  _et al._ (2010); Royo and Stengel (2019); Giannozzi  _et al._
(2020). When the atomic nature of electron density is taken into account,
numerical basis and periodic Gaussian basis are also widely usedDovesi  _et
al._ (2018); Kühne  _et al._ (2020); Sun  _et al._ (2020). However, we can
still find the footprint of planewaves in the intermediates. If the quantity
of interests, e.g. the electron density, changes smoothly in the space, they
can be efficiently evaluated with the assistance of planewaves and discrete
Fourier transform (FT). As such, crystal simulations with density functional
theory (DFT) were rapidly evolved. Large scale DFT calculations have reached
$\sim$100k atoms and millions of basis functionsBowler and Miyazaki (2010);
VandeVondele, Borštnik, and Hutter (2012).

Along with the choice of planewave basis, electronic structure simulations for
crystal face two challenges: the treatment of exact exchange and the quality
of pseudo-potential.

Evaluating Hartree-Fock (HF) exchange in crystal simulations is not trival.
Treatments of divergence and finite size errors have been addressed by many
pioneer worksPisani and Dovesi (1980); Dovesi  _et al._ (1983); Gygi and
Baldereschi (1986); Heyd and Scuseria (2004); Izmaylov, Scuseria, and Frisch
(2006); Spencer and Alavi (2008); Guidon, Hutter, and VandeVondele (2009);
Sundararaman and Arias (2013). Besides, various efforts have been made to
improve the efficiency of building HF exchange in the context of planewave
basis and the associated discrete FTMarzari and Vanderbilt (1997); Gygi
(2009); Marzari  _et al._ (2012); Gygi and Duchemin (2013); DiStasio  _et al._
(2014); Hu  _et al._ (2017); Hu, Lin, and Yang (2017). Regardless of the basis
type or the algorithm or the approximations, the computational cost for HF
exchange in terms of discrete FT scales as
$\mathcal{O}(N_{A}^{2}N_{k}^{2}N_{G}\log(N_{G}))$ for a simulation of $N_{A}$
atoms, $N_{k}$ $\mathbf{k}$-points and $N_{G}$ planewaves (or integration mesh
grids in real space). It becomes unaffordable for systems which requires a
relatively large number of planewaves. In the context of Gaussian type crystal
basis, one would use auxiliary Gaussian basis to expand electron density and
invoke density fitting (DF) technique to evaluate the Coulomb repulsion
integralsTowler, Zupan, and Causà (1996); Schwegler and Challacombe (1999);
Heyd, Scuseria, and Ernzerhof (2003); Tymczak  _et al._ (2005); Izmaylov,
Scuseria, and Frisch (2006); Maschio  _et al._ (2007); Varga (2008); Maschio
and Usvyat (2008); Guidon, Hutter, and VandeVondele (2009); Ben, Hutter, and
VandeVondele (2013); Yamada, Brewster, and Hirata (2013); Sun  _et al._
(2017); Irmler, Burow, and Pauly (2018); Wang, Lewis, and Valeev (2020). It
largely removes the demands to the planewaves (or integration grids) and
reduces the formal scaling to $\mathcal{O}(N_{AO}^{4}N_{k}^{2})$ in
computational costs and $\mathcal{O}(N_{AO}^{3}N_{k}^{2})$ in memory (or disk)
footprint for $N_{AO}$ Gaussian basis functions. Although DF can deliver great
performance for simulations with moderate system size and a moderate number of
$\mathbf{k}$ points, occasionally, memory and disk usage can become a
bottleneck when hundreds of $\mathbf{k}$ points are requested. Another issue
of DF is that it is picky on basis sets and auxiliary basis functions. Most
auxiliary basis sets were developed and tested for specific atomic basis sets
in the context of molecules. Their performance are yet extensively validated
for systems with periodic boundaries.

The contradiction of efficiency and accuracy is rooted in the development of
pseudo potentialHartwigsen, Goedecker, and Hutter (1998); Willand  _et al._
(2013); Gao and Chelikowsky (2018); Prandini  _et al._ (2018); Lu  _et al._
(2019). Because of the limited choices of available soft pseudo potentials,
one often would have to sacrifice the reliability in a real application, not
to mention the limitations of pseudo potential if one would like to study the
problems beyond closed-shell valence properties, such as relativistic effects,
core-valence correlations, core electron spectroscopy, spin polarization and
so forth. All-electron simulations would be the key of these challenges.

Developing an algorithm suitable for arbitrary basis sets and arbitrary core
treatment is the goal of the current work we would like to achieve. As shown
in Irmler’s workIrmler, Burow, and Pauly (2018), it is efficient to evaluate
exchange matrix in the two-passes algorithm: First compute the exchange matrix
based on Born-von Karman supercell in real space then transform to the
representations adapted to $\mathbf{k}$-point symmetry in reciprocal space. In
their work, they demonstrated several HF calculations with thousands of
$\mathbf{k}$ points. Nonetheless, the accuracy of Irmler’s work is affected by
the artifacts of truncated Coulomb kernel near the truncation boundary. We
adopt the similar 2-passes strategy to evaluate Coulomb and exchange matrices
in our algorithm. In contrast to their treatment of Coulomb metric, we handle
the full Coulomb metric without truncation.

The formulaism of the algorithm are established in Section II. The algorithm
in the present work requires various integral screen thresholds and cutoffs.
They are discussed in Section II.4. Results and performance benchmark are
presented in Section III.

##  II Theory

###  II.1 Fock matrix computed in the range-separation of Coulomb metric

We consider here a periodic system which consists of repeated unit cells with
translation vectors $\mathbf{m}$ and lattice vectors $\mathbf{a}$

$\mathbf{m}=m_{x}\mathbf{a}_{x}+m_{y}\mathbf{a}_{y}+m_{z}\mathbf{a}_{z}.$ (1)

Crystal orbitals for a crystal of $N\rightarrow\infty$ unit cells can be
expanded with a set of Bloch orbitals

$\phi_{\mu}(\mathbf{r})^{\mathbf{k}}=\frac{1}{\sqrt{N}}\sum_{\mathbf{m}}e^{i\mathbf{k}\cdot\mathbf{m}}\mu^{\mathbf{m}}(\mathbf{r})$
(2)

where $\mu^{\mathbf{m}}$ is an real Gaussian-type orbital centered at
$\mathbf{R}_{\mu}$ in the unit cell regarding to the translation vector
$\mathbf{m}$

$\mu^{\mathbf{m}}(\mathbf{r})=\mu(\mathbf{r}-\mathbf{R}_{\mu}-\mathbf{m})$ (3)

Solving the periodic HF model can lead us to the $\mathbf{k}$-point symmetry
adapted Roothaan equations

$\displaystyle\mathbf{F}^{\mathbf{k}}\mathbf{C}^{\mathbf{k}}=\mathbf{S}^{\mathbf{k}}\mathbf{C}^{\mathbf{k}}\epsilon^{\mathbf{k}},$
(4)

and the corresponding $\mathbf{k}$-point symmetry adapted Fock matrix and
overlap matrix

$\displaystyle\mathbf{F}^{\mathbf{k}}=\mathbf{T}^{\mathbf{k}}+\mathbf{V}_{\mathrm{Nuc}}^{\mathbf{k}}+\mathbf{J}^{\mathbf{k}}-\mathbf{K}^{\mathbf{k}},$
(5) $\displaystyle
S_{\mu\nu}^{\mathbf{k}}=\langle\phi_{\mu}^{\mathbf{k}}|\phi_{\nu}^{\mathbf{k}}\rangle$
(6)

The overlap matrix can be computed using discrete FT of the overlap integrals
in real space

$\displaystyle
S_{\mu\nu}^{\mathbf{k}}=\sum_{\mathbf{m}}e^{i\mathbf{k}\cdot\mathbf{m}}S_{\mu^{0}\nu^{m}},$
(7) $\displaystyle
S_{\mu^{0}\nu^{m}}=\langle\mu^{\mathbf{0}}|\nu^{\mathbf{m}}\rangle.$ (8)

In a similarly way, we can get the kinetic matrix in $\mathbf{F}^{\mathbf{k}}$

$\displaystyle
T_{\mu\nu}^{\mathbf{k}}=\sum_{\mathbf{m}}e^{i\mathbf{k}\cdot\mathbf{m}}T_{\mu^{0}\nu^{m}}.$
(9)

When computing the nuclear attraction matrix
$\mathbf{V}_{\mathrm{Nuc}}^{\mathbf{k}}$, we use model charge to screen the
long-range interactions then compute separately the long-range and short-range
parts of nuclear attraction integrals. Relevant techniques have been
documented in our earlier work Sun  _et al._ , 2017.

The HF Coulomb and exchange matrices in $\mathbf{F}^{\mathbf{k}}$ are

$\displaystyle
J_{\mu\nu}^{\mathbf{k}}=\sum_{\mathbf{k}^{\prime}}\sum_{\kappa\lambda}D_{\lambda\kappa}^{\mathbf{k}^{\prime}}(\phi_{\mu}^{\mathbf{k}}\phi_{\nu}^{\mathbf{k}}|\phi_{\kappa}^{\mathbf{k}^{\prime}}\phi_{\lambda}^{\mathbf{k}^{\prime}}),$
(10) $\displaystyle
K_{\mu\lambda}^{\mathbf{k}}=\sum_{\mathbf{k}^{\prime}}\sum_{\kappa\nu}D_{\nu\kappa}^{\mathbf{k}^{\prime}}(\phi_{\mu}^{\mathbf{k}}\phi_{\nu}^{\mathbf{k}^{\prime}}|\phi_{\kappa}^{\mathbf{k}^{\prime}}\phi_{\lambda}^{\mathbf{k}}).$
(11)

The four-center two-electron repulsion integrals in equation (11) follow
Mulliken notation

$\displaystyle(\phi_{\mu}^{\mathbf{k}_{\mu}}\phi_{\nu}^{\mathbf{k}_{\nu}}|\phi_{\kappa}^{\mathbf{k}_{\kappa}}\phi_{\lambda}^{\mathbf{k}_{\lambda}})=\iint\phi_{\mu}^{\mathbf{k}_{\mu}}(\mathbf{r}_{1})\phi_{\nu}^{\mathbf{k}_{\nu}}(\mathbf{r}_{1})g(\mathbf{r}_{12})\phi_{\kappa}^{\mathbf{k}_{\kappa}}(\mathbf{r}_{2})\phi_{\lambda}^{\mathbf{k}_{\lambda}}(\mathbf{r}_{2})d^{3}\mathbf{r}_{1}d^{3}\mathbf{r}_{2}.$
(12)

By splitting the Coulomb metric $g(\mathbf{r}_{12})$ into long-range (LR) plus
short-range (SR) contributions

$g(\mathbf{r}_{12})=\frac{1}{r_{12}}=\frac{\mathrm{erf}(\omega
r_{12})}{r_{12}}+\frac{\mathrm{erfc}(\omega r_{12})}{r_{12}},$ (13)

we can evaluate the electron repulsion integrals of the two components
separately in reciprocal space and real space.

We will focus on the computation of exchange matrix in the following context.
The Coulomb matrix can be computed in a similar fashion. The LR component of
the ERI (12) can be evaluated in reciprocal space, yielding a summation over
planewaves
$\mathbf{G}=\frac{(2\pi)^{3/2}}{\sqrt{N}}e^{i\mathbf{G}\cdot\mathbf{r}}$
inside the volume $\Omega$ of the unit cellDominik Marx (2009)

$\displaystyle(\phi_{\mu}^{\mathbf{k}}\phi_{\nu}^{\mathbf{k}^{\prime}}|\phi_{\kappa}^{\mathbf{k}^{\prime}}\phi_{\lambda}^{\mathbf{k}})^{\mathrm{LR}}=\frac{1}{\Omega}\sum_{\mathbf{G}}\frac{4\pi\exp(-\frac{|\mathbf{G}+\mathbf{k}^{\prime}-\mathbf{k}|^{2}}{4\omega})}{|\mathbf{G}+\mathbf{k}^{\prime}-\mathbf{k}|^{2}}\rho_{\mu\nu}^{\mathbf{k}\mathbf{k}^{\prime}}(\mathbf{G})\rho_{\kappa\lambda}^{\mathbf{k}\mathbf{k}^{\prime}}(-\mathbf{G}),$
(14) $\displaystyle\begin{aligned}
\rho_{\mu\nu}^{\mathbf{k}\mathbf{k}^{\prime}}(\mathbf{G})&=\int
e^{-i(\mathbf{G}+\mathbf{k}^{\prime}-\mathbf{k})\cdot\mathbf{r}}\phi_{\mu}^{\mathbf{k}*}(\mathbf{r})\phi_{\nu}^{\mathbf{k}^{\prime}}(\mathbf{r})d^{3}\mathbf{G}\\\
&=\sum_{\mathbf{n}}e^{i\mathbf{k}^{\prime}\cdot\mathbf{n}}\int
e^{-i(\mathbf{G}+\mathbf{k}^{\prime}-\mathbf{k})\cdot\mathbf{r}}\mu(\mathbf{r})\nu(\mathbf{r}-\mathbf{n})d^{3}\mathbf{G}.\end{aligned}$
(15)

In Eq. (14), a limited number of planewaves is required to converge the
integrals thanks to the exponential decay quantity
$\exp(-\frac{G^{2}}{4\omega})$. We then ready to obtain the LR contributions
of the exchange matrix using Eq. (10) and Eq. (11).

SR Coulomb integrals can be evaluated in real space. Following the schemes
developed in Irmler, Burow, and Pauly, 2018, we first perform the SR part of
the exchange matrix in the Born-von Karman (BvK) supercell

$\displaystyle
K_{\mu^{0}\lambda^{\mathbf{l}}}^{\mathrm{SR}}=\sum_{\nu\kappa}\sum_{\mathbf{L}}\sum_{\mathbf{m}+\mathbf{M}}\sum_{\mathbf{n}+\mathbf{N}}D_{\nu^{\mathbf{n}}\kappa^{\mathbf{m}}}(\mu^{\mathbf{0}}\nu^{\mathbf{n}+\mathbf{N}}|\kappa^{\mathbf{m}+\mathbf{M}}\lambda^{\mathbf{l}+\mathbf{L}+\mathbf{M}})^{\mathrm{SR}}.$
(16)

Here involves two types of translation vectors, the super-lattice translation
vectors for BvK supercell ($\mathbf{L},\mathbf{M},\mathbf{N}$) and the
translation inside a single BvK supercell
($\mathbf{l},\mathbf{m},\mathbf{n}$). SR Coulomb integrals as well as the
exchange matrix can be converged with a finite number of super lattice
translation vectors. Similar to the treatment for overlap matrix in Eq. (7), a
discrete FT can transform Eq. (16) to the $\mathbf{k}$-point symmetry adapted
exchange matrix. The quantity $D_{\nu^{\mathbf{n}}\kappa^{\mathbf{m}}}$ in Eq.
(16) is density matrix of BvK supercell. The translation symmetry ensures that
the density matrix is invariant against the super-lattice translations

$D_{\nu^{\mathbf{n}}\kappa^{\mathbf{m}}}=D_{\nu^{\mathbf{n}+\mathbf{N}}\kappa^{\mathbf{m}+\mathbf{M}}}.$
(17)

For the BvK supercell being filled with $N_{\mathrm{BvK}}$ primitive cells,
$D_{\nu^{\mathbf{l}}\kappa^{\mathbf{m}}}$ is obtained from back-FT of the
$\mathbf{k}$-adapted density matrices

$D_{\nu^{\mathbf{n}}\kappa^{\mathbf{m}}}=D_{\nu^{\mathbf{0}}\kappa^{\mathbf{m}-\mathbf{n}}}=\frac{1}{N_{\mathrm{BvK}}}\sum_{\mathbf{k}}e^{-i\mathbf{k}\cdot(\mathbf{m}-\mathbf{n})}D_{\nu\kappa}^{\mathbf{k}}.$
(18)

When comparing the so-obtained exchange matrix to the one based on discrete
FT, we noticed that they are different at $\mathbf{G}=0$, whose contribution
is

$K_{\mu\lambda}^{\mathbf{k}}\Big{|}_{\mathbf{G}=0}=\frac{\pi}{\Omega\omega^{2}}\sum_{\nu\kappa}S_{\mu\nu}^{\mathbf{k}}D_{\nu\kappa}^{\mathbf{k}}S_{\kappa\lambda}^{\mathbf{k}}.$
(19)

###  II.2 SR integrals in real space

When a system employs diffused basis functions, the integration scheme above
may still require a huge number repeated images of BvK supercell in the
crystal to converge the SR integrals. It is efficient to compute the integrals
in reciprocal space if the relevant basis functions are all smooth functions.
We thus divided the entire basis sets into two groups based on their locality,
using notation $D$ and $C$ to label the diffused functions and the compact
functions. The products of two basis functions $\rho_{\mu\nu}$ have four
possible combinations: $(\rho_{D})_{\mu\nu}\sim DD$ is the product of two
diffused functions; $DC$ and $CD$ are the product of one diffused function and
one compact function; $CC$ is the product of two compact functions

$\displaystyle\\{\rho_{\mu\nu}\\}=\\{(\rho_{D})_{\mu\nu}\\}\oplus\\{(\overline{\rho_{D}})_{\mu\nu}\\},$
$\displaystyle(\rho_{D})_{\mu\nu}\sim DD,$
$\displaystyle{(\overline{\rho_{D}})_{\mu\nu}}\sim CD,DC,CC.$

When computing the two-electron SR integrals, the products
$\overline{\rho_{D}}\cdot\overline{\rho_{D}}$ (in total 9 components) are
evaluated analytically in real space. All other integrals are evaluated in
reciprocal space, including: the long-range part of
$\overline{\rho_{D}}\cdot\overline{\rho_{D}}$ and the full Coulomb integrals
of $\rho_{D}\cdot\overline{\rho_{D}}$, $\overline{\rho_{D}}\cdot\rho_{D}$, and
$\rho_{D}\cdot\rho_{D}$.

###  II.3 Integrals in reciprocal space

Integrals in the reciprocal space are computed

$(\phi_{\mu}^{\mathbf{k}}\phi_{\nu}^{\mathbf{k}^{\prime}}|\phi_{\kappa}^{\mathbf{k}^{\prime}}\phi_{\lambda}^{\mathbf{k}})=\frac{1}{\Omega}\sum_{\mathbf{G}}g(\mathbf{G},\mathbf{k},\mathbf{k}^{\prime})\rho_{\mu\nu}^{\mathbf{k}\mathbf{k}^{\prime}}(\mathbf{G})\rho_{\kappa\lambda}^{\mathbf{k}\mathbf{k}^{\prime}}(-\mathbf{G}),$
(20)

where the Coulomb kernel $g(\mathbf{G},\mathbf{k},\mathbf{k}^{\prime})$ can be

$g(\mathbf{G},\mathbf{k},\mathbf{k}^{\prime})=\frac{4\pi}{|\mathbf{G}+\mathbf{k}^{\prime}-\mathbf{k}|^{2}}$
(21)

for the regular full Coulomb interactions and

$g_{LR}(\mathbf{G},\mathbf{k},\mathbf{k}^{\prime})=\frac{4\pi\exp(-\frac{|\mathbf{G}+\mathbf{k}^{\prime}-\mathbf{k}|^{2}}{4\omega})}{|\mathbf{G}+\mathbf{k}^{\prime}-\mathbf{k}|^{2}}$
(22)

for the long-range Coulomb attenuation. By using the notations $D$ and $C$ for
diffused and compact basis functions, we get 7 types of two-electron integrals
computed with the full Coulomb kernel (21). They are $(DD|DD)$, $(DD|CD)$,
$(DD|DC)$, $(DD|CC)$, $(CD|DD)$, $(DC|DD)$, $(CC|DD)$. There are 9 types of
integrals $(CC|CC)$, $(CC|CD)$, $(CC|DC)$, $(CD|CD)$, $(CD|DC)$, $(CD|CC)$,
$(DC|CD)$, $(DC|DC)$, $(DC|CC)$ associated with the LR Coulomb attenuation
kernel (22). For the 9 types of LR integrals, the contributions of
$\mathbf{G}=0$ as shown in (19) which was brought by the real-space analytical
integration need to be subtracted.

###  II.4 Truncation criteria

####  II.4.1 Coulomb attenuation parameter

Our first parameter to determine is the range-separation parameter $\omega$ in
the attenuated Coulomb metric. The optimal choice would be the one that can
balance the cost of reciprocal space and real space computation. It can be
influenced by many factors, e.g. the lattice parameters, the ratio of diffused
functions in the entire basis sets, the desired numerical precision, the
performance of underlying integral engine and FT engine. In our current
implementation, the reciprocal space computation is dominated by the
analytical FT (15), in which we need to invoke $N_{k}$ times the analytical FT
for $N_{G}$ planewaves. In each call of the FT, the complexity is
$\mathcal{O}(N_{G}N_{A}^{2}N_{I})$ regarding to the $N_{AO}$ basis functions
and $N_{I}$ the number translation vectors in lattice sum. We thus target
that, with certain value of $\omega$ and its associated $N_{G}$, the
complexity of the dominant step in the real space computation
$\mathcal{O}(N_{A}^{4}N_{I}^{3})$ is comparable to the FT complexity. Under
this guideline, we adjust certain factors and utilize the estimation
expression for $\omega$

$\omega\approx
0.01\Big{(}\frac{N_{AO}N_{I}\Omega}{N_{k}}\Big{)}^{\frac{1}{3}}.$ (23)

####  II.4.2 Integral prescreening

We employ two types of Schwarz inequalities

$\displaystyle(\mu\nu|\kappa\lambda)<\sqrt{(\mu\nu|\mu\nu)(\kappa\lambda|\kappa\lambda)},$
(24)
$\displaystyle(\mu\nu|\kappa\lambda)<\min(\sqrt{(\mu\mu|\kappa\kappa)(\nu\nu|\lambda\lambda)},\sqrt{(\mu\mu|\lambda\lambda)(\nu\nu|\kappa\kappa)}).$
(25)

The second condition is not quite useful for regular Coulomb metric due to the
LR character of Coulomb interactions. Neither
$\sqrt{(\mu\mu|\kappa\kappa)(\nu\nu|\lambda\lambda)}$ nor
$\sqrt{(\mu\mu|\lambda\lambda)(\nu\nu|\kappa\kappa)}$ will vanish regardless
how the underlying basis functions are separated. In contrast to regular
Coulomb, this condition filters about 50% terms for SR integrals. These cases
are corresponding to the interactions between the remotely separated bra and
ket functions which are quite common to the lattice sum integration.

####  II.4.3 Energy cutoff

Given Coulomb attenuation parameter $\omega$ and the desired numerical
precision $\varepsilon$, the energy cutoff
$E_{cut}=\frac{1}{2}|\mathbf{G}_{max}|^{2}$ to converge LR integrals in
reciprocal space is ready to obtain

$\frac{1}{\Omega}\sum_{\mathbf{G}=0}^{\infty}\rho(\mathbf{G})\frac{4\pi\exp(-\frac{G^{2}}{4\omega^{2}})}{G^{2}}\rho(-\mathbf{G})=\frac{1}{\Omega}\sum_{\mathbf{G}=0}^{E_{cut}}\rho(\mathbf{G})\frac{4\pi\exp(-\frac{G^{2}}{4\omega^{2}})}{G^{2}}\rho(-\mathbf{G})+\delta.$
(26)

Our estimation for the error $\delta$ is

$\displaystyle\delta$
$\displaystyle=\frac{1}{\Omega}\sum_{|\mathbf{G}|>E_{cut}}\rho(\mathbf{G})\frac{4\pi\exp(-\frac{G^{2}}{4\omega^{2}})}{G^{2}}\rho(-\mathbf{G})$
(27)
$\displaystyle\approx\int_{G_{max}-\frac{1}{2}G_{min}}^{\infty}\rho(\mathbf{G})\frac{4\pi\exp(-\frac{G^{2}}{4\omega})}{G^{2}}\rho(-\mathbf{G})d^{3}\mathbf{G}.$

$G_{min}$ is the smallest non-zero $\mathbf{G}$. The Fourier transformed
Gaussian product $\rho(\mathbf{G})$ is a Gaussian distribution.
$\rho(\mathbf{G})$ is always smaller than 1. The error estimation can be
reduction to

$\displaystyle\delta$
$\displaystyle<\int_{G_{max}-\frac{1}{2}G_{min}}^{\infty}\frac{4\pi\exp(-\frac{G^{2}}{4\omega^{2}})}{G^{2}}d^{3}\mathbf{G}$
$\displaystyle=\int_{G_{max}-\frac{1}{2}G_{min}}^{\infty}16\pi^{2}\exp(-\frac{G^{2}}{4\omega^{2}})dG$
$\displaystyle<\frac{1}{G_{max}}\int_{G_{max}-\frac{1}{2}G_{min}}^{\infty}16\pi^{2}G\exp(-\frac{G^{2}}{4\omega^{2}})dG$
$\displaystyle<\frac{32\pi^{2}\omega^{2}\exp(-\frac{(G_{max}-\frac{1}{2}G_{min})^{2}}{4\omega^{2}})}{G_{max}}$
(28)

We can set the energy cutoff

$E_{cut}=\frac{1}{2}|\mathbf{G}_{max}|^{2}>2\omega^{2}\log(\frac{32\pi^{2}\omega^{2}}{\varepsilon})$
(29)

which ensures the integral error $\delta$ of LR integrals smaller than the
required precision $\varepsilon$ for all basis functions.

####  II.4.4 Basis locality

We are now targeting a criteria to classify the diffused functions and the
compact basis functions. Based on this partitioning criteria $\rho_{D}$ the
product of two diffused functions needs to be smooth enough to converge the
reciprocal integration (20) under the full Coulomb kernel.

This error is estimated

$\delta=\frac{1}{\Omega}\sum_{\mathbf{G}\sim\sqrt{2E_{cut}}}^{\infty}\rho_{D}(\mathbf{G})\frac{4\pi}{G^{2}}\rho(-\mathbf{G})\approx\int_{\sqrt{2E_{cut}}}^{\infty}\rho_{D}(\mathbf{G})\frac{4\pi}{G^{2}}\rho(-\mathbf{G})d^{3}\mathbf{G}.$
(30)

$\rho(-\mathbf{G})$ is smaller than 1 in any circumstance.

$\delta<\int_{\sqrt{2E_{cut}}}^{\infty}\rho_{D}(\mathbf{G})\frac{4\pi}{G^{2}}d^{3}\mathbf{G}.$
(31)

If $\rho_{D}$ is formed from the same primitive $s$-type Gaussians at the same
center, we have

$\displaystyle(\rho_{D})_{\mu\nu}(\mathbf{r})=\Big{(}\frac{2\alpha}{\pi}\Big{)}^{\frac{3}{2}}\exp({-\alpha(\mathbf{r}-\mathbf{R}_{\mu})^{2}})\exp({-\alpha(\mathbf{r}-\mathbf{R}_{\nu})^{2}}).$
(32)

The product under FT appears

$\displaystyle(\rho_{D})_{\mu\nu}(\mathbf{G})=\exp(-\frac{G^{2}}{8\alpha})\exp(-\frac{i}{2}\mathbf{G}\cdot(\mathbf{R}_{\mu}+\mathbf{R}_{\nu}))$
$\displaystyle|(\rho_{D})_{\mu\nu}(\mathbf{G})|<\exp(-\frac{G^{2}}{8\alpha})$

We then derive the integration error

$\delta<\int_{\sqrt{2E_{cut}}}^{\infty}16\pi^{2}\exp(-\frac{G^{2}}{8\alpha})dG<\frac{64\pi^{2}\alpha\exp(-\frac{E_{cut}}{4\alpha})}{\sqrt{2E_{cut}}}$
(33)

which yields the boundary of the Gaussian exponent

$\displaystyle\alpha_{cut}<E_{cut}\Big{/}\Big{(}4\log(\frac{64\pi^{2}}{\sqrt{2E_{cut}}\varepsilon})\Big{)}.$
(34)

For basis functions of angular momentum $l$, we analyze the leading term of
the Fourier transformed orbital products in the similar manner and obtain an
estimation of the boundary

$\displaystyle\alpha_{cut}<E_{cut}\Big{/}\Big{(}4\log(\frac{64\pi^{2}E_{cut}^{l-\frac{1}{2}}}{2^{l+\frac{1}{2}}(2l+1)!!\varepsilon})\Big{)}.$

We use the threshold $\alpha_{cut}$ to classify the characters of basis
functions. The basis functions in the smooth set are the primitive Gaussians
with exponents smaller than $\alpha_{cut}$. The rest basis functions are
grouped in the compact set. For contracted basis functions, we first de-
contract the bases then carry out the partitioning for the primitive
functions.

####  II.4.5 Lattice sum boundary

Lattice sum in real space is required in two places. First is the FT of the
orbital product (15). Second is the SR integrals (16) computed in the real
space.

For simplicity, we ignore the $\mathbf{k}$ points in Eq. (15)

$\rho_{\mu\nu}(\mathbf{G})=\sum_{\mathbf{N}}\int
e^{-i\mathbf{G}\cdot\mathbf{r}}\mu(\mathbf{r})\nu(\mathbf{r}-\mathbf{N})d^{3}\mathbf{r}.$
(35)

The real space truncation error is

$\displaystyle\delta=\sum_{|\mathbf{N}|>R_{cut}}\int
e^{-i\mathbf{G}\cdot\mathbf{r}}\mu(\mathbf{r})\nu(\mathbf{r}-\mathbf{N})d^{3}\mathbf{r}<\sum_{|\mathbf{N}|>R_{cut}}S_{\mu\nu}^{~{}\mathbf{N}},$
(36) $\displaystyle
S_{\mu\nu}^{~{}\mathbf{N}}=\int\mu(\mathbf{r})\nu(\mathbf{r}-\mathbf{N})d^{3}\mathbf{r}.$
(37)

The overlap $S_{\mu\nu}^{~{}\mathbf{N}}$ in terms of $s$-type Gaussians is

$S_{\mu\nu}^{~{}\mathbf{N}}=\Big{(}\frac{4\alpha_{\mu}\alpha_{\nu}}{\alpha_{\mu\nu}^{2}}\Big{)}^{3/4}e^{-\eta_{\mu\nu}^{~{}\mathbf{N}}}$
(38)

where

$\displaystyle\alpha_{\mu\nu}=\alpha_{\mu}+\alpha_{\nu},$ (39)
$\displaystyle\eta_{\mu\nu}^{~{}\mathbf{N}}=\frac{\alpha_{\mu}\alpha_{\nu}}{\alpha_{\mu\nu}}|\mathbf{R}_{\mu}-\mathbf{R}_{\nu}-\mathbf{N}|^{2}.$
(40)

The error $\delta$ would be maximized when $\mu$ and $\nu$ are located at the
same atomic center. The upper bound of the error for large distance
$|\mathbf{N}|$ is

$\displaystyle\frac{2\pi\alpha_{\mu\nu}|\mathbf{N}|}{\alpha_{\mu}\alpha_{\nu}}\Big{(}\frac{4\alpha_{\mu}\alpha_{\nu}}{\alpha_{\mu\nu}^{2}}\Big{)}^{3/4}\exp\Big{(}-\frac{\alpha_{\mu}\alpha_{\nu}(|\mathbf{N}|-\frac{\sqrt[3]{\Omega}}{2})^{2}}{\alpha_{\mu\nu}}\Big{)}.$
(41)

We use this equation to filter the translation vectors in analytical FT (15).

The SR integral of four primitive $s$-type Gaussians on different centers has
analytical form

$\displaystyle(\mu\nu|\kappa\lambda)^{\mathrm{SR},\omega}$
$\displaystyle=\int\rho_{\mu\nu}(\mathbf{r}_{1})\frac{\mathrm{erfc}(\omega
r_{12})}{r_{12}}\rho_{\kappa\lambda}(\mathbf{r}_{2})d^{3}\mathbf{r}_{1}d^{3}\mathbf{r}_{2}$
$\displaystyle=\frac{8(\alpha_{\mu}\alpha_{\nu}\alpha_{\kappa}\alpha_{\lambda})^{3/4}}{\pi^{3}}\int
e^{\alpha_{\mu\nu}(\mathbf{r}_{1}-\mathbf{R}_{\mu\nu})^{2}}\frac{\mathrm{erfc}(\omega
r_{12})}{r_{12}}e^{\alpha_{\kappa\lambda}(\mathbf{r}_{2}-\mathbf{R}_{\kappa\lambda})^{2}}d\mathbf{r}_{1}^{3}d\mathbf{r}_{2}^{3}$
$\displaystyle=8\Big{(}\frac{\alpha_{\mu}\alpha_{\nu}\alpha_{\kappa}\alpha_{\lambda}}{\alpha_{\mu\nu}^{2}\alpha_{\kappa\lambda}^{2}}\Big{)}^{3/4}e^{-\eta_{\mu\nu}-\eta_{\kappa\lambda}}\frac{\mathrm{erf}(\sqrt{\beta}s)-\mathrm{erf}(\frac{\omega\sqrt{\beta}}{\sqrt{\omega^{2}+\beta}}s)}{s}$
(42)

where

$\displaystyle\beta=\frac{\alpha_{\mu\nu}\alpha_{\kappa\lambda}}{\alpha_{\mu\nu}+\alpha_{\kappa\lambda}},$
(43) $\displaystyle s=|\mathbf{R}_{\mu\nu}-\mathbf{R}_{\kappa\lambda}|,$ (44)
$\displaystyle\mathbf{R}_{\mu\nu}=\frac{\alpha_{\mu}\mathbf{R}_{\mu}+\alpha_{\nu}\mathbf{R}_{\nu}}{\alpha_{\mu}+\alpha_{\nu}}.$
(45)

In the $(DC|DC)$-type integrals, both $\rho_{\mu\nu}(\mathbf{r})$ and
$\rho_{\kappa\lambda}(\mathbf{r})$ are the product of a diffused Gaussian and
a compact Gaussian. $\mathbf{R}_{\mu\nu}$, the center of the product
$\rho_{\mu\nu}$, should be close to the center of the compact function
$\nu(\mathbf{r})$. So is the center $\mathbf{R}_{\kappa\lambda}$ in
$\rho_{\kappa\lambda}$. If function $\lambda(\mathbf{r})$ is close to function
$\nu(\mathbf{r})$, the distance $s$ between $\rho_{\mu\nu}(\mathbf{r})$ and
$\rho_{\kappa\lambda}(\mathbf{r})$ is relatively a small value otherwise the
SR Coulomb attenuation would rapidly decay to zero. Employing the formula

$\frac{\mathrm{erf}(ts)}{s}<\lim_{s\rightarrow
0}\frac{\mathrm{erf}(ts)}{s}=\frac{2t}{\sqrt{\pi}},$ (46)

the SR ERI is approximately

$\displaystyle(\mu\nu|\kappa\lambda)^{\mathrm{SR},\omega}$
$\displaystyle\approx
8\Big{(}\frac{\alpha_{\mu}\alpha_{\nu}\alpha_{\kappa}\alpha_{\lambda}}{\alpha_{\mu\nu}^{2}\alpha_{\kappa\lambda}^{2}}\Big{)}^{3/4}e^{-\eta_{\mu\nu}-\eta_{\kappa\lambda}}\frac{2}{\sqrt{\pi}}(\sqrt{\beta}-\frac{\omega\sqrt{\beta}}{\sqrt{\omega^{2}+\beta}})$
$\displaystyle<16\sqrt{\frac{\beta}{\pi}}\Big{(}\frac{\alpha_{\mu}\alpha_{\nu}\alpha_{\kappa}\alpha_{\lambda}}{\alpha_{\mu\nu}^{2}\alpha_{\kappa\lambda}^{2}}\Big{)}^{3/4}\exp\Big{(}-\frac{\alpha_{\mu}\alpha_{\nu}}{\alpha_{\mu\nu}}|\mathbf{R}_{\mu}-\mathbf{R}_{\mu\nu}|^{2}-\frac{\alpha_{\kappa}\alpha_{\lambda}}{\alpha_{\kappa\lambda}}|\mathbf{R}_{\kappa}-\mathbf{R}_{\kappa\lambda}|^{2})\Big{)}.$
(47)

Without loss of generality, we can assume $\alpha_{\mu}=\alpha_{\kappa}$ and
$\alpha_{\nu}=\alpha_{\lambda}$ and $\mu<\nu$

$(\mu\nu|\kappa\lambda)^{\mathrm{SR},\omega}<16\sqrt{\frac{\alpha_{\mu\nu}}{2\pi}}\Big{(}\frac{\alpha_{\mu}\alpha_{\nu}}{\alpha_{\mu\nu}^{2}}\Big{)}^{3/2}\exp\Big{(}-\frac{\alpha_{\mu}\alpha_{\nu}}{\alpha_{\mu\nu}}(|\mathbf{R}_{\mu}-\mathbf{R}_{\mu\nu}|^{2}+|\mathbf{R}_{\kappa}-\mathbf{R}_{\kappa\lambda}|^{2})\Big{)}.$
(48)

Due to the inequality (46), we would get the maximum value of (48) at
$s=\mathbf{R}_{\mu\nu}-\mathbf{R}_{\kappa\lambda}=0$. It indicates the
condition to filter the translation vector $\mathbf{M}$ in SR integration (16)

$16\sqrt{\frac{\alpha_{\mu\nu}}{2\pi}}\Big{(}\frac{\alpha_{\mu}\alpha_{\nu}}{\alpha_{\mu\nu}^{2}}\Big{)}^{3/2}\exp\Big{(}-\frac{\alpha_{\mu}\alpha_{\nu}|\mathbf{R}_{\mu}-(\mathbf{R}_{\kappa}+\mathbf{M})|^{2}}{2\alpha_{\mu\nu}}\Big{)}<\varepsilon.$
(49)

For $(CC|CC)$-type integrals, the orbital product $\rho$ decays rapidly. The
lattice sum in (16) is dominated by the terms of which $\mathbf{N}=0$ and
$\mathbf{L}=0$. For simplicity, we can analyze the special case $\mu=\nu$ and
$\kappa=\lambda$, and in addition the products $\rho_{\mu\nu}$ and
$\rho_{\kappa\lambda}$ being separated by the translation vector $\mathbf{M}$.
We then transform Eq. (42) to

$\displaystyle\frac{\mathrm{erfc}(\frac{\omega\sqrt{\beta}}{\sqrt{\omega^{2}+\beta}}s)-\mathrm{erfc}(\sqrt{\beta}s)}{s}<\frac{\sqrt{\omega^{2}+\alpha_{\mu}}}{\omega\sqrt{\alpha_{\mu}}}\frac{\exp(-\frac{\omega^{2}\alpha_{\mu}}{\omega^{2}+\alpha_{\mu}}s)}{2s^{2}}.$

Applying lattice sum to this equation, we get the filter condition

$\frac{\sqrt{\omega^{2}+\alpha_{\mu}}}{\omega\sqrt{\alpha_{\mu}}}\frac{\exp\Big{(}-\frac{\omega^{2}\alpha_{\mu}}{\omega^{2}+\alpha_{\mu}}(|\mathbf{M}|-\frac{\sqrt[3]{\Omega}}{2})^{2}\Big{)}}{|\mathbf{M}|-\frac{\sqrt[3]{\Omega}}{2}}.$
(50)

Depending on whether the system has diffused functions, we adopted either
condition (49) or condition (50) for the lattice sum index $\mathbf{M}$ in
(16).

In the $(DC|CC)$-type integrals, the maximum value of
$(\mu\nu|\kappa\lambda)^{\mathrm{SR},\omega}$ can be approximated with Eq.
(48). The maximum value of this integral appears around
$\mathbf{R}_{\kappa}=\mathbf{R}_{\lambda}\approx\mathbf{R}_{\nu}\approx\mathbf{R}_{\mu\nu}$
and $\alpha_{\kappa}=\alpha_{\lambda}$, i.e., the three compact basis
functions in the 4 centers are standing on the same atomic center and two of
which are the same basis functions, leading to the error estimation
proportional to the regular overlap of two basis functions

$\frac{4(\alpha_{\mu}\alpha_{\nu})^{3/4}}{\sqrt{\pi}\alpha_{\mu\nu}}\langle\mu|\nu^{\mathbf{N}}\rangle.$
(51)

By applying the effects of lattice sum (Appendix A), we would get the
condition to filter translation vector $\mathbf{N}$ in equation (16) for
$(DC|CC)$-type integrals (or vector $\mathbf{L}$ for $(CC|DC)$-type integrals)

$\frac{2\sqrt{\pi}|\mathbf{N}|}{(\alpha_{\mu}\alpha_{\nu})^{1/4}}\langle\mu|\nu^{\mathbf{N}}\rangle.$

##  III Benchmark

In this section, we benchmark the $\mathbf{k}$-point HF calculations. Our
program is available in PySCF-1.7.6Sun  _et al._ (2020). In these benchmark
calculations, the ERI prescreening parameter for J/K builder was set to
$10^{-13}$. The precision parameter $\varepsilon$ was set to $10^{-8}$ when
estimating energy cutoff (29), and the basis locality partitioning parameter
(II.4.4) and the lattice sum boundary of FT (41). All timings were obtained on
a desktop equipped with Intel I9-7940X CPU and 64 GB memory.

Table 1 lists the KHF (HF with $\mathbf{k}$-point sampling) results for
diamond crystal in the cubic unit cell with GTH-SZV basis and pseudo potential
GTH-pade. In each unit cell, this calculation contains 8 Carbon atoms and 32
basis functions. The largest calculation involves $15\times 15\times 15=3375$
samples in the first Brillouin zone, corresponding to in total 27000 atoms or
108000 basis functions. For comparison, evaluating the Coulomb integrals in
reciprocal space would require at least $N_{G}=49^{3}$ planewaves to converge
the integrals to the desired accuracy. The discrete FT engine in PySCF
requires at least $\frac{1}{2}\times 3\times 8\times
N_{G}^{4/3}N_{AO}^{2}N_{k}^{2}$ FLOPs to process the FT in K-matrix
constructor (factor 8 is obtained from the complex matrix multiplication; 3 is
due to the 3 passes in the 3D-FT; $\frac{1}{2}$ is due to the hermiticity in
the intermediate tensors). For $4\times 4\times 4$ $\mathbf{k}$-mesh, the cost
is roughly $2.9\times 10^{14}$ FLOPs.

When moving to all-electron calculations, as shown in Table 2, it is
inappropriate to call the discrete FT engine because of the steep basis
functions for core electrons. The STO-3G basis would require about $183^{3}$
planewaves for Coulomb repulsion integrals. This number becomes about
$954^{3}$ for cc-pVDZ basis due to the more compact functions presented in the
basis set. Even the minimal basis with minimal $2\times 2\times 2$
$\mathbf{k}$-mesh can easily reach $10^{15}$ FLOPs to build K-matrix, not to
mention the TB level memory footprint to hold the intermediates. The present
new algorithm would be able to finish the calculation of cc-pVDZ basis with
$4\times 4\times 4$ $\mathbf{k}$-mesh in a couple of hours on a dedicated
multi-core CPU machine. The largest all-electron calculation which consists of
112000 basis functions only requires 784 CPU hours and memory footprint less
than 6 GB.

Table 1: KHF for diamond crystal with GTH-SZV basis and pseudo potential GTH-pade. The unit cell is cubic structure with lattice parameter 3.5668 Å. $\mathbf{k}$ mesh |  $E$/$E_{\text{h}}$ per cell | $\omega$ | CPU hours | Error$\times 10^{9}$/$E_{\text{h}}$111Difference w.r.t. results of discrete FT formulaism
---|---|---|---|---
$1\times 1\times 1$ | -43.8468002862 | 2.512 | 2.76 | +0.13
$2\times 2\times 2$ | -44.0506288397 | 0.5 | 11.7 | +2.89
$3\times 3\times 3$ | -44.0240243649 | 0.5 | 13.4 | -0.75
$4\times 4\times 4$ | -44.0136700424 | 0.5 | 20.0 | -
$6\times 6\times 6$ | -44.0077526217 | 0.4 | 70.7 | -
$8\times 8\times 8$ | -44.0062794104 | 0.4 | 125 | -
$10\times 10\times 10$ | -44.0057532951 | 0.237 | 294 | -
$15\times 15\times 15$ |  |  |  | -
Table 2: All-electron KHF for diamond crystal with STO-3G and cc-pVDZ. The unit cell is cubic structure with lattice parameter 3.5668 Å. | STO-3G |  | cc-pVDZ
---|---|---|---
$\mathbf{k}$ mesh |  $E$/$E_{\text{h}}$ per cell | $\omega$ | CPU/h |  |  $E$/$E_{\text{h}}$ per cell | $\omega$ | CPU/h
$1\times 1\times 1$ | -299.328100921 | 0.6 | 2.82 |  | -302.870240306 | 0.5 | 54.3
$2\times 2\times 2$ | -299.551273686 | 0.6 | 2.58 |  | -303.076048519 | 0.609 | 19.3
$3\times 3\times 3$ | -299.525889782 | 0.579 | 6.17 |  | -303.054159127 | 0.609 | 30.1
$4\times 4\times 4$ | -299.516149849 | 0.579 | 16.1 |  | -303.044553455 | 0.616 | 79.5
$6\times 6\times 6$ | -299.510659684 | 0.332 | 25.3 |  | -303.038977867 | 0.474 | 260
$10\times 10\times 10$ | -299.508816002 | 0.237 | 95.2 |  | -303.037088495 | 0.237 | 784

##  IV Conclusions

We presented an algorithm to evaluate Coulomb repulsion integrals and
construct Fock matrix using periodic GTOs for all-electron crystals
simulations. In the algorithm, we split Coulomb integrals into smooth and
compact components based on the locality of underlying Coulomb operator and
basis functions. We then computed the intermediates either in real space or in
the reciprocal space and gain great efficiency in both space. The algorithm
does not have restrictions on the basis sets, the $\mathbf{k}$-point mesh, or
other preliminary knowledges of the system. HF self-consistency for systems of
10k atoms and 100k basis functions can be handled on regular desktops without
difficulty.

The performance of this algorithm relies on the thresholds and cutoffs in the
integrator. We discussed the prescreening thresholds and cutoffs in our
implementation based on the estimation of Coulomb integral values and lattice
sum formulaism. The estimations in this work are quite conservative. There are
plenty rooms to optimize the thresholds for better performance.

##  Appendix A Appendix

We demonstrate here how the lattice sum is applied to the estimation of
integral values.

  * •

Lattice sum for overlap

$\displaystyle\sum_{|\mathbf{L}|>R_{cut}}\exp\Big{(}-\frac{\alpha_{\mu}\alpha_{\nu}L^{2}}{\alpha_{\mu\nu}}\Big{)}$
$\displaystyle\approx\frac{1}{\Omega}\int_{R_{cut}-\frac{\sqrt[3]{\Omega}}{2}}^{\infty}\exp\Big{(}-\frac{\alpha_{\mu}\alpha_{\nu}L^{2}}{\alpha_{\mu\nu}}\Big{)}d^{3}\mathbf{L}$
$\displaystyle=\frac{4\pi}{\Omega}\int_{R_{cut}-\frac{\sqrt[3]{\Omega}}{2}}^{\infty}\exp\Big{(}-\frac{\alpha_{\mu}\alpha_{\nu}L^{2}}{\alpha_{\mu\nu}}\Big{)}L^{2}dL$
$\displaystyle=\frac{4\pi}{\Omega}\frac{\alpha_{\mu\nu}}{2\alpha_{\mu}\alpha_{\nu}}\Big{(}\big{(}{R_{cut}-\frac{\sqrt[3]{\Omega}}{2}}\big{)}\exp\Big{(}-\frac{\alpha_{\mu}\alpha_{\nu}(R_{cut}-\frac{\sqrt[3]{\Omega}}{2})^{2}}{\alpha_{\mu\nu}}\Big{)}+\int_{R_{cut}-\frac{\sqrt[3]{\Omega}}{2}}^{\infty}\exp\Big{(}-\frac{\alpha_{\mu}\alpha_{\nu}R^{2}}{\alpha_{\mu\nu}}\Big{)}dR\Big{)}$
$\displaystyle\approx\frac{2\pi(R_{cut}-\frac{\sqrt[3]{\Omega}}{2})\alpha_{\mu\nu}}{\Omega\alpha_{\mu}\alpha_{\nu}}\exp\Big{(}-\frac{\alpha_{\mu}\alpha_{\nu}(R_{cut}-\frac{\sqrt[3]{\Omega}}{2})^{2}}{\alpha_{\mu\nu}}\Big{)}\Big{(}1+\frac{\alpha_{\mu\nu}}{2\alpha_{\mu}\alpha_{\nu}(R_{cut}-\frac{\sqrt[3]{\Omega}}{2})^{2}}\Big{)}$
$\displaystyle<\frac{2\pi
R_{cut}\alpha_{\mu\nu}}{\Omega\alpha_{\mu}\alpha_{\nu}}\exp\Big{(}-\frac{\alpha_{\mu}\alpha_{\nu}(R_{cut}-\frac{\sqrt[3]{\Omega}}{2})^{2}}{\alpha_{\mu\nu}}\Big{)}$

  * •

Lattice sum for Coulomb integrals

$\displaystyle\int_{R_{cut}-\frac{\sqrt[3]{\Omega}}{2}}^{\infty}\frac{\sqrt{\omega^{2}+\alpha_{\mu}}}{\omega\sqrt{\alpha_{\mu}}}\frac{\exp(-\frac{\omega^{2}\alpha_{\mu}}{\omega^{2}+\alpha_{\mu}}s)}{2s^{2}}d^{3}\mathbf{s}$
$\displaystyle=\frac{2\pi\sqrt{\omega^{2}+\alpha_{\mu}}}{\omega\sqrt{\alpha_{\mu}}}\int_{R_{cut}-\frac{\sqrt[3]{\Omega}}{2}}^{\infty}\exp(-\frac{\omega^{2}\alpha_{\mu}}{\omega^{2}+\alpha_{\mu}}s)ds$
$\displaystyle=\frac{\omega\sqrt{\alpha_{\mu}}}{\sqrt{\pi}\sqrt{\omega^{2}+\alpha_{\mu}}}\mathrm{erfc}\Big{(}\frac{\omega\sqrt{\alpha_{\mu}}}{\sqrt{\omega^{2}+\alpha_{\mu}}}(R_{cut}-\frac{\sqrt[3]{\Omega}}{2})\Big{)}$
$\displaystyle\approx\frac{\sqrt{\omega^{2}+\alpha_{\mu}}}{\omega\sqrt{\alpha_{\mu}}}\frac{\exp\Big{(}-\frac{\omega^{2}\alpha_{\mu}}{\omega^{2}+\alpha_{\mu}}(R_{cut}-\frac{\sqrt[3]{\Omega}}{2})^{2}\Big{)}}{R_{cut}-\frac{\sqrt[3]{\Omega}}{2}}$

## References

  * Kresse and Joubert (1999) G. Kresse and D. Joubert, “From ultrasoft pseudopotentials to the projector augmented-wave method,” Phys. Rev. B 59, 1758–1775 (1999).
  * Enkovaara  _et al._ (2010) J. Enkovaara, C. Rostgaard, J. J. Mortensen, J. Chen, M. Dułak, L. Ferrighi, J. Gavnholt, C. Glinsvad, V. Haikola, H. A. Hansen, H. H. Kristoffersen, M. Kuisma, A. H. Larsen, L. Lehtovaara, M. Ljungberg, O. Lopez-Acevedo, P. G. Moses, J. Ojanen, T. Olsen, V. Petzold, N. A. Romero, J. Stausholm-Møller, M. Strange, G. A. Tritsaris, M. Vanin, M. Walter, B. Hammer, H. Häkkinen, G. K. H. Madsen, R. M. Nieminen, J. K. Nørskov, M. Puska, T. T. Rantala, J. Schiøtz, K. S. Thygesen, and K. W. Jacobsen, “Electronic structure calculations with GPAW: a real-space implementation of the projector augmented-wave method,” J. Phys.: Condens. Matter 22, 253202 (2010).
  * Royo and Stengel (2019) M. Royo and M. Stengel, “First-principles theory of spatial dispersion: Dynamical quadrupoles and flexoelectricity,” Phys. Rev. X 9, 021050 (2019).
  * Giannozzi  _et al._ (2020) P. Giannozzi, O. Baseggio, P. Bonfà, D. Brunato, R. Car, I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas, F. Ferrari Ruffino, A. Ferretti, N. Marzari, I. Timrov, A. Urru, and S. Baroni, “Quantum espresso toward the exascale,” J. Chem. Phys. 152, 154105 (2020), https://doi.org/10.1063/5.0005082 .
  * Dovesi  _et al._ (2018) R. Dovesi, A. Erba, R. Orlando, C. M. Zicovich-Wilson, B. Civalleri, L. Maschio, M. Rérat, S. Casassa, J. Baima, S. Salustro, and B. Kirtman, “Quantum-mechanical condensed matter simulations with crystal,” WIREs Comput. Mol. Sci. 8, e1360 (2018), https://onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1360 .
  * Kühne  _et al._ (2020) T. D. Kühne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöß, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack, and J. Hutter, “Cp2k: An electronic structure and molecular dynamics software package - quickstep: Efficient and accurate electronic structure calculations,” J. Chem. Phys. 152, 194103 (2020), https://doi.org/10.1063/5.0007045 .
  * Sun  _et al._ (2020) Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, “Recent developments in the pyscf program package,” J. Chem. Phys. 153, 024109 (2020), https://doi.org/10.1063/5.0006074 .
  * Bowler and Miyazaki (2010) D. R. Bowler and T. Miyazaki, “Calculations for millions of atoms with density functional theory: linear scaling shows its potential,” J. Phys.: Condens. Matter 22, 074207 (2010).
  * VandeVondele, Borštnik, and Hutter (2012) J. VandeVondele, U. Borštnik, and J. Hutter, “Linear scaling self-consistent field calculations with millions of atoms in the condensed phase,” J. Chem. Theory Comput. 8, 3565–3573 (2012), pMID: 26593003, https://doi.org/10.1021/ct200897x .
  * Pisani and Dovesi (1980) C. Pisani and R. Dovesi, “Exact-exchange hartree-fock calculations for periodic systems. i. illustration of the method,” Int. J. Quantum Chem. 17, 501–516 (1980), https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.560170311 .
  * Dovesi  _et al._ (1983) R. Dovesi, C. Pisani, C. Roetti, and V. R. Saunders, “Treatment of coulomb interactions in hartree-fock calculations of periodic systems,” Phys. Rev. B 28, 5781–5792 (1983).
  * Gygi and Baldereschi (1986) F. Gygi and A. Baldereschi, “Self-consistent hartree-fock and screened-exchange calculations in solids: Application to silicon,” Phys. Rev. B 34, 4405–4408 (1986).
  * Heyd and Scuseria (2004) J. Heyd and G. E. Scuseria, “Efficient hybrid density functional calculations in solids: Assessment of the heyd-scuseria-ernzerhof screened coulomb hybrid functional,” J. Chem. Phys. 121, 1187–1192 (2004), https://doi.org/10.1063/1.1760074 .
  * Izmaylov, Scuseria, and Frisch (2006) A. F. Izmaylov, G. E. Scuseria, and M. J. Frisch, “Efficient evaluation of short-range hartree-fock exchange in large molecules and periodic systems,” J. Chem. Phys. 125, 104103 (2006), https://doi.org/10.1063/1.2347713 .
  * Spencer and Alavi (2008) J. Spencer and A. Alavi, “Efficient calculation of the exact exchange energy in periodic systems using a truncated coulomb potential,” Phys. Rev. B 77, 193110 (2008).
  * Guidon, Hutter, and VandeVondele (2009) M. Guidon, J. Hutter, and J. VandeVondele, “Robust periodic hartree-fock exchange for large-scale simulations using gaussian basis sets,” J. Chem. Theory Comput. 5, 3010–3021 (2009), pMID: 26609981, https://doi.org/10.1021/ct900494g .
  * Sundararaman and Arias (2013) R. Sundararaman and T. A. Arias, “Regularization of the coulomb singularity in exact exchange by wigner-seitz truncated interactions: Towards chemical accuracy in nontrivial systems,” Phys. Rev. B 87, 165122 (2013).
  * Marzari and Vanderbilt (1997) N. Marzari and D. Vanderbilt, “Maximally localized generalized wannier functions for composite energy bands,” Phys. Rev. B 56, 12847–12865 (1997).
  * Gygi (2009) F. Gygi, “Compact representations of kohn-sham invariant subspaces,” Phys. Rev. Lett. 102, 166406 (2009).
  * Marzari  _et al._ (2012) N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, “Maximally localized wannier functions: Theory and applications,” Rev. Mod. Phys. 84, 1419–1475 (2012).
  * Gygi and Duchemin (2013) F. Gygi and I. Duchemin, “Efficient computation of hartree-fock exchange using recursive subspace bisection,” J. Chem. Theory Comput. 9, 582–587 (2013).
  * DiStasio  _et al._ (2014) R. A. DiStasio, B. Santra, Z. Li, X. Wu, and R. Car, “The individual and collective effects of exact exchange and dispersion interactions on the ab initio structure of liquid water,” J. Chem. Phys. 141, 084502 (2014).
  * Hu  _et al._ (2017) W. Hu, L. Lin, A. S. Banerjee, E. Vecharynski, and C. Yang, “Adaptively compressed exchange operator for large-scale hybrid density functional calculations with applications to the adsorption of water on silicene,” J. Chem. Theory Comput. 13, 1188–1198 (2017), pMID: 28177229, https://doi.org/10.1021/acs.jctc.6b01184 .
  * Hu, Lin, and Yang (2017) W. Hu, L. Lin, and C. Yang, “Interpolative separable density fitting decomposition for accelerating hybrid density functional calculations with applications to defects in silicon,” J. Chem. Theory Comput. 13, 5420–5431 (2017).
  * Towler, Zupan, and Causà (1996) M. D. Towler, A. Zupan, and M. Causà, “Density functional theory in periodic systems using local gaussian basis sets,” Comput. Phys. Commun. 98, 181–205 (1996).
  * Schwegler and Challacombe (1999) E. Schwegler and M. Challacombe, “Linear scaling computation of the fock matrix. iv. multipole accelerated formation of the exchange matrix,” J. Chem. Phys. 111, 6223–6229 (1999).
  * Heyd, Scuseria, and Ernzerhof (2003) J. Heyd, G. E. Scuseria, and M. Ernzerhof, “Hybrid functionals based on a screened coulomb potential,” J. Chem. Phys. 118, 8207–8215 (2003).
  * Tymczak  _et al._ (2005) C. J. Tymczak, V. T. Weber, E. Schwegler, and M. Challacombe, “Linear scaling computation of the fock matrix. viii. periodic boundaries for exact exchange at the $\gamma$ point,” J. Chem. Phys. 122, 124105 (2005).
  * Maschio  _et al._ (2007) L. Maschio, D. Usvyat, F. R. Manby, S. Casassa, C. Pisani, and M. Schütz, “Fast local-mp2 method with density-fitting for crystals. i. theory and algorithms,” Phys. Rev. B 76, 075101 (2007).
  * Varga (2008) v. Varga, “Long-range analysis of density fitting in extended systems,” Int. J. Quantum Chem. 108, 1518–1527 (2008), https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.21682 .
  * Maschio and Usvyat (2008) L. Maschio and D. Usvyat, “Fitting of local densities in periodic systems,” Phys. Rev. B 78, 073102 (2008).
  * Ben, Hutter, and VandeVondele (2013) M. D. Ben, J. Hutter, and J. VandeVondele, “Electron correlation in the condensed phase from a resolution of identity approach based on the gaussian and plane waves scheme,” J. Chem. Theory Comput. 9, 2654–2671 (2013), http://dx.doi.org/10.1021/ct4002202 .
  * Yamada, Brewster, and Hirata (2013) T. Yamada, R. P. Brewster, and S. Hirata, “Asymptotic expansion of two-electron integrals and its application to coulomb and exchange lattice sums in metallic, semimetallic, and nonmetallic crystals,” J. Chem. Phys. 139, 184107 (2013).
  * Sun  _et al._ (2017) Q. Sun, T. C. Berkelbach, J. D. McClain, and G. K.-L. Chan, “Gaussian and plane-wave mixed density fitting for periodic systems,” J. Chem. Phys. 147, 164119 (2017).
  * Irmler, Burow, and Pauly (2018) A. Irmler, A. Burow, and F. Pauly, “Robust periodic fock exchange with atom-centered gaussian basis sets,” J. Chem. Theory Comput. 14, 4567–4580 (2018).
  * Wang, Lewis, and Valeev (2020) X. Wang, C. A. Lewis, and E. F. Valeev, “Efficient evaluation of exact exchange for periodic systems via concentric atomic density fitting,” J. Chem. Phys. 153, 124116 (2020).
  * Hartwigsen, Goedecker, and Hutter (1998) C. Hartwigsen, S. Goedecker, and J. Hutter, “Relativistic separable dual-space gaussian pseudopotentials from h to rn,” Phys. Rev. B 58, 3641–3662 (1998).
  * Willand  _et al._ (2013) A. Willand, Y. O. Kvashnin, L. Genovese, A. Vázquez-Mayagoitia, A. K. Deb, A. Sadeghi, T. Deutsch, and S. Goedecker, “Norm-conserving pseudopotentials with chemical accuracy compared to all-electron calculations,” J. Chem. Phys. 138, 104109 (2013).
  * Gao and Chelikowsky (2018) W. Gao and J. R. Chelikowsky, “Accuracy of partial core corrections using fourier transforms in pseudopotential-density functional theory,” J. Chem. Theory Comput. 14, 6515–6520 (2018).
  * Prandini  _et al._ (2018) G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet, and N. Marzari, “Precision and efficiency in solid-state pseudopotential calculations,” npj Comput. Mater. 4, 72 (2018).
  * Lu  _et al._ (2019) J.-B. Lu, D. C. Cantu, M.-T. Nguyen, J. Li, V.-A. Glezakou, and R. Rousseau, “Norm-conserving pseudopotentials and basis sets to explore lanthanide chemistry in complex environments,” J. Chem. Theory Comput. 15, 5987–5997 (2019).
  * Dominik Marx (2009) J. H. Dominik Marx, “Ab initio molecular dynamics. basic theory and advanced methods. by dominik marx and jürg hutter.” Angew. Chem. Int. Ed. 48, 9404–9405 (2009).