An adaptive preconditioning scheme for the self-consistent field iteration and generalized stacking-fault energy calculations

Vitek defined the SFE as Vitek (1968):

Where $\Delta E$ is the energy difference between the original structure and the stacking-fault structure, $s$ is the area of the slip plane.

We first focus on the modeling methods of the single-shift and triple-shift supercell models for GSFE curve calculations. In Al, Cu, and Si structures, the stacking sequence of the close-packed atomic layers along the z-$\langle 111\rangle$ direction can be denoted by ABCABC. Here, we define a 6-layer stacking configuration as the supercell unit block to mitigate the interaction impact between two stacking-faults along the z-axis direction. The $\langle{1}\bar{1}0\rangle$ and $\langle{1}0\bar{1}\rangle$ directions on the (111) plane are selected as the basis vectors to make the supercell unit block as small as possible to reduce computational cost.

The single-shift supercell model consists of two interacting supercell unit blocks, where one remains stationary, and the other moves along the $\langle\bar{1}\bar{1}2\rangle$ direction. However, due to periodic boundary conditions, a relative displacement in the $\langle\bar{1}\bar{1}2\rangle$ direction between the blocks is induced. We add a 10 Å vacuum layer in the z-axis direction to eliminate the interaction between the blocks and their periodic images. The triple-shift supercell model consists of three interacting supercell blocks. They simultaneously slips in the $\langle{1}\bar{1}0\rangle$, $\langle 01\bar{1}\rangle$ and $\langle\bar{1}01\rangle$ directions, respectively. As a result, the displacements between adjacent blocks are $\langle 2\bar{1}\bar{1}\rangle$, $\langle\bar{1}2\bar{1}\rangle$, and $\langle\bar{1}\bar{1}2\rangle$. The sum of these three displacement vectors is zero, no need to add a vacuum layer and preserving the periodic boundary conditions along the z-axis direction. The modeling methods of the two models are shown in the Fig. 1. We took 18 equidistant points to approximate the unit (111) plane $\langle\bar{1}\bar{1}2\rangle$ direction GSFE curve of Al, Cu, and Si (the (111) plane of Si includes the glide plane and the shuffle plane, and we only compute the glide plane).

We utilized ELK-7.2.42 packages. It is based on the FLAPW method, which divides the space of the crystal into two regions: muffin-tin spheres and interstitial regions. The wavefunctions are expanded using a combination of atomic orbitals in the muffin-tin region and planewave in the interstitial region. The wave functions $\psi^{\sigma}_{\nu\textbf{k}}$ for spin-$\sigma$ valence orbital of band index $\nu$ with Bloch vector k is:

Where $c_{\nu\textbf{k}}^{\textbf{G}\sigma}$ is the expansion coefficients of the wavefunction, G are all reciprocal lattice vectors dual to the lattice vectors defining the periodic domain up to the $\textbf{G}_{\rm{max}}$, and the basis functions $\varphi^{\sigma}_{\textbf{{G}}}(\textbf{k},\textbf{r})$ can be written as:

Where $L$ abbreviates the quantum numbers $l$ and $m$, $\alpha$ is the muffin-tin sphere, $Y_{L}(\hat{\bf{r}})$ is the spherical harmonics, $A^{\sigma\alpha\textbf{G}}_{jL}(\textbf{{k}})$ is the matching coefficient, $j$ is the order of the basis function, for $j=0$ and 1, $u^{\sigma\alpha}_{jl}(r_{\alpha})$ is the solution of the radial Schrödinger equation at a linearization energy and its energy derivative evaluated at this same energy. In the interstitial region, the function is the planewave. As the basis sets in different regions of the crystal are different, it results in different representations of the density and potential in the two regions:

Based on the above, an eigenvalue problem can be derived:

where Hamiltonian $\hat{H}$ is a sum of corresponding terms:

${T}_{0}$ is kinetic energy, ${V}_{\rm ext}$ is external-potential, ${V}_{\rm H}$ is Hartree potential, ${V}_{\rm xc}$ is exchange-correlation potential. Of which, ${V}_{\rm xc}$ and ${V}_{\rm H}$ are the local potentials and explicitly density dependent, and the charge density depends on the solution of the Eq.(6), leading to the nonlinearity of the Kohn-Sham equation. One need to resort to a fixed-point iteration scheme:

The residual $\textbf{f}(\textbf{n}_{m})$ is defined as: $\textbf{f}(\textbf{n}_{m})=\textbf{Q}({\textbf{n}_{m}})-\textbf{n}_{m}.$ Based on this idea, the simple mixing scheme for the vector $\textbf{n}_{m}$ (represents charge density $\rho_{m}$ or potential ${v}_{m}$, in this work, we employ potential mixing) can be derived:

where $P$ is the preconditioner. Through SCF iterations and continuously updating the charge density and potential until the desired convergence accuracy is achieved, the ground state-energy can be obtained. From this, one can further derive various physical properties.

