{'fraction_non_alphanumeric': 0.03958723070424235, 'fraction_numerical': 0.01281354977571259, 'mean_word_length': 4.20129734254028, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical stability and efficiency of response property calculations in density functional theory

Response calculations in density functional theory aim at computing the change in ground-state density induced by an external perturbation. At finite 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-field calculations, to improve the stability and efficiency 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.


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. It offers a favourable compromise between accuracy and computational efficiency 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 coefficients, correspond physically to the response of GS quantities to a change in external parameters (e.g. nuclear positions, electromagnetic fields). As such their mathematical expressions involve derivatives of the obtained GS solution with respect to these parameters. For example interatomic forces are first-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 n th -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 defining the exchange-correlation functional [27,29,34].
Efficient 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 first-and higher-order derivatives computed by DFPT in electronic structure calculation software can be greatly simplified by Automatic Differentiation techniques [19], the efficiency of the resulting code crucially depends on the efficiency of a key building block: the computation of the linear response δρ of the GS density to an infinitesimal variation δV of the total Kohn-Sham potential.
For reasons that will be detailed below, the numerical evaluation of the linear map δV → δρ is not straightforward, 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 first consider the case of a periodic simulation cell Ω containing an even number N el 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 finite temperature T > 0 is then described by an L 2 (Ω)-orthonormal set of orbitals (φ n ) n∈N * with energies (ε n ) n∈N * , which are the eigenmodes of the Kohn-Sham Hamiltonian H associated with the GS density: Hφ n = ε n φ n ,ˆΩ φ * m (r)φ n (r)dr = δ mn , ε 1 ε 2 ε 3 · · · , together with periodic boundary conditions. The GS density in turn reads where f is a smooth occupation function converging to 2 at −∞ and to zero at +∞, e.g. the Fermi-Dirac smearing function f (x) = 2 1+e x (see Figure 1). The Fermi level ε F is the Lagrange multiplier of the neutrality charge constraint: it is the unique real number such that It follows from perturbation theory that the linear response δρ of the density to an infinitesimal variation δV of the total Kohn-Sham potential is given by where χ 0 is 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 infinitesimal perturbation δV . It holds (see Section 3) where δV mn := φ m , δV φ n , δε F is the induced variation of the Fermi level ε F , and δ mn is the Kronecker symbol. We also use the convention In practice, these equations are discretized on a finite basis set, so that the sums in (1) and (2) become finite. Since the number of basis functions N b is 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 N b of orbitals. These orbitals are then computed using efficient iterative methods [38].
For insulating systems, there is a (possibly) large band gap between ε Np and ε 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 f becomes a step function (see Figure 1). The jump from 2 to 0 in the occupations occurs exactly when the lowest N p = N el /2 energy levels ε 1 · · · ε Np are occupied with an electron pair (two electrons of opposite spin). Thus, f n = 2 for 1 n N p and f n = 0 for n > N p . As a result, N can be chosen equal to the number of electron pairs N p without 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 infinite simulation cell), causing the denominators in the right-hand side of formula (2) to formally blow up. Calculations on metallic systems are thus done at finite temperature T > 0, in which case every orbital has a fractional occupancy f n ∈ (0, 2). However, since from a classical semiclassical approximation ε n tends to infinity as n 2/3 as n → ∞, and f decays very quickly, one can safely assume that only a finite number N of orbitals have nonnegligible occupancies. This allows one to avoid computing φ n for n > N . Under this approximation, a formal differentiation of (1) gives However, while the response δρ to a given δV is well-defined by (2), the set (δφ n , δf n ) 1 n N is not. Indeed, the Kohn-Sham energy functional being in fact a function of the density matrix γ = N n=1 f n |φ n φ n |, any transformation of (δφ n , δf n ) 1 n N leaving invariant the first-order variation 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 P the orthogonal projector onto Span(φ n ) 1 n N , the space spanned by the orbitals considered as (partially) occupied, and by Q = 1 − P the 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 n where: • δφ P n = P δφ n ∈ Ran(P ) can be directly computed via a sum-over-state formula (explicit decomposition on the basis of (φ n ) n N ). This contribution can be chosen to vanish in the zero-temperature limit, as in that case P δγ P = 0. At finite 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 where H is the Kohn-Sham Hamiltonian of the system. This equation is possibly very ill-conditioned for n = N if ε N +1 − ε N is very small. This paper addresses these two issues. First, we review and analyse the different gauge choices for δφ P n proposed 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 ) n N , 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 efficiency 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 δφ n and 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 different 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.

