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 simulations were chosen in order to model 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} an appropriate MD code called {\em posic ({\bf P}recipitation {\bf o}f {\bf SiC})}\footnote{The source code will be made available for download at 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 In addition to the classical potential approach ...
  11 Highly accurate technique DFT, VASP code, VASP mods, tools in appendix ...
  12
  13 Determining the formation energies if defects is ...
  14
  15 Migration pathways were investigated using the , which is explained in cloder detail in ...
  16
  17
  18 \section{Molecular dynamics simulations}
  19 \label{section:md}
  20
  21 \begin{quotation}
  22 \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{}
  23 \begin{flushright}
  24 {\em Marquis Pierre Simon de Laplace, 1814, \cite{laplace}.}
  25 \end{flushright}
  26 \end{quotation}
  27
  28 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.
  29 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.
  30 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.
  31
  32 \subsection{Introduction to molecular dynamics simulations}
  33
  34 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.
  35 The MD method was first introduced by Alder and Wainwright in 1957 \cite{alder57,alder59} to study the interactions of hard spheres.
  36 The basis of the approach are Newton's equations of motion to describe classicaly the many-body system.
  37 MD simulation is the numerical way of solving the $N$-body problem which cannot be solved analytically ($N>3$).
  38 Quantum mechanical effects are taken into account by an analytical interaction potential between the nuclei.
  39
  40 By MD simulation techniques a complete description of the system in the sense of classical mechanics on the microscopic level is obtained.
  41 This microscopic information has to be translated to macroscopic observables by means of statistical mechanics.
  42
  43 The basic idea is to integrate Newton's equations numerically.
  44 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
  45 \begin{equation}
  46 m_i \frac{d^2}{dt^2} {\bf r}_i = {\bf F}_i \, \textrm{.}
  47 \end{equation}
  48 The forces ${\bf F}_i$ are obtained from the potential energy $U(\{{\bf r}\})$:
  49 \begin{equation}
  50 {\bf F}_i = - \nabla_{{\bf r}_i} U({\{\bf r}\}) \, \textrm{.}
  51 \end{equation}
  52
  53 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.
  54 The solution of these equations provides the complete information of a system evolving in time.
  55
  56 The following chapters cover the tools of the trade necessary for the MD simulation technique.
  57 Three ingredients are required for an MD simulation:
  58 \begin{enumerate}
  59 \item A model for the interaction between system constituents is needed.
  60       Interaction potentials and their accuracy for describing certain systems of elements will be outlined in chapter \ref{subsection:interact_pot}.
  61 \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.
  62       In chapter \ref{subsection:integrate_algo} a detailed overview of the available integration algorithms is given, including their advantages and disadvantages.
  63 \item A statistical ensemble has to be chosen, which allows certain thermodynamic quantities to be controlled or to stay constant.
  64       This is discussed in chapter \ref{subsection:statistical_ensembles}
  65 \end{enumerate}
  66
  67 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.
  68
  69 \subsection{Integration algorithms}
  70 \label{subsection:integrate_algo}
  71
  72 \subsection{Interaction potentials}
  73 \label{subsection:interact_pot}
  74
  75 The potential energy of $N$ interacting atoms can be written in the form
  76 \begin{equation}
  77 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
  78 \end{equation}
  79 where $U$ is the total potential energy.
  80 $U_1$ is a single particle potential describing external forces.
  81 Examples of single particle potentials are the gravitational force or an electric field.
  82 $U_2$ is a two body pair potential which only depends on the distance $r_{ij}$ between the two atoms $i$ and $j$.
  83 If not only pair potentials are considered, three body potentials $U_3$ or multi body potentials $U_n$ can be included.
  84 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.
  85
  86 Ordinary pair potentials have a close-packed structure like face-centered cubic (FCC) or hexagonal close-packed (HCP) as a ground state.
  87 A pair potential is thus unable to describe properly elements with other structures than FCC or HCP.
  88 Silicon and carbon for instance, have a diamand/zincblende structure with four covalent bonded neighbours, which is far from a close-packed structure.
  89 A three body potential has to be included for these types of elements.
  90
  91 In the following, relevant potentials for this work are discussed.
  92
  93 %\subsubsection{The Lennard-Jones potential}
  94 %
  95 %The L-J potential is a realistic two body pair potential and is of the form
  96 %\begin{equation}
  97 %U^{LJ}(r) = 4 \epsilon \Big[ \Big( \frac{\sigma}{r} \Big)^{12} - \Big( \frac{\sigma}{r} \Big)^6 \Big] \, \textrm{,}
  98 %\label{eq:lj-p}
  99 %\end{equation}
 100 %where $r$ denotes the distance between the two atoms.
 101 %
 102 %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.
 103 %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.
 104 %
 105 %The repulsive term $(\sim r^{-12})$ captures the non-bonded overlap of the electron clouds.
 106 %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.
 107 %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.
 108 %
 109 %The constants $\epsilon$ and $\sigma$ are usually determined by fitting to experimental data.
 110 %$\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.
 111 %
 112 %Writing down the derivative of the Lennard-Jones potential in respect to $x_i$ (the $i$th component of the distance vector ${\bf r}$)
 113 %\begin{equation}
 114 %\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)
 115 %\label{eq:lj-d}
 116 %\end{equation}
 117 %one can easily identify $\sigma$ by the equilibrium distance of the atoms $r_e=\sqrt[6]{2} \sigma$.
 118 %Applying the equilibrium distance into \eqref{eq:lj-p} $\epsilon$ turns out to be the negative well depth.
 119 %The $i$th component of the force is given by
 120 %\begin{equation}
 121 %F_i = - \frac{\partial}{\partial x_i} U^{LJ}(r) \, \textrm{.}
 122 %\label{eq:lj-f}
 123 %\end{equation}
 124 %
 125 %\subsubsection{The Stillinger Weber potential}
 126 %
 127 %The Stillinger Weber potential \cite{stillinger_weber} \ldots
 128 %
 129 %\begin{equation}
 130 %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)
 131 %\end{equation}
 132 %
 133 %\begin{equation}
 134 %U_2(r_{ij}) = \left\{
 135 %  \begin{array}{ll}
 136 %    \epsilon A \Big( B (r_{ij} / \sigma)^{-p} - 1\Big) \exp \Big[ (r_{ij} / \sigma - 1)^{-1}  \Big]  & r_{ij} < a \sigma \\
 137 %    0 & r_{ij} \ge a \sigma
 138 %  \end{array} \right.
 139 %\end{equation}
 140 %
 141 %\begin{equation}
 142 %U_3({\bf r}_i,{\bf r}_j,{\bf r}_k) =
 143 %\epsilon \Big[ h(r_{ij},r_{ik},\theta_{jik}) + h(r_{ji},r_{jk},\theta_{ijk}) + h(r_{ki},r_{kj},\theta_{ikj}) \Big]
 144 %\end{equation}
 145 %
 146 %\begin{equation}
 147 %h(r_{ij},r_{ik},\theta_{jik}) =
 148 %\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
 149 %\end{equation}
 150
 151 \subsubsection{The Tersoff potential}
 152
 153 Tersoff proposed an empirical interatomic potential for covalent systems.
 154 The Tersoff potential explicitly incorporates the dependence of bond order on local envirenments, permitting an improved description of covalent materials.
 155 Tersoff applied the potential to silicon \cite{tersoff_si1,tersoff_si2,tersoff_si3}, carbon \cite{tersoff_c} and also to multicomponent systems like $SiC$ \cite{tersoff_m}.
 156 The basic idea is that, in real systems, the bond order depends upon the local environment.
 157 An atom with many neighbours forms weaker bonds than an atom with few neighbours.
 158
 159 The interatomic potential is taken to have the form
 160 \begin{eqnarray}
 161 E & = & \sum_i E_i = \frac{1}{2} \sum_{i \ne j} V_{ij} \textrm{ ,} \\
 162 V_{ij} & = & f_C(r_{ij}) [ f_R(r_{ij}) + b_{ij} f_A(r_{ij}) ] \textrm{ .}
 163 \end{eqnarray}
 164 $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}$.
 165 The indices $i$ and $j$ correspond to the atoms of the system with $r_{ij}$ being the distance from atom $i$ to atom $j$.
 166
 167 The functions $f_R$ and $f_A$ represent a repulsive and an attractive pair potential.
 168 The repulsive part is due to the orthogonalization energy of overlapped atomic wave functions.
 169 The attractive part is associated with the bonding.
 170 \begin{eqnarray}
 171 f_R(r_{ij}) & = & A_{ij} \exp (- \lambda_{ij} r_{ij} ) \\
 172 f_A(r_{ij}) & = & -B_{ij} \exp (- \mu_{ij} r_{ij} )
 173 \end{eqnarray}
 174 The function $f_C$ is the potential cutoff function to limit the range of the potential.
 175 It is designed to have a smooth transition of the potential at distances $R_{ij}$ and $S_{ij}$.
 176 \begin{equation}
 177 f_C(r_{ij}) = \left\{
 178   \begin{array}{ll}
 179     1, & r_{ij} < R_{ij} \\
 180     \frac{1}{2} + \frac{1}{2} \cos \Big[ \pi (r_{ij} - R_{ij})/(S_{ij} - R_{ij}) \Big], & R_{ij} < r_{ij} < S_{ij} \\
 181     0, & r_{ij} > S_{ij}
 182   \end{array} \right.
 183 \label{eq:basics:fc}
 184 \end{equation}
 185
 186 The function $b_{ij}$ represents a measure of the bond order, monotonically decreasing with the coordination of atoms $i$ and $j$.
 187 It is of the form:
 188 \begin{eqnarray}
 189 b_{ij} & = & \chi_{ij} (1 + \beta_i^{n_i} \zeta^{n_i}_{ij})^{-1/2n_i} \\
 190 \zeta_{ij} & = & \sum_{k \ne i,j} f_C (r_{ik}) \omega_{ik} g(\theta_{ijk}) \\
 191 g(\theta_{ijk}) & = & 1 + c_i^2/d_i^2 - c_i^2/[d_i^2 + (h_i - \cos \theta_{ijk})^2] \\
 192 \end{eqnarray}
 193 where $\theta_{ijk}$ is the bond angle between bonds $ij$ and $ik$.
 194 This is illustrated in Figure \ref{img:tersoff_angle}.
 195
 196 \printimg{!h}{width=8cm}{tersoff_angle.eps}{Angle between bonds of atoms $i,j$ and $i,k$.}{img:tersoff_angle}
 197
 198 Here comes an explanation, energy per bond monotonically decreasing with the amount of bonds and so on and so on \ldots
 199
 200 The force acting on atom $i$ is given by the derivative of the potential energy.
 201 For a three body potential ($V_{ij} \neq V{ji}$) the derivation is of the form
 202 \begin{equation}
 203 \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{ .}
 204 \end{equation}
 205 The force is then given by
 206 \begin{equation}
 207 F^i = - \nabla_{{\bf r}_i} E \textrm{ .}
 208 \end{equation}
 209 The details of the Tersoff potential derivative can be seen in appendix \ref{app:d_tersoff}.
 210
 211 \subsubsection{A reparametrized Tersoff-like bond order potential}
 212
 213 Erhart-Albe potential ...
 214
 215 \subsection{Statistical ensembles}
 216 \label{subsection:statistical_ensembles}
 217
 218 \section{Denstiy functional theory}
 219 \label{section:dft}
 220
 221 \subsection{Hohenberg-Kohn theorem}
 222
 223 \subsection{Born-Oppenheimer (adiabatic) approximation}
 224
 225 \subsection{Effective potential}
 226
 227 \subsection{Kohn-Sham system}
 228
 229 \subsection{Approximations for exchange and correlation}
 230
 231 \subsection{Pseudopotentials}
 232
 233 \section{Modeling of defects}
 234 \label{section:basics:defects}
 235
 236 \section{Migration paths and diffusion barriers}
 237 \label{section:basics:migration}
 238