NUMERICAL STABILITY AND EFFICIENCY OF RESPONSE PROPERTY CALCULATIONS IN DENSITY FUNCTIONAL THEORY ERIC CANCÈS12 MICHAEL F. HERBST3 GASPARD KEMLIN12 ANTOINE LEVITT12 AND BENJAMIN STAMM4

2025-05-02 0 0 1.29MB 21 页 10玖币

侵权投诉

NUMERICAL STABILITY AND EFFICIENCY OF RESPONSE PROPERTY

CALCULATIONS IN DENSITY FUNCTIONAL THEORY

ERIC CANCÈS1,2, MICHAEL F. HERBST3, GASPARD KEMLIN1,2, ANTOINE LEVITT1,2, AND BENJAMIN STAMM4

Abstract. Response calculations in density functional theory aim at computing the change in ground-state density

induced by an external perturbation. At ﬁnite temperature these are usually performed by computing variations

of orbitals, which involve the iterative solution of potentially badly-conditioned linear systems, the Sternheimer

equations. Since many sets of variations of orbitals yield the same variation of density matrix this involves a choice

of gauge. Taking a numerical analysis point of view we present the various gauge choices proposed in the literature

in a common framework and study their stability. Beyond existing methods we propose a new approach, based

on a Schur complement using extra orbitals from the self-consistent-ﬁeld calculations, to improve the stability and

eﬃciency of the iterative solution of Sternheimer equations. We show the success of this strategy on nontrivial

examples of practical interest, such as Heusler transition metal alloy compounds, where savings of around 40% in

the number of required cost-determining Hamiltonian applications have been achieved.

1. Introduction

Kohn-Sham (KS) density-functional theory (DFT) [26,32] is the most popular approximation to the electronic

many-body problem in quantum chemistry and materials science. While not perfect, it oﬀers a favourable com-

promise between accuracy and computational eﬃciency for a vast majority of molecular systems and materials. In

this work, we focus on KS-DFT approaches aiming at computing electronic ground-state (GS) properties. Having

solved the minimization problem underlying DFT directly yields the ground-state density and corresponding energy.

However, many quantities of interest, such as interatomic forces, (hyper)polarizabilities, magnetic susceptibilities,

phonons spectra, or transport coeﬃcients, correspond physically to the response of GS quantities to a change in

external parameters (e.g. nuclear positions, electromagnetic ﬁelds). As such their mathematical expressions in-

volve derivatives of the obtained GS solution with respect to these parameters. For example interatomic forces are

ﬁrst-order derivatives of the GS energy with respect to the atomic positions, and can actually be obtained without

computing the derivatives of the GS density, thanks to the Hellmann-Feynman theorem [21]. On the other hand the

computation of any property corresponding to second- or higher-order derivatives of the GS energy does require the

computation of derivatives of the density. More precisely, it follows from Wigner’s (2n+ 1) theorem that nth-order

derivatives of the GS density are required to compute properties corresponding to (2n)th- and (2n+ 1)st-derivatives

of the KS energy functional. More recent applications, such as the design of machine-learned exchange-correlation

energy functionals, also require the computation of derivatives of the ground-state with respect to parameters, such

as the ones deﬁning the exchange-correlation functional [27,29,34].

Eﬃcient numerical methods for evaluating these derivatives are therefore needed. The application of generic

perturbation theory to the special case of DFT is known as density-functional perturbation theory (DFPT) [2,

15,16,18]. See also [37] for applications to quantum chemistry, [1] for applications to phonon calculations, and

[9] for a mathematical analysis of DFPT within the reduced Hartree-Fock (rHF) approximation (also called the

Hartree approximation in the physics literature). Although the practical implementation of ﬁrst- and higher-order

derivatives computed by DFPT in electronic structure calculation software can be greatly simpliﬁed by Automatic

Diﬀerentiation techniques [19], the eﬃciency of the resulting code crucially depends on the eﬃciency of a key

building block: the computation of the linear response δρ of the GS density to an inﬁnitesimal variation δV of the

total Kohn-Sham potential.

(1) CERMICS, ENPC, France

(2) MATHERIALS team, Inria Paris, France

(3) Applied and Computational Mathematics, RWTH Aachen University, Germany

(4) IANS, University of Stuttgart, Germany

E-mail addresses:eric.cances@enpc.fr, herbst@acom.rwth-aachen.de, gaspard.kemlin@enpc.fr, antoine.levitt@inria.fr, best@ians.uni-stuttgart.de.

arXiv:2210.04512v2 [math.NA] 16 Feb 2023

2 NUMERICAL STABILITY AND EFFICIENCY OF RESPONSE PROPERTY CALCULATIONS IN DFT