Mathematical framework
2.1. Periodic Kohn-Sham equations. We consider here a simulation cell Ω = [0, 1)a 1 + [0, 1)a 2 + [0, 1)a 3 with periodic boundary conditions, where (a 1 , a 2 , a 3 ) is a nonnecessarily orthonormal basis of R 3 . We denote by R = Za 1 + Za 2 + Za 3 the periodic lattice in the position space and by R * = Zb 1 + Zb 2 + Zb 3 with a i · b j = 2πδ ij the reciprocal lattice. Let us denote by the Hilbert space of complex-valued R-periodic locally square integrable functions on R 3 , endowed with its usual inner product ·, · and by H s where e G (r) = e iG·r / |Ω| is the Fourier mode with wave-vector G.
In atomic units, the KS equations for a system of N el = 2N p spin-unpolarized electrons at finite temperature read where H ρ is the Kohn-Sham Hamiltonian. It is given by where V is the potential generated by the nuclei (or the ionic cores if pseudopotentials are used) of the system, and V Hxc is an R-periodic real-valued function depending on ρ. The Hartree potential V H ρ is the unique zero-mean solution to the periodic Poisson equation −∆V H ρ (r) = 4π ρ(r) − 1 |Ω|´Ω ρ and the function V xc ρ is the exchange-correlation potential. H ρ is a self-adjoint operator on L 2 # (R 3 , 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 field (SCF) algorithms. These algorithms are based on successive partial diagonalizations of the Hamiltonian H ρn built 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 ε n and occupation number f n . At finite temperature T > 0, f n is a real number in the interval [0, 2) and we have where f is a fixed analytic smearing function, which we choose here equal to twice the Fermi-Dirac function: f (x) = 2/(1 + e x ). The Fermi level ε F is then uniquely defined by the charge constraint When T → 0, f ((· − ε F )/T ) → 2 × 1 {·<εF} almost everywhere, and only the lowest N p = N el /2 energy levels for which ε n < ε F are occupied by two electrons of opposite spins (see Figure 1): f n = 2 for n N p and f n = 0 for n > N p .
Remark 1 (The case of perfect crystals). Using a finite 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 simulations are limited to, say 10 3 − 10 4 electrons, on currently available standard computer architectures. Simulating in vacuo a small sample of the material containing, say 10 3 atoms, would lead to completely wrong results, polluted by surface effects since about half of the atoms would lay on the sample surface. Periodic boundary conditions are a trick to get rid of surface effects, at the price of artificial interactions between the sample and its periodic image. In the case of a perfect crystal with Bravais lattice L and unit cell ω, it is natural to choose a periodic simulation (super)cell Ω = Lω consisting of L 3 unit 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.
Using Bloch theory, the supercell eigenstates φ n can then be relabelled as φ n (r) = e ik·r u jk (r), where u jk now has cell periodicity, and equations (7)-(9) can be rewritten as where Here L * is the dual lattice of L and ω = R 3 /L * the first Brillouin zone of the crystal. In the thermodynamic limit L → ∞, we obtain the periodic Kohn-Sham equations at finite temperature 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 > 0 allows 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 first rewrite the Kohn-Sham equations (7) as the fixed-point equation for the density ρ where F is the potential-to-density mapping defined by with (ε n , φ n ) n∈N * an orthonormal basis of eigenmodes of − 1 2 ∆ + V and ε F defined by (10). The solution of (16) defines a mapping from V to ρ: the purpose of DFPT is to compute its derivative. Let δV 0 be a local infinitesimal perturbation, in the sense that it can be represented by a multiplication operator by a periodic function r → δV 0 (r). By taking the derivative of (16) with the chain rule, we obtain the implicit equation for δρ: where the Hartree-exchange-correlation kernel K Hxc ρ is the derivative of the map ρ → V Hxc In the field of DFT calculations, the latter operator is known as the independent-particle susceptibility operator and is denoted by χ 0 . This yields the Dyson equation This equation is commonly solved by iterative methods, which require efficient 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 χ 0 maps δV to the first-order variation δρ of the ground-state density of a noninteracting system of electrons (K Hxc = 0). Denoting A mn := φ m , Aφ n for a given operator A, it holds where δ mn is the Kronecker delta, δε F is the induced variation in the Fermi level and we use the following convention Charge conservation leads toˆΩ where δε n := δV nn . We refer to [1] for a physical discussion of this formula, and to [8,22,33], where it is proven rigorously using contour integrals.
Remark 2. Similar to the discussion above on the computation of perfect crystal employing Bloch theory, response computations of perfect crystals can be performed by decomposing δV 0 in its Bloch modes. This allows for the efficient computation of phonon spectral or dielectric functions for instance.
Remark 3. We restricted our discussion for simplicity to local potentials, but the formalism can easily be extended to nonlocal perturbations (such as the ones created by pseudopotentials in the Kleinman-Bylander form [30]).

