# The relationship between average eigenvalue and convergence performance in VASP?

I meet some convergence problems when doing self-consistent field calculation in VASP, I find some parameters on VASP wiki that might be able to tune the convergence performance. But I don't quite understand the relationship between convergence and average eigenvalue, as explained below:

In VASP the eigenvalue spectrum of the charge dielectric matrix is calculated and written to the OUTCAR file at each electronic step. This allows a rather easy optimization of the mixing parameters, if required. Search in the OUTCAR file for

eigenvalues of (default mixing * dielectric matrix)

The parameters for the mixing are optimal if the mean eigenvalue Γmean=1, and if the width of the eigenvalue spectrum is minimal. For an initial linear mixing (BMIX≈0) an optimal setting for AMIX can be found easily by setting AMIXoptimal=AMIXcurrent*Γmean. For the Kerker scheme (IMIX=1) either AMIX or BMIX can be optimized, but we recommend to change only BMIX and keep AMIX fixed (you must decrease BMIX if the mean eigenvalue is larger than one, and increase BMIX if the mean eigenvalue Γmean<1). However, the optimal AMIX depends very much on the system, for metals this parameter usually has to be rather small, e.g. AMIX= 0.02.

I checked my result, the average eigenvalue at GAMMA is 0.2025.


eigenvalues of (default mixing * dielectric matrix)
average eigenvalue GAMMA=   0.2025
1.2461  0.6193  0.4021  0.4021  0.2087  0.2087  0.0809  0.0809  0.0725  0.0725
0.0815  0.0815  0.1056  0.1056  0.1091  0.0921  0.0654  0.0654  0.0566  0.0566
0.0386



My question is:
Why is the model's convergence performance better when the average eigenvalue at GAMMA is 1?

In general, the difficulty of an optimisation problem depends on how widely spread the eigenvalues of the Hessian are. To see why this is, consider a minimisation problem where we wish to find the minimum value of the function $$f(\mathbf{r})$$, given an initial "guess" set of inputs, $$\mathbf{r}_0$$.

A common approach is to compute the gradient of the function at $$\mathbf{r}_0$$, i.e. $$\nabla f(\mathbf{r}_0)$$, and use it to determine an improved guess $$\mathbf{r}_1$$. Since the gradient is the direction in which the function increases quickest, and we wish to minimise the function, we use $$-\nabla f(\mathbf{r}_0)$$ as the direction in which to move. We then write $$\mathbf{r}_1 = \mathbf{r}_0 - \alpha \nabla f(\mathbf{r}_0),\tag{1}\label{eq:step}$$ where $$\alpha$$ is the step we take in the search direction $$-\nabla f(\mathbf{r}_0)$$. This step could be a fixed guess, or we could try various different values and try to find the optimum (often called "line minimisation"). (NB In the context of machine learning, $$\alpha$$ is sometimes called the "learning rate".)

### Finding the optimal step length

What determines what $$\alpha$$ should be? Let's Taylor-expand the function around our trial inputs, and write $$f(\mathbf{r}) \approx f(\mathbf{r}_0) + (\mathbf{r}-\mathbf{r}_0)^\dagger\nabla f(\mathbf{r}_0)+\frac{1}{2}(\mathbf{r}-\mathbf{r}_0)^\dagger\mathrm{B}.(\mathbf{r}-\mathbf{r}_0),\tag{2}$$ where $$\mathrm{B}$$ is the Hessian matrix, which is the matrix of second derivatives, $$\mathrm{B}_{ij}=\left.\frac{\partial^2 f}{\partial r_i \partial r_j}\right\vert_{\mathbf{r}=\mathbf{r}_0}.\tag{3}$$ At the minimum, $$\mathbf{r}=\mathbf{r}_\mathrm{opt}$$, we know that the gradient should be zero - so let's differentiate the Taylor expansion to get the gradient expression: $$\nabla f(\mathbf{r}) \approx \nabla f(\mathbf{r}_0)+\mathrm{B}.(\mathbf{r}-\mathbf{r}_0),\tag{4}$$ Now we set this to zero, and rearrange

$$\begin{eqnarray}{} & \nabla f(\mathbf{r}_0)+\mathrm{B}.(\mathbf{r}_\mathrm{opt}-\mathbf{r}_0) &= 0\tag{5}\\ \Rightarrow & \mathrm{B}.(\mathbf{r}_\mathrm{opt}-\mathbf{r}_0) &= -\nabla f(\mathbf{r}_0)\tag{6}\\ \Rightarrow & (\mathbf{r}_\mathrm{opt}-\mathbf{r}_0) &= -\mathrm{B}^{-1}\nabla f(\mathbf{r}_0)\tag{7}\\ \Rightarrow & \mathbf{r}_\mathrm{opt} &= \mathbf{r}_0 -\mathrm{B}^{-1}\nabla f(\mathbf{r}_0)\tag{8} \end{eqnarray}$$ Comparing this with equation \eqref{eq:step}, we see that the basic form is the same, except that the scalar $$\alpha$$ has been replaced by the matrix $$\mathrm{B}^{-1}$$.

