How does charge mixing work?

Question

When doing DFT calculations on certain magnetic materials, convergence can be difficult. This problem can be solved by changing the mixing parameters (e.g. AMIN and BMIX in the INCAR file of VASP). My question is how charge mixing occurs and what is the underlying reason for convergence issues in magnetic systems ?

Answer 1

Density mixing is a type of self-consistent field (SCF) method, which tries to find the closest density to the Kohn-Sham (or Hartree-Fock etc) ground state density by mixing previous densities and density changes.

The self-consistent field method

In the density mixing SCF method, a trial density $n_\mathrm{in}(\vec{r})$ is used to construct the Kohn-Sham Hamiltonian, and then a set of trial states $\left\{\psi_{bk}\right\}$ is computed by solving the Kohn-Sham equations,

$$ \hat{H}\left[n_\mathrm{in}\right]\psi_{bk}=E_{bk}\psi_{bk}. \tag{1} $$ These equations may be solved either by conventional matrix diagonalisation (if the basis set is compact, e.g. atom-centred Gaussians) or by iterative diagonalisation (if the basis set is large, e.g. plane-waves, such as VASP uses).

However equation (1) is solved, an important thing to note is that the density $n_\mathrm{in}$ is not changed as the states are updated; yet the density should be related to the states by,

$$ n(\vec{r}) = \sum_{bk}f_{bk}\left\vert\psi_{bk}(\vec{r})\right\vert^2, \tag{2} $$ where $f_{bk}$ are the band occupancies. We define this density to be $n_\mathrm{out}$.

Only when $n_\mathrm{out}=n_\mathrm{in}$ have we computed the correct Kohn-Sham ground state, where the density obtained from our states is the same as the density used to define the Hamiltonian; the solution in this case is said to be self-consistent.

In general, $n_\mathrm{out}\neq n_\mathrm{in}$, in which case the $\psi_{bk}$ we computed are not for the ground state Hamiltonian. The question thus arises: given an initial input density $n^{(0)}_\mathrm{in}$ and resultant output density $n^{(0)}_\mathrm{out}$, how can we find an improved guess for the input density, $n^{(1)}_\mathrm{in}$?

Let us define the change in $n^{(0)}_\mathrm{in}$ to be $\delta n^{(0)}_\mathrm{in}$, and write $$ n^{(1)}_\mathrm{in} = n^{(0)}_\mathrm{in} + \delta n^{(0)}_\mathrm{in}. \tag{3} $$ Similarly, we write the new output density as $$ n^{(1)}_\mathrm{out} = n^{(0)}_\mathrm{out} + \delta n^{(0)}_\mathrm{out}. \tag{4} $$ We wish these two densities to be the same, since that corresponds to the self-consistent solution of the Kohn-Sham equations. Thus, $$ n^{(0)}_\mathrm{in} + \delta n^{(0)}_\mathrm{in} = n^{(0)}_\mathrm{out} + \delta n^{(0)}_\mathrm{out}. \tag{5} $$ If we are near to the solution, then we can approximate $\delta n_\mathrm{out}$ as, $$ \delta n_\mathrm{out} \approx \frac{\delta n_\mathrm{out}}{\delta n_\mathrm{in}}\delta n_\mathrm{in}, \tag{6} $$ which we substitute into equation (5) to obtain, $$ \begin{array}{lrcl} & n^{(0)}_\mathrm{in} + \delta n^{(0)}_\mathrm{in} &\approx& n^{(0)}_\mathrm{out} + \frac{\delta n_\mathrm{out}}{\delta n_\mathrm{in}}\delta n_\mathrm{in}\\ \Rightarrow & \delta n^{(0)}_\mathrm{in} - \frac{\delta n_\mathrm{out}}{\delta n_\mathrm{in}}\delta n_\mathrm{in} &\approx& n^{(0)}_\mathrm{out} - n^{(0)}_\mathrm{in}\\ \Rightarrow & \left(1 - \frac{\delta n_\mathrm{out}}{\delta n_\mathrm{in}}\right)\delta n_\mathrm{in} &\approx& n^{(0)}_\mathrm{out} - n^{(0)}_\mathrm{in}\\ \Rightarrow & \delta n_\mathrm{in} &\approx& \left(1 - \frac{\delta n_\mathrm{out}}{\delta n_\mathrm{in}}\right)^{-1} \left(n_\mathrm{out}-n_\mathrm{in}\right). \tag{7} \end{array} $$ The quantity $\left(n_\mathrm{out}-n_\mathrm{in}\right)$ is the density residual, which is the difference between the input and output densities; at the solution, this should be zero.