2.3.
Planewave discretization and numerical resolution. In this paper we are interested in plane-wave DFT calculations of metallic systems. This corresponds to a specific Galerkin approximation of the Kohn-Sham model using as variational approximation space where N b denotes the dimension of the discretization space, linked to the cut-off energy E cut . Denoting by Π N b the orthogonal projection onto X N b for the L 2 # inner product, we then solve the discrete problem: where H ρ is the Kohn-Sham Hamiltonian (or one of its Bloch fibers). This discretization method for Kohn-Sham equations has been analysed for instance in [3].
We emphasize again the point that not all N b eigenpairs need to be computed. At zero temperature, only the N = N el /2 lowest energy Kohn-Sham orbitals need to be fully converged as they are the only occupied ones. At finite temperature, the number of bands with meaningful occupation numbers is usually higher than the number of electrons, but the fast decay of the occupation numbers allows to avoid computing all N b eigenpairs. Determining the number of bands to compute is not easy as, at finite temperature, we do not know a priori the number of bands that are significantly occupied. A standard heuristic is to fully converge 20% more bands than the number of electrons pairs during the SCF. For the response calculation we then select the number N of bands that have occupation numbers above some numerical threshold. On top of these bands, it is common in DFT calculations to add additional bands that are not fully converged by the successive eigensolvers. The main advantages of introducing these bands are: (i) they enhance the diagonalization procedure by increasing the gap between converged and uncomputed bands and (ii) adding extra bands is not very expensive when the diagonalization is performed with block-based methods, such as the LOBPCG algorithm [31].
3. Computing the response 3.1. Practical implementation. Using (20) as it stands is not possible because of the large sums. One possibility is to represent δρ through a collection of occupied orbital variations (δφ n ) 1 n N and occupation number variations (δf n ) 1 n N . One then has to make appropriate ansatz and gauge choices on the links between δρ and its representation. Differentiating the formula ρ(r) = N n=1 f n |φ n (r)| 2 , one gets Then, for n N , we expand f n δφ n into the basis (φ m ) m∈N . Defining yields where δφ Q n := Qδφ n and Q is the orthogonal projector onto Span(φ m ) N <m , the space spanned by the unoccupied orbitals. Plugging (27) into (25), we obtain, using symmetry between n and m, A first gauge choice can be made here. Using again the charge conservation, we get Given that we adapt δf n accordingly we can thus assume Γ nn = 0 for any 1 n N . We will make this gauge choice from this point, leaving the constraint N n=1 δf n = 0 to restrict possible choices of δf n . We now derive conditions on (Γ mn ) 1 n,m N , (δf n ) 1 n N and (δφ Q n ) 1 n N so that the ansatz we made is a valid representation of δρ, that is to say (28) coincides with (20). To this end, we rewrite (20) as where the terms f n , f m for which n, m > N + 1 have been neglected because of their small occupation numbers and we used the symmetry between n and m for the terms with 1 n N , m > N . From a term by term comparison between (28) and (30), we infer first from the n = m term and the gauge choice Γ nn = 0 that δf n = f n (δV nn − δε F ) = f n (δε n − δε F ). Note that, thanks to the definition (22) of δε F , charge conservation is indeed satisfied. Next, for the first sum to coincide between (28) and (30), we see that the Γ mn 's have to satisfy Finally, since δφ Q n ∈ Span(φ m ) N <m , we deduce from the last sum in (28) and (30) that δφ Q n can be computed as the unique solution of the linear system sometimes known in DFT as the Sternheimer equation [44]. Note that δφ Q N can be arbitrarily large, since ε N +1 − ε N may be arbitrarily small. However, this does not pose a problem in practice as δφ Q N is multiplied by f N (cf. (28)), which is very small.
Here, δf n = f n (δε n − δε F ), and δφ n is separated into two contributions: where (δφ P n , δφ Q n ) ∈ Ran(P ) × Ran(Q) with P the orthogonal projector onto Span(φ m ) 1 m N and Q = 1 − P . These two contributions are computed as follows: • δφ P n is computed via a sum-over-states m = n: where the Γ mn 's satisfy Γ mn + Γ nm = ∆ mn . An additional gauge choice has to be made as these constraints do not yet define Γ nm uniquely. We refer to this term as the occupied-occupied contribution. • δφ Q n is obtained as the solution of the Sternheimer equation (32). However, this linear system is possibly very ill-conditioned if ε N +1 − ε N is small. We refer to this term as the unoccupied-occupied contribution. Note that, at zero temperature, δφ P n vanishes so that δφ n = δφ Q n ∈ Ran(Q) and only the Sternheimer equation (32) needs to be solved. In the next two sections, we detail the practical computation of these two contributions.

