Direct optimization of the orbitals

With "direct optimization of the orbitals" we denote a category of electronic minimization 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 _{m=1}^{N}{\hat {S}}\vert \psi _{m}\rangle \langle \psi _{m}\vert {\Big )}{\hat {H}}\vert \psi _{n}\rangle +\sum _{m=1}^{N}{\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 _{m=1}^{N}{\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. Note that taking a step along this search direction amounts to a rotation of the orbitals (rotation on account of $U$ being unitary):

$\vert \psi _{n}\rangle =\sum _{m=1}^{N}U_{nm}\vert \psi _{m}\rangle$

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.

When the three contributions to the "search direction" (out-of-subspace, subspace rotational, and change in the partial occupancies) have been determined, they are used to update the orbitals and partial occupancies, either by steepest descent, by means of the conjugate-gradient approximation, or using a damped molecular dynamics scheme.

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 optimization of the orbitals stops when the change of the total energy drops below EDIFF.

Note that, when starting from scratch (ISTART = 0), the direct optimization 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 optimize the orbitals. This ensures that the orbitals, that are initialised with random numbers, have converged to reasonable starting point for the subsequent direct optimization.

---