next up previous contents
Next: Geometry Optimization Up: Practical Electronic Structure Calculations Previous: Practical Electronic Structure Calculations   Contents


Single Point Energy Calculations

A single point energy (SPE) calculation calculates the wavefunction and charge density, and hence the energy, of a particular (arbitrary) arrangement of nuclei. The total energy functional of a system of electrons and nuclei can be written as

\begin{displaymath}
E[\rho(\mathbf{r_i}),\mathbf{R}_I]=
T_e[\rho(\mathbf{r_i})...
...thbf{R}_I]+
E_{ee}[\rho(\mathbf{r_i})]+
E_{nn}(\mathbf{R}_I)
\end{displaymath} (2.32)

where $T_e[\rho(\mathbf{r_i})]$ is the electronic kinetic energy given in (2.12), $T_n$ is the nuclear kinetic energy, $E_{en}$ is the electron-nuclei interaction given in (2.9), $E_{ee}$ is the electron-electron interaction and $E_{nn}$ is the nuclei-nuclei interaction. $E_{ne}$ is described using pseudopotentials §2.5 while for a fixed set of nuclear positions $T_n$ is zero and $E_{nn}$ is a constant. Thus a SPE calculation is reduced to finding the charge density that minimizes the total energy functional, i.e. solving the Kohn-Sham equations.

The minimization of the total energy functional can be done in a number of ways. One method proceeds by direct diagonalization of the matrix equation (2.30). Starting from an initial trial density $E_{ee}$ is calculated and inserted into (2.30). A new density is then calculated by inverting the matrix equation (2.30). If this new density gives an energy that is consistent with the old density (if the change in energy between iterations is smaller than a given tolerance) it is then inserted into the total energy functional and the energy calculated. Otherwise this new density is used to calculate a new $E_{ee}$. This is repeated until the density and potential are consistent with each other, within a given tolerance.

The matrix diagonalization method has the disadvantage that the computational cost of matrix diagonalization scales as the number of plane waves cubed. An alternative method is to minimize the energy functional directly [46]. The energy functional is a functional of the density, which is determined by the expansion coefficients $c_{i,\mathbf{k}+\mathbf{G}}$. The ground state density is found from the set of $c_{i,\mathbf{k}+\mathbf{G}}$ that minimize the energy functional. Standard functional optimization methods [47] can then be used to find the minima of the total energy functional. The simplest of these techniques is the method of Steepest Descents (SD) [47]. The SD method produces a series of points $\{\mathbf{P}_i\}$

\begin{displaymath}
\mathbf{P}_{i+1}=\mathbf{P}_i+\lambda_i\mathbf{h}_i
\end{displaymath} (2.33)

where
\begin{displaymath}
\mathbf{h}_i=-\nabla f(\mathbf{P}_i).
\end{displaymath} (2.34)

In this case $\mathbf{P}_i$ are sets of plane wave expansion coefficients and $f(\mathbf{P})$ is the energy functional. The scalar $\lambda_i$ in (2.33) is the distance along the direction $\mathbf{h}_i$ from $\mathbf{P}_i$ that a minima is located. The SD method proceeds by moving in the steepest downhill direction from a point $\mathbf{P}_i$ until a minima along that direction is located at point $\mathbf{P}_{i+1}$. The steepest downhill direction from $\mathbf{P}_{i+1}$ is then determined ( $-\nabla
f(\mathbf{P}_{i+1})$) and the minima located along that direction is found. This is repeated until the change in the function is lower than a preset tolerance. The speed at which the SD method will find a minima is limited as at each step only the information at that point is taken into account. It is easy to think of examples where this will lead to slow convergence. A better method is the Conjugate Gradients (CG) method [46,48,47]. It differs from the SD method in that each search direction is conjugate to the last one. The CG method assumes that the function can be approximated by a quadratic form
\begin{displaymath}
f(\mathbf{x})\approx\frac{1}{2}\mathbf{x}.\mathbf{A}.\mathbf{x}
\end{displaymath} (2.35)

where $\mathbf{A}$ is the Hessian given by
\begin{displaymath}
A_{ij}=
\left.\frac{\partial^2f}{\partial x_i\partial x_j}
\right\vert _{\mathbf{x}=\mathbf{x}_0}.
\end{displaymath} (2.36)

The CG method produces a series of points as in (2.33) but $\mathbf{h}_i$ is given by
\begin{displaymath}
\mathbf{h}_i=
\left\{\begin{array}{c}
\mathbf{g}_i,\:i=0 ...
...\mathbf{g}_{i-1}}
\mathbf{h}_{i-1},\:i>0
\end{array} \right.
\end{displaymath} (2.37)

where $\mathbf{g}_i=-\nabla f(\mathbf{P}_i)$. Using the quadratic form (2.35) the step needed to get to the minima in the direction $\mathbf{h}_i$ is given by [47]
\begin{displaymath}
\lambda_i=\frac{\mathbf{h}_i.\mathbf{g}_i}
{\mathbf{h}_i.\mathbf{A}.\mathbf{h}_i}.
\end{displaymath} (2.38)

This however can be found from a one-dimensional minimization along $\mathbf{h}_i$, i.e. without explicitly calculating the Hessian. As the Hessian is an $N\times N$ matrix, for DFT calculations with $N=10^5-10^6$ plane waves the CG method has a large advantage over methods that explicitly use the Hessian. The CG method will find a minimum of an $N$ dimensional function in $N$ iterations. When using the CG and SD methods to minimize the total energy functional , the requirement that the wavefunctions are to be orthonormal places an additional constraint on the minimzation. For a DFT calculation with 10$^5$ plane waves this is still a monumental task. The number of iterations required can be substantially reduced by preconditioning the function [46].


next up previous contents
Next: Geometry Optimization Up: Practical Electronic Structure Calculations Previous: Practical Electronic Structure Calculations   Contents
Dr S J Clark
2003-01-30