\chapter{Basic principles of utilized simulation techniques}
\label{chapter:basics}
-In the following the simulation methods used within the scope of this study are introduced.
+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}.
-An overview of the most important tools within the MD package is given in appendix \ref{app:code}.
+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}.
+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.
+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 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}.
+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}.
\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.
-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.
+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, which is 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, i.e.\ 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.
\subsection{Introduction to molecular dynamics simulations}
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}.
+ Interaction potentials and their accuracy for describing certain systems of elements will be outlined in section~\ref{subsection:interact_pot}.
\item An integrator is needed, which propagates the particle positions and velocities from time $t$ to $t+\delta t$, realized by a finite difference scheme which moves trajectories discretely in time.
- This is explained in section \ref{subsection:integrate_algo}.
+ This is explained in section~\ref{subsection:integrate_algo}.
\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}.
+ This is discussed in section~\ref{subsection:statistical_ensembles}.
\end{enumerate}
These ingredients will be outlined in the following.
The discussion is restricted to methods employed within this study.
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 increase 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.
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 a cutoff function to limit the range of interaction to nearest neighbors.
+The function $f_C$ is a cut-off 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\{
\label{eq:basics:fc}
\end{equation}
As discussed above, $b_{ij}$ represents a measure of the bond order, monotonously decreasing with the coordination of atoms $i$ and $j$.
-It is of the form:
+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}) \\
This is illustrated in Fig.~\ref{img:tersoff_angle}.
\begin{figure}[t]
\begin{center}
-\includegraphics[width=8cm]{tersoff_angle.eps}
+\includegraphics[width=6cm]{tersoff_angle.eps}
\end{center}
\caption{Angle between bonds of atoms $i,j$ and $i,k$.}
\label{img:tersoff_angle}
The angular dependence does not give a fixed minimum angle between bonds since the expression is embedded inside the bond order term.
The relation to the above discussed bond order potential becomes obvious if $\chi=1, \beta=1, n=1, \omega=1$ and $c=0$.
Parameters with a single subscript correspond to the parameters of the elemental system~\cite{tersoff_si3,tersoff_c} while the mixed parameters are obtained by interpolation from the elemental parameters by the arithmetic or geometric mean.
-The elemental parameters were obtained by fit with respect to the cohesive energies of real and hypothetical bulk structures and the bulk modulus and bond length of the diamond structure.
+The elemental parameters were obtained by a fit with respect to the cohesive energies of real and hypothetical bulk structures and the bulk modulus and bond length of the diamond structure.
New parameters for the mixed system are $\chi$, which is used to finetune the strength of heteropolar bonds, and $\omega$, which is set to one for the C-Si interaction but is available as a feature to permit the application of the potential to more drastically different types of atoms in the future.
The force acting on atom $i$ is given by the derivative of the potential energy.
\begin{equation}
F^i = - \nabla_{{\bf r}_i} E \textrm{ .}
\end{equation}
-Details of the Tersoff potential derivative are presented in appendix \ref{app:d_tersoff}.
+Details of the Tersoff potential derivative are presented in appendix~\ref{app:d_tersoff}.
\subsubsection{Improved analytical bond order potential}
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}.
+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.
However, T3, which is used in the Si/C potential, suffers from an underestimation of the dimer binding energy.
Similar behavior is found for the C-C interaction.
For this reason, Erhart and Albe provide a reparametrization of the Tersoff potential based on three independently fitted parameter sets for the Si-Si, C-C and Si-C interaction~\cite{albe_sic_pot}.
The functional form is similar to the one proposed by Tersoff.
-Differences in the energy functional and the force evaluation routine are pointed out in appendix \ref{app:d_tersoff}.
-Concerning Si the elastic properties of the diamond phase as well as the structure and energetics of the dimer are reproduced very well.
+Differences in the energy functional and the force evaluation routine are pointed out in appendix~\ref{app:d_tersoff}.
+Concerning Si, the elastic properties of the diamond phase as well as the structure and energetics of the dimer are reproduced very well.
The new parameter set for the C-C interaction yields improved dimer properties while at the same time delivers a description of the bulk phase similar to the Tersoff potential.
The potential succeeds in the description of the low as well as high coordinated structures.
The description of elastic properties of SiC is improved with respect to the potentials available in literature.
\subsection{Statistical ensembles}
\label{subsection:statistical_ensembles}
-Using the above mentioned algorithms the most basic type of MD is realized by simply integrating the equations of motion of a fixed number of particles ($N$) in a closed volume $V$ realized by periodic boundary conditions (PBC).
+Using the above mentioned algorithms, the most basic type of MD is realized by simply integrating the equations of motion of a fixed number of particles ($N$) in a closed volume $V$ realized by periodic boundary conditions (PBC).
Providing a stable integration algorithm the total energy $E$, i.e.\ the kinetic and configurational energy of the particles, is conserved.
This is known as the $NVE$, or microcanonical ensemble, describing an isolated system composed of microstates, among which the number of particles, volume and energy are held constant.
\end{equation}
where $\beta$ is the isothermal compressibility and $p$ corresponds to the current pressure, which is determined by equation \eqref{eq:basics:ps}.
-Using this method the system does not behave like a true $NpT$ ensemble.
+Using this method, the system does not behave like a true $NpT$ ensemble.
On average $T$ and $p$ correspond to the expected values.
For large enough time constants, i.e.\ $\tau > 100 \delta t$, the method shows realistic fluctuations in $T$ and $p$.
The advantage of the approach is that the coupling can be decreased to minimize the disturbance of the system and likewise be adjusted to suit the needs of a given application.
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.
In contrast, instead of using the description by the many-body wave function, the key point in density functional theory (DFT) is to recast the problem to a description utilizing the charge density $n(\vec{r})$, which constitutes a quantity in real space depending only on the three spatial coordinates.
-In the following sections the basic idea of DFT will be outlined.
+In the following sections, the basic idea of DFT will be outlined.
As will be shown, DFT can formally be regarded as an exactification of the Thomas Fermi theory~\cite{thomas27,fermi27} and the self-consistent Hartree equations~\cite{hartree28}.
A nice review is given in the Nobel lecture of Kohn~\cite{kohn99}, one of the inventors of DFT.
The challenging problem of determining the exact ground-state is now formally reduced to the determination of the $3$-dimensional function $n(\vec{r})$, which minimizes the energy functional.
However, the complexity associated with the many-electron problem is now relocated in the task of finding the well-defined but, in contrast to the potential energy, not explicitly known functional $F[n(\vec{r})]$.
-It is worth to note, that this minimal principle may be regarded as exactification of the TF theory, which is rederived by the approximations
+It is worth to note that this minimal principle may be regarded as exactification of the TF theory, which is rederived by the approximations
\begin{equation}
T=\int n(\vec{r})\frac{3}{10}k_{\text{F}}^2(n(\vec{r}))d\vec{r}
\text{ ,}
\subsection{Pseudopotentials}
As discussed in the last part of the previous section, an extremely large basis set of plane waves would be required to perform an all-electron calculation and a vast amount of computational time would be required to calculate the electronic wave functions.
-It is worth to stress out one more time, that this is mainly due to the orthogonalization wiggles of the wave functions of valence electrons near the nuclei.
+It is worth to stress out one more time that this is mainly due to the orthogonalization wiggles of the wave functions of valence electrons near the nuclei.
Thus, existing core states practically prevent the use of a PW basis set.
However, the core electrons, which are tightly bound to the nuclei, do not contribute significantly to chemical bonding or other physical properties of the solid.
This fact is exploited in the pseudopotential (PP) approach~\cite{cohen70} by removing the core electrons and replacing the atom and the associated strong ionic potential by a pseudoatom and a weaker PP that acts on a set of pseudo wave functions rather than the true valance wave functions.
The standard generation procedure of pseudopotentials proceeds by varying its parameters until the pseudo eigenvalues are equal to the all-electron valence eigenvalues and the pseudo wave functions match the all-electron valence wave functions beyond a certain cut-off radius determining the core region.
Modified methods to generate ultra-soft pseudopotentials were proposed, which address the rapid convergence with respect to the size of the plane wave basis set~\cite{vanderbilt90,troullier91}.
-Using PPs the rapid oscillations of the wave functions near the core of the atoms are removed considerably reducing the number of plane waves necessary to appropriately expand the wave functions.
+Using PPs, the rapid oscillations of the wave functions near the core of the atoms are removed considerably reducing the number of plane waves necessary to appropriately expand the wave functions.
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}
\subsection{Structural relaxation and Hellmann-Feynman theorem}
-Up to this point, the system is in the ground state with respect to the electronic subsystem, while the positions of the ions as well as size and shape of the supercell are fixed.
+Up to this point, the system is in the ground state with respect to the electronic subsystem while the positions of the ions as well as size and shape of the supercell are fixed.
To investigate equilibrium structures, however, the ionic subsystem must also be allowed to relax into a minimum energy configuration.
Local minimum configurations can be easily obtained in a MD-like way by moving the nuclei over small distances along the directions of the forces, as discussed in the MD chapter above.
Clearly, the conjugate gradient method constitutes a more sophisticated scheme, which will locate the equilibrium positions of the ions more rapidly.
indeed reveals a contribution to the change in total energy due to the change of the wave functions $\Phi_j$.
However, provided that the $\Phi_j$ are eigenstates of $H$, it is easy to show that the last two terms cancel each other and in the special case of $H=T+V$ the force is given by
\begin{equation}
-\vec{F}_i=-\sum_j \langle \Phi_j | \Phi_j\frac{\partial V}{\partial \vec{R}_i} \rangle
+\vec{F}_i=-\sum_j \langle \Phi_j | \frac{\partial V}{\partial \vec{R}_i} \Phi_j \rangle
\text{ .}
\end{equation}
This is called the Hellmann-Feynman theorem~\cite{feynman39}, which enables the calculation of forces, called the Hellmann-Feynman forces, acting on the nuclei for a given configuration, without the need for evaluating computationally costly energy maps.
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.
+Therefor, a supercell containing the perfect crystal is generated in an initial process.
If not by construction, the system should be fully relaxed.
The substitutional or vacancy defect is realized by replacing or removing one atom contained in the supercell.
Interstitial defects are created by adding an atom at positions located in the space between regular lattice sites.
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