Accurate approximations of density functional theory
for large systems with applications to defects in
crystalline solids
Kaushik Bhattacharya
1
, Vikram Gavini
2
, Michael Ortiz
1
, Mauricio Ponga
3
,
and Phanish Suryanarayana
4
1
California Institute of Technology, Pasadena,USA
2
University of Michigan, Ann Arbor,USA
3
University of British Columbia, Vancouver, Canada
4
Georgia Institute of Technology, Atlanta, USA
December 14, 2021
Abstract
This chapter presents controlled approximations of Kohn-Sham density functional
theory (DFT) that enable very large scale simulations. The work is motivated by the
study of defects in crystalline solids, though the ideas can be used in other applica-
tions. The key idea is to formulate DFT as a minimization problem over the density
operator, and to cast spatial and spectral discretization as systematically convergent
approximations. This enables efficient and adaptive algorithms that solve the equations
of DFT with no additional modeling, and up to desired accuracy, for very large systems,
with linear and sublinear scaling. Various approaches based on such approximations
are presented, and their numerical performance demonstrated through selected exam-
ples. These examples also provide important insight about the mechanics and physics
of defects in crystalline solids.
Contents
1 Introduction
3
2 Variational formulation of density functional theory
5
2.1 Kohn-Sham density functional theory as a minimum problem . . . . . . . . . 5
2.2 Approximations resulting from spatial discretization . . . . . . . . . . . . . . 8
2.3 Spectral reformulation of the discrete Kohn-Sham problem . . . . . . . . . . 9
2.4 Approximation by numerical quadrature . . . . . . . . . . . . . . . . . . . . 10
2.5 Rayleigh-Ritz interpretation . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1
arXiv:2112.06016v1 [cond-mat.mtrl-sci] 11 Dec 2021
2.6 Convexification and thermalization . . . . . . . . . . . . . . . . . . . . . . . 11
2.7 Spatial densities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.8 Eigenvalue problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3 Filtering, spectrum splitting and pseudopotentials
13
3.1 Filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2 Spectrum splitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.3 Frozen core approximation and pseudopotentials . . . . . . . . . . . . . . . . 16
4 Spatial coarse-graining: Finite-element discretization
17
4.1 Higher-order spectral finite-elements . . . . . . . . . . . . . . . . . . . . . . 18
4.2 Spatial adaptivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.3 DFT-FE: A massively parallel code for real-space finite-element DFT calcu-
lations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.4 Enriched finite-element basis . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
5 Spectral coarse-graining: Spectral Quadrature method
24
5.1 Spectral integrals and quadrature . . . . . . . . . . . . . . . . . . . . . . . . 25
5.2 Spectral integrals and quadrature with spatial localization . . . . . . . . . . 27
5.3 Gauss Spectral Quadrature . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
5.4 Clenshaw-Curtis Spectral Quadrature . . . . . . . . . . . . . . . . . . . . . . 31
5.5 Convergence rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5.6 Scaling estimates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5.7 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
6 Spatial and spectral coarse-graining
35
6.1 Periodic systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
6.2 Coarse-grained representation . . . . . . . . . . . . . . . . . . . . . . . . . . 36
6.3 Selected results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
A Crystalline solids and the Cauchy-Born Rule
42
2
1 Introduction
Defects are common in crystalline solids [63, 35]. These include point defects (vacancies
with missing atoms, substitutional elements where an atom of an impurity (solute) replaces
the atom of the actual material, or interstitial atom where an extra atom is inserted into
the solid), cluster defects (vacancy cluster or prismatic dislocation loop with a missing disc
of atoms), line defects (dislocations where an extra plane of atoms terminates along a line)
or planar defects (twin or grain boundaries across which the crystal orientation changes or
phase boundaries across which the crystal structure changes).
Defects play a critical role in determining important properties of materials. Vacancies
mediate creep, solutes strengthen solids, vacancy clusters lead to void nucleation and dis-
locations mediate plasticity. Remarkably, they do so at extremely dilute concentrations.
Vacancies affect creep at parts per million, and dislocations densities are of the order of one
in a million amongst atomic columns during plastic deformation.
The reason that defects can have such a profound effect on properties at dilute concen-
trations is because they trigger physics at multiple length and time scales [63, 78]. In this
review, we are interested in the equilibrium structure, and therefore focus only on the length
scales. The defects cause an imbalance of forces on the neighboring atoms which in turn
lead to deformations. Even though electronic interactions decay quickly, displacement of the
atoms from their periodic equilibrium positions lead to imbalanced forces on their neighbors,
and so on and on, leading to extremely slow decay of the displacement field. The complex
quantum mechanical or chemical interactions at the defect core lead to a complex atomistic
and electronic structure that need an electronic structure theory for its description. As we
move away, the displacements from the periodic structure are less complex and may be un-
derstood through atomistic interactions. Even farther away, the displacements are smaller
and may be described by continuum elasticity theory.
Crucially, these scales interact intimately and one scale does not dominate over others.
We consider two examples. The first example is a vacancy, where a single atom is missing
from an otherwise perfect lattice. In the far field, elasticity theory tells us that the displace-
ment decays as
r
−
2
so that the stress and the strain decay as
r
−
3
[57]. This is relatively fast
decay and the energy is summable in the far field. The divergence at the origin in regularized
by atomistic and electronic interactions. One can estimate the energy due to the chemistry
of the core by the cohesive energy of a solid (the energy difference between an isolated atom
and an atom in the crystal), and this is typically few electron volts in metals. One can also
estimate the energy due to the elastic field far away, and this is of the order of a few tenths
of electron volts for a typical metal. While one is smaller than the other, it is not negligibly
so; therefore these fields do interact with each other even in this simple defect. Further, the
elastic fields generated by the vacancy are large enough to interact with macroscopic stress
due to boundary conditions resulting in stress-induced driving force on the vacancy. The
second example is that of a dislocation. The stress and strain decay as
r
−
1
away from the
line which means that the elastic energy density is logarithmic, and thus divergent at both
the origin and the far field [35, 57]. In other words, dislocations provide a direct link between
continuum scale boundary conditions and electronic scale interactions at the core. In short,
defects connect the far field to the electronic scale and it is this ability to bridge scales that
result in defects having a profound effect on macroscopic properties. This also makes defects
3
extremely difficult to study.
Kohn-Sham DFT [36, 44] has emerged as the method of choice for the study of electronic
structure in condensed matter [11]. It converts the many particle Schödinger equation to a
single particle problem with an effective single-electron potential. While one can show the
existence of such a theory, the functional that gives rise to the potential (or even the locality
or lack thereof) remains unknown, and is modelled. Thus, while DFT is nominally
ab initio
in that it is agnostic about the material (other than the atomic number), it does require
a constitutive model of the universal exchange-correlation functional. It has proved to be
an extremely useful compromise between practical application of quantum mechanics and
fidelity.
Given the importance of electronic structure in determining the structure and energetics
of defects, it is desirable to study defects using DFT. Such studies require large domains
including large numbers of electrons due to the slow decay and low concentrations that are
typical in materials; unfortunately such large domains are beyond the capability of brute
force (full-resolution) DFT calculations using existing widely-used methods with reasonable
computational resources. However, the complexity and details of the electronic structure are
important near the core, and less so in the large regions of slow decay. This has motivated
multiscale modeling of materials, where one builds a cascade of models (DFT, atomistic,
continuum) to study the phenomena at different scales [1, 89, 6, 94]. However, the interaction
between the scales means that one has to link them, which requires
further modeling
. Much
of this modeling is empirical, taking us away from the
ab initio
point of view. Moreover,
such a cascade of models linked empirically do not have an inherent or quantifiable notion
of error. A notion of error is important since the study of defects requires comparing the
energies of different configurations, and therefore evaluating a small difference between two
large numbers. Therefore an estimate of the accuracy and an ability to control the accuracy
is important.
This chapter presents a line of work that seeks to
solve the equations of DFT, and only
the equations of DFT, with no further modeling
on large domains relevant to defects by in-
troducing approximations where the error can be controlled. The idea is to formulate DFT
as a well-defined minimum problem over spaces of operators, and then introduce system-
atically convergent approximations. Specifically, there are two approximations: i) Spatial
discretization, resulting in finite-dimensional approximate problems obtained by constrained
minimization of the DFT energy functional; and ii) Spectral discretization, based on ap-
proximations of the DFT energy functional itself, or
variational crimes
in the parlance of
approximation theory.
A related issue is the fact that widely used DFT methods for condensed matter are limited
to periodic systems. This is motivated by crystals that involve a periodic arrangement of
atoms. The periodicity enables one to work in Fourier space using a plane wave basis, and
this has proven to be extremely efficient on moderate computational resources. One can
extend this approach to defects using “super-cells”, i.e., studying a periodic arrangement
of defects and the resulting unit cell. However, since defects interact over long distances,
these can lead to artifacts. Further, dislocations are topological defects and therefore not
amenable to periodic arrangement unless one studies defect complexes (dislocation dipoles or
quadrapoles) with zero topological content, further leading to potential artifacts. Therefore,
it is desirable to move away from periodic arrangements. Finally, defects are interesting
4
because they interact with far-field stimuli. Therefore, it is desirable to study defects under
arbitrary boundary conditions. These motivate the need to solve DFT in real space.
The systematically convergent approximations lead to a various algorithms that enable
the solution of DFT with controlled error on large systems in real space. This chapter de-
scribes three in some detail. The first involves variable spatial discretization by exploiting
adaptive higher finite elements. The second involves spectral discretization using quadra-
tures. The final method combines spectral and spatial discretization.
The chapter is organized as follows. Section 2 provides a variational formulation of DFT.
An important result is the reformulation (13) which presents DFT as a nested variational
problem. We then use spectral theory to rewrite the inner variational problem. This reformu-
lation enables us to introduce spatial discretization and spectral quadratures as convergent
(Rayleigh-Ritz) approximations. Section 3 introduces three ideas that are useful for the ef-
ficient practical implementation of the methods. Section 4 introduces spatial discretization
using (higher order) finite-elements, and describes how this can be used for spatial adaptivity.
This section also presents a series of examples that describe the efficacy of such an approach
in studying defects, and the overall performance of the method. We turn to spectral quadra-
tures in Section 5. We discuss the relation of this method to other approaches including the
recursion method (widely used in tight binding), Padé approximations and Fermi-operator
expansion. We discuss convergence and demonstrate the performance of the method using
various examples. We turn to combined spatial and spectral coarse-graining in Section 6.
We describe a sub-linear scaling method for the study of defects, and its application to study
vacancy clusters and dislocation cores in magnesium.
2 Variational formulation of density functional theory
In this section, we proceed formally following the notation and presentation of Refs. [3, 88]
to formulate Kohn-Sham DFT as a well-defined minimum problem over spaces of operators,
and the approximation schemes – spatial and spectral – that it suggests.
2.1 Kohn-Sham density functional theory as a minimum problem
We consider a closed shell spin-unpolarized system in an insulated, bounded, open and
Lipschitz domain
Ω
⊂
R
3
for simplicity. The presentation may be extended to spin-polarized
systems and unbounded domains [3]. Let
V
=
H
1
0
(Ω)
,
H
=
L
2
(Ω;
C
)
,
D
(
h
) =
H
1
0
(Ω;
C
)
,
and
X
=
{
γ
∈ S
(
H
)
,
R
(
γ
)
⊂ D
(
h
)
}
, where
h
denotes the single-particle Hamiltonian,
γ
represents the single-particle density matrix, and
S
(
H
)
denotes the vector space of bounded
self-adjoint operators on
H
. Let
K
=
{
γ
∈X
: 0
≤
γ
≤
1
,
2tr(
γ
) =
N
}
(1)
be a constraint set defining the admissible density operators. Define the Kohn-Sham energy
functional
F
:
X →
̄
R
as
F
(
γ
) =
{
2tr
(
−
1
2
∆
γ
)
+
G
(
ρ
)
,
if
γ
∈K
,
+
∞
,
otherwise
,
(2)
5
which can be written as
F
(
γ
) =
I
K
+ 2tr
(
−
1
2
∆
γ
)
+
G
(
ρ
)
,
(3)
where
ρ
(
r
) = 2
γ
(
r
,
r
′
)
(4)
is the electron density, with
σ
denoting the spin. In addition,
G
(
ρ
) =
∫
Ω
ρv d
r
+
J
(
ρ
) +
E
xc
(
ρ
)
,
(5)
where
v
is an external potential, and
J
(
ρ
) = sup
φ
∈V
{
∫
Ω
ρφd
r
−
1
4
π
∫
Ω
1
2
|∇
φ
|
2
d
r
}
(6)
is the classical electrostatic energy [39], A formal connection between (6) and the oft-used
equivalent expression (up to an inessential constant) based on the Coulombic interaction
formula
1
2
∫
Ω
∫
Ω
ρ
(
r
)
ρ
(
r
′
)
|
r
−
r
′
|
d
r
d
r
′
can be established simply by writing out the Euler-Lagrange
equation of (6) and solving for
φ
using the Green’s function for the Laplacian. The expression
(6) simply recognizes the fact that
J
(
ρ
)
is the dual of the Dirichlet functional. Representation
(6) is advantageous over the Coulombic representation from the standpoint of approximation,
which only requires local conforming interpolation of the electrostatic field
φ
. Finally,
E
xc
(
ρ
)
is the exchange-correlation energy functional. It must necessarily be modeled. Here, for
simplicity, we choose the local density approximation (LDA) [59]. The
Kohn-Sham DFT
problem
is to find the ground state energy
E
KS
= inf
γ
∈X
F
(
γ
)
,
(7)
and attendant energy-minimizing states. For subsequent purposes, we use duality to refor-
mulate the energy in trace form. To this end, assume for definiteness that the exchange
correlation
E
xc
:
V →
̄
R
is convex with dual
E
∗
xc
:
V →
̄
R
, so that [88]
E
xc
(
ρ
) = sup
v
xc
∈V
{
ρu
−
E
∗
xc
(
v
xc
)
}
.
(8)
Now, define the Hamiltonian
h
(
φ,v
xc
) =
−
1
2
∆ +
Φ
+
V
xc
+
V,
(9)
where the electrostatic potential operator
Φ
, the exchange-correlation potential operator
V
xc
,
and the external potential operator
V
are bounded self-adjoint operators over
H
defined by
the properties
2tr(
Φγ
) =
∫
Ω
ρφd
r
,
2tr(
V
xc
γ
) =
∫
Ω
ρv
xc
d
r
,
2tr(
V γ
) =
∫
Ω
ρv d
r
.
(10)
6
Then,
F
(
γ
) =
I
K
(
γ
) + sup
φ
∈V
sup
v
xc
∈V
(
2tr(
h
(
φ,v
xc
)
γ
)
−
1
4
π
∫
Ω
1
2
|∇
φ
|
2
d
r
−
E
∗
xc
(
v
xc
)
)
(11)
and the Kohn-Sham DFT problem (7) becomes
E
KS
= inf
γ
∈X
sup
φ
∈V
sup
v
xc
∈V
(
I
K
(
γ
) +
(
2tr(
h
(
φ,v
xc
)
γ
)
−
1
4
π
∫
Ω
1
2
|∇
φ
|
2
d
r
−
E
∗
xc
(
v
xc
)
)
.
(12)
It is possible to exchange the order of the
inf
and
sup
operations in the above equation [88]
to arrive at the
reformulated Kohn-Sham DFT problem
,
E
KS
= sup
φ
∈V
sup
v
xc
∈V
[
inf
γ
∈X
(
I
K
(
γ
) + 2tr(
h
(
φ,v
xc
)
γ
)
)
−
1
4
π
∫
Ω
1
2
|∇
φ
|
2
d
r
−
E
∗
xc
(
v
xc
)
]
.
(13)
This reformulation offers various advantages and serves as the basis for the approximations
to follow. First, in the same spirit as in (6), the representation (13) only involves local
operators and requires local or conforming interpolation of
φ
and
v
xc
. Second, the functional
be expressed in terms of linear operators acting on
γ
only, thus paving the way for a spectral
treatment of the problem as we do presently.
We now focus on the inner
inf
operation that yields the energy-minimizing density matrix
for fixed
(
φ,u
)
:
U
= inf
γ
∈X
(
I
K
(
γ
) + 2tr(
hγ
)
)
,
(14)
where, here and subsequently, we omit the dependence of
h
on the fixed fields
(
φ,v
xc
)
for
simplicity of notation. Note that the quantity
U
is commonly referred to as the
band structure
energy
in the physics literature.
It follows from spectral theory (cf., e. g., [68]) that the minimizing density matrix operator
γ
of (14) shares the same spectral measure as the Hamiltonian
h
, i.e., we may write
h
=
∫
R
εd
E
(
ε
) =
h
T
, γ
=
∫
R
f
(
ε
)
d
E
(
ε
) =
γ
T
,
(15)
for
0
≤
f
≤
1
, where
E
is a resolution of the identity over the Borel sets of the real line. In
addition,
γ
and
h
have the same spectral measure if and only if they commute, i. e.,
γh
=
hγ.
(16)
Finally, we can show that there is a minimizer
f
∈ {
0
,
1
}
so that
γ
2
=
γ
. Therefore, the
minimum problem is
minimize:
U
(
γ
) := 2tr(
hγ
)
,
(17a)
subject to:
γ
T
=
γ, γh
=
hγ,γ
2
=
γ,
2tr(
γ
) =
N.
(17b)
The variational problem (13) is often solved by a fixed point iteration of
self-consistent
field
(SCF) iteration where (17) is solved for
γ
with a fixed
φ,v
xc
and the outer sup problem
in (13) is used to update
φ
and
v
xc
for fixed
γ
.
7
2.2 Approximations resulting from spatial discretization
We proceed to discretize problem (13)
à la
Rayleigh-Ritz, i.e., by restriction to finite-
dimensional spaces. To this end, let
V
h
be a nested sequence of finite dimensional spaces
of
V
spanned by orthonormal bases
{
e
1
,...,e
N
g
}
, e.g., corresponding to a finite element
discretization
1
. Let
H
h
=
V
h
be the corresponding sequence of subspaces of
H
. Then, the
discrete wave function, electrostatic field and exchange-correlation potential field are of the
form
φ
h
(
r
) =
N
g
∑
a
=1
φ
a
e
a
(
r
)
, φ
h
(
r
) =
N
g
∑
a
=1
φ
a
e
a
(
r
)
, v
xc,h
(
r
) =
N
g
∑
a
=1
v
xc,a
e
a
(
r
)
.
(18)
Likewise, the discrete Hamiltonian is
h
h
(
r
,
r
′
) =
N
g
∑
a
1
=1
N
g
∑
a
2
=1
H
a
1
a
2
e
a
1
(
r
)
e
a
2
(
r
′
)
,
(19)
with
H
=
A
h
+
Φ
h
+
V
xc,h
+
V
h
,
(20)
A
h
;
a
1
a
2
=
∫
Ω
1
2
∇
e
a
1
(
r
)
·∇
e
a
2
(
r
)
d
r
=
B
h
;
a
1
a
2
,
(21)
Φ
h
;
a
1
a
2
=
∫
Ω
(
N
g
∑
a
=1
φ
a
e
a
(
r
)
)
e
a
1
(
r
)
e
a
2
(
r
)
d
r
,
(22)
V
xc,h
;
a
1
a
2
=
∫
Ω
(
N
g
∑
a
=1
v
xc,a
e
a
(
r
)
)
e
a
1
(
r
)
e
a
2
(
r
)
d
r
,
(23)
and
V
h
;
a
1
a
2
=
∫
Ω
v
(
r
)
e
a
1
(
r
)
e
a
2
(
r
)
d
r
.
(24)
We note the additional linear structure
Φ
h
=
T
h
φ
h
, V
xc,h
=
T
h
u
h
,
(25)
where
T
a
1
a
2
a
3
=
∫
Ω
e
a
1
(
r
)
e
a
2
(
r
)
e
a
3
(
r
)
d
r
.
(26)
Finally, the discrete density matrix is of the form
γ
h
(
r
,
r
′
) =
N
g
∑
a
1
=1
N
g
∑
a
2
=1
P
a
1
a
2
e
a
1
(
r
)
e
a
2
(
r
′
)
,
(27)
and the discrete electron density follows as
ρ
h
(
r
) = 2
γ
h
(
r
,
r
)
.
(28)
1
Note that we use the subscript
h
to index the nested spaces following typical notation in computational
science, and not to signify a relationship with the Hamiltonian which we have denoted using by the letter
h
.
8
This sequence
{V
k
}
of finite-dimensional subspaces of
V
defines a nested sequence of
subspaces of
X
of
density matrices
X
h
=
{
γ
h
∈ S
(
H
h
)
}
, where
S
(
H
h
)
denotes the vector
space of symmetric linear operators on
H
h
. This in turn defines a sequence of discrete
constraint sets
K
h
=
{
γ
h
∈ X
h
: 0
≤
γ
h
≤
1
,
2tr(
γ
h
) =
N
}
, where
0
≤
γ
h
≤
1
expresses
the requirement that
0
≤
(
φ
h
|
γ
h
|
φ
h
)
≤
1
for all
φ
h
∈ H
h
. We note that, if the spaces
H
h
are nested, then
K
h
defines a decreasing sequence of sets in
X
and that
K⊂K
h
. Then, the
corresponding sequence of discrete energies
F
h
:
X
h
→
̄
R
follows as
F
h
(
γ
h
) =
I
K
h
(
γ
h
) + sup
φ
h
∈V
h
sup
v
xc,h
∈V
h
(
2tr(
h
h
γ
h
)
−
1
8
π
(
φ
h
|
B
h
|
φ
h
)
−
E
∗
xc
(
v
xc,h
)
)
,
(29)
where
I
K
h
is the indicator function of
K
h
, and
B
h
=
D
∗
h
D
h
with
D
h
the discrete gradient
operator. The discrete Kohn-Sham DFT problem becomes
E
KS
,h
= inf
γ
h
∈X
h
F
h
(
γ
h
)
= sup
φ
h
∈V
h
sup
v
xc,h
∈V
h
[
inf
γ
h
∈X
h
(
I
K
h
(
γ
h
) + 2tr(
h
h
γ
h
)
)
−
1
8
π
(
φ
h
|
B
h
|
φ
h
)
−
E
∗
xc
(
v
xc,h
)
]
,
(30)
where we have again exchanged the order of the inf and sup operations [88].
As before, we may rewrite the inner inf problem as
minimize:
U
h
(
γ
h
) := 2tr(
h
h
γ
h
)
,
(31a)
subject to:
γ
T
h
=
γ
h
, γ
h
h
h
=
h
h
γ
h
,γ
2
h
=
γ
h
,
2tr(
γ
h
) =
N.
(31b)
2.3 Spectral reformulation of the discrete Kohn-Sham problem
By the spectral decomposition theorem (cf., e. g., [68]), we can write
H
=
∫
R
εd
E
h
(
ε
)
, P
=
∫
R
f
h
(
ε
)
d
E
h
(
ε
)
(32)
where
E
h
is an operator valued measure. In this representation, we have
tr(
HP
) =
∫
R
εf
h
(
ε
)
d
M
h
(
ε
)
≡
F
h
(
f
h
)
,
(33)
and
tr(
P
) =
∫
R
f
h
(
ε
)
d
M
h
(
ε
)
≡
M
h
(
f
h
)
,
(34)
where
M
h
= tr(
E
h
)
(35)
is a
spectral measure
with:
d
M
h
=
N
g
∑
i
=1
δ
ε
i
,
(36)
where
δ
is the Dirac delta function. Given the spectral measure
M
h
, the calculation of the
energy-minimizing discrete density matrix
γ
h
at fixed
(
φ
h
,u
h
)
reduces to the scalar problem
inf
f
h
∈
X
h
{
F
h
(
f
h
)
,
0
≤
f
≤
1
,
2
M
h
(
f
h
) =
N
}
,
(37)
where
X
h
denotes the space of bounded real-valued Borel functions over the real line.
9
2.4 Approximation by numerical quadrature
We proceed to reduce problem (37) by recourse to numerical quadrature. Let
∫
R
g
(
ε
)
d
M
h
(
ε
)
≈
k
∑
j
=0
A
j
g
(
ε
j
)
(38)
be a sequence of quadrature rules, parameterized by
k
∈
N
, with weights
A
j
and nodes
ε
j
.
Here,
A
j
=
∫
R
l
j
(
ε
)
d
M
h
(
ε
)
,
(39)
where
l
j
(
ε
) =
k
∏
i
=0
i
6
=
j
ε
−
ε
i
ε
j
−
ε
i
,
(40)
for
j
= 0
,...,k
are the Lagrange polynomials.
Define the sequence of approximate energies
F
k
(
f
h
) =
k
∑
j
=0
A
j
εf
h
(
ε
j
)
,
(41)
and the sequence of approximate masses
M
k
(
f
h
) =
k
∑
j
=0
A
j
f
h
(
ε
j
)
.
(42)
Then, we have a corresponding sequence of
discretized
problems
inf
f
h
∈
X
h
{
F
k
(
f
h
)
,
0
≤
f
h
≤
1
,
2
M
k
(
f
h
) =
N
}
.
(43)
The solution of these approximate problems then follows from the algorithm:
i) Set
f
0
(
ε
) = 0
,
i
= 0
,...,k
,
I
0
=
{
0
,...,k
}
,
N
0
= 0
,
n
= 1
.
ii) Let
i
n
∈
argmin
{
ε, i
∈
I
n
−
1
}
,
N
n
=
N
n
−
1
+
A
i
n
.
iii) If
N
n
< N
, set
f
n
(
ε
i
n
) = 1
,
I
n
=
I
n
−
1
\{
i
n
}
,
n
←
n
+ 1
, go to (ii).
iv) Otherwise, set
f
n
(
ε
i
n
) = (
N
−
N
n
−
1
)
/A
i
n
,
f
h
=
f
n
, exit.
2.5 Rayleigh-Ritz interpretation
The numerical-quadrature reduction can again be given an appealing Rayleigh-Ritz inter-
pretation. Begin by noting the identity
∫
R
l
i
(
ε
)
εd
M
h
(
ε
) =
k
∑
j
=0
A
j
l
i
(
ε
j
)
ε
j
=
A
i
ε
i
=
(
∫
R
l
i
(
ε
)
d
M
h
(
ε
)
)
ε
i
.
(44)
10
From this identity and (39) we have
F
k
(
f
h
) =
k
∑
i
=0
(
∫
R
l
i
(
ε
)
d
M
h
(
ε
)
)
ε
i
f
(
ε
i
) =
∫
R
(
k
∑
i
=0
l
i
(
ε
)
ε
i
f
(
ε
i
)
)
d
M
h
(
x
)
=
∫
R
(
k
∑
i
=0
l
i
(
ε
)
f
(
ε
i
)
)
εd
M
h
(
ε
) =
F
h
(
f
k
)
,
(45)
where
f
k
=
k
∑
i
=0
l
i
(
ε
)
f
h
(
ε
i
)
.
(46)
Likewise,
M
k
(
f
h
) =
k
∑
i
=0
(
∫
R
l
i
(
ε
)
d
M
h
(
ε
)
)
f
h
(
ε
i
) =
∫
R
(
k
∑
i
=0
l
i
(
ε
)
f
h
(
ε
i
)
)
d
M
h
(
ε
) =
M
h
(
f
k
)
.
(47)
Define now the sequence of spaces
X
k
= span
{
l
i
(
H
)
, i
= 0
,...,k
}
=
{
P
=
∑
k
i
=0
f
i
l
i
(
H
)
, f
i
∈
R
, i
= 0
,...,k
}
, where
{
l
i
, i
= 0
,...,k
}
are the Lagrange polynomials defined by the roots
of the orthogonal polynomial
p
k
+1
generated by
H
. Define, in addition, the sequence of
relaxed constraint sets
K
k
=
{
P
=
∑
k
i
=0
f
i
l
i
(
H
)
,
0
≤
f
i
≤
1
, i
= 0
,...,k
}
. Then, the
reduced problem (43) is equivalent to solving
inf
γ
h
∈X
k
(
I
K
k
(
γ
h
) + 2tr(
h
h
γ
h
)
)
,
(48)
which corresponds to a Rayleigh-Ritz reduction of problem (14) to the subspaces of density
matrices
X
k
generated by numerical quadrature, and to the corresponding relaxed constraint
sets
K
k
.
2.6 Convexification and thermalization
The DFT problem (13) and the inner minimization problem (17) is not convex due to the
constraint that
f
(
ε
)
take values in
{
0
,
1
}
. We convexify the problem by allowing the function
to take values in the entire interval
[0
,
1]
, the resulting function henceforth referred to as
f
β
(
ε
)
. We expect the minimizers to take extreme values only and thus the convexified and
original problems to yield the same minimizers and the same minimum energy. We can now
enforce the constraint
0
≤
f
β
(
ε
)
≤
1
by entropic penalization. We present it for the infinite-
dimensional version (13), though it can readily be extended to the versions with spatial and
spectral discretization. Introduce the entropy
S
(
γ
) = tr[
γ
log(
γ
) + (
I−
γ
) log(
I−
γ
)]
,
(49)
and the thermalized problem
minimize:
U
β
(
γ
) =
U
(
γ
) +
1
β
S
(
γ
) = 2tr(
hγ
) +
2
β
tr[
γ
log(
γ
) + (
I−
γ
) log(
I−
γ
)]
,
(50a)
subject to:
γ
T
=
γ, γh
=
hγ,
2tr(
γ
) =
N,
(50b)
11
where
β
is an inverse temperature. The minimizer of
U
β
(
γ
)
is
γ
β
=
f
β
(
h
) = (
I
+ e
β
(
h
−
μ
I
)
)
−
1
,
(51)
where
μ
is a chemical potential introduced to enforce the number constraint.
μ
and
f
β
are
commonly referred to as the
Fermi level
and
Fermi-Dirac distribution
, respectively, with
β
→∞
representing the zero-temperature limit.
The corresponding minimum value of
U
β
is
−
2
β
tr
(
log(
I
+ e
β
(
h
−
μ
I
)
)
)
=
−
2
β
log det(
I
+ e
β
(
h
−
μ
I
)
)
.
(52)
leading to a thermalized total energy
E
KS
,β
= sup
φ
∈V
sup
v
xc
∈V
[
U
β
(
γ
β
)
−
1
4
π
∫
Ω
1
2
|∇
φ
|
2
d
r
−
E
∗
xc
(
v
xc
)
]
(53)
= sup
φ
∈V
sup
v
xc
∈V
[
2
β
log det(
I
+ e
β
(
h
−
μ
I
)
)
−
1
4
π
∫
Ω
1
2
|∇
φ
|
2
d
r
−
E
∗
xc
(
v
xc
)
]
.
(54)
Finally, we may estimate the zero temperature ground state energy as
E
KS
,
0
≈
sup
φ
∈V
sup
v
xc
∈V
[
U
β
(
γ
β
)
−
1
2
β
S
(
γ
β
)
−
1
4
π
∫
Ω
1
2
|∇
φ
|
2
d
r
−
E
∗
xc
(
v
xc
)
]
.
(55)
2.7 Spatial densities
For later use, we note that the quantities in (13) and (55) have associated spatial densities
and can be rewritten in terms of volume integrals. Following (15), and introducing explicitly
the spatial variables
h
(
r
,
r
′
) =
∫
R
εd
E
r
,
r
′
(
ε
) =
h
T
, γ
β
(
r
,
r
′
) =
∫
R
f
β
(
ε
)
d
E
r
,
r
′
(
ε
) =
γ
T
β
.
(56)
Therefore, the number of electrons, the band structure energy and the entropy may be
written as
N
= 2
tr
(
γ
β
) =
∫
Ω
(
2
∫
R
f
β
(
ε
)
d
E
r
,
r
(
ε
)
)
d
r
=
∫
Ω
ρ
(
r
)
d
r
,
(57)
U
β
= 2
tr
(
hγ
β
) =
∫
Ω
(
2
∫
R
εf
β
(
ε
)
d
E
r
,
r
(
ε
)
)
d
r
=
∫
Ω
u
(
r
)
d
r
,
(58)
S
= 2
tr
(
γ
β
log
γ
β
−
(
I−
γ
β
) log(
I−
γ
β
))
=
∫
Ω
(
2
∫
R
[
f
β
(
ε
) log
f
β
(
ε
) + (1
−
f
β
(
ε
)) log(1
−
f
β
(
ε
))]
d
E
r
,
r
(
ε
)
)
d
r
=
∫
Ω
s
(
r
)
d
r
(59)
in terms of the
charge or number density
ρ
,
band structure energy density
u
and
entropy
density
. Indeed, recall that
ρ
(
r
) =
γ
β
(
r
,
r
)
. These densities play a key role in later sections.
12
2.8 Eigenvalue problem
We close the formulation by connecting our formulation to the way DFT is usually presented
as an eigenvalue problem. The direct solution of problem (30) entails the computation of
N/
2
eigenvalues and eigenvectors. To see this, consider the inner
inf
operation in (30). Write
P
=
Q
h
Q
T
h
,
(60)
Q
h
∈L
(
H
h
,
R
N/
2
)
, with
Q
T
h
Q
h
=
I
N/
2
,
(61)
where
I
N/
2
denotes identity in
L
(
R
N/
2
)
. Here and subsequently,
L
(
A
,
B
)
denotes the space
of linear transformations between two linear spaces
A
and
B
, and
L
(
A
)
the space of linear
transformations from a linear space
A
to itself. Then,
P
T
=
P
,
0
≤
P
≤
1
and
tr(
P
) =
N/
2
,
hence
γ
h
∈K
h
. The problem under consideration thus becomes
Q
h
∈
argmin
{
2tr(
Q
T
h
HQ
h
)
, Q
T
h
Q
h
=
I
N/
2
}
,
(62)
The Euler-Lagrange equations of this problem are
HQ
h
=
Q
h
Λ
h
,
(63)
where
Λ
h
∈ L
(
R
N/
2
)
,
Λ
T
h
= Λ
h
, is a Lagrange multiplier. Clearly, these Euler-Lagrange
equations are solved if the columns of
Q
h
consist of eigenvectors of
H
and
Λ
h
stores the
corresponding eigenvalues in its diagonal. In addition, if
{
ε
1
,...,ε
N
g
}
are the ordered eigen-
values of
H
is ascending order and
{
φ
1
,...,φ
N
g
}
are the corresponding eigenvectors, then the
minimum problem is solved by
Q
h
=
{
φ
1
,...,φ
N/
2
}
and
Λ
h
= diag
{
ε
1
,...,ε
N/
2
}
. Finally,
the energy follows as
2tr(
Q
T
h
HQ
h
) = 2tr(
Q
T
h
Q
h
Λ
h
) = 2tr(Λ
h
)
.
(64)
Clearly, this computation becomes intractable for large material samples containing a large
number electrons
N
. Therefore, computational tractability of large samples requires an
additional reduction (beyond spatial discretization) that we refer to as
spectral reduction
above.
3 Filtering, spectrum splitting and pseudopotentials
This section introduces three ideas that enable faster calculations. The first two, filtering and
spectrum splitting, are convergent approaches and take advantage of the spectral formulation.
The third, pseudopotentials, involves modeling.
3.1 Filtering
The discrete DFT problem (30) is posed as a problem in
N
g
-dimensional subspace
V
h
of
V
=
H
1
0
(Ω)
. In practice, the accurate solution of the equations requires that
N
g
>> N
.
However, the solution to our problem, the density matrix
γ
h
, has rank
N
(in the thermalized
problem, the thermalized density matrix
γ
β,h
has rank larger than but close to
N
). Therefore,
13
one can obtain significant savings in computational effort if one could identify
a priori
a sub-
space
V
f
h
such that range
(
γ
h
)
⊂V
f
h
⊂⊂V
h
, and restrict the problem (31) and specifically the
Hamiltonian
h
h
to the sub-subspace
V
f
h
. This can be achieved using filtering. While many
approaches have been proposed based on filtering such as purification (cf. e.g. [47, 33, 58, 74])
and approximations to the Fermi-Dirac functions (cf. e.g. [26, 4, 48]), the Chebyshev filtering
technique [95, 96] is being adopted in many recent DFT codes [24, 23, 50, 52]. The main
idea in Chebyshev filtering is to approximate the subspace
V
f
h
as
V
f
h
≈
T
m
(
g
(
h
h
))
X
h
(65)
where
X
h
⊂V
h
with dim
(
X
h
) =
dim
(
V
f
h
)
,
T
m
is a Chebyshev polynomial of order
m
, and
g
(
x
) =
2
b
−
a
(
x
−
b
+
a
2
)
, b > a
(66)
with
b
= max
σ
(
h
h
)
>>
1
(
σ
denoting the spectrum) and
a
= max
σ
(
γ
h
h
h
) +
O
(1)
. In
particular,
a
≈
μ
+
O
(1)
is a reasonable choice. We note that
g
transforms the spectrum of
h
h
such that
σ
(
g
(
h
h
))
<
1
and
σ
(
g
(
γ
h
h
h
)))
∈
(
−∞
,
−
1)
. Thus, as
T
m
(
x
)
>
1
for
x
∈
(
−∞
,
−
1)
,
T
m
(
g
(
h
h
))
X
h
provides a good approximation to
V
f
h
. We note that the suitable choice of
m
depends on the
(
b
−
a
)
, with a larger
m
that would be needed for larger values of
(
b
−
a
)
.
For instance, based on numerical studies, if
(
b
−
a
)
∼ O
(10
2
)
, values of
m
∼
10
−
30
are
sufficient to construct a good approximation to
V
f
h
[96, 52]. However, if
(
b
−
a
)
∼ O
(10
6
)
,
values of
m
∼
1
,
000
are needed [69].
If
P
f
:
V
h
→ V
f
h
denotes the projection operator onto the filtered subspace, then the
solution to the DFT problem can be obtained by replacing
h
h
in (31) with
P
f
h
h
P
f
. As
the spectral width
Σ(
P
f
h
h
P
f
)
<<
Σ(
h
h
)
, it enables faster numerical solution of the DFT
problem, and has been the basis for subspace projection methods (cf. e.g. [18, 53]).
3.2 Spectrum splitting
The next idea combines filtering with a feature of the solution of typical problem. Here, we
assume that the DFT problem has already been projected onto
V
f
h
, and denote
h
f
=
P
f
h
h
P
f
.
We denote
σ
h
=
σ
(
h
f
)
as the spectrum of
h
f
, and assume in the following that
max
σ
h
<
0
(i.e.,
h
f
is appropriately shifted such that this condition is satisfied). It has long been
recognized that the spectrum of
h
f
has a gap that separates the so-called
core
, or deeply
bound states at the lower end, from the rest. In other words, the spectrum
σ
h
=
σ
c
h
∪
σ
r
h
with
ε
′
+
E
g
≤
ε
′′
∀
ε
′
∈
σ
c
h
,ε
′′
∈
σ
r
h
for a gap
E
g
>
0
. We can therefore split the Hamiltonian
h
f
and the density operator
γ
h
(corresponding to
h
f
) into
h
f
=
h
c
f
+
h
r
f
, γ
h
=
γ
c
h
+
γ
r
h
,
(67)
where the spectrum of
h
c
f
is
σ
c
h
∪{
0
}
. It follows that we can divide
V
f
h
into two orthogonal
subspaces,
V
f
h
=
V
c
h
⊕V
r
h
(68)
where
V
c,r
h
is the range of
h
c,r
f
. Further, since
σ
c
h
is the lower end of the spectrum, it follows
that
γ
c
h
=
P
c
h
(69)
14
is the projection operator from
V
f
h
to
V
c
h
.
Now, in light of the spectral gap, we can again use filtering on
h
f
, and then readily
identify
V
c
h
as the range of
h
c
f
. Therefore, we can use (69) to easily compute
γ
c
h
. Further,
using the orthogonality of the subspaces,
h
r
f
= (
I−P
c
h
)
h
f
(
I−P
c
h
)
.
(70)
Since the spectra
σ
c
h
and
σ
r
h
are disjoint, it follows
tr(
h
f
γ
h
) = tr(
h
c
f
γ
c
h
) + tr(
h
r
f
γ
r
h
)
.
(71)
We may now reduce (31) as
minimize:
U
r
h
(
γ
r
h
) := 2tr(
h
r
f
γ
r
h
)
,
(72a)
subject to:
(
γ
r
h
)
T
=
γ
r
h
, γ
r
h
h
r
f
=
h
r
f
γ
r
h
,
(
γ
r
h
)
2
=
γ
r
h
,
2tr(
γ
r
h
) =
N
−
N
c
,
(72b)
where
N
c
denotes the number of core electrons. This approach of spectrum splitting provides
a number of advantages. First, the computation of
γ
c
h
, the core part of the density matrix,
is relatively simple as described above. Second, in practice, the width of spectrum of
h
r
f
(
Σ(
h
r
f
)
) is significantly smaller than that of
h
f
(
Σ(
h
f
)
), and this allows for a more efficient
numerical solution. Finally, the core subspace
V
c
h
consists of functions which have a compact
support close to the nuclei. In other words, this is the subspace spanned by the orbitals
of the core electrons. This can be further exploited to gain numerical efficiency. Further,
its complement,
V
r
h
, that contains so-called valance and conduction electrons, consists of
functions that vary smoothly outside a core region around the nucleus. Therefore, we can
use a spatially adaptive resolution to discretize it.
We may proceed similarly in the thermalized problem to find that (69) and (71) still
holds, and
γ
r
h,β
=
f
β
(
h
r
f
)
.
(73)
It is common to compute this by expanding this in a polynomial basis (Fermi operator ex-
pansion [26, 27]), which we shall show later in Section 5 is related to the spectral quadratures.
Therefore, the advantages of spectrum splitting carry over to the thermalized setting.
The accuracy and efficacy of this approach for large-scale all-electron DFT calculations
has been demonstrated in [55]. Here, we present some representative results on Si and Au
nanoclusters. Figure 1 shows the results from ground-state energies computed using two ap-
proaches: (i) SubPJ-FE: A subspace projection approach via filtering (Sec 3.1) implemented
in finite-element basis, where
γ
h,β
=
f
β
(
h
f
)
is computed via Fermi-operator expansion using
Chebyshev polynomials for various orders; (ii) Spectrum-splitting method: In addition to
the subspace projection via filtering, spectrum splitting is used, where
γ
h,β
=
γ
c
h
+
f
β
(
h
r
f
)
and
f
β
(
h
r
f
)
is evaluated via Fermi-operator expansion using Chebyshev polynomials for var-
ious orders. The results for Si
95
are provided for two values of
β
corresponding to
T
= 500
and
1000
K, and results for Au
6
cluster are shown for
T
= 500
K. As is evident, spectrum
splitting not only provides computationally efficiency—due to a substantial reduction in the
polynomial order required in Fermi operator expansion—it is indispensable to obtain the
desired accuracy for systems with large atomic numbers, like Au.
15
0
500
1000
1500
2000
10
−
8
10
−
6
10
−
4
10
−
2
10
0
Fermi
−
operator Expansion Degree
Abs. Error in Energy (Ha/atom)
Spectrum splitting method: T=500 K
SubPJ
−
FE: T=500 K
Spectrum splitting method: T=1000 K
SubPJ
−
FE: T=1000 K
0
500
1000
1500
2000
10
−
6
10
−
4
10
−
2
10
0
10
2
Fermi
−
operator Expansion Degree
Abs. Error in Energy (Ha/atom)
Spectrum splitting method
SubPJ
−
FE
Figure 1:
Accuracy and computational efficiency afforded by spectrum splitting in all-electron
calculations via Fermi operator expansion.
(Left)
Si
95
nanocluster;
(Right)
Au
6
nanocluster.
Adapted from [55].
We conclude this subsection by noting that spectrum splitting is also closely related to
the so-called
enrichment
methods. Note that the identity (69) means that we can use any
basis set to represent
V
c
h
. Therefore, picking
∼
N
c
/
2
functions that are computationally
convenient and approximate the span of
V
c
h
provides a good starting point. Subsequently,
choosing a spatial discretization sufficient to span
V
f
h
provides the desired accuracy. This
is computationally effective since the spatial discretization does not have to be so fine as
to represent the core electrons. The basis set approximately spanning
V
c
h
can be iteratively
updated as the calculation proceeds. These ideas lead to augmented plane wave (APW),
linearized augmented plane wave (LAPW) [71], and enriched finite basis [92, 41, 69]. We
refer the reader to the chapter by Chen and Schneider [13] for a detailed discussion of these
methods.
3.3 Frozen core approximation and pseudopotentials
The formulations discussed till now have consider all electrons in the system. However, it is
a long-held observation in the field that core electrons play a minimal role in the bonding
between atoms. Specifically, it is observed that the
γ
c
h
is relatively independent of the external
potentials
v
that arise in molecules and crystals. This motivates the desire to exclude these
electrons from the calculations, and to focus on the valance and conduction electrons.
One approach to doing so is the so-called
frozen core approximation
. Here, a high resolu-
tion all-electron calculation for a single atom is conducted to obtain the core density matrix,
γ
c
h,Z
for a single atom (the subscript
Z
here refers to the single atom of atomic number
Z
with the nucleus located at the origin). Subsequently, this is used as an ansatz for the
core electrons for any given problem. Specifically, for a problem with
N
a
atoms with atomic
numbers
{
Z
i
}
located at
{
r
i
}
,
̄
γ
c
h
(
r
,
r
′
) =
N
a
∑
i
=1
γ
c
h,Z
i
(
r
−
r
i
,
r
′
−
r
i
)
(74)
16
Table 1:
Error incurred from the frozen core approximation for various systems [51].
System
||
ρ
c
0
−
̄
ρ
c
||
L
2
/
||
ρ
c
0
||
L
2
||
ρ
0
−
̄
ρ
0
||
L
2
/
||
ρ
0
||
L
2
Li
2
0.00703
0.00787
O
2
0.00102
0.00128
CO
0.00181
0.00129
Si
18
0.01272
0.0130
Si
31
0.01273
0.0134
is used as an ansatz for
γ
h
= ̄
γ
c
h
+
γ
r
h
(75)
in (31) to solve for
γ
r
h
. Note that the computational complexity of the problem is now reduced
from
N
electrons to
N
−
N
c
electrons. Further, as noted above, the range of
γ
r
h
is spanned by
relatively smooth functions outside the core, and therefore one can use a spatially adaptive
discretization to represent this problem.
Note that this is an uncontrolled approximation since it is based on an ansatz. Table 1
from Ref. [51] shows the errors from the frozen core approximation for a range of systems.
In particular, the two metrics used to measure the approximation are: (i) the relative error
in the core electron density at the ground-state
||
ρ
c
0
−
̄
ρ
c
||
L
2
/
||
ρ
c
0
||
L
2
, where
ρ
c
0
is the core
electron density at the ground-state from the all-electron calculation and
̄
ρ
c
= ̄
γ
c
h
(
r
,
r
)
; (ii)
the relative error in the total electron density at the ground-state
||
ρ
0
−
̄
ρ
0
||
L
2
/
||
ρ
0
||
L
2
, where
ρ
0
is the total electron density at the ground-state from the all-electron calculation, and
̄
ρ
0
is the total electron density at the ground-state from the frozen core approximation. As
evident, while the approximation is good for some systems, it can incur larger errors for
others (such as Si nanoclusters).
A closely related idea is that of a
pseudopotential
. Here, the objective is to fully exclude
the core states by using a fictitious potentials, namely pseudopotentials, thus replacing
h
h
with
h
PS
. The pseudopotentials are generated such that
γ
r
PS
=
f
β
(
h
PS
)
closely approximates
γ
r
h
outside a core radius around each atom, but the range of
γ
r
PS
is smooth all through the
simulation domain. Thus, this alleviates the need for a spatially refined basis to resolve
the core states. Various pseudopotentials have been proposed and are widely used (cf.
e.g. [86, 8, 31]). Despite the errors and the uncontrolled nature of these approximations, it
is often the only practical route to proceed in large systems of interest.
4 Spatial coarse-graining: Finite-element discretization
Spatial discretization (cf. Sec 2.2) plays a central role in the practical aspects of computing
the solution to the Kohn-Sham DFT problem in an efficient manner. Many discretization
schemes have been adopted by the scientific community in solving the Kohn-Sham problem,
and besides the algorithms employed, the discretization schemes have been the main differ-
entiator for the various DFT codes and their performance based on computational efficiency
and scalability. The widely used discretization methods include the plane-wave basis (cf.
17
e.g. [45, 29, 25]) and atomic orbital type basis functions (cf. e.g. [34, 38, 9, 84]). While
the plane-wave basis offers spectral convergence, it is primarily efficient for periodic problem
owing to lack of spatial adaptivity, and is constrained by limited parallel scalability of nu-
merical implementations. The atomic orbital type basis functions present a reduced order
basis, but in practice may not guarantee a robust and systematically convergent solution,
especially for metallic systems. Also, they suffer from limited parallel scalability owing to
the global nature of the basis functions. The finite-element and finite difference discretiza-
tion schemes, while have been explored over two decades ago [81, 82, 61, 62, 46], are only
recently gaining traction as efficient and scalable approaches for solving the Kohn-Sham
problem [56, 52, 24, 23].
The finite-element discretization in particular offers many attractive features including
the following: (i) Systematic convergence. Piecewise polynomials of a fixed degree
p
are dense
in
H
1
(Ω)
as the finite-element mesh-size
h
becomes small. Further, polynomials of increasing
p
are dense for a fixed
h
. (ii) Flexibility. Ability to easily handle complex geometries and
mixed boundary conditions that is especially important to treat defects where periodicity
may not be appropriate. (iii) Spatial adaptivity. The discretization can be exploited to
provide desired basis resolution in regions of interest and coarse-graining elsewhere. (iv)
Parallel scalability. The locality of the FE basis provides for efficient parallel scalability of
numerical implementation. We also refer the reader to the chapter by Dai and Zhou [14] for
a broad discussion of the application of finite element discretization to DFT.
4.1 Higher-order spectral finite-elements
Despite the aforementioned advantages of the finite-element basis, and many prior efforts
that explored the use of finite-element basis for electronic structure calculations, they have
not been competitive with widely used plane-wave and atomic orbital basis sets until re-
cently. The two main issues limiting the performance of finite-element basis in Kohn-Sham
DFT had been: (i) the significant degree of freedom disadvantage of commonly used linear
finite-elements in comparison to plane-wave basis that affects the computational efficiency in
practical DFT calculations; (ii) the non-orthogonality of the finite-element basis that either
limits the available solution schemes or requires an additional evaluation of the inverse of
the overlap matrix.
Figure 2 provides insights into the lack of computational efficiency of linear finite-elements
observed in prior studies. The figure shows the error in the ground-state energy for various
finite element discretizations of different finite-element orders for two materials systems The
higher order finite-elements employed in the study are hexahedral finite-elements, where the
finite-element basis functions are constructed as a tensor product of basis functions in each
dimension. The hexahedral finite-element basis functions in the isoparametric formulation
are constructed from polynomial basis functions in the reference domain
[
−
1
,
1]
3
as
P
i,j,k
(
ξ,η,κ
) =
l
i
(
ξ
)
l
j
(
η
)
l
k
(
κ
)
, l
i
(
ξ
) :=
∏
0
≤
m<p
m
6
=
i
ξ
−
ξ
m
ξ
i
−
ξ
m
, i,j,k
= 0
,
1
,...p
(76)
where
l
i
(
ξ
)
is a Lagrange polynomial of degree
p
constructed based on the
p
+ 1
nodes
of the finite-element. Conventionally, the finite-element nodes are chosen to be equidistant,
18