Index: src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp
===================================================================
--- src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp	(revision e7579ed03a2594b79fb214806a5efc0a927b5480)
+++ src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp	(revision ffc3688dd044a0c874c10a37de01139caff883c1)
@@ -43,4 +43,5 @@
 
 #include "CodePatterns/Assert.hpp"
+//#include "CodePatterns/Info.hpp"
 
 #include "Potentials/helpers.hpp"
@@ -54,13 +55,15 @@
 
 ManyBodyPotential_Tersoff::ManyBodyPotential_Tersoff(
-    const double &_cutoff_offset,
-    const double &_cutoff_halfwidth,
+    const double &_R,
+    const double &_S,
     const double &_A,
     const double &_B,
-    const double &_lambda1,
-    const double &_lambda2,
+    const double &_lambda,
+    const double &_mu,
     const double &_lambda3,
     const double &_alpha,
     const double &_beta,
+    const double &_chi,
+    const double &_omega,
     const double &_n,
     const double &_c,
@@ -71,13 +74,16 @@
   triplefunction(_triplefunction)
 {
-  params[cutoff_offset] = _cutoff_offset;
-  params[cutoff_halfwidth] = _cutoff_halfwidth;
+//  Info info(__func__);
+  params[R] = _R;
+  params[S] = _S;
   params[A] = _A;
   params[B] = _B;
-  params[lambda1] = _lambda1;
-  params[lambda2] = _lambda2;
+  params[lambda] = _lambda;
+  params[mu] = _mu;
   params[lambda3] = _lambda3;
   params[alpha] = _alpha;
   params[beta] = _beta;
+  params[chi] = _chi;
+  params[omega] = _omega;
   params[n] = _n;
   params[c] = _c;
@@ -91,28 +97,189 @@
     ) const
 {
-  const double &distance = arguments[0].distance;
-  const double cutoff = function_cutoff(distance);
-  const double result = (cutoff == 0.) ? 0. : cutoff * (
-      function_prefactor(
-          alpha,
-          boost::bind(&ManyBodyPotential_Tersoff::function_eta,
-              boost::cref(*this),
-              boost::cref(arguments[0])))
-      * function_smoother(
-          distance,
-          A,
-          lambda1)
-      +
-      function_prefactor(
-          beta,
-          boost::bind(&ManyBodyPotential_Tersoff::function_zeta,
-              boost::cref(*this),
-              boost::cref(arguments[0])))
-      * function_smoother(
-          distance,
-          -B,
-          lambda2)
-  );
+//  Info info(__func__);
+  const argument_t &r_ij = arguments[0];
+  const double cutoff = function_cutoff(r_ij.distance);
+  const double result = (cutoff == 0.) ?
+      0. :
+      cutoff * (
+          function_prefactor(
+              params[alpha],
+              function_eta(r_ij))
+          * function_smoother(
+              params[A],
+              params[lambda],
+              r_ij.distance)
+          +
+          function_prefactor(
+              params[beta],
+              function_zeta(r_ij))
+          * function_smoother(
+              -params[B],
+              params[mu],
+              r_ij.distance)
+      );
+//  LOG(2, "DEBUG: operator()(" << r_ij.distance << ") = " << result);
   return std::vector<result_t>(1, result);
+}
+
+ManyBodyPotential_Tersoff::derivative_components_t
+ManyBodyPotential_Tersoff::derivative(
+    const arguments_t &arguments
+    ) const
+{
+//  Info info(__func__);
+  return ManyBodyPotential_Tersoff::derivative_components_t();
+}
+
+ManyBodyPotential_Tersoff::results_t
+ManyBodyPotential_Tersoff::parameter_derivative(
+    const arguments_t &arguments,
+    const size_t index
+    ) const
+{
+//  Info info(__func__);
+//  ASSERT( arguments.size() == 1,
+//      "PairPotential_Harmonic::parameter_derivative() - requires exactly one argument.");
+  const argument_t &r_ij = arguments[0];
+  switch (index) {
+    case R:
+    {
+      const double result = 0.;
+      return results_t(1, result);
+      break;
+    }
+    case S:
+    {
+      const double result = 0.;
+      return results_t(1, result);
+      break;
+    }
+    case A:
+    {
+      const double result = 0.;
+      return results_t(1, result);
+      break;
+    }
+    case B:
+    {
+      const double result = 0.;
+      return results_t(1, result);
+      break;
+    }
+    case lambda:
+    {
+      const double cutoff = function_cutoff(r_ij.distance);
+      const double result = (cutoff == 0.) ?
+          0. :
+          -r_ij.distance * cutoff * params[lambda] * (
+          function_prefactor(
+              params[alpha],
+              function_eta(r_ij))
+          * function_smoother(
+              params[A],
+              params[lambda],
+              r_ij.distance)
+          );
+      return results_t(1, result);
+      break;
+    }
+    case mu:
+    {
+      const double cutoff = function_cutoff(r_ij.distance);
+      const double result = (cutoff == 0.) ?
+          0. :
+          -r_ij.distance * cutoff * params[mu] *(
+          function_prefactor(
+              params[beta],
+              function_zeta(r_ij))
+          * function_smoother(
+              -params[B],
+              params[mu],
+              r_ij.distance)
+      );
+      return results_t(1, result);
+      break;
+    }
+    case lambda3:
+    {
+      const double result = 0.;
+      return results_t(1, result);
+      break;
+    }
+    case alpha:
+    {
+      const double temp =
+          pow(params[alpha]*function_eta(r_ij), params[n]);
+      const double cutoff = function_cutoff(r_ij.distance);
+      const double result = (cutoff == 0.) || (params[alpha] == 0. )?
+          0. :
+          function_smoother(
+              -params[A],
+              params[lambda],
+              r_ij.distance)
+          * (-.5) * params[alpha] * (temp/params[alpha])
+          / (1. + temp)
+          ;
+      return results_t(1, result);
+      break;
+    }
+    case beta:
+    {
+      const double temp =
+          pow(params[beta]*function_zeta(r_ij), params[n]);
+      const double cutoff = function_cutoff(r_ij.distance);
+      const double result = (cutoff == 0.) || (params[beta] == 0. )?
+          0. : cutoff *
+          function_smoother(
+              -params[B],
+              params[mu],
+              r_ij.distance)
+          * (-.5) * params[beta] * (temp/params[beta])
+          / (1. + temp)
+          ;
+      return results_t(1, result);
+      break;
+    }
+    case n:
+    {
+      const double temp =
+          pow(params[beta]*function_zeta(r_ij), params[n]);
+      const double cutoff = function_cutoff(r_ij.distance);
+      const double result = (cutoff == 0.) ?
+          0. : cutoff *
+          function_smoother(
+              -params[B],
+              params[mu],
+              r_ij.distance)
+          * params[B]
+          * ( log(1.+temp)/(params[n]*params[n]) - temp
+              * (log(function_zeta(r_ij)) + log(params[beta]))
+              /(params[n]*(1.+temp)))
+          ;
+      return results_t(1, result);
+      break;
+    }
+    case c:
+    {
+      const double result = 0.;
+      return results_t(1, result);
+      break;
+    }
+    case d:
+    {
+      const double result = 0.;
+      return results_t(1, result);
+      break;
+    }
+    case h:
+    {
+      const double result = 0.;
+      return results_t(1, result);
+      break;
+    }
+    default:
+      break;
+  }
+  return results_t(1, 0.);
 }
 
@@ -122,12 +289,15 @@
   ) const
 {
-  const double offset = (distance - cutoff_offset);
-  if (offset < - cutoff_halfwidth)
-    return 1.;
-  else if (offset > cutoff_halfwidth)
-    return 0.;
+//  Info info(__func__);
+  double result = 0.;
+  if (distance < params[R])
+    result = 1.;
+  else if (distance > params[S])
+    result = 0.;
   else {
-    return (0.5 - 0.5 * sin( .5 * M_PI * offset/cutoff_halfwidth));
+    result = (0.5 + 0.5 * cos( M_PI * (distance - params[R])/(params[S]-params[R])));
   }
+//  LOG(2, "DEBUG: function_cutoff(" << distance << ") = " << result);
+  return result;
 }
 