The term $$ \left(1 - \frac{\delta n_\mathrm{out}}{\delta n_\mathrm{in}}\right)^{-1} $$ is a matrix inversion. This may not be immediately apparent, since I omitted the positional variables for brevity; re-inserting them, we have, $$ \left(1 - \frac{\delta n_\mathrm{out}(\vec{r})}{\delta n_\mathrm{in}(\vec{r}^\prime)}\right)^{-1}, \tag{8} $$ which reflects the fact that a change in the input density at a point $\vec{r}^\prime$ can cause a change in the output density at a different point, $\vec{r}$.

We don't know what this matrix is a priori, so how should we proceed? One physics point to note is that the response of the density ($\delta n_\mathrm{out}$) to a perturbing potential (which is what changing the input density by $\delta n_\mathrm{in}$ does) is the dielectric response, and this matrix is the inverse dielectric tensor of the system.

Linear mixing

As a first assumption, we may simply approximate the matrix by a scalar times the identity, $$ \left(1 - \frac{\delta n_\mathrm{out}(\vec{r})}{\delta n_\mathrm{in}(\vec{r}^\prime)}\right)^{-1} \approx \alpha \mathrm{I}. \tag{9} $$ Inserting this into equation (7) yields, $$ \begin{array}{lrcl} & \delta n_\mathrm{in} & \approx & \alpha \left(n_\mathrm{out}-n_\mathrm{in}\right)\\ \Rightarrow & n^{(1)}_\mathrm{in} &=& n^{(0)}_\mathrm{in} + \delta n^{(0)}_\mathrm{in} \\ & & \approx & n^{(0)}_\mathrm{in} + \alpha \left(n^{(0)}_\mathrm{out}-n^{(0)}_\mathrm{in}\right) \\ & & = & (1-\alpha) n^{(0)}_\mathrm{in} + \alpha n^{(0)}_\mathrm{out}. \tag{10} \end{array} $$ We see that the new input density,$n^{(1)}_\mathrm{in}$ is a linear mixture of the previous input and output densities, and it is usually referred to simply as linear mixing.

Linear mixing is very simple and easy to implement. For small systems it can be surprisingly effective, but it quickly becomes unstable for systems larger and/or more complex.

Dielectric preconditioning

The problems with linear mixing are largely due to the crude approximation of a matrix by a scaled identity. We could compute the dielectric tensor (e.g. with density functional perturbation theory), but since the system may be quite a long way from the ground state in these early SCF cycles, that may require a lot of computational effort which must be redone every SCF cycle.

A popular alternative is to approximate the dielectric by that of the homogeneous electron gas (also known as jellium). The dielectric response of jellium is diagonal in reciprocal-space, so it is convenient to express it not in terms of $\vec{r}$ and $\vec{r}^\prime$, but in terms of the reciprocal lattice vectors $\vec{G}$ and $\vec{G}^\prime$,

$$ \left(1 - \frac{\delta n_\mathrm{out}(\vec{G})}{\delta n_\mathrm{in}(\vec{G}^\prime)}\right)^{-1} \approx \delta_{\vec{G}\vec{G}^\prime}\frac{\vert\vec{G}\vert^2}{\vert\vec{G}\vert^2+\vert G_0\vert^2}, \tag{11} $$ where the parameter $G_0$ is related to the Thomas-Fermi screening wavevector, and $\delta_{\vec{G}\vec{G}^\prime}$ is the Kronecker delta (i.e. it means the matrix is diagonal).

