Robust mixing in self-consistent linearized augmented planewave calculations

We devise a mixing algorithm for full-potential (FP) all-electron calculations in the linearized augmented planewave (LAPW) method. Pulay's direct inversion in the iterative subspace is complemented with the Kerker preconditioner and further improvements to achieve smooth convergence, avoiding charge sloshing and noise in the exchange-correlation potential. As the Kerker preconditioner was originally designed for the planewave basis, we have adapted it to the FP-LAPW method and implemented in the exciting code. Applications to the $2\times 2$ Au(111) surface with a vacancy and to the Pd(111) surface demonstrate that this approach and our implementation work reliably with both density and potential mixing.


Introduction
Density-functional theory (DFT) [1,2] is used for a wide range of problems that cover electronic, mechanical, and vibrational properties of atoms, molecules, and solids. The diversity of applications to various materials requires that numerical algorithms employed in electronic-structure codes are robust. Otherwise, a poorly chosen or implemented algorithm may lead to unnecessary long computation times as well as to non-converging calculations [3]. Such issues hinder not only individual calculations but also are particularly harmful high-throughput studies.
One particular problem commonly recurring in electronic-structure theory is convergence of the self-consistent Kohn-Sham equation: where ik , v KS (r), and ψ ik (r) represent Kohn-Sham eigenenergies, potential, and wavefunctions, respectively. v KS (r) depends on the electron density: where ω k is the weight of the k-point and f ik the occupation factor. Eq. 1 defines a non-linear eigenvalue problem. To solve it, one considers a linear problem with a fixed potential v in KS (r). After computing ψ ik (r) and subsequently calculating ρ(r), one obtains a new potential v out KS (r) that serves, in principle, as v in KS (r) for the next iteration of the KS equation (Eq. 1). This procedure is repeated until self-consistency is reached, namely, v in KS (r) ≈ v out KS (r). In practice, at every step, the potential is constructed as a mixture (linear combination) of v out KS (r) and history of v in KS (r). One can likewise target self-consistency in the electron density instead of the potential, applying a corresponding mixing scheme. Regardless of the choice, the success and efficiency of this procedure depends on how exactly either of these two quantities is updated, and the problem resembles a multivariate optimization. As a result, a number of mixing methods [4,5,6,7,8] relies on standard optimization techniques or is intimately related to them. Despite successful applications to numerous systems, there are challenging cases that require further improvement of these methods. In particular, metallic systems with large unit cells suffer from an instability known as charge sloshing [3,9,10,11,12,13,14,15]. This issue originates from (i) a constant non-zero susceptibility, χ, of metallic systems at small wave-vectors and / or (ii) a factor of G −2 in the Hartree potential at small G components. It can, e.g., be seen in the error of the output potential that indicates that a small error in the input potential for short G is amplified due to these two factors. Several studies have offered a solution to this problem [9, 10, 11, 12, 16, 17, 18]. Among them, the Kerker preconditioner is widely used to reduce the charge-sloshing instability induced at small G (long-wavelength) [10], thus making the self-consistency procedure stable for metallic systems with large unit cells.
Mixing schemes are most frequently developed with planewave / pseudopotential methods in mind. Nevertheless, one can adopt the same techniques for all-electron implementations. For instance, the multisecant Broyden method has been successfully applied in FP calculations with the LAPW basis set [19]. Its implementation has become the default choice in two FP-LAPW codes, i.e. exciting [20] and Wien2k [21]. Our experience shows that this method, especially a multisecant Broyden type-1 method, which estimates a Jacobian [19], performs well in most cases, but is unstable with respect to charge-sloshing. Unfortunately, also the aforementioned solutions to this problem target planewave implementations. Mixing potentials rather than densities, pose additional numerical difficulties. For instance, the generalized-gradient approximation (GGA) for the exchange-correlation functional may lead to noisy potentials in the lowdensity regions, making the problem of finding a self-consistent potential ill-conditioned.
In this paper, we tackle the issues of charge-sloshing and noisy potentials in FP-LAPW calculations. We implement the Kerker scheme and a modified Pulay mixer in the exciting code and discuss their performance with two benchmark systems: (i) the 2 × 2 Au(111) surface with a vacancy and (ii) the Pd(111) surface. We show that the modified Pulay mixer with the Kerker preconditioner is robust for both density and potential mixing.

