Guaranteed Convergence of a Regularized Kohn-Sham Iteration in Finite Dimensions

TL;DR

This paper proves that a regularized Kohn-Sham iteration with adaptive damping reliably converges to the true ground-state density in finite-dimensional density-functional theory.

Contribution

It demonstrates convergence of a regularized Kohn-Sham iteration with adaptive damping in finite dimensions, ensuring correct ground-state density calculation.

Findings

Convergence is guaranteed with Moreau-Yosida regularization.

Adaptive damping ensures stability and convergence.

The method accurately finds the ground-state density.

Abstract

The exact Kohn-Sham iteration of generalized density-functional theory in finite dimensions witha Moreau-Yosida regularized universal Lieb functional and an adaptive damping step is shown toconverge to the correct ground-state density.

Equations

$$\tilde{F}(\rho )=\underset{\Gamma \mapsto \rho }{\inf }\left\{\mathrm{Tr}((H_{\mathrm{kin}}+H_{\mathrm{int}})\Gamma )\right\}.$$

$$E(v)=\underset{\rho \in \tilde{X}}{\inf }\{\tilde{F}(\rho )+⟨v,\rho ⟩\}.$$

$$F(x)=\underset{v\in X^{\ast }}{\sup }\{E(v)-⟨v,x⟩\},$$

$$E(v)=\underset{x\in X}{\inf }\{F(x)+⟨v,x⟩\}.$$

$$F_{\epsilon }(x)=\underset{y\in X}{\inf }\left\{F(y)+\frac{1}{2\epsilon }\Vert x-y{\Vert }^2\right\}.$$

$$E_{\epsilon }(v)=\underset{x\in X}{\inf }\{F_{\epsilon }(x)+⟨v,x⟩\}.$$

$$\nabla F_{\epsilon }(z)+v=0.$$

$$E(v)$$

$$\overline{\partial }E(v)$$

$${\tilde{F}}^0(\rho )=\underset{\Gamma \mapsto \rho }{\inf }\left\{\mathrm{Tr}(H_{\mathrm{kin}}\Gamma )\right\}$$

$$\nabla F_{\epsilon }(z)+v=\nabla F_{\epsilon }^0(z)+v_{\mathrm{KS}}^{\epsilon },$$

$$v_{i+1}=v+\nabla F_{\epsilon }(x_i)-\nabla F_{\epsilon }^0(x_i),$$

$$⟨\nabla F_{\epsilon }(x_i)+v,x_{i+1}^{\prime }-x_i⟩=⟨v_{i+1}+\nabla F_{\epsilon }^0(x_i),x_{i+1}^{\prime }-x_i⟩.$$

$${\tilde{x}}_{i+1}^{\prime }=x_{i+1}^{\prime }+\epsilon v_{i+1}\in \overline{\partial }{E}^0(v_{i+1}),$$

$${\tilde{x}}_i=x_i+\epsilon \nabla F_{\epsilon }^0(x_i)\in \overline{\partial }{E}^0(\nabla F_{\epsilon }^0(x_i)),$$

$$⟨v_{i+1}+\nabla F_{\epsilon }^0(x_i),{\tilde{x}}_{i+1}^{\prime }-{\tilde{x}}_i⟩-\epsilon \Vert v_{i+1}+\nabla F_{\epsilon }^0(x_i){\Vert }^2.$$

$$⟨\nabla F_{\epsilon }(x_i)+v,x_{i+1}^{\prime }-x_i⟩$$

$$\le -\epsilon \Vert v_{i+1}+\nabla F_{\epsilon }^0(x_i){\Vert }^2$$

$$=-\epsilon \Vert \nabla F_{\epsilon }(x_i)+v{\Vert }^2.$$

$$⟨\nabla F_{\epsilon }(x_i)+v,y_i⟩$$

$$=-\frac{1}{\epsilon }\Vert x_{i+1}-x_i\Vert =-\frac{{\tau }_i}{\epsilon }.$$

$$\frac{{\tau }_i^2}{2\epsilon }=\frac{1}{2\epsilon }\Vert x_{i+1}-x_i{\Vert }^2=e_i-m_i\le e_i-e_{i+1}\to 0,$$