3.2.
Occupied-occupied contributions. In this section we discuss possible gauge choices for Γ mn to obtain a unique solution to (31). Throughout this section we assume m = n and Γ nn = 0.
3.2.1. Orthogonal gauge. The orthogonal gauge choice is motivated from the zero temperature setting, where δφ P n = 0 allows to trivially preserve the orthogonality amongst the computed orbitals φ n under the perturbation. For the case involving temperature, we additionally impose and therefore yielding As a result f n δφ n features a possibly large contribution Γ orth mn , which is going to be almost compensated in (33) by the large contribution Γ orth nm to f n δφ * n due the requirement to sum to the moderate-size contribution Γ orth mn +Γ orth nm = ∆ mn . This can lead to numerical instabilities because small errors, e.g. due to the fact that the φ n 's in (33) are eigenvectors only up to the solver tolerance, will get amplified by the Γ mn . The next gauge choices provide solutions to this issue.

Simple gauge choice. Possibly the simplest gauge choice is
The DFPT framework presented in [1] and implemented in Quantum Espresso [12] suggests choosing where f FD = 1 2 f is the Fermi-Dirac functional. Since f FD (x) + f FD (−x) = 1, we have Γ QE mn + Γ QE mn = ∆ mn .

Abinit gauge.
In the Abinit software [17,41], the choice is 3.2.5. Minimal gauge. Motivated by our analysis of the instabilities we suggest minimizing δφ n , that is to ensure Γ mn /f n to stay as small as possible. This leads to the minimization problem min N n,m=1,m =n As the constraint (31) only couples (n, m) and (m, n), this translates into an uncoupled system of constrained minimization problems: for 1 n, m N , m = n, solve whose solution is This gauge choice is implemented by default in the DFTK software [25]. Another gauge choice inspired from this one would be to directly minimize |Γ mn | 2 + Γ nm 2 but it can be shown that this leads to the simple gauge choice 3.2.6. Comparison of gauge choices. From (20) we can see that the growth of δρ with respect to δV can not be higher than the growth of ∆ mn with respect to δV . The latter is of the order of max x∈R 1 T |f (x)| = 1 2T , which thus provides an intrinsic limit to the conditioning of the problem. For all gauge choices but the orthogonal one easily verifies |Γ mn | |∆ mn | max If we make an error on δV it is thus at most amplified by a factor of 1 2T . All choices but the orthogonal one thus manage to stay within the intrinsic conditioning limit, see Figure 2 3.3. Computation of unoccupied-occupied contributions employing a Schur complement. Since the φ m for m > N are not exactly known, a different approach is needed for obtaining the contribution δφ Q n . Usually one resorts to solving the Sternheimer equation using iterative schemes restricted to Ran(Q). However, for n = N the difference ε N +1 − ε N can become small, which deteriorates conditioning and increases the number of iterations required for convergence.
We overcome this issue by making use of the N ex extra bands, which are anyway available after the SCF algorithm has completed. Following Section 2.3 the N ex extra bands can be divided into two categories: (1) Some (usually the lower-energy ones) have been discarded during the response calculation because they have a too small occupation. Up to the eigensolver tolerance these are exact eigenvectors. (2) The remaining ones have served to accelerate the successive diagonalization steps during the SCF. These have not yet been fully converged. In any case these extra bands thus offer (at least) approximate information about some φ m for m > N , which is the underlying reason why the following approach accelerates the computation of δφ Q n . For the sake of clarity, we place ourselves here in the discrete setting: We assume that the number of computed bands N + N ex is larger than the number of occupied states N and that we trust Φ = (φ 1 , . . . , φ N ) but not Φ = ( φ N +1 , . . . , φ N +Nex ) to be eigenvectors. These N ex extra bands consist of both contributions (1) and (2) described at the beginning of this section. We assume in addition that (Φ, Φ) forms an orthonormal family and that Φ * H ρ Φ is a diagonal matrix whose elements, denoted by ( ε n ) n=N +1,...,N +Nex , are not necessarily all exact eigenvalues. Note that Rayleigh-Ritz based iterative methods such as the LOBPCG algorithm fit exactly in this framework. We decompose where T is the orthogonal projector onto Span( φ m ) N <m N +Nex and R = Q − T is the projector onto the remaining (uncomputed) states, see Figure 3. Then, as δφ Q n ∈ Ran(Q), we can decompose δφ Q n = Φα n + δφ R n , where α n ∈ C Nex and δφ R n ∈ Ran(R). Plugging this into (45) we get Using a Schur complement we deduce Inserting (49) into (48) and projecting on Ran(R) yields an equation in δφ R n : This equation can then be solved for δφ R n with a Conjugate Gradient (CG) method which is enforced to stay in Ran(R) at each iteration. Afterwards we compute α n from (49), which yields δφ Q n from (47). This scheme has been implemented as the default solver for the Sternheimer equation in DFTK.

