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}} 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}.
+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/posic}.
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.
\section{Molecular dynamics simulations}
\label{section:md}
-% todo
-% rewrite!
+% todo - rewrite md intro chapter
\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{}
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 a cutoff function to limit the range of interaction to nearest neighbors.
+The function $f_C$ is 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\{
\subsection{Verlet integration}
\label{subsection:integrate_algo}
-A numerical method to integrate Newton's equation of motion was presented by Verlet in 1967 \cite{verlet67}.
+A numerical method to integrate Newton's equations 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}
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.
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$.
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|$.
-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
+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
\begin{equation}
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}
\text{ .}
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.
\subsection{Brillouin zone sampling}
+\label{subsection:basics:bzs}
Following Bloch's theorem only a finite number of electronic wave functions need to be calculated for a periodic system.
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.
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.
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.
+However, there is hardly ever only one defect in a crystal, not even only one kind of defect.
+Again, energetic considerations can be used to investigate the existing interaction of two defects.
+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}}}$.
+This can be expressed by
+\begin{equation}
+E_{\text{b}}=
+E_{\text{f}}^{\text{comb}}-
+E_{\text{f}}^{1^{\text{st}}}-
+E_{\text{f}}^{2^{\text{nd}}}
+\label{eq:basics:e_bind}
+\end{equation}
+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.
+Accordingly, energetically favorable configurations result in binding energies below zero while unfavorable configurations show positive values for the binding energy.
+The interaction strength, i.e. the absolute value of the binding energy, approaches zero for increasingly non-interacting isolated defects.
+Thus, $E_{\text{b}}$ indeed can be best thought of a binding energy, which is required to bring the defects to infinite separation.
+
The methods presented in the last two chapters can be used to investigate defect structures and energetics.
Therefore, a supercell containing the perfect crystal is generated in an initial process.
If not by construction, the system should be fully relaxed.
\section{Migration paths and diffusion barriers}
\label{section:basics:migration}
-Investigating diffusion mechanisms is based on determining migration paths inbetween two local minimum configurations of an atom at different locations in the lattice.
+Investigating diffusion mechanisms is based on determining migration paths in between two local minimum configurations of an atom at different locations in the lattice.
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.
The maximum difference in energy is the barrier necessary for the respective migration process.
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.
\begin{figure}[t]
\begin{center}
-\subfigure[]{\label{fig:basics:crto}
-\includegraphics[width=0.5\textwidth]{crt_orig.eps}}
-\subfigure[]{\label{fig:basics:crtm}
-\includegraphics[width=0.5\textwidth]{crt_mod.eps}}
+\subfigure[]{\label{fig:basics:crto}\includegraphics[width=0.45\textwidth]{crt_orig.eps}}
+\subfigure[]{\label{fig:basics:crtm}\includegraphics[width=0.45\textwidth]{crt_mod.eps}}
\end{center}
\caption{Schematic of the constrained relaxation technique (a) and of a modified version (b) used to obtain migration pathways and corresponding configurational energies.}
\label{fig:basics:crt}
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.
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.
-
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}.
-Due to these constraints obtained activation energies can effectively be higher.
-% todo
-% advantages of pw basis with respect to hellmann feynman forces / pulay forces
+% todo - advantages of pw basis concenring hf forces + inc font in crt sketch
projects / lectures / latex.git / blobdiff