The important thing to note about this expression is that for small wavevectors it is approximately $\frac{\vert\vec{G}\vert^2}{\vert G_0\vert^2}$, which is small. Thus, it has the effect of damping the small $\vec{G}$ components of the residual and, since a small reciprocal-space vector corresponds to large real-space changes, this means it suppresses the long-range density changes which are especially problematic when simulating metals.

It is common practice to also include a mixing parameter $\alpha$, so that the actual density is given by,

$$ \begin{array}{lrcl} & n^{(1)}_\mathrm{in}(\vec{G}) &=& n^{(0)}_\mathrm{in}(\vec{G}) + \delta n^{(0)}_\mathrm{in}(\vec{G}) \\ & & \approx & n^{(0)}_\mathrm{in}(\vec{G}) + \frac{\alpha\vert\vec{G}\vert^2}{\vert\vec{G}\vert^2+\vert G_0\vert^2} \left(n^{(0)}_\mathrm{out}(\vec{G})-n^{(0)}_\mathrm{in}(\vec{G})\right). \tag{12} \end{array} $$

To my knowledge, this approach was first used by J Arponen, P Hautojärvi, R Nieminen and E Pajanne in their studies of electron and positron defects in metals, but it was Kerker who proposed it as a general method for self-consistent electronic-structure calculations (in a comment on a later M Manninen, R Nieminen, P Hautojärvi and J Arponen paper), and this approach is usually known as Kerker mixing.

Quasi-Newton methods

As the calculation proceeds, we see how the system actually responded to the new input densities, and so we gain more and more information about the actual dielectric response of the system. There is a whole field of optimisation theory to draw on when approximating a matrix by its application to sample vectors (in our case, the density residuals). The most widespread methods are probably those of Broyden and Pulay (or, equivalently, Anderson).

All of these methods involve constructing a matrix in the basis of the densities and residuals your system has used, and so the result is that your new input density is usually a mixture of all previous densities.

For plane-wave based methods, these methods are almost always combined with Kerker preconditioning so that effectively the quasi-Newton method builds up a model of the difference between the actual dielectric response, and the jellium model.

Spin and magnetism

Now we come to the heart of your question: why are magnetic systems harder to converge? There are several reasons, but probably the two most generally problematic are:

• Magnetic systems typically have many local energy minima
  These minima correspond to different atomic magnetic states; quasi-Newton methods often perform poorly for multiple minima problems

• The Kerker model is less appropriate
  The Kerker model is based on jellium's charge response, but it is usually used to precondition the spin mixing as well (albeit with different parameters for the mixing amplitude and screening wavevector). Using the charge dielectric model for the spin mixing is generally found to be an improvement over linear mixing, but it is not clear that it is actually an appropriate model.

The Kerker model may also be derived from the Kohn-Sham equations directly, assuming that the change in potential is dominated by the change in the Hartree potential; however, the Hartree potential has no contribution from the spin at all so it is not at all clear that the same model will be a good approximation to the spin response.

Other approaches

I think I have covered all of the key machinery VASP uses, but these aren't the only possible approaches. For a fuller discussion of the different density mixing methods, and the causes of instabilities, you might find this paper from my group useful, and the discussion of "difficult cases" here. This focuses on the quasi-Newton part of the problem, not the dielectric preconditioning.

Many quantum chemistry packages use a density matrix mixing method, which has some additional advantages, since the energy can be expressed directly as the trace of the density matrix times the Hamiltonian. These were not covered in the above paper, as we did that in the CASTEP plane-wave program and density matrix methods are very awkward to implement in a plane-wave basis.

Some packages also use more sophisticated models for the dielectric, for example computing the dielectric matrix (perhaps in a reduced basis, or implicitly via perturbation theory), or the auxiliary Thomas-Fermi model of Raczkowski and Canning.