Numerical tests
For all the numerical tests, we use the DFTK software [25], a recent plane-wave DFT package in Julia. All the codes to run the simulation of this paper are available online 1 . The Brillouin zone is discretized using a uniform Monkhorst-Pack grid [36]. We use the PBE exchange-correlation functional [39] and GTH pseudopotentials [13,20]. The other parameters of the calculation will be specified for each example. In all the tests, we generate a perturbation δV from atomic displacements, with local and nonlocal contributions. Then, we perform two response calculations: one with the standard approach to solve directly the Sternheimer equation (45) to compute δφ Q n , the other with the (new) Schur complement approach (50). Both linear systems are solved using the conjugate gradient (CG) algorithm, with kinetic energy preconditioning (the linear solver is preconditioned with the inverse Laplacian, which is diagonal in Fourier representation), and we compare the number of iterations required to converge the norm of the residual below 10 −9 . Note that the Sternheimer equation is solved for all N occupied orbitals and for each k-point.
If T > 0 the contribution δφ P n is nonzero and has been computed using the sum-over-states formula with the minimal gauge choice (43). In terms of runtime we expect only negligible differences between the gauge choices. Moreover, since the time for this contribution is much smaller compared to the time required to solve the Sternheimer equation, we do not report a detailed performance comparison on this step in the following.

Insulators and semiconductors. For insulators and semiconductors the gap between occupied and virtual
states is usually large. One would therefore not expect a large gain from using the Schur complement (50) when computing δφ n = δφ Q n . However, for distorted semiconductor structures or semiconductors with defects the gap can be made arbitrarily small, such that one would expect to see the Schur-complement approach to be in the advantage. We test this using an FCC Silicon crystal for which we increase the lattice constant from 10 bohrs to 11.4 bohrs to artificially decrease and eventually close the gap. All calculations have been performed using a cut-off energy of E cut = 50 Ha and a single k-point (the Γ-point). In Figure 4 we plot the number of iterations required for the linear solver of the Sternheimer equation to converge, for n = N = 4. Using the Schur complement the number of iterations stays almost constant even when the gap decreases. In contrast, with a direct approach, the linear solver requires about 30% more iterations near the closing gap.

4.2.
Metals. The real advantage of using the Schur complement (50) instead of directly solving the Sternheimer equation (45) becomes apparent when computing response properties for metals at finite temperature. We use a standard heuristic which suggests to fully converge 20% more bands than the number of electrons pairs of the system, with 3 additional extra bands that are not converged by the successive eigensolvers of the SCF. We then select the "occupied" orbitals with an occupation threshold of 10 −8 . In addition to the number of iterations, we also compare the cost of the response calculations with and without the Schur complement (50). For this we consider the total number of Hamiltonian applications which was required to compute the response δρ. For the small to mediumsized systems we consider here, the Hamiltonian-vector-product is the most expensive step in an DFT calculation and thus provides a representative cost indicator. Notice that both the implementation of the Schur complement and the direct method require exactly one Hamiltonian application per iteration of the CG. Additionally the Schur approach requires the computation of H ρ φ, which is only a negligible additional cost as this is only needed once per k-point.

Aluminium.
We start by considering an elongated aluminium supercell with 40 atoms. We use a cut-off energy E cut = 40 Ha, a temperature T = 10 −3 Ha with Fermi-Dirac smearing and a 3 × 3 × 1 discretization of the Brillouin zone. To ensure convergence of the SCF iterations we employ the Kerker preconditioner [28]. Since the system has 120 electrons per unit cell our usual heuristic converges 72 bands up to the tolerance of the eigensolver accompanied by 3 bands, which are not fully converged.
The convergence behaviour when solving the Sternheimer equation for k-points of particular interest is shown in Table 1 and Figure 5. As expected, for k-points with a small difference ε N +1 −ε N , the Schur complement (50) brings a noteworthy improvement with roughly 50% fewer iterations required to achieve convergence. Since the system we consider here has numerous occupied bands -between 60 and 70 depending on the k-point -most bands already feature a well-conditioned Sternheimer equation. Considering the cost for computing the total response, the Schur approach therefore overall only achieves a reduction by 17% in the number of Hamiltonian applications, from about 17, 800 (direct) to 14, 800 (with Schur). However, it should be noted that this improvement essentially comes for free as the extra bands are anyway provided by the SCF computation as a byproduct.  Table 1 -Convergence data for k-points 1, 2 and 5 for Al 40 . Other k-points are not displayed but they all behave as one of these by symmetry. N is the number of occupied bands, for an occupation threshold of 10 −8 .