@@ -135,10 +305,26 @@
 ManyBodyPotential_Tersoff::function_prefactor(
     const double &alpha,
-    boost::function<result_t()> etafunction
+    const double &eta
   ) const
 {
-  return pow(
-      (1. + Helpers::pow(alpha * etafunction(), n)),
-      -1./(2.*n));
+//  Info info(__func__);
+  const double result = params[chi] * pow(
+      (1. + pow(alpha * eta, params[n])),
+      -1./(2.*params[n]));
+//  LOG(2, "DEBUG: function_prefactor(" << alpha << "," << eta << ") = " << result);
+  return result;
+}
+
+ManyBodyPotential_Tersoff::result_t
+ManyBodyPotential_Tersoff::function_smoother(
+    const double &prefactor,
+    const double &lambda,
+    const double &distance
+    ) const
+{
+//  Info info(__func__);
+  const double result = prefactor * exp(-lambda * distance);
+//  LOG(2, "DEBUG: function_smoother(" << prefactor << "," << lambda << "," << distance << ") = " << result);
+  return result;
 }
 
@@ -148,8 +334,9 @@
   ) const
 {
+//  Info info(__func__);
   result_t result = 0.;
 
   // get all triples within the cutoff
-  std::vector<arguments_t> triples = triplefunction(r_ij, cutoff_offset+cutoff_halfwidth);
+  std::vector<arguments_t> triples = triplefunction(r_ij, params[S]);
   for (std::vector<arguments_t>::const_iterator iter = triples.begin();
       iter != triples.end(); ++iter) {
@@ -158,7 +345,8 @@
     const argument_t &r_ik = (*iter)[0];
     result += function_cutoff(r_ik.distance)
-        * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3));
+        * exp( Helpers::pow(params[lambda3] * (r_ij.distance - r_ik.distance) ,3));
   }
 
