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 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{ math\gt {\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 in 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 each band has been optimized a Rayleigh-Ritz optimization in the complete subspace [math]\displaystyle{ \{ \phi_k| k=1,..,N_{\rm bands}\} }[/math] 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.