×× ×××

ε1×

εNp

εF× ××× ××××

×× ××× ×+

εF× ××× ××××

fε−εF

T

Figure 1 – The occupation numbers fnfor T= 0 (left) and T > 0(right).

For reasons that will be detailed below, the numerical evaluation of the linear map δV 7→ δρ is not straight-

forward, especially for periodic metallic systems. Indeed, DFT calculations for metallic systems usually require

the introduction of a smearing temperature T, a numerical parameter which has nothing to do with the physical

temperature (in practice, its value is often higher than the melting temperature of the crystal). For the sake of

simplicity, let us ﬁrst consider the case of a periodic simulation cell Ωcontaining an even number Nel of electrons in

a spin-unpolarized state (see Remark 1 for details on how this formalism allows for the computation of properties of

perfect crystals). The Kohn-Sham GS at ﬁnite temperature T > 0is then described by an L2(Ω)-orthonormal set of

orbitals (φn)n∈N∗with energies (εn)n∈N∗, which are the eigenmodes of the Kohn-Sham Hamiltonian Hassociated

with the GS density:

Hφn=εnφn,ˆΩ

φ∗

m(r)φn(r)dr=δmn, ε16ε26ε36· · · ,

together with periodic boundary conditions. The GS density in turn reads

ρ(r) =

+∞

n=1

fn|φn(r)|2with fn:=fεn−εF

T,(1)

where fis a smooth occupation function converging to 2at −∞ and to zero at +∞,e.g. the Fermi-Dirac smearing

function f(x) = 2

1+ex(see Figure 1). The Fermi level εFis the Lagrange multiplier of the neutrality charge

constraint: it is the unique real number such that

ˆΩ

ρ(r)dr=

+∞

n=1

fn=

+∞

n=1

fεn−εF

T=Nel.

It follows from perturbation theory that the linear response δρ of the density to an inﬁnitesimal variation δV of the

total Kohn-Sham potential is given by

δρ =χ0δV,

where χ0is the independent-particle susceptibility operator (also called noninteracting density response function).

Equivalently, this operator describes the linear response of a system of noninteracting electrons of density ρsubject

to an inﬁnitesimal perturbation δV . It holds (see Section 3)

δρ(r):= (χ0δV )(r) =

+∞

n=1

+∞

m=1

fn−fm

εn−εm

φ∗

n(r)φm(r)(δVmn −δεFδmn),(2)

where δVmn :=hφm, δV φni,δεFis the induced variation of the Fermi level εF, and δmn is the Kronecker symbol.

We also use the convention (fn−fn)/(εn−εn) = 1

Tf0εn−εF

T.

In practice, these equations are discretized on a ﬁnite basis set, so that the sums in (1) and (2) become ﬁnite.

Since the number of basis functions Nbis often very large compared to the number of electrons in the system, it is

very expensive to compute the sums as such. However, in practice it is possible to restrict to the computation of a

number NNbof orbitals. These orbitals are then computed using eﬃcient iterative methods [38].

For insulating systems, there is a (possibly) large band gap between εNpand εNp+1 which remains non-zero in

the thermodynamic limit of a growing simulation cell. As a result, the calculation can be done at zero temperature,

such that the occupation function fbecomes a step function (see Figure 1). The jump from 2to 0in the occupations

occurs exactly when the lowest Np=Nel/2energy levels ε16· · · 6εNpare occupied with an electron pair (two

electrons of opposite spin). Thus, fn= 2 for 16n6Npand fn= 0 for n>Np. As a result, Ncan be chosen equal

to the number of electron pairs Npwithout any approximation. In contrast, for metallic systems εNp=εNp+1 =εF

in the zero-temperature thermodynamic limit (more precisely there is a positive density of states at the Fermi level

in the limit of an inﬁnite simulation cell), causing the denominators in the right-hand side of formula (2) to formally

NUMERICAL STABILITY AND EFFICIENCY OF RESPONSE PROPERTY CALCULATIONS IN DFT 3

blow up. Calculations on metallic systems are thus done at ﬁnite temperature T > 0, in which case every orbital

has a fractional occupancy fn∈(0,2). However, since from a classical semiclassical approximation εntends to

inﬁnity as n2/3as n→ ∞, and fdecays very quickly, one can safely assume that only a ﬁnite number Nof orbitals

have nonnegligible occupancies. This allows one to avoid computing φnfor n > N . Under this approximation, a

formal diﬀerentiation of (1) gives

δρ(r) =

n=1

fn(φ∗

n(r)δφn(r) + δφ∗

n(r)φn(r)) + δfn|φn(r)|2.(3)

