nearly finished pseudopotentials
[lectures/latex.git] / posic / thesis / basics.tex
index 46f871c..231345a 100644 (file)
@@ -277,7 +277,7 @@ It provides a stable algorithm that allows smooth changes of the system to new v
 \section{Denstiy functional theory}
 \label{section:dft}
 
-Dirac declared that chemistry has come to an end, its content being entirely contained in the powerul equation published by Schr\"odinger in 1926 \cite{schroeder26} marking the beginning of wave mechanics.
+Dirac declared that chemistry has come to an end, its content being entirely contained in the powerul equation published by Schr\"odinger in 1926 \cite{schroedinger26} marking the beginning of wave mechanics.
 Following the path of Schr\"odinger the problem in quantum-mechanical modeling of describing the many-body problem, i.e. a system of a large amount of interacting particles, is manifested in the high-dimensional Schr\"odinger equation for the wave function $\Psi({\vec{R}},{\vec{r}})$ that depends on the coordinates of all nuclei and electrons.
 The Schr\"odinger equation contains the kinetic energy of the ions and electrons as well as the electron-ion, ion-ion and electron-electron interaction.
 This cannot be solved exactly and finding approximate solutions requires several layers of simplification in order to reduce the number of free parameters.
@@ -329,13 +329,15 @@ E_0=\min_{n(\vec{r})}
   F[n(\vec{r})] + \int n(\vec{r}) V(\vec{r}) d\vec{r}
  \right\}
  \text{ ,}
+\label{eq:basics:hkm}
 \end{equation}
 where $F[n(\vec{r})]$ is a universal functional of the charge density $n(\vec{r})$, which is composed of the kinetic energy functional $T[n(\vec{r})]$ and the interaction energy functional $U[n(\vec{r})]$.
-The challenging problem of determining the exact ground-state is now formally reduced to the determination of the $3$-dimensional function $n(\vec{r})$ via a well-defined but not explicitly known functional of the charge density.
+The challenging problem of determining the exact ground-state is now formally reduced to the determination of the $3$-dimensional function $n(\vec{r})$, which minimizes the energy functional.
+However, the complexity associated with the many-electron problem is now relocated in the task of finding the well-defined but, in contrast to the potential energy, not explicitly known functional $F[n(\vec{r})]$.
 
-It is worth to note, that this minimal principle may be regarded as exactification of TF theory, which is rederived by the approximations
+It is worth to note, that this minimal principle may be regarded as exactification of the TF theory, which is rederived by the approximations
 \begin{equation}
-T=\int n(\vec{r})\frac{3}{10}k_{\text{F}}^2[n(\vec{r})]d\vec{r}
+T=\int n(\vec{r})\frac{3}{10}k_{\text{F}}^2(n(\vec{r}))d\vec{r}
 \text{ ,}
 \end{equation}
 \begin{equation}