If $$\mathrm{B}^{-1}$$ is a diagonal matrix, and all the eigenvalues are the same, then it has the form, $$\tag{9} \mathrm{B}^{-1}=\left(\begin{array}{ccc} \lambda & 0 & 0 & \ldots\\ 0 & \lambda & 0 & \ldots\\ 0 & 0 &\lambda & \\ \vdots & \vdots & & \ddots \end{array}\right) =\lambda\mathrm{I},$$ where $$\mathrm{I}$$ is the identity matrix. In this case, setting $$\alpha=\lambda$$ will give the ideal step length and, in fact, will jump straight to the minimum of the function in a single step.

If the eigenvalues of $$\mathrm{B}^{-1}$$ are different, then there is no single value of $$\alpha$$ which can mimic $$\mathrm{B}^{-1}$$ perfectly. The best approximation is for $$\alpha$$ to match $$\mathrm{B}^{-1}$$ as closely as possible. Since we used $$\mathrm{B}^{-1}$$ because of the relation, $$\mathrm{B}^{-1}\mathrm{B}=I,\tag{10}$$ this leads naturally to the conclusion that we want $$\alpha\mathrm{B} \approx I.\tag{11}$$ The identity matrix $$\mathrm{I}$$ has eigenvalues of 1, and the best approximation for $$\alpha$$ is when the average eigenvalue of $$\alpha\mathrm{B}$$ is one.

### Density mixing

In the case of density mixing, the function we wish to minimise is the difference between the input and output densities in the self-consistent field cycle. The relevant Hessian is the static dielectric matrix for the system, which is why VASP is using that for its reporting -- although note that VASP is not actually computing the dielectric matrix, because that is computationally quite expensive. Instead, it builds up an approximate form for the dielectric matrix over the course of the calculation (see below).

NB it is actually possible to use the dielectric matrix, either by computing it directly (perhaps only for the lowest wavevectors; ABINIT can do this), or using it implicitly by solving the perturbation theory equations iteratively for the dielectric response.

There are two additional methods which are commonly used: preconditioning; and quasi-Newton methods. I will not go into details here, but briefly preconditioning is a technique for transforming your problem into one where the eigenvalues are closer together ("compressing the eigenspectrum"), and quasi-Newton methods are a way to build up a model of $$\mathrm{B}^{-1}$$ in order to capture the change of search direction as well.

In the case of density-mixing methods, preconditioning is commonly done using an approximate inverse dielectric matrix proposed by Manninen et al, although it is usually named after Kerker, who noted its wider applicability. It is based on the asymptotic limit of the dielectric matrix at long wavelengths, where the response of the material is dominated by the electron-electron interaction (the Hartree potential, in density functional theory).

Quasi-Newton methods start with an initial approximation for the dielectric matrix (the Hessian), which is often the Kerker form, and successively improve it using the actual change in the density during the self-consistent field cycles. There are different ways to do this, and the most common are those of Broyden and Pulay (Anderson). In the case of VASP, it is actually this approximate dielectric matrix which is used for the eigenvalue decomposition, since this is much faster to compute.

### References

1. "Electrons and positrons in metal vacancies", M. Manninen, R. Nieminen, P. Hautojärvi, and J. Arponen, Phys. Rev. B 12, 4012 (1975).

2. "Efficient iteration scheme for self-consistent pseudopotential calculations", G. P. Kerker, Phys. Rev. B 23, 3082 (1981).

3. "A class of methods for solving nonlinear simultaneous equations", C. G. Broyden, Math. Comp. 19, 577-593 (1965).

4. "Convergence acceleration of iterative sequences. the case of scf iteration", P. Pulay, Chem. Phys. Lett. 73, 393 (1980).

5. "Improved SCF convergence acceleration", P. Pulay, J. Comput. Chem. 3, 556(1982).

6. "Iterative Procedures for Nonlinear Integral Equations", D. G. Anderson, J. ACM 12 4 (1965).

• +100. Thanks for writing such an excellent and thorough answer again! I also appreciate you staying up late to write it :) Mar 17, 2022 at 2:39
• @NikeDattani my pleasure & thanks for sorting out the numbering & references - I couldn't quite face that by the time I got to the end! Mar 17, 2022 at 16:39
• Instead of letting the average eigenvalue of $\alpha\mathrm{B}$ be 1, there are many seemingly equally justified way of choosing $\alpha$, for example one may let the average eigenvalue of $\alpha^{-1}\mathrm{B}^{-1}$ be 1, or essentially use the average eigenvalue of $\alpha^{n}\mathrm{B}^{n}$ for any non-zero $n$. The resulting $\alpha$ clearly depends on $n$ if $\mathrm{B}$ is not diagonal. Could you elaborate on the significance of choosing $n=1$ here? Mar 20, 2022 at 11:35
• @wzkchem5 I was trying not to stray too far into optimisation theory, but I'll try to add something at the end when I have time. It's actually more complicated, since the choice of 1 is "optimal" in the sense of reducing the mean error, but the whole argument is predicated on remaining in the quadratic well; if eigenvalues are large, you sometimes need a smaller alpha for stability, even though that brings the average below 1. Mar 21, 2022 at 10:21