Difficult cases for converging Kohn-Sham SCF calculations

Question

In a recent article by Woods et. al.^[1] a couple of methods for converging ground-state Kohn-Sham DFT calculations are reviewed and compared. In their test suite they give plenty of examples for badly convergent cases for self-consistent field calculations. Taking a look at their examples it seems that mainly isolated atoms, large cells, slabs and unusual spin systems that cause trouble with convergence.

I am wondering what the experience of the community are with respect to cases where standard approaches (e.g. Kerker mixing, some temperature) fail. Did you encounter systems where you had a tough time? What type / class of system were they? What are common approaches to circumvent the problem or improve convergence?

My interest mainly concerns with numerical convergence problems in the self-consistent field procedure itself. So when one would expect Kohn-Sham DFT to be a sensible model for a particular system, making the SCF problem well-posed, but still the SCF fails to converge (like e.g. the slab of Gold they mention in the cited paper).

References:

1. Woods, N. D.; Payne, M. C.; Hasnip, P. J. Computing the self-consistent field in Kohn–Sham density functional theory. J. Phys.: Condens. Matter 2019, 31 (45), 453001. DOI: 10.1088/1361-648X/ab31c0.

Answer 1

A few materials/simulation boxes I've had some proper trouble with:

1. HSE06 + noncollinear magnetism + antiferromagnetism, Vasp noncollinear:

   • This was a strongly antiferro material (4 Fe atoms, in an up-down-up-down configuration).
   • HSE06 is apparently difficult to converge anyway.
   • Noncollinear magnetism/antiferromagnetism apparently creates problems for any charge density/spin density mixer (or so I've been told/seen). But this was quite bad.
   • Solved with: AMIX = 0.01, BMIX=1e-5, AMIX_MAG=0.01, BMIX_MAG=1e-5, Methfessel-Paxton order 1 smearing of 0.2 eV. Davidson solver (ALGO=Fast). It took ~160 SCF steps, but it did converge.

2. Cell with really different a, b, c, GPAW, PW mode:

   • This was a much simpler, spin-paired, metallic system. However, the cell was 5.8 x 5.0 x ~70 angstroms.
   • This happens in general when a cell is really elongated along a particular axis/in general has a very 'non-cubic' shape.
   • This paper, which I believe is the reference for Quantum Espresso's 'local-TF' charge-density mixing explains (I think) why this happens. Turns out that really large lattice vector (or in general a very 'non-cubic' cell) ill-conditions the charge-mixing problem. This is the precise problem that is addressed in the paper.
   • Since GPAW does not implement the 'local-TF' mixing yet, I solved this by using mixer=Mixer(beta=0.01). The convergence is very slow, but again, it did converge.

3. A colleague and I discussed the convergence of a (possibly antiferromagnetic) nickel compound (that's about as much I know), Vasp. I don't have many details on this other than it was a similar pain as case #1 above. It took a lot of work and mostly just turning down AMIX and AMIX_MAG.

4. I recently saw another person attempting to converge a single Ni atom in a box (GPAW, LCAO mode), with magnetic moment set by Hund's rule, and I think the thing actually converged density to a log10-error of < -2.4 (you want to aim for < -4.0). To converge further, they proceeded with a Fermi-Dirac smearing of... 0.5 eV. This was one of the more extreme cases I've seen.

I will keep updating this answer as I see/find more pathological cases, but these are what I have so far.

Answer 2

This will be a long answer, so I will divide it in parts.

---

Woods paper

A significant limitation of the Woods et al paper is that it excludes atomic-basis set calculations where convergence acceleration is much more powerful than in plane wave codes. Namely, the update schemes discussed in the article talk about just the input and output densities, whereas if you can store and diagonalize the Kohn-Sham-Fock matrix, you can formulate much faster converging methods for the solution. Typical quantum chemistry codes extrapolate the Fock matrix, not the density. This method typically achieves convergence in a few dozen iterations.

---

Which spin state?

As far as I know, some solid-state codes determine the spin state on the fly. If you don't fix the spin multiplicity, this may contribute to convergence problems. A study of convergence problems should be run for a fixed spin state; one can always carry out separate calculations for each spin state.

---

What is "Kohn-Sham"?

I also have to point out that the notion of "Kohn-Sham calculations" is a bit ill-defined, since typical solid-state calculations are run at a finite temperature; I think this is typically referred to as Mermin-Kohn-Sham theory. Kohn-Sham to me means integer occupations. If you have a finite temperature, you get fractional occupations.

---

Fractional occupations

Now, you often get convergence problems when you have solutions of different symmetries close together (which is why atoms and diatomic molecules are often challenging). The reason for the lack of convergence is that the occupations switch between the SCF cycles. In some cases you can even find that LUMO and HOMO swap places when you optimize the orbitals: you find that the LUMO is below the HOMO, you reoptimize the orbitals at this symmetry, and now you realize the new LUMO is below the new HOMO.

But, fractional occupations at the Fermi level are in principle allowed by the Aufbau scheme. Allowing fractional occupations helps in this case, and you get much better convergence.

However, variational minimization of the energy with respect to both the orbitals and the fractional occupation numbers is very hard (which is why AFAIK almost nobody does it).

---

Finite temperature

The alternative is to use a smearing function, e.g. Fermi-Dirac occupations. Also in this case the occupation numbers depend on the orbital energies, which depend on the orbitals, which depend on the occupation numbers. Solving the coupling between these might also make calculations slowly convergent; I am not sure how tightly these are converged in solid state codes.

Answer 3

It is well-established that it can be quite difficult to converge the SCF for certain (hybrid) meta-GGA functionals compared to their (hybrid) GGA counterparts. This is especially true for plane-wave periodic DFT and is most often the case with many of the popular Minnesota functionals. In large part because of this reason, a revised version of M06-L functional was developed. Anecdotally, I have studied several systems in VASP where it is near-impossible to get the SCF to converge with M06-L unless starting from the PBE wavefunction. Generally, there are also algorithmic tricks that can be used to help the convergence of meta-GGA functionals though. In VASP, using a pre-conditioned conjugate gradient algorithm (algo=All) is found to greatly improve convergence compared to the other commonly used algorithms. This is discussed in the VASP wiki.

This may seem to silly, but it is also worth mentioning that a very common case of the SCF converging poorly is when your system isn't appropriately set up! I have personally encountered many scenarios where the SCF is extremely difficult to converge, only to find out the crystal structure I took from a database was missing an atom or ion such that the system was no longer charge-neutral. Of course, this is somewhat of an artificial scenario, but if you run into such troubles, it is one worth considering.

For the purposes of this discussion, I'd also like to add cases where the SCF is converged but not necessarily to the ground state solution.

One very tricky case is related to spin. For a given spin multiplicity (i.e. net magnetic moment), there are some systems where there can be multiple relatively low-lying arrangements of electrons. You can fix the net spin multiplicity, but there's no guarantee you converge to the lowest energy arrangement of electrons at that spin multiplicity. In these cases, the user must provide an initial guess for the individual magnetic moments and hope the SCF converges to the desired electronic state. Occasionally, it can be very difficult to converge to the lowest energy state even if the next lowest electron configuration is several tens of kJ/mol higher in energy.

Another scenario is with DFT+$U$, which can have multiple self-consistent solutions of different energies. There are several cases regarding bulk oxides where the SCF converges to a solution that is not the ground state, such as for $\mathrm{UO_{2}}$ in the fluorite structure. Here is a discussion on how to get around that.

Answer 4

Here's some calculations I had some problems with:

LaFeO$_3$ on LaAlO$_3$ with an adsorbed O atom. $\sqrt{2}\times\sqrt{2}$ perovskite cell in the x-y direction. 5 layers of LaAlO$_3$ (alternating LaO - AlO$_2$ - LaO etc.) with 3 layers of LaFeO$_3$ on top. The bottom 3 layers of LaAlO$_3$ are fixed. Anti-ferromagnetic order on the Fe atoms. Around 450 electrons or so. I used Quantum Espresso 6.1.

Things I played with to get convergence:

• trying a dipole correction OR the effective screening medium method
• changing the amount of vacuum in the unit cell z direction (usually 16-20 Å)
• mixing_mode = 'local-TF' with varying number of iterations (mixing_dim of 8, 10, 12)
• lowered mixing_beta (~0.2)
• removing random atomic orbitals from the initial charge density guess (startingwfc = 'atomic' instead of 'atomic+random')
• different pseudopotential (sometimes pseudopotentials that behaved extremely similarly during early testing could result in very subtle and fickle changes in the tendency to converge)
• slightly perturbing atomic positions (during relaxation, sometimes things would converge for many steps, then you arrive at a certain geometry that throws a wrench in scf convergence). I used scf_must_converge = .false. here.
• I also tried using a fork of QE by Satomichi Nishihara, who implemented RMM-DIIS diagonalization and SR1-BFGS relaxation algorithms

Eventually some combination of the above worked. What a pain!