Direct optimization of the orbitals

With "direct optimisation of the orbitals" we denote a category of electronic minimisation algorithms that use the gradient of the free energy with respect to the orbitals to move towards the ground state of the system: the orbitals are changed such that the total energy is lowered, using, e.g. the conjugate gradient approximation, or damped molecular dynamics.

This gradient of the free energy with respect to an orbital $\psi_n$ is given by:

$| g_n \rangle = f_n \Big(1-\sum^N_{m=1} \hat{S} \vert \psi_m \rangle \langle \psi_m \vert\Big) \hat{H} \vert \psi_n \rangle + \sum^N_{m=1} \frac{1}{2} {\bf H}_{nm} (f_n - f_m) \hat{S} \vert \psi_m \rangle$

where $\{ f_i | i=1,..,N \}$ are the partial occupancies, and

${\bf H}_{nm}=\langle \psi_m \vert \hat{H} \vert \psi_n \rangle$

is the Hamiltonian expressed within the subspace spanned by the current orbitals $\{ \psi_i | i=1,..,N \}$.

The structure of the gradient may be understood as follows: the first part on the right-hand side describes the change of the free energy with respect to changes in the orbitals that are outside (orthogonal) the subspace spanned by the current set of orbitals, whereas the second part describes the changes of the free energy due to a unitary transformation between the orbitals within this subspace.

To derive a search direction, i.e., actual change in the orbitals from the gradient these aforementioned parts are treated separately. The search direction related to the out-of-subspace part of the gradient is:

$\vert p_n \rangle = f_n K \Big(1-\sum^N_{m=1} \hat{S} \vert \psi_m \rangle \langle \psi_m \vert\Big) \hat{H} \vert \psi_n \rangle$

where $K$ is a preconditioning function.

The search direction associated with the subspace rotational part of the gradient may be constructed using Loewdin perturbation theory:

$U_{nm} = \delta_{nm} - \Delta \frac{H_{nm}}{H_{mm}-H_{nn}}$

where $\Delta$ denotes the stepsize.

Per default, however, VASP constructs a search direction for the subspace rotational part of the gradient in the manner proposed by Freysoldt et al. Changes in the partial occupancies are computed in accordance with the aforementioned work as well.

After every change of the orbitals and partial occupancies, the total energy and electronic density are recomputed. Per default, the electronic density is constructed directly from the orbitals and partial occupancies at each step along the way, without any density mixing. Optionally, though, density mixing may be used to stabilise these optimisation procedures when charge sloshing occurs.

The direct optimisation of the orbitals stops when the change of the total energy drops below EDIFF.

Note that, when starting from scratch (ISTART = 0), the direct optimisation procedures in VASP always begin with several (NELMDL) self-consistency cycles where the density is kept fixed at the initial approximation (overlapping atomic charge densities), using the blocked-Davidson algorithm to optimise the orbitals. This ensures that the orbitals, that are initialised with random numbers, have converged to reasonable starting point for the subsequent direct optimisation.

---