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^1_k| k=1,..,n_1\}$.

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

$\left \{ \psi^1_k \, / \, g^1_k = \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^1_k \, | \, 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^2_k| k=1,..,n_1\}$

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

$\left \{ \psi^2_k \,/ \, g^1_k \, / \, g^2_k = \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^2_k \, | \, 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^3_k| 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:

$\{\psi^1/g^1/g^2/../g^{d-1}\}$

Per default VASP will not iterate deeper than $d=4$.

• Continue iteration by adding a fourth set of preconditioned vectors if required, etc etc.
• When the iteration is finished, store the optimized block of orbitals back in the set:

$\{ \psi^d_k| k=1,..,n_1\} \Rightarrow \{ \psi_k| k=1,..,N_{\rm bands}\}$.

• Continue with the next block $\{ \psi^1_k| k=n_1+1,..,2 n_1\}$.
• After each band has been optimized a Rayleigh-Ritz optimization in the complete subspace $\{ \phi_k| k=1,..,N_{\rm bands}\}$ is performed.

This method is approximately a factor of 1.5-2 slower than RMM-DIIS, but always stable. It is available in parallel for any data distribution.