We employed the FLAPW combined with the local orbital basis set in ELK, the muffin-tin radii were set to $R_{\rm MT}^{\rm Al}=2$ bohr, $R_{\rm MT}^{\rm Cu}=2.2$ bohr, and $R_{\rm MT}^{\rm Si}=2$ bohr for Al, Cu, and Si, respectively. The cut-off radii for these systems were all set to $R^{\rm MT}_{\rm mix}K_{\rm max}=10$. The Brillouin zone sampling was set as $11\times 11\times 1$ for k-points. We adopt the generalized gradient approximation in the Perdew–Burke–Ernzerhof parametrization for the exchange-correlation functional Perdew et al. (1998). We employed the original Pulay mixing scheme Johnson (1988), and set the convergence criterion as energy between two consecutive iteration steps being smaller than $10^{-6}$ Hartree. All atoms relax along the z-axis direction during the calculation, and atomic forces were converged to less than $0.001$ Hatree/bohr.

To calculate the GSFE curve in the $\langle\bar{1}\bar{1}2\rangle$ direction, we took 18 equidistant points along the unit $\langle\bar{1}\bar{1}2\rangle$ direction, and adopted the 1/18$\langle\bar{1}\bar{1}2\rangle$ point as an example to show the convergence of the single-shift and triple-shift supercell models calculations in Fig. 2. It is shown that only the single-shift supercell model of Al converges in 50 steps; others do not converge in 50 steps. Obviously, calculating the GSFE curve under current scheme is inefficient.

Firstly, we define the Jacobian as:

The residual propagation of Eq.(9) is given by:

The convergence condition can be written as:

If there are small eigenvalues in $(JP)^{-1}$, the conditions of Eq.(12) are difficult to satisfy, which implies the long-wavelength divergence behavior of the Jacobian. But it is hard to calculate the eigenvalues of $(JP)^{-1}$. Practically, we can use the subspace of Pulay scheme to approximate $(JP)^{-1}$.

In the Broyden’s second method, a sequence of low-rank modifications are made to modify the initial guess of the inverse Jacobian, the recursive formula can be derived from the following constrained optimization problem Marks and Luke (2008); Fang and Saad (2009):

where $H_{m-1}$ is the approximation to the inverse Jacobian in the $(m-1)$th Broyden update. $N_{m-1}$ and $F_{m-1}$ are:

where $\ell$ is the size of the subspace. Eq.(13) has an analytical solution as Zhou et al. (2018):

And we can replace the $H_{m-1}$ in Eq.(16) to the initial guess $H_{1}$ as the inverse of Jacobian:

Then we follow the quasi Newton approach to solve the fixed-point iteration, this leads to the Pulay mixing scheme, which is also the method used in Elk codes:

The $H_{m}$ of Eq.(17) satisfies the constraint condition:

and for the Jacobian $J$, it also satisfies:

From Eq.(19) and Eq.(20), we can calculate the eigenvalues of $P^{-1}H_{m}$ replacing the eigenvalues of $(JP)^{-1}$, because $H_{m}$ provides an effective approximation to $J$ within the subspace defined by $F_{m}$. The eigenvalues of $P^{-1}H_{m}$ can be obtained by solving a generalized eigenvalue problem:

Eq.(21) can be shift as:

We only need to calculate the eigenvalues of the left-hand side of the Eq.(22), select the smallest eigenvalue as a posterior indicator. When the indicator is too small, it signifies the long-wavelength divergence behavior of the Jacobian. Note that solving Eq.(22) requires a little computational overhead because $(H_{m}-P)F_{m-1}$ has been calculated in Pulay’s update:

and the size of the subspace $\ell$ is generally smaller than 50. The details of relevant derivation can be found in Ref Zhou et al. (2018).

In the general scheme (also in our tests), $H_{1}$ of Eq.(17) is set to $\alpha I$, where $\alpha$ is a scalar parameter usually set to 0.4. The indicators of Al, Cu, and Si systems for the single-shift and triple-shift supercell models tests are shown in Fig. 3. We found that they all have small indicators during the SCF iterations. It means that the convergence condition of Eq.(12) cannot be satisfied, and the residual may be magnified in the corresponding step, leading to the SCF iteration is hard to converge.

When we find the long-wavelength divergence behavior of the Jacobian through a posterior indicator, an effective method to suppress the long-wavelength divergence behavior of the Jacobian is the Kerker preconditioner. The Kerker preconditioner is based on the Thomas-Fermi screening model. The dielectric function of the Thomas-Fermi screening model in reciprocal space can be expressed as:

where $k_{\rm TF}$ is the Thomas-Fermi vector Zhou et al. (2018). The Kerker preconditioner can be derived by inverting the dielectric matrix Kresse and Furthmüller (1996b); Lin and Yang (2013), and taking the Kerker preconditioner into simple mixing scheme can be expressed as:

In the FLAPW method, due to the violent oscillation of the charge density near the nucleus, the Fourier transform is hard, so the Kerker preconditioner must be implemented in real space Wimmer et al. (1981):

