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/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}.
+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.
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 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}.
+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}
\end{quotation}
\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.
+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 may be regarded as the basic principle of molecular dynamics, considering an isolated system of particles, the behavior of which is fully determined by the solution of the classical equations of motion.
$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 many body potentials $U_n$ can be included.
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.
+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 diamond 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.
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.
+Due to the covalent character Tersoff restricted the interaction to nearest neighbor atoms accompanied by an increase 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.
+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.
\subsubsection{Improved analytical bond order potential}
-Although the Tersoff potential is one of the most widely used potentials there are some shortcomings.
+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.
\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.
+where $\vec{v}_i=\frac{d}{dt}\vec{r}_i$ are the velocities, $\vec{f}_i=m\frac{d^2}{dt^2}\vec{r}_i$ are the forces and $\vec{b}_i=\frac{d^3}{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 subtracting equations \eqref{basics:verlet:taylor1} and \eqref{basics:verlet:taylor2}
\begin{equation}
\vec{r}_i(t+\delta t)=
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$
+The so-called velocity Verlet algorithm is obtained by combining equation \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)]
++\mathcal{O}(\delta t^3)
\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{ .}
+\vec{r}_i(t)+\delta t\vec{v}_i(t)+\frac{\delta t^2}{2m_i}\vec{f}_i(t)
++\mathcal{O}(\delta t^3)
+ \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.
\label{section:dft}
Dirac declared that chemistry has come to an end, its content being entirely contained in the powerful equation published by Schr\"odinger in 1926 \cite{schroedinger26} 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.
+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 determinants to be evaluated massively increases.
\subsection{Hohenberg-Kohn theorem and variational principle}
Investigating the energetics of Cu$_x$Zn$_{1-x}$ alloys, which for different compositions exhibit different transfers of charge between the Cu and Zn unit cells due to their chemical difference and, thus, varying electrostatic interactions contributing to the total energy, the work of Hohenberg and Kohn had a natural focus on the distribution of charge.
-Although it was clear that the Thomas Fermi (TF) theory only provides a rough approximation to the exact solution of the many-electron Schr\"odinger equation the theory was of high interest since it provides an implicit relation of the potential and the electron density distribution.
+Although it was clear that the Thomas Fermi (TF) theory only provides a rough approximation to the exact solution of the many-electron Schr\"odinger equation, the theory was of high interest since it provides an implicit relation of the potential and the electron density distribution.
This raised the question how to establish a connection between TF expressed in terms of $n(\vec{r})$ and the exact Schr\"odinger equation expressed in terms of the many-electron wave function $\Psi({\vec{r}})$ and whether a complete description in terms of the charge density is possible in principle.
The answer to this question, whether the charge density completely characterizes a system, became the starting point of modern DFT.
\subsection{Kohn-Sham system}
Inspired by the Hartree equations, i.e. a set of self-consistent single-particle equations for the approximate solution of the many-electron problem \cite{hartree28}, which describe atomic ground states much better than the TF theory, Kohn and Sham presented a Hartree-like formulation of the Hohenberg and Kohn minimal principle \eqref{eq:basics:hkm} \cite{kohn65}.
-However, due to a more general approach, the new formulation is formally exact by introducing the energy functional $E_{\text{xc}}[n(vec{r})]$, which accounts for the exchange and correlation energy of the electron interaction $U$ and possible corrections due to electron interaction to the kinetic energy $T$.
+However, due to a more general approach, the new formulation is formally exact by introducing the energy functional $E_{\text{xc}}[n(\vec{r})]$, which accounts for the exchange and correlation energy of the electron interaction $U$ and possible corrections due to electron interaction to the kinetic energy $T$.
The respective Kohn-Sham equations for the effective single-particle wave functions $\Phi_i(\vec{r})$ take the form
\begin{equation}
\left[
\end{equation}
\begin{equation}
V_{\text{eff}}(\vec{r})=V(\vec{r})+\int\frac{e^2n(\vec{r}')}{|\vec{r}-\vec{r}'|}d\vec{r}'
- + V_{\text{xc}(\vec{r})}
+ + V_{\text{xc}}(\vec{r})
\text{ ,}
\label{eq:basics:kse2}
\end{equation}
Although, even in such a simple case, no exact form of the correlation part of $\epsilon_{\text{xc}}(n)$ is known, highly accurate numerical estimates using Monte Carlo methods \cite{ceperley80} and corresponding parametrizations exist \cite{perdew81}.
Obviously exact for the homogeneous electron gas, the LDA was {\em a priori} expected to be useful only for densities varying slowly on scales of the local Fermi or TF wavelength.
Nevertheless, LDA turned out to be extremely successful in describing some properties of highly inhomogeneous systems accurately within a few percent.
-Although LDA is known to overestimate the cohesive energy in solids by \unit[10-20]{\%}, the lattice parameters are typically determined with an astonishing accuracy of about \unit[1]{\%}.
+Although LDA is known to overestimate the cohesive energy in solids by \unit[10--20]{\%}, the lattice parameters are typically determined with an astonishing accuracy of about \unit[1]{\%}.
More accurate approximations of the exchange-correlation functional are realized by the introduction of gradient corrections with respect to the density \cite{kohn65}.
Respective considerations are based on the concept of an exchange-correlation hole density describing the depletion of the electron density in the vicinity of an electron.
However, these methods rely on the fact that the wave functions are localized and exhibit an exponential decay resulting in a sparse Hamiltonian.
Another approach is to represent the KS wave functions by plane waves.
-In fact, the employed {\textsc vasp} software is solving the KS equations within a plane-wave (PW) basis set.
+In fact, the employed \textsc{vasp} software is solving the KS equations within a plane-wave (PW) basis set.
The idea is based on the Bloch theorem \cite{bloch29}, which states that in a periodic crystal each electronic wave function $\Phi_i(\vec{r})$ can be written as the product of a wave-like envelope function $\exp(i\vec{kr})$ and a function that has the same periodicity as the lattice.
The latter one can be expressed by a Fourier series, i.e. a discrete set of plane waves whose wave vectors just correspond to reciprocal lattice vectors $\vec{G}$ of the crystal.
Thus, the one-electron wave function $\Phi_i(\vec{r})$ associated with the wave vector $\vec{k}$ can be expanded in terms of a discrete PW basis set
In order to achieve these properties different PPs are required for the different shapes of the orbitals, which are determined by their angular momentum.
Mathematically, a non-local PP, which depends on the angular momentum, has the form
\begin{equation}
-V_{\text{nl}}(\vec{r}) = \sum_{lm} \mid lm \rangle V_l(\vec{r}) \langle lm \mid
+V_{\text{nl}}(\vec{r}) = \sum_{lm} | lm \rangle V_l(\vec{r}) \langle lm |
\text{ .}
\end{equation}
Applying of the operator $V_{\text{nl}}(\vec{r})$ decomposes the electronic wave functions into spherical harmonics $\mid lm \rangle$, i.e. the orbitals with azimuthal angular momentum $l$ and magnetic number $m$, which are then multiplied by the respective pseudopotential $V_l(\vec{r})$ for angular momentum $l$.
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 calculating the respective vibrational modes.
-Modifications used to add the CRT feature to the {\textsc vasp} code and a short instruction on how to use it can be found in appendix \ref{app:patch_vasp}.
+Modifications used to add the CRT feature to the \textsc{vasp} code and a short instruction on how to use it can be found in appendix \ref{app:patch_vasp}.
% todo - advantages of pw basis concenring hf forces
projects / lectures / latex.git / blobdiff