However, while the response δρ to a given δV is well-deﬁned by (2), the set (δφn, δfn)16n6Nis not. Indeed, the Kohn-

Sham energy functional being in fact a function of the density matrix γ=PN

n=1 fn|φnihφn|, any transformation of

(δφn, δfn)16n6Nleaving invariant the ﬁrst-order variation

δγ :=

n=1

δfn|φnihφn|+

n=1

fn(|φnihδφn|+|δφnihφn|)(4)

of the density matrix is admissible. This gauge freedom can be used to stabilize linear response calculations or,

in the contrary, may lead to numerical instabilities. Denote by Pthe orthogonal projector onto Span(φn)16n6N,

the space spanned by the orbitals considered as (partially) occupied, and by Q= 1 −Pthe orthogonal projector

onto the space Span(φn)n>N spanned by the orbitals considered as unoccupied. Then, the linear response of any

occupied orbital can be decomposed as δφn=δφP

n+δφQ

nwhere:

•δφP

n=P δφn∈Ran(P)can be directly computed via a sum-over-state formula (explicit decomposition

on the basis of (φn)n6N). This contribution can be chosen to vanish in the zero-temperature limit, as in

that case P δγ P = 0. At ﬁnite temperature, a gauge choice has to be made and several options have been

proposed in the literature;

•δφQ

n=Qδφn∈Ran(Q)is the unique solution of the so-called Sternheimer equation [44]

Q(H−εn)QδφQ

n=−QδV φn,(5)

where His the Kohn-Sham Hamiltonian of the system. This equation is possibly very ill-conditioned for

n=Nif εN+1 −εNis very small.

This paper addresses these two issues. First, we review and analyse the diﬀerent gauge choices for δφP

nproposed

in the literature and introduce a new one. We bring all these various gauge choices together in a new common

framework and analyse their performance in terms of numerical stability. Second, for the contribution δφQ

n, we

investigate how to improve the conditioning of the linear system (5), which is usually solved with iterative solvers

and we propose a new approach. This new approach is based on the fact that, as a byproduct of the iterative

computation of the ground state orbitals (φn)n6N, one usually obtains relatively good approximations of the

following eigenvectors. This information is often discarded for response calculations; we use them in a Schur

complement approach to improve the conditioning of the iterative solve of the Sternheimer equation. We quantify

the improvement of the conditioning obtained by this new approach and illustrate its eﬃciency on several metallic

systems, from aluminium to transition metal alloys. We observe a reduction of typically 40% of the number of

Hamiltonian applications (the most costly step of the calculation). The numerical tests have been performed

with the DFTK software [25], a recently developed plane-wave DFT package in Julia allowing for both easy

implementation of novel algorithms and numerical simulations of challenging systems of physical interest. The

improvements suggested in this work are now the default choice in DFTK to solve response problems.

This paper is organized as follows. In Section 2, we review the periodic KS-DFT equations and the associated

approximations. We also present the mathematical formulation of DFPT and we detail the links between the

orbitals’ response δφnand the ground-state density response δρ for a given external perturbation, as well as the

derivation of the Sternheimer equation (5). In Section 3, we propose a common framework for diﬀerent natural

gauge choices. Then, with focus on the Sternheimer equation and the Schur complement, we present the improved

resolution to obtain δφQ

n. Finally, in Section 4, we perform numerical simulations on relevant physical systems. In

the appendix, we propose a strategy for choosing the number of extra orbitals motivated by a rough convergence

analysis of the Sternheimer equation.

4 NUMERICAL STABILITY AND EFFICIENCY OF RESPONSE PROPERTY CALCULATIONS IN DFT

2. Mathematical framework

2.1. Periodic Kohn-Sham equations. We consider here a simulation cell Ω = [0,1)a1+ [0,1)a2+ [0,1)a3

with periodic boundary conditions, where (a1,a2,a3)is a nonnecessarily orthonormal basis of R3. We denote by

R=Za1+Za2+Za3the periodic lattice in the position space and by R∗=Zb1+Zb2+Zb3with ai·bj= 2πδij

the reciprocal lattice. Let us denote by

#(R3,C):={u∈L2

loc(R3,C)|uis R-periodic}(6)

the Hilbert space of complex-valued R-periodic locally square integrable functions on R3, endowed with its usual

inner product h·,·i and by Hs

#(R3,C)the R-periodic Sobolev space of order s∈R

#(R3,C):=(u=X

G∈R∗buGeG,X

G∈R∗

(1 + |G|2)s|buG|2<∞)

where eG(r)=eiG·r/p|Ω|is the Fourier mode with wave-vector G.