+//  LOG(2, "DEBUG: function_eta(" << r_ij.distance << ") = " << result);
   return result;
 }
@@ -169,8 +357,9 @@
   ) const
 {
+//  Info info(__func__);
   result_t result = 0.;
 
   // get all triples within the cutoff
-  std::vector<arguments_t> triples = triplefunction(r_ij, cutoff_offset+cutoff_halfwidth);
+  std::vector<arguments_t> triples = triplefunction(r_ij, params[S]);
   for (std::vector<arguments_t>::const_iterator iter = triples.begin();
       iter != triples.end(); ++iter) {
@@ -181,8 +370,10 @@
     result +=
         function_cutoff(r_ik.distance)
+        * params[omega]
         * function_angle(r_ij.distance, r_ik.distance, r_jk.distance)
-        * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3));
+        * exp( Helpers::pow(params[lambda3] * (r_ij.distance - r_ik.distance) ,3));
   }
 
+//  LOG(2, "DEBUG: function_zeta(" << r_ij.distance << ") = " << result);
   return result;
 }
@@ -195,11 +386,19 @@
   ) const
 {
+//  Info info(__func__);
   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(c, 2)/Helpers::pow(d, 2))
-      - Helpers::pow(c, 2)/(Helpers::pow(d, 2) +
-          Helpers::pow(h - cos(angle/divisor),2));
-  return result;
-}
+  const double divisor = 2.* r_ij * r_ik;
+//  LOG(2, "DEBUG: cos(theta)= " << angle/divisor);
+  if (divisor != 0.) {
+    const double result =
+        1.
+        + (Helpers::pow(params[c]/params[d], 2))
+        - Helpers::pow(params[c], 2)/(Helpers::pow(params[d], 2) +
+            Helpers::pow(params[h] - angle/divisor,2));
+
+//    LOG(2, "DEBUG: function_angle(" << r_ij << "," << r_ik << "," << r_jk << ") = " << result);
+    return result;
+  } else
+    return 0.;
+}
+
Index: src/Potentials/Specifics/ManyBodyPotential_Tersoff.hpp
===================================================================
--- src/Potentials/Specifics/ManyBodyPotential_Tersoff.hpp	(revision e7579ed03a2594b79fb214806a5efc0a927b5480)
+++ src/Potentials/Specifics/ManyBodyPotential_Tersoff.hpp	(revision ffc3688dd044a0c874c10a37de01139caff883c1)
@@ -17,8 +17,8 @@
 #include <cmath>
 
+#include "CodePatterns/Assert.hpp"
+
 #include "Potentials/EmpiricalPotential.hpp"
 #include "FunctionApproximation/FunctionModel.hpp"