Kerker mixing
Considering the i-th step of the self-consistent field procedure, the (i + 1)-st guess of the Kohn-Sham potential can be obtained via a simple linear relation: where α is some prefactor. This simple approach may work in many cases, however, it fails in calculations of metallic systems with a large unit cell. The linear-response properties in these systems are such that a small change in v in i triggers a large change in v out i at short wavevectors G as outlined above. Such an instability causes massive charge fluctuations over the self-consistency cycle which is known as charge sloshing. These short-G fluctuations can be suppressed using the Kerker scheme [10]. In a planewave basis, the mixing relation in Eq.
where λ is a parameter that determines the range of G over which changes of the potential are suppressed. Eq. 5 is straightforward to implement in the planewave methods, but not as simple in others. To describe its implementation in the FP-LAPW method, we briefly sketch its characteristics. In this method, the unit cell is divided into the interstitial region (I ) and atomic spheres (commonly known as muffin-tin spheres), that are centered at the atomic positions. The potential (as is the density) is expanded in terms of planewaves in I and in spherical harmonics in the atomic spheres (labelled α), respectivley: Since performing a Fourier transform is not feasible in this case, calculating v in i+1 directly from Eq. 5 is not practical.
To make Eq. 5 applicable in the FP-LAPW case, we transform it to the real space and obtain v in The operatorK = (∇ 2 − λ 2 ) −1 cannot be applied to a function directly. It rather means that V (r) =Kf (r) with some function f (r) is a solution of the screened Poisson equation: If we set f (r) = −4πρ(r), the function V (r) corresponds to the electrostatic potential due to the charge density ρ(r) in a system with the Thomas-Fermi screening. Eq. 8 can be solved in the spirit of the pseudo-charge method suggested by Weinert [23]. It was originally proposed for calculating the Hartree potential in FP-LAPW calculations, but the same idea can be exploited for solving the screened Poisson equation as suggested in Ref. [24]. The original purpose of the algorithm presented in that study was an implementation of the screened exchange in a hybrid exchange-correlation functional but it also meets our needs. The crucial steps of the method are described below. At first, we calculate the screened potential in the interstitial region. To do so, the charge density ρ(r) is replaced by a smooth pseudo-charge densityρ(r), for which it is easy to perform a Fourier transform. Note that the input charge density is given in the same form as the potential in Eq. 6. The pseudo-charge density is constructed as where the first term is the given interstitial density. Note that unlike in Eq. 6 planewaves are allowed to enter the muffin-tin spheres in Eq. 9. The second term is a smooth function that is non-zero only in the muffin-tin spheres.ρ α (r) are chosen such thatρ(r) yields the correct screened potential in the interstitial region. To ensure it, we expand the second term as follows: where Q α lm are constants chosen in such a way thatρ(r) has the same charge multipoles as ρ(r). Since we consider the screened Coulomb interaction, the multipoles should be calculated as where i l (r) is modified spherical Bessel function of the first kind. Finally, the radial functions in Eq. 10 are defined as with n = R α G max /4. Due to this choice of σ α lm (r), the Fourier transform ofρ α (r) can be performed analytically, yielding With the Fourier transforms of both parts in Eq. 9, the potential in the interstitial reads To obtain the potential in the muffin-tin region, we solve the Dirichlet boundaryvalue problem with the original density ρ(r). Employing the Green-function method [25], the potential reads where the Green function and its normal derivative are expressed as and respectively. k l (r) is the modified spherical Bessel function of the second kind, and r < (r > ) is the maximum (minimum) value between r and r . ∂G ∂n is the normal derivative at the atomic sphere boundary.
The method described here is not limited to mixing potentials. The input and output potentials in Eq. 7 can be replaced by respective densities, which does not change the procedure of solving the screened Poisson equation.