In atomic units, the KS equations for a system of Nel = 2Npspin-unpolarized electrons at ﬁnite temperature

read

Hρφn=εnφn, ε16ε26· · · ,hφn, φmi=δnm, ρ(r) =

+∞

n=1

fn|φn(r)|2,

+∞

n=1

fn=Nel,(7)

where Hρis the Kohn-Sham Hamiltonian. It is given by

Hρ=−1

2∆ + V+VHxc

ρ(8)

where Vis the potential generated by the nuclei (or the ionic cores if pseudopotentials are used) of the system, and

VHxc

ρ(r) = VH

ρ(r) + Vxc

ρ(r)is an R-periodic real-valued function depending on ρ. The Hartree potential VH

ρis the

unique zero-mean solution to the periodic Poisson equation −∆VH

ρ(r)=4πρ(r)−1

|Ω|´Ωρand the function Vxc

is the exchange-correlation potential. Hρis a self-adjoint operator on L2

#(R3,C), bounded below and with compact

resolvent. Its spectrum is therefore composed of a nondecreasing sequence of eigenvalues (εn)n∈N∗converging to

+∞. Since Hρdepends on the electronic density ρ, which in turn depends on the eigenfunctions φn, (7) is a

nonlinear eigenproblem, usually solved with self-consistent ﬁeld (SCF) algorithms. These algorithms are based on

successive partial diagonalizations of the Hamiltonian Hρnbuilt from the current iterate ρn. See [6,7,35] and

references therein for a mathematical presentation of SCF algorithms.

In (7), the φn’s are the Kohn-Sham orbitals, with energy εnand occupation number fn. At ﬁnite temperature

T > 0,fnis a real number in the interval [0,2) and we have

fn=fεn−εF

T,(9)

where fis a ﬁxed analytic smearing function, which we choose here equal to twice the Fermi-Dirac function:

f(x)=2/(1 + ex). The Fermi level εFis then uniquely deﬁned by the charge constraint

+∞

n=1

fn=Nel.(10)

When T→0,f((· − εF)/T )→2×1{·<εF}almost everywhere, and only the lowest Np=Nel/2energy levels for

which εn< εFare occupied by two electrons of opposite spins (see Figure 1): fn= 2 for n6Npand fn= 0 for

n>Np.

Remark 1 (The case of perfect crystals).Using a ﬁnite simulation cell Ωwith periodic boundary conditions is

usually the best way to compute the bulk properties of a material in the condensed phase. Indeed, KS-DFT simu-

lations are limited to, say 103−104electrons, on currently available standard computer architectures. Simulating

in vacuo a small sample of the material containing, say 103atoms, would lead to completely wrong results, polluted

by surface eﬀects since about half of the atoms would lay on the sample surface. Periodic boundary conditions

are a trick to get rid of surface eﬀects, at the price of artiﬁcial interactions between the sample and its periodic

image. In the case of a perfect crystal with Bravais lattice Land unit cell ω, it is natural to choose a periodic

simulation (super)cell Ω = Lω consisting of L3unit cells (we then have R=LL). In the absence of spontaneous

symmetry breaking, the KS ground-state density has the same L-translational invariance as the nuclear potential.

NUMERICAL STABILITY AND EFFICIENCY OF RESPONSE PROPERTY CALCULATIONS IN DFT 5

Using Bloch theory, the supercell eigenstates φncan then be relabelled as φn(r) = eik·rujk(r), where ujknow has

cell periodicity, and equations (7)–(9) can be rewritten as

Hρ,kujk=εjkujk, ε1k6ε2k6· · · ,hujk, uj0ki=δjj0,(11)

ρ(r) = 1

L3X

k∈GL

+∞

j=1

fjk|ujk(r)|2,1

L3X

k∈GL

+∞

j=1

fjk=Nel, fjk=fεjk−εF

T(12)

Hρ,k=1

2(−i∇+k)2+V+VHxc

ρ,(13)

where GL=L−1L∗∩ω∗. Here L∗is the dual lattice of Land ω=R3/L∗the ﬁrst Brillouin zone of the crystal. In

the thermodynamic limit L→ ∞, we obtain the periodic Kohn-Sham equations at ﬁnite temperature

Hρ,kujk=εjkujk, ε1k6ε2k6· · · ,hujk, uj0ki=δjj0,(14)

ρ(r) = ω∗

+∞

j=1

fjk|ujk(r)|2dk, ω∗

+∞

j=1

fjkdk=Nel, fjk=fεjk−εF

T.(15)

This is a massive reduction in complexity, as now only computations on the unit cell have to be performed. For

metals, the integrand on the Brillouin zone is discontinuous at zero temperature, which makes standard quadrature

