\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.
-This understanding may be regarded as the basic principle of molecular dynamics, considering an isolated system of particles, the behaviour of which is fully determined by the solution of the classical equations of motion.
+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}
-Molecular dynamics (MD) simulation is a technique to compute a system of particles, referred to as molecules, with their positions, volocities and forces among each other evolving in time.
+Molecular dynamics (MD) simulation is a technique to compute a system of particles, referred to as molecules, with their positions, velocities and forces among each other evolving in time.
The MD method was first introduced by Alder and Wainwright in 1957 \cite{alder57,alder59} to study the interactions of hard spheres.
-The basis of the approach are Newton's equations of motion to describe classicaly the many-body system.
+The basis of the approach are Newton's equations of motion to describe classically the many-body system.
MD is the numerical way of solving the $N$-body problem which cannot be solved analytically for $N>3$.
A potential is necessary to describe the interaction of the particles.
By MD, a complete description of the system in the sense of classical mechanics on the microscopic level is obtained.
\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 An integrator is needed, which propagtes the particle positions and velocities from time $t$ to $t+\delta t$, realised by a finite difference scheme which moves trajectories discretely in time.
+\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}.
\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}
-These ingredients will be outlined in the follwoing.
+These ingredients will be outlined in the following.
The discussion is restricted to methods employed within this study.
\subsection{Interaction potentials for silicon and carbon}
$U_1$ is a single particle potential describing external forces.
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.
+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.
-Silicon and carbon for instance, have a diamand and zincblende structure with four covalently bonded neighbors, which is far from a close-packed structure.
+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.
\subsubsection{The Tersoff bond order potential}
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.
-Based on pseudopotential theory the bond order term $b_{ijk}$ limitting the attractive pair interaction is of the form $b_{ijk}\propto Z^{-\delta}$ where $Z$ is the coordination number and $\delta$ a constant \cite{abell85}, which is $\frac{1}{2}$ in the seond-moment approximation within the tight binding scheme \cite{horsfield96}.
+Based on pseudopotential theory, the bond order term $b_{ijk}$ limiting the attractive pair interaction is of the form $b_{ijk}\propto Z^{-\delta}$ where $Z$ is the coordination number and $\delta$ a constant \cite{abell85}, which is $\frac{1}{2}$ in the second-moment approximation within the tight binding scheme \cite{horsfield96}.
Tersoff incorporated the concept of bond order in a three-body potential formalism.
The interatomic potential is taken to have the form
\label{subsection:integrate_algo}
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.
+The idea of the so-called Verlet and a variant, the velocity Verlet algorithm, which additionally 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)=
\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}
+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)=
2\vec{r}_i(t)-\vec{r}_i(t-\delta t)+\frac{\delta t^2}{m_i}\vec{f}_i(t)+
\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).
-Providing a stable integration algorithm the total energy $E$, i.e. the kinetic and configurational energy of the paticles, is conserved.
+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.
However, the successful formation of SiC dictates precise control of temperature by external heating.
While the temperature of such a system is well defined, the energy is no longer conserved.
The microscopic states of a system, which is in thermal equilibrium with an external thermal heat bath, are represented by the $NVT$ ensemble.
-In the so-called canonical ensemble the temperature $T$ is related to the expactation value of the kinetic energy of the particles, i.e.
+In the so-called canonical ensemble the temperature $T$ is related to the expectation value of the kinetic energy of the particles, i.e.
\begin{equation}
\langle E_{\text{kin}}\rangle = \frac{3}{2}Nk_{\text{B}}T \text{, }
E_{\text{kin}}=\sum_i \frac{\vec{p}^2_i}{2m_i} \text{ .}
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.
It provides a stable algorithm that allows smooth changes of the system to new values of temperature or pressure, which is ideal for the investigated problem.
-\section{Denstiy functional theory}
+\section{Density functional theory}
\label{section:dft}
-Dirac declared that chemistry has come to an end, its content being entirely contained in the powerul equation published by Schr\"odinger in 1926 \cite{schroedinger26} marking the beginning of wave mechanics.
+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.
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 determinats to be evaluated massively increases.
+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.
\text{ ,}
\end{equation}
where $Z_l$ are the atomic numbers of the nuclei and $\Psi$ is the many-electron wave function, which depends on the positions and spins of the electrons.
-Accordingly, there is only a parametrical dependence on the ionic coordinates $\vec{R}_l$.
+Accordingly, there is only a parametric dependence on the ionic coordinates $\vec{R}_l$.
However, the remaining number of free parameters is still too high and need to be further decreased.
\subsection{Hohenberg-Kohn theorem and variational principle}
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.
-Considering a system with a nondegenerate ground state there is obviously only one ground-state charge density $n_0(\vec{r})$ that correpsonds to a given potential $V(\vec{r})$.
+Considering a system with a nondegenerate ground state there is obviously only one ground-state charge density $n_0(\vec{r})$ that corresponds to a given potential $V(\vec{r})$.
In 1964 Hohenberg and Kohn showed the opposite and far less obvious result \cite{hohenberg64}.
Employing no more than the Rayleigh-Ritz minimal principle it is concluded by {\em reductio ad absurdum} that for a nondegenerate ground state the same charge density cannot be generated by different potentials.
Thus, the charge density of the ground state $n_0(\vec{r})$ uniquely determines the potential $V(\vec{r})$ and, consequently, the full Hamiltonian and ground-state energy $E_0$.
In mathematical terms the full many-electron ground state is a unique functional of the charge density.
-Im particular, $E$ is a functional $E[n(\vec{r})]$ of $n(\vec{r})$.
+In particular, $E$ is a functional $E[n(\vec{r})]$ of $n(\vec{r})$.
The ground-state charge density $n_0(\vec{r})$ minimizes the energy functional $E[n(\vec{r})]$, its value corresponding to the ground-state energy $E_0$, which enables the formulation of a minimal principle in terms of the charge density \cite{hohenberg64,levy82}
\begin{equation}
%\end{equation}
The system of interacting electrons is mapped to an auxiliary system, the Kohn-Sham (KS) system, of non-interacting electrons in an effective potential.
-The exact effective potential $V_{\text{eff}}(\vec{r})$ may be regarded as a fictious external potential yielding a gound-state density for non-interacting electrons, which is equal to that for interacting electrons in the external potential $V(\vec{r})$.
+The exact effective potential $V_{\text{eff}}(\vec{r})$ may be regarded as a fictitious external potential yielding a ground-state density for non-interacting electrons, which is equal to that for interacting electrons in the external potential $V(\vec{r})$.
The one-electron KS orbitals $\Phi_i(\vec{r})$ as well as the respective KS energies $\epsilon_i$ are not directly attached to any physical observable except for the ground-state density, which is determined by equation \eqref{eq:basics:kse3} and the ionization energy, which is equal to the highest occupied relative to the vacuum level.
The KS equations may be considered the formal exactification of the Hartree theory, which it is reduced to if the exchange-correlation potential and functional are neglected.
In addition to the Hartree-Fock (HF) method, KS theory includes the difference of the kinetic energy of interacting and non-interacting electrons as well as the remaining contributions to the correlation energy that are not part of the HF correlation.
\end{equation}
which is, thus, called local density approximation (LDA).
Here, the exchange-correlation energy per particle of the uniform electron gas of constant density $n$ is used for $\epsilon_{\text{xc}}(n(\vec{r}))$.
-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 paramterizations exist \cite{perdew81}.
+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]{\%}.
Finally, a set of basis functions is required to represent the one-electron KS wave functions.
With respect to the numerical treatment it is favorable to approximate the wave functions by linear combinations of a finite number of such basis functions.
-Covergence of the basis set, i.e. convergence of the wave functions with respect to the amount of basis functions, is most crucial for the accuracy of the numerical calulations.
+Convergence of the basis set, i.e. convergence of the wave functions with respect to the amount of basis functions, is most crucial for the accuracy of the numerical calculations.
Two classes of basis sets, the plane-wave and local basis sets, exist.
Local basis set functions usually are atomic orbitals, i.e. mathematical functions that describe the wave-like behavior of electrons, which are localized, i.e. centered on atoms or bonds.
The basis set is orthonormal by construction and, as mentioned above, it is simple to check for convergence.
The biggest advantage, however, is the ability to perform exact calculations by a discrete sum over a numerical grid.
This is due to the related construction of the grid and the PW basis.
-Ofcourse, exactness is restricted by the fact that the PW basis set is finite.
+Of course, exactness is restricted by the fact that the PW basis set is finite.
The simple form of the PW representation of the KS equations
\begin{equation}
\sum_{\vec{G}'} \left[
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.
This way, the pseudo wave functions become smooth near the nuclei.
-Most PPs statisfy four general conditions.
+Most PPs satisfy four general conditions.
The pseudo wave functions generated by the PP should not contain nodes, i.e. the pseudo wave functions should be smooth and free of wiggles in the core region.
Outside the core region, the pseudo and real valence wave functions as well as the generated charge densities need to be identical.
The charge enclosed within the core region must be equal for both wave functions.
V_{\text{nl}}(\vec{r}) = \sum_{lm} \mid lm \rangle V_l(\vec{r}) \langle lm \mid
\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 psuedopotential $V_l(\vec{r})$ for angular momentum $l$.
-The standard generation procedure of pseudopotentials proceeds by varying its parameters until the pseudo eigenvalues are eual to the all-electron valence eigenvalues and the pseudo wave functions match the all-electron valence wave functions beyond a certain cut-off radius detrmining the core region.
+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$.
+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.
Since the values of the wave function within a small interval around $\vec{k}$ are almost identical, it is possible to approximate the infinite sum by a sum over an affordable number of $k$ points, each representing the respective region of the wave function in $\vec{k}$ space.
Methods have been derived for obtaining very accurate approximations by a summation over special sets of $\vec{k}$ points with distinct, associated weights \cite{baldereschi73,chadi73,monkhorst76}.
If present, symmetries in reciprocal space may further reduce the number of calculations.
-For supercells, i.e. repeating unit cells that contain several primitive cells, restricting the sampling of the Brillouin zone (BZ) to the $\Gamma$ point can yield quite accurat results.
+For supercells, i.e. repeating unit cells that contain several primitive cells, restricting the sampling of the Brillouin zone (BZ) to the $\Gamma$ point can yield quite accurate results.
In fact, with respect to BZ sampling, calculating wave functions of a supercell containing $n$ primitive cells for only one $\vec{k}$ point is equivalent to the scenario of a single primitive cell and the summation over $n$ points in $\vec{k}$ space.
-In general, finer $\vec{k}$ point meshes better account for the periodicity of a system, which in some cases, however, might be fictious anyway.
+In general, finer $\vec{k}$ point meshes better account for the periodicity of a system, which in some cases, however, might be fictitious anyway.
\subsection{Structural relaxation and Hellmann-Feynman theorem}
+\sum_j \Phi_j^* H \frac{\partial \Phi_j}{\partial \vec{R}_i}
\text{ ,}
\end{equation}
-indeed reveals a contributon to the chnage in total energy due to the change of the wave functions $\Phi_j$.
+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}
\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}
\end{figure}
-One possibility to compute the migration path from one stable cofiguration into another one is provided by the constrained relaxation technique (CRT) \cite{kaukonen98}.
+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}.
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, ohowever, some atoms far away from the defect region should be fixed in all three coordinate directions.
+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.
However, for the present study, the method tremendously failed.
-Abrupt changes in structure and configurational energy occured among relaxed structures of two successive displacement steps.
+Abrupt changes in structure and configurational energy occurred among relaxed structures of two successive displacement steps.
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 repsect to energy and structure.
+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 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 copnfiguration is reached.
+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.
-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.
+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}.
% todo - advantages of pw basis concenring hf forces
projects / lectures / latex.git / blobdiff