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}.
+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.
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 Fig. \ref{img:tersoff_angle}.
+This is illustrated in Fig.~\ref{img:tersoff_angle}.
\begin{figure}[t]
\begin{center}
\includegraphics[width=8cm]{tersoff_angle.eps}
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$.
+Applying the operator $V_{\text{nl}}(\vec{r})$ decomposes the electronic wave functions into spherical harmonics $| 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$.
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}.
Writing down the derivative of the total energy $E$ with respect to the position $\vec{R}_i$ of ion $i$
\begin{equation}
\frac{dE}{d\vec{R_i}}=
- \sum_j \Phi_j^* \frac{\partial H}{\partial \vec{R}_i} \Phi_j
-+\sum_j \frac{\partial \Phi_j^*}{\partial \vec{R}_i} H \Phi_j
-+\sum_j \Phi_j^* H \frac{\partial \Phi_j}{\partial \vec{R}_i}
+ \sum_j \langle \Phi_j | \frac{\partial H}{\partial\vec{R}_i} | \Phi_j \rangle
++\sum_j \langle \frac{\partial \Phi_j}{\partial\vec{R}_i} | H \Phi_j \rangle
++\sum_j \langle \Phi_j H | \frac{\partial \Phi_j}{\partial \vec{R}_i} \rangle
\text{ ,}
\end{equation}
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 \Phi_j^*\Phi_j\frac{\partial V}{\partial \vec{R}_i}
+\vec{F}_i=-\sum_j \langle \Phi_j | \Phi_j\frac{\partial V}{\partial \vec{R}_i} \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.
\caption[Insertion positions for interstitial defect atoms in the diamond lattice.]{Insertion positions for the tetrahedral ({\color{red}$\bullet$}), hexagonal ({\color{green}$\bullet$}), \hkl<1 0 0> dumbbell ({\color{yellow}$\bullet$}), \hkl<1 1 0> dumbbell ({\color{magenta}$\bullet$}) and bond-centered ({\color{cyan}$\bullet$}) interstitial defect atom in the diamond lattice. The black dots correspond to the lattice atoms and the blue lines indicate the covalent bonds of the perfect diamond structure.}
\label{fig:basics:ins_pos}
\end{figure}
-The respective estimated interstitial insertion positions for various interstitial structures in a diamond lattice are displayed in Fig. \ref{fig:basics:ins_pos}.
+The respective estimated interstitial insertion positions for various interstitial structures in a diamond lattice are displayed in Fig.~\ref{fig:basics:ins_pos}.
The labels of the interstitial types indicate their positions in the interstitial lattice.
In a dumbbell (DB) configuration two atoms share a single lattice site along a certain direction that is also comprehended in the label of the defect.
For the DB configurations the nearest atom of the bulk lattice is slightly displaced by $(0,0,-1/8)$ and $(-1/8,-1/8,0)$ of the unit cell length respectively.
\end{figure}
One possibility to compute the migration path from one stable configuration into another one is provided by the constrained relaxation technique (CRT) \cite{kaukonen98}.
The atoms involving great structural changes in the diffusion process are moved stepwise from the starting to the final position and relaxation after each step is only allowed in the plane perpendicular to the direction of the vector connecting its starting and final position.
-This is illustrated in Fig. \ref{fig:basics:crto}.
+This is illustrated in Fig.~\ref{fig:basics:crto}.
The number of steps required for smooth transitions depends on the shape of the potential energy surface.
No constraints are applied to the remaining atoms to allow for the relaxation of the surrounding lattice.
To prevent the remaining lattice to shift according to the displacement of the defect, however, some atoms far away from the defect region should be fixed in all three coordinate directions.
For some structures even the expected final configurations are not obtained.
Thus, the method mentioned above is adjusted adding further constraints in order to obtain smooth transitions with respect to energy and structure.
In the modified method all atoms are stepwise displaced towards their final positions.
-In addition to this, relaxation of each individual atom is only allowed in the plane perpendicular to the last individual displacement vector, as displayed in Fig. \ref{fig:basics:crtm}.
+In addition to this, relaxation of each individual atom is only allowed in the plane perpendicular to the last individual displacement vector, as displayed in Fig.~\ref{fig:basics:crtm}.
In the modified version respective energies could be higher than the real ones due to the additional constraints that prevent a more adequate relaxation until the final configuration is reached.
-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.
+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}.
projects / lectures / latex.git / blobdiff