$$\underset{i\to \infty }{\lim }{v}_{i+1}=v+\nabla F_{\epsilon }^1(x)-\nabla F_{\epsilon }^0(x),$$

$$\Vert x_{i+1}^{\prime }-x_i\Vert ⟨\nabla F_{\epsilon }(x_i)+v,y_i⟩\le -\epsilon \Vert \nabla F_{\epsilon }(x_i)+v{\Vert }^2,$$

$$\Vert x_{i+1}^{\prime }-x_i\Vert \frac{{\tau }_i}{\epsilon }\ge \epsilon \Vert \nabla F_{\epsilon }(x_i)+v{\Vert }^2.$$

$$x_{i+1}^{\prime }\in \overline{\partial }{E}_{\epsilon }^0(v_{i+1})=\overline{\partial }{E}^0(v_{i+1})-\epsilon v_{i+1},$$

$$\rho =z+\epsilon v,{\rho }_{\mathrm{KS}}^{\epsilon }=z+\epsilon v_{\mathrm{KS}}^{\epsilon },$$

Full text

The Kohn–Sham (KS) scheme Kohn and Sham (1965) of ground-state density-functional theory (DFT) is the cornerstone of electronic structure calculations in quantum chemistry and solid-state physics Burke (2012). It maps a complicated system of interacting electrons onto an auxiliary, non-interacting KS system. This yields a set of coupled one-particle equations that need to be solved self-consistently. Since a direct solution is unfeasible, practical approaches are variations of self-consistent field methods taking the form of fixed-point iterations or energy minimization algorithms Karlström (1979); Rabuck and Scuseria (1999); Cancès and Le Bris (2000); Helgaker et al. (2000); Kudin and Scuseria (2002); Høst et al. (2008). To date, no method has been rigorously shown to converge to the correct ground-state density. Convergence results for approximate schemes are available for auxiliary assumptions Liu et al. (2015), and reliably achieving convergence in systems with small band gaps or for transition metals remains a hard practical challenge Kudin and Scuseria (2007). Approximation techniques face the problem of an exponential growth of local minima with increasing number of particles Fukutome (1971). Such local minima appear as ‘false’ solutions in the energy landscape and distract from the global, absolute minimum Stillinger (1999). Hence, a method with mathematically guaranteed convergence to the correct minimizer is of central importance and has been listed as one of twelve outstanding problems in DFT Ruzsinszky and Perdew (2011).

An early insight is that iterations commonly fail unless oscillations between trial states are damped Karlström (1979). Work by Cancès and Le Bris Cancès and Le Bris (2000); Cancès (2001) led to the optimal damping algorithm (ODA) based on energy minimization by line search along the descent direction. Wagner et al. (2013, 2014) presented a similar scheme and claimed to have proven convergence in the setting of exact DFT, while only the strict descent of energies was secured. In such efforts, functional differentiability is almost always tacitly assumed or wrongly claimed, prominently in Ref. (Engel and Dreizler, 2011, Eq. (2.105)), while the underlying universal functionals are known to be non-differentiable Lammert (2007). This means the usual presentations of DFT already assume some form of regularization of the functionals. Other special forms of DFT like with internal magnetic fields Tellgren (2018) or finite temperatures Chayes et al. (1985); Giesbertz and Ruggenthaler (2019) automatically include regularization effects.

This issue was addressed in Laestadius et al. (2018), where a similar iterative scheme was proposed that proved a weak type of convergence after Moreau–Yosida (MY) regularization to ensure differentiability of the universal Lieb functional Lieb (1983). Weak-type convergence here means that the energy converges to either the correct energy or an upper bound. MY regularization has been introduced to DFT by Kvaal et al. (2014).

A rich study of possible strategies for self-consistent field iteration was recently put forward by Lammert . Yet in all those works the question of a limit density and corresponding KS potential was left open. On the other hand, the result in Laestadius et al. (2018) is applicable to not only standard DFT, but to all DFT flavors that fit into the given framework of reflexive Banach spaces. It has already been successfully applied to paramagnetic current DFT (CDFT) Laestadius et al. (2019). This general approach is also pursued in this Letter.

