X-Git-Url: https://hackdaworld.org/gitweb/?p=lectures%2Flatex.git;a=blobdiff_plain;f=posic%2Fthesis%2Fbasics.tex;h=ba62d82ca738adf57a3e6e1a7750cb8da6d57d94;hp=bd259d0653eb5aba4254d8b5899cc9e6f436be08;hb=9eca74b850bc1fa7c26fe65218b7d89d55c546c1;hpb=4686b1709127ca6616d82ca7ed97307f79598433

diff --git a/posic/thesis/basics.tex b/posic/thesis/basics.tex
index bd259d0..ba62d82 100644
--- a/posic/thesis/basics.tex
+++ b/posic/thesis/basics.tex
@@ -4,7 +4,7 @@
 In the following the simulation methods used within the scope of this study are introduced.
 Enabling the investigation of the evolution of structure on the atomic scale, molecular dynamics (MD) simulations are chosen for modeling the behavior and precipitation of C introduced into an initially crystalline Si environment.
 To be able to model systems with a large amount of atoms computational efficient classical potentials to describe the interaction of the atoms are most often used in MD studies.
-For reasons of flexibility in executing this non-standard task and in order to be able to use a novel interaction potential \cite{albe_sic_pot} an appropriate MD code called {\textsc posic}\footnote{{\textsc posic} is an abbreviation for {\bf p}recipitation {\bf o}f {\bf SiC}}\footnote{Source code: http://www.physik.uni-augsburg.de/\~{}zirkelfr/posic/posic.tar.bz2} including a library collecting respective MD subroutines was developed from scratch.
+For reasons of flexibility in executing this non-standard task and in order to be able to use a novel interaction potential \cite{albe_sic_pot} an appropriate MD code called {\textsc posic}\footnote{{\textsc posic} is an abbreviation for {\bf p}recipitation {\bf o}f {\bf SiC}} including a library collecting respective MD subroutines was developed from scratch\footnote{Source code: http://www.physik.uni-augsburg.de/\~{}zirkelfr/posic}.
 The basic ideas of MD in general and the adopted techniques as implemented in {\textsc posic} in particular are outlined in section \ref{section:md}, while the functional form and derivative of the employed classical potential is presented in appendix \ref{app:d_tersoff}.
 An overview of the most important tools within the MD package is given in appendix \ref{app:code}.
 Although classical potentials are often most successful and at the same time computationally efficient in calculating some physical properties of a particular system, not all of its properties might be described correctly due to the lack of quantum-mechanical effects.
@@ -19,8 +19,7 @@ The method used to investigate migration pathways to identify the prevalent diff
 \section{Molecular dynamics simulations}
 \label{section:md}
 
-% todo
-% rewrite!
+% todo - rewrite md intro chapter
 
 \begin{quotation}
 \dq We may regard the present state of the universe as the effect of the past and the cause of the future. An intellect which at any given moment knew all of the forces that animate nature and the mutual positions of the beings that compose it, if this intellect were vast enough to submit the data to analysis, could condense into a single formula the movement of the greatest bodies of the universe and that of the lightest atom; for such an intellect nothing could be uncertain and the future just like the past would be present before its eyes.\dq{}
@@ -32,7 +31,7 @@ The method used to investigate migration pathways to identify the prevalent diff
 \noindent
 Pierre Simon de Laplace phrased this vision in terms of a controlling, omniscient instance - the {\em Laplace demon} - which would be able to look into the future as well as into the past due to the deterministic nature of processes, governed by the solution of differential equations.
 Although Laplace's vision is nowadays corrected by chaos theory and quantum mechanics, it expresses two main features of classical mechanics, the determinism of processes and time reversibility of the fundamental equations.
-This understanding was one of the first ideas for doing molecular dynamics simulations, considering an isolated system of particles, the behaviour of which is fully determined by the solution of the classical equations of motion.
+This understanding may be regarded as the basic principle of molecular dynamics, considering an isolated system of particles, the behaviour of which is fully determined by the solution of the classical equations of motion.
 
 \subsection{Introduction to molecular dynamics simulations}
 
@@ -40,8 +39,8 @@ Molecular dynamics (MD) simulation is a technique to compute a system of particl
 The MD method was first introduced by Alder and Wainwright in 1957 \cite{alder57,alder59} to study the interactions of hard spheres.
 The basis of the approach are Newton's equations of motion to describe classicaly the many-body system.
 MD is the numerical way of solving the $N$-body problem which cannot be solved analytically for $N>3$.
-A potential is necessary describing the interaction of the particles.
-By MD a complete description of the system in the sense of classical mechanics on the microscopic level is obtained.
+A potential is necessary to describe the interaction of the particles.
+By MD, a complete description of the system in the sense of classical mechanics on the microscopic level is obtained.
 The microscopic information can then be translated to macroscopic observables by means of statistical mechanics.
 
 The basic idea is to assume that the particles can be described classically by Newton's equations of motion, which are integrated numerically.
@@ -50,15 +49,17 @@ A system of $N$ particles of masses $m_i$ ($i=1,\ldots,N$) at positions ${\bf r}
 %m_i \frac{d^2}{dt^2} {\bf r}_i = {\bf F}_i \Leftrightarrow
 %m_i \frac{d}{dt} {\bf r}_i = {\bf p}_i\textrm{ , } \quad
 %\frac{d}{dt} {\bf p}_i = {\bf F}_i\textrm{ .}
-m_i \ddot{{\bf r}_i} = {\bf F}_i \Leftrightarrow
-m_i \dot{{\bf r}_i} = {\bf p}_i\textrm{, }
-\dot{{\bf p}_i} = {\bf F}_i\textrm{ .}
+m_i \ddot{\bf r}_i = {\bf F}_i \Leftrightarrow
+m_i \dot{\bf r}_i = {\bf p}_i\textrm{, }
+\dot{\bf p}_i = {\bf F}_i\textrm{ .}
+\label{eq:basics:newton}
 \end{equation}
 The forces ${\bf F}_i$ are obtained from the potential energy $U(\{{\bf r}\})$:
 \begin{equation}
 {\bf F}_i = - \nabla_{{\bf r}_i} U({\{\bf r}\}) \, \textrm{.}
+\label{eq:basics:force}
 \end{equation}
-Given the initial conditions ${\bf r}_i(t_0)$ and $\dot{{\bf r}}_i(t_0)$ the equations can be integrated by a certain integration algorithm.
+Given the initial conditions ${\bf r}_i(t_0)$ and $\dot{\bf r}_i(t_0)$, the equations can be integrated by a certain integration algorithm.
 The solution of these equations provides the complete information of a system evolving in time.
 The following sections cover the tools of the trade necessary for the MD simulation technique.
 Three ingredients are required for a MD simulation:
@@ -120,7 +121,7 @@ The attractive part is associated with the bonding.
 f_R(r_{ij}) & = & A_{ij} \exp (- \lambda_{ij} r_{ij} ) \\
 f_A(r_{ij}) & = & -B_{ij} \exp (- \mu_{ij} r_{ij} )
 \end{eqnarray}
-The function $f_C$ is the a cutoff function to limit the range of interaction to nearest neighbors.
+The function $f_C$ is a cutoff function to limit the range of interaction to nearest neighbors.
 It is designed to have a smooth transition of the potential at distances $R_{ij}$ and $S_{ij}$.
 \begin{equation}
 f_C(r_{ij}) = \left\{
@@ -139,7 +140,7 @@ b_{ij} & = & \chi_{ij} (1 + \beta_i^{n_i} \zeta^{n_i}_{ij})^{-1/2n_i} \\
 g(\theta_{ijk}) & = & 1 + c_i^2/d_i^2 - c_i^2/[d_i^2 + (h_i - \cos \theta_{ijk})^2]
 \end{eqnarray}
 where $\theta_{ijk}$ is the bond angle between bonds $ij$ and $ik$.
-This is illustrated in Figure \ref{img:tersoff_angle}.
+This is illustrated in Fig. \ref{img:tersoff_angle}.
 \begin{figure}[t]
 \begin{center}
 \includegraphics[width=8cm]{tersoff_angle.eps}
@@ -182,12 +183,12 @@ The potential succeeds in the description of the low as well as high coordinated
 The description of elastic properties of SiC is improved with respect to the potentials available in literature.
 Defect properties are only fairly reproduced but the description is comparable to previously published potentials.
 It is claimed that the potential enables modeling of widely different configurations and transitions among these and has recently been used to simulate the inert gas condensation of Si-C nanoparticles \cite{erhart04}.
-Therefore the Erhart/Albe (EA) potential is considered the superior analytical bond order potential to study the SiC precipitation and associated processes in Si.
+Therefore, the Erhart/Albe (EA) potential is considered the superior analytical bond order potential to study the SiC precipitation and associated processes in Si.
 
 \subsection{Verlet integration}
 \label{subsection:integrate_algo}
 
-A numerical method to integrate Newton's equation of motion was presented by Verlet in 1967 \cite{verlet67}.
+A numerical method to integrate Newton's equations of motion was presented by Verlet in 1967 \cite{verlet67}.
 The idea of the so-called Verlet and a variant, the velocity Verlet algorithm, which additionaly generates directly the velocities, is explained in the following.
 Starting point is the Taylor series for the particle positions at time $t+\delta t$ and $t-\delta t$
 \begin{equation}
@@ -358,7 +359,7 @@ The respective Kohn-Sham equations for the effective single-particle wave functi
 \text{ ,}
 \end{equation}
 \begin{equation}
-V_{\text{eff}}=V(\vec{r})+\int\frac{e^2n(\vec{r}')}{|\vec{r}-\vec{r}'|}d\vec{r}'
+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}
@@ -369,7 +370,7 @@ n(\vec{r})=\sum_{i=1}^N |\Phi_i(\vec{r})|^2
 \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 to the interaction energy.
+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})}
@@ -417,9 +418,9 @@ E_{\text{xc}}[n(\vec{r})]=\int\epsilon_{\text{xc}}(\vec{r};[n(\tilde{\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$.
+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
+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{ .}
@@ -441,10 +442,10 @@ Thus, local basis sets enable the implementation of methods that scale linearly
 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 basis set.
+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 plane-wave basis set
+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}}}
@@ -456,32 +457,219 @@ Thus, the one-electron wave function $\Phi_i(\vec{r})$ associated with the wave
 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.
-matrix elements of the Hamiltonian have a simple form (pw rep of ks equations)
-As mentioned above ... simple to check for convergence.
 
-Disadvantage ... periodic system required, but escapable by respective choice of the supercell.
-size of matrix to diagonalize determined by cut-off energy, severe 
+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}
 
