Index: src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp
===================================================================
--- src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp	(revision f904d5f531c1e72c1dd49182c9acbc071c10c0ac)
+++ src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp	(revision ca8d8233e18d2fac0946021d984e388e215e9212)
@@ -170,5 +170,22 @@
     case A:
     {
-      const double result = 0.;
+      const double cutoff = function_cutoff(r_ij.distance);
+      const double result = (cutoff == 0.) ?
+          0. :
+          cutoff *
+              function_prefactor(
+                  alpha,
+                  function_eta(r_ij))
+              * function_smoother(
+                  1.,
+                  params[lambda],
+                  r_ij.distance);
+//          cutoff * function_prefactor(
+//              alpha,
+//              function_eta(r_ij))
+//          * function_smoother(
+//              1.,
+//              params[mu],
+//              r_ij.distance);
       return results_t(1, result);
       break;
@@ -176,5 +193,21 @@
     case B:
     {
-      const double result = 0.;
+      const double cutoff = function_cutoff(r_ij.distance);
+      const double result = (cutoff == 0.) ?
+          0. :
+          cutoff * function_prefactor(
+              params[beta],
+              function_zeta(r_ij))
+          * function_smoother(
+              -1.,
+              params[mu],
+              r_ij.distance);
+//          cutoff * function_prefactor(
+//              beta,
+//              function_zeta(r_ij))
+//          * function_smoother(
+//              -params[B],
+//              params[mu],
+//              r_ij.distance)/params[B];
       return results_t(1, result);
       break;
@@ -185,13 +218,12 @@
       const double result = (cutoff == 0.) ?
           0. :
-          -r_ij.distance * cutoff * params[lambda] * (
-          function_prefactor(
-              alpha,
-              function_eta(r_ij))
-          * function_smoother(
-              params[A],
-              params[lambda],
-              r_ij.distance)
-          );
+          -r_ij.distance * cutoff *
+              function_prefactor(
+                  alpha,
+                  function_eta(r_ij))
+              * function_smoother(
+                  params[A],
+                  params[lambda],
+                  r_ij.distance);
       return results_t(1, result);
       break;
@@ -228,5 +260,5 @@
 //          0. :
 //          function_smoother(
-//              -params[A],
+//              params[A],
 //              params[lambda],
 //              r_ij.distance)
@@ -294,9 +326,11 @@
     {
       const double zeta = function_zeta(r_ij);
+      if (zeta == 0.)
+        break;
       const double temp =
-          pow(params[beta]*zeta, params[n]);
-      const double cutoff = function_cutoff(r_ij.distance);
-      const double result = (cutoff == 0.) ?
-          0. : - cutoff *
+          pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
+      const double cutoff = function_cutoff(r_ij.distance);
+      const double result = (cutoff == 0.) ?
+          0. : cutoff *
           function_smoother(
               -params[B],
@@ -306,18 +340,6 @@
               params[beta],
               zeta)
-          * (temp/zeta)/(1+(temp/zeta));
-      double factor = 0.;
-      std::vector<arguments_t> triples = triplefunction(r_ij, S);
-      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];
-        factor +=
-            function_cutoff(r_ik.distance)
-            * (params[c]/Helpers::pow(params[c],2))
-            * (1. - 1./(1.+(params[h]-function_theta(r_ij.distance, r_ik.distance, r_jk.distance))/Helpers::pow(params[c],2)));
-      }
+           * (-1.) * temp / (1.+temp*zeta);
+      double factor = function_derivative_c(r_ij);
       return results_t(1, result*factor);
       break;
@@ -327,8 +349,8 @@
       const double zeta = function_zeta(r_ij);
       const double temp =
-          pow(params[beta]*zeta, params[n]);
-      const double cutoff = function_cutoff(r_ij.distance);
-      const double result = (cutoff == 0.) ?
-          0. : - cutoff *
+          pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
+      const double cutoff = function_cutoff(r_ij.distance);
+      const double result = (cutoff == 0.) ?
+          0. : cutoff *
           function_smoother(
               -params[B],
@@ -338,20 +360,6 @@
               params[beta],
               zeta)
-          * (temp/zeta)/(1+(temp/zeta));
-      double factor = 0.;
-      std::vector<arguments_t> triples = triplefunction(r_ij, S);
-      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];
-        factor +=
-            function_cutoff(r_ik.distance)
-            * (-Helpers::pow(params[c],2)/Helpers::pow(params[d],3)
-                - Helpers::pow(params[c],2) * params[d]
-                 / Helpers::pow(Helpers::pow(params[d],2) + Helpers::pow(params[h] - function_theta(r_ij.distance, r_ik.distance, r_jk.distance),2),2)
-            );
-      }
+           * (-1.) * temp / (1.+temp*zeta);
+      double factor = function_derivative_d(r_ij);
       return results_t(1, result*factor);
       break;