In what follows we give a fully rigorous proof of convergence for the KS scheme in a finite-dimensional state space. Because of the techniques involved, the new iteration scheme was baptized “MYKSODA” in Laestadius et al. (2019). The employed damping critically depends on MY regularization that bounds the curvature of the universal Lieb functional from above.

We will now present the mathematical framework. For a much more detailed discussion of generalized KS schemes in Banach spaces that can also be infinite dimensional, we refer to Laestadius et al. (2018). The spaces for densities and potentials are chosen to be the Hilbert space $X=X^{*}=\ell^{2}(M)$, $M\in\mathbb{N}$, which corresponds to a finite one-particle basis, a lattice system with $M$ sites, or many other possible settings. The reason for this dual choice of spaces is how densities and potentials couple in the energy expression. What is denoted $\int_{\mathbb{R}^{3}}v\rho\,\mathrm{d}x$ in standard DFT is a finite sum in the given setting and will further be written $\langle v,\rho\rangle$ with $\rho\in X,v\in X^{*}$. For the internal energy of the full system, a universal functional $\tilde{F}$, like the one defined by constrained search Levy (1979); Lieb (1983) over all $N$-particle density matrices $\Gamma$ that yield a given density $\rho\in X$, is introduced,

[TABLE]

Here $H_{\mathrm{kin}}$ stands for the kinetic energy and $H_{\mathrm{int}}$ for interactions. Consequently, the functional $\tilde{F}$ is defined on a set $\tilde{X}\subset X$ of physical densities that come from an $N$-particle density matrix (ensemble $N$-representability). This set $\tilde{X}$ will be assumed bounded in $X$. Since all physical densities are normalized in the $\ell^{1}$ norm and all norms are equivalent in finite dimensions, this follows naturally. It also holds for CDFT on a finite lattice, since the current density is bounded by the hopping parameter (Farzanehpour and Tokatly, 2012, Eq. (25)), and for one-body reduced density matrix functional theory (RDMFT) in finite basis sets, since the off-diagonal elements of the reduced density matrix are bounded by the diagonal ones that give the usual density (Giesbertz and Ruggenthaler, 2019, Eq. (3.49)).

On the other hand, elements in $X$ will in general not constitute physical densities. In standard DFT this means that an arbitrary $x\in X$ does not have to be normalized or even positive. Such an $x\in X$ will thus be called a quasi-density. We reserve the notation $\rho$ for physical densities.

The total energy is the infimum of $\tilde{F}(\rho)$ plus the the potential energy coming from a given external potential $v\in X^{*}$, taken over all physical densities,

[TABLE]

It is linked to a functional $F$ on $X$ by the Legendre–Fenchel transformation (convex conjugate). Then $F$ can be transformed back to the same $E$ as

[TABLE]

The functional $F$ is by construction convex and lower semi-continuous (Lieb, 1983, Th. 3.6) and has $F(x)=+\infty$ whenever $x$ is not in the domain $\tilde{X}$ of $\tilde{F}$. Minimizers of (2) are the ground-state densities, which establishes a link to the Schrödinger equation. They stay the same if one switches from $\tilde{F}$ to $F$, and thus minimizers of (4) are always in $\tilde{X}$. Finding such minimizers $\rho$ of (4) is equivalent to determining the superdifferential of $E$, i.e., the set of functionals in $X^{**}=X$ that yield a graph completely above $E$, written $\rho\in\overline{\partial}E(v)\subset\tilde{X}$.

The MY regularization of the functional $F$ on $X$ is defined as

[TABLE]

The visual understanding of this is the following. As the vertex of the regularization parabola $\frac{1}{2\varepsilon}\|x\|^{2}$ moves along the graph of $F$, the regularized $F_{\varepsilon}$ is given by the traced out lower envelope (the “Moreau envelope” (Rockafellar and Wets, 2009, Def. 1.22)). This is visualized in Fig. 1. It also means the regularization puts an upper bound of $\varepsilon^{-1}$ on the (positive) curvature of $F_{\varepsilon}$. This will be an important ingredient in the convergence proof: A bound on the curvature means the convex functional $F_{\varepsilon}$ cannot change from falling to rising too quickly, yielding a secure bound on the possible step length for descent.

