Index: src/Potentials/Makefile.am
===================================================================
--- src/Potentials/Makefile.am	(revision 6bb72a1d333eb13b241f2296056b61a3f431aa72)
+++ src/Potentials/Makefile.am	(revision 610c11a5323c7e86a4ffb55afc832f78510a0136)
@@ -3,4 +3,5 @@
 
 POTENTIALSSOURCE = \
+  Potentials/Specifics/ManyBodyPotential_Tersoff.cpp \
   Potentials/Specifics/PairPotential_Harmonic.cpp
 
@@ -8,6 +9,6 @@
   Potentials/EmpiricalPotential.hpp \
   Potentials/helpers.hpp \
+  Potentials/Specifics/ManyBodyPotential_Tersoff.hpp \
   Potentials/Specifics/PairPotential_Harmonic.hpp
-
 
 noinst_LTLIBRARIES += libMolecuilderPotentials.la
Index: src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp
===================================================================
--- src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp	(revision 610c11a5323c7e86a4ffb55afc832f78510a0136)
+++ src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp	(revision 610c11a5323c7e86a4ffb55afc832f78510a0136)
@@ -0,0 +1,165 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2012 University of Bonn. All rights reserved.
+ * Please see the COPYING file or "Copyright notice" in builder.cpp for details.
+ * 
+ *
+ *   This file is part of MoleCuilder.
+ *
+ *    MoleCuilder is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 2 of the License, or
+ *    (at your option) any later version.
+ *
+ *    MoleCuilder is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>. 
+ */
+
+/*
+ * ManyBodyPotential_Tersoff.cpp
+ *
+ *  Created on: Sep 26, 2012
+ *      Author: heber
+ */
+
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include "ManyBodyPotential_Tersoff.hpp"
+
+#include <boost/bind.hpp>
+#include <cmath>
+
+#include "CodePatterns/Assert.hpp"
+
+#include "Potentials/helpers.hpp"
+
+EmpiricalPotential::result_t
+ManyBodyPotential_Tersoff::operator()(
+    const arguments_t &arguments
+    ) const
+{
+  const double &distance = arguments[0].distance;
+  const double cutoff = function_cutoff(distance);
+  const double result = (cutoff == 0.) ? 0. : cutoff * (
+      function_prefactor(
+          manybodyparameter_alpha,
+          boost::bind(&ManyBodyPotential_Tersoff::function_eta,
+              boost::cref(*this),
+              boost::cref(arguments[0])))
+      * function_smoother(
+          distance,
+          manybodyparameter_A,
+          manybodyparameter_lambda1)
+      +
+      function_prefactor(
+          manybodyparameter_beta,
+          boost::bind(&ManyBodyPotential_Tersoff::function_zeta,
+              boost::cref(*this),
+              boost::cref(arguments[0])))
+      * function_smoother(
+          distance,
+          -manybodyparameter_B,
+          manybodyparameter_lambda2)
+  );
+  return result;
+}
+
+EmpiricalPotential::result_t
+ManyBodyPotential_Tersoff::function_cutoff(
+    const double &distance
+  ) const
+{
+  const double offset = (distance - cutoff_offset);
+  if (offset < - cutoff_halfwidth)
+    return 1.;
+  else if (offset > cutoff_halfwidth)
+    return 0.;
+  else {
+    return (0.5 - 0.5 * sin( .5 * M_PI * offset/cutoff_halfwidth));
+  }
+}
+
+EmpiricalPotential::result_t
+ManyBodyPotential_Tersoff::function_prefactor(
+    const double &alpha,
+    boost::function<result_t()> etafunction
+  ) const
+{
+  return pow(
+      (1. + Helpers::pow(alpha * etafunction(), manybodyparameter_n)),
+      -1./(2.*manybodyparameter_n));
+}
+
+EmpiricalPotential::result_t
+ManyBodyPotential_Tersoff::function_eta(
+    const argument_t &r_ij
+  ) const
+{
+  result_t result = 0.;
+
+  // get all triples within the cutoff
+  std::vector<arguments_t> triples = triplefunction(r_ij, cutoff_offset+cutoff_halfwidth);
+  for (std::vector<arguments_t>::const_iterator iter = triples.begin();
+      iter != triples.end(); ++iter) {
+    ASSERT( iter->size() == 2,
+        "ManyBodyPotential_Tersoff::function_zeta() - the triples result must contain of exactly two distances.");
+    const argument_t &r_ik = (*iter)[0];
+    result += function_cutoff(r_ik.distance)
+        * exp( Helpers::pow(manybodyparameter_lambda3 * (r_ij.distance - r_ik.distance) ,3));
+  }
+
+  return result;
+}
+
+EmpiricalPotential::result_t
+ManyBodyPotential_Tersoff::function_zeta(
+    const argument_t &r_ij
+  ) const
+{
+  result_t result = 0.;
+
+  // get all triples within the cutoff
+  std::vector<arguments_t> triples = triplefunction(r_ij, cutoff_offset+cutoff_halfwidth);
+  for (std::vector<arguments_t>::const_iterator iter = triples.begin();
+      iter != triples.end(); ++iter) {
+    ASSERT( iter->size() == 2,
+        "ManyBodyPotential_Tersoff::function_zeta() - the triples result must contain exactly two distances.");
+    const argument_t &r_ik = (*iter)[0];
+    const argument_t &r_jk = (*iter)[1];
+    result +=
+        function_cutoff(r_ik.distance)
+        * function_angle(r_ij.distance, r_ik.distance, r_jk.distance)
+        * exp( Helpers::pow(manybodyparameter_lambda3 * (r_ij.distance - r_ik.distance) ,3));
+  }
+
+  return result;
+}
+
+EmpiricalPotential::result_t
+ManyBodyPotential_Tersoff::function_angle(
+    const double &r_ij,
+    const double &r_ik,
+    const double &r_jk
+  ) const
+{
+  const double angle = Helpers::pow(r_ij,2) + Helpers::pow(r_ik,2) - Helpers::pow(r_jk,2);
+  const double divisor = r_ij * r_ik;
+  const double result =
+      1.
+      + (Helpers::pow(manybodyparameter_c, 2)/Helpers::pow(manybodyparameter_d, 2))
+      - Helpers::pow(manybodyparameter_c, 2)/(Helpers::pow(manybodyparameter_d, 2) +
+          Helpers::pow(manybodyparameter_h - cos(angle/divisor),2));
+  return result;
+}
Index: src/Potentials/Specifics/ManyBodyPotential_Tersoff.hpp
===================================================================
--- src/Potentials/Specifics/ManyBodyPotential_Tersoff.hpp	(revision 610c11a5323c7e86a4ffb55afc832f78510a0136)
+++ src/Potentials/Specifics/ManyBodyPotential_Tersoff.hpp	(revision 610c11a5323c7e86a4ffb55afc832f78510a0136)
@@ -0,0 +1,174 @@
+/*
+ * ManyBodyPotential_Tersoff.hpp
+ *
+ *  Created on: Sep 26, 2012
+ *      Author: heber
+ */
+
+#ifndef MANYBODYPOTENTIAL_TERSOFF_HPP_
+#define MANYBODYPOTENTIAL_TERSOFF_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "Potentials/EmpiricalPotential.hpp"
+
+#include <boost/function.hpp>
+#include <cmath>
+
+/** This class is the implementation of the Tersoff potential function.
+ *
+ * \note The arguments_t argument list is here in the following order:
+ * -# first \f$ r_{ij}$ \f$,
+ * -# then all \f$ r_{ik}$ \f$ that are within the cutoff, i.e. \f$ r_{ik}$ < R + D\f$
+ *
+ */
+class ManyBodyPotential_Tersoff : public EmpiricalPotential
+{
+private:
+  //!> cutoff_offset offset for cutoff parameter
+  const double cutoff_offset;
+  //!> cutoff_halfwidth half-width for cutoff parameter
+  const double cutoff_halfwidth;
+  //!> many body prefactor parameter for attractive part
+  const double manybodyparameter_A;
+  //!> many body prefactor parameter for attractive part
+  const double manybodyparameter_B;
+  //!> many body scale parameter for attractive part
+  const double manybodyparameter_lambda1;
+  //!> many body scale parameter for attractive part
+  const double manybodyparameter_lambda2;
+  //!> length scale for coordination influence
+  const double manybodyparameter_lambda3;
+  //!> many body parameter for attractive part
+  const double manybodyparameter_alpha;
+  //!> many body parameter for repulsive part
+  const double manybodyparameter_beta;
+  //!> many body type-dependent parameter giving the power
+  const double manybodyparameter_n;
+  //!> many-body type-dependent parameter in angular dependence
+  const double manybodyparameter_c;
+  //!> many-body type-dependent parameter in angular dependence
+  const double manybodyparameter_d;
+  //!> many-body type-dependent parameter in angular dependence
+  const double manybodyparameter_h;
+public:
+  /** Constructor for class ManyBodyPotential_Tersoff.
+   *
+   * @param _triplefunction function that returns a list of triples (i.e. the
+   *        two remaining distances) to a given pair of points (contained as
+   *        indices within the argument)
+   */
+  ManyBodyPotential_Tersoff(
+      const double _cutoff_offset,
+      const double _cutoff_halfwidth,
+      const double A,
+      const double B,
+      const double lambda1,
+      const double lambda2,
+      const double lambda3,
+      const double alpha,
+      const double beta,
+      const double n,
+      const double c,
+      const double d,
+      const double h,
+      boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction) :
+    cutoff_offset(_cutoff_offset),
+    cutoff_halfwidth(_cutoff_offset),
+    manybodyparameter_A(A),
+    manybodyparameter_B(B),
+    manybodyparameter_lambda1(lambda1),
+    manybodyparameter_lambda2(lambda2),
+    manybodyparameter_lambda3(lambda3),
+    manybodyparameter_alpha(alpha),
+    manybodyparameter_beta(beta),
+    manybodyparameter_n(n),
+    manybodyparameter_c(d),
+    manybodyparameter_d(d),
+    manybodyparameter_h(h),
+    triplefunction(_triplefunction)
+  {}
+  /** Destructor of class ManyBodyPotential_Tersoff.
+   *
+   */
+  virtual ~ManyBodyPotential_Tersoff() {}
+
+  /** Evaluates the Tersoff potential for the given arguments.
+   *
+   * @param arguments single distance
+   * @return value of the potential function
+   */
+  result_t operator()(const arguments_t &arguments) const;
+
+  /** Evaluates the derivative of the Tersoff potential with respect to the
+   * input variables.
+   *
+   * @param arguments single distance
+   * @return vector with components of the derivative
+   */
+  derivative_components_t derivative(const arguments_t &arguments) const;
+
+private:
+  /** This function represents the cutoff \f$ f_C \f$.
+   *
+   * @param distance variable of the function
+   * @return a value in [0,1].
+   */
+  result_t function_cutoff(
+      const double &distance
+      ) const;
+  /** This  function has the exponential feature from the Morse potential.
+   *
+   * @param distance variable of the function
+   * @param prefactor prefactor parameter to exp function
+   * @param lambda scale parameter of exp function's argument
+   * @return
+   */
+  result_t function_smoother(
+      const double &distance,
+      const double &prefactor,
+      const double &lambda
+      ) const
+  {
+    return prefactor * exp(-lambda * distance);
+  }
+
+  /** This function represents \f$ (1 + \alpha^n \eta^n)^{-1/2n} \f$.
+   *
+   * @param alpha prefactor to eta function
+   * @param r_ij distance argument
+   * @param etafunction eta or zeta
+   * @return \f$ (1 + \alpha^n \eta^n)^{-1/2n} \f$
+   */
+  result_t function_prefactor(
+      const double &alpha,
+      boost::function<result_t()> etafunction
+        ) const;
+
+  result_t
+  function_eta(
+      const argument_t &r_ij
+    ) const;
+
+  result_t
+  function_zeta(
+      const argument_t &r_ij
+    ) const;
+
+  result_t
+  function_angle(
+      const double &r_ij,
+      const double &r_ik,
+      const double &r_jk
+    ) const;
+
+private:
+  //!> bound function that obtains the triples for the internal coordinationb summation.
+  const boost::function< std::vector< arguments_t >(const argument_t &, const double)> &triplefunction;
+};
+
+
+#endif /* MANYBODYPOTENTIAL_TERSOFF_HPP_ */
