1 Introduction↩︎

Kohn-Sham density functional theory (KS-DFT) is a widely used method in solid-state physics and chemistry for simulating the electronic properties of materials [1], [2]. The Kohn-Sham equations are typically solved using self-consistent field (SCF) iterations, and improving the convergence of these iterations has been an ongoing challenge since the first numerical resolution of these equations.

The ground-state density in Kohn-Sham DFT is the solution of a nonlinear equation \(\mathcal{F}(\rho) = 0\), where \(\mathcal{F}\) is the fixed-point residual. The Jacobian of \(\mathcal{F}\) is given by the adjoint dielectric matrix \(\varepsilon^\dagger = 1 - \chi_0 K_\mathrm{HXC}\), with \(\chi_0\) denoting the non-interacting susceptibility and \(K_\mathrm{HXC}\) the Hartree-exchange-correlation kernel. Methods to accelerate the convergence of SCF iterations generally involve constructing approximations of this Jacobian, to apply a quasi-Newton method. Two distinct approaches exist to build these approximations, and in practice, both are combined.

The first approach consists of generic methods that do not rely on prior knowledge of \(\mathcal{F}\). Instead, they use the iteration history to construct approximate Jacobians. In the context of SCF iterations in Kohn-Sham DFT, these methods are referred to as mixing. Some widely used mixing schemes include Broyden’s method [3], [4], Anderson acceleration [5], or Pulay mixing [6]. The development and refinement of these algorithms have been an active area of research, with numerous variants and improvements proposed in the literature (e.g., [7]–[15]).

The second approach is KS-DFT-specific methods that use simple physical models or knowledge of the structure of \(\mathcal{F}\) to construct estimates of the adjoint dielectric matrix. These methods are known as preconditioning techniques. The earliest proposition of a preconditioner for Kohn-Sham DFT is the Kerker preconditioner [16], which approximates the dielectric response with the Thomas-Fermi screening function. Another widely used preconditioner is the Resta preconditioner that uses the model dielectric function for semiconductors from [17], also based on the Thomas-Fermi model.

Many improvements to these two preconditioners have been proposed since. The Kerker and Resta preconditioners, initially designed for plane-wave calculations, have been adapted for real-space implementations or other more complex frameworks [18]–[20]. Significant efforts have also been directed toward developing preconditioners suitable for inhomogeneous systems. For instance, Lin and Yang [21] proposed to use an elliptic partial differential equation with inhomogeneous coefficients, while Zhou et al. [22] proposed a generalization of the Kerker preconditioner for mixed systems. In [23], Herbst and Levitt developed a parameter-free preconditioner based on the local density of states. Alternative strategies include the use of the full dielectric matrix computed on reduced systems [24], [25]. Beyond linear preconditioners, nonlinear approaches also exist. These involve solving a simplified problem at each SCF iteration, such as the Thomas-Fermi-von Weizsäcker equation for inhomogeneous systems [26], a second-variational prediction operator [27], or an auxiliary bosonic system [28].

Dederichs and Zeller [29] and later Herbst and Levitt [23] identified three sources of slow convergence related to the dielectric operator of the system: charge sloshing in large simulation cells caused by the long range of the Coulomb interaction, localized orbitals that yield large contributions in the non-interacting susceptibility, and lastly, phase transitions associated with a non-invertible dielectric matrix, typically in magnetic systems. Most of the available preconditioners are designed to mitigate charge sloshing issues. In this work, we focus on the convergence issue associated with magnetic phase transitions. We propose a hybrid preconditioner \(P^\mathrm{hybrid}= 1 - \chi_0^\mathrm{LDOS}K_H - \chi_0^\mathrm{diag}K_\mathrm{XC}\), which extends the LDOS preconditioner of [23] by incorporating an additional term \(\chi_0^\mathrm{diag}\) in which we neglect the orbital variations. This approach is inspired by the Stoner model, to address convergence challenges near ferromagnetic phase transitions.

This paper is organized as follows: section 2 reviews Kohn-Sham DFT, specifically its magnetic formulation in the collinear-spin DFT formalism, and discusses the convergence properties of SCF iterations. Section 3 details the specific convergence issues related to magnetism. Section 4 introduces our proposition for a preconditioning scheme, inspired by the Stoner model, to reduce convergence issues specific to magnetic systems. Lastly, section 5 presents numerical results for our preconditioner.

2 The self-consistent field iterations in spin-polarized Kohn-Sham DFT↩︎

Consider a system of \(N\) electrons in a domain \(\Omega\) evolving in a potential \(V_\mathrm{ext}\) created by the nuclei and the core electrons. In the Kohn-Sham density functional theory (KS-DFT) framework, these electrons are treated as non-interacting particles in a density-dependent effective potential \(V_\mathrm{HXC}= V_H + V_\mathrm{XC}\). To allow for spin-polarization, we work within spin-DFT, where the exchange-correlation potential \(V_\mathrm{XC}\) depends not only on the total electron density but also on its spin-resolved components. A common restriction is the collinear spin approximation, in which the exchange-correlation potential acts independently on each spin channel. Under this assumption, the spin, \(\alpha \in \{+, -\}\), becomes a good quantum number, leading to the following simplified formulation of Kohn-Sham equations for collinear-spin-DFT [30], [31]: \[\begin{align} \begin{aligned} & \left[ -\frac{1}{2} \nabla^2 + V_\mathrm{ext} + V_\mathrm{HXC}^\alpha(\rho^+, \rho^-) \right] \psi_i^\alpha = \epsilon_i^\alpha \psi_i^\alpha \\ &\int_\Omega \overline{\psi_i^\alpha} \psi_j^\alpha = \delta_{ij} \\ & \rho^\alpha = \sum_{i=1}^\infty f(\epsilon_i^\alpha - \epsilon_F) |\psi_i^\alpha|^2 \\ & V_\mathrm{HXC}^\alpha(\rho^+, \rho^-) = v_c(\rho^++\rho^-) + V_\mathrm{XC}^\alpha(\rho^+, \rho^-) \\ & \int_\Omega \rho^+ + \rho^- = N \end{aligned} \label{eq:KS} \end{align}\tag{1}\] for \(i \geq 1\) and \(\alpha \in \{+, -\}\). Here \(f(x) = 1/(1+\exp(x/T))\) is the Fermi-Dirac function at smearing temperature \(T\), \(v_c\) is the Coulomb kernel, and \(V_\mathrm{XC}^\alpha\) is the exchange-correlation potential for spin channel \(\alpha \in \{+, -\}\). We will denote \(\rho = \begin{pmatrix} \rho^+\\ \rho^- \end{pmatrix}\) and \(V = \begin{pmatrix} V^+\\ V^- \end{pmatrix}\) the spin-resolved densities and potentials.