@@ -361,8 +369,8 @@
       const double zeta = function_zeta(r_ij);
       const double temp =
-          pow(params[beta]*zeta, params[n]);
-      const double cutoff = function_cutoff(r_ij.distance);
-      const double result = (cutoff == 0.) ?
-          0. : - cutoff *
+          pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
+      const double cutoff = function_cutoff(r_ij.distance);
+      const double result = (cutoff == 0.) ?
+          0. : cutoff *
           function_smoother(
               -params[B],
@@ -372,20 +380,6 @@
               params[beta],
               zeta)
-          * (temp/zeta)/(1+(temp/zeta));
-      double factor = 0.;
-      std::vector<arguments_t> triples = triplefunction(r_ij, S);
-      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];
-        const double tempangle = params[h] - function_theta(r_ij.distance, r_ik.distance, r_jk.distance);
-        factor +=
-            function_cutoff(r_ik.distance)
-            * (-1.*Helpers::pow(params[c],2)*tempangle)
-            / Helpers::pow(Helpers::pow(params[d],2)+Helpers::pow(tempangle,2),2)
-            ;
-      }
+           * (-1.) * temp / (1.+temp*zeta);
+      double factor = function_derivative_h(r_ij);
       return results_t(1, result*factor);
       break;
@@ -395,4 +389,77 @@
   }
   return results_t(1, 0.);