4.2.2.
Heusler system. Next we study the response calculation of Heusler-type transition-metal alloys. We focus mainly on the Fe 2 MnAl system but other compounds, such as the Fe 2 CrGa and CoFeMnGa alloy systems, have been tested and similar results were obtained. Heusler alloys are of considerable practical interest due to their rich and unusual magnetic and electronic properties. For instance, Fe 2 MnAl shows halfmetallic behaviour: the majority spin channel (denoted by ↑) behaves like a metal whereas the minority spin channel (denoted by ↓) behaves like an insulator as it has a vanishing density of states at the Fermi level. See [23], and reference therein, for more details as well as an analysis of the SCF convergence on such systems. For these systems we use a cut-off E cut = 45 Ha, a temperature T = 10 −2 Ha with Gaussian smearing and a 13 × 13 × 13 discretization of the Brillouin zone. The SCF was converged using a Kerker preconditioner [28]. Moreover, as we deal with a spin-polarized system, the numerical simulation slightly differs. The orbitals φ σ (n,k) and the occupation numbers f σ (n,k) depend on the spin orientation σ ∈ {↑, ↓} and the f σ (n,k) 's belong to [0, 1) instead of [0, 2). Furthermore we modify the heuristic to determine the number of bands to be computed: Fe 2 MnAl has N el = 50 electrons per unit cell and we use 25 + 0.2 × 50 = 35 fully converged bands per k-point, complemented by 3 additional bands, which are not checked for convergence.
We show in Table 2 and Figure 6 the results for the two spin channels of the k-point with reduced coordinates [0.385, 0.231, 0.077]. The other k-points behave similarly. Since both channels feature a small difference ε N +1 − ε N using the Schur complement (50) to solve the Sternheimer equation has a significant impact: for the orbitals with highest energy it reduces the number of iterations by half. For the direct approach we notice a plateau where the solver encounters difficulties to converge the Sternheimer equation for the N -th orbital due to the small gap. Unlike the aluminium case the improvements observed for the Heusler alloys are not restricted to a small number of bands. In Figure 7 we contrast the number of iterations required to solve the Sternheimer equation for every band at every k-point with and without using the Schur complement. Notice that lattice symmetries allow to reduce the number     Figure 8 -Resolution of the Sternheimer equation for both spin channels of one specific k-point for the Fe 2 MnAl system, with the Schur approach (50) and the shifted approach (52). Note that for this particular k-point, the spin ↓ channel has a starting point with already small residual for the highest occupied energy level. This is due to the term fn appearing in (52), and the convergence is still slow.
of explicitly treated k-points to 140 albeit we are using a 13 × 13 × 13 k-point grid. In terms of the total number of Hamiltonian applications required for the response calculation, the Schur complement achieves a reduction by roughly 40%, from around 344, 000 (without Schur) to 208, 000 (with Schur). It should be noted that in this system the standard heuristic caused a large portion of the available extra bands to be fully converged, thus providing an ideal setting for the Schur complement approach to be effective. For example for the k-point discussed in Table 2 seven extra bands have been fully converged and an additional three partially. Given the enormous importance of Heusler systems and the known numerical difficulties for computing response properties in these systems, our result is encouraging and motivates the development of a more economical heuristic for choosing the number of converged bands in future work.

4.3.
Comparison to shifted Sternheimer approaches. In the literature other strategies for computing δρ have been reported. We briefly consider the approach proposed in [1], where the response is computed as Instead of splitting δφ n into two contributions, the full δφ n is computed for all n N by solving the shifted Sternheimer equation (52) Here S : Ran(P ) → Ran(P ) is a shift operator acting on the space of occupied orbitals, chosen so that the linear system is nonsingular (for any n N , H ρ − ε n is not invertible). Then, S n is chosen for every n N such that δρ from (51) satisfies (20). However, as S only acts on Ran(P ), equation still becomes badly conditioned if ε N +1 − ε N is too small. This becomes apparent when solving the shifted Sternheimer equation (52) for the Fe 2 MnAl system, see Figure 8. For the orbital responses of the highest-energy occupied bands the CG iterations on the shifted Sternheimer equation converge very slowly -in contrast to the Schur complement approach (50) we proposed in this work. In terms of the number of Hamiltonian applications, the shifted Sternheimer strategy required around 492, 000 applications versus 208, 000 for the Schur complement approach.

