\chapter{Basic principles of utilized simulation techniques}
+\label{chapter:basics}
+
+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.
+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.
+Thus, in order to obtain more accurate results quantum-mechanical calculations from first principles based on density functional theory (DFT) were performed.
+The Vienna {\em ab initio} simulation package ({\textsc vasp}) \cite{kresse96} is used for this purpose.
+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}.
+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.
+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.
+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}.
+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}.
+
+\section{Molecular dynamics simulations}
+\label{section:md}
+
+% todo
+% rewrite!
\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{}
-\cite{laplace}
\begin{flushright}
-{\em Marquis Pierre Simon de Laplace, 1814.}
+{\em Marquis Pierre Simon de Laplace, 1814.} \cite{laplace}
\end{flushright}
\end{quotation}
-\section{Molecular dynamics simulations}
-
+\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.
\subsection{Introduction to molecular dynamics simulations}
-Basically, 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.
-The MD method was first introduced by Alder and Wainwright in 1957 \cite{alder1,alder2} to study the interactions of hard spheres.
+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.
+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 simulation is the numerical way of solving the $N$-body problem which cannot be solved analytically ($N>3$).
-Quantum mechanical effects are taken into account by an analytical interaction potential between the nuclei.
-
-By MD simulation techniques a complete description of the system in the sense of classical mechanics on the microscopic level is obtained.
-This microscopic information has to be translated to macroscopic observables by means of statistical mechanics.
+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.
+The microscopic information can then be translated to macroscopic observables by means of statistical mechanics.
-The basic idea is to integrate Newton's equations numerically.
+The basic idea is to assume that the particles can be described classically by Newton's equations of motion, which are integrated numerically.
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
\begin{equation}
-m_i \frac{d^2}{dt^2} {\bf r}_i = {\bf F}_i \, \textrm{.}
+%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{ .}
\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{.}
\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.
The solution of these equations provides the complete information of a system evolving in time.
-
-The following chapters cover the tools of the trade necessary for the MD simulation technique.
-Three ingredients are required for an MD simulation:
+The following sections cover the tools of the trade necessary for the MD simulation technique.
+Three ingredients are required for a MD simulation:
\begin{enumerate}
\item A model for the interaction between system constituents is needed.
- Interaction potentials and their accuracy for describing certain systems of elements will be outlined in chapter \ref{subsection:interact_pot}.
+ Interaction potentials and their accuracy for describing certain systems of elements will be outlined in section \ref{subsection:interact_pot}.
\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.
- In chapter \ref{subsection:integrate_algo} a detailed overview of the available integration algorithms is given, including their advantages and disadvantages.
+ This is explained in section \ref{subsection:integrate_algo}.
\item A statistical ensemble has to be chosen, which allows certain thermodynamic quantities to be controlled or to stay constant.
- This is discussed in chapter \ref{subsection:statistical_ensembles}
+ This is discussed in section \ref{subsection:statistical_ensembles}.
\end{enumerate}
+These ingredients will be outlined in the follwoing.
+The discussion is restricted to methods employed within this study.
-In addition special techniques will be outlined which reduce the complexity of the MD algorithm, though the force/energy evaluation almost inevitably dictates the overall speed.
-
-\subsection{Integration algorithms}
-\label{subsection:integrate_algo}
-
-\subsection{Interaction potentials}
+\subsection{Interaction potentials for silicon and carbon}
\label{subsection:interact_pot}
The potential energy of $N$ interacting atoms can be written in the form
Examples of single particle potentials are the gravitational force or an electric field.
$U_2$ is a two body pair potential which only depends on the distance $r_{ij}$ between the two atoms $i$ and $j$.
If not only pair potentials are considered, three body potentials $U_3$ or multi body potentials $U_n$ can be included.
-Mainly 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.
-
+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.
Ordinary pair potentials have a close-packed structure like face-centered cubic (FCC) or hexagonal close-packed (HCP) as a ground state.
-A pair potential is thus unable to describe properly elements with other structures than FCC or HCP.
-Silicon and carbon for instance, have a diamand/zincblende structure with four covalent bonded neighbours, which is far from a close-packed structure.
+A pair potential is, thus, unable to describe properly elements with other structures than FCC or HCP.
+Silicon and carbon for instance, have a diamand and zincblende structure with four covalently bonded neighbors, which is far from a close-packed structure.
A three body potential has to be included for these types of elements.
-In the following, relevant potentials for this work are discussed.
-
-%\subsubsection{The Lennard-Jones potential}
-%
-%The L-J potential is a realistic two body pair potential and is of the form
-%\begin{equation}
-%U^{LJ}(r) = 4 \epsilon \Big[ \Big( \frac{\sigma}{r} \Big)^{12} - \Big( \frac{\sigma}{r} \Big)^6 \Big] \, \textrm{,}
-%\label{eq:lj-p}
-%\end{equation}
-%where $r$ denotes the distance between the two atoms.
-%
-%The attractive tail for large separations $(\sim r^{-6})$ is essentially due to correlations between electron clouds surrounding the atoms. The attractive part is also known as {\em van der Waals} or {\em London} interaction.
-%It can be derived classically by considering how two charged spheres induce dipol-dipol interactions into each other, or by considering the interaction between two oscillators in a quantum mechanical way.
-%
-%The repulsive term $(\sim r^{-12})$ captures the non-bonded overlap of the electron clouds.
-%It does not have a true physical motivation, other than the exponent being larger than $6$ to get a steep rising repulsive potential wall at short distances.
-%Chosing $12$ as the exponent of the repulsive term it is just the square of the attractive term which makes the potential evaluable in a very efficient way.
-%
-%The constants $\epsilon$ and $\sigma$ are usually determined by fitting to experimental data.
-%$\epsilon$ accounts to the depth of the potential well, where $\sigma$ is regarded as the radius of the particle, also known as the van der Waals radius.
-%
-%Writing down the derivative of the Lennard-Jones potential in respect to $x_i$ (the $i$th component of the distance vector ${\bf r}$)
-%\begin{equation}
-%\frac{\partial}{\partial x_i} U^{LJ}(r) = 4 \epsilon x_i \Big( -12 \frac{\sigma^{12}}{r^{14}} + 6 \frac{\sigma^6}{r^8} \Big)
-%\label{eq:lj-d}
-%\end{equation}
-%one can easily identify $\sigma$ by the equilibrium distance of the atoms $r_e=\sqrt[6]{2} \sigma$.
-%Applying the equilibrium distance into \eqref{eq:lj-p} $\epsilon$ turns out to be the negative well depth.
-%The $i$th component of the force is given by
-%\begin{equation}
-%F_i = - \frac{\partial}{\partial x_i} U^{LJ}(r) \, \textrm{.}
-%\label{eq:lj-f}
-%\end{equation}
-%
-%\subsubsection{The Stillinger Weber potential}
-%
-%The Stillinger Weber potential \cite{stillinger_weber} \ldots
-%
-%\begin{equation}
-%U = \sum_{i,j} U_2({\bf r}_i,{\bf r}_j) + \sum_{i,j,k} U_3({\bf r}_i,{\bf r}_j,{\bf r}_k)
-%\end{equation}
-%
-%\begin{equation}
-%U_2(r_{ij}) = \left\{
-% \begin{array}{ll}
-% \epsilon A \Big( B (r_{ij} / \sigma)^{-p} - 1\Big) \exp \Big[ (r_{ij} / \sigma - 1)^{-1} \Big] & r_{ij} < a \sigma \\
-% 0 & r_{ij} \ge a \sigma
-% \end{array} \right.
-%\end{equation}
-%
-%\begin{equation}
-%U_3({\bf r}_i,{\bf r}_j,{\bf r}_k) =
-%\epsilon \Big[ h(r_{ij},r_{ik},\theta_{jik}) + h(r_{ji},r_{jk},\theta_{ijk}) + h(r_{ki},r_{kj},\theta_{ikj}) \Big]
-%\end{equation}
-%
-%\begin{equation}
-%h(r_{ij},r_{ik},\theta_{jik}) =
-%\lambda \exp \Big[ \gamma (r_{ij}/\sigma -a)^{-1} + \gamma (r_{ik}/\sigma - a)^{-1} \Big] \Big( \cos \theta_{jik} + \frac{1}{3} \Big)^2
-%\end{equation}
-
-\subsubsection{The Tersoff potential}
-
-Tersoff proposed an empirical interatomic potential for covalent systems.
-The Tersoff potential explicitly incorporates the dependence of bond order on local envirenments, permitting an improved description of covalent materials.
-Tersoff applied the potential to silicon \cite{tersoff_silicon1,tersoff_silicon2,tersoff_silicon3}, carbon \cite{tersoff_carbon} and also to multicomponent systems like $SiC$ \cite{tersoff_multi}.
-The basic idea is that, in real systems, the bond order depends upon the local environment.
-An atom with many neighbours forms weaker bonds than an atom with few neighbours.
-
+\subsubsection{The Tersoff bond order potential}
+
+Tersoff proposed an empirical interatomic potential for covalent systems \cite{tersoff_si1,tersoff_si2}.
+The Tersoff potential explicitly incorporates the dependence of bond order on local environments, permitting an improved description of covalent materials.
+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.
+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}.
+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}.
+Atoms with many neighbors form weaker bonds than atoms with only a few neighbors.
+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.
+If the energy per bond decreases rapidly enough with increasing coordination the most stable structure will be the dimer.
+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.
+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.
+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}.
+
+Tersoff incorporated the concept of bond order in a three-body potential formalism.
The interatomic potential is taken to have the form
\begin{eqnarray}
E & = & \sum_i E_i = \frac{1}{2} \sum_{i \ne j} V_{ij} \textrm{ ,} \\
\end{eqnarray}
$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}$.
The indices $i$ and $j$ correspond to the atoms of the system with $r_{ij}$ being the distance from atom $i$ to atom $j$.
-
The functions $f_R$ and $f_A$ represent a repulsive and an attractive pair potential.
The repulsive part is due to the orthogonalization energy of overlapped atomic wave functions.
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 potential cutoff function to limit the range of the potential.
+The function $f_C$ is the 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\{
\frac{1}{2} + \frac{1}{2} \cos \Big[ \pi (r_{ij} - R_{ij})/(S_{ij} - R_{ij}) \Big], & R_{ij} < r_{ij} < S_{ij} \\
0, & r_{ij} > S_{ij}
\end{array} \right.
+\label{eq:basics:fc}
\end{equation}
-
-The function $b_{ij}$ represents a measure of the bond order, monotonically decreasing with the coordination of atoms $i$ and $j$.
+As discussed above, $b_{ij}$ represents a measure of the bond order, monotonously decreasing with the coordination of atoms $i$ and $j$.
It is of the form:
\begin{eqnarray}
b_{ij} & = & \chi_{ij} (1 + \beta_i^{n_i} \zeta^{n_i}_{ij})^{-1/2n_i} \\
\zeta_{ij} & = & \sum_{k \ne i,j} f_C (r_{ik}) \omega_{ik} g(\theta_{ijk}) \\
-g(\theta_{ijk}) & = & 1 + c_i^2/d_i^2 - c_i^2/[d_i^2 + (h_i - \cos \theta_{ijk})^2] \\
+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}.
-
-\printimg{!h}{width=8cm}{tersoff_angle.eps}{Angle between bonds of atoms $i,j$ and $i,k$.}{img:tersoff_angle}
-
-Here comes an explanation, energy per bond monotonically decreasing with the amount of bonds and so on and so on \ldots
+\begin{figure}[t]
+\begin{center}
+\includegraphics[width=8cm]{tersoff_angle.eps}
+\end{center}
+\caption{Angle between bonds of atoms $i,j$ and $i,k$.}
+\label{img:tersoff_angle}
+\end{figure}
+The angular dependence does not give a fixed minimum angle between bonds since the expression is embedded inside the bond order term.
+The relation to the above discussed bond order potential becomes obvious if $\chi=1, \beta=1, n=1, \omega=1$ and $c=0$.
+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.
+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.
+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.
The force acting on atom $i$ is given by the derivative of the potential energy.
For a three body potential ($V_{ij} \neq V{ji}$) the derivation is of the form
\begin{equation}
F^i = - \nabla_{{\bf r}_i} E \textrm{ .}
\end{equation}
-The details of the Tersoff potential derivative can be seen in appendix \ref{app:d_tersoff}.
-
-\subsubsection{A reparametrized Tersoff-like bond order potential}
+Details of the Tersoff potential derivative are presented in appendix \ref{app:d_tersoff}.
+
+\subsubsection{Improved analytical bond order potential}
+
+Although the Tersoff potential is one of the most widely used potentials there are some shortcomings.
+Describing the Si-Si interaction Tersoff was unable to find a single parameter set to describe well both, bulk and surface properties.
+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}.
+While T2 describes well surface properties, T3 yields improved elastic constants and should be used for describing bulk properties.
+However, T3, which is used in the Si/C potential, suffers from an underestimation of the dimer binding energy.
+Similar behavior is found for the C-C interaction.
+
+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}.
+The functional form is similar to the one proposed by Tersoff.
+Differences in the energy functional and the force evaluation routine are pointed out in appendix \ref{app:d_tersoff}.
+Concerning Si the elastic properties of the diamond phase as well as the structure and energetics of the dimer are reproduced very well.
+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.
+The potential succeeds in the description of the low as well as high coordinated structures.
+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.
+
+\subsection{Verlet integration}
+\label{subsection:integrate_algo}
-Erhart-Albe potential ...
+A numerical method to integrate Newton's equation 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}
+\vec{r}_i(t+\delta t)=
+\vec{r}_i(t)+\delta t\vec{v}_i(t)+\frac{\delta t^2}{2m_i}\vec{f}_i(t)+
+\frac{\delta t^3}{6}\vec{b}_i(t) + \mathcal{O}(\delta t^4)
+\label{basics:verlet:taylor1}
+\end{equation}
+\begin{equation}
+\vec{r}_i(t-\delta t)=
+\vec{r}_i(t)-\delta t\vec{v}_i(t)+\frac{\delta t^2}{2m_i}\vec{f}_i(t)-
+\frac{\delta t^3}{6}\vec{b}_i(t) + \mathcal{O}(\delta t^4)
+\label{basics:verlet:taylor2}
+\end{equation}
+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.
+The Verlet algorithm is obtained by summarizing and substracting equations \eqref{basics:verlet:taylor1} and \eqref{basics:verlet:taylor2}
+\begin{equation}
+\vec{r}_i(t+\delta t)=
+2\vec{r}_i(t)-\vec{r}_i(t-\delta t)+\frac{\delta t^2}{m_i}\vec{f}_i(t)+
+\mathcal{O}(\delta t^4)
+\end{equation}
+\begin{equation}
+\vec{v}_i(t)=\frac{1}{2\delta t}[\vec{r}_i(t+\delta t)-\vec{r}_i(t-\delta t)]+
+\mathcal{O}(\delta t^3)
+\end{equation}
+the truncation error of which is of order $\delta t^4$ for the positions and $\delta t^3$ for the velocities.
+The velocities, although not used to update the particle positions, are not synchronously determined with the positions but drag behind one step of discretization.
+The Verlet algorithm can be rewritten into an equivalent form, which updates the velocities and positions in the same step.
+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$
+\begin{equation}
+\vec{v}_i(t+\delta t)=
+\vec{v}_i(t)+\frac{\delta t}{2m_i}[\vec{f}_i(t)+\vec{f}_i(t+\delta t)]
+\end{equation}
+\begin{equation}
+\vec{r}_i(t+\delta t)=
+\vec{r}_i(t)+\delta t\vec{v}_i(t)+\frac{\delta t^2}{2m_i}\vec{f}_i(t) \text{ .}
+\end{equation}
+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.
\subsection{Statistical ensembles}
\label{subsection:statistical_ensembles}
+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).
+Providing a stable integration algorithm the total energy $E$, i.e. the kinetic and configurational energy of the paticles, is conserved.
+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.
+
+However, the successful formation of SiC dictates precise control of temperature by external heating.
+While the temperature of such a system is well defined, the energy is no longer conserved.
+The microscopic states of a system, which is in thermal equilibrium with an external thermal heat bath, are represented by the $NVT$ ensemble.
+In the so-called canonical ensemble the temperature $T$ is related to the expactation value of the kinetic energy of the particles, i.e.
+\begin{equation}
+\langle E_{\text{kin}}\rangle = \frac{3}{2}Nk_{\text{B}}T \text{, }
+E_{\text{kin}}=\sum_i \frac{\vec{p}^2_i}{2m_i} \text{ .}
+\label{eq:basics:ts}
+\end{equation}
+The volume of the synthesized material can hardly be controlled in experiment.
+Instead the pressure can be adjusted.
+Holding constant the pressure in addition to the temperature of the system its states are represented by the isothermal-isobaric $NpT$ ensemble.
+The expression for the pressure of a system derived from the equipartition theorem is given by
+\begin{equation}
+pV=Nk_{\text{B}}T+\langle W\rangle\text{, }W=-\frac{1}{3}\sum_i\vec{r}_i\nabla_{\vec{r}_i}U
+\text{, }
+\label{eq:basics:ps}
+\end{equation}
+where $W$ is the virial and $U$ is the configurational energy.
+
+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.
+Control of the respective variable is based on the relations given in equations \eqref{eq:basics:ts} and \eqref{eq:basics:ps}.
+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
+\begin{equation}
+\lambda=\left[1+\frac{\delta t}{\tau_T}(\frac{T_0}{T}-1)\right]^\frac{1}{2}
+\text{ ,}
+\end{equation}
+where $T$ is the current temperature according to equation \eqref{eq:basics:ts}.
+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
+\begin{equation}
+\mu=\left[1-\frac{\beta\delta t}{\tau_p}(p_0-p)\right]^\frac{1}{3}\text{ ,}
+\end{equation}
+where $\beta$ is the isothermal compressibility and $p$ corresponds to the current pressure, which is determined by equation \eqref{eq:basics:ps}.
+
+Using this method the system does not behave like a true $NpT$ ensemble.
+On average $T$ and $p$ correspond to the expected values.
+For large enough time constants, i.e. $\tau > 100 \delta t$, the method shows realistic fluctuations in $T$ and $p$.
+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.
+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.
+
\section{Denstiy functional theory}
\label{section:dft}
-\subsection{Born-Oppenheimer (adiabatic) approximation}
+Dirac declared that chemistry has come to an end, its content being entirely contained in the powerul equation published by Schr\"odinger in 1926 \cite{schroeder26} marking the beginning of wave mechanics.
+Following the path of Schr\"odinger the problem in quantum-mechanical modeling of describing the many-body problem, i.e. a system of a large amount of interacting particles, is manifested in the high-dimensional Schr\"odinger equation for the wave function $\Psi({\vec{R}},{\vec{r}})$ that depends on the coordinates of all nuclei and electrons.
+The Schr\"odinger equation contains the kinetic energy of the ions and electrons as well as the electron-ion, ion-ion and electron-electron interaction.
+This cannot be solved exactly and finding approximate solutions requires several layers of simplification in order to reduce the number of free parameters.
+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.
+
+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.
+In the following sections the basic idea of DFT will be outlined.
+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}.
+
+\subsection{Born-Oppenheimer approximation}
+
+Born and Oppenheimer proposed a simplification enabling the effective decoupling of the electronic and ionic degrees of freedom \cite{born27}.
+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.
+Thus, on the timescale of electronic motion the ions appear at fixed positions and, on the other way round, for the nuclei the electrons appear blurred in space adding an extra term to the ion-ion potential.
+The simplified Schr\"odinger equation is rewritten without the kinetic energy of the ions and its positions enter as fixed parameters.
+
+\subsection{Bloch theorem}
\subsection{Hohenberg-Kohn theorem}
-\subsection{Exchange correlation}
+\subsection{Effective potential}
+
+\subsection{Kohn-Sham system}
+
+\subsection{Approximations for exchange and correlation}
\subsection{Pseudopotentials}
+\section{Modeling of defects}
+\label{section:basics:defects}
+
+\section{Migration paths and diffusion barriers}
+\label{section:basics:migration}
+
projects / lectures / latex.git / blobdiff