+}
+
+ManyBodyPotential_Tersoff::result_t
+ManyBodyPotential_Tersoff::function_derivative_c(
+    const argument_t &r_ij
+  ) const
+{
+  double result = 0.;
+  std::vector<arguments_t> triples = triplefunction(r_ij, S);
+  for (std::vector<arguments_t>::const_iterator iter = triples.begin();
+      iter != triples.end(); ++iter) {
+    ASSERT( iter->size() == 2,
+        "ManyBodyPotential_Tersoff::function_derivative_c() - the triples result must contain exactly two distances.");
+    const argument_t &r_ik = (*iter)[0];
+    const argument_t &r_jk = (*iter)[1];
+    const double tempangle = params[h] - function_theta(r_ij.distance, r_ik.distance, r_jk.distance);
+    const double cutoff = function_cutoff(r_ik.distance);
+    result += (cutoff == 0.) ?
+        0. : cutoff * omega * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3)) * (
+            params[c]/Helpers::pow(params[d],2)
+            - params[c] / ( Helpers::pow(params[d],2) + Helpers::pow(tempangle,2) )
+        );
+  }
+  return result;
+}
+
+ManyBodyPotential_Tersoff::result_t
+ManyBodyPotential_Tersoff::function_derivative_d(
+    const argument_t &r_ij
+  ) const
+{
+  double result = 0.;
+  std::vector<arguments_t> triples = triplefunction(r_ij, S);
+  for (std::vector<arguments_t>::const_iterator iter = triples.begin();
+      iter != triples.end(); ++iter) {
+    ASSERT( iter->size() == 2,
+        "ManyBodyPotential_Tersoff::function_derivative_d() - the triples result must contain exactly two distances.");
+    const argument_t &r_ik = (*iter)[0];
+    const argument_t &r_jk = (*iter)[1];
+    const double tempangle = params[h] - function_theta(r_ij.distance, r_ik.distance, r_jk.distance);
+    const double cutoff = function_cutoff(r_ik.distance);
+    result += (cutoff == 0.) ?
+        0. : cutoff * omega  * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3)) * (
+          - Helpers::pow(params[c],2)/Helpers::pow(params[d],3)
+          + Helpers::pow(params[c],2) * params[d]
+            / Helpers::pow(Helpers::pow(params[d],2) + Helpers::pow(tempangle,2),2)
+        );
+  }
+  return result;
+}
+
+ManyBodyPotential_Tersoff::result_t
+ManyBodyPotential_Tersoff::function_derivative_h(
+    const argument_t &r_ij
+  ) const
+{
+  double result = 0.;
+  std::vector<arguments_t> triples = triplefunction(r_ij, S);
+  for (std::vector<arguments_t>::const_iterator iter = triples.begin();
+      iter != triples.end(); ++iter) {
+    ASSERT( iter->size() == 2,
+        "ManyBodyPotential_Tersoff::function_derivative_h() - the triples result must contain exactly two distances.");
+    const argument_t &r_ik = (*iter)[0];
+    const argument_t &r_jk = (*iter)[1];
+    const double tempangle = params[h] - function_theta(r_ij.distance, r_ik.distance, r_jk.distance);
+    const double cutoff = function_cutoff(r_ik.distance);
+    result += (cutoff == 0.) ?
+        0. : cutoff * omega  * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3)) * (
+          ( Helpers::pow(params[c],2)*tempangle )
+            / Helpers::pow(Helpers::pow(params[d],2) + Helpers::pow(tempangle,2),2)
+        );
+  }
+  return result;
 }
 
@@ -516,14 +583,12 @@
 {
 //  Info info(__func__);
-  const double angle = function_theta(r_ij, r_ik, r_jk);
-//  LOG(2, "DEBUG: cos(theta)= " << angle/divisor);
   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,2));
-
-//    LOG(2, "DEBUG: function_angle(" << r_ij << "," << r_ik << "," << r_jk << ") = " << result);
-  return result;
-}
-
+          Helpers::pow(params[h] - function_theta(r_ij, r_ik, r_jk),2));
+
+//  LOG(2, "DEBUG: function_angle(" << r_ij << "," << r_ik << "," << r_jk << ") = " << result);
+  return result;
+}
+
Index: src/Potentials/Specifics/ManyBodyPotential_Tersoff.hpp
===================================================================
--- src/Potentials/Specifics/ManyBodyPotential_Tersoff.hpp	(revision f904d5f531c1e72c1dd49182c9acbc071c10c0ac)
+++ src/Potentials/Specifics/ManyBodyPotential_Tersoff.hpp	(revision ca8d8233e18d2fac0946021d984e388e215e9212)
@@ -186,4 +186,20 @@
 
 private:
+  result_t
+  function_derivative_c(
+      const argument_t &r_ij
+    ) const;
+
+  result_t
+  function_derivative_d(
+      const argument_t &r_ij
+    ) const;
+
+  result_t
+  function_derivative_h(
+      const argument_t &r_ij
+    ) const;
+
+public:
   enum parameter_enum_t {
     A,
@@ -223,8 +239,9 @@
   void setParameters(const parameters_t &_params)
   {
-    ASSERT( _params.size() == getParameterDimension(),
-        "ManyBodyPotential_Tersoff::setParameters() - we need exactly "
+    ASSERT( _params.size() <= getParameterDimension(),
+        "ManyBodyPotential_Tersoff::setParameters() - we need not more than "
         +toString(getParameterDimension())+" parameters.");
-    params = _params;
+    for (size_t i=0; i< _params.size(); ++i)
+      params[i] = _params[i];
   }
 