The regularized $F_{\varepsilon}$ is then differentiable and even has a continuous gradient $\nabla F_{\varepsilon}$ (Fréchet differentiability) (Laestadius et al., 2018, Th. 9), something that will also become important in the convergence proof. We define the associated energy functional

[TABLE]

The functional in (6) is not the MY regularization of $E$ but the Legendre–Fenchel transformation of $F_{\varepsilon}$. If $z\in X$ is a minimizer in (6), called the ground-state quasi-density, then the gradient of $F_{\varepsilon}+v$ at $z$ must be zero,

[TABLE]

Since the regularized functional is differentiable everywhere, the usual problem of $v$-representability is avoided. For $E_{\varepsilon}$, which is still not differentiable, we can resort to the superdifferential. Since any such element $z\in\overline{\partial}E_{\varepsilon}(v)$ automatically solves (7), it is the ground-state quasi-density of the regularized problem with potential $v$. Two important properties of $E_{\varepsilon}$ are (Laestadius et al., 2018, Th. 10 and Cor. 11)

[TABLE]

which relate the regularized problem back to the unregularized one. Since $E$ is already concave, the subtraction of a parabola in (8) makes $E_{\varepsilon}$ strongly concave. Note that the $\varepsilon$ in (9) takes a role comparable to that of permittivity, linking potentials to densities.

To set up a KS scheme we define a reference system that is non-interacting by

[TABLE]

on the same $\tilde{X}\subset X$ and define $E^{0},F_{\varepsilon}^{0},$ and $E_{\varepsilon}^{0}$ analogously. The analogue of (7) for $F_{\varepsilon}^{0}$ at the same quasi-density $z\in X$ is $\nabla F^{0}_{\varepsilon}(z)+v_{\mathrm{KS}}^{\varepsilon}=0$ and defines the KS potential $v_{\mathrm{KS}}^{\varepsilon}$. Simply equating this equation and (7) gives

[TABLE]

where the ground-state quasi-density $z$ and the auxiliary $v_{\mathrm{KS}}^{\varepsilon}$ for the reference system are still unknown and neither $F_{\varepsilon}$ nor $F^{0}_{\varepsilon}$ have a simple, explicit expression. The trick is to determine $z$ and $v_{\mathrm{KS}}^{\varepsilon}$ in an iterative algorithm by replacing them with sequences $x_{i}\to z$, $v_{i}\to v_{\mathrm{KS}}^{\varepsilon}$. The indicated convergence is our major concern in the following proof. We get an update rule for the potential sequence $(v_{i})_{i}$ directly from (11),

[TABLE]

and determine the next quasi-density by solving for the ground-state of the regularized reference system with $v_{i+1}$. This iteration has the stopping condition $v_{i+1}=-\nabla F_{\varepsilon}^{0}(x_{i})$, which means $v=-\nabla F_{\varepsilon}(x_{i})$ by (12). Then $x_{i}$ is already the sought-after ground-state quasi-density $z$ and thus also $v_{i+1}=v_{\mathrm{KS}}^{\varepsilon}$ is the respective KS potential that yields the same quasi-density for the reference system.

The most important ingredient of practical KS calculations enters by giving suitable approximations for the expression $\nabla F_{\varepsilon}-\nabla F^{0}_{\varepsilon}$ (Hartree-exchange-correlation (Hxc) potential including the correlated kinetic energy). For the purpose of showing convergence it is not crucial that this object comes from the exact functional or that it is the result of an approximation, as long as $F_{\varepsilon},F^{0}_{\varepsilon}$ have the stated properties.

The MYKSODA algorithm is then the following. In step (a) get the new potential by (12) above. In step (b) solve the (simpler) ground-state problem for the reference system by choosing the next quasi-density from $\overline{\partial}E_{\varepsilon}^{0}(v_{i+1})$. From (9) it follows that the set of quasi-densities $\overline{\partial}E_{\varepsilon}^{0}(v_{i+1})$ can be determined from the set of ground-state densities of the reference system, which means solving the non-interacting Schrödinger equation. Finally, to ensure a strictly descending energy and to show convergence of $(x_{i})_{i},(v_{i})_{i}$, we include a damping step (c) with an adaptively chosen step length.