Pulay mixing
The Pulay mixer [6,7] also known as direct inversion in the iterative subspace represents an improvement over the linear mixer. It uses a sequence of input potentials v in i (r) and residuals R i (r) = v out i (r) − v in i (r) of previous iterations to estimate the optimal potential and the corresponding residual in the following manner: with the weights ω i subject to the constraint Minimizing the residual R in opt , the weights ω i are determined by [11,12] with According to the conventional Pulay mixing scheme, the next trial potential is given as where α is the mixing parameter. The optimal value of α depends on the system [7]. For example, it is typically bigger in the case of semiconductors and insulators than in the case of metallic systems. Equation 21 contains an integral over the entire unit cell consisting of contributions from the interstitial region and the muffin-tin spheres. When v in i (r) in the interstitial region is noisy (e.g. in case the GGA is employed), this noise propagates into the calculated weights resulting in a numerically problematic potential in Eq. 18. We solve this issue by calculating Eq. 21 as In other words, we ignore the interstitial part of the residuals at this stage and, thus, R in opt is minimized strictly within the atomic spheres only. The summation in Eq. 18 is, however, still performed over entire unit cell. The level of noise in the charge density is typically noticeably lower compared to that in the potential, and, hence, employing Eq. 23 is not useful in case of density mixing. Therefore, we use different equations, i.e. Eqs. 21 and 23, for density and potential mixing, respectively. In both cases, we refer to the method as simple Pulay mixing.
Previous studies [11,12] have shown that introducing the inverse Kerker metric in Eq. 21 helps to prevent charge sloshing. The matrix elements in this approach read Similarly to Eq. 5, this equation is written in the planewave basis and is thus not directly applicable in FP-LAPW calculations. We solve this issue analogously to the case of the Kerker mixer: R j (G)/G 2 corresponds to the bare Coulomb potential due to the charge density ρ(r) = −R j (G)/4π. Finding such a potential is a standard problem in an FP-LAPW code, which is routinely solved using Weinert's method [23]. The expression for the matrix elements can be formally written in real space: The inner integral is evaluated over the entire space. If the range of the outer integral is the whole unit cell, Eq. 25 is equivalent to Eq. 24. We use it for mixing densities. The Kerker transformation can also be used as a preconditioner in the Pulay method [11,12,14,15]. In this case, the Kohn-Sham potential is updated as The expression on the right-hand side exactly matches Eq. 5, and it can be evaluated using the same procedure as described in Sec. 2. We use this preconditioner only in the first n iterations (typically, n = 5). Afterwards, Eq. 22 is employed. This modification of the Pulay method with the inverse Kerker metric and the Kerker preconditionning (termed Pulay-KP) is the most successful method as we will show in the example below. Both methods, the simple Pulay method and the one with the Kerker preconditioner and the inverse Kerker metric, depend on the mixing parameter α that we set to 0.4 for all considered cases. The latter method contains also the screening parameter λ, which strongly influences the convergence properties of the method. Previous studies [31,32] have suggested to use the wavevector of the Thomas-Fermi screening, k T F . Following Ref. [32], we estimate it as where N (ε F ) corresponds to the density of states at the Fermi energy.

Computational details
The first system studied here is the 2 × 2 Au(111) surface with a vacancy. It consists of 5 layers of Au and is termed 5L-(2 × 2)Au(111)-V. The second example is the Pd(111) surface with 15 atomic layers, termed 15L-Pd(111). Both structures are displayed in Fig. 1. The corresponding slabs are constructed based on the bulk geometries with the lattice constants of 4.19Å and 3.95Å for Au and Pd, respectively. These structural parameters are obtained using the PBE exchange-correlation functional. Further structural optimizations of the slabs are not performed. In order to eliminate spurious interactions between the periodic images of the metal slabs, the vacuum spacing in the z direction is set to 30Å and 20Å for both surfaces. The above described methods have been implemented in the all-electron fullpotential code exciting [20], which is then used for all calculations. It employs the linearized augmented planewave plus local orbitals basis-set (LAPW+lo) [26,27,28]. The chosen muffin-tin radii are R Au M T = 2.4 bohr and R P d M T = 1.9 bohr. An LAPW cutoff of R M T G max = 7 is adopted in all considered systems. Exchange and correlation effects are described within the GGA using the Perdew-Berke-Ernzerhof (PBE) parametrization [22]. The Perdew-Wang parametrization [29] of the local-density approximation (LDA) is applied for the sake of comparison. The Brillioun zone (BZ) is sampled on a 4 × 4 × 1 and a 10 × 10 × 1 k-mesh for 5L-(2 × 2)Au(111)-V and 15L-Pd(111), respectively. The self-consistency cycle is considered converged once the totalenergy difference between two consecutive steps is smaller than 10 −6 Ha. Simple Pulay, Pulay with Kerker preconditioner and inverse Kerker metric, and multisecant Broyden methods employ stored potentials (densities) and residuals obtained from m = 12 previous iterations for the 5L-(2 × 2)Au(111)-V system. In 15L-Pd(111) calculations, we have set m = 30, since it is required for a stable performance.