We replace the $k_{\rm TF}$ by $\lambda$, because Thomas-Fermi wave number of the homogeneous electron gas and the electronic properties of the actual system are hard to equal Winkelmann et al. (2020). Then, the key of the problem is the treatment of $(\nabla^{2}-\lambda^{2})^{-1}$, which can be related to the solution of screened Poisson’s equation Lin and Yang (2013); Hinzen et al. (2020):

For the potential mixing, where ${\rho}(\textbf{r})=-{v}_{m}(\textbf{r})/4\pi$. In order to solving the screened Poisson’s equation in the FLAPW basis, we need to divide the solution into the muffin-tin and the interstitial regions, and firstly solve the screened Poisson’s equation for the interstitial regions in reciprocal space, the solution of interstitial regions provides the boundary conditions for muffin-tin regions. Then, the solution of the muffin-tin regions can be solved by Green’s function method. This is based on the pseudo-charge method proposed by Weinert Weinert (1981).

For the solution of the interstitial region, we constructed the pseudo-charge $\bar{{\rho}}(\textbf{r})$ that can be expanded in terms of the finite set of G-vectors. The charge density in the interstitial region is extended to the entire space, and the charge density in the muffin-tin region is the original charge density subtracted from the extended interstitial charge density. We now replace the real muffin-tin charge density with a smooth density $\tilde{{\rho}}^{\alpha}(\textbf{r})$:

Concurrently, the pseudo-charge must preserve equivalent multipole moments (the same as the real charge distribution) within the muffin-tin region. This implies that the multipole moments of smooth density $\tilde{q}^{\alpha}_{lm}$ are determined by the difference between the muffin-tin and interstitial multipole moments:

Where the multipole moments for the muffin-tin $q^{\alpha,\rm{MT}}_{lm}$ and interstitial regions $q_{lm}^{\rm I}$ can be expressed as:

$i_{l}(r)$ and $k_{l}(r)$ are the modified spherical Bessel function of the first and second kind. Then the pseudo-charge density in reciprocal space can be written as:

The pseudo-charge has a more rapidly converging Fourier expansion than the original charge density, and the $V_{I}({\bf G})$ in the interstitial can be get:

The potential of the interstitial region provides the value at the boundary of the muffin-tin sphere. The solution of potential within the muffin-tin sphere is a Dirichlet boundary-value problem, which can be solved by the Green’s function method Jackson (1999):

One characteristic of Kerker preconditioner is that its convergence is not affected by an increase in the number of atoms on the long-axis for metal systems. Based on the implementation of the Kerker preconditioner, we tested the convergence for the triple-shift supercell model of Al, Cu with different thicknesses on the z-axis. We set $\lambda_{\rm Al}=0.6$, $\lambda_{\rm Cu}=0.8$ Winkelmann et al. (2020) , the convergence results are shown in Fig. 4. We found that the convergence speed is significantly improved by using the Kerker preconditioner. At the same time, we found that the step of convergence almost unchanged for different atomic layer thicknesses. Obviously, the Kerker preconditioner successfully suppressed the long-wavelength divergence behavior of the Jacobian.

In the SCF iteration, we establish an initiation criterion for the Kerker preconditioner by setting a threshold of less than 0.2 for a posterior indicator. This leads to an adaptive preconditioning algorithm. The algorithm process is shown in Algorithm. 1.

At the beginning of the program, we set $H_{1}$ to $\alpha I$. A posterior indicator is calculated in every SCF iteration. When the indicator is less than 0.2, break the current step, and the current subspace will also be discarded. Subsequently, a new SCF iteration will be started with the Kerker preconditioner initiating, continue iteration until convergence accuracy is reached.

We tested the convergence and posterior indicators of Al, Cu, and Si using the adaptive preconditioning scheme in the two supercell models. We set $\lambda_{\rm Al}=0.6$, $\lambda_{\rm Cu}=0.8$ and $\lambda_{\rm Si}=0.6$. As shown in Fig. 5, small indicator appears around the 5 steps. Then the Kerker preconditioner is initiated, the indicators are approximately 1 and remains stable. It leads to a significantly improved convergence speed compared with the original algorithm. This scheme allows the system to adjust the choice of preconditioning method more reasonably without prior information.

It is worth noting that the Kerker preconditioner is aimed at metallic systems, but in out tests, it is also effective for Si, which is a semiconductor. We calculated the density of states (DOS) for the original structure, the single-shift supercell model, and the triple-shift supercell model of Si, as shown in Fig. 6. Although the original structure of Si is insulated, in the single-shift and triple-shift supercell models, planar defects in stacking-fault structures generate free electrons, giving the system metallic properties. Moreover, in the single-shift supercell model, the planar defects introduced by a vacuum layer enhanced the metallicity. Our programs identified the long-wavelength divergence behavior of the Jacobian and initiated Kerker preconditioner improving convergence in these models.