Computing the self-consistent field in Kohn–Sham density functional theory

The concept of self-consistency has been prevalent across many domains of physics, typically as a characteristic requirement when one invokes a mean-field approximation. For example, Hartree theory replaces the two-body Coulomb interaction between electrically-charged quantum particles with a mean-field, the Hartree potential, generated by the distribution of the electric charge in the system. Each particle is influenced by the Hartree potential, which in turn alters the distribution of charge in the system. This charge distribution can then be used to construct a new Hartree potential. The Hartree potential is self-consistent when these two fields are the same, i.e. the potential leads to a charge distribution which gives rise to the same potential. In fact, this was the context in which self-consistency was first introduced,

'If the final field is the same as the initial field, the field will be called 'self-consistent', and the determination of self-consistent fields for various atoms is the main object of this paper'.

–D.R. Hartree (1927) [20].

Later refined by Fock [21] and Slater [22], Hartree–Fock theory became widely adopted in computational quantum chemistry to compute ground state properties of molecules [3]. Whilst Hartree and Hartree–Fock theory are mean-field approximations, Hohenberg et al [5, 23] showed that a mean-field exists which reproduces the ground-state energy and particle density exactly. This 'density functional theory' allows, in principle, the computation of the exact electronic structure of any quantum system; however the exact density functional is not known, and must be approximated in any practical application of DFT. For a more detailed examination of the origins and physical foundations of Kohn–Sham theory, the reader is directed to the following resources [1, 4, 24], and references therein.