Conclusion
In density-functional theory the simulation of many physical properties requires the computation of the response of the ground-state density to an external perturbation. In this work we have reviewed the standard formalism of such response calculations from the point of view of numerical analysis. We provided an overview of the possible gauge choices for representing the density response, summarizing and contrasting the approaches employed by state-of-the-art codes such as Quantum Espresso [12] or Abinit [41] in a common framework.
Based on our analysis we furthermore suggested two novel approaches for DFT response calculations. For the occupied-occupied part of the response we developed a gauge choice based on the idea to maximize numerical stability in the involved sums by minimizing the numerical range of the individual orbital contributions. For the occupied-unoccupied part of the response we suggested a novel approach to solving the Sternheimer equation based on a Schur complement. Key idea of this approach is to make use of the additional (partially) converged bands, which are available as a byproduct from the preceding self-consistent field (SCF) procedure (which yields the groundstate density). Without additional computational effort this allows to improve the conditioning of the Sternheimer equation and thus accelerate its convergence. We demonstrated this numerically on a number of practically relevant problems, including response calculations on small-gapped semiconductors, elongated metallic slabs or numerically challenging Heusler alloy systems. Overall the Schur complement approach allowed to obtain a converged response saving up to 40% in the required Hamiltonian applications -the cost-dominating step in small to medium-sized DFT problems. For larger systems we similarly expect savings from introducing a Schur complement technique, even though algorithms commonly employ different trade-offs.
In this work we followed standard heuristics for selecting the number of extra bands to employ in the SCF calculations and thus the number of additional bands available when solving the response problem. However, our results emphasize the need for a more robust understanding between the computed number of bands and the observed rate of convergence. We have provided some initial ideas for such an analysis in the appendix, but leave a more exhaustive discussion for future work.
CG is indeed led by the square root of the condition number of the system matrix (see [43,Section 9]) when using the Schur complement. Key idea will thus be to estimate the condition number of the linear system (50).
A.1. Numerical analysis. To analyse the condition number of the Schur complement, we consider the following specific setting where H ρ ∈ C N b ×N b is typically the discretized self-consistent Hamiltonian of the system, at some k-point. We assume that we have N occupied orbitals that have an occupation number higher than the threshold we fixed and that we have N ex extra bands, as explained in Section 3.3. In summary, we have at our disposal N + N ex bands in total: Φ = (φ 1 , . . . , φ N ) are occupied, fully converged bands and the extra bands Φ ex = (φ N +1 , . . . , φ N +Nex ) are not necessarily all converged. We added here the exponent as we make the following assumptions: • for any ∈ N, (Φ, Φ ex ) is an orthonormal family; • for any ∈ N, (Φ ex ) * H ρ Φ ex ∈ C Nex×Nex is a diagonal matrix whose elements are labelled ε m := φ m , H ρ φ m for N + 1 m N + N ex ; • as → +∞, (φ m , ε m ) → (φ m , ε m ). All these assumptions are satisfied for instance if the sequence (Φ, Φ ex ) l∈N is generated by any Rayleigh-Ritz based eigensolver (for instance the LOBPCG eigensolver [31]), which is the case by default in DFTK. For every , we can thus decompose the plane-wave approximation space H = X N b (with N b N + N ex ) in two different ways: where P := N n=1 φ n φ * n and are all orthogonal projectors. In these two decompositions, H ρ has the associated block representations: where E := Diag(ε 1 , . . . , ε n ), E ex := Diag(ε N +1 , . . . , ε N +Nex ) and E ex := Diag(ε N +1 , . . . , ε N +Nex ) are diagonal matrices. Moreover, note that as Φ ex → Φ ex , the residuals R H ρ T converge to 0. Now, we fix n N and we compute the condition number of the linear system (50). Enforcing the CG to stay at each iteration in Ran(R ), this condition number is given by the ratio of the largest and smallest nonzero eigenvalues of where Here E ex − ε n is diagonal and thus explicitly invertible if is large enough as ε N +1 → ε N +1 > ε N ε n . We focus for the moment on the smallest nonzero eigenvalue, that is ε N +Nex+1 − ε n . The condition number being proportional to the inverse of the smallest eigenvalue, we now derive a lower bound of ε N +Nex+1 − ε n in order to get an upper bound on the condition number of (57). When → +∞, we have X n → 0 (as RH ρ Φ ex = 0) and H n → H n := R(H ρ − ε n )R whose smallest nonzero eigenvalue is ε N +Nex+1 − ε n . We use next a perturbative approach to effectively approximate the condition number of (57).
We use a standard eigenvalue perturbation result, whose proof is recalled for the sake of completeness. It is directly adapted from the general case of self-adjoint bounded below operators with symmetric perturbations studied for instance in [11]. Proposition 1. Let N ∈ N, H 0 , W ∈ C N ×N be Hermitian matrices and α 0 such that H 0 + α > 0. Then, the eigenvalues of H := H 0 + W and H 0 satisfy where W H0,α is the operator norm of (H 0 + α) −1/2 W (H 0 + α) −1/2 and ν i (A) is the i-th lowest eigenvalue of the matrix A.
In our case, we can apply this result to with H 0 = H n , W = W n := (H n − H n ) + X n and α = ε N +Nex+1 − ε n > 0. Proposition 1 applied to the (N + N ex + 1)-th eigenvalues then yields where we assume that 2 W n Hn,ε N +Nex +1 −εn is small enough to be negligible with respect to 1, which is the case if the extra states are sufficiently converged. Now, if this bound is valid in theory, in practice we do not have access to ε N +Nex+1 as we work with N + N ex bands only. However, up to loosing sharpness, we can use that ε N +Nex+1 ε N +Nex where ε N +Nex can be estimated using the last extra band. Indeed, using for instance the Bauer-Fike bound ([24, Theorem 1] or [42]), we obtain where r N +Nex is the residual associated to the last extra band. Of course, this estimate is not sharp as we expect the error on the eigenvalue to behave as the square of the residual, but this requires to estimate the gap to the rest of the spectrum, see for instance the Kato-Temple bound [24,Theorem 2]. In the end, we have the following lower bound for ε N +Nex+1 − ε n : where we assume again that r N +Nex is small enough with respect to ε N +Nex − ε n .
We can now derive an upper bound on κ n , the condition number of (57). It is given by the ratio of its highest eigenvalue and ε N +Nex+1 − ε n . Since the Laplace operator is the higher-order term in the Kohn-Sham Hamiltonian, the highest eigenvalue is, as usually in plane-wave simulations, of order E cut . With proper kinetic preconditioning, we can assume that its contribution to the condition number of the linear system is constant with respect to E cut and n so that, finally, Therefore the condition number is bounded from above by C/(ε N +Nex − ε n ) to first-order. The number of CG iterations to solve the linear system (50) with a given accuracy is then proportional to the square root of the condition number of the matrix (57) (see [43]): Note that this upper bound is valid provided that the extra bands are converged enough, not necessarily fully, and proper kinetic preconditioning is employed.
Estimate (68) leads, as expected, to the qualitative conclusion that the more extra bands we use, the higher the difference ε N +Nex − ε n and the faster the convergence. However, note that it is not possible to evaluate directly the convergence speed as the constant C is a priori unknown, in particular if we use preconditioners.
A.2. An adaptive strategy to choose the number of extra bands. The main bottleneck of (68) is the estimation of the constant C. However, one can reasonably assume that this constant does not depend too much on n, so that the ratio between the number of iterations to reach convergence between the last occupied band (n = N ) and the first band (n = 1) can be estimated by This ratio can be of interest as (68) suggests that the Sternheimer solver converges the fastest for n = 1 and the slowest for n = N .
We plot in Figure 9 the upper bound ξ Nex as well as the computed ratios between the number of iterations of the first and last bands for the systems we considered in Section 4. These plots show that ξ Nex is indeed an upper bound of the actual ratio. This bound does not seem to be sharp however. This is due to the successive approximations we made to obtain this estimate. Plots in Figure 9 also confirm that if, for every k-point, the ratio of the number of iterations between the first and last occupied bands is assumed to be an accurate indicator of the efficiency of the Sternheimer solver, then using the Schur complement (50) always make this ratio smaller.
If we want the ratio of the number of iterations between the first and the last occupied bands to be lower than some target ratio ξ T (for instance 3), Figure 9 suggests that the computable ratio ξ Nex can help in choosing the number of extra bands to reach this target ratio. We propose in Algorithm 1 an adaptive algorithm to select the number of extra bands as a post-processing step after termination of the SCF. The basic idea is that, given the initial output (Φ, Φ ex ) with = 0 of an SCF calculation, one iterates Φ ex → Φ +1 ex where Φ ex gathers the extra bands. At each iteration , we compute ξ Nex and check if it is below the target ratio. If not, we compute more approximated eigenvectors, that we converge until the residual r N +Nex is negligible with respect to ε N +Nex − ε N , and so on. To generate such a residual, after adding a random extra band properly orthonormalized, we update the extra bands using a LOBPCG with tolerance tol = (ε N +Nex−1 − ε N )/50. It appears that Algorithm 1 can be used to choose the number of extra bands in order to reach a given ratio ξ T . However, using the algorithm as such is not useful in practice as it requires a too high number of Hamiltonian applications, making this strategy less interesting than the Schur approach we proposed with the default choice of extra bands. Strategies to reduce the number of Hamiltonian applications in order to choose an appropriate number of extra bands will be subject of future work.