-Since core electrons tend to be concentrated very close to the atomic nuclei, resulting in large wavefunction and density gradients near the nuclei which are not easily described by a plane-wave basis set unless a very high energy cutoff, and therefore small wavelength, is used.
+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 mainly 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 (PP) approach \cite{cohen70} by removing the core electrons and replacing the atom and the associated strong ionic potential by a pseudoatom and a weaker PP that acts on a set of pseudo wave functions rather than the true valance wave functions.
+This way, the pseudo wave functions become smooth near the nuclei.
+
+Most PPs statisfy four general conditions.
+The pseudo wave functions generated by the PP should not contain nodes, i.e. the pseudo wave functions should be smooth and free of wiggles in the core region.
+Outside the core region, the pseudo and real valence wave functions as well as the generated charge densities need to be identical.
+The charge enclosed within the core region must be equal for both wave functions.
+Last, almost redundantly, the valence all-electron and pseudopotential eigenvalues must be equal.
+Pseudopotentials that meet the conditions outlined above are referred to as norm-conserving pseudopotentials \cite{hamann79}.
+
+%Certain properties need to be fulfilled by PPs and the resulting pseudo wave functions.
+%The pseudo wave functions should yield the same energy eigenvalues than the true valence wave functions.
+%The PP is called norm-conserving if the pseudo and real charge contained within the core region matches.
+%To guarantee transferability of the PP the logarithmic derivatives of the real and pseudo wave functions and their first energy derivatives need to agree outside of the core region.
+%A simple procedure was proposed to extract norm-conserving PPs obyeing the above-mentioned conditions from {\em ab initio} atomic calculations \cite{hamann79}.
+
+In order to achieve these properties different PPs are required for the different shapes of the orbitals, which are determined by their angular momentum.
+Mathematically, a non-local PP, which depends on the angular momentum, has the form
+\begin{equation}
+V_{\text{nl}}(\vec{r}) = \sum_{lm} \mid lm \rangle V_l(\vec{r}) \langle lm \mid
+\text{ .}
+\end{equation}
+Applying of the operator $V_{\text{nl}}(\vec{r})$ decomposes the electronic wave functions into spherical harmonics $\mid lm \rangle$, i.e. the orbitals with azimuthal angular momentum $l$ and magnetic number $m$, which are then multiplied by the respective psuedopotential $V_l(\vec{r})$ for angular momentum $l$.
+The standard generation procedure of pseudopotentials proceeds by varying its parameters until the pseudo eigenvalues are eual to the all-electron valence eigenvalues and the pseudo wave functions match the all-electron valence wave functions beyond a certain cut-off radius detrmining the core region.
+Modified methods to generate ultra-soft pseudopotentials were proposed, which address the rapid convergence with respect to the size of the plane wave basis set \cite{vanderbilt90,troullier91}.
+
+Using PPs the rapid oscillations of the wave functions near the core of the atoms are removed considerably reducing the number of plane waves necessary to appropriately expand the wave functions.
+More importantly, less accuracy is required compared to all-electron calculations to determine energy differences among ionic configurations, which almost totally appear in the energy of the valence electrons that are typically a factor $10^3$ smaller than the energy of the core electrons.
 
 \subsection{Brillouin zone sampling}