MYKSODA iteration scheme.—Assume $\tilde{X}$ bounded and $E^{0}$ finite everywhere. For $v\in X^{*}$ fixed, set $v_{1}=v$ and select $x_{1}\in\overline{\partial}E_{\varepsilon}^{0}(v)$. Iterate $i=1,2,\dots$ according to

1. (a)

set $v_{i+1}=v+\nabla F_{\varepsilon}(x_{i})-\nabla F_{\varepsilon}^{0}(x_{i})$ and stop if $v_{i+1}=-\nabla F_{\varepsilon}^{0}(x_{i})=v_{\mathrm{KS}}^{\varepsilon}$, 2. (b)

select $x_{i+1}^{\prime}\in\overline{\partial}E_{\varepsilon}^{0}(v_{i+1})$ and get the step direction $y_{i}=(x_{i+1}^{\prime}-x_{i})/\|x_{i+1}^{\prime}-x_{i}\|$, 3. (c)

choose the step length $\tau_{i}=-\varepsilon\langle\nabla F_{\varepsilon}(x_{i})+v,y_{i}\rangle>0$ and set $x_{i+1}=x_{i}+\tau_{i}y_{i}$.

We prove below that this algorithm guarantees convergence to the correct KS potential, $v_{i}\to v_{\mathrm{KS}}^{\varepsilon}$, and to the ground-state quasi-density, $x_{i}\to z$, of both the full system with $v$ and the reference system with $v_{\mathrm{KS}}^{\varepsilon}$. The corresponding energy is then determined by $E_{\varepsilon}(v)=F_{\varepsilon}(z)+\langle v,z\rangle>-\infty$. These are still solutions of the regularized problem, but with (8) and (9) a transformation back to the unregularized setting is easily achieved. This, unlike the usually assumed unregularized KS scheme, gives different ground-state densities for the non-interacting and the interacting system, while circumventing all problems of differentiability and thus of $v$-representability. The assumption that $E^{0}$ is finite everywhere is trivially fulfilled in a finite-dimensional setting because $E^{0}$ is a finite sum. It is still kept here to connect more closely to standard DFT and CDFT, where $E^{0}$ is finite even in the infinite-dimensional setting, see Lieb (1983, Th. 3.1(iii)) and Laestadius et al. (2019, Lem. 20), respectively.

Convergence proof.—We refer to the the first part of the proof of Theorem 12 in Laestadius et al. (2018) to show that the superdifferential $\overline{\partial}E_{\varepsilon}^{0}(v_{i+1})$ is everywhere non-empty because of $E^{0}$ finite, guaranteeing that the (regularized) ground-state problem in (b) always has at least one solution. The directional derivative of $F_{\varepsilon}+v$ at $x_{i}$ in direction $x^{\prime}_{i+1}-x_{i}$ can be rewritten by (a),

[TABLE]

Realizing that $x^{\prime}_{i+1}\in\overline{\partial}E_{\varepsilon}^{0}(v_{i+1})$ from (b) and $x_{i}\in\overline{\partial}E_{\varepsilon}^{0}(\nabla F_{\varepsilon}^{0}(x_{i}))$ from invertibility (Laestadius et al., 2018, Lem. 4), we rewrite the right-hand side of (13) with the help of (9), substituting

[TABLE]

which gives

[TABLE]

Now, since $\tilde{x}^{\prime}_{i+1},\tilde{x}_{i}$ are selected from the superdifferential of $E^{0}$ for the respective potentials $v_{i+1},\nabla F_{\varepsilon}^{0}(x_{i})$, the inner product is always smaller or equal to zero (Laestadius et al., 2018, Lem. 5). This property is called monotonicity of $\overline{\partial}E^{0}$ and directly follows from concavity of $E^{0}$. What follows is strong monotonicity of $\overline{\partial}E^{0}_{\varepsilon}$, i.e.,

[TABLE]

