From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Tue, 10 May 2011 15:14:52 +0000 (+0200)
Subject: finished classical simulation methods, go on with dft
X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=561ce1a0f44231abc16b975d4528f6c332a96c86;p=lectures%2Flatex.git

finished classical simulation methods, go on with dft
---

diff --git a/posic/thesis/basics.tex b/posic/thesis/basics.tex
index 0e735ee..92b402c 100644
--- a/posic/thesis/basics.tex
+++ b/posic/thesis/basics.tex
@@ -44,10 +44,15 @@ A potential is necessary describing 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.
 The microscopic information can then be translated to macroscopic observables by means of statistical mechanics.
 
-The basic idea is to integrate Newton's equations numerically.
+The basic idea is to assume that the particles can be described classically by Newton's equations of motion, which are integrated numerically.
 A system of $N$ particles of masses $m_i$ ($i=1,\ldots,N$) at positions ${\bf r}_i$ and velocities $\dot{{\bf r}}_i$ is given by
 \begin{equation}
-m_i \frac{d^2}{dt^2} {\bf r}_i = {\bf F}_i \, \textrm{.}
+%m_i \frac{d^2}{dt^2} {\bf r}_i = {\bf F}_i \Leftrightarrow
+%m_i \frac{d}{dt} {\bf r}_i = {\bf p}_i\textrm{ , } \quad
+%\frac{d}{dt} {\bf p}_i = {\bf F}_i\textrm{ .}
+m_i \ddot{{\bf r}_i} = {\bf F}_i \Leftrightarrow
+m_i \dot{{\bf r}_i} = {\bf p}_i\textrm{, }
+\dot{{\bf p}_i} = {\bf F}_i\textrm{ .}
 \end{equation}
 The forces ${\bf F}_i$ are obtained from the potential energy $U(\{{\bf r}\})$:
 \begin{equation}
@@ -65,50 +70,8 @@ Three ingredients are required for a MD simulation:
 \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}
-Furthermore special techniques will be outlined, which reduce the complexity of the MD algorithm, though the evaluation of energy and force almost inevitably dictates the overall speed.
-
-\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}.
-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}
-\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)+
-\frac{\delta t^3}{6}\vec{b}_i(t) + \mathcal{O}(\delta t^4)
-\label{basics:verlet:taylor1}
-\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)-
-\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.
-The Verlet algorithm is obtained by summarizing and substracting 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)+
-\mathcal{O}(\delta t^4)
-\end{equation}
-\begin{equation}
-\vec{v}_i(t)=\frac{1}{2\delta t}[\vec{r}_i(t+\delta t)-\vec{r}_i(t-\delta t)]+
-\mathcal{O}(\delta t^3)
-\end{equation}
-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$
-\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)]
-\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{ .}
-\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.
+These ingredients will be outlined in the follwoing.
+The discussion is restricted to methods employed within this study.
 
 \subsection{Interaction potentials for silicon and carbon}
 \label{subsection:interact_pot}
@@ -188,7 +151,7 @@ The angular dependence does not give a fixed minimum angle between bonds since t
 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.
-New parameters for the mixed system are $\chi$, which is used to finetune the strength of heteropolar bonds, and $\omeag$, 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.
+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.
 For a three body potential ($V_{ij} \neq V{ji}$) the derivation is of the form
@@ -215,20 +178,101 @@ 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.
 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 provides succeeds in the description of the low as well as high coordinated structures is provided.
+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.
 Defect properties are only fairly reproduced but the description is comparable to previously published potentials.
-It is claimed that the potential enables the modeling of widely different configurations and transitions among these and has recently been used to simulate the inert gas condensation of Si-C nanoparticles \cite{erhart04}.
+It is claimed that the potential enables modeling of widely different configurations and transitions among these and has recently been used to simulate the inert gas condensation of Si-C nanoparticles \cite{erhart04}.
 Therefore the Erhart/Albe (EA) potential is considered the superior analytical bond order potential to study the SiC precipitation and associated processes in Si.
 
+\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}.
+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}
+\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)+
+\frac{\delta t^3}{6}\vec{b}_i(t) + \mathcal{O}(\delta t^4)
+\label{basics:verlet:taylor1}
+\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)-
+\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.
+The Verlet algorithm is obtained by summarizing and substracting 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)+
+\mathcal{O}(\delta t^4)
+\end{equation}
+\begin{equation}
+\vec{v}_i(t)=\frac{1}{2\delta t}[\vec{r}_i(t+\delta t)-\vec{r}_i(t-\delta t)]+
+\mathcal{O}(\delta t^3)
+\end{equation}
+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$
+\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)]
+\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{ .}
+\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.
+
 \subsection{Statistical ensembles}
 \label{subsection:statistical_ensembles}
 
-By default ... NVE
+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.
+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, we need to control T -> NVT
-.. and p -> NpT ...
+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.
+\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{ .}
+\label{eq:basics:ts}
+\end{equation}
+The volume of the synthesized material can hardly be controlled in experiment.
+Instead the pressure can be adjusted.
+Holding constant the pressure in addition to the temperature of the system its states are represented by the isothermal-isobaric $NpT$ ensemble.
+The expression for the pressure of a system derived from the equipartition theorem is given by
+\begin{equation}
+pV=Nk_{\text{B}}T+\langle W\rangle\text{, }W=-\frac{1}{3}\sum_i\vec{r}_i\nabla_{\vec{r}_i}U
+\text{, }
+\label{eq:basics:ps}
+\end{equation}
+where $W$ is the virial and $U$ is the configurational energy.
+
+Berendsen~et~al.~\cite{berendsen84} proposed a method, which is easy to implement, to couple the system to an external bath with constant temperature $T_0$ or pressure $p_0$ with adjustable time constants $\tau_T$ and $\tau_p$ determining the strength of the coupling.
+Control of the respective variable is based on the relations given in equations \eqref{eq:basics:ts} and \eqref{eq:basics:ps}.
+The thermostat is achieved by scaling the velocities of all atoms in every time step $\delta t$ from $\vec{v}_i$ to $\lambda \vec{v}_i$, with
+\begin{equation}
+\lambda=\left[1+\frac{\delta t}{\tau_T}(\frac{T_0}{T}-1)\right]^\frac{1}{2}
+\text{ ,}
+\end{equation}
+where $T$ is the current temperature according to equation \eqref{eq:basics:ts}.
+The barostat adjusts the pressure by changing the virial through scaling of the particle positions $\vec{r}_i$ to $\mu \vec{r}_i$ and the volume $V$ to $\mu^3 V$, with
+\begin{equation}
+\mu=\left[1-\frac{\beta\delta t}{\tau_p}(p_0-p)\right]^\frac{1}{3}\text{ ,}
+\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.
+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.
+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}
 \label{section:dft}