Index: src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp
===================================================================
--- src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp	(revision 550f2a98ea5244ef1fc747f6debbc4635e941f7c)
@@ -184,29 +184,153 @@
 ManyBodyPotential_Tersoff::results_t
 ManyBodyPotential_Tersoff::operator()(
-    const arguments_t &arguments
+    const list_of_arguments_t &listarguments
     ) const
 {
 //  Info info(__func__);
-  double result = 0.;
-  for(arguments_t::const_iterator argiter = arguments.begin();
-      argiter != arguments.end();
-      ++argiter) {
-    const argument_t &r_ij = *argiter;
-    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-        arguments_t(1, r_ij), getParticleTypes()),
-        "ManyBodyPotential_Tersoff::operator() - types don't match with ones in arguments.");
- 
-    const double cutoff = function_cutoff(r_ij.distance);
-    const double temp = (cutoff == 0.) ?
-        0. :
-        cutoff * (
-            function_prefactor(
-                alpha,
-                function_eta(r_ij))
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    for(arguments_t::const_iterator argiter = arguments.begin();
+        argiter != arguments.end();
+        ++argiter) {
+      const argument_t &r_ij = *argiter;
+      ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+          arguments_t(1, r_ij), getParticleTypes()),
+          "ManyBodyPotential_Tersoff::operator() - types don't match with ones in arguments.");
+
+      const double cutoff = function_cutoff(r_ij.distance);
+      const double temp = (cutoff == 0.) ?
+          0. :
+          cutoff * (
+              function_prefactor(
+                  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)
+          );
+      result += temp;
+    }
+  }
+//  LOG(2, "DEBUG: operator()(" << r_ij.distance << ") = " << result);
+  return results_t(1, result);
+}
+
+ManyBodyPotential_Tersoff::derivative_components_t
+ManyBodyPotential_Tersoff::derivative(
+    const list_of_arguments_t &listarguments
+    ) const
+{
+//  Info info(__func__);
+  return derivative_components_t();
+}
+
+ManyBodyPotential_Tersoff::results_t
+ManyBodyPotential_Tersoff::parameter_derivative(
+    const list_of_arguments_t &listarguments,
+    const size_t index
+    ) const
+{
+//  Info info(__func__);
+//  ASSERT( arguments.size() == 1,
+//      "ManyBodyPotential_Tersoff::parameter_derivative() - requires exactly one argument.");
+
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    for(arguments_t::const_iterator argiter = arguments.begin();
+        argiter != arguments.end();
+        ++argiter) {
+      const argument_t &r_ij = *argiter;
+      ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+          arguments_t(1, r_ij), getParticleTypes()),
+          "ManyBodyPotential_Tersoff::operator() - types don't match with ones in arguments.");
+
+    switch (index) {
+//      case R:
+//      {
+//        result += 0.;
+//        break;
+//      }
+//      case S:
+//      {
+//        result += 0.;
+//        break;
+//      }
+      case A:
+      {
+        const double cutoff = function_cutoff(r_ij.distance);
+        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);
+        break;
+      }
+      case B:
+      {
+        const double cutoff = function_cutoff(r_ij.distance);
+        result += (cutoff == 0.) ?
+            0. :
+            cutoff * function_prefactor(
+                params[beta],
+                function_zeta(r_ij))
             * function_smoother(
-                params[A],
-                params[lambda],
-                r_ij.distance)
-            +
+                -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];
+        break;
+      }
+      case lambda:
+      {
+        const double cutoff = function_cutoff(r_ij.distance);
+        result += (cutoff == 0.) ?
+            0. :
+            -r_ij.distance * cutoff *
+                function_prefactor(
+                    alpha,
+                    function_eta(r_ij))
+                * function_smoother(
+                    params[A],
+                    params[lambda],
+                    r_ij.distance);
+        break;
+      }
+      case mu:
+      {
+        const double cutoff = function_cutoff(r_ij.distance);
+        result += (cutoff == 0.) ?
+            0. :
+            -r_ij.distance * cutoff *(
             function_prefactor(
                 params[beta],
@@ -217,272 +341,156 @@
                 r_ij.distance)
         );
-    result += temp;
-  }
-//  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,
-//      "ManyBodyPotential_Tersoff::parameter_derivative() - requires exactly one argument.");
-
-  double result = 0.;
-  for(arguments_t::const_iterator argiter = arguments.begin();
-      argiter != arguments.end();
-      ++argiter) {
-    const argument_t &r_ij = *argiter;
-    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-        arguments_t(1, r_ij), getParticleTypes()),
-        "ManyBodyPotential_Tersoff::operator() - types don't match with ones in arguments.");
-
-  switch (index) {
-//    case R:
-//    {
-//      result += 0.;
-//      break;
-//    }
-//    case S:
-//    {
-//      result += 0.;
-//      break;
-//    }
-    case A:
-    {
-      const double cutoff = function_cutoff(r_ij.distance);
-      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);
-      break;
+        break;
+      }
+//      case lambda3:
+//      {
+//        result += 0.;
+//        break;
+//      }
+//      case alpha:
+//      {
+//        const double temp =
+//            pow(alpha*function_eta(r_ij), params[n]);
+//        const double cutoff = function_cutoff(r_ij.distance);
+//        result += (cutoff == 0.) || (alpha == 0. )?
+//            0. :
+//            function_smoother(
+//                params[A],
+//                params[lambda],
+//                r_ij.distance)
+//            * (-.5) * alpha * (temp/alpha)
+//            / (1. + temp)
+//            ;
+//        break;
+//      }
+//      case chi:
+//      {
+//        result += 0.;
+//        break;
+//      }
+//      case omega:
+//      {
+//        result += 0.;
+//        break;
+//      }
+      case beta:
+      {
+        const double temp =
+            pow(params[beta]*function_zeta(r_ij), params[n]);
+        const double cutoff = function_cutoff(r_ij.distance);
+        result += (cutoff == 0.) || (params[beta] == 0. )?
+            0. : cutoff *
+            function_smoother(
+                -params[B],
+                params[mu],
+                r_ij.distance)
+            * (-.5)
+            * function_prefactor(
+                params[beta],
+                function_zeta(r_ij))
+            * (temp/params[beta])
+            / (1. + temp)
+            ;
+        break;
+      }
+      case n:
+      {
+        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 tempres = ((cutoff == 0.) || (zeta == 0.)) ? // zeta must be caught if zero due to log
+            0. : .5 * cutoff *
+            function_smoother(
+                -params[B],
+                params[mu],
+                r_ij.distance)
+            * function_prefactor(
+                params[beta],
+                function_zeta(r_ij))
+            * ( log(1.+temp)/(params[n]*params[n]) - temp
+                * (log(function_zeta(r_ij)) + log(params[beta]))
+                /(params[n]*(1.+temp)))
+            ;
+//        if (tempres != tempres)
+//    LOG(2, "DEBUG: tempres is NaN.");
+//        LOG(2, "DEBUG: Adding " << tempres << " for p.d. w.r.t n, temp=" << temp << ", cutoff=" << cutoff);
+        result += tempres;
+        break;
+      }
+      case c:
+      {
+        const double zeta = function_zeta(r_ij);
+        if (zeta == 0.)
+          break;
+        const double temp =
+            pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
+        const double cutoff = function_cutoff(r_ij.distance);
+        const double tempres = (cutoff == 0.) ?
+            0. : cutoff *
+            function_smoother(
+                -params[B],
+                params[mu],
+                r_ij.distance)
+            * function_prefactor(
+                params[beta],
+                zeta)
+             * (-1.) * temp / (1.+temp*zeta);
+        double factor = function_derivative_c(r_ij);
+        result += tempres*factor;
+        if (result != result)
+          ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
+        break;
+      }
+      case d:
+      {
+        const double zeta = function_zeta(r_ij);
+        const double temp =
+            pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
+        const double cutoff = function_cutoff(r_ij.distance);
+        const double tempres = (cutoff == 0.) ?
+            0. : cutoff *
+            function_smoother(
+                -params[B],
+                params[mu],
+                r_ij.distance)
+            * function_prefactor(
+                params[beta],
+                zeta)
+             * (-1.) * temp / (1.+temp*zeta);
+        double factor = function_derivative_d(r_ij);
+        result += tempres*factor;
+        if (result != result)
+          ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
+        break;
+      }
+      case h:
+      {
+        const double zeta = function_zeta(r_ij);
+        const double temp =
+            pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
+        const double cutoff = function_cutoff(r_ij.distance);
+        const double tempres = (cutoff == 0.) ?
+            0. : cutoff *
+            function_smoother(
+                -params[B],
+                params[mu],
+                r_ij.distance)
+            * function_prefactor(
+                params[beta],
+                zeta)
+             * (-1.) * temp / (1.+temp*zeta);
+        double factor = function_derivative_h(r_ij);
+        result += tempres*factor;
+        if (result != result)
+          ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
+        break;
+      }
+      default:
+        ASSERT(0, "ManyBodyPotential_Tersoff::parameter_derivative() - derivative to unknown parameter desired.");
+        break;
     }
