SCF energy keeps on fluctuating between two values, instead of converging, why?

Question

Sometimes I have observed this phenomenon where my calculations' SCF iterations keep on fluctuating between two values, and not converging. I have observed it in all sort of systems, ranging from molecules and materials.

What is the reason for this see-saw behavior and how to alleviate it? My google search was bit inconclusive

I have had similar problems once before in some gold clusters as well, but that project got scrapped so I never properly followed up that time.

CP2K input and output files for replication (I am not the most well versed guy in DFT so please highlight any obvious error as well!).

&GLOBAL
  PROJECT "Sb8 | Energy | SZV-GTH"
  RUN_TYPE ENERGY
  PRINT_LEVEL LOW
&END GLOBAL

&FORCE_EVAL
  METHOD QUICKSTEP

    &DFT
    BASIS_SET_FILE_NAME  BASIS_MOLOPT_UCL
    POTENTIAL_FILE_NAME  GTH_POTENTIALS
    &SCF
      MAX_SCF 300
      EPS_SCF 1e-05
      SCF_GUESS ATOMIC
       &SMEAR ON
        METHOD FERMI_DIRAC
        ELECTRONIC_TEMPERATURE [K] 300
      &END SMEAR
      !&DIAGONALIZATION
      !  ALGORITHM STANDARD
      !&END DIAGONALIZATION
      ADDED_MOS 80
    &END SCF
    &XC
      &XC_FUNCTIONAL PBE
      &END XC_FUNCTIONAL
    &END XC
     &MGRID
      NGRIDS 5
      CUTOFF 650
      REL_CUTOFF 65
    &END MGRID
    &QS
      EPS_RHO      1.0E-12
      EPS_PGF_ORB  1.0E-07
    &END QS
  &END DFT

  &SUBSYS
    &KIND Sb
      ELEMENT Sb
      BASIS_SET TZVP-MOLOPT-SR-GTH-q5
      POTENTIAL GTH-PBE-q5
    &END KIND
    &CELL
      A       8.6168003       0.0000000       0.0000000
      B      -4.3084006       7.4623675       0.0000000
      C      -0.0000010      -0.0000017      22.5480003
    &END CELL 
    &COORD
      Sb      -0.00000     -0.00000      2.63237
      Sb      -0.00000     -0.00000     13.90637
      Sb      -2.15420      3.73118      2.63237
      Sb       2.15420      3.73118     13.90637
    &END COORD
  &END SUBSYS
&END FORCE_EVAL

 Step     Update method      Time    Convergence         Total energy    Change
  ------------------------------------------------------------------------------
...
    46 P_Mix/Diag. 0.40E+00    4.1     0.83462647       -21.3544184161  1.19E-07
    47 P_Mix/Diag. 0.40E+00    4.1     0.83462647       -21.3544185344 -1.18E-07
    48 P_Mix/Diag. 0.40E+00    4.1     0.83462647       -21.3544184158  1.19E-07
    49 P_Mix/Diag. 0.40E+00    4.2     0.83462647       -21.3544185342 -1.18E-07
    50 P_Mix/Diag. 0.40E+00    4.2     0.83462647       -21.3544184156  1.19E-07
    51 P_Mix/Diag. 0.40E+00    4.3     0.83462647       -21.3544185342 -1.19E-07
    52 P_Mix/Diag. 0.40E+00    4.2     0.83462647       -21.3544184156  1.19E-07
    53 P_Mix/Diag. 0.40E+00    4.1     0.83462647       -21.3544185342 -1.19E-07
    54 P_Mix/Diag. 0.40E+00    4.2     0.83462647       -21.3544184156  1.19E-07
    55 P_Mix/Diag. 0.40E+00    4.2     0.83462647       -21.3544185341 -1.19E-07
    56 P_Mix/Diag. 0.40E+00    4.2     0.83462647       -21.3544184156  1.19E-07
    57 P_Mix/Diag. 0.40E+00    4.2     0.83462647       -21.3544185341 -1.19E-07
    58 P_Mix/Diag. 0.40E+00    4.2     0.83462647       -21.3544184156  1.19E-07
    59 P_Mix/Diag. 0.40E+00    4.2     0.83462647       -21.3544185341 -1.19E-07
    60 P_Mix/Diag. 0.40E+00    4.2     0.83462647       -21.3544184156  1.19E-07
    61 P_Mix/Diag. 0.40E+00    4.2     0.83462647       -21.3544185341 -1.19E-07