The last line follows from (a) and is strictly smaller than zero if not $\|\nabla F_{\varepsilon}(x_{i})+v\|=0$, which would mean that we have already converged to the ground-state quasi-density. We thus infer that, unless converged, we always have a negative directional derivative of $F_{\varepsilon}+v$ at $x_{i}$ in the step direction $y_{i}$ which is parallel to $x^{\prime}_{i+1}-x_{i}$, i.e., $\langle\nabla F_{\varepsilon}(x_{i})+v,y_{i}\rangle<0.$ Such a negative directional derivative means the left leg of the regularization parabola is aligned tangentially to the (differentiable) energy functional $F_{\varepsilon}+v$, like depicted in Fig. 2. The next quasi-density $x_{i+1}=x_{i}+\tau_{i}y_{i}$ is then chosen at the vertex of this regularization parabola.

This corresponds to a choice of step length $\tau_{i}$ where the directional derivatives at $x_{i}$ in direction $y_{i}$ of the regularization parabola $\frac{1}{2\varepsilon}\|\cdot-x_{i+1}\|^{2}$ and of $F_{\varepsilon}+v$ are equal,

[TABLE]

This construction yields a $x_{i+1}\neq x_{i}$, where the energy $e_{i}=F_{\varepsilon}(x_{i})+\langle v,x_{i}\rangle$ is always larger than the energy value $m_{i}$ at the vertex, see Fig. 2. Since the regularization parabola lies fully above the energy functional $F_{\varepsilon}+v$ by construction, the energy $e_{i+1}$ at $x_{i+1}$ must obey $e_{i+1}\leq m_{i}<e_{i}$. The strictly decreasing $e_{i}$ is now by definition bounded below by $E_{\varepsilon}(v)$ from (6) and thus converges. By determining $e_{i}-m_{i}$ from the regularization parabola and then combining it with (18),

[TABLE]

we can infer convergence of $(x_{i})_{i}$. Step (a) then defines an associated potential

[TABLE]

since the gradients are both continuous. After having proved that the densities and potentials converge, it shall be demonstrated that they converge to the expected ground-state quasi-density $z$ and KS potential $v_{\mathrm{KS}}^{\varepsilon}$. We come back to (17), where substituting $x_{i+1}^{\prime}-x_{i}=y_{i}\,\|x_{i+1}^{\prime}-x_{i}\|$ gives

[TABLE]

which together with (18) results in

[TABLE]

We already know from the convergence of densities that $(x_{i})_{i}$ is bounded, further

[TABLE]

by (b) and (9). But $\overline{\partial}E^{0}(v_{i+1})\subset\tilde{X}$, which is bounded, and $(v_{i})_{i}$ converges as well. Thus, $\|x_{i+1}^{\prime}-x_{i}\|$ is bounded, and since $\tau_{i}\to 0$, it follows $\|\nabla F_{\varepsilon}(x_{i})+v\|\to 0$ and $\|\nabla F_{\varepsilon}^{0}(x_{i})+v_{i+1}\|\to 0$. This in turn means $v=-\nabla F_{\varepsilon}(\lim x_{i})$, so $\lim x_{i}=z$ is the ground-state quasi-density for the potential $v$ in the full, regularized problem. Finally, $\lim v_{i+1}=-\lim\nabla F_{\varepsilon}^{0}(x_{i})=-\nabla F_{\varepsilon}^{0}(z)=v_{\mathrm{KS}}^{\varepsilon}$ is the KS potential. $\Box$

As noted above, the reference system reproduces the quasi-density $z$ of the full system and they link back to the real densities by (9),

[TABLE]

where typically $\rho_{\mathrm{KS}}^{\varepsilon}\neq\rho$. Then $v_{\mathrm{KS}}^{\varepsilon}-v=\varepsilon^{-1}(\rho_{\mathrm{KS}}^{\varepsilon}-\rho)$ is precisely the Hxc potential that depends on the regularization parameter $\varepsilon$ here. This means every choice of $\varepsilon$ defines a different reference system. A limit $\varepsilon\to 0$ in the algorithm is unfeasible because of its relation to the step length.