This work concerns the need to achieve self-consistency in the context of DFT simulations of atoms, molecules and materials. Namely, we focus on computing the particle density $\rho(x)$ for a set of atomic species and positions within the framework of Kohn–Sham DFT. Each of the N particles in the system are influenced by an external potential $v_{{\rm ext}}$ which is uniquely defined by the species and positions of the atoms, the level of approximation employed, and more. For the purposes of this article, finding the ground state energy E in Kohn–Sham theory is viewed as a constrained minimisation problem,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{kohnsham_min} E[v_{{\rm ext}}] = \inf_{\{\phi_i \}} \Big\{E_{\tiny{\rm KS}}[\{\phi_i \} ] \ \Big| \ \phi_i \in \mathcal{H}^1(\mathbb{R}^3) & \nonumber \\ \int_{\mathbb{R}^3} \phi^\ast_i(x) & \phi_j(x) = \delta_{ij} \nonumber \\ \nonumber & i,j \in [1,N] \Big\}, \nonumber \end{align} \tag{ 1 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E_{\tiny{\rm KS}}[\{\phi_i \} ] = \sum_{i=1}^N & \frac{1}{2} \int_{\mathbb{R}^3} \ |\nabla \phi_i(x) |^2 + \frac{1}{2}\int_{\mathbb{R}^3 \times \mathbb{R}^3} \frac{\rho(x)\rho(x')}{|x-x'|} \nonumber \\ &+ \int_{\mathbb{R}^3} \rho(x) v_{{\rm ext}}(x) + E_{{\rm xc}}[\rho], \nonumber \end{align} \tag{ 2 }$$

where atomic units are used, and, for now, spin degrees of freedom are omitted. The particle density, $\rho(x)$, is defined in terms of the single-particle orbitals, $\{\phi_i \}$, via

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho(x) = \sum_{i=1}^N |\phi_i(x)|^2. \label{purestate_density} \nonumber \end{align} \tag{ 3 }$$

That is, one must minimise the Kohn–Sham objective functional equation (1) over a set of N orthogonal, normalisable functions $\{\phi_i \}$ whose first derivative is also normalisable, i.e. they exist in the Sobolev space $\mathcal{H}^1(\mathbb{R}^3)$. The exchange-correlation functional $E_{{\rm xc}}$ is a yet undetermined functional of the density designed to capture the effects of exchange and correlation missing from the remainder of the functional. In principle, the Hohenberg–Kohn theorems guarantee that the Kohn–Sham objective functional is a functional of the density alone [23]. However, in the case of Kohn–Sham theory, recourse to an orbital-dependent functional is necessitated by the definition of the single-particle kinetic energy.

Explicit constrained variation of the orbitals allows one to approach the optimisation problem in equation (1) directly. This can be done, for example, with a series of line searches in the direction of steepest descent of $E_{\tiny{\rm KS}}$ with respect to the orbitals [25, 26]. Alternatively, assuming differentiability [27], the associated Lagrangian problem can be formulated, and the functional derivative of the Lagrangian set to zero. This yields the Euler–Lagrange equations for the problem, the solution of which is a stationary point of the functional. In the present context, the Euler–Lagrange equations constitute a nonlinear eigenvalue problem,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H_{\tiny{\rm KS}}[\rho] \phi_i(x) = \epsilon_i \phi_i(x), \nonumber \end{align} \tag{ 4 }$$

where the Hamiltonian operator $H_{\tiny{\rm KS}}$ depends on its eigenvectors via

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H_{\tiny{\rm KS}}[\rho] = -\frac{1}{2}\nabla^2 + v_{{\rm ext}} + v_{\rm h}[\rho] + v_{{\rm xc}}[\rho], \nonumber \end{align} \tag{ 5 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v_{\rm h}[\rho](x) = \int_{\mathbb{R}^3} \ \frac{\rho(x')}{|x-x'|}, \nonumber \end{align} \tag{ 6 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v_{{\rm xc}}[\rho] = \frac{\delta E_{{\rm xc}}}{\delta \rho}. \nonumber \end{align} \tag{ 7 }$$

These are the Kohn–Sham equations. The eigenvalues (quasi-particle energies) $\newcommand{\e}{{\rm e}} \epsilon_i$ are the Lagrange multipliers associated with the orbital orthonormality constraint. Solving the Kohn–Sham equations to find a stationary point of the Kohn–Sham functional is a necessary but not sufficient condition for (local) optimality. A sufficient condition would require the second derivative (curvature) about the stationary point to be everywhere positive. Furthermore, in general, the Kohn–Sham functional for some approximate $E_{{\rm xc}}$ is not a convex functional of the orbitals, meaning that verifying global optimality is a difficult task. In practice, solving the Kohn–Sham equations with certain methods of biasing the solution toward a (possibly local) minimum are often chosen rather than direct minimisation methods [28]. The advantages and drawbacks of each approach will be examined in section 4.

It is now possible to formally define what is meant by self-consistency. In order to construct the Kohn–Sham Hamiltonian, one requires a density as input $\rho^{{\rm in}}$ to compute the Hartree and exchange-correlation potentials. An output density $\rho^{{\rm out}}$ is then calculated (non-linearly) from the eigenfunctions of the Kohn–Sham Hamiltonian

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H_{\tiny{\rm KS}}[\rho^{{\rm in}}]\phi_i(x) = \epsilon_i \phi_i(x), \label{rho_in} \nonumber \end{align} \tag{ 8 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{{\rm out}}(x) = \sum_{i=1}^N |\phi_i(x)|^2. \label{rho_out} \nonumber \end{align} \tag{ 9 }$$

In general, the input density is not equal to the output density. For a given external potential and exchange-correlation functional, the density $\rho_*$ is self-consistent when $\rho_* = \rho^{{\rm in}} = \rho^{{\rm out}}$, and hence the non-linear eigenvalue problem of equations (8) and (9) is solved. The non-linearity in equation (9) necessitates an iterative procedure that takes an initial estimate of the density as input and iterates this density toward a self-consistent solution of the Kohn–Sham equations: the self-consistent field procedure, figure 1. As one might expect, an infinity of self-consistent densities exist for a given external potential and exchange-correlation functional [29]. However, we are interested primarily in the subset of these densities that are local minima of the Kohn–Sham objective functional.

Zoom In Zoom Out Reset image size

Modern computational implementations of Kohn–Sham theory can vary significantly due to various factors. The key distinguishing factor is the choice of basis set, which leads to the related problem of whether one treats all the electrons in the computation explicitly, or treats core electrons with a pseudopotential [30]. Despite these differences, perhaps with the exception of linear scaling methods [14], the self-consistent field techniques to be discussed here are adaptable to most implementations. Indeed, some of the most popular software, such as vasp [15, 31], abinit [10, 32], quantum espresso [33], and castep [34], use similar default methods to achieve self-consistency: preconditioned multisecant methods, which are discussed in section 4.

The overarching goal of this work is to quantify the utility of a given algorithm for reaching self-consistency in Kohn–Sham theory. In turn, this allows us to compare and analyse the performance of a sample of existing algorithms from the literature. Assessing these algorithms requires the creation of a test suite of Kohn–Sham inputs, representative of a range of numerical issues. This test suite generates a standard which can be used to test, improve, and present new algorithms designed by method developers. Furthermore, the test suite allows DFT developers to more effectively assess which algorithms they wish to implement. With these aims in mind, this article is structured in two partitions, as follows.

The first part constitutes a review of self-consistency in Kohn–Sham theory. As such, the relevant sections are ideal for an interested party who is not actively involved in development to gain a more in-depth understanding of self-consistency from an algorithmic perspective. In particular, this review collates decades of past literature on self-consistency in Kohn–Sham theory, thus elucidating conclusions that have become conventional wisdom. Section 3 examines the Kohn–Sham framework abstractly from a mathematical and computational perspective in order to study where and why algorithms encounter difficulties. This involves, for example, a discussion on the nature of so-called 'charge-sloshing', the initial guess density, sources of ill-conditioning, and more. Section 4 then examines and categorises the range of available algorithms in present literature. A focus will be placed on detailing the algorithms which have proven to be particularly successful.

The second part then utilises the analysis presented in the prior sections to perform a study akin to recent benchmarking efforts such as GW100 [35] and the $\Delta$-project [36] for assessing reproducibility in GW and DFT codes respectively. Whilst the study presented here will take on a similar structure to these examples, it differs in the following way. The aim of the $\Delta$-project is to assess an 'error' for each DFT software in a given computable property compared to a reference software over a set of test systems. Here, we instead aim to assess the utility of an algorithm, rather than an error, which we do with two competing measures: efficiency and robustness, defined in section 5. A test suite of Kohn–Sham inputs is then constructed to target weaknesses in contemporary algorithms and exploit the difficulties discussed in section 3. This test suite is designed to be representative of the range of systems practitioners may encounter and with which they may have difficulties reaching convergence. Each algorithm is then assigned a robustness and efficiency score when tested over the full suite. The methods of Pareto analysis then provide a prescription for the definition of optimal when there exist two or more competing measures of utility. Section 6 demonstrates these concepts by using this workflow on a selection of algorithms described in section 4, implemented in the plane-wave, pseudopotential software castep. This study allows conclusions to be drawn about the current state of self-consistency algorithms in Kohn–Sham codes. Finally, we discuss how one might utilise the test suite and workflow demonstrated here to present and assess future methods and algorithms.

The central approximation involved in converting the framework of Kohn–Sham theory into a form suitable for computation is called the finite-basis approximation, or the Galerkin approximation [28]. The orbitals $\phi_i$ are continuous functions of a continuous three dimensional variable, x. These functions are equivalent to vectors existing in an infinite dimensional vector space, spanned by a complete basis $\chi_\mu$. Provided this basis does indeed span the space, the orbitals can be expressed exactly as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \phi_i(x) = \sum_{\mu=1}^\infty \alpha_{i \mu} \chi_\mu (x). \nonumber \end{align} \tag{ 10 }$$

Note that the basis set is required to be complete, but not necessarily orthogonal,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int_{\mathbb{R}^3} \chi_\mu(x) \chi_\nu (x) = S_{\mu \nu} = \delta_{\mu \nu}. \nonumber \end{align} \tag{ 11 }$$

If the basis is non-orthogonal, the finite-basis Kohn–Sham equations generalise to include the overlap matrix, $S_{\mu \nu} \neq \delta_{\mu \nu}$, see [28] for more details. This distinction does not affect the ensuing analysis, and so the basis set is hereafter taken to be orthogonal for clarity. Once the basis is specified, the equations can be rearranged and solved for the infinity of coefficients to the basis $\alpha_{\mu i}$. In practice, one must truncate the basis such that it is no longer complete and instead captures only the most relevant regions of the formally infinite Hilbert space,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \phi_i(x) \approx \sum_{\mu=1}^{N_b} \alpha_{i \mu} \chi_\mu (x). \nonumber \end{align} \tag{ 12 }$$

The characteristic size of the basis N_b will depend primarily on the choice of basis functions. Within the finite-basis approximation, the Kohn–Sham Hamiltonian becomes an $N_b \times N_b$ matrix, of which a subset of the eigenvalues and eigenvectors is required to progress toward a solution of the non-linear eigenvalue problem, equations (8) and (9). Basis functions which are localised about the atomic cores [37] are a popular choice. These tend to be particularly accurate per basis function, meaning N_b is typically the same order of magnitude as the number of electrons, N. Methods utilising local basis functions are often able to form and diagonalise the Kohn–Sham Hamiltonian matrix explicitly. In such implementations, the Kohn–Sham energy functional can be expressed in terms of the density matrix,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle D(x,x') &= \sum_{i=1}^N \phi_i^*(x) \phi_i(x') \nonumber \\ &= \sum_{\mu \nu=0}^{N_b} D_{\mu \nu} \chi^*_{\mu}(x) \chi_{\nu}(x'), \label{density_matrix} \nonumber \end{align} \tag{ 13 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho(x) = D(x,x), \nonumber \end{align} \tag{ 14 }$$

rather than the orbitals, where the density matrix is also of dimension $N_b \times N_b$. The Kohn–Sham energy functional has a closed-form expression in terms of the density matrix (see [28] and section 4.2), and therefore the constrained optimisation in equation (1) becomes an optimisation over allowed variations in the density matrix. From the point of view of the work to follow, the ability to construct, store, and optimise the density matrix directly is the distinguishing characteristic of localised basis sets with respect to the basis set considered in the following work: namely, the set of N_b plane-waves,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi_{G}(x) = {\rm e}^{{\rm i} G.x}, \nonumber \end{align} \tag{ 15 }$$

with the same periodicity as the unit cell [15, 31], labelled by the frequency of the plane-wave G. This basis set is delocalised, meaning the functions $\chi_G$ are non-zero across the whole unit cell. The introduction of a delocalised basis results in a reduction in accuracy per basis function, which in turn necessitates a much larger value of N_b to reproduce the same accuracy as a computation using localised basis sets. The advantage of a plane-wave, or similar, basis set lies elsewhere [15, 31]. This will become relevant in section 4, as certain algorithms exploit the ability to construct $D(x,x')$ explicitly. Nevertheless, much of the analysis to follow in this section will remain largely independent of basis set. The discussion will, however, be framed in the language of an entirely plane-wave basis set.

As already stated, Kohn–Sham theory is a constrained optimisation problem, equation (1). The associated Euler–Lagrange equations provide a method for transforming the optimisation problem into a fixed-point problem: the Kohn–Sham equations. That is, we seek the density $\rho_*$ such that it is a fixed-point of the discretised Kohn–Sham map,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K: \mathbb{R}^{N_b} \rightarrow \mathbb{R}^{N_b}, \nonumber \end{align} \tag{ 16 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K[\rho_*] = \rho_*. \nonumber \end{align} \tag{ 17 }$$

In general, $K[\rho^{{\rm in}}] = \rho^{{\rm out}}$, where K is defined using equations (8) and (9). That is, K takes an input density which is used to construct the Hartree and exchange-correlation potentials, then the associated Kohn–Sham Hamiltonian is diagonalised, and an output density is constructed as the sum of the square of N eigenfunctions. Formally, the Kohn–Sham map is a map from the set of non-interacting $v$-representable densities onto itself. Here, a non-interacting $v$-representable density is a density that can be constructed via equation (9) for a given Kohn–Sham Hamiltonian. The 'size' of this set, as a subset of $\mathbb{R}^{N_b}$, is an open problem [27]. Hence, it is entirely possible that algorithms generate input densities that are not non-interacting $v$-representable; however this appears to not be an issue in practice⁴. The aim now is to generate a converging sequence of densities $\{\rho_0^{{\rm in}}, \rho_1^{{\rm in}},..., \rho_n^{{\rm in}} \}$ starting from an initial guess density $\rho_0^{{\rm in}}$, where $\rho_n^{{\rm in}} \approx \rho_*$ to within some desired tolerance. The ease with which this sequence can be generated in practice depends on the functional properties of K, which are examined later in this section.

The Kohn–Sham map K, can be used to define a new map R, the residual

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle R[\rho_*] = K[\rho_*] - \rho_* = 0, \nonumber \end{align} \tag{ 18 }$$

which transforms the fixed-point problem into a root-finding problem. An absolute scalar measure of convergence is thus provided by the norm of the residual $||R[\rho^{{\rm in}}]||_2$, where ||.||₂ is used to denote the vector L²-norm. However, ||R||₂ is a quantity which lacks transparent physical interpretation, making it difficult to assess just how converged a calculation is by consideration of ||R||₂ alone. Hence, convergence is conventionally defined in terms of fluctuations in the total energy, a more tractable measure. When fluctuations in the total energy are sufficiently low to satisfy the accuracy requirements of the users' calculation, the iterations are terminated and the calculation is converged. In practice, the total energy is often not calculated by evaluating the Kohn–Sham energy functional $E_{\tiny{\rm KS}}[\rho_n^{{\rm in}}]$. Instead, the Harris-Foulkes functional $\widetilde{E}_{\tiny{\rm KS}}$ is defined [38],

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \widetilde{E}_{\tiny{\rm KS}} = \sum_{i=1}^N & \frac{1}{2} \int_{\mathbb{R}^3} \ |\nabla \phi_i |^2 + \frac{1}{2}\int_{\mathbb{R}^3} \rho^\mathrm{out}(x) v_h[\rho^{{\rm in}}] \nonumber \\ & +\frac{1}{2}\int_{\mathbb{R}^3} \left(\rho^\mathrm{out}(x) - \rho^\mathrm{in}(x)\right) v_h[\rho^{{\rm in}}] \nonumber \\ & + \int_{\mathbb{R}^3} \left(\rho^\mathrm{out}(x) - \rho^\mathrm{in}(x)\right) v_{xc}[\rho^{{\rm in}}] +E_{{\rm xc}}[\rho^{{\rm in}}]\nonumber \\ & + \int_{\mathbb{R}^3} \rho^{{\rm out}}(x) v_{{\rm ext}}(x),\nonumber \end{align} \tag{ 19 }$$

which can be shown to give the exact ground state energy correct to quadratic order in the density error about the fixed-point density $\rho_*$—i.e. it is correct to $\mathcal{O}((\rho_* - \delta \rho){}^2)$. Note that it is not the Harris-Foulkes functional that is minimised during the computation, as it possesses incorrect behaviour away from $\rho_*$ [39, 40]. However, evaluating the energy using this functional when near $\rho_*$ allows one to terminate the iterations at a desired accuracy earlier than if one evaluates the energy using the Kohn–Sham functional, which is correct to linear order in the density. Finally, recall that $||R||_2 \rightarrow 0$ is the criterion for solving the Kohn–Sham equations, not for finding a minimum of the Kohn–Sham functional. Indeed, to verify that a local minimiser of the Kohn–Sham functional is obtained, one would need to ensure all eigenvalues of the Hessian were positive. Such a procedure is not practical in plane-wave codes, and hence the exit criterion for algorithms in section 4 is based solely on fluctuations in the total energy.

Identifying properties unique to K can help guide and narrow the choice of algorithms in section 4. Firstly, we note that it is computationally expensive to 'query the oracle', meaning evaluate K for a given input density to generate the pair $\{\rho^{{\rm in}}_i, \rho^{{\rm out}}_i \}$ on the ith iteration. This is because, when one has specified $\rho^{{\rm in}}$, finding the corresponding $\rho^{{\rm out}}$ requires one to construct and diagonalise the Kohn–Sham Hamiltonian. In plane-wave codes, this diagonalisation is done iteratively, and only the relevant N eigenfunctions and eigenvalues are computed. This procedure scales as approximately $\mathcal{O}(N^3)$, and is (in a sense) the bottleneck of the computation [15]. Hence, an algorithm that uses all past iterative data optimally so as to reduce evaluations of K is desirable. Here, the past iterative data constitutes the set of n iterative density pairs $\{(\rho^{{\rm in}}_i, \rho^{{\rm out}}_i) \ | \ i \in [0,n] \}$. In order to utilise this set to generate the subsequent density $\rho^{{\rm in}}_{n+1}$ from some algorithm, one is required to store the history of iterative densities in memory. Each density is represented by an array that scales with N_b⁵, meaning as the iteration number n grows large, so does the memory requirement of storing the entire history. Therefore, a limited memory algorithm is also desirable here, meaning no more than m of the most recent density pairs are stored. The final feature of computational Kohn–Sham theory that we will mention here is the accuracy of the initial guess, $\rho_0^{{\rm in}}$. A discussion on the generation of the initial guess is left to later in this section, but it suffices to note that the initial guess is typically 'close' to the converged density $\rho_*$. By 'close' we mean that a linear response approximation can be employed effectively, see section 4. As perhaps would be expected when this is the case, some of the most successful algorithms are able to utilise the past iterations cleverly with limited memory requirements, and employ some form of linearising approximation.

As mentioned previously, convergence of the self-consistent field iterations depends on the functional properties K, where we recall that each K is specified by the framework of Kohn–Sham theory plus an exchange-correlation approximation and external potential. Despite little being known about the precise functional properties of K [41–43], empirical wisdom allows us to make certain broad statements about it. For the sake of analysis, we now introduce the fixed-point iteration,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{{\rm in}}_{n+1} = \rho^{{\rm out}}_n = K[\rho^{{\rm in}}_n]. \nonumber \end{align} \tag{ 20 }$$

This is perhaps the most simple iterative scheme one could envisage, yet it remains profoundly important from the point of view of functional analysis [44]. An example algorithm that makes use of the fixed-point iteration scheme is given in figure 2. This algorithm, on iteration n, constructs and diagonalises the Kohn–Sham Hamiltonian for a given $\rho^{{\rm in}}_n$, and computes the output density $\rho^{{\rm out}}_n$ from the N eigenvectors corresponding to the lowest N eigenvalues, otherwise known as the aufbau principle. The fixed-point iteration is then used as one sets $\rho^{{\rm in}}_{n+1} = \rho^{{\rm out}}_n$, and the procedure is repeated. For the algorithm in figure 2 to converge, K must be so-called locally k-contractive in the region of the initial guess. For the Kohn–Sham map to be k-contractive under the L²-norm, it must satisfy

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle ||K(\rho_1) - K(\rho_2)||_2 \leqslant k||\rho_1 - \rho_2||_2 \label{lipschitz} \nonumber \end{align} \tag{ 21 }$$

for some real number 0  <  k  <  1. The intuition here is that, for any two points in the 'contractive region', the map K brings these points closer in the L²-norm. A contractive region is simply defined as a domain in density space where equation (21) is true for all densities within this domain. Successive application of K—the fixed-point iteration scheme—thus continues to bring these points closer toward a locally unique fixed-point, $\rho_*$. (See the Banach fixed-point theorem [45] or its generalisations [46] for k-contractive maps). Unfortunately, as section 6 shows, the Kohn–Sham map is not locally k-contractive for the vast majority of Kohn–Sham inputs. However, perhaps surprisingly, certain calculations do lead to a k-contractive Kohn–Sham map, such as spin-independent fcc aluminium at the PBE [47] level of theory, figure 3. In these cases, sophisticated acceleration algorithms tend to do little-to-nothing to assist convergence. The fixed-point iteration is also referred to as the Roothaan iteration in the physics and quantum chemistry communities [48]. It has been demonstrated that, in the context of Hartree–Fock theory, the Roothaan algorithm either converges linearly toward a solution or oscillates between two densities about the solution [49]. It is expected that this behaviour will carry over to Kohn–Sham theory [50].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We now define a new iterative scheme, the damped iteration (or one of its many other aliases, such as Krasnosel'skii–Mann or averaged iteration [51, 52]) such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{{\rm in}}_{n+1} &= \rho^{{\rm in}}_n + \alpha (K[\rho^{{\rm in}}_n] - \rho^{{\rm in}}_n) \nonumber \\ &= \rho^{{\rm in}}_n + \alpha R[\rho^{{\rm in}}_n]. \label{linear_mixing} \nonumber \end{align} \tag{ 22 }$$

Hereafter, we refer to a scheme utilising the damped iteration as linear mixing. This scheme constitutes a series of steps in the residual L²-norm steepest descent direction R weighted by the parameter $\alpha \in (0,1)$. It can be shown that provided K is non-expansive, there always exists some $\alpha$ such that the damped iteration converges [53–55]. Here, non-expansive refers to an instance whereby k  =  1 in equation (21), i.e. densities do not get further apart upon successive application of K. This property is typically assumed, just as we also assume differentiability of $E_{\tiny{\rm KS}}$, in theorems relating to convergence features of algorithms discussed in section 4 (e.g. [56]). Indeed, the past few decades of computation using Kohn–Sham theory has lead to the wisdom that one can always find some $\alpha$ such that one's calculation converges [50], albeit often impractically slowly. Fortunately, rather large damping parameters of $\alpha \sim 0.5$ are sometimes able to significantly improve convergence, as demonstrated in figure 4 [57]. In this sense, the Kohn–Sham map is relatively well-behaved, although many problems of physical interest are not so well-behaved. In these cases, sophisticated algorithms are required in order to accelerate and stabilise convergence. However, as section 6 demonstrates, even when recourse to a sophisticated algorithm is required, most inputs excluding those belonging to certain problematic classes are able to converge effectively. This is a testament to the Kohn–Sham map often being dominated by its linear response within some relatively large region about the current iterate, a property which is examined further in section 3.7.

Zoom In Zoom Out Reset image size

The behaviour of K discussed here could be interpreted as arising due to the lack of convexity of the underlying functional $E_{\tiny{\rm KS}}$ used to generate it. Convexity is defined formally in section 4, but for now it suffices to note that it can be taken to mean $E_{\tiny{\rm KS}}$ has a unique minimum, which is the unique fixed-point of K, and moreover this minimum is global [58]. In other words, solving the Euler–Lagrange equations is a necessary and sufficient condition to verify global optimality. While this is clearly an attractive quality for an energy functional, not least because only the global minimum has direct physical meaning in Kohn–Sham theory, it is not the case here (in general). The lack of convexity of $E_{\tiny{\rm KS}}$ is particularly pronounced in spin-dependent Kohn–Sham theory, where it is not uncommon for many minima to exist, which are interpreted as representing different meta-stable spin states of the system [59]. In this case, one could, for example, employ some form of global optimisation in an attempt to explore the landscape of local minima with hopes of finding the global minimum.

In summary, while a large class of Kohn–Sham inputs are well-behaved and convergent for relatively high values of the damping parameter, many inputs, especially the increasingly complex ones involved in modern technologies, are not. The remainder of this section explores the precise characteristics of K that lead to ill-behaved convergence.

The question remains of how one might go about choosing which N eigenfunctions of the Kohn–Sham Hamiltonian are used to iteratively construct the output densities toward convergence. For N_b  ≫  N, there is of course a large number of permutations of N eigenfunctions from which to choose. While it is perhaps taken for granted, [29] demonstrates that, in the case of Hartree–Fock theory, the lowest energy solution to the Hartree–Fock equations will necessarily be one which corresponds to the N eigenvectors with the lowest eigenvalues of $H_{\tiny{\rm KS}}$. This is otherwise known as the aufbau principle, and appears in the algorithm presented in figure 2. These eigenfunctions $\phi_i$ are termed 'occupied' orbitals, with associated quasi-particle energies $\newcommand{\e}{{\rm e}} \epsilon_i$. However, just because the exact ground state solution satisfies the aufbau principle does not guarantee that doing so at each iteration is optimal [43, 49, 60]. Furthermore, iteratively satisfying the aufbau principle does not guarantee a global, or even local, minimum of $E_{\tiny{\rm KS}}$ will be obtained as a solution to the Kohn–Sham equations [28]. Nevertheless, iteratively satisfying the aufbau pricinple has proven a successful heuristic for finding minima of $E_{\tiny{\rm KS}}$ via the Kohn–Sham equations. Here, the aufbau principle serves to bias our solution of the Kohn–Sham equations toward a minimum of $E_{\tiny{\rm KS}}$, rather than an inflection point or maximum.

Iterative procedures utilising the aufbau principle are well-defined and work best primarily when the input possesses a Kohn–Sham gap, i.e. when it is not a (Kohn–Sham) metal. The Kohn–Sham gap is defined in the limit of large system size as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E_{{\rm gap}} = \epsilon_{N+1} - \epsilon_N, \nonumber \end{align} \tag{ 23 }$$

otherwise known as the HOMO-LUMO gap—the difference in energy between the highest energy occupied and lowest energy unoccupied (molecular) orbitals. When this gap disappears, meaning there exists a non-zero density of states at the Fermi energy, convergence becomes increasingly difficult [25, 26]. Here, the Fermi energy $\mu$ is defined as the energy of the highest occupied orbital, which is the Lagrange multiplier corresponding to the constraint that the particle number N remain fixed in the minimisation of the Kohn–Sham energy functional. Such cases are prone to the phenomenon of occupancy sloshing: iterations become hindered by a continual iterative switching of binary occupation of orbitals whose energies are close to the Fermi energy. In some circumstances, an aufbau solution to the Kohn–Sham equations does not exist for binary occupation of orbitals [27, 61–63]. For example, [61] demonstrates that, in the case of the C2 molecule, the Kohn–Sham solution possesses a 'hole' below the highest occupied orbital. In the context of self-consistent field iterations, this would mean any algorithm would continue to switch orbital occupancies at each iteration ad infinitum. This occurrence is a consequence of degeneracy in the highest occupied Kohn–Sham orbitals, which can occur even in the absence of symmetry and degeneracy in the exact many-body system. Here, and in other cases like this, the density should be constructed from a density matrix

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle D = \sum_{i=1}^q \lambda_{i} \Psi^*_i(x_1,...,x_N) \Psi_i(x'_1,...,x'_N) \label{degen_density_matrix} \nonumber \end{align} \tag{ 24 }$$

via

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho(x) = {{\rm Tr}}D \hat{\rho}(x). \label{ensemble_density} \nonumber \end{align} \tag{ 25 }$$

The wavefunctions $\{\Psi_i \}$ are Slater determinants of Kohn–Sham orbitals corresponding to each degenerate solution within some q-fold degenerate subspace. After rearrangement, we find that the density can now be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho(x) = \sum_{\epsilon_i < \mu} | \phi_i(x) |^2 + \sum_{\epsilon_i = \mu} f_i | \phi_i(x) |^2, \nonumber \end{align} \tag{ 26 }$$

where the fractional occupancies f_i are determined as some combination of the weights $\lambda_i$ in equation (24). This form of the density allows one to see more transparently that we have now introduced fractional occupancy of the orbitals whose energy is degenerate at the Fermi energy. In the example of C2 in [61], the degenerate subspace is first identified, and then the occupancies f _i are varied smoothly until the energies of the identified orbitals are equal. This procedure, termed evaporation of the hole, yielded accurate energy predictions when compared to configuration interaction calculations. In this case, the Kohn–Sham degeneracy is interpreted as being due to the presence of strong electron correlation. These degeneracies lead to densities that are so-called ensemble non-interacting v-representable. That is, the exact Kohn–Sham density can no longer be constructed from a pure state via the sum of the square of orbitals as in equation (3), but instead must be constructed from some ensemble of states via equations (24) and (25). The extension of Kohn–Sham theory to include fractional occupancy is described well in [64, 65].

This so-called ensemble extension to Kohn–Sham theory is also utilised when constructing a non-interacting theory of Mermin's finite temperature formulation of DFT [66]. It is this version of Kohn–Sham theory that is usually used in modern Kohn–Sham codes that include fractional occupancy. As we are interested primarily in how this extension mitigates convergence issues, the reader interested in an in-depth discussion of finite temperature Kohn–Sham theory is referred to [25, 64], and references therein. Here, it suffices to observe that we now seek to minimise the following free energy functional

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E[\{\phi_i \}, \{\,f_i \}, T] = \sum_{i=1}^\infty & \frac{1}{2} f_i \int_{\mathbb{R}^3} \ |\nabla \phi_i |^2 \nonumber \\ & + \frac{1}{2} \int_{\mathbb{R}^3 \times \mathbb{R}^3} \frac{\rho(x)\rho(x')}{|x-x'|} \nonumber \\ &+ \int_{\mathbb{R}^3} \rho(x) v_{{\rm ext}}(x) + E_{{\rm xc}}[\rho] \nonumber \\ &- T S[\{\,f_i \}],\nonumber \end{align} \tag{ 27 }$$

where the entropy functional and density are defined respectively as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S = \sum_{i=1}^\infty f_i \ln (\,f_i) + (1-f_i) \ln (1-f_i), \nonumber \end{align} \tag{ 28 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho(x) = \sum_{i=1}^\infty f_i |\phi_i(x)|^2. \nonumber \end{align} \tag{ 29 }$$

The real-valued fractional occupancies $f_i \in [0,1]$ now constitute additional variational parameters alongside the orbitals. Minimisation of the finite temperature Kohn–Sham functional can be tackled directly as in [26], which is discussed in section 4 and tested in section 6. Alternatively, the associated fixed-point problem can be formulated, whereby the occupancies are given a fixed functional form dependent on both T and the Kohn–Sham Hamiltonian eigenenergies $\newcommand{\e}{{\rm e}} \epsilon_i$. This is otherwise known as the smearing scheme, an example of which is the Fermi–Dirac function,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_i = \frac{1}{{\rm e}^{(\epsilon_i - \mu)/T} + 1}. \nonumber \end{align} \tag{ 30 }$$

The electronic temperature T is now an input parameter which determines the degree of broadening of occupancies about the Fermi energy, figure 5. At each iteration, the occupancies are updated with new values of $\newcommand{\e}{{\rm e}} \epsilon_i$, and this process is continued toward convergence. This procedure demonstrably mitigates occupancy sloshing for Kohn–Sham metals with large density of states at the Fermi energy [67, 68]⁶. Furthermore, introducing finite temperature also assists with sampling of the Brillouin zone in periodic Kohn–Sham codes. That is, interpolation techniques for evaluating integrals across the Brillouin zone are inaccurate when many band crossings (discontinuous changes of occupancy) exist, i.e. in Kohn–Sham metals. This necessitates a fine sampling of k-space in order to accurately evaluate the integrals. As discussed, fractional occupancies negate these discontinuities, allowing for a coarser sampling of the Brillouin zone, meaning interpolation techniques become increasingly accurate—see [67, 69] for more details. In any case, finite electronic temperatures are a valuable numerical tool to assist convergence of the self-consistent field iterations in the event of inputs with large density of states at the Fermi energy. Hence, the test suite in section 5 includes many such systems, and in particular a variety of electronic temperatures are considered.

Zoom In Zoom Out Reset image size

As one might expect, a more accurate initial guess of the variable to be optimised leads monotonically to more efficient and stable convergence rates [70]. In the case of self-consistent field methodology, and for plane-wave and similar codes, the initial guess charge density is often computed as a sum of pseudoatomic densities [15, 31]. That is, once the exchange-correlation and pseudopotential for the atomic species in the computation has been specified, the charge density for these atoms in vacuum is calculated. Then, each individual density is overlaid at positions centered on the atomic cores in order to construct the initial guess density, figure 6. This figure demonstrates visually the accuracy of this prescription for generating initial guess charge densities. Note that different considerations are required in order to generate an initial guess for the density matrix or orbitals. The accuracy of the initial guess is, in part, responsible for the relative success of methods that employ linearising approximations, such as quasi-Newton methods, see section 4. Notable cases in which the initial guess density is relatively poor include polar materials such as magnesium oxide. The initial guess is charge neutral by construction, meaning the charge is required to shift onto the electro-negative species for convergence. Furthermore, inputs whereby the atomic species are subject to large inter-atomic forces can also lead to inaccurate initial guesses. This is partly due to the fact that the initial guess becomes exact in the limit of large atomic separation, and since large inter-atomic forces imply low inter-atomic separation, this can result in potentially inaccurate initial guess densities. Such inputs are generated routinely during structure searching applications [71]. The test suite includes various examples of these 'far-from-equilibrium' systems.

Zoom In Zoom Out Reset image size

Spin polarised Kohn–Sham theory presents more serious issues: there is no widely successful method for generating initial guess spin densities. In the spin polarised or 'unrestricted' formalism, the following spin densities are introduced (see [2]),

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{\uparrow}(x) = \sum_{i=1}^\infty f^{\uparrow}_i|\phi_i^{\uparrow}(x)|^2 \nonumber \end{align} \tag{ 31 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{\downarrow}(x) = \sum_{i=1}^\infty f^{\downarrow}_i|\phi_i^{\downarrow}(x)|^2, \nonumber \end{align} \tag{ 32 }$$

generated from spin up and down particles occupying separate spin orbitals $\{\phi_i^{\uparrow}, \phi_i^{\downarrow} \}$. This leads now to two coupled non-linear eigenvalue problems, one for each spin. A method for generating the initial guess spin densities is thus required, rather than just the initial guess charge density. As one, in general, has no knowledge of the spin state a priori, this initial guess can be relatively far away from the ground state. In practice, one conventionally deals with charge and spin densities, rather than spin up and spin down densities,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{{\rm charge}}(x) = \rho^{\uparrow}(x) + \rho^{\downarrow}(x), \label{charge_density} \nonumber \end{align} \tag{ 33 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{{\rm spin}}(x) = \rho^{\uparrow}(x) - \rho^{\downarrow}(x). \label{spin_density} \nonumber \end{align} \tag{ 34 }$$

The charge density can be initialised similarly to the spin-independent case, with a sum of independent pseudoatomic charge densities. The spin density can be initialised to zero, or be scaled by specifying some magnetic character on the atoms, e.g. ferromagnetic. Such a prescription typically leads to initial guess densities that are further away in the residual L²-norm than spin-independent initial guess charge densities. This observation at least partially accounts for the reason that spin polarised systems tend to be much harder to converge than spin unpolarised systems. For this reason, and others cited in the following section, many spin polarised inputs are included in the test suite. Recently, various schemes have been proposed that aim to better predict self-consistent densities to use as the initial guess [70, 72]. In particular, [70] considers a data-derived approach to predicting and assessing uncertainty in a guess density away from the ground state.

The condition of a problem, loosely speaking, can be taken as characteristic of the difficulty a black-box algorithm will have in solving the problem. Due to the complexity of the Kohn–Sham map, evaluating its condition number directly is impossible in practice. However, within the context of linear response theory, it is possible to explore certain causes of ill-conditioning generic to either all, or certain broad classes, of inputs. Hence, we begin by linearising the map K about a fixed-point⁷,

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \rho_* + \delta \rho^{{\rm out}}_n =& K[\rho_* + \delta \rho^{{\rm in}}_{n}], \nonumber \\ \approx& K[\rho_*] + \frac{\delta K[\rho_n^{{\rm in}}]}{\delta \rho_n^{{\rm in}}}\bigg\vert_{\rho_*} \delta \rho^{{\rm in}}_{n}, \nonumber \\ \Longrightarrow \delta \rho^{{\rm out}}_n = &\frac{\delta K[\rho_n^{{\rm in}}]}{\delta \rho_n^{{\rm in}}}\bigg\vert_{\rho_*} \delta \rho^{{\rm in}}_n. \label{linear_response_approx} \nonumber \end{align} \tag{ 35 }$$

This is the definition of linearisation in the present context, i.e. a small change in the input density yields a change in the output density proportional to the initial change, where the constant of proportionality is shown by the components of the Jacobian of the map K,

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle J_K(x,x') &= \frac{\delta K[\rho_n^{{\rm in}}]}{\delta \rho_n^{{\rm in}}}\bigg\vert_{\rho_*}(x,x') \nonumber \\ &= \frac{\delta \rho^{{\rm out}}_n}{\delta \rho^{{\rm in}}_n}\bigg\vert_{\rho_*}(x,x'). \label{residual_jacobian} \nonumber \end{align} \tag{ 36 }$$

Within the language of linear response theory, the Jacobian can be identified with the non-interacting charge dielectric via

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \varepsilon_0(x,x') = I - J_K(x,x'), \nonumber \end{align} \tag{ 37 }$$

which is the linear response function of the residual map, rather than the Kohn–Sham map. The dielectric can be expanded as such

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \varepsilon_0(x,x') = I - \int_{\mathbb{R}^3} {\rm d}x'' \ \frac{\delta v^{{\rm in}}_{{\rm hxc}}(x'')}{\delta \rho^{{\rm in}}(x')}\frac{\delta \rho^{{\rm out}}(x)}{\delta v^{{\rm in}}_{{\rm hxc}}(x'')}, \nonumber \end{align} \tag{ 38 }$$

where $v_{{\rm hxc}} = v_{\rm h} + v_{\rm xc}$, which are the only two potentials which have a dependence on the density. Hence, the dielectric is given is terms of the non-interacting susceptibility $\chi_0$ as

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \varepsilon_0(x,x') = I - \int {\rm d}x'' \ \big(\,f_{{\rm h}}(x', x'') + f_{{\rm xc}}(x', x'') \big) \chi_0(x, x''), \label{dielectric} \nonumber \end{align} \tag{ 39 }$$

where $f_{{\rm h}}$ and $f_{{\rm xc}}$ are the kernels of the Hartree (Coulomb) and exchange-correlation integrals. Therefore, the linear response of a system to a density perturbation is given by the interplay between the exchange-correlation and Coulomb kernels, and the susceptibility

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi_0 = \frac{\delta \rho^{{\rm out}}}{\delta v^{{\rm in}}_{{\rm hxc}}}, \nonumber \end{align} \tag{ 40 }$$

which is highly system dependent [57, 73]. As the non-interacting susceptibility plays a central role in the description of many physical phenomena, such as absorption spectra, it is a relatively well-studied object [74, 75]. The remainder of this section classifies certain generic behaviours of $\varepsilon_0$ so that causes of divergence in the self-consistency iterations can be studied. First, in order to see why the linear response function is important for self-consistency iterations, note that one may consider each iteration as a perturbation in the density about the current iterate. Knowledge of the exact response function $\chi_0$, and subsequently $\varepsilon_0$, would thus allow one to take a controlled step toward the fixed-point density, depending on how well-behaved⁸ the map is about the current iterate. An iterative scheme utilising the exact response function is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{{\rm in}}_{n+1} = \rho^{{\rm in}}_{n} + \varepsilon_0^{-1} R[\rho^{{\rm in}}_{n}], \nonumber \end{align} \tag{ 41 }$$

which one may recognise as Newton's method. While Newton's method is not global, it has many attractive features, see section 4.1. However, one is rarely privileged with knowledge of the exact dielectric as it is vastly expensive to compute and store [78–81]. In practice, one is left to estimate, or iteratively build, this response function. Cases in which the input is very sensitive to density perturbations, characterised by large eigenvalues of the discretised $\varepsilon_0^{-1}$, tend to amplify errors in iterates, and thus potentially move one away from the fixed-point.

Consider now completely neglecting higher order terms in the Taylor expansion of the Kohn–Sham map, and let us examine the map as if it were linear. This allows us to borrow results from numerical analysis of linear systems, and apply these results as well-motivated heuristics to convergence in the non-linear case. In particular, assuming linearity, absolute convergence can be identified as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &\delta \rho^{{\rm in}}_{n+1} \rightarrow 0 ~{\rm as}~ n \rightarrow \infty, \nonumber \\ \Longrightarrow & (\varepsilon_0^{-1})^n \delta \rho^{{\rm in}}_0 \rightarrow 0 \nonumber \end{align} \tag{ 42 }$$

using equation (35), meaning $\lambda_i < 1$ for all i, where $\lambda_i$ are the eigenvalues of the inverse dielectic matrix, which have been shown to be real and positive for some appropriate $f_{{\rm xc}}$ [57]. Hence, simply multiplying the dielectric by a scalar $\alpha$ such that $\lambda_{{\rm max}}$ is below unity can ensure convergence. This comes at the cost of reducing the efficiency of convergence for components of the density corresponding to low eigenvalues of the dielectric matrix. Defining the condition number of the dielectric as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \kappa = \frac{\lambda_{{\rm max}}(\varepsilon_0^{-1})}{\lambda_{{\rm min}}(\varepsilon_0^{-1})}, \label{cond_number} \nonumber \end{align} \tag{ 43 }$$

it can be seen that the efficiency of the linear mixing procedure is limited by how close this ratio is to unity. One ansatz for the scalar premultiplying the dielectric is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha = \frac{2}{\lambda_{{\rm max}} + \lambda_{{\rm min}}}, \nonumber \end{align} \tag{ 44 }$$

which ensures, as much as the linear approximation is valid, that components of the density corresponding to the maximal and minimal eigenvalues of $\varepsilon_0^{-1}$ converge at the same rate [57, 73, 82]. However, this form of $\alpha$ ignores the distribution (e.g. clustering) of eigenvalues [83], and is not commonly used in conjunction with more sophisticated schemes such as those in section 4. An additional strategy to improve convergence would be to construct a matrix, the preconditioner, such that when the preconditioner is applied to $\varepsilon_0$, the eigenspectrum of the product is compressed toward unity. This is done in practice, see [15] for example, and is the core idea behind the Kerker preconditioner [84, 85], as discussed shortly.

It is clear from equation (43) that the convergence depends critically on the spectrum of the inverse dielectric matrix. The minimum eigenvalue is unity, and the large eigenvalues are dominated by the contributions from the Coulomb kernel, rather than the exchange-correlation kernel [73, 78, 86]. To see why this is, it is first asserted that the discretised Coulomb kernel possesses large (divergent) eigenvalues due to its $|x-x'|^{-1}$ dependence, which will be demonstrated in the work to follow shortly. These large eigenvalues act as a generic amplification factor when multiplied by the eigenvalues of $\chi_0$ in equation (39), thus resulting in large (divergent) eigenvalues in the dielectric. Conversely, the exchange-correlation kernel does not generically possess such large eigenvalues. In semi-local Kohn–Sham theory, the exchange-correlation kernel is a polynomial of the density and potentially its higher order derivatives. Crucially, it has a simple, local (explicit) dependence on position given by $\delta(x-x')$,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{{\rm xc}}(x,x') = \delta(x-x') \frac{\delta v_{\rm xc}}{\delta \rho}, \nonumber \end{align} \tag{ 45 }$$

that cannot, by definition, introduce x-dependent ill-conditioning. As such, no generic amplification of the eigenvalues of $\chi_0$ occurs, and hence the exchange-correlation kernel can be ignored relative to the Coulomb kernel from the perspective of ill-conditioning. In other words, the following analysis utilises in the random phase approximation (RPA) by setting $f_{{\rm xc}} = 0$ in equation (39). As [78] notes, even in situations whereby the density vanishes in some region, meaning that negative powers of the density are divergent, the linear response function tempers this divergence, and the exchange-correlation contribution remains well-conditioned.

The principle categorisation one can make when analysing generic behaviour of the response function is the distinction between Kohn–Sham metals and insulators. Consider a homogeneous and isotropic system, i.e. the homogeneous electron gas, such that $\chi_0(x,x') \rightarrow \chi_0(|x-x'|)$, which satisfies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \delta \rho^{{\rm out}} (x) = \int {\rm d}x' {{\rm ~}} \chi_0(|x-x'|) \delta v^{{\rm in}}_{{\rm h}}(x'). \label{nonlocal_suscept} \nonumber \end{align} \tag{ 46 }$$

This is a convolution in real space, and hence a product in reciprocal space

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \delta \tilde{\rho}^{{\rm out}} (G) = \tilde{\chi}_0(|G|) \delta \tilde{v}^{{\rm in}}_{{\rm h}}(G), \nonumber \end{align} \tag{ 47 }$$

where we label the Fourier components G by convention. This susceptibility is local and homogeneous in reciprocal space, and relates perturbations in the input density to a response by the output density (within the RPA) via

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \delta \tilde{\rho}^{{\rm out}} (G) &= f_{{\rm h}}(|G|) \tilde{\chi}_0(|G|) \delta \tilde{\rho}^{{\rm in}} (G) \nonumber \\ &= \frac{4 \pi \tilde{\chi}_0(|G|)}{|G|^2} \delta \tilde{\rho}^{{\rm in}} (G). \nonumber \end{align} \tag{ 48 }$$

The susceptibility of the homogeneous electron gas, which constitutes a simple metal in the present context, is derived from Thomas–Fermi theory as the Thomas–Fermi wavevector $\tilde{\chi}_0 \sim k^2_{\tiny{\rm TF}}$, which is constant⁹ [87]. It can therefore be seen that if there is any error in a trial input density, generated by an iterative algorithm, away from the optimal update, then this error is amplified by a factor of |G|⁻² for $|G|<1$, where $|G|=0$ does not contribute. This sensitivity to error in iterates is identified as the source of charge sloshing, and is a somewhat generic feature of Kohn–Sham metals. Whilst the above derivation utilises Thomas–Fermi theory of the homogeneous electron gas to demonstrate constant susceptibility, it can be shown that all Kohn–Sham metals display this behaviour in the small $|G|$ limit [75, 88]. A demonstration of charge sloshing is illustrated in figure 7, whereby a linear mixing algorithm purposefully takes slightly too large steps in the density. This leads to vast over-corrections in each iteration, giving the appearance that charge is 'sloshing' about the unit cell. This is not the only source of large eigenvalues of the dielectric in Kohn–Sham metals, as the susceptibility possesses inherently divergent eigenvalues independent from the amplification by the Coulomb kernel. To see this, consider the Adler–Wiser equation which is defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi_0(x,x') = \sum_{n=1}^{N} \sum_{m=N + 1}^\infty \frac{\phi_n(x) \phi^*_m(x) \phi^*_n(x') \phi_m(x')}{\epsilon_n - \epsilon_m}, \label{Adler--Wiser} \nonumber \end{align} \tag{ 49 }$$

which is an expression from perturbation theory for the exact Kohn–Sham susceptibility [74, 75]. As [82] originally noted, the denominator $\newcommand{\e}{{\rm e}} \epsilon_n - \epsilon_m$ approaches zero when the input is gapless, i.e. it has a large density of states about the Fermi energy. If left untreated, this observation, in conjunction with the amplifying factor from the low $|G|$ components of the Coulomb kernel, lead to significant ill-conditioning. The largest condition numbers arise when $|G|$ is extremely small; since G is a reciprocal lattice vector, this will occur for unit cells that are large in any (or all) of the three real-space dimensions. Whilst the dependence of the eigenvalues of the dielectric on unit cell size is in practice complicated [73], it suffices to note that increased unit cell size is a significant source of ill-conditioning. As compute power continues to grow, larger and larger systems are being tackled using Kohn–Sham theory, and the increase in required number of self-consistency iterations as a result of this instability poses serious issues for Kohn–Sham calculations. Inefficiencies of this kind are best dealt with using preconditioners, as section 4 demonstrates. On the other hand, insulators possess no such divergences in the eigenvalues of the dielectric matrix. It can be shown that in the low $|G|$ limit the behaviour of the susceptibility for gapped materials is [75, 88]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{\chi}_0 \propto |G|^2. \nonumber \end{align} \tag{ 50 }$$

Zoom In Zoom Out Reset image size

This functional dependence cancels the |G|⁻² dependence from the Coulomb kernel, and thus the eigenvalues of the dielectric become constant for all G. This constant is unknown in general, and guaranteed convergence for simple insulators amounts to finding the damping parameter $\alpha$ such that this constant is below unity. This is in line with the empirical wisdom that insulators are much easier to converge than metals, provided that the insulator does not artificially assume a metallic character during the self-consistency iterations.

In this vein, inputs that are increasingly complex, i.e. deviating from simple metals or insulators, are likely to exhibit problematic behaviour. As discussed in section 4, preconditioners are able to alleviate charge sloshing in simple metals. However, when a metal-insulator interface is used as input, there are regions with starkly different behaviour in the response function, which is difficult to capture analytically. Hence, preconditioning techniques may fail to assist, and even hinder, iterations in calculations on interfaces of this kind [73]. Furthermore, it is possible that artificial phase transitions between gapped and gapless phases occur during the self-consistency iterations. Many algorithms function by building up an approximation to the dielectric using past iterates. The discontinuous change in behaviour of the dielectric in differing phases causes parts of the iterative history to actively interfere in correctly modelling the dielectric. Hence, iterations become hindered or divergent. An artificial phase change of this kind is demonstrated to occur in [89] for an isolated iron atom. Various examples of the aforementioned problematic classes of inputs are included in the test suite.

Finally, a brief comment is provided on how the above analysis translates to spin-dependent Kohn–Sham theory. As discussed, in the spin-dependent case one solves two non-linear eigenvalue problems that independently look very similar to equations (8) and (9), but crucially are coupled through the Hartree and exchange-correlation potential. That is, an algorithm that perturbs the spin up (spin down) density will lead to a response by the spin up (spin down) density given by the prior analysis. However, one must now also consider how a perturbation in the spin up density affects the spin down density, and vice versa, which is again through the Coulomb and exchange-correlation kernels. Hence, all of the above sources of ill-conditioning translate directly to the spin-dependent case, with the added difficulty that the number of optimisation parameters has doubled, and these parameters are coupled in such a way that potentially introduces further ill-conditioning. In other words, the response function now contains four components,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \varepsilon_0 = \left[\begin{array}{@{}cc@{}} \varepsilon^{\uparrow \uparrow}_0 & \varepsilon^{\uparrow \downarrow}_0 \nonumber \\ \varepsilon^{\downarrow \uparrow}_0 & \varepsilon^{\downarrow \downarrow}_0 \end{array}\right], \nonumber \end{align} \tag{ 51 }$$

where the diagonal elements describe the response of the spin up (spin down) density to perturbations in the spin up (spin down) density, and the off-diagonal elements describe the response of the spin up (spin down) density to perturbations in the spin down (spin up) density. Implementations of spin-dependent Kohn–Sham theory often utilise the charge and spin densities, defined in equations (33) and (34), rather than the spin up and spin down densities. In this representation, the behaviour of the off-diagonal (coupling) elements is determined via the exchange-correlation potential, as the Hartree contribution cancels, see [2, 57]. To the authors' knowledge, there is less literature on the manifestation of this coupling in the self-consistency iterations than on the spin-independent counterpart. Dederichs and Zeller [57] uses self-consistency in the Stoner model to demonstrate that the condition of the system is indeed worsened in the presence of magnetism due to the coupling. However, it is noted that the charge and spin densities decouple near self-consistency. In any case, for these reasons, and perhaps for reasons yet unexplored, empiricism demonstrates that spin polarised calculations are, in general, more difficult to converge than spin unpolarised calculations. Analysing this coupling in terms of self-consistency and ill-conditioning is an avenue for further research, which could be used to assist self-consistency iterations through the use of improved preconditioners, see section 4.3.

First, we make some general comments about the Newton update. The Newton update is the optimal first order update in the density at the current iteration. In other words, if the current iterate is within the linear response radius of the root, then the exact Newton update would lead to convergence in one iteration by definition. That is, we seek the update $\delta \rho_n$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle R[\rho^{{\rm in}}_n + \delta \rho_n] \approx R[\rho^{{\rm in}}_n] + J|_{\rho^{{\rm in}}_n} \delta \rho_n = 0, \nonumber \end{align} \tag{ 58 }$$

where J:  =  J_R is the Jacobian of the residual, defined similarly to equation (36), evaluated at the current iterate. Rearranging for $\delta \rho_n$, the update is given as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{{\rm in}}_{n+1} = \rho^{{\rm in}}_n - J|_{\rho^{{\rm in}}_n}^{-1} R[\rho^{{\rm in}}_n]. \nonumber \end{align} \tag{ 59 }$$

Assuming the Jacobian exists and is Lipschitz continuous¹⁰, this update is shown to have quadratic convergence in some region about the root [83]. The Jacobian must be computed numerically, which can be done with either the Adler–Wiser equation equation (49) [74, 75], or with finite-difference numerical differentiation [94]. As section 4.3 will explore in more depth, in the absence of further approximation, both of these techniques are inadequate for modern calculations due to the computational complexity and the size of the basis set. The former strategy is an $O(N_b^4)$ process that requires the computation and storage of all eigenvectors of the Kohn–Sham Hamiltonian [78, 79]. The latter strategy requires excessively many evaluations of K [95]. We now examine the class of methods that can be cast as a Newton step with some iteratively updated approximation to the Jacobian: quasi-Newton methods.

4.1.1. Broyden's methods.

Consider having knowledge of an approximate Jacobian at the previous iteration, J_n−1. We seek a prescription for generating an approximate Jacobian at the current iteration, J_n, such that the following quasi-Newton update can be performed,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{{\rm in}}_{n+1} = \rho^{{\rm in}}_n - J^{-1}_n R_n, \label{QN_step} \nonumber \end{align} \tag{ 60 }$$

where $R_n := R[\rho^{{\rm in}}_n]$. First, note that all methods of this kind must begin from some initial estimate of the Jacobian, J₀. For lack of a better option, this can be taken as a scaled identity, $J_0 = \alpha I$. Although, in the present context, the Kerker matrix is used, which is defined in section 4.3. We begin with a description of Broyden's two methods [96]. These methods, as they are about to be presented, are not commonly used in modern Kohn–Sham software. However, the conceptual foundation of Broyden's methods, that is, low rank updates to a Jacobian that satisfies secant conditions, remain foundational to contemporary methodology. First, the meaning of a secant condition is defined. For illustrative purposes, a finite-difference approximation for the derivative of a one-dimensional function f  at the current iterate x_n is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f'(x_n) = \frac{f(x_n) - f(x_{n-1})}{x_n - x_{n-1}}, \nonumber \end{align} \tag{ 61 }$$

which is increasingly accurate as the iterates become closer. Since the Jacobian is the derivative of the residual map, the N_b-dimensional equivalent of this finite-difference equation is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J_n (\rho^{{\rm in}}_n - \rho^{{\rm in}}_{n-1}) = R_n - R_{n-1}, \label{secant_condition} \nonumber \end{align} \tag{ 62 }$$

where hereafter we define $\Delta \rho^{{\rm in}}_n = \rho^{{\rm in}}_n - \rho^{{\rm in}}_{n-1}$ and $\Delta R_n = R_n - R_{n-1}$. If the nth Jacobian satisfies equation (62), it is said to satisfy the secant condition of the current iteration, and thus belongs to Broyden's family of methods. Since J_n is an $N_b \times N_b$ matrix, and the secant condition only specifies how J_n acts on the vector $\Delta \rho^{{\rm in}}_n$, there are a remaining $N_b^2 - N_b$ components of the Jacobian that are yet unspecified. Broyden fixes these remaining components by requiring J_n acts on all vectors orthogonal to $\Delta \rho^{{\rm in}}_n$ similarly to J_n−1. This is equivalent to requiring that the Jacobian of the current iteration solves the following constrained quadratic programming problem,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm minimise}~|| J_n - J_{n-1} ||_f \label{broyden_1} \nonumber \end{align} \tag{ 63 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm subject~to}~ J_n \Delta \rho_n^{{\rm in}} = \Delta R_n, \nonumber \end{align} \tag{ 64 }$$

as demonstrated by [97, 98], which defines Broyden's first method. The Frobenius norm ||.||_f of a square matrix A is defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle ||A||_f = \sqrt{\sum_{i,j=1}^N |a_{ij}|^2}. \nonumber \end{align} \tag{ 65 }$$

In other words, the current Jacobian J_n is required to satisfy the current secant condition, and otherwise minimise the difference between itself and the previous Jacobian J_n−1 in the sense of the Frobenius norm. Note that the nth Jacobian satisfies all of the previous n secant equations provided the past iterates are mutually orthogonal, $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}(\Delta \rho^{{\rm in}}_i){}^{\dagger} \Delta \rho^{{\rm in}}_j = 0$ for $i \neq j$. However, the space of past iterates is often linearly independent, but not mutually orthogonal. Therefore, if one requires J_n to satisfy only the most recent secant equation, one loses information about past secant equations, i.e. J_n no longer satisfies the past secant equations. Schemes that ensure J_n satisfies multiple previous secant equations are studied in the next section.

The constrained optimsiation problem of equation (63) has a unique analytic solution, which is obtained in [97, 98] by means of unconstrained optimisation using the method of Lagrange multipliers,

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle J_n = J_{n-1} + \frac{\Delta R_n - J_{n-1} \Delta \rho_n^{{\rm in}}}{||\Delta \rho^{{\rm in}}_n||_2^2} (\Delta \rho^{{\rm in}}_n)^{\dagger}. \label{b1_update} \nonumber \end{align} \tag{ 66 }$$

The notation $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}uv^{\dagger}$ defines the outer product of the vectors $u,v$. One can now observe from equation (66) that this prescription has lead transparently to a rank-one update of the Jacobian at each iteration. The full quasi-Newton update for Broyden's first method involves subsequently inverting equation (66), applying it to residual vector, and performing the quasi-Newton step equation (60). The apparent excessive cost of inverting equation (66) is negated as the inverse of a rank-one matrix can be computed analytically using the Sherman–Morrison–Woodbury formula [99]. Furthermore, as matrix-vector multiplication is associative, one can compute the vector $J_n^{-1}R_n$ without constructing or storing J_n explicitly, and instead using a series of vector–vector products. This was originally demonstrated in [100], so that at a given instance Broyden's first method only requires the storage of two N_b-length vectors, and the computation of a few vector–vector products. Broyden's second method optimises the components of the matrix $H_n := J_n^{-1}$ directly via

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm minimise}~ || H_n - H_{n-1} ||_f \label{broyden_2} \nonumber \end{align} \tag{ 67 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm subject~to} ~ \Delta \rho_n^{{\rm in}} = H_n \Delta R_n, \label{b2_constraint} \nonumber \end{align} \tag{ 68 }$$

instead of optimising the Jacobian, then subsequently inverting. Hereafter, methods that optimise the Jacobian are referred to as 'type-I' methods, and methods that optimise the inverse Jacobian are referred to as 'type-II' methods, see [91]. Note that the constraint in equation (68) is simply the inverse secant condition. Similarly to Broyden's first method, this has the analytic solution,

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle H_n = H_{n-1} + \frac{\Delta \rho_n^{{\rm in}} - H_{n-1} \Delta R_n }{||\Delta R_n||_2^2} (\Delta R_n)^{\dagger}, \label{b2_update} \nonumber \end{align} \tag{ 69 }$$

which can be substituted directly into the quasi-Newton step¹¹. The conventional wisdom has emerged that Broyden's second method tends to provide more robust and efficient convergence than Broyden's first method. However, both methods are shown to be q  −  superlinearly convergent [83, 97] in the sense that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \lim_{n \rightarrow \infty} \frac{||\rho^{{\rm in}}_{n+1} - \rho_*||_2}{||\rho^{{\rm in}}_{n} - \rho_*||_2} = 0, \nonumber \end{align} \tag{ 70 }$$

which is a necessary condition for some q  >  1 in equation (56). Broyden's second method is implemented and tested in section 6.

4.1.2. Multisecant Broyden's methods.

A natural extension to Broyden's methods is to consider all prior secant conditions at each iteration, rather than just the most recent secant condition. This leads to a so-called generalised or multisecant version of Broyden's methods, which are examined extensively in both optimisation and electronic structure literature [83, 91, 100, 102]. The ensuing summary follows a similar structure to that of [91]. A multisecant method is defined as a method that generates an iterative Jacobian J_n such that this Jacobian satisfies the most recent m secant conditions. That is, the following $N_b \times m$ matrices are defined

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta \mathcal{R}_n := (\Delta R_{n-m+1}, \Delta R_{n-m+2}, ..., \Delta R_n) \nonumber \end{align} \tag{ 71 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta \mathcal{P}_n := (\Delta \rho^{{\rm in}}_{n-m+1}, \Delta \rho^{{\rm in}}_{n-m+2}, ..., \Delta \rho^{{\rm in}}_n) \nonumber \end{align} \tag{ 72 }$$

such that a Jacobian satisfying the previous m secant conditions must satisfy the matrix equation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J_n \Delta \mathcal{P}_n = \Delta \mathcal{R}_n. \label{multisecant_condition} \nonumber \end{align} \tag{ 73 }$$

The parameter m introduced here defines the history length, i.e. the number of iterates that are stored and used for secant conditions. If m is less than the full history size n then the method takes on its modified limited memory form. If m  =  n, then the method satisfies all prior secant conditions. The generalisation of Broyden's two methods is now readily established: alter the constraints in the optimisation problems equations (63) and (67) to reflect the multisecant condition equation (73). The multisecant version of Broyden's first and second method respectively are

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm minimise}~|| J_n - J_{n-1} ||_f \label{MSB1} \nonumber \\ {\rm subject~to} ~J_n \Delta \mathcal{P}_n = \Delta \mathcal{R}_n, \nonumber \end{align} \tag{ 74 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm minimise}~|| H_n - H_{n-1} ||_f \label{MSB2} \nonumber \\ {\rm subject~to}~\Delta \mathcal{P}_n = H_n\Delta \mathcal{R}_n, \nonumber \end{align} \tag{ 75 }$$

which are of type-I and type-II respectively. These both have a unique analytic solution in the form of a rank-m update,

$$\begin{align*} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle &J_n = J_{n-1} + (\Delta \mathcal{R}_n - J_{n-1} \Delta \mathcal{P}_n)(\Delta \mathcal{P}^{\dagger}_n \Delta \mathcal{P}_n)^{-1} \Delta \mathcal{P}_n^{\dagger}, \nonumber \\ &H_n = H_{n-1} + (\Delta \mathcal{P}_n - H_{n-1} \Delta \mathcal{R}_n)(\Delta \mathcal{R}_n^{\dagger} \Delta \mathcal{R}_n)^{-1} \Delta \mathcal{R}_n^{\dagger}, \nonumber \end{align*}$$

which are found by solving the associated Lagrangian problems. The former Jacobian update can be inverted similarly to Broyden's first method with the Sherman–Morrison–Woodbury formula. As [89, 92] conclude, and section 6 also examines, the type-II variant tends to outperform the type-I variant in the context of multisecant Broyden's methods, in line with the conventional wisdom from Broyden's original methods. As stated previously, if the space of past iterates is mutually orthogonal, this method is equivalent to Broyden's original methods.

Finally, we remark on the connection between the above methods and the method examined by Eyert, Vanderbilt & Louie and Johnson in [101–103]. First, the following unconstrained minimisation problem for variations in H_n is defined,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm minimise}~ &w_0 || H_n - H_{n-1} ||_f \nonumber \\ &+ \sum_{i=n-m+1}^n w_i || \Delta \mathcal{P}_n - H_n\Delta \mathcal{R}_n ||_2^2, \label{VL} \nonumber \end{align} \tag{ 76 }$$

where we choose to update the inverse Jacobian H_n−1, although a similar method can be formulated in terms of Jacobian updates. The weights $\{w_0,w_i \}$ are introduced as free parameters that act as penalty coefficients. That is, the weights are chosen to signify how 'important' it is to satisfy the corresponding constraint. In this sense, inspection of equation (76) shows that w₀ controls the degree to which the inverse Jacobian can change iteration-to-iteration, and w_i controls the degree to which the ith secant equation should be satisfied by H_n. Therefore, this method also constitutes a multisecant method, but the multisecant conditions are allowed to be weighted according to relative importance. Various common fixed-point methods can be recovered as special cases of these weights. Notably, as [15, 31] demonstrate, the choice w_i  =  0 for i  <  n, and $w_0 << w_n$, leads to Broyden's second method. This can be intuited from equation (76): the weights w_i now favour exclusively the most recent secant condition, and in directions orthogonal to that secant condition, the minimum norm condition on H_n is applied. In the original work of [102, 103], the weights $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}w_i = (R_i^{\dagger} R_i){}^{-1}$ are considered, which favour secant conditions closer to convergence. This was used in the context of electronic structure calculations with success in [101–103]. However, as [101] demonstrates, the optimal set of weights require $w_0 \rightarrow 0$, and if {w_i} are to be non-zero, these weights in fact cancel in the update formula. Hence, w_i  =  1 can be set without loss of generality, and the method can be identified with a standard multisecant method; see [101] for additional detail. An interesting aspect of the multisecant methods discussed here are their relationship Pulay's or Anderson's method—a ubiquitous method in electronic structure theory software—which is now examined.

4.1.3. Pulay's method.

Pulay's method [104, 105], or the discrete inversion in the iterative subspace (DIIS), as it is known in electronic structure literature, or Anderson's method, as it is known in optimisation literature [106], has proven extremely effective at converging Kohn–Sham calculations. The simplicity of its formulation combined with its impressive efficiency and robustness has lead to Pulay's method becoming the default algorithm in a range of Kohn–Sham codes [10, 15, 33, 34, 107, 108]. The past few decades of wisdom suggest that Pulay's method systematically outperforms the unmodified Broyden's methods in both the single and multisecant formulation. This conclusion will be tested in section 6. First, a brief review of Pulay's method as it was originally formulated is given.

Consider constructing a so-called 'optimum' residual—a residual whose argument is an optimum density—as a linear combination of past residuals in the m-dimensional iterative subspace,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle R[\rho^{{\rm in}}_{\rm opt}] = \sum_{i=n-m+1}^n c_i R[\rho^{{\rm in}}_i]. \label{opt_residual_1} \nonumber \end{align} \tag{ 77 }$$

Here, optimum is defined by the method one chooses to fix the coefficients c_i. In Pulay's method, these coefficients are fixed by requiring that the L²-norm of the residual is minimal, i.e. solve

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm minimise}~& ||R[\rho^{{\rm in}}_{\rm opt}]||_2^2 \label{pulay_statement} \nonumber \\ {\rm subject~to} ~&\sum_i c_i = 1, \nonumber \end{align} \tag{ 78 }$$

where the constraint that the coefficients must sum to unity is an exact requirement at convergence. Substitution of equations (77) into (78), and use of Lagrange multipliers, allows the optimisation problem to be cast as an $(n+1)$-dimensional linear system,

$$\begin{align*} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \left(\begin{array}{@{}ccccc@{}} R_1^{\dagger}R_1 & R_1^{\dagger}R_2 & \dots & R_1^{\dagger}R_n & 1 \nonumber \\ R_2^{\dagger}R_1 & \ddots & & & 1 \nonumber \\ \vdots & & & & \vdots \nonumber \\ R_n^{\dagger}R_1 & & & & \nonumber \\ 1 & 1 & \dots & \dots & 0 \end{array}\right) \left(\begin{array}{@{}c@{}} c_1 \nonumber \\ c_2 \nonumber \\ c_3 \nonumber \\ \vdots \nonumber \\ c_n \nonumber \\ \lambda \end{array}\right) = \left(\begin{array}{@{}c@{}} 0 \nonumber \\ 0 \nonumber \\ 0 \nonumber \\ \vdots \nonumber \\ 0 \nonumber \\ 1 \end{array}\right), \nonumber \end{align*}$$

for $n \leqslant m$, which is readily generalised for n  >  m. Assuming the space of past iterates is of full rank (comprised of linearly independent vectors), solution of this linear system provides the set of coefficients c_i. Given these coefficients, the density update remains to be defined. Following [56, 90, 91, 101, 109], the optimum residual can be first be expanded as such,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle R[\rho^{{\rm in}}_{\rm opt}] = \sum_{i=n-m+1}^n c_i K[\rho^{{\rm in}}_i] - \sum_{i=n-m+1}^n c_i \rho^{{\rm in}}_i. \label{opt_residual} \nonumber \end{align} \tag{ 79 }$$

If K is assumed to be linear, the rightmost term in equation (79) can be interpreted as the optimal input density,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{{\rm in}}_{{\rm opt}} = \sum_{i=n-m+1}^n c_i \rho^{{\rm in}}_i. \nonumber \end{align} \tag{ 80 }$$

Hence, the optimal update can take the standard undamped form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{{\rm in}}_{n+1} = \rho^{{\rm in}}_{{\rm opt}} + R[\rho^{{\rm in}}_{\rm opt}] \nonumber \end{align} \tag{ 81 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle = \sum_{i=n-m+1}^n c_i K[\rho^{{\rm in}}_i]. \nonumber \end{align} \tag{ 82 }$$

This update is favoured over $\rho^{{\rm in}}_{n+1} = \rho^{{\rm in}}_{{\rm opt}}$ so that the algorithm does not stagnate in the subspace of past input densities. Alternatively, as originally studied in [106], a damped step can be taken,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{{\rm in}}_{n+1} = \rho^{{\rm in}}_{{\rm opt}} + \alpha R[\rho^{{\rm in}}_{{\rm opt}}], \label{damped_pulay_update} \nonumber \end{align} \tag{ 83 }$$

for $\alpha \in (0,1)$. An example algorithm that implements this formulation of Pulay's method is given in algorithm 1.

Algorithm 1. Pulay's algorithm.

1: Input: tol, $\rho_0^{{\rm in}}, m, \alpha$

2: for $n=0,1,2$,...

3:   Compute $R[\rho^{{\rm in}}_n]$ and store the pair $\{\rho^{{\rm in}}_n, R[\rho^{{\rm in}}_n] \}$

4:   Solve equation (78) for {ci}

5:   if $n \leqslant m$

6:     Set $\rho^{{\rm in}}_{n+1} = \sum_{i=1}^n c_i (\rho^{{\rm in}}_i + \alpha R[\rho^{{\rm in}}_i])$

7:   else

8:     Set $\rho^{{\rm in}}_{n+1} = \sum_{i=n-m+1}^n c_i (\rho^{{\rm in}}_i + \alpha R[\rho^{{\rm in}}_i])$

9:   if $||R[\rho^{{\rm in}}_n]||_2 <$ tol, exit

At a first glance, Pulay’s method bears little resemblance to the secant-based methods discussed in the previous section. However, as described in [56, 91, 109], a rearrangement of the optimisation problem in equation (78) reveals a close relationship between Pulay’s method and type-II multisecant methods. A more detailed treatment of this correspondence is found in [56, 91, 109]; here, we simply state that the following unconstrained optimisation problem

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm minimise}~ ||R_n - \Delta \mathcal{R}_n \gamma ||_2, \nonumber \end{align} \tag{ 84 }$$

for variations in $\gamma = (\gamma_{n-m+1},\gamma_{n-m+2},...,\gamma_n)$ is equivalent to Pulay's optimisation problem equation (78). The new coefficients $\{\gamma_i \}$ are related to the old coefficients {c_i} such that the update in equation (83) now takes the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{{\rm in}}_{n+1} = \rho^{{\rm in}}_n + \alpha R_n - (\Delta \mathcal{P}_n + \alpha \Delta \mathcal{R}_n) \gamma, \nonumber \end{align} \tag{ 85 }$$

where $\gamma$ on iteration n is solved by

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \gamma = (\Delta \mathcal{R}_n^{\dagger} \Delta \mathcal{R}_n)^{-1}\Delta \mathcal{R}_n^{\dagger} R_n. \nonumber \end{align} \tag{ 86 }$$

The parallel between Pulay's method and multisecant methods becomes apparent when these equations are combined to give the final update,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{{\rm in}}_{n+1} = \rho^{{\rm in}}_n + H_n R_n, \nonumber \end{align} \tag{ 87 }$$

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle H_n = \alpha I - (\Delta \mathcal{P}_n - \alpha \Delta \mathcal{R}_n)(\Delta \mathcal{R}_n^{\dagger} \Delta \mathcal{R}_n)^{-1} \Delta \mathcal{R}_n^{\dagger}. \label{pulay_update_v2} \nonumber \end{align} \tag{ 88 }$$

By comparison with the updated inverse Jacobian in equation (75), we can observe that Pulay's method is a type-II quasi-Newton step where the iterative inverse Jacobian is updated according to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm minimise} ~&|| H_n - H_0 ||_f \label{pulay_opt_II}\nonumber \\ {\rm subject~to} ~&\Delta \mathcal{P}_n = H_n\Delta \mathcal{R}_n. \nonumber \end{align} \tag{ 89 }$$

In other words, this optimisation problem minimises the difference between the components of the inverse Jacobian H_n and the initial guess inverse Jacobian H₀, while also requiring the previous m secant conditions to be fulfilled. Note that $H_0= \alpha I$ is required in order to recover the update in equation (88). This reformulation not only connects Pulay's method to the type-II variant of multisecant Broyden's method, but also uncovers another flavour of Pulay's method,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm minimise} ~&|| J_n - J_0 ||_f \label{pulay_opt_I}\nonumber \\ {\rm subject~to} ~& J_n\Delta \mathcal{P}_n = \Delta \mathcal{R}_n, \nonumber \end{align} \tag{ 90 }$$

which is of type-I, i.e. the Jacobian is optimised, rather than the inverse Jacobian. This form of Pulay's method, as originally described in [91], has seen comparatively less application and testing in the context of Kohn–Sham codes [56]. These methods differ from multisecant Broyden methods precisely when m  >  n, in which case the multisecant Broyden methods retain information from all prior secant equations implicitly, whereas Pulay's method(s) ignore completely secant equations not in the (size m) history.

A few modifications to Pulay's method are now examined; although we note that these modifications are adaptable to all the secant-based methods discussed previously. First, the work in [110], based on [111], suggests that alternating between Pulay and linear mixing steps can improve the robustnesss of the iterations over standard Pulay—the 'Periodic Pulay' method. Each new Pulay step utilises the history from the linear mixing and past Pulay steps to solve the optimisation subproblem equation (78). This is demonstrated to have a stabilising effect as the linear mixing history data is used well by the Pulay extrapolation. In [110], an input parameter k determines the number of linear mixing steps performed between each Pulay step, i.e. k linear steps per Pulay step. As suggested in the original work, the values k  =  2 are tested in section 6 with a damping parameter of $\alpha=0.2$ for both the Pulay and linear mixing steps.

Second, [112] considers occasionally flushing the history every time a certain criterion is met, rather than iteratively overriding the history—'Restarted Pulay'. This criterion is chosen to be whenever the current iteration number is an integer multiple of the maximum history size, i.e. n  =  am for some $a \in \mathbb{Z}^+$. For inputs with a considerable degree of non-linearity, either due to a poor initial guess, or inherent to the Kohn–Sham map, the history can actively interfere with modelling an accurate iterative Jacobian at the current iteration. Restarted Pulay thus represents a strategy for dealing with this issue by periodically removing the history.

The final technique we discuss here is the 'Guaranteed Reduction Pulay' algorithm of [113]. The approach of Guaranteed Reduction Pulay involves ensuring that the Pulay predicted optimal residual $||R_{{\rm opt}}||_2$ decreases each iteration. This is achieved by rearranging the stored history of residuals {R_i} such that, at a given iteration, the Pulay predicted optimal residual is added to the history, rather than the residual obtained from evaluating the Kohn–Sham map. The subsequent iteration then involves a linear mixing step, which generates a new exact residual that is added to the history. The coefficients {c_i} of the now current iteration are determined by solving Pulay's optimisation problem equation (78). However, note that the previous Pulay predicted optimal residual is an element of the set of residuals that are used to determine {c_i}. Hence, the addition of the residual from the linear mixing step can only lower the Pulay predicted optimal residual, or at worst leave it the same. This new reduced Pulay predicted optimal residual replaces the exact linear mixing residual in the history, and the process repeats. As expected, this method performs best when the predicted optimal residual accurately models what the residual would have been were the optimal density evaluated with the Kohn–Sham map. Pulay's method predicts the residual increasingly well the closer it is to the linear response regime from the root. Therefore, when the behaviour of Kohn–Sham map is highly non-linear, the guaranteed reductions in the predicted residual tend to stagnate, while the exact residual does not decrease. All three of these techniques are benchmarked in section 6.

4.1.4. Modern multisecant-based algorithms.

Here, we highlight one modern use of multisecant methods in particular, the methods outlined in [89, 92], which are now default self-consistency methods in wien2k [114]. These algorithms can be considered a sophisticated modern variant of the standard multisecant already methods discussed. Furthermore, they are designed with the aim of converging Kohn–Sham calculations. The range of strategies utilised make [89, 92] an interesting case to isolate and examine here. The most recent published form of these algorithms is that given in [89] titled 'multisecant rank one' (MSR1). However, this algorithm is designed to converge both the atomic (geometry optimisation) and electronic degrees of freedom. Hence, we focus on the techniques that are relevant to the self-consistent field iterations, and the reader is referred to [89, 92] for a more in-depth treatment.

First, the updates considered in [89, 92] are defined in equations (90) and (89), which are of the form

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle H_n = \alpha I + (\Delta \mathcal{P}_n - \alpha \Delta \mathcal{R}_n)(\Delta \mathcal{R}_n^{\dagger} W_n)^{-1} W^{\dagger}. \label{MSR1_update} \nonumber \end{align} \tag{ 91 }$$

The initial guess inverse Jacobian is $H_0 = \alpha I$, and $W_n=\Delta \mathcal{P}_n$, $W_n=\Delta \mathcal{R}_n$ define a type-I and type-II update respectively. It is demonstrated in [92], and further verified in section 6, that type-II methods are superior for the self-consistency problem than type-I methods. However, if atomic degrees of freedom are included, it is advantageous to consider a linear combination of updates,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle W_n = Y_n + \beta S_n , \nonumber \end{align} \tag{ 92 }$$

for $\beta \in \mathbb{R}_{\geqslant 0}$. The parameter $\beta$ controls the degree to which the method takes a type-II step, $\beta = 0$, or a type-I step, $\beta \rightarrow \infty$. As noted in [89], this is similar to a technique used in [115] whereby a criterion is defined to assess whether a type-I or type-II step will be optimal, and then the corresponding step is taken. In MSR1, the parameter $\beta$ is determined based on an ansatz that seeks to ensure the eigenvalues of the Jacobian are positive, as they should be in the case that the fixed-point corresponds to a variational minimum. However, as stated, in the context of density mixing type-II methods consistently outperform type-I methods, meaning we hereafter consider $\beta=0$.

Second, a core feature of the methods in [89, 92] involve partitioning of the full update into a predicted and unpredicted component, which are now defined. Consider the update generated from equation (91),

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{{\rm in}}_{n+1} = \rho^{{\rm in}}_n + H_n R_n, \nonumber \end{align} \tag{ 93 }$$

which is now split in two,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{{\rm in}}_{n+1} = \rho^{{\rm in}}_n + (H^{\rm p}_n + H^{\rm u}_n) R_n. \nonumber \end{align} \tag{ 94 }$$

Since we are now considering the type-II variant, $W_n=\Delta \mathcal{R}_n$, the unpredicted component of the update, $H^{\rm u}_n R_n$, is defined as the orthogonal projection of the current residual R_n onto the past residual differences,

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle (H^{\rm u}_n R_n)^{\dagger} \Delta \mathcal{R}_{n} = 0. \label{ortho_proj} \nonumber \end{align} \tag{ 95 }$$

In other words, the unpredicted vector is the part of the full update that is not described within the iterative subspace of residuals. In this sense, the remaining update can be considered to be the part of the update that the iterative subspace does describe. Equation (95) is shown in [92] to have the solution,

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle H_n^u R_n = - (I - \Delta \mathcal{R}_n (\Delta \mathcal{R}_n^{\dagger}\Delta \mathcal{R}_n)^{-1} \Delta \mathcal{R}^{\dagger}_n)R_n. \nonumber \end{align} \tag{ 96 }$$

The update now takes the rearranged form

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle H_n^u R_n = (I - \Delta \mathcal{R}_n (\Delta \mathcal{R}_n^{\dagger} \Delta \mathcal{R}_n)^{-1} \Delta \mathcal{R}_n^{\dagger})R_n, \nonumber \end{align} \tag{ 97 }$$

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle H_n^{\,p} R_n = - S_n (\Delta \mathcal{R}_n^{\dagger} \Delta \mathcal{R}_n)^{-1} W_n^{\dagger} R_n. \nonumber \end{align} \tag{ 98 }$$

This partitioning is used to introduce the concept of algorithmic greed, which is quantified with an iterative damping parameter $\alpha_n$ that multiplies the unpredicted update, rather than the full update,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{{\rm in}}_{n+1} = \rho^{{\rm in}}_n + (H^{\rm p}_n + \alpha_n H^{\rm u}_n) R_n. \nonumber \end{align} \tag{ 99 }$$

In the implementation tested in section 6, the unpredicted direction is also Kerker preconditioned (see section 4.3). References [89, 92] refer to the updating of the parameter $\alpha_n$ as an implicit trust region. That is, $\alpha_n$ is allowed to increase, within some bounds, provided the algorithm is performing well, by some definition of 'well' which defines the greed controls. If the algorithm is performing poorly, the damping parameter is decreased accordingly, thus reducing the need for user intervention. The precise greed controls for MSR1 are found in [89]; the controls used for section 6 are slightly modified for performance in the plane-wave pseudopotential setting. Furthermore, the matrix inverse in equation (91) is Tikhonov regularised [116] to prevent spurious behaviour due to rank-deficiencies, and the matrices involved are scaled appropriately. The method that has been described thus far is similar to 'multisecant Broyden 2' (MSB2) of [92]. The type-I variant, MSB1, can be derived in a similar fashion, and both are tested in section 6.

Both of these methods remain prone to charge sloshing due to the primary form of step length control being an implicit trust region, the greed controls. Hence, particularly ill-conditioned simulations can still lead to divergence of the self-consistency iterations. In this context, a further stabilising measure is taken in the form of an explicit trust region. Given R_n as a descent direction on iteration n, the trust region subproblem can take a standard form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm minimise}~ || R_n - H_n^{-1} X ||_2^2 \nonumber \end{align} \tag{ 100 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm subject~to}~ ||X||^2 - \delta^2 = 0 \nonumber \end{align} \tag{ 101 }$$

for variations in X with some trust region radius $\delta$. The variable X is the new trial step generated from the trust-region subproblem, and the scalar ||X||₂ is the trail total step length, where the exact step length is $||H_n R_n ||_2$. If the exact step length exceeds the trust region radius $\delta$, then the trust region problem is required to be solved to generate this new step. Details on how the trust-region subproblem is solved are given in [83, 89]. As [89] also states, the trust region problem does not need to be defined in terms of the full update, $||H_n R_n ||_2$. Instead, the trust region problem can be solved for the predicted component of the update, $||H^{\,p}_n R_n ||_2$. This is motivated by the fact that the size of the step in the unpredicted direction is already restricted by the greed controls. Note that $\delta$ is also iteratively updated based on algorithmic progress. The inclusion of an explicit trust region for the predicted update, in conjunction with the strategies described earlier, are able to significantly stabalise the self-consistency iterations [89, 92].

A second method we highlight is a global variant of the type-I Pulay update [56], which has not yet been tested in the context of self-consistency iterations. (However, note that this method utilises a selection of the techniques from section 4.1.3 that do perform well for self-consistency iterations). First, similarly to the previous method, the Jacobian is required to be non-singular, which is achieved here through a form of Powell regularisation [117]. The regularisation parameter that is introduced, in some sense, scales the update between an unregularised Pulay step, and a fixed-point step, and thus must be chosen appropriately as not to negatively impact the efficiency of the method. The remaining modifications are aimed at stabilising the iterations and preventing stagnation in the reduction of the residual norm. Type-I methods, more so than type-II methods, tend to suffer from stagnation due to rank-deficiency in the iterative subspace of density differences $\Delta \mathcal{P}_n$, see [56]. Rank-deficiency is avoided using a technique derived from the Restarted Pulay method [112] in the previous section, i.e. occasionally restart the history of stored iterates. Here, the restart condition can be triggered based on two separate criteria, the first of which is whenever the number of iterates stored reaches a maximum value m, similarly to [112]. The second of the restart criteria is based on ensuring the iterative subspace remains approximately of full rank. This is done by first using Gram–Schmit orthonormalisation on set of vectors used to construct $\Delta \mathcal{P}_n$ in order to generate a new set of vectors with the same span, $\Delta \hat{\mathcal{P}}_n$. Elements of this new set of vectors are identically zero if the original set had linear dependencies. The condition

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle ||(\Delta \hat{\mathcal{P}}_n)_i||_2 < \tau ||(\Delta \mathcal{P}_n)_i||_2 \nonumber \end{align} \tag{ 102 }$$

is therefore used to quantify the degree of linear independence we require from the iterative subspace, parametised by $\tau$. If this condition is triggered, for some sensible value of $\tau$, the iterative history is deemed too linearly dependent, and the history is reset.

Given that update Jacobian is now bounded, and the iterative history is linearly independent, global convergence is further guaranteed by using techniques similar to the Periodic Pulay method [110, 118, 119]. That is, linear mixing steps are included in between Pulay steps, and added to the history. However, instead of performing the linear mixing steps with a fixed period, the linear mixing steps are performed based on a criterion that ensures progress is made. Namely, the following inequality is defined

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle ||R_n||_2 \leqslant D ||R_0||_2 (n' + 1)^{-(1 + \varepsilon)} \nonumber \end{align} \tag{ 103 }$$

for parameters D and $\varepsilon$, where $n'$ is the number of Pulay steps performed so far in the simulation. If the inequality is not satisfied, linear mixing steps are performed until the inequality is satisfied again, and then Pulay steps are resumed. In other words, if the Pulay steps are not making sufficient progress, then linear steps are performed until a certain amount of progress, defined by D and $\varepsilon$, has been made, and then Pulay steps are continued. The proof of global convergence for this algorithm is given in [56].

The methods introduced previously define density mixing schemes, by which we mean the one-particle density is iteratively updated in order to find a fixed-point of the Kohn–Sham equations. Density mixing is common in implementations that utilise some form of delocalised basis set. In these implementations, both the Kohn–Sham Hamiltonian and the density matrix can be prohibitive to compute and store due to the size of the basis N_b. However, in the context of localised basis sets, it is common to formulate Kohn–Sham theory so that the density matrix equation (13) is the optimised variable, rather than the density¹². The Kohn–Sham energy functional has a closed-form expression in terms of the density matrix,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E_{\tiny{\rm KS}}[D] = &-\frac{1}{2} \int_{\mathbb{R}^3} \nabla^2_{x'} D(x,x')|_{x=x'} \nonumber \\ &+ \int_{\mathbb{R}^3 \times \mathbb{R}^3} \frac{D(x,x)D(x',x')}{|x-x'|} \nonumber \\ &+ \int_{\mathbb{R}^3} v_{{\rm ext}}D(x,x) + E_{{\rm xc}}[D(x,x)] \label{ks_objective_dmatrix} \nonumber \end{align} \tag{ 104 }$$

which is now minimised over allowed variations in D,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int_{\mathbb{R}^3} D(x,x) = N, \nonumber \end{align} \tag{ 105 }$$

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle D^{\dagger} = D, \nonumber \end{align} \tag{ 106 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int_{\mathbb{R}^3} D(x,x'') D(x'',x') = D(x,x'), \label{binary_occ_D} \nonumber \end{align} \tag{ 107 }$$

for binary occupancy of the Kohn–Sham orbitals. These conditions define the set of density matrices $\mathcal{D}$ which form the domain of $E_{\tiny{\rm KS}}[D]$, i.e. the set of density matrices for which the Kohn–Sham functional is defined. These are necessary but not sufficient conditions for a density matrix to be so-called N-representable, meaning it can be constructed from an antisymmetric wavefunction [120]. The first condition forces the density matrix to correspond to an N-particle system. The second condition forces the density matrix to be Hermition, which is true by construction. The third condition is called idempotency, which leads to the requirement that the eigenvalues of discretised density matrix, the orbital occupancies, are f _i  =  {0,1} [28]. The Kohn–Sham Hamiltonan, sometimes referred to as the Fock matrix, can also be constructed and used to solve the Kohn–Sham equations for a self-consistent density matrix. As expected, a variety of self-consistent field techniques exist that are well suited to the this formulation, such as leveling shifting and its modern variants [121–124], methods that minimise a local model energy functional [122, 125–127], and more [19, 60, 128–131]. As many of these methods are not readily adaptable to the plane-wave setting, we keep the discussion here brief, and instead refer readers to the following review articles [19, 130, 132, 133] and references therein. In fact, we highlight one class of methods in particular, the relaxed constraints algorithms given in [49, 60, 129, 131, 134]. A member of this class, 'Energy DIIS', is now the default self-consistency method in gaussian09 [18].

Relaxed constraints algorithms are introduced in the context of Hartree–Fock theory, where, unlike in extended Kohn–Sham theory discussed in section 3.6, the notion of fractional occupancy is not a part of the framework. Relaxed constraints algorithms operate by permitting fractional occupation of the Hartree–Fock orbitals as a tool to reach convergence. The binary occupation fixed-point solution is recovered at the end of the calculation. First, consider the set of allowed discretised density matrices,

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \mathcal{D} = \left\{D \ | \ {{\rm Tr}}D = N, \ D^{\dagger}=D, \ D^2 = D \right\}. \nonumber \end{align} \tag{ 108 }$$

It is now apparent that the idempotency condition equation (107) requires that the eigenvalues of the density matrix, the orbital occupancies, are binary, $f_i \in \{0,1\}$. The extension to fractional occupancy, as described in section 3.6, thus alters the set of allowed density matrices,

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \tilde{\mathcal{D}} = \left\{D \ | \ {{\rm Tr}}D = N, \ D^{\dagger}=D, \ D^2 \leqslant D \right\}; \nonumber \end{align} \tag{ 109 }$$

i.e. the eigenvalues of the density matrix must now satisfy $f_i^2 \leqslant f_i$, which is the case for $0 \leqslant f_i \leqslant 1$. Relaxed constraint algorithms are now founded based on two theorems, the proofs of which are given in [135]. The first theorem states that the Hartree–Fock functional varied over $\mathcal{D}$ has the same stationary points as the Hartree–Fock functional varied over $\tilde{\mathcal{D}}$. In other words, minimising the Hartree–Fock functional over $\tilde{\mathcal{D}}$ will always lead to physical (binary occupation) solutions. The second theorem states that the set $\tilde{\mathcal{D}}$ is convex. This means that, given a convex combination of density matrices $D_1, D_2 \in \tilde{\mathcal{D}}$,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle D = (1 - \gamma) D_1 + \gamma D_2, \nonumber \end{align} \tag{ 110 }$$

for $\gamma \in [0,1]$, then $D \in \tilde{\mathcal{D}}$. This is not the case for elements of $\mathcal{D}$. Algorithms that utilise these theorems in conjunction with the set $\tilde{\mathcal{D}}$ belong to the class of relaxed constraints algorithms.

The question remains of how these theorems translate to the Kohn–Sham functional equation (104). In fact, the former theorem no longer holds [134], meaning that the functional varied over $\tilde{\mathcal{D}}$ does lead to different stationary points than in the case of binary occupation. This is expected, as these solutions correspond to solutions of the extended Kohn–Sham theory developed in section 3.6. Therefore, the absence of the former theorem does not pose a problem, as solutions that are obtained that are members of $\tilde{\mathcal{D}}$ retain meaning within extended Kohn–Sham theory. Furthermore, the set $\tilde{\mathcal{D}}$ remains convex. The optimal damping algorithm was the first of the relaxed constraints algorithms, as examined in [60]. This algorithm seeks, at each iteration, to find an 'optimal' damping parameter $\alpha_n$ in the linear mixing scheme equation (22). Consider the pair $\{D^{{\rm in}}_n, D^{{\rm out}}_n \}$ on iteration n. A trail density matrix is now constructed as a convex combination of this pair,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle D^{{\rm trial}} = (1 - \alpha)D^{{\rm in}}_n + \alpha D^{{\rm out}}_n, \label{line_seg} \nonumber \end{align} \tag{ 111 }$$

for $\alpha \in [0,1]$, where we note that $D^{{\rm out}}_n - D^{{\rm in}}_n$ is a descent direction. The Kohn–Sham energy functional $E_{{\rm ks}}[D^{{\rm trial}}]$ is then minimised along this line-segment by varying $\alpha$. The value of the damping parameter that leads to the minimum energy is the so-called optimal damping value $\alpha^{{\rm opt}}$. The subsequent input density matrix is thus

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle D^{{\rm in}}_{n+1} = D^{{\rm trial}} = D^{{\rm in}}_n + \alpha^{{\rm opt}} (D^{{\rm out}}_n - D^{{\rm in}}_n). \nonumber \end{align} \tag{ 112 }$$

Note that this method works precisely because the elements along the line-segment equation (111) remain in the domain of $E_{\tiny{\rm KS}}$. Furthermore, it is much cheaper to evaluate the energy for a given density matrix, rather than construct and diagonalise the Fock matrix for a new output density matrix.

This method is improved by instead considering an (at most) m-dimensional iterative subspace of past density matrices, rather than just the most recent pair. In other words, similar to Pulay's method, a trail density matrix is constructed

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle D^{{\rm trial}} = \sum_{i=n-m+1}^m c_i D^{{\rm in}}_i \nonumber \end{align} \tag{ 113 }$$

for some unknown coefficients c_i. These coefficients are then fixed as the coefficients that minimise $E_{\tiny{\rm KS}}[D^{{\rm trial}}]$. An algorithm to accomplish this is provided in [131]. However, note that only convex combinations of past density matrices are permissible, otherwise the resulting density matrix may not be a valid according to the conditions imposed in $\tilde{\mathcal{D}}$. This observation necessitates the restriction $0 \leqslant c_i \leqslant 1$, which leads to the following constrained minimisation problem,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm minimise}~E_{\tiny{\rm KS}} \left[\sum_{i=n-m+1}^m c_i D^{{\rm in}}_i \right] \label{EDIIS} \nonumber \end{align} \tag{ 114 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm subject~to}~ 0 \leqslant c_i \leqslant 1. \nonumber \end{align} \tag{ 115 }$$

The coefficients that solve this problem form the optimal density matrix that is set equal to the subsequent input density matrix $D^{{\rm in}}_{n+1}$. This method is titled 'Energy DIIS' (EDIIS) and has been demonstrated to perform well in a variety of cases [130, 131]. Furthermore, this method is global, as the energy functional is required to be minimised at each iteration in equation (114). The minimisation of the energy functional is an interpolation step in the space of past density matrices due to the constraint on the coefficients, and can be slow when the iterates are near convergence [131]. For this reason, EDIIS is commonly augmented with Pulay iterations (DIIS) when close to convergence, which demonstrably improves efficiency [130]. These methods, and similar methods, can be impractical in plane-wave codes as one is required to construct and store the density matrix.

Preconditioning refers to the modification of an optimisation problem such that the condition number of the problem is improved. Crucially, the modified problem is required to have the same minimum, and minimiser, as the original problem. Algorithms applied to the modified problem thus have more stable and accelerated convergence. A preconditioner is most transparently defined for linear systems as being the matrix P such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P^{-1} A x = P^{-1} b, \nonumber \end{align} \tag{ 116 }$$

where the P⁻¹A has a lower condition number equation (43) than A. In the case P  =  A the linear system is solved, and P is the exact preconditioner. The definition of a preconditioner for non-linear systems is less transparent. Consider the optimisation of the Kohn–Sham residual L²-norm, which now takes the preconditioned form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle ||P(R[\rho_*])||_2 = 0. \label{precond_residual} \nonumber \end{align} \tag{ 117 }$$

The preconditioned residual $P(R)$ is required to have the same solution as R, but has, in some sense, improved convergence properties. The perfect preconditioner here would modify R such that only one step of appropriate size in the steepest descent direction is required for convergence. A successfully preconditioned problem therefore represents a problem whose landscape is easier to traverse toward a minimum using, for example, Newton's algorithm. It is known that the Jacobian (Hessian) eigenvalue spectrum of the objective function determines the rate of convergence of Newton, and quasi-Newton, methods [58, 83, 97]. Therefore, a preconditioner should accomplish one or more of the following: reduce the number of eigenvalue clusters; reduce the width of the eignevalue clusters; or compress the spectrum as a whole. A discussion on the importance of the clustering of the eigenvalues, rather than just the condition number, can be found in [83, 89, 136]. To simplify matters, we consider P to be a matrix constant with respect to the optimised variable, the density, such that $P(R) = P R$. Hence, in the present context, preconditioning amounts to finding the matrix P such that the spectrum of the dielectric is more suitable for quasi-Newton algorithms. Note that the preconditioning matrix is permitted to change iteration-to-iteration. The strategy used in practice is to identify the source of divergent eigenvalues of the dielectric, as examined in section 3, and temper this divergence in such a fashion that is generally applicable to all, or large classes, of Kohn–Sham inputs. Depending on the implementation, the preconditioning approach can differ. For example, in augmented plane-wave implementations [114], the unit cell is partitioned into the regions surrounding atomic cores, represented by local basis functions, and an interstitial region, represented by plane-waves. Naturally, due to the differing number of basis functions involved in each region, among other properties, the preconditioning for each region is separate [89, 92]. The following work assumes an entirely plane-wave basis set, although the preconditioners can be adapted for a variety of implementations.

Recall from section 3.8 that the Coulomb kernel, in combination with the susceptibility, is principally responsible for the large eigenvalues of the Kohn–Sham residual linear response function. As discussed, in the general case, the susceptibility is a complicated object about which it is difficult to make sweeping statements. However, when the system is homogeneous and isotropic, the dielectric eigenvalues of gapped and gapless phases are approximately determined by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \epsilon_0^{-1} = \left(1 + \frac{4\pi \gamma}{|G|^2} \right)^{-1}, \label{response_metal} \nonumber \end{align} \tag{ 118 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \epsilon_0^{-1} = \left(1 + 4\pi \gamma \right)^{-1} \label{response_insulator}, \nonumber \end{align} \tag{ 119 }$$

for some a priori unknown system-dependent constant $\gamma$. In the case of the homogeneous electron gas, this constant in equation (118) is identified with the square of the Thomas–Fermi screening wavevector $k_{\tiny{\rm TF}}^2$. Modest departures from homogeneity and isotropy remain accurately modelled by equations (118) and (119), particularly in the low $|G|$ limit [75, 88]. Therefore, these model dielectrics can be used to improve the condition of the residual map by allowing $\newcommand{\e}{{\rm e}} P=\epsilon_0^{-1}$ for either equation (118) or (119) depending on whether one suspects the input to be metallic or insulating. Collecting and relabelling the unknown constants, the preconditioner becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P = \alpha \frac{|G|^2}{|G|^2 + |G_0|^2}, \label{kerker} \nonumber \end{align} \tag{ 120 }$$

where |G₀| and $\alpha$ are parameters that are determined by the linear response of the input system; e.g. |G₀|  =  0 for Kohn–Sham insulators, and |G₀| is related to the Wigner–Seitz radius for the homogeneous electron gas. The values of |G₀| and $\alpha$ naturally differ depending on the input, although fixing |G₀|  =  1.5 $\mathring{\rm A}$ and $\alpha=0.8$ [15, 31] demonstrably improves convergence, see section 6. The modified Kohn–Sham problem is now solved

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Bigg| \Bigg| \alpha \frac{|G|^2}{|G|^2 + |G_0|^2} R[\rho_*(G)] \Bigg| \Bigg|_2 = 0, \nonumber \end{align} \tag{ 121 }$$

which is referred to as Kerker preconditioned equation (120) [84, 85]. The Kerker preconditioner suppresses charge sloshing, as defined in section 3.8, by damping eigenvalues of the dielectric corresponding to low $|G|$ components of the density, see figure 8. It is these components that have a generically amplified response due to the Coulomb kernel. Note that the optimal damping algorithm detailed in section 4.2 can be regarded as a preconditioner that updates the value of $\alpha_n$ at each iteration. The adaptation of the Kerker preconditioner to real space implementations of Kohn–Sham theory more difficult. The dielectric response function of the homogeneous electron gas is non-local in real space, meaning the integral in equation (39) leads to a dense $N_b \times N_b$ matrix that must be computed and stored. The susceptibility takes the Yukawa screening form in real space [76]. An efficient real space implementation of the Kerker preconditioner is given in [137].

Zoom In Zoom Out Reset image size

One aspect of the preconditioning presented here is that, for increasingly homogeneous and isotropic inputs, the dependence of the dielectric condition number on unit cell dimension $L \sim G^{-1}$ cancels [73]. This is by construction, and is identified by examining the eigenvalues of the Kerker preconditioned dielectric. For inputs that are not so accommodating, while the Kerker preconditioner does help, the scaling of iterations with unit cell dimension persists [73, 138]. Removing this scaling can be considered one of the primary goals of preconditioners in Kohn–Sham theory as large simulation cells are required for many modern applications. Moderate extensions to the Kerker preconditioner have been proposed [139], which involve more accurately modelling the dielectric response, e.g. [140]. However, in these examples, the main drawback of the Kerker preconditioner remains: there is no scope for systematically including anisotropy and inhomogeneity. Furthermore, the exchange-correlation kernel is ignored. This is a reasonable practice in non spin polarised systems, but in spin polarised systems, the spin density interacts entirely through the exchange-correlation kernel, and is thus not preconditioned.

There have been various efforts to construct a preconditioner that provides an improved description of inhomogeneous and anisotropic inputs. In the most extreme examples, computation of the exact linear response function is considered [78–81, 86, 141]. As discussed, the exact linear response function does not represent a preconditioning scheme, rather it is the exact Newton method. A density-dependent preconditioner necessarily alters the Jacobian, equation (117), in a non-trivial manner. Hence, it is not obvious the extent to which the condition of the modified problem will improve. Nevertheless, even approximate attempts at computing the exact response function, when treated as a Newton step, are able to improve the iterations over accelerated fixed-point algorithms such as Pulay's method [78, 141]. The central difficulty in computing the exact susceptability, and subsequently the exact dielectric, is that it requires a summation over all unoccupied-occupied eigenfunction pairs equation (49). In plane-wave codes, both the computation and storage of all eigenfunctions of the Kohn–Sham Hamiltonian is infeasible.

Past attempts [79, 80] are able to implement and examine computation of the exact susceptibility with success, albeit with a basis set size now unsuited to modern computation. The problem of having to compute the full set of eigenfunctions, and having too large a basis, can be remedied with a few differing approaches. First, the size of the effective basis can be reduced by recalling that the low $|G|$ components of the density are those responsible for divergent eigenvalues of the dielectirc. Therefore, the susceptibility need only be computed for a reduced set of plane-waves—those with low frequency G. Working in this reduced space for the purposes of density mixing leads to a significant reduction in compute and memory overhead¹³ [15, 31, 78]. Second, we highlight two strategies to remove the infinite summation over unoccupied eigenfunctions in the Adler–Wiser equation equation (49). The first strategy utilises density functional pertubation theory, and in particular the Sternheimer equation [142]. Solution of the Sternheimer equation allows one to obtain the first-order response in a perturbed quantity—here, the density in the direction of the residual—without requiring the full eigendecomposition of the Kohn–Sham Hamiltonian. This is utilised in, for example, the implementation of the GW approximation [143], and hence is available functionality within many Kohn–Sham codes. The density functional pertubation theory approach is used to construct a Newton step in [141] and is shown to significantly reduce the number of iterations taken to converge when compared to Pulay's method. The second strategy involves exploiting the completeness of eigenfunctions,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \delta(x - x') = \sum_{j=1}^{\infty} \phi^*_i(x) \phi_j(x'). \label{completeness} \nonumber \end{align} \tag{ 122 }$$

This identity cannot be applied to the Adler–Wiser equation equation (49) without further approximation due to the denominator depending on the j th eigenenergy, $\newcommand{\e}{{\rm e}} \epsilon_j$. However, these eigenenergies can be replaced with some approximate constant value $\bar{E}$ above a certain cut-off number of eigenfunctions $N_{{\rm cut\mbox{-}off}}$. This allows equation (122) to transform the Adler–Wiser equation into a sum over $N_{\rm cut\mbox{-}off} << N_b$ eigenfunctions, which is used [78, 81] to successfully reduce the number of iterations taken to converge. However, note that these methods retain a poor scaling with the number of electrons; $\mathcal{O}(N^4)$ in the most recent examples. Whilst the prefactor of the scaling is much reduced compared to earlier efforts [79], such methods, without further approximation or development, are precluded for larger system sizes. A major advantage of the methods discussed here though is the ability to take into account the spin linear response function. To the authors knoweldge, no other methods based around model dielectrics attempt to include the spin response function, and thus the spin density is often not preconditioned.

Finally, methods are examined that attempt to include inhomogeneity and anisotropy through model response functions, rather than with the exact methods discussed prior. An extension to the Kerker preconditioner, which is based on Thomas–Fermi theory of the homogeneous electron gas, is considered [144]. Here, Thomas–Fermi–Von Weizäcker¹⁴ theory is used, in combination with Pulay's method, to generate the subsequent density. Crucially, this work does not restrict to the case of the homogeneous electron gas, and instead numerically minimises the relevant functional. That is, the modified orbital-free functional

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle E_{\tiny{\rm TFW}}[\rho] = &\int_{\mathbb{R}^3} \big| \nabla \rho(x)^{\frac{1}{2}} \big|^2 + \frac{3}{5} \int_{\mathbb{R}^3} \ \rho(x)^{\frac{5}{3}} \nonumber \\ &- \int_{\mathbb{R}^3} \ v_{{\rm ef\,\!f}}[\rho](x)\rho(x) + \delta E[\rho^{{\rm in}}_n,\rho] \label{tfwfunctional} \nonumber \end{align} \tag{ 123 }$$

is minimised with respect to variations in $\rho$, where $v_{\rm ef\,\!f} = v_{\rm h} + v_{{\rm xc}} + v_{\rm ext}$, the effective potential used to solve the Kohn–Sham equations. The details of the modification term $\delta E[\rho^{{\rm in}}_n,\rho]$ are given in [144], and is derived such that the minimiser of equation (123) can be used as the subsequent input density. Namely, the method calculates the repsonse that is required to bring the current input density to self-consistency within the framework of Thomas–Fermi–Von Weizäcker theory, and uses this as a model of the exact Kohn–Sham response. Minimisation of the Thomas–Fermi–Von Weizäcker functional, which is done using the conjugate gradient method [144], is vastly more efficient than minimisation of the Kohn–Sham functional due to that fact it is orbital-free. In certain test cases this method is demonstrated to reduce the time taken to converge by up to a factor of three, and is implemented in the software package abinit [10, 32].

Alternatively, given an input that can be transparently partitioned into metallic and insulating regions, such as an interface, inhomogeneity can be included explicitly by varying |G₀| in each region. That is, set |G₀|  =  0 for the insulating region, and have finite |G₀| in the metallic region, which is considered in [73, 145]. This is non-trivial as the dielectric equation (39) becomes non-diagonal in both Fourier and real space, and hence becomes unfavourable to construct and store. However, note that one requires the dielectric applied to the residual vector R_n, rather than the dielectric itself. Therefore, inhomogeneity can be included using an algorithm that successively switches between Fourier and real space as to avoid constructing the non-diagonal dielectric; see [145]. Furthermore, if the potential is treated as the optimisation variable rather than the density, a modified Poission equation can be solved for the updated potential residual, see [73]. Here, the inhomogeneity is specified a priori with two functions that are now inputs to the calculation. These methods thus have the drawback that they are not black-box, as one is required to specify the inhomogeneity using prior knowledge of ones input system. Nonetheless, for specific systems, these frameworks provide an expert user with an additional degree of freedom for optimising convergence.

Whilst self-consistent field methods are widespread, alternative techniques are available based on direct minimisation of the energy functional. These methods exploit the variational principle, and are thus global, varying $\{\phi_i \}$ to minimise $E_{\tiny{\rm KS}}[\{\phi_i \} ]$. The density in these schemes is a dependent quantity, always derived directly from $\{\phi_i \}$ with no history of previous densities. From an initial guess set of orbitals, $\{\phi^{(0)}_i \}$, a density

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{(0)}(x)=\sum_i \left\vert\phi^{(0)}_i(x)\right\vert^2 \nonumber \end{align} \tag{ 124 }$$

is constructed, and the energy $E_{\tiny{\rm KS}}[\{\phi^{(0)}_i \} ]$ and Kohn–Sham eigenvalue estimates $\newcommand{\e}{{\rm e}} \epsilon^{(0)}_i$ are evaluated, along with the energy gradient with respect to changes in $\{\phi^{(0)}_i \}$,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left.\frac{\delta E_{\tiny{\rm KS}}}{\delta \phi_i }\right\vert_{\phi_i=\phi^{(0)}_i} = H_{\tiny{\rm KS}}[\rho]\phi^{(0)}_i - \epsilon^{(0)}_i\phi^{(0)}_i . \nonumber \end{align} \tag{ 125 }$$

Since the energy gradients are the steepest ascent directions, the steepest descent direction is the negative of this. This steepest descent direction may be interpreted as a force acting on the degrees of freedom of the trial states $\{\phi^{(0)}_i \}$. If masses are assigned to these degrees of freedom, then the states may be evolved forward in time according to a suitable equation of motion, and this forms the foundation of the Car–Parrinello method [17, 146]. By damping the motion appropriately, the system evolves towards the ground state.

An alternative approach is to consider the search for the ground state as a minimisation problem. A candidate search direction $D^{(0)}_i$ may be constructed to minimise the energy as,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle D^{(0)}_i = -\left.\frac{\delta E_{\tiny{\rm KS}}}{\delta \phi_i }\right\vert_{\phi_i=\phi^{(0)}_i} \nonumber \end{align} \tag{ 126 }$$

i.e. the steepest descent direction. In practice more sophisticated methods are used to construct a search direction, usually based on preconditioned quasi-Newton methods such as conjugate gradients.

Once a search direction has been obtained, an improved set of trial orbitals are constructed, e.g.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \phi^{(1)}_i = \phi^{(0)}_i + \alpha D^{(0)}_i , \nonumber \end{align} \tag{ 127 }$$

where $\alpha$ is a scalar parameter, chosen to minimise $E_{\tiny{\rm KS}}[\{\phi^{(1)}_i \} ]$. Note that, in general, $\{\phi^{(1)}_i\}$ will not be orthonormal and must be orthonormalised explicitly first. The search for an optimal value of $\alpha$ is known as the line-minimisation step. It is also important to note that, in the evaluation of $E_{\tiny{\rm KS}}[\{\phi^{(1)}_i \} ]$, $v_{\rm h}$ and $v_{\rm xc}$ are always constructed from the density

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^{(1)}(x)=\sum_i \left\vert\phi^{(1)}_i(x)\right\vert^2. \nonumber \end{align} \tag{ 128 }$$

This is the critical difference between the self-consistent field methods and the direct energy minimisation methods. In self-consistent field methods the corresponding optimisation of the orbitals is carried out using the original $v_{\rm h}$ and $v_{\rm xc}$. Figure 9 shows a direct comparison between these two approaches for a simulation of silicon, using an 8-atom conventional unit cell. The effect of updating the Kohn–Sham potential at each step along the line minimisation is to increase the curvature of the energy with respect to the step-length $\alpha$ leading to higher energies and an energy minimum at a smaller value of $\alpha$. In contrast, the over-estimation of $\alpha$ when the energy curve along the steepest descent direction is computed with the fixed initial potential (from $\alpha=0$) can be considered one of the root causes of charge-sloshing instabilities in self-consistent field methods, see section 3.8.

Zoom In Zoom Out Reset image size

The direct energy minimisation method discussed thus far is suitable for simulations with a band-gap. For metals and finite temperature insulators, however, it is not sufficient to consider only the lowest N eigenstates as occupied and there is an additional dependence on the partial occupancies f _i. The ground state must now be found by searching over both the Kohn–Sham states and their occupancies. One of the most robust methods of this form is the ensemble density functional theory (EDFT) method of Marzari et al [26]. In EDFT every update of the trial states $\{\phi_i \}$ is followed by a direct energy-minimisation over the density matrix in the basis of the trial states occupancies, {f _ij}. This density matrix is nothing more than a generalisation of the occupancies to the case when the trial states $\{\phi_i \}$ do not diagonalise the subspace Hamiltonian h directly, where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_{ij} = \int {\rm d}x \ \phi^\ast_i(x)H_{\tiny{\rm KS}}[\rho] \phi_j(x). \nonumber \end{align} \tag{ 129 }$$