-
-#include "CodePatterns/Assert.hpp"
 
 /** This class is the implementation of the Tersoff potential function.
@@ -31,4 +31,6 @@
 class ManyBodyPotential_Tersoff : virtual public EmpiricalPotential, virtual public FunctionModel
 {
+  //!> grant unit test access to internal parts
+  friend class ManyBodyPotential_TersoffTest;
   // some repeated typedefs to avoid ambiguities
   typedef FunctionModel::arguments_t arguments_t;
@@ -50,13 +52,15 @@
   /** Constructor for class ManyBodyPotential_Tersoff.
    *
-   * @param _cutoff_offset offset for cutoff
-   * @param _cutoff_halfwidth halfwidth for cutoff relative to \a _cutoff_offset
+   * @param _R offset for cutoff
+   * @param _S halfwidth for cutoff relative to \a _R
    * @param A
    * @param B
-   * @param lambda1
-   * @param lambda2
+   * @param lambda
+   * @param mu
    * @param lambda3
    * @param alpha
    * @param beta
+   * @param chi
+   * @param omega
    * @param n
    * @param c
@@ -68,17 +72,19 @@
    */
   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,
+      const double &_R,
+      const double &_S,
+      const double &_A,
+      const double &_B,
+      const double &_lambda,
+      const double &_mu,
+      const double &_lambda3,
+      const double &_alpha,
+      const double &_beta,
+      const double &_chi,
+      const double &_omega,
+      const double &_n,
+      const double &_c,
+      const double &_d,
+      const double &_h,
       boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction);
 