+\label{subsection:basics:bzs}
 
-Due to the Bloch theorem only a finite number of electronic wave functions need to be calculated for a periodic system.
+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 an intergration over special sets of $\vec{k}$ points \cite{}.
+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, sampling of the Brillouin zone restricted to the $\Gamma$ point can be quite accurately used, which is equivalent to calculating a single primitive cell using multiple $\vec{k}$ points.
+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}
+\subsection{Structural relaxation and Hellmann-Feynman theorem}
 
-\section{Modeling of defects}
+Up to this point, the system is in the ground state with respect to the electronic subsystem, while the positions of the ions as well as size and shape of the supercell are fixed.
+To investigate equilibrium structures, however, the ionic subsystem must also be allowed to relax into a minimum energy configuration.
+Local minimum configurations can be easily obtained in a MD-like way by moving the nuclei over small distances along the directions of the forces, as discussed in the MD chapter above.
+Clearly, the conjugate gradient method constitutes a more sophisticated scheme, which will locate the equilibrium positions of the ions more rapidly.
+To find the global minimum, i.e. the absolute ground state, methods like simulated annealing or the Monte Carlo technique, which allow the system to escape local minima, have to be used for the search.
+
+The force on an ion is given by the negative derivative of the total energy with respect to the position of the ion.
+However, moving an ion, i.e. altering its position, changes the wave functions to the KS eigenstates corresponding to the new ionic configuration.
+Writing down the derivative of the total energy $E$ with respect to the position $\vec{R}_i$ of ion $i$
+\begin{equation}
+\frac{dE}{d\vec{R_i}}=
+ \sum_j \Phi_j^* \frac{\partial H}{\partial \vec{R}_i} \Phi_j
++\sum_j \frac{\partial \Phi_j^*}{\partial \vec{R}_i} H \Phi_j
++\sum_j \Phi_j^* H \frac{\partial \Phi_j}{\partial \vec{R}_i}
+\text{ ,}
+\end{equation}
+indeed reveals a contributon to the chnage in total energy due to the change of the wave functions $\Phi_j$.
+However, provided that the $\Phi_j$ are eigenstates of $H$, it is easy to show that the last two terms cancel each other and in the special case of $H=T+V$ the force is given by
+\begin{equation}
+\vec{F}_i=-\sum_j \Phi_j^*\Phi_j\frac{\partial V}{\partial \vec{R}_i}
+\text{ .}
+\end{equation}
+This is called the Hellmann-Feynman theorem \cite{feynman39}, which enables the calculation of forces, called the Hellmann-Feynman forces, acting on the nuclei for a given configuration, without the need for evaluating computationally costly energy maps.
+
+\section{Modeling of point defects}
 \label{section:basics:defects}
 
+Point defects are defects that affect a single lattice site.
+At this site the crystalline periodicity is interrupted.
+An empty lattice site, which would be occupied in the perfect crystal structure, is called a vacancy defect.
+If an additional atom is incorporated into the perfect crystal, this is called interstitial defect.
+A substitutional defect exists, if an atom belonging to the perfect crystal is replaced with an atom of another species.
+The disturbance caused by these defects may result in the distortion of the surrounding atomic structure and is accompanied by an increase in configurational energy.
+Thus, next to the structure of the defect, the energy needed to create such a defect, i.e. the defect formation energy, is an important value characterizing the defect and likewise determining its relative stability.
+
+The formation energy of a defect is defined by
+\begin{equation}
+E_{\text{f}}=E-\sum_i N_i\mu_i
+\text{ ,}
+\label{eq:basics:ef2}
+\end{equation}
+where $E$ is the total energy of the interstitial structure involving $N_i$ atoms of type $i$ with chemical potential $\mu_i$.
+Here, the chemical potentials are determined by the chemical potential of the respective equilibrium bulk structure, i.e. the cohesive energy per atom for the fully relaxed structure at zero temperature and pressure.
+Considering C interstitial defects in Si, the chemical potential for C could also be determined by the cohesive energies of Si and SiC according to the relation $\mu_{\text{C}}=\mu_{\text{SiC}}-\mu_{\text{Si}}$ of the chemical potentials.
+In this way, SiC is chosen as a reservoir for the C impurity.
+For defect configurations consisting of a single atom species the formation energy reduces to
+\begin{equation}
+E_{\text{f}}=\left(E_{\text{coh}}^{\text{defect}}
+                  -E_{\text{coh}}^{\text{defect-free}}\right)N
+\text{ ,}
+\label{eq:basics:ef1}
+\end{equation}
+where $N$ and $E_{\text{coh}}^{\text{defect}}$ are the number of atoms and the cohesive energy per atom in the defect configuration and $E_{\text{coh}}^{\text{defect-free}}$ is the cohesive energy per atom of the defect-free structure.
+Clearly, for a single atom species equation \eqref{eq:basics:ef2} is equivalent to equation \eqref{eq:basics:ef1} since $NE_{\text{coh}}^{\text{defect}}$ is equal to the total energy of the defect structure and $NE_{\text{coh}}^{\text{defect-free}}$ corresponds to $N\mu$, provided the structure is fully relaxed at zero temperature.
+
+However, there is hardly ever only one defect in a crystal, not even only one kind of defect.
+Again, energetic considerations can be used to investigate the existing interaction of two defects.
+The binding energy $E_{\text{b}}$ of a defect pair is given by the difference of the formation energy of the defect combination $E_{\text{f}}^{\text{comb}} $ and the sum of the two separated defect configurations $E_{\text{f}}^{1^{\text{st}}}$ and $E_{\text{f}}^{2^{\text{nd}}}$.
+This can be expressed by
+\begin{equation}
+E_{\text{b}}=
+E_{\text{f}}^{\text{comb}}-
+E_{\text{f}}^{1^{\text{st}}}-
+E_{\text{f}}^{2^{\text{nd}}}
+\label{eq:basics:e_bind}
+\end{equation}
+where the formation energies $E_{\text{f}}^{\text{comb}}$, $E_{\text{f}}^{1^{\text{st}}}$ and $E_{\text{f}}^{2^{\text{nd}}}$ are determined as discussed above.
+Accordingly, energetically favorable configurations result in binding energies below zero while unfavorable configurations show positive values for the binding energy.
+The interaction strength, i.e. the absolute value of the binding energy, approaches zero for increasingly non-interacting isolated defects.
+Thus, $E_{\text{b}}$ indeed can be best thought of a binding energy, which is required to bring the defects to infinite separation.
+
+The methods presented in the last two chapters can be used to investigate defect structures and energetics.
+Therefore, a supercell containing the perfect crystal is generated in an initial process.
+If not by construction, the system should be fully relaxed.
+The substitutional or vacancy defect is realized by replacing or removing one atom contained in the supercell.
+Interstitial defects are created by adding an atom at positions located in the space between regular lattice sites.
+Once the intuitively created defect structure is generated structural relaxation methods will yield the respective local minimum configuration.
+Since the supercell approach applies periodic boundary conditions enough bulk material surrounding the defect is required to exclude possible interaction of the defect with its periodic image.
+
+\begin{figure}[t]
+\begin{center}
+\includegraphics[width=9cm]{unit_cell_e.eps}
+\end{center}
+\caption[Insertion positions for interstitial defect atoms in the diamond lattice.]{Insertion positions for the tetrahedral ({\color{red}$\bullet$}), hexagonal  ({\color{green}$\bullet$}), \hkl<1 0 0> dumbbell ({\color{yellow}$\bullet$}), \hkl<1 1 0> dumbbell ({\color{magenta}$\bullet$}) and bond-centered ({\color{cyan}$\bullet$}) interstitial defect atom in the diamond lattice. The black dots correspond to the lattice atoms and the blue lines indicate the covalent bonds of the perfect diamond structure.}
+\label{fig:basics:ins_pos}
+\end{figure}
+The respective estimated interstitial insertion positions for various interstitial structures in a diamond lattice are displayed in Fig. \ref{fig:basics:ins_pos}.
+The labels of the interstitial types indicate their positions in the interstitial lattice.
+In a dumbbell (DB) configuration two atoms share a single lattice site along a certain direction that is also comprehended in the label of the defect.
+For the DB configurations the nearest atom of the bulk lattice is slightly displaced by $(0,0,-1/8)$ and $(-1/8,-1/8,0)$ of the unit cell length respectively.
+This is indicated by the dashed, unfilled circles.
+By this, high forces, which might enable the system to overcome barriers of the local minimum configuration and, thus, result in a different structure, are avoided.
+
 \section{Migration paths and diffusion barriers}
 \label{section:basics:migration}
 
+Investigating diffusion mechanisms is based on determining migration paths in between two local minimum configurations of an atom at different locations in the lattice.
+During migration, the total energy of the system increases, traverses at least one maximum of the configurational energy and finally decreases to a local minimum value.
+The maximum difference in energy is the barrier necessary for the respective migration process.
+The path exhibiting the minimal energy difference determines the diffusion path and associated diffusion barrier and the maximum configuration turns into a saddle point configuration.
+
+\begin{figure}[t]
+\begin{center}
+\subfigure[]{\label{fig:basics:crto}\includegraphics[width=0.45\textwidth]{crt_orig.eps}}
+\subfigure[]{\label{fig:basics:crtm}\includegraphics[width=0.45\textwidth]{crt_mod.eps}}
+\end{center}
+\caption{Schematic of the constrained relaxation technique (a) and of a modified version (b) used to obtain migration pathways and corresponding configurational energies.}
+\label{fig:basics:crt}
+\end{figure}
+One possibility to compute the migration path from one stable cofiguration into another one is provided by the constrained relaxation technique (CRT) \cite{kaukonen98}.
+The atoms involving great structural changes in the diffusion process are moved stepwise from the starting to the final position and relaxation after each step is only allowed in the plane perpendicular to the direction of the vector connecting its starting and final position.
+This is illustrated in Fig. \ref{fig:basics:crto}.
+The number of steps required for smooth transitions depends on the shape of the potential energy surface.
+No constraints are applied to the remaining atoms to allow for the relaxation of the surrounding lattice.
+To prevent the remaining lattice to shift according to the displacement of the defect, ohowever, some atoms far away from the defect region should be fixed in all three coordinate directions.
+However, for the present study, the method tremendously failed.
+Abrupt changes in structure and configurational energy occured among relaxed structures of two successive displacement steps.
+For some structures even the expected final configurations are not obtained.
+Thus, the method mentioned above is adjusted adding further constraints in order to obtain smooth transitions with repsect to energy and structure.
+In the modified method all atoms are stepwise displaced towards their final positions.
+In addition to this, relaxation of each individual atom is only allowed in the plane perpendicular to the last individual displacement vector, as displayed in Fig. \ref{fig:basics:crtm}.
+In the modified version respective energies could be higher than the real ones due to the additional constraints that prevent a more adequate relaxation until the final copnfiguration is reached.
+
+Structures of maximum configurational energy do not necessarily constitute saddle point configurations, i .e. the method does not guarantee to find the true minimum energy path.
+Whether a saddle point configuration and, thus, the minimum energy path is obtained by the CRT method, needs to be verified by caculating the respective vibrational modes.
+Modifications used to add the CRT feature to the VASP code and a short instruction on how to use it can be found in appendix \ref{app:patch_vasp}.
+
+% todo - advantages of pw basis concenring hf forces
+