...

Answer 1

The cause of the oscillations in the SCF energies is probably a "sloshing instability"; the two most common types are "charge sloshing" and "occupancy sloshing".

Sloshing instabilities

Sloshing instabilities arise in the self-consistent field (SCF) method because the trial solutions to the Kohn-Sham equations are optimised for a fixed Kohn-Sham potential, and then the potential is updated. There are many ways to think about why this happens, but for a physical picture consider this scenario with a very simple SCF method:

At the start of some iteration $i$ of the SCF optimisation, one region (A) of the system has too high an electron density, and another region (B) has too low an electron density.

We construct the Kohn-Sham potential for this iteration, $V^{(i)}(r)$, from the density. The electrostatic repulsion (Hartree term) means that $V^{(i)}(r)$ is very large in region A, whereas it is low in region B.

The Kohn-Sham states are optimised by minimising the energy of the system, keeping the potential $V^{(i)}(r)$ fixed. Electronic density is moved from region A (high potential) to region B (low potential).

Note: as electronic density is removed from A, the potential in that region should be reduced; conversely, the potential in region B should increase as electron density is moved into that region. However, since the optimisation is performed for a fixed potential, this doesn't happen in the simulation. The result is that the optimisation moves more electron density from A to B than it should.

Once the new optimised states have been computed, the new electron density is calculated. Since too much density was moved from A to B, the electron density is now too low in region A and too high in region B.

The new Kohn-Sham potential is calculated from the density, and used as the input to the next iteration, $V^{(i+1)}(r)$. The electrostatic repulsion means that the $V^{(i+1)}(r)$ is low in region A, and very large in region B.

The next optimisation is performed for fixed $V^{(i+1)}(r)$, and now electron density is moved from region B (where the potential is now high) back to region A (because its potential is now low). We are now back where we started!

Stabilising convergence

Sloshing instabilities are extremely common, and in fact the simple SCF method in the previous section almost never converges. Real-life software doesn't just take the output of one SCF iteration as the input to the next, it "damps" the changes by mixing together the outputs from two or more iterations. These methods are called "density mixing" or "potential mixing", depending on whether it is the density or potential which is mixed. I'll assume density mixing here for brevity, but you can just replace "density" with "potential" for potential mixing.

Most mixing methods are based on the behaviour of a parameterised model system (e.g. the homogeneous electron gas), with a method to dynamically update the model (e.g. Pulay/Anderson or Broyden). The key user parameter is the total amount of the density change to "mix in" to the previous input, usually defined such that 0 means don't change the density at all, and 1 means just use the new density (the simple SCF method from the previous section). In CP2K the relevant parameter is DFT/SCF/MIXING/ALPHA.

The smaller you set this mixing parameter, the smaller the changes you allow the software to make each SCF iteration. This has the effect of making the SCF more stable, but could slow the convergence. Conversely, a large mixing parameter makes much larger changes to the density, so the SCF could converge very rapidly, but could also become unstable.

Summary

What you're seeing is a common instability in SCF methods, and you can often suppress it by reducing the mixing amplitude parameter. For a fuller discussion of density mixing, see:

How does charge mixing work?

and a related discussion on the mixing amplitude and stability:

The relationship between average eigenvalue and convergence performance in VASP?

One of my former students, Nick Woods, looked into this in some detail, so if you'd like to know more you might find his Masters thesis and paper useful:

"On the Nature of Self-Consistency in Density Functional Theory", N. Woods; https://arxiv.org/abs/1803.01763

"Computing the Self-Consistent Field in Kohn-Sham Density Functional Theory", N. Woods, P. Hasnip, M. C. Payne, J. Phys.: Condens. Matter, 31, 453001 (2019); https://doi.org/10.1088/1361-648X/ab31c0