+If not only pair potentials are considered, three body potentials $U_3$ or multi 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.
+A three body potential has to be included for these types of elements.
+
+\subsubsection{The Tersoff bond order potential}
+
+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 increases 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 depends upon the local environment.
+An atom with many neighbors forms weaker bonds than an atom with few neighbors since .
+
+Here comes an explanation, energy per bond monotonically decreasing with the amount of bonds and so on and so on \ldots
+
+The interatomic potential is taken to have 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}
+$E$ is the total energy of the system, constituted either by the sum over the site energies $E_i$ or by the bond energies $V_{ij}$.
+The indices $i$ and $j$ correspond to the atoms of the system with $r_{ij}$ being the distance from atom $i$ to atom $j$.
+The functions $f_R$ and $f_A$ represent a repulsive and an attractive pair potential.
+The repulsive part is due to the orthogonalization energy of overlapped atomic wave functions.
+The attractive part is associated with the bonding.
+\begin{eqnarray}
+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 the a cutoff 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\{
+ \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.
+\label{eq:basics:fc}
+\end{equation}
+The function $b_{ij}$ represents a measure of the bond order, monotonically decreasing with the coordination of atoms $i$ and $j$.
+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}) \\
+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 Figure \ref{img:tersoff_angle}.
+\begin{figure}[t]
+\begin{center}
+\includegraphics[width=8cm]{tersoff_angle.eps}
+\end{center}
+\caption{Angle between bonds of atoms $i,j$ and $i,k$.}
+\label{img:tersoff_angle}
+\end{figure}
+
+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
+\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}
+The force is then given by
+\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}.