@@ -345,12 +347,179 @@ U=\frac{1}{2}\int\frac{n(\vec{r})n(\vec{r}')}{|\vec{r}-\vec{r}'|}d\vec{r}d\vec{r
 
 \subsection{Kohn-Sham system}
 
-Now find $F[n]$ ...
+Inspired by the Hartree equations, i.e. a set of self-consistent single-particle equations for the approximate solution of the many-electron problem \cite{hartree28}, which describe atomic ground states much better than the TF theory, Kohn and Sham presented a Hartree-like formulation of the Hohenberg and Kohn minimal principle \eqref{eq:basics:hkm} \cite{kohn65}.
+However, due to a more general approach, the new formulation is formally exact by introducing the energy functional $E_{\text{xc}}[n(vec{r})]$, which accounts for the exchange and correlation energy of the electron interaction $U$ and possible corrections due to electron interaction to the kinetic energy $T$.
+The respective Kohn-Sham equations for the effective single-particle wave functions $\Phi_i(\vec{r})$ take the form
+\begin{equation}
+\left[
+ -\frac{\hbar^2}{2m}\nabla^2 + V_{\text{eff}}(\vec{r})
+\right] \Phi_i(\vec{r})=\epsilon_i\Phi_i(\vec{r})
+\label{eq:basics:kse1}
+\text{ ,}
+\end{equation}
+\begin{equation}
+V_{\text{eff}}(\vec{r})=V(\vec{r})+\int\frac{e^2n(\vec{r}')}{|\vec{r}-\vec{r}'|}d\vec{r}'
+ + V_{\text{xc}(\vec{r})}
+\text{ ,}
+\label{eq:basics:kse2}
+\end{equation}
+\begin{equation}
+n(\vec{r})=\sum_{i=1}^N |\Phi_i(\vec{r})|^2
+\text{ ,}
+\label{eq:basics:kse3}
+\end{equation}
+where the local exchange-correlation potential $V_{\text{xc}}(\vec{r})$ is the partial derivative of the exchange-correlation functional $E_{\text{xc}}[n(\vec{r})]$ with respect to the charge density $n(\vec{r})$ for the ground-state $n_0(\vec{r})$.
+The first term in equation \eqref{eq:basics:kse1} corresponds to the kinetic energy of non-interacting electrons and the second term of equation \eqref{eq:basics:kse2} is just the Hartree contribution $V_{\text{H}}(\vec{r})$ to the interaction energy.
+%\begin{equation}
+%V_{\text{xc}}(\vec{r})=\frac{\partial}{\partial n(\vec{r})}
+% E_{\text{xc}}[n(\vec{r})] |_{n(\vec{r})=n_0(\vec{r})}
+%\end{equation}
+
+The system of interacting electrons is mapped to an auxiliary system, the Kohn-Sham (KS) system, of non-interacting electrons in an effective potential.
+The exact effective potential $V_{\text{eff}}(\vec{r})$ may be regarded as a fictious external potential yielding a gound-state density for non-interacting electrons, which is equal to that for interacting electrons in the external potential $V(\vec{r})$.
+The one-electron KS orbitals $\Phi_i(\vec{r})$ as well as the respective KS energies $\epsilon_i$ are not directly attached to any physical observable except for the ground-state density, which is determined by equation \eqref{eq:basics:kse3} and the ionization energy, which is equal to the highest occupied relative to the vacuum level.
+The KS equations may be considered the formal exactification of the Hartree theory, which it is reduced to if the exchange-correlation potential and functional are neglected.
+In addition to the Hartree-Fock (HF) method, KS theory includes the difference of the kinetic energy of interacting and non-interacting electrons as well as the remaining contributions to the correlation energy that are not part of the HF correlation.
+
+The self-consistent KS equations \eqref{eq:basics:kse1}, \eqref{eq:basics:kse2} and \eqref{eq:basics:kse3} are non-linear partial differential equations, which may be solved numerically by an iterative process.
+Starting from a first approximation for $n(\vec{r})$ the effective potential $V_{\text{eff}}(\vec{r})$ can be constructed followed by determining the one-electron orbitals $\Phi_i(\vec{r})$, which solve the single-particle Schr\"odinger equation for the respective potential.
+The $\Phi_i(\vec{r})$ are used to obtain a new expression for $n(\vec{r})$.
+These steps are repeated until the initial and new density are equal or reasonably converged.
+
+Again, it is worth to note that the KS equations are formally exact.
+Assuming exact functionals $E_{\text{xc}}[n(\vec{r})]$ and potentials $V_{\text{xc}}(\vec{r})$ all many-body effects are included.
+Clearly, this directs attention to the functional, which now contains the costs involved with the many-electron problem.
 
 \subsection{Approximations for exchange and correlation}
+\label{subsection:ldagga}
+
+As discussed above, the HK and KS formulations are exact and so far no approximations except the adiabatic approximation entered the theory.
+However, to make concrete use of the theory, effective approximations for the exchange and correlation energy functional $E_{\text{xc}}[n(\vec{r})]$ are required.
+
+Most simple and at the same time remarkably useful is the approximation of $E_{\text{xc}}[n(\vec{r})]$ by a function of the local density \cite{kohn65}
+\begin{equation}
+E^{\text{LDA}}_{\text{xc}}[n(\vec{r})]=\int\epsilon_{\text{xc}}(n(\vec{r}))n(\vec{r}) d\vec{r}
+\text{ ,}
+\label{eq:basics:xca}
+\end{equation}
+which is, thus, called local density approximation (LDA).
+Here, the exchange-correlation energy per particle of the uniform electron gas of constant density $n$ is used for $\epsilon_{\text{xc}}(n(\vec{r}))$.
+Although, even in such a simple case, no exact form of the correlation part of $\epsilon_{\text{xc}}(n)$ is known, highly accurate numerical estimates using Monte Carlo methods \cite{ceperley80} and corresponding paramterizations exist \cite{perdew81}.
+Obviously exact for the homogeneous electron gas, the LDA was {\em a priori} expected to be useful only for densities varying slowly on scales of the local Fermi or TF wavelength.
+Nevertheless, LDA turned out to be extremely successful in describing some properties of highly inhomogeneous systems accurately within a few percent.
+Although LDA is known to overestimate the cohesive energy in solids by \unit[10-20]{\%}, the lattice parameters are typically determined with an astonishing accuracy of about \unit[1]{\%}.
+
+More accurate approximations of the exchange-correlation functional are realized by the introduction of gradient corrections with respect to the density \cite{kohn65}.
+Respective considerations are based on the concept of an exchange-correlation hole density describing the depletion of the electron density in the vicinity of an electron.
+The averaged hole density can be used to give a formally exact expression for $E_{\text{xc}}[n(\vec{r})]$ and an equivalent expression \cite{kohn96,kohn98}, which makes use of the electron density distribution $n(\tilde{\vec{r}})$ at positions $\tilde{\vec{r}}$ near $\vec{r}$, yielding the form
+\begin{equation}
+E_{\text{xc}}[n(\vec{r})]=\int\epsilon_{\text{xc}}(\vec{r};[n(\tilde{\vec{r}})])n(\vec{r}) d\vec{r}
+\end{equation}
+for the exchange-correlation energy, where $\epsilon_{\text{xc}}(\vec{r};[n(\tilde{\vec{r}})])$ becomes a nearsighted functional of $n(\tilde{\vec{r}})$.
+Expressing $n(\tilde{\vec{r}})$ in a Taylor series, $\epsilon_{\text{xc}}$ can be thought of as a function of coefficients, which correspond to the respective terms of the expansion.
+Neglecting all terms of order $\mathcal{O}(\nabla n(\vec{r})$ results in the functional equal to LDA, which requires the function of variable $n$.
+Including the next element of the Taylor series introduces the gradient correction to the functional, which requires the function of variables $n$ and $|\nabla n|$.
+This is called the generalized gradient approximation (GGA), which expresses the exchange-correlation energy density as a function of the local density and the local gradient of the density
+\begin{equation}
+E^{\text{GGA}}_{\text{xc}}[n(\vec{r})]=\int\epsilon_{\text{xc}}(n(\vec{r}),|\nabla n(\vec{r})|)n(\vec{r}) d\vec{r}
+\text{ .}
+\end{equation}
+These functionals constitute the simplest extensions of LDA for inhomogeneous systems.
+At modest computational costs gradient-corrected functionals very often yield much better results than the LDA with respect to cohesive and binding energies.
+
+\subsection{Plane-wave basis set}
+
+Finally, a set of basis functions is required to represent the one-electron KS wave functions.
+With respect to the numerical treatment it is favorable to approximate the wave functions by linear combinations of a finite number of such basis functions.
+Covergence of the basis set, i.e. convergence of the wave functions with respect to the amount of basis functions, is most crucial for the accuracy of the numerical calulations.
+Two classes of basis sets, the plane-wave and local basis sets, exist.
+
+Local basis set functions usually are atomic orbitals, i.e. mathematical functions that describe the wave-like behavior of electrons, which are localized, i.e. centered on atoms or bonds.
+Molecular orbitals can be represented by linear combinations of atomic orbitals (LCAO).
+By construction, only a small number of basis functions is required to represent all of the electrons of each atom within reasonable accuracy.
+Thus, local basis sets enable the implementation of methods that scale linearly with the number of atoms.
+However, these methods rely on the fact that the wave functions are localized and exhibit an exponential decay resulting in a sparse Hamiltonian.
+
+Another approach is to represent the KS wave functions by plane waves.
+In fact, the employed {\textsc vasp} software is solving the KS equations within a plane-wave (PW) basis set.
+The idea is based on the Bloch theorem \cite{bloch29}, which states that in a periodic crystal each electronic wave function $\Phi_i(\vec{r})$ can be written as the product of a wave-like envelope function $\exp(i\vec{kr})$ and a function that has the same periodicity as the lattice.
+The latter one can be expressed by a Fourier series, i.e. a discrete set of plane waves whose wave vectors just correspond to reciprocal lattice vectors $\vec{G}$ of the crystal.
+Thus, the one-electron wave function $\Phi_i(\vec{r})$ associated with the wave vector $\vec{k}$ can be expanded in terms of a discrete PW basis set
+\begin{equation}
+\Phi_i(\vec{r})=\sum_{\vec{G}
+%, |\vec{G}+\vec{k}|<G_{\text{cut}}}
+}c_{i,\vec{k}+\vec{G}} \exp\left(i(\vec{k}+\vec{G})\vec{r}\right)
+\text{ .}
+%E_{\text{cut}}=\frac{\hbar^2 G^2_{\text{cut}}}{2m}
+%\text{, }
+\end{equation}
+The basis set, which in principle should be infinite, can be truncated to include only plane waves that have kinetic energies $\hbar^2|\vec{k}+\vec{G}|^2/2m$ less than a particular cut-off energy $E_{\text{cut}}$.
+Although coefficients $c_{i,\vec{k}+\vec{G}}$ corresponding to small kinetic energies are typically more important, convergence with respect to the cut-off energy is crucial for the accuracy of the calculations.
+Convergence with respect to the basis set, however, is easily achieved by increasing $E_{\text{cut}}$ until the respective differences in total energy approximate zero.
+
+Next to their simplicity, plane waves have several advantages.
+The basis set is orthonormal by construction and, as mentioned above, it is simple to check for convergence.
+The biggest advantage, however, is the ability to perform exact calculations by a discrete sum over a numerical grid.
+This is due to the related construction of the grid and the PW basis.
+Ofcourse, exactness is restricted by the fact that the PW basis set is finite.
+The simple form of the PW representation of the KS equations
+\begin{equation}
+\sum_{\vec{G}'} \left[
+ \frac{\hbar^2}{2m}|\vec{k}+\vec{G}|^2 \delta_{\vec{GG}'}
+ + \tilde{V}(\vec{G}-\vec{G}')
+ + \tilde{V}_{\text{H}}(\vec{G}-\vec{G}')
+ + \tilde{V}_{\text{xc}}(\vec{G}-\vec{G}')
+\right] c_{i,\vec{k}+\vec{G}} = \epsilon_i c_{i,\vec{k}+\vec{G}}
+\label{eq:basics:pwks}
+\end{equation}
+reveals further advantages.
+The various potentials are described in terms of their Fourier transforms.
+Equation \eqref{eq:basics:pwks} is solved by diagonalization of the Hamiltonian matrix $H_{\vec{k}+\vec{G},\vec{k}+\vec{G}'}$ given by the terms in the brackets.
+The gradient operator is diagonal in reciprocal space whereas the exchange-correlation potential has a diagonal representation in real space.
+This suggests to carry out different operations in real and reciprocal space, which requires frequent Fourier transformations.
+These, however, can be efficiently achieved by the fast Fourier transformation (FFT) algorithm.
+
+There are likewise disadvantages associated with the PW representation.
+By construction, PW calculations require a periodic system.
+This does not pose a severe problem since non-periodic systems can still be described by a suitable choice of the simulation cell.
+Describing a defect, for instance, requires the inclusion of enough bulk material in the simulation to prevent or reduce the interaction with its periodic, artificial images.
+As a consequence the number of atoms involved in the calculations are increased.
+To describe surfaces, sufficiently thick vacuum layers need to be included to avoid interaction of adjacent crystal slabs.
+Clearly, to appropriately approximate the wave functions and the respective charge density of a system composed of vacuum in addition to the solid in a PW basis, an increase of the cut-off energy is required.
+According to equation \eqref{eq:basics:pwks} the size of the Hamiltonian depends on the cut-off energy and, therefore, the computational effort is likewise increased.
+For the same reason, the description of tightly bound core electrons and the respective, highly localized charge density is hindered.
+However, a much more profound problem exists whenever wave functions for the core as well as the valence electrons need to be calculated within a PW basis set.
+Wave functions of the valence electrons exhibit rapid oscillations in the region occupied by the core electrons near the nuclei.
+The oscillations maintain the orthogonality between the wave functions of the core and valence electrons, which is compulsory due to the exclusion principle.
+Accurately approximating these oscillations demands for an extremely large PW basis set, which is too large for practical use.
+Fortunately, the impossibility to model the core in addition to the valence electrons is eliminated in the pseudopotential approach discussed in the next section.
 
 \subsection{Pseudopotentials}
 
+As discussed in the last part of the previous section, an extremely large basis set of plane waves would be required to perform an all-electron calculation and a vast amount of computational time would be required to calculate the electronic wave functions.
+It is worth to stress out one more time, that this is due to the orthogonalization wiggles of the wave functions of valence electrons near the nuclei.
+Thus, existing core states practically prevent the use of a PW basis set.
+However, the core electrons, which are tightly bound to the nuclei, do not contribute significantly to chemical bonding or other physical properties of the solid.
+This fact is exploited in the pseudopotential approach \cite{} by removing the core electrons and replacing the atom and the associated strong ionic potential by a pseudoatom and a weaker pseudopotential that acts on a set of pseudo wave functions rather than the true valance wave functions.
+Certain conditions need to be fulfilled by the constructed pseudopotentials and the resulting pseudo wave functions.
+Outside the core region, the pseudo and real wafe functions as well as the generated charge densities need to be identical.
+...
+A pseudopotential is called norm-conserving if the pseudo and real charge contained within the core region match.
+...
+
+\subsection{Brillouin zone sampling}
+
+Following Bloch's theorem only a finite number of electronic wave functions need to be calculated for a periodic system.
+However, to calculate quantities like the total energy or charge density, these have to be evaluated in a sum over an infinite number of $\vec{k}$ points.
+Since the values of the wave function within a small interval around $\vec{k}$ are almost identical, it is possible to approximate the infinite sum by a sum over an affordable number of $k$ points, each representing the respective region of the wave function in $\vec{k}$ space. 
+Methods have been derived for obtaining very accurate approximations by a summation over special sets of $\vec{k}$ points with distinct, associated weights \cite{baldereschi73,chadi73,monkhorst76}.
+If present, symmetries in reciprocal space may further reduce the number of calculations.
+For supercells, i.e. repeating unit cells that contain several primitive cells, restricting the sampling of the Brillouin zone (BZ) to the $\Gamma$ point can yield quite accurat results.
+In fact, with respect to BZ sampling, calculating wave functions of a supercell containing $n$ primitive cells for only one $\vec{k}$ point is equivalent to the scenario of a single primitive cell and the summation over $n$ points in $\vec{k}$ space.
+In general, finer $\vec{k}$ point meshes better account for the periodicity of a system, which in some cases, however, might be fictious anyway.
+
+\subsection{Hellmann-Feynman forces}
+
 \section{Modeling of defects}
 \label{section:basics:defects}