methods fail. Introducing a smearing temperature T > 0allows one to smooth out the integrand, see [5,33] for a

numerical analysis of the smearing technique. We also refer for instance to [40, Section XIII.16] for more details

on Bloch theory, to [4,10] for a proof of the thermodynamic limit for perfect crystals in the rHF setting for both

insulators and metals, and to [14] for the numerical analysis for insulators.

2.2. Density-functional perturbation theory. We detail in this section the mathematical framework of DFPT.

We ﬁrst rewrite the Kohn-Sham equations (7) as the ﬁxed-point equation for the density ρ

FV+VHxc

ρ=ρ, (16)

where Fis the potential-to-density mapping deﬁned by

F(V)(r) =

+∞

n=1

fεn−εF

T|φn(r)|2(17)

with (εn, φn)n∈N∗an orthonormal basis of eigenmodes of −1

2∆ + Vand εFdeﬁned by (10). The solution of (16)

deﬁnes a mapping from Vto ρ: the purpose of DFPT is to compute its derivative. Let δV0be a local inﬁnitesimal

perturbation, in the sense that it can be represented by a multiplication operator by a periodic function r7→ δV0(r).

By taking the derivative of (16) with the chain rule, we obtain the implicit equation for δρ:

δρ =F0V+VHxc

ρ·(δV0+KHxc

ρδρ),(18)

where the Hartree-exchange-correlation kernel KHxc

ρis the derivative of the map ρ7→ VHxc

ρand F0V+VHxc

ρis

the derivative of Fcomputed at V+VHxc

ρ. In the ﬁeld of DFT calculations, the latter operator is known as the

independent-particle susceptibility operator and is denoted by χ0. This yields the Dyson equation

δρ =χ0(δV0+KHxc

ρδρ)⇔δρ =1−χ0KHxc

ρ−1χ0δV0.(19)

This equation is commonly solved by iterative methods, which require eﬃcient and robust computations of χ0δV

for various right-hand sides δV ’s. In the rest of this article, we forget about the solution of (19) and focus on the

computation of the noninteracting response δρ :=χ0δV for a given δV .

The operator χ0maps δV to the ﬁrst-order variation δρ of the ground-state density of a noninteracting system

of electrons (KHxc = 0). Denoting Amn :=hφm, Aφnifor a given operator A, it holds

δρ(r) =

+∞

n=1

+∞

m=1

fn−fm

εn−εm

φ∗

n(r)φm(r)(δVmn −δεFδmn),(20)

where δmn is the Kronecker delta, δεFis the induced variation in the Fermi level and we use the following convention

fn−fn

εn−εn

Tf0εn−εF

T=:f0

n.(21)

文档加载中……请稍候！
如果长时间未打开，您也可以点击刷新试试。

下载文档到电脑，查找使用更方便

10 玖币 0人已下载

立即下载

摘要：

NUMERICALSTABILITYANDEFFICIENCYOFRESPONSEPROPERTYCALCULATIONSINDENSITYFUNCTIONALTHEORYERICCANCÈS1;2,MICHAELF.HERBST3,GASPARDKEMLIN1;2,ANTOINELEVITT1;2,ANDBENJAMINSTAMM4Abstract.Responsecalculationsindensityfunctionaltheoryaimatcomputingthechangeinground-statedensityinducedbyanexternalperturbation.At...

展开>> 收起<<

NUMERICAL STABILITY AND EFFICIENCY OF RESPONSE PROPERTY CALCULATIONS IN DENSITY FUNCTIONAL THEORY ERIC CANCÈS12 MICHAEL F. HERBST3 GASPARD KEMLIN12 ANTOINE LEVITT12 AND BENJAMIN STAMM4.pdf

共21页,预览5页

还剩页未读，继续阅读

声明：本站为文档C2C交易模式，即用户上传的文档直接被用户下载，本站只是中间服务平台，本站所有文档下载所得的收益归上传人(含作者)所有。玖贝云文库仅提供信息存储空间，仅对用户上传内容的表现方式做保护处理，对上载内容本身不做任何修改或编辑。若文档所含内容侵犯了您的版权或隐私，请立即通知玖贝云文库，我们立即给予删除！

NUMERICAL STABILITY AND EFFICIENCY OF RESPONSE PROPERTY CALCULATIONS IN DENSITY FUNCTIONAL THEORY ERIC CANCÈS12 MICHAEL F. HERBST3 GASPARD KEMLIN12 ANTOINE LEVITT12 AND BENJAMIN STAMM4

相关推荐

开通VIP享超值会员特权

作者详情

相关内容

热门标签

举报选择: