1 \chapter{Basic principles of utilized simulation techniques}
4 In the following the simulation methods used within the scope of this study are introduced.
5 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.
6 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.
7 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/download/posic/posic.tar.bz2}.
8 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}.
9 An overview of the most important tools within the MD package is given in appendix \ref{app:code}.
10 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.
11 Thus, in order to obtain more accurate results quantum-mechanical calculations from first principles based on density functional theory (DFT) were performed.
12 The Vienna {\em ab initio} simulation package ({\textsc vasp}) \cite{kresse96} is used for this purpose.
13 The relevant basics of DFT are described in section \ref{section:dft} while an overview of utilities mainly used to create input or parse output data of {\textsc vasp} is given in appendix \ref{app:code}.
14 The gain in accuracy achieved by this method, however, is accompanied by an increase in computational effort constraining the simulated system to be much smaller in size.
15 Thus, investigations based on DFT are restricted to single defects or combinations of two defects in a rather small Si supercell, their structural relaxation as well as some selected diffusion processes.
16 Next to the structure, defects can be characterized by the defect formation energy, a scalar indicating the costs necessary for the formation of the defect, which is explained in section \ref{section:basics:defects}.
17 The method used to investigate migration pathways to identify the prevalent diffusion mechanism is introduced in section \ref{section:basics:migration} and modifications to the {\textsc vasp} code implementing this method are presented in appendix \ref{app:patch_vasp}.
19 \section{Molecular dynamics simulations}
22 % todo - rewrite md intro chapter
25 \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{}
27 {\em Marquis Pierre Simon de Laplace, 1814.} \cite{laplace}
32 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.
33 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.
34 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.
36 \subsection{Introduction to molecular dynamics simulations}
38 Molecular dynamics (MD) simulation is a technique to compute a system of particles, referred to as molecules, with their positions, volocities and forces among each other evolving in time.
39 The MD method was first introduced by Alder and Wainwright in 1957 \cite{alder57,alder59} to study the interactions of hard spheres.
40 The basis of the approach are Newton's equations of motion to describe classicaly the many-body system.
41 MD is the numerical way of solving the $N$-body problem which cannot be solved analytically for $N>3$.
42 A potential is necessary describing the interaction of the particles.
43 By MD a complete description of the system in the sense of classical mechanics on the microscopic level is obtained.
44 The microscopic information can then be translated to macroscopic observables by means of statistical mechanics.
46 The basic idea is to assume that the particles can be described classically by Newton's equations of motion, which are integrated numerically.
47 A system of $N$ particles of masses $m_i$ ($i=1,\ldots,N$) at positions ${\bf r}_i$ and velocities $\dot{{\bf r}}_i$ is given by
49 %m_i \frac{d^2}{dt^2} {\bf r}_i = {\bf F}_i \Leftrightarrow
50 %m_i \frac{d}{dt} {\bf r}_i = {\bf p}_i\textrm{ , } \quad
51 %\frac{d}{dt} {\bf p}_i = {\bf F}_i\textrm{ .}
52 m_i \ddot{\bf r}_i = {\bf F}_i \Leftrightarrow
53 m_i \dot{\bf r}_i = {\bf p}_i\textrm{, }
54 \dot{\bf p}_i = {\bf F}_i\textrm{ .}
55 \label{eq:basics:newton}
57 The forces ${\bf F}_i$ are obtained from the potential energy $U(\{{\bf r}\})$:
59 {\bf F}_i = - \nabla_{{\bf r}_i} U({\{\bf r}\}) \, \textrm{.}
60 \label{eq:basics:force}
62 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.
63 The solution of these equations provides the complete information of a system evolving in time.
64 The following sections cover the tools of the trade necessary for the MD simulation technique.
65 Three ingredients are required for a MD simulation:
67 \item A model for the interaction between system constituents is needed.
68 Interaction potentials and their accuracy for describing certain systems of elements will be outlined in section \ref{subsection:interact_pot}.
69 \item An integrator is needed, which propagtes the particle positions and velocities from time $t$ to $t+\delta t$, realised by a finite difference scheme which moves trajectories discretely in time.
70 This is explained in section \ref{subsection:integrate_algo}.
71 \item A statistical ensemble has to be chosen, which allows certain thermodynamic quantities to be controlled or to stay constant.
72 This is discussed in section \ref{subsection:statistical_ensembles}.
74 These ingredients will be outlined in the follwoing.
75 The discussion is restricted to methods employed within this study.
77 \subsection{Interaction potentials for silicon and carbon}
78 \label{subsection:interact_pot}
80 The potential energy of $N$ interacting atoms can be written in the form
82 U(\{{\bf r}\}) = \sum_i U_1({\bf r}_i) + \sum_i \sum_{j>i} U_2({\bf r}_i,{\bf r}_j) + \sum_i \sum_{j>i} \sum_{k>j>i} U_3({\bf r}_i,{\bf r}_j,{\bf r}_k) \ldots
84 where $U$ is the total potential energy.
85 $U_1$ is a single particle potential describing external forces.
86 Examples of single particle potentials are the gravitational force or an electric field.
87 $U_2$ is a two body pair potential which only depends on the distance $r_{ij}$ between the two atoms $i$ and $j$.
88 If not only pair potentials are considered, three body potentials $U_3$ or multi body potentials $U_n$ can be included.
89 Usually these higher order terms are avoided since they are not easy to model and it is rather time consuming to evaluate potentials and forces originating from these many body terms.
90 Ordinary pair potentials have a close-packed structure like face-centered cubic (FCC) or hexagonal close-packed (HCP) as a ground state.
91 A pair potential is, thus, unable to describe properly elements with other structures than FCC or HCP.
92 Silicon and carbon for instance, have a diamand and zincblende structure with four covalently bonded neighbors, which is far from a close-packed structure.
93 A three body potential has to be included for these types of elements.
95 \subsubsection{The Tersoff bond order potential}
97 Tersoff proposed an empirical interatomic potential for covalent systems \cite{tersoff_si1,tersoff_si2}.
98 The Tersoff potential explicitly incorporates the dependence of bond order on local environments, permitting an improved description of covalent materials.
99 Due to the covalent character Tersoff restricted the interaction to nearest neighbor atoms accompanied by an increases in computational efficiency for the evaluation of forces and energy based on the short-range potential.
100 Tersoff applied the potential to silicon \cite{tersoff_si1,tersoff_si2,tersoff_si3}, carbon \cite{tersoff_c} and also to multicomponent systems like silicon carbide \cite{tersoff_m}.
101 The basic idea is that, in real systems, the bond order, i.e. the strength of the bond, depends upon the local environment \cite{abell85}.
102 Atoms with many neighbors form weaker bonds than atoms with only a few neighbors.
103 Although the bond strength intricately depends on geometry the focus on coordination, i.e. the number of neighbors forming bonds, is well motivated qualitatively from basic chemistry since for every additional formed bond the amount of electron pairs per bond and, thus, the strength of the bonds is decreased.
104 If the energy per bond decreases rapidly enough with increasing coordination the most stable structure will be the dimer.
105 In the other extreme, if the dependence is weak, the material system will end up in a close-packed structure in order to maximize the number of bonds and likewise minimize the cohesive energy.
106 This suggests the bond order to be a monotonously decreasing function with respect to coordination and the equilibrium coordination being determined by the balance of bond strength and number of bonds.
107 Based on pseudopotential theory the bond order term $b_{ijk}$ limitting the attractive pair interaction is of the form $b_{ijk}\propto Z^{-\delta}$ where $Z$ is the coordination number and $\delta$ a constant \cite{abell85}, which is $\frac{1}{2}$ in the seond-moment approximation within the tight binding scheme \cite{horsfield96}.
109 Tersoff incorporated the concept of bond order in a three-body potential formalism.
110 The interatomic potential is taken to have the form
112 E & = & \sum_i E_i = \frac{1}{2} \sum_{i \ne j} V_{ij} \textrm{ ,} \\
113 V_{ij} & = & f_C(r_{ij}) [ f_R(r_{ij}) + b_{ij} f_A(r_{ij}) ] \textrm{ .}
115 $E$ is the total energy of the system, constituted either by the sum over the site energies $E_i$ or by the bond energies $V_{ij}$.
116 The indices $i$ and $j$ correspond to the atoms of the system with $r_{ij}$ being the distance from atom $i$ to atom $j$.
117 The functions $f_R$ and $f_A$ represent a repulsive and an attractive pair potential.
118 The repulsive part is due to the orthogonalization energy of overlapped atomic wave functions.
119 The attractive part is associated with the bonding.
121 f_R(r_{ij}) & = & A_{ij} \exp (- \lambda_{ij} r_{ij} ) \\
122 f_A(r_{ij}) & = & -B_{ij} \exp (- \mu_{ij} r_{ij} )
124 The function $f_C$ is the a cutoff function to limit the range of interaction to nearest neighbors.
125 It is designed to have a smooth transition of the potential at distances $R_{ij}$ and $S_{ij}$.
127 f_C(r_{ij}) = \left\{
129 1, & r_{ij} < R_{ij} \\
130 \frac{1}{2} + \frac{1}{2} \cos \Big[ \pi (r_{ij} - R_{ij})/(S_{ij} - R_{ij}) \Big], & R_{ij} < r_{ij} < S_{ij} \\
135 As discussed above, $b_{ij}$ represents a measure of the bond order, monotonously decreasing with the coordination of atoms $i$ and $j$.
138 b_{ij} & = & \chi_{ij} (1 + \beta_i^{n_i} \zeta^{n_i}_{ij})^{-1/2n_i} \\
139 \zeta_{ij} & = & \sum_{k \ne i,j} f_C (r_{ik}) \omega_{ik} g(\theta_{ijk}) \\
140 g(\theta_{ijk}) & = & 1 + c_i^2/d_i^2 - c_i^2/[d_i^2 + (h_i - \cos \theta_{ijk})^2]
142 where $\theta_{ijk}$ is the bond angle between bonds $ij$ and $ik$.
143 This is illustrated in Fig. \ref{img:tersoff_angle}.
146 \includegraphics[width=8cm]{tersoff_angle.eps}
148 \caption{Angle between bonds of atoms $i,j$ and $i,k$.}
149 \label{img:tersoff_angle}
151 The angular dependence does not give a fixed minimum angle between bonds since the expression is embedded inside the bond order term.
152 The relation to the above discussed bond order potential becomes obvious if $\chi=1, \beta=1, n=1, \omega=1$ and $c=0$.
153 Parameters with a single subscript correspond to the parameters of the elemental system \cite{tersoff_si3,tersoff_c} while the mixed parameters are obtained by interpolation from the elemental parameters by the arithmetic or geometric mean.
154 The elemental parameters were obtained by fit with respect to the cohesive energies of real and hypothetical bulk structures and the bulk modulus and bond length of the diamond structure.
155 New parameters for the mixed system are $\chi$, which is used to finetune the strength of heteropolar bonds, and $\omega$, which is set to one for the C-Si interaction but is available as a feature to permit the application of the potential to more drastically different types of atoms in the future.
157 The force acting on atom $i$ is given by the derivative of the potential energy.
158 For a three body potential ($V_{ij} \neq V{ji}$) the derivation is of the form
160 \nabla_{{\bf r}_i} E = \frac{1}{2} \big[ \sum_j ( \nabla_{{\bf r}_i} V_{ij} + \nabla_{{\bf r}_i} V_{ji} ) + \sum_k \sum_j \nabla_{{\bf r}_i} V_{jk} \big] \textrm{ .}
162 The force is then given by
164 F^i = - \nabla_{{\bf r}_i} E \textrm{ .}
166 Details of the Tersoff potential derivative are presented in appendix \ref{app:d_tersoff}.
168 \subsubsection{Improved analytical bond order potential}
170 Although the Tersoff potential is one of the most widely used potentials there are some shortcomings.
171 Describing the Si-Si interaction Tersoff was unable to find a single parameter set to describe well both, bulk and surface properties.
172 Due to this and since the first approach labeled T1 \cite{tersoff_si1} turned out to be unstable \cite{dodson87}, two further parametrizations exist, T2 \cite{tersoff_si2} and T3 \cite{tersoff_si3}.
173 While T2 describes well surface properties, T3 yields improved elastic constants and should be used for describing bulk properties.
174 However, T3, which is used in the Si/C potential, suffers from an underestimation of the dimer binding energy.
175 Similar behavior is found for the C-C interaction.
177 For this reason, Erhart and Albe provide a reparametrization of the Tersoff potential based on three independently fitted parameter sets for the Si-Si, C-C and Si-C interaction \cite{albe_sic_pot}.
178 The functional form is similar to the one proposed by Tersoff.
179 Differences in the energy functional and the force evaluation routine are pointed out in appendix \ref{app:d_tersoff}.
180 Concerning Si the elastic properties of the diamond phase as well as the structure and energetics of the dimer are reproduced very well.
181 The new parameter set for the C-C interaction yields improved dimer properties while at the same time delivers a description of the bulk phase similar to the Tersoff potential.
182 The potential succeeds in the description of the low as well as high coordinated structures.
183 The description of elastic properties of SiC is improved with respect to the potentials available in literature.
184 Defect properties are only fairly reproduced but the description is comparable to previously published potentials.
185 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}.
186 Therefore, the Erhart/Albe (EA) potential is considered the superior analytical bond order potential to study the SiC precipitation and associated processes in Si.
188 \subsection{Verlet integration}
189 \label{subsection:integrate_algo}
191 A numerical method to integrate Newton's equation of motion was presented by Verlet in 1967 \cite{verlet67}.
192 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.
193 Starting point is the Taylor series for the particle positions at time $t+\delta t$ and $t-\delta t$
195 \vec{r}_i(t+\delta t)=
196 \vec{r}_i(t)+\delta t\vec{v}_i(t)+\frac{\delta t^2}{2m_i}\vec{f}_i(t)+
197 \frac{\delta t^3}{6}\vec{b}_i(t) + \mathcal{O}(\delta t^4)
198 \label{basics:verlet:taylor1}
201 \vec{r}_i(t-\delta t)=
202 \vec{r}_i(t)-\delta t\vec{v}_i(t)+\frac{\delta t^2}{2m_i}\vec{f}_i(t)-
203 \frac{\delta t^3}{6}\vec{b}_i(t) + \mathcal{O}(\delta t^4)
204 \label{basics:verlet:taylor2}
206 where $\vec{v}_i=\frac{d}{dt}\vec{r}_i$ are the velocities, $\vec{f}_i=m\frac{d}{dt^2}\vec{r}_i$ are the forces and $\vec{b}_i=\frac{d}{dt^3}\vec{r}_i$ are the third derivatives of the particle positions with respect to time.
207 The Verlet algorithm is obtained by summarizing and substracting equations \eqref{basics:verlet:taylor1} and \eqref{basics:verlet:taylor2}
209 \vec{r}_i(t+\delta t)=
210 2\vec{r}_i(t)-\vec{r}_i(t-\delta t)+\frac{\delta t^2}{m_i}\vec{f}_i(t)+
211 \mathcal{O}(\delta t^4)
214 \vec{v}_i(t)=\frac{1}{2\delta t}[\vec{r}_i(t+\delta t)-\vec{r}_i(t-\delta t)]+
215 \mathcal{O}(\delta t^3)
217 the truncation error of which is of order $\delta t^4$ for the positions and $\delta t^3$ for the velocities.
218 The velocities, although not used to update the particle positions, are not synchronously determined with the positions but drag behind one step of discretization.
219 The Verlet algorithm can be rewritten into an equivalent form, which updates the velocities and positions in the same step.
220 The so-called velocity Verlet algorithm is obtained by combining \eqref{basics:verlet:taylor1} with equation \eqref{basics:verlet:taylor2} displaced in time by $+\delta t$
222 \vec{v}_i(t+\delta t)=
223 \vec{v}_i(t)+\frac{\delta t}{2m_i}[\vec{f}_i(t)+\vec{f}_i(t+\delta t)]
226 \vec{r}_i(t+\delta t)=
227 \vec{r}_i(t)+\delta t\vec{v}_i(t)+\frac{\delta t^2}{2m_i}\vec{f}_i(t) \text{ .}
229 Since the forces for the new positions are required to update the velocity the determination of the forces has to be carried out within the integration algorithm.
231 \subsection{Statistical ensembles}
232 \label{subsection:statistical_ensembles}
234 Using the above mentioned algorithms the most basic type of MD is realized by simply integrating the equations of motion of a fixed number of particles ($N$) in a closed volume $V$ realized by periodic boundary conditions (PBC).
235 Providing a stable integration algorithm the total energy $E$, i.e. the kinetic and configurational energy of the paticles, is conserved.
236 This is known as the $NVE$, or microcanonical ensemble, describing an isolated system composed of microstates, among which the number of particles, volume and energy are held constant.
238 However, the successful formation of SiC dictates precise control of temperature by external heating.
239 While the temperature of such a system is well defined, the energy is no longer conserved.
240 The microscopic states of a system, which is in thermal equilibrium with an external thermal heat bath, are represented by the $NVT$ ensemble.
241 In the so-called canonical ensemble the temperature $T$ is related to the expactation value of the kinetic energy of the particles, i.e.
243 \langle E_{\text{kin}}\rangle = \frac{3}{2}Nk_{\text{B}}T \text{, }
244 E_{\text{kin}}=\sum_i \frac{\vec{p}^2_i}{2m_i} \text{ .}
247 The volume of the synthesized material can hardly be controlled in experiment.
248 Instead the pressure can be adjusted.
249 Holding constant the pressure in addition to the temperature of the system its states are represented by the isothermal-isobaric $NpT$ ensemble.
250 The expression for the pressure of a system derived from the equipartition theorem is given by
252 pV=Nk_{\text{B}}T+\langle W\rangle\text{, }W=-\frac{1}{3}\sum_i\vec{r}_i\nabla_{\vec{r}_i}U
256 where $W$ is the virial and $U$ is the configurational energy.
258 Berendsen~et~al.~\cite{berendsen84} proposed a method, which is easy to implement, to couple the system to an external bath with constant temperature $T_0$ or pressure $p_0$ with adjustable time constants $\tau_T$ and $\tau_p$ determining the strength of the coupling.
259 Control of the respective variable is based on the relations given in equations \eqref{eq:basics:ts} and \eqref{eq:basics:ps}.
260 The thermostat is achieved by scaling the velocities of all atoms in every time step $\delta t$ from $\vec{v}_i$ to $\lambda \vec{v}_i$, with
262 \lambda=\left[1+\frac{\delta t}{\tau_T}(\frac{T_0}{T}-1)\right]^\frac{1}{2}
265 where $T$ is the current temperature according to equation \eqref{eq:basics:ts}.
266 The barostat adjusts the pressure by changing the virial through scaling of the particle positions $\vec{r}_i$ to $\mu \vec{r}_i$ and the volume $V$ to $\mu^3 V$, with
268 \mu=\left[1-\frac{\beta\delta t}{\tau_p}(p_0-p)\right]^\frac{1}{3}\text{ ,}
270 where $\beta$ is the isothermal compressibility and $p$ corresponds to the current pressure, which is determined by equation \eqref{eq:basics:ps}.
272 Using this method the system does not behave like a true $NpT$ ensemble.
273 On average $T$ and $p$ correspond to the expected values.
274 For large enough time constants, i.e. $\tau > 100 \delta t$, the method shows realistic fluctuations in $T$ and $p$.
275 The advantage of the approach is that the coupling can be decreased to minimize the disturbance of the system and likewise be adjusted to suit the needs of a given application.
276 It provides a stable algorithm that allows smooth changes of the system to new values of temperature or pressure, which is ideal for the investigated problem.
278 \section{Denstiy functional theory}
281 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.
282 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.
283 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.
284 This cannot be solved exactly and finding approximate solutions requires several layers of simplification in order to reduce the number of free parameters.
285 Approximations that consider a truncated Hilbert space of single-particle orbitals yield promising results, however, with increasing complexity and demand for high accuracy the amount of Slater determinats to be evaluated massively increases.
287 In contrast, instead of using the description by the many-body wave function, the key point in density functional theory (DFT) is to recast the problem to a description utilizing the charge density $n(\vec{r})$, which constitutes a quantity in real space depending only on the three spatial coordinates.
288 In the following sections the basic idea of DFT will be outlined.
289 As will be shown, DFT can formally be regarded as an exactification of the Thomas Fermi theory \cite{thomas27,fermi27} and the self-consistent Hartree equations \cite{hartree28}.
290 A nice review is given in the Nobel lecture of Kohn \cite{kohn99}, one of the inventors of DFT.
292 \subsection{Born-Oppenheimer approximation}
294 Born and Oppenheimer proposed a simplification enabling the effective decoupling of the electronic and ionic degrees of freedom \cite{born27}.
295 Within the Born-Oppenheimer (BO) approximation the light electrons are assumed to move much faster and, thus, follow adiabatically to the motion of the heavy nuclei, if the latter are only slightly deflected from their equilibrium positions.
296 Thus, on the timescale of electronic motion the ions appear at fixed positions.
297 On the other way round, on the timescale of nuclear motion the electrons appear blurred in space adding an extra term to the ion-ion potential.
298 The simplified Schr\"odinger equation no longer contains the kinetic energy of the ions.
299 The momentary positions of the ions enter as fixed parameters and, therefore, the ion-ion interaction may be regarded as a constant added to the electronic energies.
300 The Schr\"odinger equation describing the remaining electronic problem reads
302 \left[-\frac{\hbar^2}{2m}\sum_j\nabla^2_j-
303 \sum_{j,l} \frac{Z_le^2}{|\vec{r}_j-\vec{R}_l|}+
304 \frac{1}{2}\sum_{j\neq j'}\frac{e^2}{|\vec{r}_j-\vec{r}_{j'}|}
305 \right] \Psi = E \Psi
308 where $Z_l$ are the atomic numbers of the nuclei and $\Psi$ is the many-electron wave function, which depends on the positions and spins of the electrons.
309 Accordingly, there is only a parametrical dependence on the ionic coordinates $\vec{R}_l$.
310 However, the remaining number of free parameters is still too high and need to be further decreased.
312 \subsection{Hohenberg-Kohn theorem and variational principle}
314 Investigating the energetics of Cu$_x$Zn$_{1-x}$ alloys, which for different compositions exhibit different transfers of charge between the Cu and Zn unit cells due to their chemical difference and, thus, varying electrostatic interactions contributing to the total energy, the work of Hohenberg and Kohn had a natural focus on the distribution of charge.
315 Although it was clear that the Thomas Fermi (TF) theory only provides a rough approximation to the exact solution of the many-electron Schr\"odinger equation the theory was of high interest since it provides an implicit relation of the potential and the electron density distribution.
316 This raised the question how to establish a connection between TF expressed in terms of $n(\vec{r})$ and the exact Schr\"odinger equation expressed in terms of the many-electron wave function $\Psi({\vec{r}})$ and whether a complete description in terms of the charge density is possible in principle.
317 The answer to this question, whether the charge density completely characterizes a system, became the starting point of modern DFT.
319 Considering a system with a nondegenerate ground state there is obviously only one ground-state charge density $n_0(\vec{r})$ that correpsonds to a given potential $V(\vec{r})$.
320 In 1964 Hohenberg and Kohn showed the opposite and far less obvious result \cite{hohenberg64}.
321 Employing no more than the Rayleigh-Ritz minimal principle it is concluded by {\em reductio ad absurdum} that for a nondegenerate ground state the same charge density cannot be generated by different potentials.
322 Thus, the charge density of the ground state $n_0(\vec{r})$ uniquely determines the potential $V(\vec{r})$ and, consequently, the full Hamiltonian and ground-state energy $E_0$.
323 In mathematical terms the full many-electron ground state is a unique functional of the charge density.
324 Im particular, $E$ is a functional $E[n(\vec{r})]$ of $n(\vec{r})$.
326 The ground-state charge density $n_0(\vec{r})$ minimizes the energy functional $E[n(\vec{r})]$, its value corresponding to the ground-state energy $E_0$, which enables the formulation of a minimal principle in terms of the charge density \cite{hohenberg64,levy82}
328 E_0=\min_{n(\vec{r})}
330 F[n(\vec{r})] + \int n(\vec{r}) V(\vec{r}) d\vec{r}
333 \label{eq:basics:hkm}
335 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})]$.
336 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.
337 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})]$.
339 It is worth to note, that this minimal principle may be regarded as exactification of the TF theory, which is rederived by the approximations
341 T=\int n(\vec{r})\frac{3}{10}k_{\text{F}}^2(n(\vec{r}))d\vec{r}
345 U=\frac{1}{2}\int\frac{n(\vec{r})n(\vec{r}')}{|\vec{r}-\vec{r}'|}d\vec{r}d\vec{r}'
349 \subsection{Kohn-Sham system}
351 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}.
352 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$.
353 The respective Kohn-Sham equations for the effective single-particle wave functions $\Phi_i(\vec{r})$ take the form
356 -\frac{\hbar^2}{2m}\nabla^2 + V_{\text{eff}}(\vec{r})
357 \right] \Phi_i(\vec{r})=\epsilon_i\Phi_i(\vec{r})
358 \label{eq:basics:kse1}
362 V_{\text{eff}}(\vec{r})=V(\vec{r})+\int\frac{e^2n(\vec{r}')}{|\vec{r}-\vec{r}'|}d\vec{r}'
363 + V_{\text{xc}(\vec{r})}
365 \label{eq:basics:kse2}
368 n(\vec{r})=\sum_{i=1}^N |\Phi_i(\vec{r})|^2
370 \label{eq:basics:kse3}
372 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})$.
373 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.
375 %V_{\text{xc}}(\vec{r})=\frac{\partial}{\partial n(\vec{r})}
376 % E_{\text{xc}}[n(\vec{r})] |_{n(\vec{r})=n_0(\vec{r})}
379 The system of interacting electrons is mapped to an auxiliary system, the Kohn-Sham (KS) system, of non-interacting electrons in an effective potential.
380 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})$.
381 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.
382 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.
383 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.
385 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.
386 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.
387 The $\Phi_i(\vec{r})$ are used to obtain a new expression for $n(\vec{r})$.
388 These steps are repeated until the initial and new density are equal or reasonably converged.
390 Again, it is worth to note that the KS equations are formally exact.
391 Assuming exact functionals $E_{\text{xc}}[n(\vec{r})]$ and potentials $V_{\text{xc}}(\vec{r})$ all many-body effects are included.
392 Clearly, this directs attention to the functional, which now contains the costs involved with the many-electron problem.
394 \subsection{Approximations for exchange and correlation}
395 \label{subsection:ldagga}
397 As discussed above, the HK and KS formulations are exact and so far no approximations except the adiabatic approximation entered the theory.
398 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.
400 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}
402 E^{\text{LDA}}_{\text{xc}}[n(\vec{r})]=\int\epsilon_{\text{xc}}(n(\vec{r}))n(\vec{r}) d\vec{r}
404 \label{eq:basics:xca}
406 which is, thus, called local density approximation (LDA).
407 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}))$.
408 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}.
409 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.
410 Nevertheless, LDA turned out to be extremely successful in describing some properties of highly inhomogeneous systems accurately within a few percent.
411 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]{\%}.
413 More accurate approximations of the exchange-correlation functional are realized by the introduction of gradient corrections with respect to the density \cite{kohn65}.
414 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.
415 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
417 E_{\text{xc}}[n(\vec{r})]=\int\epsilon_{\text{xc}}(\vec{r};[n(\tilde{\vec{r}})])n(\vec{r}) d\vec{r}
419 for the exchange-correlation energy, where $\epsilon_{\text{xc}}(\vec{r};[n(\tilde{\vec{r}})])$ becomes a nearsighted functional of $n(\tilde{\vec{r}})$.
420 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.
421 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$.
422 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|$.
423 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
425 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}
428 These functionals constitute the simplest extensions of LDA for inhomogeneous systems.
429 At modest computational costs gradient-corrected functionals very often yield much better results than the LDA with respect to cohesive and binding energies.
431 \subsection{Plane-wave basis set}
433 Finally, a set of basis functions is required to represent the one-electron KS wave functions.
434 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.
435 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.
436 Two classes of basis sets, the plane-wave and local basis sets, exist.
438 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.
439 Molecular orbitals can be represented by linear combinations of atomic orbitals (LCAO).
440 By construction, only a small number of basis functions is required to represent all of the electrons of each atom within reasonable accuracy.
441 Thus, local basis sets enable the implementation of methods that scale linearly with the number of atoms.
442 However, these methods rely on the fact that the wave functions are localized and exhibit an exponential decay resulting in a sparse Hamiltonian.
444 Another approach is to represent the KS wave functions by plane waves.
445 In fact, the employed {\textsc vasp} software is solving the KS equations within a plane-wave (PW) basis set.
446 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.
447 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.
448 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
450 \Phi_i(\vec{r})=\sum_{\vec{G}
451 %, |\vec{G}+\vec{k}|<G_{\text{cut}}}
452 }c_{i,\vec{k}+\vec{G}} \exp\left(i(\vec{k}+\vec{G})\vec{r}\right)
454 %E_{\text{cut}}=\frac{\hbar^2 G^2_{\text{cut}}}{2m}
457 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}}$.
458 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.
459 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.
461 Next to their simplicity, plane waves have several advantages.
462 The basis set is orthonormal by construction and, as mentioned above, it is simple to check for convergence.
463 The biggest advantage, however, is the ability to perform exact calculations by a discrete sum over a numerical grid.
464 This is due to the related construction of the grid and the PW basis.
465 Ofcourse, exactness is restricted by the fact that the PW basis set is finite.
466 The simple form of the PW representation of the KS equations
468 \sum_{\vec{G}'} \left[
469 \frac{\hbar^2}{2m}|\vec{k}+\vec{G}|^2 \delta_{\vec{GG}'}
470 + \tilde{V}(\vec{G}-\vec{G}')
471 + \tilde{V}_{\text{H}}(\vec{G}-\vec{G}')
472 + \tilde{V}_{\text{xc}}(\vec{G}-\vec{G}')
473 \right] c_{i,\vec{k}+\vec{G}} = \epsilon_i c_{i,\vec{k}+\vec{G}}
474 \label{eq:basics:pwks}
476 reveals further advantages.
477 The various potentials are described in terms of their Fourier transforms.
478 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.
479 The gradient operator is diagonal in reciprocal space whereas the exchange-correlation potential has a diagonal representation in real space.
480 This suggests to carry out different operations in real and reciprocal space, which requires frequent Fourier transformations.
481 These, however, can be efficiently achieved by the fast Fourier transformation (FFT) algorithm.
483 There are likewise disadvantages associated with the PW representation.
484 By construction, PW calculations require a periodic system.
485 This does not pose a severe problem since non-periodic systems can still be described by a suitable choice of the simulation cell.
486 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.
487 As a consequence the number of atoms involved in the calculations are increased.
488 To describe surfaces, sufficiently thick vacuum layers need to be included to avoid interaction of adjacent crystal slabs.
489 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.
490 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.
491 For the same reason, the description of tightly bound core electrons and the respective, highly localized charge density is hindered.
492 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.
493 Wave functions of the valence electrons exhibit rapid oscillations in the region occupied by the core electrons near the nuclei.
494 The oscillations maintain the orthogonality between the wave functions of the core and valence electrons, which is compulsory due to the exclusion principle.
495 Accurately approximating these oscillations demands for an extremely large PW basis set, which is too large for practical use.
496 Fortunately, the impossibility to model the core in addition to the valence electrons is eliminated in the pseudopotential approach discussed in the next section.
498 \subsection{Pseudopotentials}
500 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.
501 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.
502 Thus, existing core states practically prevent the use of a PW basis set.
503 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.
504 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.
505 This way, the pseudo wave functions become smooth near the nuclei.
507 Most PPs statisfy four general conditions.
508 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.
509 Outside the core region, the pseudo and real valence wave functions as well as the generated charge densities need to be identical.
510 The charge enclosed within the core region must be equal for both wave functions.
511 Last, almost redundantly, the valence all-electron and pseudopotential eigenvalues must be equal.
512 Pseudopotentials that meet the conditions outlined above are referred to as norm-conserving pseudopotentials \cite{hamann79}.
514 %Certain properties need to be fulfilled by PPs and the resulting pseudo wave functions.
515 %The pseudo wave functions should yield the same energy eigenvalues than the true valence wave functions.
516 %The PP is called norm-conserving if the pseudo and real charge contained within the core region matches.
517 %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.
518 %A simple procedure was proposed to extract norm-conserving PPs obyeing the above-mentioned conditions from {\em ab initio} atomic calculations \cite{hamann79}.
520 In order to achieve these properties different PPs are required for the different shapes of the orbitals, which are determined by their angular momentum.
521 Mathematically, a non-local PP, which depends on the angular momentum, has the form
523 V_{\text{nl}}(\vec{r}) = \sum_{lm} \mid lm \rangle V_l(\vec{r}) \langle lm \mid
526 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$.
527 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.
528 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}.
530 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.
531 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.
533 \subsection{Brillouin zone sampling}
534 \label{subsection:basics:bzs}
536 Following Bloch's theorem only a finite number of electronic wave functions need to be calculated for a periodic system.
537 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.
538 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.
539 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}.
540 If present, symmetries in reciprocal space may further reduce the number of calculations.
541 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.
542 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.
543 In general, finer $\vec{k}$ point meshes better account for the periodicity of a system, which in some cases, however, might be fictious anyway.
545 \subsection{Structural relaxation and Hellmann-Feynman theorem}
547 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.
548 To investigate equilibrium structures, however, the ionic subsystem must also be allowed to relax into a minimum energy configuration.
549 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.
550 Clearly, the conjugate gradient method constitutes a more sophisticated scheme, which will locate the equilibrium positions of the ions more rapidly.
551 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.
553 The force on an ion is given by the negative derivative of the total energy with respect to the position of the ion.
554 However, moving an ion, i.e. altering its position, changes the wave functions to the KS eigenstates corresponding to the new ionic configuration.
555 Writing down the derivative of the total energy $E$ with respect to the position $\vec{R}_i$ of ion $i$
557 \frac{dE}{d\vec{R_i}}=
558 \sum_j \Phi_j^* \frac{\partial H}{\partial \vec{R}_i} \Phi_j
559 +\sum_j \frac{\partial \Phi_j^*}{\partial \vec{R}_i} H \Phi_j
560 +\sum_j \Phi_j^* H \frac{\partial \Phi_j}{\partial \vec{R}_i}
563 indeed reveals a contributon to the chnage in total energy due to the change of the wave functions $\Phi_j$.
564 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
566 \vec{F}_i=-\sum_j \Phi_j^*\Phi_j\frac{\partial V}{\partial \vec{R}_i}
569 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.
571 \section{Modeling of point defects}
572 \label{section:basics:defects}
574 Point defects are defects that affect a single lattice site.
575 At this site the crystalline periodicity is interrupted.
576 An empty lattice site, which would be occupied in the perfect crystal structure, is called a vacancy defect.
577 If an additional atom is incorporated into the perfect crystal, this is called interstitial defect.
578 A substitutional defect exists, if an atom belonging to the perfect crystal is replaced with an atom of another species.
579 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.
580 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.
582 The formation energy of a defect is defined by
584 E_{\text{f}}=E-\sum_i N_i\mu_i
586 \label{eq:basics:ef2}
588 where $E$ is the total energy of the interstitial structure involving $N_i$ atoms of type $i$ with chemical potential $\mu_i$.
589 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.
590 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.
591 In this way, SiC is chosen as a reservoir for the C impurity.
592 For defect configurations consisting of a single atom species the formation energy reduces to
594 E_{\text{f}}=\left(E_{\text{coh}}^{\text{defect}}
595 -E_{\text{coh}}^{\text{defect-free}}\right)N
597 \label{eq:basics:ef1}
599 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.
600 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.
602 However, there is hardly ever only one defect in a crystal, not even only one kind of defect.
603 Again, energetic considerations can be used to investigate the existing interaction of two defects.
604 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}}}$.
605 This can be expressed by
608 E_{\text{f}}^{\text{comb}}-
609 E_{\text{f}}^{1^{\text{st}}}-
610 E_{\text{f}}^{2^{\text{nd}}}
611 \label{eq:basics:e_bind}
613 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.
614 Accordingly, energetically favorable configurations result in binding energies below zero while unfavorable configurations show positive values for the binding energy.
615 The interaction strength, i.e. the absolute value of the binding energy, approaches zero for increasingly non-interacting isolated defects.
616 Thus, $E_{\text{b}}$ indeed can be best thought of a binding energy, which is required to bring the defects to infinite separation.
618 The methods presented in the last two chapters can be used to investigate defect structures and energetics.
619 Therefore, a supercell containing the perfect crystal is generated in an initial process.
620 If not by construction, the system should be fully relaxed.
621 The substitutional or vacancy defect is realized by replacing or removing one atom contained in the supercell.
622 Interstitial defects are created by adding an atom at positions located in the space between regular lattice sites.
623 Once the intuitively created defect structure is generated structural relaxation methods will yield the respective local minimum configuration.
624 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.
628 \includegraphics[width=9cm]{unit_cell_e.eps}
630 \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.}
631 \label{fig:basics:ins_pos}
633 The respective estimated interstitial insertion positions for various interstitial structures in a diamond lattice are displayed in Fig. \ref{fig:basics:ins_pos}.
634 The labels of the interstitial types indicate their positions in the interstitial lattice.
635 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.
636 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.
637 This is indicated by the dashed, unfilled circles.
638 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.
640 \section{Migration paths and diffusion barriers}
641 \label{section:basics:migration}
643 Investigating diffusion mechanisms is based on determining migration paths in between two local minimum configurations of an atom at different locations in the lattice.
644 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.
645 The maximum difference in energy is the barrier necessary for the respective migration process.
646 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.
650 \subfigure[]{\label{fig:basics:crto}\includegraphics[width=0.45\textwidth]{crt_orig.eps}}
651 \subfigure[]{\label{fig:basics:crtm}\includegraphics[width=0.45\textwidth]{crt_mod.eps}}
653 \caption{Schematic of the constrained relaxation technique (a) and of a modified version (b) used to obtain migration pathways and corresponding configurational energies.}
654 \label{fig:basics:crt}
656 One possibility to compute the migration path from one stable cofiguration into another one is provided by the constrained relaxation technique (CRT) \cite{kaukonen98}.
657 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.
658 This is illustrated in Fig. \ref{fig:basics:crto}.
659 The number of steps required for smooth transitions depends on the shape of the potential energy surface.
660 No constraints are applied to the remaining atoms to allow for the relaxation of the surrounding lattice.
661 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.
662 However, for the present study, the method tremendously failed.
663 Abrupt changes in structure and configurational energy occured among relaxed structures of two successive displacement steps.
664 For some structures even the expected final configurations are not obtained.
665 Thus, the method mentioned above is adjusted adding further constraints in order to obtain smooth transitions with repsect to energy and structure.
666 In the modified method all atoms are stepwise displaced towards their final positions.
667 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}.
668 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.
670 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.
671 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.
672 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}.
674 % todo - advantages of pw basis concenring hf forces + inc font in crt sketch