Following the formalism used in [23], we define the potential-to-density mapping \(F\) by \[\begin{align} F^\alpha\left(V\right) = \sum_{i=1}^\infty f(\epsilon_i^\alpha - \epsilon_F) |\psi_i^\alpha|^2 \end{align}\] where, for \(\alpha \in \{+, -\}\), \((\epsilon_i^\alpha, \psi_i^\alpha)\) are the normalized eigenpairs of \(-\frac{1}{2} \nabla^2 + V_\mathrm{ext} + V^\alpha\) and the Fermi level \(\epsilon_F\) is adjusted so that \(\sum_{\alpha \in \{+, -\}} \sum_{i=1}^\infty f(\epsilon_i^\alpha - \epsilon_F) = N\). This allows us to rewrite the Kohn-Sham equations 1 as the following fixed-point, or self-consistent, equation: \[\begin{align} F(V_\mathrm{HXC}(\rho)) = \rho . \label{eq:KSFP} \end{align}\tag{2}\]

This equation is often solved using self-consistent field (SCF) iterations. One of the simplest type of SCF iterations is damped iterations, where the density is updated with \(\rho_{n+1} = \rho_n + \tau [ F(V_\mathrm{HXC}(\rho_n)) - \rho_n ]\) for some damping parameter \(0<\tau \leq 1\). As it is standard (see [23], [32]), the linearization of these iterations, around a solution \(\rho^*\), with an optimal \(\tau\) yields a convergence rate \(R = (\kappa_s(\varepsilon^\dagger)-1)/(\kappa_s(\varepsilon^\dagger)+1)\) where \(\kappa_s(\varepsilon^\dagger) = \lambda_{\max}(\varepsilon^\dagger)/\lambda_{\min}(\varepsilon^\dagger)\) is the spectral condition number, i.e., the ratio of the largest to the smallest eigenvalues, of the adjoint dielectric matrix \[\begin{align} \varepsilon^\dagger = 1 - \chi_0 K_\mathrm{HXC}. \end{align}\] Here, \(\chi_0 = F'(V_\mathrm{HXC}(\rho^*))\) is the non-interacting susceptibility, and \(K_\mathrm{HXC}= V_\mathrm{HXC}'(\rho^*)\) denotes the Hartree-exchange-correlation kernel. Because \(\chi_0\) is self-adjoint and non-positive, \(\varepsilon^\dagger\) has a real spectrum and its spectral condition numbers is well-defined (see [23]).

Within the collinear-spin formalism, these derivatives can be expressed as \(2 \times 2\) Jacobian matrices of operators: \[\begin{align} \chi_0 = \begin{pmatrix} \chi_0^{++} & \chi_0^{+-} \\ \chi_0^{-+} & \chi_0^{--} \end{pmatrix} \end{align}\] where the matrix elements are derived within perturbation theory [23], [33]–[35].

The Hartree-exchange-correlation kernel can be decomposed into two contributions, \(K_\mathrm{HXC}= K_H + K_\mathrm{XC}\), with \[\begin{align} K_H = \begin{pmatrix} v_c & v_c \\ v_c & v_c \end{pmatrix} . \end{align}\]

3 Convergence issues in magnetic systems↩︎

The adjoint dielectric matrix \(\varepsilon^\dagger = I - \chi_0 (K_H + K_\mathrm{XC})\) is ill-conditioned when it has either large or small eigenvalues. Since the non-interacting susceptibility \(\chi_0\) is negative definite and the Hartree kernel \(K_H\) is positive semi-definite, the operator \(I-\chi_0 K_H\) has all its eigenvalues larger than or equal to 1. Hence, convergence issues associated with small eigenvalues of \(\varepsilon^\dagger\) originate from strong exchange–correlation effects [23]. This is typically the case in magnetic systems, especially near phase transitions.

Indeed, when a non-zero magnetization appears continuously as a function of one of the system parameters, \(\gamma\), this transition must be associated with a non-invertible dielectric matrix. Due to symmetries, spin-polarized solutions of the Kohn-Sham equation are degenerate. However, if the dielectric matrix were invertible, the implicit function theorem would guarantee a locally unique solution to the Kohn-Sham equation with respect to \(\gamma\), thereby preventing the emergence of multiple distinct solutions at this point.

In this section, we first discuss quantitatively how a few small eigenvalues can affect the convergence of SCF iterations, especially with advanced mixing strategies. Then, we study magnetic phase transitions, first with the Stoner model and then in Kohn-Sham DFT.

3.1 Convergence issues associated with small eigenvalues↩︎

When the Kohn-Sham fixed-point equation is solved with optimally damped iterations, the convergence rate is given by \(R = (\kappa_s-1)/(\kappa_s+1)\). For large values of \(\kappa_s\), this rate can be approximated by \(R \approx 1-2/\kappa_s\) and the number of iterations required to reach a tolerance \(e\) is \(n_\mathrm{ite} \approx \kappa_s \frac{1}{2} \log(1/e)\). This implies that if the smallest eigenvalue of \(\varepsilon^\dagger\) is reduced by a given factor, the number of iterations increases by this same factor. Thus, small outlier eigenvalues of \(\varepsilon^\dagger\) severely degrade the convergence. However, damped iterations are rarely used in practice, as more sophisticated and more efficient fixed-point solvers are available. One such method is Anderson acceleration, whose convergence properties and limitations are discussed in [36].

Anderson acceleration, when applied to a linear problem with an untruncated history, is equivalent to the Generalized Minimal Residual (GMRES) method [37], a Krylov subspace technique (see, e.g., [38]). This section therefore focuses on the impact of small eigenvalues on the convergence of GMRES, for a matrix \(A = \varepsilon^\dagger\).

Krylov subspace methods [39], [40] construct iterative approximate solutions \(x_n\) to the linear problem \(Ax=b\), within the subspace \(x_0 + \mathcal{K}_n(A, r_0)\), where \(\mathcal{K}_n(A, r_0)\) is the order-\(n\) Krylov subspace generated by \(A\) and the initial residual \(r_0\). Equivalently, Krylov methods can be viewed as constructing approximations whose residuals can be expressed as \[\begin{align} r_n = P_n(A) r_0 \end{align}\] where \(P_n\) is a polynomial of degree \(n\) that satisfies \(P_n(0)=1\). We call the space of such polynomials \(\mathcal{P}_n\).

Within this framework, the GMRES method minimizes the \(L_2\)-norm of the residual, \(\|r_n\|\), over all polynomials in \(\mathcal{P}_n\) [39], that is: \[\begin{align} \|r_n\| = \min_{P \in \mathcal{P}_n} \|P(A) r_0\| . \end{align}\] In the linearization of the SCF procedure, the matrix \(A\) is the adjoint dielectric matrix \(\varepsilon^\dagger\), which is diagonalizable with a non-negative real spectrum. This gives the following bound: \[\begin{align} \|r_n\| \leq \kappa(V) \|r_0\| \inf_{P \in \mathcal{P}_n} \max_{\lambda \in \Lambda(A)} |P(\lambda)| \end{align}\] where \(\kappa(V) = \|V\|\|V^{-1}\|\) is the condition number of \(V\), the eigenvector matrix of \(A\).

Chebyshev polynomials provide a classical bound for the above minimization problem when the spectrum of \(A\) is contained in an interval \([\lambda_{\min} , \lambda_{\max}]\) with \(0 < \lambda_{\min}\) [39]. In this case, the minimum \[\begin{align} \min_{P \in \mathcal{P}_n} \max_{\lambda \in [\lambda_{\min}, \lambda_{\max}]} |P(\lambda)| \end{align}\] is achieved by a scaled and shifted Chebyshev polynomial of degree \(n\). It yields the following bound for the residual: \[\begin{align} \|r_n\| \leq 2 \kappa(V) \|r_0\| \left( \frac{\sqrt{\kappa_s} - 1}{\sqrt{\kappa_s} + 1} \right)^n \end{align}\] with \(\kappa_s = \lambda_{\max}/\lambda_{\min}\) being the spectral condition number of \(A\). This bound is relevant when the spectrum of \(A\) is uniformly distributed over the interval \([\lambda_{\min}, \lambda_{\max}]\). However, if \(A\) has outlier eigenvalues, this bound becomes overly pessimistic, as polynomials do not need to attenuate the entire interval \([\lambda_{\min}, \lambda_{\max}]\).

The impact of spectral outliers on the convergence of Krylov methods is nicely explained in [40] from a polynomial approximation perspective.

Consider a symmetric matrix \(A\) with a real spectrum including a small outlier eigenvalue \(\lambda_\mathrm{out}\):

We denote \(\tilde{A}\) the restriction of \(A\) to the eigensubspace associated with the spectrum of \(A\) excluding \(\lambda_\mathrm{out}\). Suppose the GMRES method is applied to solve a linear system associated with \(\tilde{A}\), yielding a sequence of polynomials \((P_n)\) and residuals \(\tilde{r}_n = P_n(A) \tilde{r}_0\) that achieve a convergence rate \(\gamma\). We consider here that it means \(\| P_n(A) \tilde{r}_0 \| \approx \gamma^n \|\tilde{r}_0\|\). Now consider the two polynomials \(P_n\) and \(Q_n(X) = (1-X/\lambda_\mathrm{out})P_{n-1}(X)\). These give the following bound for the residual \(r_n\) of the GMRES method applied to \(A\) with initial residual \(r_0 = s_0 + \tilde{r}_0\), where \(s_0\) is along the eigenvector associated with \(\lambda_\mathrm{out}\): \[\begin{align} &\|r_n\| \leq |P(\lambda_\mathrm{out})| \|s_0\| + \| P_n(A) \tilde{r}_0\| \approx \|s_0\| + \gamma^{n} \|\tilde{r}_0\| \tag{3} \\ &\|r_n\| \leq \|Q_n(A) \tilde{r}_0\| \approx \frac{\lambda_{\max}}{\lambda_\mathrm{out}} \gamma^{n-1} \|\tilde{r}_0\| \tag{4} . \end{align}\]

These two bounds are illustrated in Figure 1. They suggest that the convergence curves for this type of system initially follows the convergence rate of the matrix without outliers, until the residual along the mode associated with \(\lambda_\mathrm{out}\) dominates, and then exhibits a plateau of \(n_\mathrm{plateau}\) iterations, before aligning again with the convergence of the matrix without outliers. Here \(n_\mathrm{plateau}\) is the number of additional iterations needed for the residual on the remaining spectrum in the second bound (equation 4 ) to match the residual on the remaining spectrum in the first bound (equation 3 ), that is \[\begin{align} \gamma^{n_\mathrm{plateau}-1} \approx \frac{\lambda_\mathrm{out}}{\lambda_{\max}} . \end{align}\] This means that the number of additional steps introduced by the small outlier \(\lambda_\mathrm{out}\) corresponds to the number of iterations needed to improve the convergence by a factor \(\frac{\lambda_\mathrm{out}}{\lambda_{\max}}\) for the matrix \(A\).

This argument can be readily extended to the case of a matrix with multiple outliers, where the length of the plateau will correspond to the sum of the above estimates for each (effective) outlier.

Thus, while Anderson acceleration reduces the impact of small outlier eigenvalues of the dielectric matrix, it does not eliminate it entirely. The number of required additional iterations scales logarithmically with the outlier; with a convergence that is all the more affected when outliers are numerous and when the underlying matrix (without outliers) is ill-conditioned. This analysis also provides a diagnostic for convergence issues associated with a plateau, linking them to the presence of small eigenvalues in the adjoint dielectric matrix.

3.2 Stoner model↩︎

The Stoner model provides a simple framework to understand convergence issues associated with magnetic transitions, as discussed in [29].

The Stoner model for ferromagnetism [41], [42] gives a scalar description for the total magnetization \(M\) of a system. In this model, we consider that the total magnetization \(M\) rigidly shifts a prescribed density of states by \(IM\) between the \(+\) and \(-\) spin channels. The Stoner parameter \(I\) characterizes the strength of the exchange interaction in the system. This leads to the following fixed-point equation for \(M\): \[\begin{align} G(M) = M \label{eq:stoner} \end{align}\tag{5}\] with \[\begin{align} \begin{aligned} G(M) = & \int_{-\infty}^{\epsilon_F(M)} D\left(\epsilon + \frac{1}{2} IM\right) d\epsilon \\ & - \int_{-\infty}^{\epsilon_F(M)} D\left(\epsilon - \frac{1}{2} IM\right) d\epsilon . \end{aligned} \end{align}\] where \(D\) is the density of states of the non-spin-polarized system, and \(\epsilon_F(M)\) is determined by the constraints on the total number of electrons, \(\int_{-\infty}^{\epsilon_F(M)} [ D\left(\epsilon + 1/2 IM\right) + D\left(\epsilon - 1/2 IM\right) ] d\epsilon = N\).

