\chapter{Force evaluation for the three body Tersoff potential}
\label{app:d_tersoff}

  \section{Form of the Tersoff potential and its derivative}

The Tersoff potential is of the form
\begin{eqnarray}
E & = & \sum_i E_i = \frac{1}{2} \sum_{i \ne j} V_{ij} \textrm{ ,} \\
V_{ij} & = & f_C(r_{ij}) [ f_R(r_{ij}) + b_{ij} f_A(r_{ij}) ] \textrm{ .}
\end{eqnarray}
The repulsive $f_R$ and attractive $f_A$ part is given by
\begin{eqnarray}
f_R(r_{ij}) & = & A_{ij} \exp (- \lambda_{ij} r_{ij} ) \textrm{ ,} \\
f_A(r_{ij}) & = & -B_{ij} \exp (- \mu_{ij} r_{ij} ) \textrm{ .}
\end{eqnarray}
The bond order function $b_{ij}$ is
\begin{equation}
b_{ij} = \chi_{ij} (1 + \beta_i^{n_i} \zeta^{n_i}_{ij})^{-1/2n_i}
\end{equation}
with
\begin{eqnarray}
\zeta_{ij} & = & \sum_{k \ne i,j} f_C (r_{ik}) \omega_{ik} g(\theta_{ijk}) \textrm{ ,}\\
g(\theta_{ijk}) & = & 1 + c_i^2/d_i^2 - c_i^2/[d_i^2 + (h_i - \cos \theta_{ijk})^2] \textrm{ .}
\end{eqnarray}
The cutoff function $f_C$ is taken to be
\begin{equation}
f_C(r_{ij}) = \left\{
  \begin{array}{ll}
    1, & r_{ij} < R_{ij} \\
    \frac{1}{2} + \frac{1}{2} \cos \Big[ \pi (r_{ij} - R_{ij})/(S_{ij} - R_{ij}) \Big], &
 R_{ij} < r_{ij} < S_{ij} \\
    0, & r_{ij} > S_{ij}
  \end{array} \right.
\end{equation}
with $\theta_{ijk}$ being the bond angle between bonds $ij$ and $ik$ as shown in Figure \ref{img:tersoff_angle}.\\
\\
For a three body potential, if $V_{ij}$ is not equal to $V_{ji}$, the derivative is of the form
\begin{equation}
\nabla_{{\bf r}_i} E = \frac{1}{2} \big[ \sum_j ( \nabla_{{\bf r}_i} V_{ij} + \nabla_{{\bf r}_i} V_{ji} ) + \sum_k \sum_j \nabla_{{\bf r}_i} V_{jk} \big] \textrm{ .}
\end{equation}
In the following all the necessary derivatives to calculate $\nabla_{{\bf r}_i} E$ are done.

  \section{Derivative of $V_{ij}$ with respect to ${\bf r}_i$}

\begin{eqnarray}
\nabla_{{\bf r}_i} V_{ij} & = & \nabla_{{\bf r}_i} f_C(r_{ij}) \big[ f_R(r_{ij}) + b_{ij} f_A(r_{ij}) \big] + \nonumber \\
 & & + f_C(r_{ij}) \big[ \nabla_{{\bf r}_i} f_R(r_{ij}) + b_{ij} \nabla_{{\bf r}_i} f_A(r_{ij}) + f_A(r_{ij}) \nabla_{{\bf r}_i} b_{ij} \big]
\end{eqnarray}
\begin{eqnarray}
\nabla_{{\bf r}_i} f_R(r_{ij}) & = & A_{ij} \lambda_{ij} \frac{{\bf r}_{ij}}{r_{ij}} \exp(-\lambda_{ij} r_{ij}) \\
\nabla_{{\bf r}_i} f_A(r_{ij}) & = & - B_{ij} \mu_{ij} \frac{{\bf r}_{ij}}{r_{ij}} \exp(-\mu_{ij} r_{ij})
\end{eqnarray}
\begin{equation}
\nabla_{{\bf r}_i} f_C(r_{ij}) = \left\{
  \begin{array}{ll}
    \frac{1}{2} \sin \Big( \frac{\pi(r_{ij}-R_{ij})}{S_{ij}-R_{ij}} \Big) \frac{\pi}{S_{ij}-R_{ij}} \frac{{\bf r}_{ij}}{r_{ij}}, & R_{ij} < r_{ij} < S_{ij} \\
    0, & \textrm{else.}
  \end{array} \right.
\end{equation}
\begin{eqnarray}
\nabla_{{\bf r}_i} b_{ij} &= & - \frac{\chi_{ij}}{2} (1+\beta^{n_i} \zeta_{ij}^{n_i})^{-\frac{1}{2n_i}-1} \beta^{n_i} \zeta_{ij}^{n_i-1} \nabla_{{\bf r}_i} \zeta_{ij} \\
\nabla_{{\bf r}_i} \zeta_{ij} & = & \sum_{k \neq i,j} \big( g(\theta_{ijk}) \nabla_{{\bf r}_i} f_C(r_{ik}) + f_C(r_{ik}) \nabla_{{\bf r}_i} g(\theta_{ijk}) \big) \\
\nabla_{{\bf r}_i} g(\theta_{ijk}) & = & - \frac{2(h_i-\cos\theta_{ijk})c_i^2}{\big[d_i^2 + (h_i - \cos\theta_{ijk})^2\big]^2} \nabla_{{\bf r}_i} (\cos\theta_{ijk})
\end{eqnarray}
\begin{eqnarray}
\nabla_{{\bf r}_i} \cos \theta_{ijk} & = & \nabla_{{\bf r}_i} \Big( \frac{{\bf r}_{ij} {\bf r}_{ik}}{r_{ij} r_{ik}} \Big) \nonumber \\
 & = & \Big[ \frac{\cos\theta_{ijk}}{r_{ij}^2} - \frac{1}{r_{ij} r_{ik}} \Big] {\bf r}_{ij} + \Big[ \frac{\cos\theta_{ijk}}{r_{ik}^2} - \frac{1}{r_{ij} r_{ik}} \Big] {\bf r}_{ik}
\end{eqnarray}

  \section{Derivative of $V_{ji}$ with respect to ${\bf r}_i$}

\begin{eqnarray}
\nabla_{{\bf r}_i} V_{ji} & = & \nabla_{{\bf r}_i} f_C(r_{ji}) \big[ f_R(r_{ji}) + b_{ji} f_A(r_{ji}) \big] + \nonumber \\
 & & + f_C(r_{ji}) \big[ \nabla_{{\bf r}_i} f_R(r_{ji}) + b_{ji} \nabla_{{\bf r}_i} f_A(r_{ji}) + f_A(r_{ji}) \nabla_{{\bf r}_i} b_{ji} \big]
\end{eqnarray}
\begin{eqnarray}
\nabla_{{\bf r}_i} f_R(r_{ji}) & = & - A_{ji} \lambda_{ji} \frac{{\bf r}_{ji}}{r_{ji}} \exp(-\lambda_{ji} r_{ji}) = \nabla_{{\bf r}_i} f_R(r_{ij}) \\
\nabla_{{\bf r}_i} f_A(r_{ji}) & = & + B_{ji} \mu_{ji} \frac{{\bf r}_{ji}}{r_{ji}} \exp(-\mu_{ji} r_{ji}) = \nabla_{{\bf r}_i} f_A(r_{ij})
\end{eqnarray}
\begin{equation}
\nabla_{{\bf r}_i} f_C(r_{ij}) = \nabla_{{\bf r}_i} f_C(r_{ij}) = \left\{
  \begin{array}{ll}
    - \frac{1}{2} \sin \Big( \frac{\pi(r_{ji}-R_{ji})}{S_{ji}-R_{ji}} \Big) \frac{\pi}{S_{ji}-R_{ji}} \frac{{\bf r}_{ji}}{r_{ji}}, & R_{ji} < r_{ji} < S_{ji} \\
    0, & \textrm{else.}
  \end{array} \right.
\end{equation}
\begin{eqnarray}
\nabla_{{\bf r}_i} b_{ji} &= & - \frac{\chi_{ji}}{2} (1+\beta^{n_j} \zeta_{ji}^{n_j})^{-\frac{1}{2n_j}-1} \beta^{n_j} \zeta_{ji}^{n_j-1} \nabla_{{\bf r}_i} \zeta_{ji} \\
\nabla_{{\bf r}_i} \zeta_{ji} & = & \sum_{k \neq j,i} \big( g(\theta_{jik}) \nabla_{{\bf r}_i} f_C(r_{jk}) + f_C(r_{jk}) \nabla_{{\bf r}_i} g(\theta_{jik}) \big) \nonumber \\
 & = & \sum_{k \neq j,i} f_C(r_{jk}) \nabla_{{\bf r}_i} g(\theta_{jik}) \quad \Big(\textrm{Reason: }\nabla_{{\bf r}_i} f_C(r_{jk}) = 0\Big) \\
\nabla_{{\bf r}_i} g(\theta_{jik}) & = & - \frac{2(h_j-\cos\theta_{jik})c_j^2}{\big[d_j^2 + (h_j - \cos\theta_{jik})^2\big]^2} \nabla_{{\bf r}_i} (\cos\theta_{jik})
\end{eqnarray}
\begin{eqnarray}
\nabla_{{\bf r}_i} \cos \theta_{jik} & = & \nabla_{{\bf r}_i} \Big( \frac{{\bf r}_{ji} {\bf r}_{jk}}{r_{ji} r_{jk}} \Big) \nonumber \\
 & = & \frac{1}{r_{ji} r_{jk}} {\bf r}_{jk} - \frac{\cos\theta_{jik}}{r_{ji}^2} {\bf r}_{ji}
\end{eqnarray}

  \section{Derivative of $V_{jk}$ with respect to ${\bf r}_i$}

\begin{eqnarray}
\nabla_{{\bf r}_i} V_{jk} & = & f_C(r_{jk}) f_A(r_{jk}) \nabla_{{\bf r}_i} b_{jk} \\
\nabla_{{\bf r}_i} b_{jk} & = & - \frac{\chi_{jk}}{2} (1+\beta^{n_j} \zeta_{jk}^{n_j})^{-\frac{1}{2n_j}-1} \beta^{n_j} \zeta_{jk}^{n_j-1} \nabla_{{\bf r}_i} \zeta_{jk} \\
\nabla_{{\bf r}_i} \zeta_{jk} & = & \sum_{l \neq j,k} \big( g(\theta_{jkl}) \nabla_{{\bf r}_i} f_C(r_{jl}) + f_C(r_{jl}) \nabla_{{\bf r}_i} g(\theta_{jkl}) \big) \nonumber \\
 & = & f_C(r_{ji}) \nabla_{{\bf r}_i} g(\theta_{jki}) + g(\theta_{jki}) \nabla_{{\bf r}_i} f_C(r_{ji}) \\
\nabla_{{\bf r}_i} g(\theta_{jki}) & = & - \frac{2(h_j-\cos\theta_{jki})c_j^2}{\big[d_j^2 + (h_j - \cos\theta_{jki})^2\big]^2} \nabla_{{\bf r}_i} (\cos\theta_{jki}) \\
\nabla_{{\bf r}_i} \cos \theta_{jki} & = & \nabla_{{\bf r}_i} \Big( \frac{{\bf r}_{jk} {\bf r}_{ji}}{r_{jk} r_{ji}} \Big) \nonumber \\
 & = & \frac{1}{r_{jk} r_{ji}} {\bf r}_{jk} - \frac{\cos\theta_{jki}}{r_{ji}^2} {\bf r}_{ji}
\end{eqnarray}

  \section{Implementation issues}

As seen in the last sections the derivatives of $V_{ij}$, $V_{ji}$ and $V_{jk}$
with respect to ${\bf r}_i$ are necessary to compute the forces for atom $i$.
According to this, for every triple $(ijk)$ the derivatives of the three
potential contributions, denoted by $V_{ijk}$, $V_{jik}$ and $V_{jki}$
have to be computed.
In simulation, however, it is not practical to evaluate all three potential
derivatives for each $(ijk)$ triple.
The $V_{jik}$ and $V_{jki}$ potential and its derivatives will be calculated
in subsequent loops anyways.
To avoid multiple computation of the same potential derivatives
the force contributions for atom $j$ and $k$ due to the $V_{ijk}$ contribution
have to be considered by calculating the derivatives of $V_{ijk}$
with respect to ${\bf r}_j$ and ${\bf r}_k$
inside the loop of the $(ijk)$ triple
in addition to the derivative with respect to ${\bf r}_i$.
This poses a more convenient method to obtain the forces
keeping in mind that all the necessary force contributions for atom $i$
are calculated and added in subsequent loops.

The following symmetry considerations help to obtain the 

  \subsection{Derivative of $V_{ij}$ with respect to ${\bf r}_j$}

\begin{eqnarray}
 \nabla_{{\bf r}_j} V_{ij} & = &
 \nabla_{{\bf r}_j} f_C(r_{ij}) \big[ f_R(r_{ij}) +
 b_{ij} f_A(r_{ij}) \big] + \nonumber \\
 & & + f_C(r_{ij}) \big[ \nabla_{{\bf r}_j} f_R(r_{ij}) +
 b_{ij} \nabla_{{\bf r}_j} f_A(r_{ij}) +
 f_A(r_{ij}) \nabla_{{\bf r}_j} b_{ij} \big]
\end{eqnarray}
Using the equality $\nabla_{{\bf r}_i} r_{ij}=-\nabla_{{\bf r}_j} r_{ij}$
the following relations are valid:
\begin{eqnarray}
 \nabla_{{\bf r}_j} f_R(r_{ij}) &=& - \nabla_{{\bf r}_i} f_R(r_{ij}) \\
 \nabla_{{\bf r}_j} f_A(r_{ij}) &=& - \nabla_{{\bf r}_i} f_A(r_{ij}) \\
 \nabla_{{\bf r}_j} f_C(r_{ij}) &=& - \nabla_{{\bf r}_i} f_C(r_{ij})
\end{eqnarray}
The pair contributions .... easy.
Now having a look at $b_{ij}$.
\begin{eqnarray}
 \nabla_{{\bf r}_j}\cos\theta_{ijk} &=&
 \nabla_{{\bf r}_j}\Big(\frac{{\bf r}_{ij}{\bf }r_{ik}}{r_{ij}r_{ik}}\Big)
 \nonumber \\
 &=& \frac{1}{r_{ij}r_{ik}}{\bf r}_{ik} -
     \frac{\cos\theta_{ijk}}{r_{ij}^2}{\bf r}_{ij}
\end{eqnarray}

  \subsection{Derivative of $V_{ij}$ with respect to ${\bf r}_k$}

The derivative of $V_{ij}$ with respect to ${\bf r}_k$ just consists of the
single term
\begin{eqnarray}
 \nabla_{{\bf r}_k} V_{ij} & = &
 f_C(r_{ij})f_A(r_{ij})\nabla_{{\bf r}_{k}}b_{ij}
\end{eqnarray}
since the derivatives of the functions only depending on atom $i$ and $j$
vanish.
\begin{eqnarray}
 \nabla_{{\bf r}_k} f_R(r_{ij}) &=& 0 \\
 \nabla_{{\bf r}_k} f_A(r_{ij}) &=& 0 \\
 \nabla_{{\bf r}_k} f_C(r_{ij}) &=& 0
\end{eqnarray}
Now look at $b_{ij}$, not only angle important here!
\begin{eqnarray}
 \nabla_{{\bf r}_k}\zeta_{ij} &=&
 g(\theta_{ijk})\nabla_{{\bf r}_k}f_C(r_{ik}) +
 f_C(r_{ik})\nabla_{{\bf r}_k}g(\theta_{ijk}) \\
 \nabla_{{\bf r}_k}f_C(r_{ik}) &=& - \nabla_{{\bf r}_i}f_C(r_{ik}) \\
 \nabla_{{\bf r}_k}\cos\theta_{ijk} &=&
 \nabla_{{\bf r}_k}\Big(\frac{{\bf r}_{ij}{\bf }r_{ik}}{r_{ij}r_{ik}}\Big)
 \nonumber \\
 &=&\frac{1}{r_{ij}r_{ik}}{\bf r}_{ij} -
 \frac{\cos\theta_{ijk}}{r_{ik}^2}{\bf r}_{ik}
\end{eqnarray}

  \subsection{Code realization}

The implementation of the force evaluation shown in the following
is applied to the potential designed by Erhard and Albe.
There are slight differences comparted to the original potential by Tersoff:
\begin{itemize}
 \item Difference in sign of the attractive part.
 \item $c$, $d$ and $h$ values depend on atom $k$ in addition to atom $i$.
 \item Difference in sign of the $\cos\theta_{ijk}$ term.
 \item There are no parameters $\beta$ and $\chi$.
 \item The exponent of the $b$ term is constantly $-\frac{1}{2}$.
\end{itemize}
These differences actually slightly ease code realization.

\begin{figure}
\renewcommand\labelitemi{}
\renewcommand\labelitemii{}
\renewcommand\labelitemiii{}
LOOP i \{
\begin{itemize}
 \item // nop (only used in orig. Tersoff)
 \item LOOP j \{
       \begin{itemize}
        \item // nop (only used in orig. Tersoff)
       \end{itemize}
 \item \}
 \item LOOP j \{
       \begin{itemize}
        \item $\zeta_{ij}=0$
        \item set $S_{ij}$ (cutoff)
        \item calculate: $r_{ij}$, $r_{ij}^2$
        \item IF $r_{ij} > S_{ij}$ THEN CONTINUE
	\item
	\item LOOP k \{
	      \begin{itemize}
               \item set $ik$-depending values
               \item calculate: $r_{ik}$, $r_{ik}^2$
	       \item IF $r_{ik} > S_{ik}$ THEN CONTINUE
	       \item calculate: $\theta_{ijk}$, $\cos(\theta_{ijk})$,
	                        $dg(\theta_{ijk})$, $g(\theta)$,
	                        $f_C(r_{ik})$, $df_C(r_{ik})$
	       \item $\zeta_{ij}=\zeta_{ij}+f_C(r_{ik})g(\theta)$
	      \end{itemize}
        \item \}
	\item
	\item calculate: $f_C(r_{ij})$, $df_C(r_{ij})$, $f_A(r_{ij})$,
	                 $df_A(r_{ij})$, $f_R(r_{ij})$, $df_R(r_{ij})$,
			 $b_{ij}$, $db_{ij}$
	\item calculate:
$
F=-\frac{1}{2}\big(
\nabla_{{\bf r}_i} f_C(r_{ij}) \big[ f_R(r_{ij}) - b_{ij} f_A(r_{ij}) \big] +
f_C(r_{ij}) \big[ \nabla_{{\bf r}_i} f_R(r_{ij}) -
                  b_{ij} \nabla_{{\bf r}_i} f_A(r_{ij}) \big]\big)
$
	\item $F_{Atom\, i} = F_{Atom\, i} + F$
	\item $F_{Atom\, j} = F_{Atom\, j} - F$
	\item $E=E+\frac{1}{2}f_C(r_{ij})[f_R(r_{ij})-b_{ij}f_A(r_{ij})]$
	\item $d\zeta_{ij}=\frac{1}{2}f_A(r_{ij})f_C(r_{ij})db_{ij}$
	\item
	\item LOOP k \{
	      \begin{itemize}
               \item calculate: $\nabla_{{\bf r}_i}\cos\theta_{ijk}$,
	                        $\nabla_{{\bf r}_j}\cos\theta_{ijk}$,
				$\nabla_{{\bf r}_k}\cos\theta_{ijk}$
               \item $
F_{Atom\, i}+= d\zeta_{ij}\big(
g(\theta_{ijk})\nabla_{{\bf r}_i}f_C(r_{ik}) +
f_C(r_{ik})dg(\theta_{ijk})\nabla_{{\bf r}_i}\cos\theta_{ijk}\big)$
               \item $
F_{Atom\, j}+= d\zeta_{ij}
f_C(r_{ik})dg(\theta_{ijk})\nabla_{{\bf r}_j}\cos\theta_{ijk}$
               \item $
F_{Atom\, k}+= d\zeta_{ij}\big(
g(\theta_{ijk})\nabla_{{\bf r}_k}f_C(r_{ik}) +
f_C(r_{ik})dg(\theta_{ijk})\nabla_{{\bf r}_k}\cos\theta_{ijk}\big)$
	      \end{itemize}
        \item \}
       \end{itemize}
 \item \}
\end{itemize}
\}
\caption{Implementation of the force evaluation for Tersoff like bond-order
         potentials using pseudocode.}
\end{figure}

