Blocked-Davidson algorithm

The workflow of the blocked-Davidson iterative matrix diagonalization scheme implemented in VASP is as follows:

• Take a subset (block) of $n_{1}$ orbitals out of the total set of NBANDS orbitals:

$\{\psi _{n}|n=1,..,N_{\rm {bands}}\}\Rightarrow \{\psi _{k}^{1}|k=1,..,n_{1}\}$.

• Extend the subspace spanned by $\{\psi ^{1}\}$ by adding the preconditioned residual vectors of $\{\psi ^{1}\}$:

$\left\{\psi _{k}^{1}\,/\,g_{k}^{1}=\left(1-\sum _{n=1}^{N_{\rm {bands}}}|\psi _{n}\rangle \langle \psi _{n}|{\bf {S}}\right){\bf {K}}\left({\bf {H}}-\epsilon _{\rm {app}}{\bf {S}}\right)\psi _{k}^{1}\,|\,k=1,..,n_{1}\right\}.$

• Rayleigh-Ritz optimization ("subspace rotation") within the $2n_{1}$-dimensional space spanned by $\{\psi ^{1}/g^{1}\}$, to determine the $n_{1}$ lowest eigenvectors:

${\rm {diag}}\{\psi ^{1}/g^{1}\}\Rightarrow \{\psi _{k}^{2}|k=1,..,n_{1}\}$

• Extend the subspace with the residuals of $\{\psi ^{2}\}$:

$\left\{\psi _{k}^{2}\,/\,g_{k}^{1}\,/\,g_{k}^{2}=\left(1-\sum _{n=1}^{N_{\rm {bands}}}|\psi _{n}\rangle \langle \psi _{n}|{\bf {S}}\right){\bf {K}}\left({\bf {H}}-\epsilon _{\rm {app}}{\bf {S}}\right)\psi _{k}^{2}\,|\,k=1,..,n_{1}\right\}.$

• Rayleigh-Ritz optimization ("subspace rotation") within the $3n_{1}$-dimensional space spanned by $\{\psi ^{1}/g^{1}/g^{2}\}$:

${\rm {diag}}\{\psi ^{1}/g^{1}/g^{2}\}\Rightarrow \{\psi _{k}^{3}|k=1,..,n_{1}\}$

• If need be the subspace may be extended by repetition of this cycle of adding residual vectors and Rayleigh-Ritz optimization of the resulting subspace:

${\rm {diag}}\{\psi ^{1}/g^{1}/g^{2}/../g^{d-1}\}\Rightarrow \{\psi _{k}^{d}|k=1,..,n_{1}\}$

Per default VASP will not iterate deeper than $d=4$, though it may break off even sooner when certain criteria that measure the convergence of the orbitals have been met.

• When the iteration is finished, store the optimized block of orbitals back into the set:

$\{\psi _{k}^{d}|k=1,..,n_{1}\}\Rightarrow \{\psi _{k}|k=1,..,N_{\rm {bands}}\}$.

• Move on to the next block $\{\psi _{k}^{1}|k=n_{1}+1,..,2n_{1}\}$.
• When LDIAG=.TRUE. (default), a Rayleigh-Ritz optimization in the complete subspace $\{\psi _{k}|k=1,..,N_{\rm {bands}}\}$ is performed after all orbitals have been optimized.

The blocksize $n_{1}$ used in the blocked-Davidson algorithm can be set by means of the NSIM tag. In principle $n_{1}=2\times$ NSIM, but for technical reasons it needs to be dividable by an integer N:

$n_{1}={\rm {int}}\left({\frac {2*{\rm {NSIM}}+N-1}{N}}\right)N$

where $N$ is the "number of band groups per k-point group":

$N={\frac {\rm {\#\;of\;MPI\;ranks}}{{\rm {IMAGES}}*{\rm {KPAR}}*{\rm {NCORE}}}}$

(see the section on parallelization basics).

As mentioned before, the optimization of a block of orbitals is stopped when either the maximum iteration depth (NRMM), or a certain convergence threshold has been reached. The latter may be fine-tuned by means of the EBREAK, DEPER, and WEIMIN tags. Note: we do not recommend you to do so! Rather rely on the defaults instead.

The blocked-Davidson algorithm is approximately a factor of 1.5-2 slower than the RMM-DIIS, but more robust.

---