@@ -132,17 +138,14 @@
   /** 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
+   * @param distance variable of the function
    * @return
    */
   result_t function_smoother(
-      const double &distance,
       const double &prefactor,
-      const double &lambda
-      ) const
-  {
-    return prefactor * exp(-lambda * distance);
-  }
+      const double &lambda,
+      const double &distance
+      ) const;
 
   /** This function represents \f$ (1 + \alpha^n \eta^n)^{-1/2n} \f$.
@@ -150,10 +153,10 @@
    * @param alpha prefactor to eta function
    * @param r_ij distance argument
-   * @param etafunction eta or zeta
+   * @param eta result value of 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 double &eta
         ) const;
 
@@ -177,18 +180,19 @@
 private:
   enum parameter_enum_t {
-    cutoff_offset=0,
-    equilibrium_distance=1,
-    cutoff_halfwidth=2,
-    A=3,
-    B=4,
-    lambda1=5,
-    lambda2=6,
-    lambda3=7,
-    alpha=8,
-    beta=9,
-    n=10,
-    c=11,
-    d=12,
-    h=13,
+    R=0,
+    S=1,
+    A=2,
+    B=3,
+    lambda=4,
+    mu=5,
+    lambda3=6,
+    alpha=7,
+    beta=8,
+    chi=9,
+    omega=10,
+    n=11,
+    c=12,
+    d=13,
+    h=14,
     MAXPARAMS
   };
Index: src/Potentials/Specifics/unittests/Makefile.am
===================================================================
--- src/Potentials/Specifics/unittests/Makefile.am	(revision e7579ed03a2594b79fb214806a5efc0a927b5480)
+++ src/Potentials/Specifics/unittests/Makefile.am	(revision ffc3688dd044a0c874c10a37de01139caff883c1)
@@ -3,12 +3,15 @@
 
 POTENTIALSSPECIFICSTESTSSOURCES = \
+	../Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.cpp \
 	../Potentials/Specifics/unittests/PairPotential_HarmonicUnitTest.cpp \
 	../Potentials/Specifics/unittests/PairPotential_MorseUnitTest.cpp
 
 POTENTIALSSPECIFICSTESTSHEADERS = \
+	../Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.hpp \
 	../Potentials/Specifics/unittests/PairPotential_HarmonicUnitTest.hpp \
 	../Potentials/Specifics/unittests/PairPotential_MorseUnitTest.hpp
 
 POTENTIALSSPECIFICSTESTS = \
+	ManyBodyPotential_TersoffUnitTest \
 	PairPotential_HarmonicUnitTest \
 	PairPotential_MorseUnitTest
@@ -20,7 +23,13 @@
 POTENTIALSSPECIFICSLIBS = \
 	../libMolecuilderPotentials.la \
+	$(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \
 	${CodePatterns_LIBS} \
   $(BOOST_SERIALIZATION_LDFLAGS) $(BOOST_SERIALIZATION_LIBS) \
 	$(BOOST_LIB)
+
+ManyBodyPotential_TersoffUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
+	../Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.cpp \
+	../Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.hpp
+ManyBodyPotential_TersoffUnitTest_LDADD = ${POTENTIALSSPECIFICSLIBS}
 
 PairPotential_HarmonicUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
Index: src/Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.cpp
===================================================================
--- src/Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.cpp	(revision ffc3688dd044a0c874c10a37de01139caff883c1)
+++ src/Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.cpp	(revision ffc3688dd044a0c874c10a37de01139caff883c1)
@@ -0,0 +1,378 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2012 University of Bonn. All rights reserved.
+ * 
+ *
+ *   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_TersoffUnitTest.cpp
+ *
+ *  Created on: Oct 04, 2012
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include "ManyBodyPotential_TersoffUnitTest.hpp"
+
+#include <boost/assign.hpp>
+#include <boost/function.hpp>
+
+#include <limits>
+
+#include "CodePatterns/Assert.hpp"
+#include "CodePatterns/Log.hpp"
+
+#include "FunctionApproximation/FunctionArgument.hpp"
+#include "Potentials/Specifics/ManyBodyPotential_Tersoff.hpp"
+#include "Potentials/helpers.hpp"
+
+using namespace boost::assign;
+
+#ifdef HAVE_TESTRUNNER
+#include "UnitTestMain.hpp"
+#endif /*HAVE_TESTRUNNER*/
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( ManyBodyPotential_TersoffTest );
+
+ManyBodyPotential_TersoffTest::configurations_t ManyBodyPotential_TersoffTest::configurations;
+ManyBodyPotential_TersoffTest::configuration_t *ManyBodyPotential_TersoffTest::CurrentConfiguration = NULL;
+
+/** This function looks up all distances ik and jk to a given ij and
+ * returns a vector of pairs.
+ */
+std::vector<FunctionModel::arguments_t>
+triplefunction(const argument_t &arguments, const double cutoff)
+{
+  std::vector<FunctionModel::arguments_t> result;
+  // go through current configuration and gather all other distances
+  ManyBodyPotential_TersoffTest::configuration_t::const_iterator firstiter =
+      ManyBodyPotential_TersoffTest::CurrentConfiguration->begin();
+  std::advance(firstiter, arguments.indices.first);
+  ManyBodyPotential_TersoffTest::configuration_t::const_iterator seconditer =
+      ManyBodyPotential_TersoffTest::CurrentConfiguration->begin();
+  std::advance(seconditer, arguments.indices.second);
+  for (ManyBodyPotential_TersoffTest::configuration_t::const_iterator iter =
+      ManyBodyPotential_TersoffTest::CurrentConfiguration->begin();
+      iter != ManyBodyPotential_TersoffTest::CurrentConfiguration->end();
+      ++iter) {
+    // skip k==i and k==j
+    if ((iter == firstiter) || (iter == seconditer))
+      continue;
+    FunctionModel::arguments_t args(2);
+    // ik
+    args[0].distance = firstiter->distance(*iter);
+    args[0].indices = std::make_pair(
+        std::distance( // enforce const_iterator return from begin()
+            const_cast<const ManyBodyPotential_TersoffTest::configuration_t *>(
+                ManyBodyPotential_TersoffTest::CurrentConfiguration
+                )->begin(), firstiter),
+        std::distance( // enforce const_iterator return from begin()
+            const_cast<const ManyBodyPotential_TersoffTest::configuration_t *>(
+                ManyBodyPotential_TersoffTest::CurrentConfiguration
+                )->begin(), iter)
+    );
+    // jk
+    args[1].distance = seconditer->distance(*iter);
+    args[1].indices = std::make_pair(
+        std::distance( // enforce const_iterator return from begin()
+            const_cast<const ManyBodyPotential_TersoffTest::configuration_t *>(
+                ManyBodyPotential_TersoffTest::CurrentConfiguration
+                )->begin(), seconditer),
+        std::distance( // enforce const_iterator return from begin()
+            const_cast<const ManyBodyPotential_TersoffTest::configuration_t *>(
+                ManyBodyPotential_TersoffTest::CurrentConfiguration
+                )->begin(), iter)
+    );
+    result.push_back(args);
+  }
+  return result;
+}
+
+void ManyBodyPotential_TersoffTest::setUp()
+{
+  // failing asserts should be thrown
+  ASSERT_DO(Assert::Throw);
+
+  setVerbosity(10);
+
+  // we use parameters from tremolo/defaults/tersoff/tersoff.potentials
+  // [Tersoff, '89]
+  params.resize(ManyBodyPotential_Tersoff::MAXPARAMS,0.);
+  params[ManyBodyPotential_Tersoff::R] = 1.800000e+00;
+  params[ManyBodyPotential_Tersoff::S] = 2.100000e+00;
+  params[ManyBodyPotential_Tersoff::A] = 1.393600e+03;
+  params[ManyBodyPotential_Tersoff::B] = 3.467000e+02;
+  params[ManyBodyPotential_Tersoff::lambda] = 3.487900e+00;
+  params[ManyBodyPotential_Tersoff::mu] = 2.211900e+00;
+  params[ManyBodyPotential_Tersoff::lambda3] = 0.;
+  params[ManyBodyPotential_Tersoff::alpha] = 0.;
+  params[ManyBodyPotential_Tersoff::beta] = 1.572400e-07;
+  params[ManyBodyPotential_Tersoff::chi] = 1.;
+  params[ManyBodyPotential_Tersoff::omega] = 1.;
+  params[ManyBodyPotential_Tersoff::n] = 7.275100e-01;
+  params[ManyBodyPotential_Tersoff::c] = 3.804900e+04;
+  params[ManyBodyPotential_Tersoff::d] = 4.384000e+00;
+  params[ManyBodyPotential_Tersoff::h] = -5.705800e-01;
+
+  // initial configuration as in tremolo/default/tersoff/tersoff.data with
+  // Si renamed to C.
+  // create test configurations of 5 C atoms, constructed via:
+  // for file in tersoff.vis.00?0.xyz; do
+  // echo -e "\t{\n\t\tconfiguration_t positions;\n\t\tpositions +="
+  // tail -n 5 $file | awk -F " " {'print "\t\t\t\tVector("$2","$3","$4")"(NR!=5 ? "," : ";")'}
+  // echo -e -n "\t}\n"
+  // done
+  {
+    configuration_t positions;
+    positions +=
+        Vector(1.000000e+01,1.100000e+01,1.100000e+01),
+        Vector(1.120000e+01,1.000000e+01,1.000000e+01),
+        Vector(8.800000e+00,1.000000e+01,1.000000e+01),
+        Vector(1.000000e+01,1.200000e+01,1.000000e+01),
+        Vector(1.000000e+01,1.100000e+01,1.240000e+01);
+    configurations.push_back(positions);
+  }
+  {
+    configuration_t positions;
+    positions +=
+        Vector(1.000000e+01,1.099315e+01,1.099482e+01),
+        Vector(1.119779e+01,1.000179e+01,1.000235e+01),
+        Vector(8.802208e+00,1.000179e+01,1.000235e+01),
+        Vector(1.000000e+01,1.200262e+01,9.999421e+00),
+        Vector(1.000000e+01,1.100066e+01,1.240107e+01);
+    configurations.push_back(positions);
+  }
+  {
+    configuration_t positions;
+    positions +=
+        Vector(1.000000e+01,1.097354e+01,1.098018e+01),
+        Vector(1.119164e+01,1.000675e+01,1.000902e+01),
+        Vector(8.808358e+00,1.000675e+01,1.000902e+01),
+        Vector(1.000000e+01,1.201036e+01,9.997816e+00),
+        Vector(1.000000e+01,1.100260e+01,1.240397e+01);
+    configurations.push_back(positions);
+  }
+  {
+    configuration_t positions;
+    positions +=
+        Vector(1.000000e+01,1.094419e+01,1.095884e+01),
+        Vector(1.118308e+01,1.001364e+01,1.001884e+01),
+        Vector(8.816924e+00,1.001364e+01,1.001884e+01),
+        Vector(1.000000e+01,1.202283e+01,9.995566e+00),
+        Vector(1.000000e+01,1.100570e+01,1.240790e+01);
+    configurations.push_back(positions);
+  }
+  {
+    configuration_t positions;
+    positions +=
+        Vector(1.000000e+01,1.090791e+01,1.093293e+01),
+        Vector(1.117318e+01,1.002158e+01,1.003102e+01),
+        Vector(8.826818e+00,1.002158e+01,1.003102e+01),
+        Vector(1.000000e+01,1.203924e+01,9.993210e+00),
+        Vector(1.000000e+01,1.100969e+01,1.241182e+01);
+    configurations.push_back(positions);
+  }
+  {
+    configuration_t positions;
+    positions +=
+        Vector(1.000000e+01,1.086664e+01,1.090321e+01),
+        Vector(1.116216e+01,1.003043e+01,1.004546e+01),
+        Vector(8.837839e+00,1.003043e+01,1.004546e+01),
+        Vector(1.000000e+01,1.205848e+01,9.991167e+00),
+        Vector(1.000000e+01,1.101403e+01,1.241470e+01);
+    configurations.push_back(positions);
+  }
+  {
+    configuration_t positions;
+    positions +=
+        Vector(1.000000e+01,1.082297e+01,1.087033e+01),
+        Vector(1.115036e+01,1.003997e+01,1.006213e+01),
+        Vector(8.849644e+00,1.003997e+01,1.006213e+01),
+        Vector(1.000000e+01,1.207923e+01,9.989652e+00),
+        Vector(1.000000e+01,1.101786e+01,1.241577e+01);
+    configurations.push_back(positions);
+  }
+  {
+    configuration_t positions;
+    positions +=
+        Vector(1.000000e+01,1.077905e+01,1.083482e+01),
+        Vector(1.113831e+01,1.004988e+01,1.008085e+01),
+        Vector(8.861694e+00,1.004988e+01,1.008085e+01),
+        Vector(1.000000e+01,1.210048e+01,9.989022e+00),
+        Vector(1.000000e+01,1.102071e+01,1.241446e+01);
+    configurations.push_back(positions);
+  }
+  {
+    configuration_t positions;
+    positions +=
+        Vector(1.000000e+01,1.073683e+01,1.079745e+01),
+        Vector(1.112678e+01,1.005973e+01,1.010129e+01),
+        Vector(8.873218e+00,1.005973e+01,1.010129e+01),
+        Vector(1.000000e+01,1.212139e+01,9.989594e+00),
+        Vector(1.000000e+01,1.102232e+01,1.241038e+01);
+    configurations.push_back(positions);
+  }
+  {
+    configuration_t positions;
+    positions +=
+        Vector(1.000000e+01,1.069834e+01,1.075920e+01),
+        Vector(1.111685e+01,1.006891e+01,1.012296e+01),
+        Vector(8.883151e+00,1.006891e+01,1.012296e+01),
+        Vector(1.000000e+01,1.214131e+01,9.991602e+00),
+        Vector(1.000000e+01,1.102252e+01,1.240327e+01);
+    configurations.push_back(positions);
+  }
+
+  // cut from tersoff.etot via:
+  // for i in `seq 0 1 9`; do
+  // grep $i.000000e-01 tersoff.etot | awk -F" " {'print "\t\t\t\t"$2","'}
+  // done
+  // (though timestep 0 is missing and is added manually)
+  output +=
+      -1.333927e+01,
+      -1.359628e+01,
+      -1.420701e+01,
+      -1.479974e+01,
+      -1.537942e+01,
+      -1.584603e+01,
+      -1.615832e+01,
+      -1.630598e+01,
+      -1.626654e+01,
+      -1.603360e+01;
+
+  CPPUNIT_ASSERT_EQUAL( output.size(), configurations.size() );
+}
+
+void ManyBodyPotential_TersoffTest::tearDown()
+{
+  configurations.clear();
+  CurrentConfiguration = NULL;
+}
+
+/** UnitTest for operator()
+ */
+void ManyBodyPotential_TersoffTest::operatorTest()
+{
+  boost::function<
+      std::vector<FunctionModel::arguments_t>(const argument_t &, const double)
+    > fct =
+      triplefunction;
+  ManyBodyPotential_Tersoff tersoff(fct);
+  tersoff.setParameters(params);
+  for (size_t index = 0; index < configurations.size(); ++index) {
+    CurrentConfiguration = &(configurations[index]);
+    double temp = 0.;
+    for (size_t i=0; i < CurrentConfiguration->size(); ++i)
+      for (size_t j=0; j < CurrentConfiguration->size(); ++j) {
+        if (i == j)
+          continue;
+        argument_t arg;
+        arg.indices = std::make_pair(i,j);
+        arg.distance = (*CurrentConfiguration)[i].distance((*CurrentConfiguration)[j]);
+        FunctionModel::arguments_t args(1,arg);
+        const ManyBodyPotential_Tersoff::results_t res = tersoff(args);
+        temp += res[0];
+      }
+    // this little precision is because tremolo does seem to handle coordinates
+    // a little differently than we do, the precise difference in the x coordinate
+    // of the first and second position for the second configuration is easy to
+    // see as 1.119779. However, tremolo obtains 1.1197792 for unknown reasons.
+    // Maybe, there is some float floating around in the code ...
+    CPPUNIT_ASSERT(
+        Helpers::isEqual(
+            output[index],
+            .5*temp,
+            1.e-4/std::numeric_limits<double>::epsilon()
+        )
+    );
+  }
+}
+
+/** UnitTest for derivative()
+ */
+void ManyBodyPotential_TersoffTest::derivativeTest()
+{
+//  boost::function<
+//      std::vector<FunctionModel::arguments_t>(const argument_t &, const double)
+//    > fct =
+//      triplefunction;
+//  ManyBodyPotential_Tersoff tersoff(fct);
+//  tersoff.setParameters(params);
+//  CPPUNIT_ASSERT(
+//      Helpers::isEqual(
+//          0.,
+//          tersoff.derivative(
+//              FunctionModel::arguments_t(1,argument_t(1.))
+//          )[0]
+//      )
+//  );
+}
+
+
+/** UnitTest for parameter_derivative()
+ */
+void ManyBodyPotential_TersoffTest::parameter_derivativeTest()
+{
+//  boost::function<
+//      std::vector<FunctionModel::arguments_t>(const argument_t &, const double)
+//    > fct =
+//      triplefunction;
+//  ManyBodyPotential_Tersoff tersoff(fct);
+//  tersoff.setParameters(params);
+//  CPPUNIT_ASSERT(
+//      Helpers::isEqual(
+//          0.,
+//          tersoff.parameter_derivative(
+//              FunctionModel::arguments_t(1,argument_t(1.)),
+//              0
+//          )[0]
+//      )
+//  );
+//  CPPUNIT_ASSERT(
+//      Helpers::isEqual(
+//          0.,
+//          tersoff.parameter_derivative(
+//              FunctionModel::arguments_t(1,argument_t(1.)),
+//              1
+//          )[0]
+//      )
+//  );
+//  CPPUNIT_ASSERT(
+//      Helpers::isEqual(
+//          1.,
+//          tersoff.parameter_derivative(
+//              FunctionModel::arguments_t(1,argument_t(1.)),
+//              2
+//          )[0]
+//      )
+//  );
+}
Index: src/Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.hpp
===================================================================
--- src/Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.hpp	(revision ffc3688dd044a0c874c10a37de01139caff883c1)
+++ src/Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.hpp	(revision ffc3688dd044a0c874c10a37de01139caff883c1)
@@ -0,0 +1,54 @@
+/*
+ * ManyBodyPotential_TersoffUnitTest.hpp
+ *
+ *  Created on: Oct 04, 2012
+ *      Author: heber
+ */
+
+#ifndef MANYBODYPOTENTIAL_TERSOFFUNITTEST_HPP_
+#define MANYBODYPOTENTIAL_TERSOFFUNITTEST_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include <vector>
+
+#include "LinearAlgebra/Vector.hpp"
+
+#include "FunctionApproximation/FunctionModel.hpp"
+
+class ManyBodyPotential_Tersoff;
+
+/********************************************** Test classes **************************************/
+
+class ManyBodyPotential_TersoffTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( ManyBodyPotential_TersoffTest) ;
+    CPPUNIT_TEST ( operatorTest );
+    CPPUNIT_TEST ( derivativeTest );
+    CPPUNIT_TEST ( parameter_derivativeTest );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+      void setUp();
+      void tearDown();
+      void operatorTest();
+      void derivativeTest();
+      void parameter_derivativeTest();
+
+      typedef std::vector<Vector> configuration_t;
+      typedef std::vector<configuration_t> configurations_t;
+      static configurations_t configurations;
+      static configuration_t *CurrentConfiguration;
+private:
+      FunctionModel::parameters_t params;
+      typedef std::vector<double> data_t;
+      data_t output;
+};
+
+#endif /* MANYBODYPOTENTIAL_TERSOFFUNITTEST_HPP_ */
