posic/thesis/basics.tex

   1 \chapter{Basic principles of utilized simulation techniques}
   2 \label{chapter:basics}
   3
   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}}\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.
   8 The basic ideas of MD in general and the adopted techniques as implemented in {\em 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 system to be described 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}.
  18
  19 \section{Molecular dynamics simulations}
  20 \label{section:md}
  21
  22 % todo
  23 % rewrite!
  24
  25 \begin{quotation}
  26 \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 \begin{flushright}
  28 {\em Marquis Pierre Simon de Laplace, 1814.} \cite{laplace}
  29 \end{flushright}
  30 \end{quotation}
  31
  32 \noindent
  33 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.
  34 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.
  35 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
  37 \subsection{Introduction to molecular dynamics simulations}
  38
  39 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.
  40 The MD method was first introduced by Alder and Wainwright in 1957 \cite{alder57,alder59} to study the interactions of hard spheres.
  41 The basis of the approach are Newton's equations of motion to describe classicaly the many-body system.
  42 MD is the numerical way of solving the $N$-body problem which cannot be solved analytically for $N>3$.
  43 A potential is necessary describing the interaction of the particles.
  44 By MD a complete description of the system in the sense of classical mechanics on the microscopic level is obtained.
  45 The microscopic information can then be translated to macroscopic observables by means of statistical mechanics.
  46
  47 The basic idea is to integrate Newton's equations numerically.
  48 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 \begin{equation}
  50 m_i \frac{d^2}{dt^2} {\bf r}_i = {\bf F}_i \, \textrm{.}
  51 \end{equation}
  52 The forces ${\bf F}_i$ are obtained from the potential energy $U(\{{\bf r}\})$:
  53 \begin{equation}
  54 {\bf F}_i = - \nabla_{{\bf r}_i} U({\{\bf r}\}) \, \textrm{.}
  55 \end{equation}
  56 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.
  57 The solution of these equations provides the complete information of a system evolving in time.
  58 The following sections cover the tools of the trade necessary for the MD simulation technique.
  59 Three ingredients are required for a MD simulation:
  60 \begin{enumerate}
  61 \item A model for the interaction between system constituents is needed.
  62       Interaction potentials and their accuracy for describing certain systems of elements will be outlined in section \ref{subsection:interact_pot}.
  63 \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.
  64       This is explained in section \ref{subsection:integrate_algo}.
  65 \item A statistical ensemble has to be chosen, which allows certain thermodynamic quantities to be controlled or to stay constant.
  66       This is discussed in section \ref{subsection:statistical_ensembles}.
  67 \end{enumerate}
  68 Furthermore special techniques will be outlined, which reduce the complexity of the MD algorithm, though the evaluation of energy and force almost inevitably dictates the overall speed.
  69
  70 \subsection{Verlet integration}
  71 \label{subsection:integrate_algo}
  72
  73 A numerical method to integrate Newton's equation of motion was presented by Verlet in 1967 \cite{verlet67}.
  74 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.
  75 Starting point is the Taylor series for the particle positions at time $t+\delta t$ and $t-\delta t$
  76 \begin{equation}
  77 \vec{r}_i(t+\delta t)=
  78 \vec{r}_i(t)+\delta t\vec{v}_i(t)+\frac{\delta t^2}{2m_i}\vec{f}_i(t)+
  79 \frac{\delta t^3}{6}\vec{b}_i(t) + \mathcal{O}(\delta t^4)
  80 \label{basics:verlet:taylor1}
  81 \end{equation}
  82 \begin{equation}
  83 \vec{r}_i(t-\delta t)=
  84 \vec{r}_i(t)-\delta t\vec{v}_i(t)+\frac{\delta t^2}{2m_i}\vec{f}_i(t)-
  85 \frac{\delta t^3}{6}\vec{b}_i(t) + \mathcal{O}(\delta t^4)
  86 \label{basics:verlet:taylor2}
  87 \end{equation}
  88 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.
  89 The Verlet algorithm is obtained by summarizing and substracting equations \eqref{basics:verlet:taylor1} and \eqref{basics:verlet:taylor2}
  90 \begin{equation}
  91 \vec{r}_i(t+\delta t)=
  92 2\vec{r}_i(t)-\vec{r}_i(t-\delta t)+\frac{\delta t^2}{m_i}\vec{f}_i(t)+
  93 \mathcal{O}(\delta t^4)
  94 \end{equation}
  95 \begin{equation}
  96 \vec{v}_i(t)=\frac{1}{2\delta t}[\vec{r}_i(t+\delta t)-\vec{r}_i(t-\delta t)]+
  97 \mathcal{O}(\delta t^3)
  98 \end{equation}
  99 the truncation error of which is of order $\delta t^4$ for the positions and $\delta t^3$ for the velocities.
 100 The velocities, although not used to update the particle positions, are not synchronously determined with the positions but drag behind one step of discretization.
 101 The Verlet algorithm can be rewritten into an equivalent form, which updates the velocities and positions in the same step.
 102 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$
 103 \begin{equation}
 104 \vec{v}_i(t+\delta t)=
 105 \vec{v}_i(t)+\frac{\delta t}{2m_i}[\vec{f}_i(t)+\vec{f}_i(t+\delta t)]
 106 \end{equation}
 107 \begin{equation}
 108 \vec{r}_i(t+\delta t)=
 109 \vec{r}_i(t)+\delta t\vec{v}_i(t)+\frac{\delta t^2}{2m_i}\vec{f}_i(t) \text{ .}
 110 \end{equation}
 111 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.
 112
 113 \subsection{Interaction potentials for silicon and carbon}
 114 \label{subsection:interact_pot}
 115
 116 % todo
 117 % rewrite!
 118
 119 The potential energy of $N$ interacting atoms can be written in the form
 120 \begin{equation}
 121 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
 122 \end{equation}
 123 where $U$ is the total potential energy.
 124 $U_1$ is a single particle potential describing external forces.
 125 Examples of single particle potentials are the gravitational force or an electric field.
 126 $U_2$ is a two body pair potential which only depends on the distance $r_{ij}$ between the two atoms $i$ and $j$.
 127 If not only pair potentials are considered, three body potentials $U_3$ or multi body potentials $U_n$ can be included.
 128 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.
 129 Ordinary pair potentials have a close-packed structure like face-centered cubic (FCC) or hexagonal close-packed (HCP) as a ground state.
 130 A pair potential is, thus, unable to describe properly elements with other structures than FCC or HCP.
 131 Silicon and carbon for instance, have a diamand and zincblende structure with four covalently bonded neighbors, which is far from a close-packed structure.
 132 A three body potential has to be included for these types of elements.
 133
 134 \subsubsection{The Tersoff bond order potential}
 135
 136 Tersoff proposed an empirical interatomic potential for covalent systems \cite{tersoff_si1,tersoff_si2}.
 137 The Tersoff potential explicitly incorporates the dependence of bond order on local environments, permitting an improved description of covalent materials.
 138 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.
 139 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}.
 140 ...
 141 The basic idea is that, in real systems, the bond order depends upon the local environment.
 142 An atom with many neighbors forms weaker bonds than an atom with few neighbors since .
 143
 144 Here comes an explanation, energy per bond monotonically decreasing with the amount of bonds and so on and so on \ldots
 145
 146 The interatomic potential is taken to have the form
 147 \begin{eqnarray}
 148 E & = & \sum_i E_i = \frac{1}{2} \sum_{i \ne j} V_{ij} \textrm{ ,} \\
 149 V_{ij} & = & f_C(r_{ij}) [ f_R(r_{ij}) + b_{ij} f_A(r_{ij}) ] \textrm{ .}
 150 \end{eqnarray}
 151 $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}$.
 152 The indices $i$ and $j$ correspond to the atoms of the system with $r_{ij}$ being the distance from atom $i$ to atom $j$.
 153 The functions $f_R$ and $f_A$ represent a repulsive and an attractive pair potential.
 154 The repulsive part is due to the orthogonalization energy of overlapped atomic wave functions.
 155 The attractive part is associated with the bonding.
 156 \begin{eqnarray}
 157 f_R(r_{ij}) & = & A_{ij} \exp (- \lambda_{ij} r_{ij} ) \\
 158 f_A(r_{ij}) & = & -B_{ij} \exp (- \mu_{ij} r_{ij} )
 159 \end{eqnarray}
 160 The function $f_C$ is the a cutoff function to limit the range of interaction to nearest neighbors.
 161 It is designed to have a smooth transition of the potential at distances $R_{ij}$ and $S_{ij}$.
 162 \begin{equation}
 163 f_C(r_{ij}) = \left\{
 164   \begin{array}{ll}
 165     1, & r_{ij} < R_{ij} \\
 166     \frac{1}{2} + \frac{1}{2} \cos \Big[ \pi (r_{ij} - R_{ij})/(S_{ij} - R_{ij}) \Big], & R_{ij} < r_{ij} < S_{ij} \\
 167     0, & r_{ij} > S_{ij}
 168   \end{array} \right.
 169 \label{eq:basics:fc}
 170 \end{equation}
 171 The function $b_{ij}$ represents a measure of the bond order, monotonically decreasing with the coordination of atoms $i$ and $j$.
 172 It is of the form:
 173 \begin{eqnarray}
 174 b_{ij} & = & \chi_{ij} (1 + \beta_i^{n_i} \zeta^{n_i}_{ij})^{-1/2n_i} \\
 175 \zeta_{ij} & = & \sum_{k \ne i,j} f_C (r_{ik}) \omega_{ik} g(\theta_{ijk}) \\
 176 g(\theta_{ijk}) & = & 1 + c_i^2/d_i^2 - c_i^2/[d_i^2 + (h_i - \cos \theta_{ijk})^2]
 177 \end{eqnarray}
 178 where $\theta_{ijk}$ is the bond angle between bonds $ij$ and $ik$.
 179 This is illustrated in Figure \ref{img:tersoff_angle}.
 180 \begin{figure}[t]
 181 \begin{center}
 182 \includegraphics[width=8cm]{tersoff_angle.eps}
 183 \end{center}
 184 \caption{Angle between bonds of atoms $i,j$ and $i,k$.}
 185 \label{img:tersoff_angle}
 186 \end{figure}
 187
 188 The force acting on atom $i$ is given by the derivative of the potential energy.
 189 For a three body potential ($V_{ij} \neq V{ji}$) the derivation is of the form
 190 \begin{equation}
 191 \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{ .}
 192 \end{equation}
 193 The force is then given by
 194 \begin{equation}
 195 F^i = - \nabla_{{\bf r}_i} E \textrm{ .}
 196 \end{equation}
 197 Details of the Tersoff potential derivative are presented in appendix \ref{app:d_tersoff}.
 198
 199 \subsubsection{Improved analytical bond order potential}
 200
 201
 202 \subsection{Statistical ensembles}
 203 \label{subsection:statistical_ensembles}
 204
 205 \section{Denstiy functional theory}
 206 \label{section:dft}
 207
 208 \subsection{Hohenberg-Kohn theorem}
 209
 210 \subsection{Born-Oppenheimer (adiabatic) approximation}
 211
 212 \subsection{Effective potential}
 213
 214 \subsection{Kohn-Sham system}
 215
 216 \subsection{Approximations for exchange and correlation}
 217
 218 \subsection{Pseudopotentials}
 219
 220 \section{Modeling of defects}
 221 \label{section:basics:defects}
 222
 223 \section{Migration paths and diffusion barriers}
 224 \label{section:basics:migration}
 225