A simulation of two electrons on a ring lattice MYr allows us to illustrate the above method. Compared to a previous implementation in a CDFT setting Laestadius et al. (2019), the version given here uses the more conservative damping step that helped to prove convergence. To distinguish the two versions, we denote them “MYKSODA-S” for shorter, conservative steps, and “-L” for the original longer steps Laestadius et al. (2018, 2019). Both versions have been adapted to a pure DFT setting, $H_{\mathrm{kin}}$, taking the form of a standard second-order finite difference. A radius of $R=1$ bohr, a uniform grid with 30 points, and the interaction energy $H_{\mathrm{int}}=3\sqrt{1+\cos(\theta_{1}-\theta_{2})}$ were used. As expected, larger $\varepsilon$ leads to faster convergence. Also, the more conservative steps taken by MYKSODA-S often lead to slower convergence in practice. Surprisingly, however, in some cases MYKSODA-S overtakes the less conservative MYKSODA-L. An example is shown in Fig. 3. Such a crossover is possible as the two algorithms follow different paths through the space of densities and potentials. Yet, when the starting point is the same, the first step by MYKSODA-L always lowers the energy more than the first step by MYKSODA-S. Although it is a plausible conjecture that also MYKSODA-L, taking maximal steps, is guaranteed to converge, the present proof does not establish this. “Maximal steps” here means taking $\tau_{i}$ maximally such that $\langle\nabla F_{\varepsilon}(x_{i+1})+v,x^{\prime}_{i+1}-x_{i}\rangle\leq 0$, which yields maximal decrease in energy in the direction chosen by step (b).

In this Letter we proved convergence of the regularized KS scheme with special adaptive damping. In short, this means that KS-DFT is a veritable method to calculate the correct ground-state density. This strong statement holds for all flavors of DFT that are defined on a finite-dimensional density space $X=\ell^{2}(M)$ and have a linear coupling to external potentials of type $\langle v,\rho\rangle$. This includes CDFT, where the potential $v$ is a combination of scalar and vector potential, and the density $\rho$ includes the paramagnetic current density. To allow for a combination of these different entities into one Banach space setting, the respective function spaces for one-particle densities and current densities have to fulfill a condition termed “compatibility” in Laestadius et al. (2019). A proof of MYKSODA convergence for infinite-dimensional Banach spaces $X$ is feasible but much more technical and will be presented elsewhere. The choice of step length in (c) is an essential part of the proof and similar choices could be of value in showing convergence of related iteration schemes and for other settings such as Hartree–Fock theory Cancès and Le Bris (2000). Next to this damping step the MY regularization is a vital part of the proof at hand, not only to have functional differentiability, but also for the strong monotonicity estimate needed to show convergence. How those findings can be transferred to realistic KS implementations will be the content of future research, but it is expected that they serve as useful guidelines for better convergence results.

Erratum et Corrigendum

Additional author: Paul E. Lammert, Department of Physics, Pennsylvania State University, USA

The central argument for convergence of the “MYKSODA iteration scheme” introduced in our article Penz et al. (2019) contains a critical mistake that we hereby want to correct. The claim was that because of $\frac{1}{2\varepsilon}\|x_{i+1}-x_{i}\|^{2}\leq e_{i}-e_{i+1}\to 0$ in Eq. (19) the sequence of quasidensities $(x_{i})_{i}$ also converges. This is wrong, since this does not show the necessary Cauchy property $\|x_{m}-x_{n}\|\to 0$ for all $m,n\to\infty$. Convergence could be secured, on the other hand, in cases where $e_{i}-e_{i+1}\leq r^{i}$ with $0<r<1$, because then $\|x_{i+1}-x_{i}\|\leq C\sqrt{r}^{i}$ and summability of the geometric series would give the Cauchy property for $(x_{i})_{i}$. Yet, such a condition for the convergent energy sequence $(e_{i})_{i}$ is in general not available. The corrigendum given below proceeds along a different path that effectively takes degeneracies of the ground state into account. First, it is shown that the sequence of quasidensities $(x_{i})_{i}$ is bounded. Since we operate in a finite-dimensional setting this means by the Bolzano–Weierstrass theorem that the sequence has at least one convergent subsequence. We then show that all such subsequences converge to ground-state quasidensities of the regularized, fully interacting problem. Thereby, the accumulation points of the sequence give the desired solutions of which there can be more than one in case of degeneracy. Although in this way the algorithm can accommodate multiple valid solutions, we cannot be sure that we map out all possible solutions. The selection of specific solutions might still depend on the initial value $x_{0}$.