Results and discussion
We start discussing the performance of the implemented mixing methods with the case of 5L-(2 × 2)Au(111)-V. This is a metallic system with a large unit cell and, therefore, it is prone to the charge-sloshing instability [9,10,11,12,13,14,15]. Indeed, 5L-(2 × 2)Au(111)-V is a too complicated problem for the linear mixing and the Kerker method. To illustrate this, we perform calculations applying these methods by mixing the potentials and set the corresponding parameters in Eqs. 4 and 5 to α = 0.4 and λ = 1.0 bohr −1 . To monitor the convergence, we display in Fig. 2 the variation of the total energy between two consecutive steps, ∆E. With the linear mixer, this quantity not only does not converge during the first 100 iterations to a desired threshold, but remains large overall, i.e. within range of 10 4 -10 5 Ha. The Kerker mixing, in turn, shows an improved performance compared to the linear mixing. ∆E decays with the number of iterations, yet this happens too slowly for the Kerker mixer to be practical. Still, we acknowledge that it suppresses charge sloshing.
A further issue that makes 5L-(2 × 2)Au(111)-V challenging is the (typically) noisy GGA potential in the vacuum region. To demonstrate this fact, we compare in Fig. 3 the self-consistent exchange-correlation potentials v xc (r) obtained by LDA and GGA, together with the corresponding electron densities, for 5L-(2 × 2)Au(111)-V. While the densities are essentially indistinguishable on the displayed scale, the differences in the exchange-correlation potentials are immediately apparent. The LDA potential is smooth throughout the unit cell and only slighty jagged in the low-density region (ρ < 10 −4 bohr). The noise obviously stems from the gradient term (∇ρ/ρ 4/3 ), which is prone to rapid variations in low-density regions [30] since a small variation in ρ(r) can cause a large variation in v xc (r). Thus, GGA calculations of low-dimensional systems, i.e. with vacuum, are less stable than those using the electron density only, in particular when mixing potentials. Our implementation of both the simple Pulay and the modified Pulay approaches is insensitive to both of these instabilities. This is demonstrated in Fig. 4 (top panel) where the convergence behavior of the total energy (using the potential mixing) is compared for these two methods and the multisecant Broyden scheme. All three methods converge to the target precision within 25-50 steps, in contrast to the linear mixing and the Kerker method (see Fig. 2). Simple Pulay and Pulay-KP perform better than the multisecant Broyden method, where the noise in the potential is not taken care of.
We also apply these approaches for mixing densities. The convergence behavior of the total energy is shown in the bottom panel of Fig. 4. The simple Pulay and multisecant Broyden methods do not reach self-consistency within 100 steps, while the Pulay-KP method converges in just a few iterations more than in the case of potential mixing. Overall, our calculations for 5L-(2 × 2)Au(111)-V show that the Pulay-KP method is the most efficient method among the considered ones.  Our second benchmark case is 15-layer Pd(111). Besides being metallic and having a large unit cell, this system has a high density of states near the Fermi level, therefore charge-sloshing is even more pronounced [3,33]. The total-energy convergence in calculations using different approaches, i.e. multisecant Broyden, simple Pulay, and Pulay-KP is presented in Fig. 5. Regardless whether the density or the potential is mixed, the only method that succeeds for this system is Pulay-KP. The other two methods do not converge for either kind of mixing. In fact, multisecant Broyden leads in case of density mixing to such an unphysical density that the self-consistency loop cannot be terminated. Pulay-KP, however, reaches the target precision for the total energy in 35 (density) and 36 (potential) steps, respectively. This further indicates that Pulay-PK successfully attenuates the long-wavelength instability which induces charge sloshing.

Summary and Conclusions
In summary, we have reformulated the Kerker preconditioner and the inverse Kerker metric to make them applicable in FP-LAPW calculations. We have implemented the Pulay mixing algorithm modified with these features in the electronic-structure code exciting. Our applications demonstrate that this method is robust and superior to the standard Pulay method and the multisecant Broyden approach.
looked up there.