Blocked-Davidson algorithm

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

Take a subset (block) of [math]\displaystyle{ n_1 }[/math] orbitals out of the total set of NBANDS orbitals:

[math]\displaystyle{ \{ \psi_n| n=1,..,N_{\rm bands}\}\Rightarrow \{ \psi^1_k| k=1,..,n_1\} }[/math].

Extend the subspace spanned by [math]\displaystyle{ \{\psi^1\} }[/math] by adding the preconditioned residual vectors of [math]\displaystyle{ \{\psi^1\} }[/math]:

[math]\displaystyle{ \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 \}. }[/math]

Rayleigh-Ritz optimization ("subspace rotation") within the [math]\displaystyle{ 2n_1 }[/math]-dimensional space spanned by [math]\displaystyle{ \{\psi^1/g^1\} }[/math], to determine the [math]\displaystyle{ n_1 }[/math] lowest eigenvectors:

[math]\displaystyle{ {\rm diag}\{\psi^1/g^1\} \Rightarrow \{ \psi^2_k| k=1,..,n_1\} }[/math]

Extend the subspace with the residuals of [math]\displaystyle{ \{\psi^2\} }[/math]:

[math]\displaystyle{ \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 \}. }[/math]

Rayleigh-Ritz optimization ("subspace rotation") within the [math]\displaystyle{ 3n_1 }[/math]-dimensional space spanned by [math]\displaystyle{ \{\psi^1/g^1/g^2\} }[/math]:

[math]\displaystyle{ {\rm diag}\{\psi^1/g^1/g^2\} \Rightarrow \{ \psi^3_k| k=1,..,n_1\} }[/math]

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:

[math]\displaystyle{ {\rm diag}\{\psi^1/g^1/g^2/../g^{d-1}\}\Rightarrow \{ \psi^d_k| k=1,..,n_1\} }[/math]

Per default VASP will not iterate deeper than [math]\displaystyle{ d=4 }[/math], 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:

[math]\displaystyle{ \{ \psi^d_k| k=1,..,n_1\} \Rightarrow \{ \psi_k| k=1,..,N_{\rm bands}\} }[/math].

Continue with the next block [math]\displaystyle{ \{ \psi^1_k| k=n_1+1,..,2 n_1\} }[/math].
After all orbitals have been optimized, a Rayleigh-Ritz optimization in the complete subspace [math]\displaystyle{ \{ \psi_k| k=1,..,N_{\rm bands}\} }[/math] is performed.

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