While the MYKSODA iteration scheme remains precisely the same, the statement following it has to be modified. We prove that for any convergent subsequence $x_{\alpha(i)}\to z$ the corresponding potential sequence $v_{\alpha(i)}$ converges to the correct KS potential that reproduces $z$ in the regularized reference problem. We will now repeat the last part of the proof which needs to be adjusted.

Corrected convergence proof. Up to Eq. (19) all statements stay intact, just the inference of convergence right after Eq. (19) is erroneous. The $F_{\varepsilon}$ increases like $\|x\|^{2}$ asymptotically by construction (cf. the regularization procedure in Eq. (5)) because $F(x)=+\infty$ whenever $x$ is outside the bounded $\tilde{X}$. Since $e_{i}=F_{\varepsilon}(x_{i})+\langle v,x_{i}\rangle\leq e_{1}$ for all $i$ this means that $(x_{i})_{i}$ must already be bounded. This translates to boundedness of $(v_{i})_{i}$ defined by step (a), because $\nabla F_{\varepsilon},\nabla F^{0}_{\varepsilon}$ are Lipschitz-bounded (with constant $\varepsilon^{-1}$; Barbu and Precupanu (2012, Corollary 2.59)). We can then argue for $\|\nabla F_{\varepsilon}(x_{i})+v\|\to 0$ and $\|\nabla F^{0}_{\varepsilon}(x_{i})+v_{i+1}\|\to 0$, exactly like in the original proof. This means $v=-\lim\nabla F_{\varepsilon}(x_{i})$ and for any accumulation point $z$ of $(x_{i})_{i}$ (Bolzano–Weierstrass guarantees that there is at least one) there must be a convergent subsequence $x_{\alpha(i)}\to z$ such that $v=-\lim\nabla F_{\varepsilon}(x_{\alpha(i)})=-\nabla F_{\varepsilon}(z)$. By inversion we get $z\in\overline{\partial}E_{\varepsilon}(v)$ which means that $z$ is the correct ground-state quasidensity of the regularized full problem. Finally, $\lim v_{\alpha(i+1)}=-\lim\nabla F^{0}_{\varepsilon}(x_{\alpha(i)})=-\nabla F^{0}_{\varepsilon}(z)=v_{\mathrm{KS}}^{\varepsilon}$ is the KS potential that exactly reproduces this $z$ in the regularized reference problem. A different accumulation point with its associated convergent subsequence $(x_{\beta(i)})_{i}$ will lead to another possible ground-state quasidensity of the regularized full problem in the case of degeneracy. Then the procedure assigns it the appropriate KS potential as the limit of $(v_{\beta(i)})_{i}$ that reproduces the same quasidensity in the regularized reference problem. In this manner the whole sequence $(x_{i})_{i}$ partitions into convergent subsequences that have ground-state quasidensities as limit points. $\Box$

Note that it is still true that the distance between successive terms goes to zero, i.e. $\|x_{i+1}-x_{i}\|\to 0$. This means that no isolated accumulation points are possible. But since the ground-state solutions in terms of density matrices $\Gamma$ form a convex set and the mapping $\Gamma\mapsto\rho$ is linear, also the set of ground-state (quasi)densities will be convex and thus corresponds to a whole ‘accumulation region’. The sequence $(x_{i})_{i}$ will thus converge to this set of solutions. If the ground-state density is non-degenerate (unique) then $\overline{\partial}E_{\varepsilon}(v)$ is single valued and the whole sequence converges to this point only.

Another small mistake should be noted: After (13), two times in (15), and after (16) there must be a minus sign in front of $\nabla F^{0}_{\varepsilon}$. The error does not influence any subsequent results. We would also like to point out a misleading feature of Fig. 1, where the impression is given that the minima of the one-dimensional sections through $F+v$ and $F_{\varepsilon}+v$ coincide, which holds true in the special case of the absolute minima but not along an arbitrary direction.