Equation 5 has a trivial solution \(M_0=0\). Since \(M\) is bounded by the total number of electrons, a sufficient condition for a magnetic solution to exist is \(G'(0) > 1\). This condition rewrites \(ID(\epsilon_F^0) > 1\) and is the well-known Stoner criterion.

The fixed-point equation 5 can be interpreted as the stationary condition for the following Stoner energy: \[\begin{align} \begin{aligned} E_\mathrm{Stoner}(M) = &\int_{-\infty}^{\epsilon_F(M)+\lambda(M)} \epsilon D\left(\epsilon+\frac{I}{2}M\right) d\epsilon \\ &+ \int_{-\infty}^{\epsilon_F(M)-\lambda(M)} \epsilon D\left(\epsilon-\frac{I}{2}M\right) d\epsilon \\ &+ \frac{1}{4} IM^2 \end{aligned} \end{align}\] in which \(\epsilon_F(M)\) and \(\lambda(M)\) are determined by the condition on the total number of electrons and total magnetization: \[\begin{align} &\begin{aligned} N = &\int_{-\infty}^{\epsilon_F(M) + \lambda(M)} D\left(\epsilon + 1/2 IM\right) d\epsilon \\ &+ \int_{-\infty}^{\epsilon_F(M) - \lambda(M)}D\left(\epsilon - 1/2 IM\right)d\epsilon \\ \end{aligned} \\ &\begin{align} M = &\int_{-\infty}^{\epsilon_F(M) + \lambda(M)} D\left(\epsilon + 1/2 IM\right) d\epsilon \\ &- \int_{-\infty}^{\epsilon_F(M) - \lambda(M)}D\left(\epsilon - 1/2 IM\right)d\epsilon \end{align} . \end{align}\] The stationary condition imposes \(\lambda(M)=0\), which recovers the Stoner fixed-point equation. This energy gives a stability condition for the non-spin-polarized solution \(M_0=0\), which is \(E_\mathrm{Stoner}''(0) = (2D(\epsilon_F^0))^{-1}(1-ID(\epsilon_F^0)) > 0\). This is not fulfilled when the Stoner criterion is satisfied. In this case, the system will spontaneously polarize in spin.

In order to understand more precisely how the magnetic solution appears, we use a perturbation expansion of \(G(M)\) around 0, which allows us to rewrite equation 5 as \[\begin{align} 0 = M(ID(\epsilon_F^0)-1 + a M^2) + \mathcal{O}(M^4) \label{eq:stoner95lin} \end{align}\tag{6}\] with \(a = 1/24 I^3 D''(\epsilon_F^0)\). The sign of \(a\) determines two distinct scenarios.

For \(a < 0\) (i.e., a concave density of states at the Fermi level), equation 6 has two non-zero solutions \(M_\pm = \pm (\sqrt{-a})^{-1}\sqrt{ID(\epsilon_F^0)-1} + \mathcal{O}({(M_\pm)}^2)\) when the Stoner criterion is met, which are actual minima of the energy, i.e., \(E_{\mathrm{Stoner}}''(M_\pm) + \mathcal{O}({(M_\pm)}^3) > 0\). The magnetic solution therefore appears continuously with respect to \(I\), so the system exhibits a second-order magnetic transition.

For \(a >0\) (i.e., a convex density of states at the Fermi level), equation 6 has two non-zero solutions \(M_\pm = \pm (\sqrt{a})^{-1}\sqrt{1 - ID(\epsilon_F^0)} + \mathcal{O}({(M_\pm)}^2)\) when the Stoner criterion is not met and those solutions are maxima of the Stoner energy (\(E_\mathrm{Stoner}''(M_\pm) + \mathcal{O}({(M_\pm)}^3) < 0\)). The linearization of equation 5 does not give a magnetic solution when the Stoner criterion is met, but in this case, the Stoner fixed-point equation has at least one magnetic solution. However, in this case, the magnetic solution doesn’t appear continuously with respect to \(I\), resulting in a first-order magnetic transition.

Figure 2 shows an example of the magnetization predicted by the Stoner model for a given density of states in the case \(a<0\), that yields a second-order magnetic transition. We have also represented a few energy surfaces that illustrate the correspondence between the local extrema and the predicted magnetization.

The scalar fixed-point equation \(M=G(M)\) can be solved using fixed-point iterations. The convergence of this method is determined by the derivative of \(G\) at the fixed point (see [29]), and the analogue of the dielectric matrix in this simple model is \(\varepsilon_\mathrm{Stoner}(M) = 1 - G'(M)\). At the Stoner criterion (\(ID(\epsilon_F^0) = 1\)), we have \(\varepsilon_\mathrm{Stoner}(0)=0\); thus, for systems exhibiting a second-order magnetic transition, convergence can be arbitrarily slow near this point. For systems exhibiting a first-order magnetic transition, the transition occurs before reaching the Stoner criterion, which mitigates convergence issues.

3.3 Magnetic transitions in Kohn-Sham DFT↩︎

In Kohn-Sham DFT, we also expect second-order magnetic phase transitions to be associated with zero eigenvalues of the dielectric matrix.

Consider a parameter-dependent Kohn-Sham equation of the form \[\begin{align} \mathcal{F}(\rho, \gamma) = 0 \end{align}\] where \(\mathcal{F}(\rho, \gamma)\) represents the fixed-point residual of \(\rho\) at the parameter \(\gamma\) (e.g. pressure, volume, electronic temperature…). The derivative of \(\mathcal{F}\) with respect to \(\rho\) at a solution \(\rho^*(\gamma)\) is the adjoint dielectric matrix \(\varepsilon^\dagger(\gamma)\). As previously discussed, a magnetic solution can only appear continuously with respect to \(\gamma\) at points where \(\varepsilon^\dagger(\gamma)\) is not invertible. Assume that \(\rho^*(0)\) is non-magnetic and that its associated adjoint dielectric matrix \(\varepsilon^\dagger(0)\) has a zero eigenvalue associated with the spin mode \(\rho_1\).

The bifurcation of the solution at \(\gamma=0\) can be analyzed using a perturbation expansion of \(\mathcal{F}\) around \((\rho^*(0), 0)\), and a Lyapunov-Schmidt reduction. We denote \(F_1 = \mathrm{Ker}(\varepsilon^\dagger(0)) = \mathrm{Span}(\rho_1)\), \(F_2 = \mathrm{Ran}(\varepsilon^\dagger(0))\) and \(P_{F_1}\), \(P_{F_2}\) the corresponding projectors. The equation \(\mathcal{F}(\rho, \gamma) = 0\) can be decomposed into: \[\begin{align} & P_{F_1} \mathcal{F}(\rho(0) + t\rho_0 + \rho_2, \gamma) = 0 \\ & P_{F_2} \mathcal{F}(\rho(0) + t\rho_0 + \rho_2, \gamma) = 0 \end{align}\] where \(\rho_2 \in \mathcal{F}_2\) and \(t \in \mathbb{R}\). For \(t\) and \(\gamma\) given, the implicit function theorem ensures that the second equation has a unique solution near \((t=0, \gamma=0)\), denoted here \(\rho_2(t, \gamma)\). Substituting this into the first equation reduces it to an equation in \(t\) and \(\gamma\) alone, \(f(t, \gamma)=0\), with \(\frac{d}{dt}f(0, 0) = \braket{\rho_1}{ \varepsilon^\dagger(0) \rho_1} = 0\). Using the spin inversion symmetries, we can show that \(\frac{d}{d\gamma}f(0, 0) = 0\) and \(\frac{d^2}{dt^2}f(0, 0) = 0\). This leads to an equation analogous to that obtained in the perturbation expansion of the Stoner model: \[\begin{align} 0 = t(a\gamma t + bt^2) + \mathcal{O}(t^4 + (\gamma t)^2) \end{align}\] with \(a=\frac{d^2}{dtd\gamma}f(0, 0)\) and \(b=\frac{d^3}{dt^3}f(0, 0)\). We recognize a pitchfork bifurcation, where the magnetization emerges as the square root of the parameter \(\gamma\), with \(t = \pm \sqrt{-(a/b) \gamma} + \mathcal{O}(\gamma)\) provided that \(-a/b > 0\).

Besides, note that \(\frac{d}{dt}\rho_2(0, 0) = 0\), implying that the magnetic solution writes \[\begin{align} \rho^*(\gamma) = \rho^*(0) \pm \sqrt{- (a/b) \gamma} \rho_1 + \mathcal{O}(\gamma) . \end{align}\] This shows that, at the transition, the spin component of the density aligns with the eigenvector associated with the zero eigenvalue of the adjoint dielectric matrix.

In order to illustrate this type of transition and the associated convergence issues, we have done some numerical experiments on the ferromagnetic transition of BCC iron. The calculations were performed with the density-functional toolkit (DFTK) [43]. We have used a plane-wave basis with an energy cutoff of 50 Hartree and a Monkhorst-Pack k-points grid of size \(20\times20\times20\), with norm-conserving pseudopotentials [44], the Perdew-Burke-Ernzerhof (PBE) [45] exchange-correlation functional, and the Fermi-Dirac occupation scheme with a temperature of 0.01 Hartree. To induce the magnetic transition, we varied the lattice constant around a reference lattice constant for BCC iron (2.46 Å). Note that this transition is not the physical transition of iron under pressure, as a phase transition to HCP iron occurs before the magnetic moment of BCC iron disappears, around 13 GPa [46], [47].

Figure 3 shows the magnetization as a function of the lattice constant for this BCC iron system, as well as the smallest eigenvalue and condition number of the adjoint dielectric matrix \(\varepsilon^\dagger\). We observe that the smallest eigenvalue of the dielectric matrix decreases as the lattice constant increases, eventually reaching zero and triggering a phase transition. At this transition, the spectral conditioning of the dielectric matrix blows up, which should cause convergence issues.

These convergence issues are illustrated in Figure 4, which shows the number of SCF iterations needed to converge the density of BCC iron for various lattice constants. An increase in the number of iterations is observed around the transition, where the adjoint dielectric matrix becomes singular. We also include the predicted number of iterations derived in section 3.1, \(n_\mathrm{prediction} = 1 + \log(2 \cdot 10^{-10}\lambda_1/\lambda_M) / \log(\gamma)\), where \(\gamma = (\sqrt{\lambda_M/\lambda_2}-1)/(\sqrt{\lambda_M/\lambda_2}+1)\) is the convergence rate that exclude the transition eigenvalue and \(\lambda_1 \leq \cdots \leq \lambda_M\) are the eigenvalues of \(\varepsilon^\dagger\). Because the number of iterations scales logarithmically with the smallest eigenvalue of \(\varepsilon^\dagger\), the convergence issues remain highly localized around the transition, even though the estimation based on GMRES appears to be optimistic. Figure 5 shows an example of a convergence plot for BCC iron at the transition, on which we see the characteristic plateau caused by the small eigenvalue of the dielectric matrix.

4 Preconditioning SCF iterations with a diagonal approximation of \(\chi_0\)↩︎

Preconditioning the SCF iterations in KS-DFT consists in replacing the residuals \(\rho_n - F(V_\mathrm{HXC}(\rho_n))\) with the preconditioned residuals \[\begin{align} P^{-1}(\rho_n - F(V_\mathrm{HXC}(\rho_n))) \end{align}\] in the SCF iterations. The matrix \(P\) is referred to as the preconditioner. The jacobian associated with this residual is now the preconditioned adjoint dielectric matrix \(P^{-1}\varepsilon^\dagger\). Thus, to improve the conditioning of the SCF iterations, \(P\) should approximate the adjoint dielectric matrix, \(\varepsilon^\dagger = 1- \chi_0 K_\mathrm{HXC}\).

For semi-local exchange and correlation functionals, the \(K_\mathrm{HXC}\) operator is relatively computationally inexpensive, whereas the operator \(\chi_0\) is too expensive to use exactly in practice. In [23], the authors propose to approximate directly the non-interacting susceptibility \(\chi_0\) instead of the inverse dielectric matrix. This shifts the focus of developing effective preconditioners to finding a computationally efficient yet relevant approximation for \(\chi_0\). The well-known Kerker preconditioner [16] falls into this as it approximates \(\chi_0\) by a constant (cf [23]), as does the LDOS preconditioner, proposed in the aforementioned article, that uses an approximation based of the local density of states.

When deriving its expression using perturbation theory [23], the non interacting susceptibility can be decomposed into three contributions, \[\begin{align} \chi_0 = \chi_0^{\delta \epsilon} + \chi_0^{\delta \epsilon_F} + \chi_0^{\delta \psi} , \end{align}\] arising respectively from the variations of the eigenvalues, the Fermi level, and the Kohn-Sham orbitals: \[\begin{align} &\left( \chi_0^{\delta \epsilon} \right)^{\alpha \beta} = \delta_{\alpha \beta} \sum_{i=1}^\infty f'(\epsilon_i-\epsilon_F) \ket{\rho_{ii}^\alpha} \bra{\rho_{ii}^\alpha} \\ &\left( \chi_0^{\delta \epsilon_F} \right)^{\alpha \beta} = \frac{1}{D} \ket{D_\mathrm{loc}^\alpha} \bra{D_\mathrm{loc}^\beta}\\ &\left( \chi_0^{\delta \psi} \right)^{\alpha \beta} = \delta_{\alpha, \beta} \sum_{\substack{i, j=1 \\i \neq j}}^\infty \frac{f(\epsilon_i^\alpha - \epsilon_F) - f(\epsilon_j^\alpha - \epsilon_F)}{\epsilon_i^\alpha - \epsilon_j^\alpha} \ket{\rho_{ij}^\alpha} \bra{\rho_{ij}^\alpha} \end{align}\] for \(\alpha\), \(\beta\) in \(\{+, -\}\). In these expressions, we use the orbital codensities \(\rho_{ij}(\mathbf{r}) = \psi_i(\mathbf{r})\overline{\psi_j}(\mathbf{r})\), the density of states \(D = \sum_\alpha \sum_i f'(\epsilon_i^\alpha - \epsilon_F)\), and the local density of states \(D_\mathrm{loc}^\alpha(\mathbf{r}) = \sum_i f'(\epsilon_i^\alpha - \epsilon_F) \rho_{ii}^\alpha(\mathbf{r})\).

Among these three terms, the most computationally intensive contribution when applying \(\chi_0\) to a vector is by far \(\chi_0^{\delta \psi}\).

In this section, we suggest a preconditioner aimed at addressing convergence issues associated with ferromagnetism. Drawing inspiration from the Stoner model in which only the eigenvalues are affected by magnetization, we propose to neglect the computationally demanding contribution \(\chi_0^{\delta \psi}\). Therefore, we propose the following model for \(\chi_0\): \[\begin{align} \chi_0^\mathrm{diag}= \chi_0^{\delta \epsilon} + \chi_0^{\delta \epsilon_F} \end{align}\] that keeps only the diagonal contributions of Kohn-Sham orbitals interactions.

This approximation generalizes the Stoner model. Consider a non-magnetic KS-DFT reference state \(\rho^\pm = \sum_{i=1}^\infty \\ f(\epsilon_i - \varepsilon_F) \rho_{ii}\) and the associated model non-interacting susceptibility \(\chi_0^\mathrm{diag}\). We study the following linear fixed-point equation for a variation of the density \(\delta \rho\) with respect to the reference state: \[\begin{align} \begin{aligned} &\delta V = K_\mathrm{Stoner} \delta \rho \\ &\delta \rho = \chi_0^\mathrm{diag}\delta V , \end{aligned} \end{align}\] where the potential to density map is given by \(\chi_0^\mathrm{diag}\) and the density to potential map is given by the Stoner model, with \(\left( K_\mathrm{Stoner} \delta \rho \right)^\alpha = -\frac{\alpha}{2} I M\) and the magnetization \(M = \int_\Omega \delta \rho^+ - \delta \rho^-\). This leads to: \[\begin{align} \delta \rho^+ - \delta \rho^- = I M \left( -\sum_{i=1}^1 f'(\epsilon_i - \epsilon_F) \rho_{ii} \right) . \end{align}\] Integrating the above equation over the unit cell yields \[\begin{align} M(1-ID^0) = 0 , \end{align}\] which recovers the first-order term of the Stoner equation \(M=G(M)\), where a phase transition is only possible at the Stoner criterion, \(1 - ID^0=0\).

4.1 \(\chi_0^{\mathrm{diag}}\) on supercells↩︎

To evaluate the range of applicability of the proposed model, \(\chi_0^\mathrm{diag}\), we examine its behavior in supercell.

Consider a supercell \(\Omega_\mathcal{N}\) with periodic boundary conditions, composed of \(\mathcal{N}\) repetitions of the unit cell \(\Omega\), in which the external potential \(V_\mathrm{ext}\) is \(\Omega\)-periodic. The Bloch theorem ensures that both the orbital densities \(\rho_{ii}^\alpha\) and the local densities of states \(D_\mathrm{loc}^\alpha\) are also \(\Omega\)-periodic. Hence, if we consider a potential perturbation \(\delta V\) over the supercell, without \(\Omega\)-periodic component, the scalar products \(\braket{\rho^\alpha_{ii}}{\delta V^\alpha}\) and \(\braket{D_\mathrm{loc}^\beta}{\delta V^\alpha}\) vanishes and, \[\begin{align} \chi^\mathrm{diag}_0 \delta V = 0 . \end{align}\]

Put differently, the model \(\chi_0^\mathrm{diag}\) acts exclusively in \(\Omega\), and will fail to capture any effect taking place across unit cells.

In particular, this diagonal model will not be able to treat charge sloshing issues, as they are caused by long-range modes that extend beyond the unit cell. Additionally, it will fail to treat convergence issues caused by small antiferromagnetic modes. Indeed, such modes are characterized by spin densities whose sign alternate between neighboring atoms, and their periodicity exceeds that of \(\Omega\).

4.2 Hybrid preconditioner \(P^\mathrm{hybrid}\)↩︎

Convergence issues caused by the long-range divergence of the Coulombic interaction happen systematically in metallic systems with large simulation cells. A preconditioner that fails to address these convergence issues is therefore of limited utility.

In order to build a robust preconditioner treating both convergence issues caused by charge sloshing and by small ferromagnetic modes, we propose to combine the model \(\chi_0^\mathrm{diag}\) with the LDOS preconditioner developed in [23]. The LDOS preconditioner is specifically designed to address charge sloshing problems. Within the collinear spin framework, the non-interacting susceptibility model \(\chi_0^\mathrm{LDOS}\) of the LDOS preconditioner is given by the integral kernels \[\begin{align} \begin{aligned} \left(\chi_0^\mathrm{LDOS}\right)^{\alpha \alpha}(\mathbf{r}, \mathbf{r}') = &-D_\mathrm{loc}^\alpha(\mathbf{r})\delta(\mathbf{r}-\mathbf{r}') \\ &+ \frac{1}{D} D_\mathrm{loc}^\alpha(\mathbf{r}) D_\mathrm{loc}^\alpha(\mathbf{r}')\\ \end{aligned} \\ \begin{align} \left(\chi_0^\mathrm{LDOS}\right)^{\alpha \beta}(\mathbf{r}, \mathbf{r}') = \frac{1}{D} D_\mathrm{loc}^\alpha(\mathbf{r}) D_\mathrm{loc}^\beta(\mathbf{r}') \quad \text{when } \alpha \neq \beta \end{align} \end{align}\] for \(\alpha, \beta \in \{+, -\}\). This model uses the local density of states (LDOS) \(D_\mathrm{loc}^\alpha\), which localizes the preconditioner, making it suitable for inhomogeneous systems. The spin channel representation shows that \(\chi_0^\mathrm{LDOS}\) is also “local” in spin as it uses the corresponding local density of states for each spin channel. This allows the LDOS preconditioner to capture different behaviors in each spin channel. For instance, it is well adapted for half-metals, that are metallic in one spin channel and insulating in the other.

Since the two models \(\chi_0^\mathrm{LDOS}\) and \(\chi_0^\mathrm{diag}\) are designed to capture different modes, respectively, stiff long-range modes arising from the Coulombic interaction and soft short-range ferromagnetic modes arising from the exchange and correlation, it is reasonable to use the following hybrid preconditioner: \[\begin{align} P^\mathrm{hybrid}= I - \chi_0^\mathrm{LDOS}K_H - \chi_0^\mathrm{diag}K_\mathrm{XC}. \end{align}\] This is the final expression of the proposed preconditioner, that we will refer to as the hybrid preconditioner in the rest of the paper. Its performances are tested in section 5.

4.3 Implementation details and computation cost↩︎

In SCF iterations, we need to apply the inverse of the preconditioner \(P^\mathrm{hybrid}\). This inversion is done iteratively with the GMRES [37] algorithm.

The \(\chi_0^\mathrm{diag}= \chi_0^{\delta\epsilon} + \chi_0^{\delta\epsilon_F}\) operator is applied with the following formula for \(\chi_0^{\delta\epsilon}\): \[\begin{align} \chi_0^{\delta\epsilon} \delta V^\alpha = \sum_{i=1}^\infty f'(\epsilon_{i}^\alpha-\epsilon_F) \braket{\rho_{ii}^\alpha}{\delta V^\alpha} \rho_{ii}^\alpha , \end{align}\] where the quantity \(f'(\epsilon_i^\alpha - \epsilon_F)\) is non-negligible only for \(\epsilon_i^\alpha\) close to the Fermi level. If \(p\) represents the proportion of computed bands sufficiently close to the Fermi level, computing the application of \(\chi_0^\mathrm{diag}\) to a vector has a cost similar to \(p\) times the cost of computing the density from the wavefunctions \(\psi_i\).

In plane wave KS-DFT codes, the computational cost of one self-consistent iteration is dominated by the diagonalization step, which scales in \(\mathcal{O}(N^3)\), where \(N\) is the number of electrons. The application of \(\chi_0^\mathrm{diag}\), on the other hand, is dominated by the FFTs and scales in \(\mathcal{O}(N^2 \log(N))\), as computing the orbital densities \(\rho_{ii}\) involves transforming \(\psi_i\) from from the reciprocal space to the direct space. The preconditioner cost becomes negligible for large systems but remain significant for smaller ones.

Since \(\chi_0^\mathrm{diag}\) must be applied at each GMRES iteration, there is a tradeoff to find between the computational cost of applying the same FFT several times and the memory cost of saving these FFTs throughout the GMRES iterations. If the FFT results are saved, the computational cost does not significantly increase compared to the use of the LDOS preconditioner alone. If the FFT results are not saved, the preconditioner will have a cost similar to \(p \cdot n_\mathrm{GMRES}\) density calculations, and we have found that limiting the number of GMRES iterations \(n_\mathrm{GMRES}\) to 20 yields satisfactory results.

We have also noticed that, like the LDOS preconditioner but to a greater extend, the hybrid preconditioner is very sensitive to the convergence in k-points. Indeed, the operator \(\chi_0\) is more challenging to converge in k-points than other KS-DFT quantities. To overcome this difficulty, we propose to artificially increase the smearing temperature in the preconditioner to ensure its smoothness. This introduces one parameter in our otherwise parameter-free preconditioner, but it appears that we can largely increase the temperature before degrading the convergence, so this parameter doesn’t need to be adjusted precisely. In practice, we used a preconditioner smearing temperature of 0.01 Hartree for all our tests.

5 Numerical results↩︎

In this section, we present the results obtained with the proposed preconditioner \(P^\mathrm{hybrid}\) on different magnetic systems. We start by analyzing in detail the behavior of the hybrid preconditioner on BCC iron, whose phase transition we studied in section 3. Then we present numerical results on a wider variety of systems and compare the performance of the hybrid preconditioner with the widely used Kerker preconditioner [16] and the LDOS preconditioner [23], of which \(P^\mathrm{hybrid}\) is an extension.

5.1 Preconditioning BCC iron with \(P^\mathrm{hybrid}\)↩︎

For the hybrid preconditioner to efficiently address convergence issues related to small magnetic modes, it must accurately reproduce the small modes of the dielectric matrix. In Figure 6, we represent the smallest eigenvalue (\(\lambda_1\)) of the preconditioner \(P^\mathrm{hybrid}\) and of the adjoint dielectric matrix. The smallest eigenvalue of \(P^\mathrm{hybrid}\) follows closely that of the adjoint dielectric matrix, with a constant discrepancy of approximately 0.1. We also compute the angle, \(\angle(x, y) = \arccos \left( \langle x, y \rangle / (\|x\| \|y\|) \right)\), between the two eigenvectors associated with the smallest eigenvalue of \(\varepsilon^\dagger\) and \(P^\mathrm{hybrid}\). We observe that this angle remains around 10° near the first magnetic phase transition.

This rather good approximation of the smallest eigenmode of the adjoint dielectric matrix by the preconditioner \(P^\mathrm{hybrid}\) translates to a reduction of the spectral condition number of the preconditioned adjoint dielectric matrix \((P^\mathrm{hybrid})^{-1} \varepsilon^\dagger\). Figure 7 shows the spectral condition number of \(\varepsilon^\dagger\) and of \(\left( P^\mathrm{hybrid}\right)^{-1} \varepsilon^\dagger\). We can see that the hybrid preconditioner significantly reduces the spectral condition number of the preconditioned dielectric matrix near the magnetic transition, where the unpreconditioned dielectric matrix is ill-conditioned. Away from the transition, the behavior of the hybrid preconditioner is governed by the LDOS preconditioner. In this specific example, with no charge sloshing issues, the LDOS preconditioner happens to worsen the condition number of the dielectric matrix.

This reduced spectral condition number around the transition translates into a significantly smaller number of SCF iterations at the transition, as illustrated in Figure 8. Figure 5 shows an example of a convergence plot at the transition, on which we see that the hybrid preconditioner eliminates the plateau. Thus, the hybrid preconditioner seems appropriate to address convergence issues in this stereotypical ferromagnetic system that is cubic iron.

5.2 \(P^\mathrm{hybrid}\) on a wider range of magnetic systems↩︎

We now present numerical results on a wider variety of magnetic systems. The DFT simulations presented here were performed with the Abinit software [48], with a plane-wave basis set and norm-conserving pseudopotentials [44]. We used the same cutoff energy of 50 Hartree for all systems and a Monkhorst-Pack k-point grid. The exchange-correlation functional used is the GGA PBE functional [45], combined with a Fermi-Dirac occupation scheme. For all systems, convergence with respect to the cutoff energy and k-point grid ensures a maximum error of \(10^{-3}\) Ha in the total energy and a maximum relative error of \(2\%\) in the magnetization.

We consider a range of ferromagnetic systems to compare the efficiency of the Kerker preconditioner, the LDOS preconditioner, and the hybrid preconditioner. We chose systems with a small total magnetization per unit cell, and we adjusted the lattice constants of these systems to bring them closer to a magnetic phase transition, where small magnetic modes are expected to appear in the dielectric matrix. All SCF procedures were accelerated using Pulay mixing with a history size of 10 and a damping parameter of 0.8. The Kohn-Sham eigenproblem was solved iteratively with a degree 16 Chebyshev polynomial filter.

When the system parameters are intentionally adjusted to place the system near a phase transition, we observe significant accelerations with the hybrid preconditioner. Figure 10 shows the convergence plots for four systems that are efficiently accelerated with the use of \(P^\mathrm{hybrid}\), BCC ¹, BCC ² , ³ and ⁴, whose structures are taken from the Materials Project [49]. For these systems, we observe a plateau in the convergence curves obtained with the Kerker and LDOS preconditioners, characteristic of the presence of small eigenvalues in the dielectric matrix. The hybrid preconditioner either eliminates this plateau or significantly reduces its size. Still, we observe that \(P^\mathrm{hybrid}\) doesn’t solve all convergence issues caused by ferromagnetic phase transitions. For instance, when the BCC Cobalt system is brought even closer to its transition, it doesn’t converge with any of the preconditioners, including the hybrid scheme.

In systems simulated away from their magnetic phase transition, which present no convergence issues related to magnetism, the hybrid preconditioner generally has no significant effect and doesn’t worsen the convergence of the SCF iterations. However, the hybrid preconditioner is more sensitive to the convergence in k-points. In our experiment, we have found a few systems, with insufficient k-point sampling, for which the use of \(P^\mathrm{hybrid}\) significantly worsens the convergence, even with the increased smearing temperature in the preconditioner mentioned in section 4.3.

The newly proposed hybrid preconditioning scheme thus appears promising to address convergence issues related to magnetism. We successfully accelerated calculations for the three typical ferromagnetic elements, iron, cobalt, and nickel, as well as several alloys. Still the applicability of this preconditioning scheme remains limited to ferromagnetic phase transitions. That said, it seems that in most cases, when \(P^\mathrm{hybrid}\) is used outside its range of applicability, it has little effect on convergence, and importantly, it doesn’t worsen it.

6 Conclusion↩︎

We have investigated convergence issues in magnetic systems and demonstrated that, although these convergence issues are localized around phase transitions, they can be significant. Inspired by the Stoner model, we derived a model that targets problematic ferromagnetic modes, combining it with the LDOS preconditioner to address both magnetic and charge sloshing issues.

Numerical results have shown overall good results in the small applicability range of the proposed hybrid preconditioner, where convergence issues caused by magnetism become significant. Outside its small range of applicability, the hybrid preconditioner does not worsen the convergence.

Overall, this study provides a deeper understanding of convergence issues in spin-KS-DFT and offers a practical solution for improving the robustness of SCF iterations in challenging systems.

Another issue commonly encountered in magnetic systems is the possible existence of multiple stable solutions to the Kohn–Sham equations. Although preconditioning can influence which solution the SCF iterations converge to, the problem of identifying the physically relevant solution is largely independent of accelerating convergence toward a given solution. A natural extension to this work would to explore the solution landscapes in magnetic systems.

Acknowledgments↩︎

Part of the work was performed using HPC resources from CCRT (CEA, France) on the Topaze supercomputer.

References↩︎

---

1. [mp-102] with a \(49\times49\times49\) k-point grid, a smearing temperature of 0.001 Ha and the lattice constants scaled by a factor 0.91 (convergence issues from 0.908 to 0.915).↩︎

2. [mp-23] with a \(49\times49\times49\) k-point grid, a smearing temperature of 0.005 Ha and the lattice constants scaled by a factor 0.830104 (convergence issues from 0.830096 to 0.830106).↩︎

3. [mp-1401] with a \(25\times25\times25\) k-point grid, a smearing temperature of 0.005 Ha and the lattice constants scaled by a factor 0.9332 (convergence issues from 0.9312 to 0.9336).↩︎

4. [mp-570428] with a \(25\times25\times25\) k-point grid, a smearing temperature of 0.005 Ha and the lattice constants scaled by a factor 0.96 (convergence issues from 0.94 to 0.97).↩︎