-    case B:
-    {
-      const double cutoff = function_cutoff(r_ij.distance);
-      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];
-      break;
+    if (result != result)
+      ELOG(1, "result is NaN.");
     }
-    case lambda:
-    {
-      const double cutoff = function_cutoff(r_ij.distance);
-      result += (cutoff == 0.) ?
-          0. :
-          -r_ij.distance * cutoff *
-              function_prefactor(
-                  alpha,
-                  function_eta(r_ij))
-              * function_smoother(
-                  params[A],
-                  params[lambda],
-                  r_ij.distance);
-      break;
-    }
-    case mu:
-    {
-      const double cutoff = function_cutoff(r_ij.distance);
-      result += (cutoff == 0.) ?
-          0. :
-          -r_ij.distance * cutoff *(
-          function_prefactor(
-              params[beta],
-              function_zeta(r_ij))
-          * function_smoother(
-              -params[B],
-              params[mu],
-              r_ij.distance)
-      );
-      break;
-    }
-//    case lambda3:
-//    {
-//      result += 0.;
-//      break;
-//    }
-//    case alpha:
-//    {
-//      const double temp =
-//          pow(alpha*function_eta(r_ij), params[n]);
-//      const double cutoff = function_cutoff(r_ij.distance);
-//      result += (cutoff == 0.) || (alpha == 0. )?
-//          0. :
-//          function_smoother(
-//              params[A],
-//              params[lambda],
-//              r_ij.distance)
-//          * (-.5) * alpha * (temp/alpha)
-//          / (1. + temp)
-//          ;
-//      break;
-//    }
-//    case chi:
-//    {
-//      result += 0.;
-//      break;
-//    }
-//    case omega:
-//    {
-//      result += 0.;
-//      break;
-//    }
-    case beta:
-    {
-      const double temp =
-          pow(params[beta]*function_zeta(r_ij), params[n]);
-      const double cutoff = function_cutoff(r_ij.distance);
-      result += (cutoff == 0.) || (params[beta] == 0. )?
-          0. : cutoff *
-          function_smoother(
-              -params[B],
-              params[mu],
-              r_ij.distance)
-          * (-.5)
-          * function_prefactor(
-              params[beta],
-              function_zeta(r_ij))
-          * (temp/params[beta])
-          / (1. + temp)
-          ;
-      break;
-    }
-    case n:
-    {
-      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 tempres = ((cutoff == 0.) || (zeta == 0.)) ? // zeta must be caught if zero due to log
-          0. : .5 * cutoff *
-          function_smoother(
-              -params[B],
-              params[mu],
-              r_ij.distance)
-          * function_prefactor(
-              params[beta],
-              function_zeta(r_ij))
-          * ( log(1.+temp)/(params[n]*params[n]) - temp
-              * (log(function_zeta(r_ij)) + log(params[beta]))
-              /(params[n]*(1.+temp)))
-          ;
-//      if (tempres != tempres)
-//	LOG(2, "DEBUG: tempres is NaN."); 
-//      LOG(2, "DEBUG: Adding " << tempres << " for p.d. w.r.t n, temp=" << temp << ", cutoff=" << cutoff);
-      result += tempres;
-      break;
-    }
-    case c:
-    {
-      const double zeta = function_zeta(r_ij);
-      if (zeta == 0.)
-        break;
-      const double temp =
-          pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
-      const double cutoff = function_cutoff(r_ij.distance);
-      const double tempres = (cutoff == 0.) ?
-          0. : cutoff *
-          function_smoother(
-              -params[B],
-              params[mu],
-              r_ij.distance)
-          * function_prefactor(
-              params[beta],
-              zeta)
-           * (-1.) * temp / (1.+temp*zeta);
-      double factor = function_derivative_c(r_ij);
-      result += tempres*factor;
-      if (result != result)
-        ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
-      break;
-    }
-    case d:
-    {
-      const double zeta = function_zeta(r_ij);
-      const double temp =
-          pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
-      const double cutoff = function_cutoff(r_ij.distance);
-      const double tempres = (cutoff == 0.) ?
-          0. : cutoff *
-          function_smoother(
-              -params[B],
-              params[mu],
-              r_ij.distance)
-          * function_prefactor(
-              params[beta],
-              zeta)
-           * (-1.) * temp / (1.+temp*zeta);
-      double factor = function_derivative_d(r_ij);
-      result += tempres*factor;
-      if (result != result)
-        ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
-      break;
-    }
-    case h:
-    {
-      const double zeta = function_zeta(r_ij);
-      const double temp =
-          pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
-      const double cutoff = function_cutoff(r_ij.distance);
-      const double tempres = (cutoff == 0.) ?
-          0. : cutoff *
-          function_smoother(
-              -params[B],
-              params[mu],
-              r_ij.distance)
-          * function_prefactor(
-              params[beta],
-              zeta)
-           * (-1.) * temp / (1.+temp*zeta);
-      double factor = function_derivative_h(r_ij);
-      result += tempres*factor;
-      if (result != result)
-        ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
-      break;
-    }
-    default:
-      ASSERT(0, "ManyBodyPotential_Tersoff::parameter_derivative() - derivative to unknown parameter desired.");
-      break;
-  }
-  if (result != result)
-    ELOG(1, "result is NaN.");
   }
   return results_t(1,-result);
