\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 were chosen in order to model the behavior and precipitation of C introduced into an initially crystalline Si environment.
+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 {\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}.
\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{}
\begin{flushright}
\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 sections cover the tools of the trade necessary for the MD simulation technique.
-Three ingredients are required for an MD simulation:
+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 section \ref{subsection:interact_pot}.
\item A statistical ensemble has to be chosen, which allows certain thermodynamic quantities to be controlled or to stay constant.
This is discussed in section \ref{subsection:statistical_ensembles}.
\end{enumerate}
-Furthermore special techniques will be outlined which reduce the complexity of the MD algorithm, though the evaluation of the energy and force almost inevitably dictates the overall speed.
+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.
-\subsection{Integration algorithms}
+\subsection{Verlet integration}
\label{subsection:integrate_algo}
-\subsection{Interaction potentials}
+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{Interaction potentials for silicon and carbon}
\label{subsection:interact_pot}
+% todo
+% rewrite!
+
The potential energy of $N$ interacting atoms can be written in the form
\begin{equation}
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
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 and 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.
-%\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_si1,tersoff_si2,tersoff_si3}, carbon \cite{tersoff_c} and also to multicomponent systems like $SiC$ \cite{tersoff_m}.
+\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 depends upon the local environment.
-An atom with many neighbours forms weaker bonds than an atom with few neighbours.
+An atom with many neighbors forms weaker bonds than an atom with few neighbors since .
+
+Here comes an explanation, energy per bond monotonically decreasing with the amount of bonds and so on and so on \ldots
The interatomic potential is taken to have the form
\begin{eqnarray}
\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\{
\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$.
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 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}.
+Details of the Tersoff potential derivative are presented in appendix \ref{app:d_tersoff}.
-\subsubsection{A reparametrized Tersoff-like bond order potential}
+\subsubsection{Improved analytical bond order potential}
-Erhart-Albe potential ...
\subsection{Statistical ensembles}
\label{subsection:statistical_ensembles}
projects / lectures / latex.git / commitdiff