Index: src/Potentials/Specifics/ConstantPotential.cpp
===================================================================
--- src/Potentials/Specifics/ConstantPotential.cpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/ConstantPotential.cpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -112,53 +112,69 @@
 ConstantPotential::results_t
 ConstantPotential::operator()(
-    const arguments_t &arguments
+    const list_of_arguments_t &listarguments
     ) const
 {
-  ASSERT( arguments.size() == 0,
-      "ConstantPotential::operator() - requires no argument.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "ConstantPotential::operator() - types don't match with ones in arguments.");
-  const result_t result = params[energy_offset];
-  return std::vector<result_t>(1, result);
+  // directly set result as energy constant as independent of arg list
+  result_t result = params[energy_offset];
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    ASSERT( arguments.size() == 0,
+        "ConstantPotential::operator() - requires no argument.");
+    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+        arguments, getParticleTypes()),
+        "ConstantPotential::operator() - types don't match with ones in arguments.");
+    result += 0.;
+  }
+  return results_t(1, result);
 }
 
 ConstantPotential::derivative_components_t
 ConstantPotential::derivative(
-    const arguments_t &arguments
+    const list_of_arguments_t &listarguments
     ) const
 {
-  ASSERT( arguments.size() == 0,
-      "ConstantPotential::operator() - requires no argument.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "ConstantPotential::operator() - types don't match with ones in arguments.");
-  derivative_components_t result(1, 0.);
-  return result;
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    ASSERT( arguments.size() == 0,
+        "ConstantPotential::operator() - requires no argument.");
+    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+        arguments, getParticleTypes()),
+        "ConstantPotential::operator() - types don't match with ones in arguments.");
+    result += 0.;
+  }
+  return derivative_components_t(1, result);
 }
 
 ConstantPotential::results_t
 ConstantPotential::parameter_derivative(
-    const arguments_t &arguments,
+    const list_of_arguments_t &listarguments,
     const size_t index
     ) const
 {
-  ASSERT( arguments.size() == 0,
-      "ConstantPotential::parameter_derivative() - requires no argument.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "ConstantPotential::operator() - types don't match with ones in arguments.");
-  switch (index) {
-    case energy_offset:
-    {
-      // Maple result: 1
-      const result_t result = +1.;
-      return std::vector<result_t>(1, result);
-      break;
-    }
-    default:
-      break;
+  // Maple result: 1
+  result_t result = 1.;   // energy_offset
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    ASSERT( arguments.size() == 0,
+        "ConstantPotential::parameter_derivative() - requires no argument.");
+    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+        arguments, getParticleTypes()),
+        "ConstantPotential::operator() - types don't match with ones in arguments.");
+//    switch (index) {
+//      case energy_offset:
+//      {
+//        // Maple result: 1
+//        result = +1.;
+//        break;
+//      }
+//      default:
+//        break;
+//    }
   }
-  return std::vector<result_t>(1, 0.);
+  return results_t(1, result);
 }
 
Index: src/Potentials/Specifics/ConstantPotential.hpp
===================================================================
--- src/Potentials/Specifics/ConstantPotential.hpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/ConstantPotential.hpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -36,4 +36,5 @@
   friend class PotentialFactory;
   // some repeated typedefs to avoid ambiguities
+  typedef FunctionModel::list_of_arguments_t list_of_arguments_t;
   typedef FunctionModel::arguments_t arguments_t;
   typedef FunctionModel::result_t result_t;
@@ -90,25 +91,25 @@
   /** Evaluates the harmonic potential function for the given arguments.
    *
-   * @param arguments single distance
+   * @param listarguments empty list
    * @return value of the potential function
    */
-  results_t operator()(const arguments_t &arguments) const;
+  results_t operator()(const list_of_arguments_t &listarguments) const;
 
   /** Evaluates the derivative of the potential function.
    *
-   * @param arguments single distance
+   * @param listarguments empty list
    * @return vector with derivative with respect to the input degrees of freedom
    */
-  derivative_components_t derivative(const arguments_t &arguments) const;
+  derivative_components_t derivative(const list_of_arguments_t &listarguments) const;
 
   /** Evaluates the derivative of the function with the given \a arguments
    * with respect to a specific parameter indicated by \a index.
    *
-   * \param arguments set of arguments as input variables to the function
+   * \param listarguments empty list
    * \param index derivative of which parameter
    * \return result vector containing the derivative with respect to the given
    *         input
    */
-  results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;
+  results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const;
 
   /** Returns the functor that converts argument_s into the
Index: src/Potentials/Specifics/FourBodyPotential_Torsion.cpp
===================================================================
--- src/Potentials/Specifics/FourBodyPotential_Torsion.cpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/FourBodyPotential_Torsion.cpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -140,85 +140,99 @@
 FourBodyPotential_Torsion::results_t
 FourBodyPotential_Torsion::operator()(
-    const arguments_t &arguments
+    const list_of_arguments_t &listarguments
     ) const
 {
-  ASSERT( arguments.size() == getSpecificArgumentCount(),
-      "FourBodyPotential_Torsion::operator() - requires exactly three arguments.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments.");
-  const argument_t &r_ij = arguments[0]; // 01
-  const argument_t &r_ik = arguments[1]; // 02
-  const argument_t &r_il = arguments[2]; // 03
-  const argument_t &r_jk = arguments[3]; // 12
-  const argument_t &r_jl = arguments[4]; // 13
-  const argument_t &r_kl = arguments[5]; // 23
-  const result_t result =
-      params[spring_constant]
-             * Helpers::pow( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance], 2 );
-  return std::vector<result_t>(1, result);
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    ASSERT( arguments.size() == getSpecificArgumentCount(),
+        "FourBodyPotential_Torsion::operator() - requires exactly three arguments.");
+    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+        arguments, getParticleTypes()),
+        "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments.");
+    const argument_t &r_ij = arguments[0]; // 01
+    const argument_t &r_ik = arguments[1]; // 02
+    const argument_t &r_il = arguments[2]; // 03
+    const argument_t &r_jk = arguments[3]; // 12
+    const argument_t &r_jl = arguments[4]; // 13
+    const argument_t &r_kl = arguments[5]; // 23
+    result +=
+        params[spring_constant]
+               * Helpers::pow( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance)
+                   - params[equilibrium_distance], 2 );
+  }
+  return results_t(1, result);
 }
 
 FourBodyPotential_Torsion::derivative_components_t
 FourBodyPotential_Torsion::derivative(
-    const arguments_t &arguments
+    const list_of_arguments_t &listarguments
     ) const
 {
-  ASSERT( arguments.size() == getSpecificArgumentCount(),
-      "FourBodyPotential_Torsion::operator() - requires exactly three arguments.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments.");
-  derivative_components_t result;
-  const argument_t &r_ij = arguments[0]; // 01
-  const argument_t &r_ik = arguments[1]; // 02
-  const argument_t &r_il = arguments[2]; // 03
-  const argument_t &r_jk = arguments[3]; // 12
-  const argument_t &r_jl = arguments[4]; // 13
-  const argument_t &r_kl = arguments[5]; // 23
-  result.push_back( 2. * params[spring_constant] * ( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance]) );
-  ASSERT( result.size() == 1,
-      "FourBodyPotential_Torsion::operator() - we did not create exactly one component.");
-  return result;
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    ASSERT( arguments.size() == getSpecificArgumentCount(),
+        "FourBodyPotential_Torsion::operator() - requires exactly three arguments.");
+    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+        arguments, getParticleTypes()),
+        "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments.");
+    const argument_t &r_ij = arguments[0]; // 01
+    const argument_t &r_ik = arguments[1]; // 02
+    const argument_t &r_il = arguments[2]; // 03
+    const argument_t &r_jk = arguments[3]; // 12
+    const argument_t &r_jl = arguments[4]; // 13
+    const argument_t &r_kl = arguments[5]; // 23
+    result +=
+        2. * params[spring_constant] *
+        ( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance)
+            - params[equilibrium_distance]);
+  }
+  return derivative_components_t(1, result);
 }
 
 FourBodyPotential_Torsion::results_t
 FourBodyPotential_Torsion::parameter_derivative(
-    const arguments_t &arguments,
+    const list_of_arguments_t &listarguments,
     const size_t index
     ) const
 {
-  ASSERT( arguments.size() == getSpecificArgumentCount(),
-      "FourBodyPotential_Torsion::parameter_derivative() - requires exactly three arguments.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments.");
-  const argument_t &r_ij = arguments[0]; // 01
-  const argument_t &r_ik = arguments[1]; // 02
-  const argument_t &r_il = arguments[2]; // 03
-  const argument_t &r_jk = arguments[3]; // 12
-  const argument_t &r_jl = arguments[4]; // 13
-  const argument_t &r_kl = arguments[5]; // 23
-  switch (index) {
-    case spring_constant:
-    {
-      const result_t result =
-                 Helpers::pow( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance], 2 );
-      return std::vector<result_t>(1, result);
-      break;
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    ASSERT( arguments.size() == getSpecificArgumentCount(),
+        "FourBodyPotential_Torsion::parameter_derivative() - requires exactly three arguments.");
+    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+        arguments, getParticleTypes()),
+        "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments.");
+    const argument_t &r_ij = arguments[0]; // 01
+    const argument_t &r_ik = arguments[1]; // 02
+    const argument_t &r_il = arguments[2]; // 03
+    const argument_t &r_jk = arguments[3]; // 12
+    const argument_t &r_jl = arguments[4]; // 13
+    const argument_t &r_kl = arguments[5]; // 23
+    switch (index) {
+      case spring_constant:
+      {
+        result +=
+            Helpers::pow( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance], 2 );
+        break;
+      }
+      case equilibrium_distance:
+      {
+        result +=
+            -2. * params[spring_constant]
+                   * ( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance]);
+        break;
+      }
+      default:
+        ASSERT(0, "FourBodyPotential_Torsion::parameter_derivative() - derivative to unknown parameter desired.");
+        break;
     }
-    case equilibrium_distance:
-    {
-      const result_t result =
-          -2. * params[spring_constant]
-                 * ( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance]);
-      return std::vector<result_t>(1, result);
-      break;
-    }
-    default:
-      ASSERT(0, "FourBodyPotential_Torsion::parameter_derivative() - derivative to unknown parameter desired.");
-      break;
   }
-  return std::vector<result_t>(1);
+  return results_t(1, result);
 }
 
Index: src/Potentials/Specifics/FourBodyPotential_Torsion.hpp
===================================================================
--- src/Potentials/Specifics/FourBodyPotential_Torsion.hpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/FourBodyPotential_Torsion.hpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -37,4 +37,5 @@
   friend class PotentialFactory;
   // some repeated typedefs to avoid ambiguities
+  typedef FunctionModel::list_of_arguments_t list_of_arguments_t;
   typedef FunctionModel::arguments_t arguments_t;
   typedef FunctionModel::result_t result_t;
@@ -92,25 +93,25 @@
   /** Evaluates the harmonic potential function for the given arguments.
    *
-   * @param arguments single distance
+   * @param listarguments list of six distances
    * @return value of the potential function
    */
-  results_t operator()(const arguments_t &arguments) const;
+  results_t operator()(const list_of_arguments_t &listarguments) const;
 
   /** Evaluates the derivative of the potential function.
    *
-   * @param arguments single distance
+   * @param listarguments list of six distances
    * @return vector with derivative with respect to the input degrees of freedom
    */
-  derivative_components_t derivative(const arguments_t &arguments) const;
+  derivative_components_t derivative(const list_of_arguments_t &listarguments) const;
 
   /** Evaluates the derivative of the function with the given \a arguments
    * with respect to a specific parameter indicated by \a index.
    *
-   * \param arguments set of arguments as input variables to the function
+   * \param listarguments list of six distances
    * \param index derivative of which parameter
    * \return result vector containing the derivative with respect to the given
    *         input
    */
-  results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;
+  results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const;
 
   /** Returns the functor that converts argument_s into the
Index: src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp
===================================================================
--- src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -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);
Index: src/Potentials/Specifics/ManyBodyPotential_Tersoff.hpp
===================================================================
--- src/Potentials/Specifics/ManyBodyPotential_Tersoff.hpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/ManyBodyPotential_Tersoff.hpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -38,4 +38,5 @@
   friend class PotentialFactory;
   // some repeated typedefs to avoid ambiguities
+  typedef FunctionModel::list_of_arguments_t list_of_arguments_t;
   typedef FunctionModel::arguments_t arguments_t;
   typedef FunctionModel::result_t result_t;
@@ -106,26 +107,26 @@
   /** Evaluates the Tersoff potential for the given arguments.
    *
-   * @param arguments single distance
+   * @param listarguments list of distances
    * @return value of the potential function
    */
-  results_t operator()(const arguments_t &arguments) const;
+  results_t operator()(const list_of_arguments_t &listarguments) const;
 
   /** Evaluates the derivative of the Tersoff potential with respect to the
    * input variables.
    *
-   * @param arguments single distance
+   * @param listarguments list of distances
    * @return vector with components of the derivative
    */
-  derivative_components_t derivative(const arguments_t &arguments) const;
+  derivative_components_t derivative(const list_of_arguments_t &listarguments) const;
 
   /** Evaluates the derivative of the function with the given \a arguments
    * with respect to a specific parameter indicated by \a index.
    *
-   * \param arguments set of arguments as input variables to the function
+   * \param listarguments list of distances
    * \param index derivative of which parameter
    * \return result vector containing the derivative with respect to the given
    *         input
    */
-  results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;
+  results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const;
 
   /** Returns the functor that converts argument_s into the
Index: src/Potentials/Specifics/PairPotential_Harmonic.cpp
===================================================================
--- src/Potentials/Specifics/PairPotential_Harmonic.cpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/PairPotential_Harmonic.cpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -117,71 +117,82 @@
 PairPotential_Harmonic::results_t
 PairPotential_Harmonic::operator()(
-    const arguments_t &arguments
+    const list_of_arguments_t &listarguments
     ) const
 {
-  ASSERT( arguments.size() == 1,
-      "PairPotential_Harmonic::operator() - requires exactly one argument.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "PairPotential_Harmonic::operator() - types don't match with ones in arguments.");
-  const argument_t &r_ij = arguments[0];
-  const result_t result =
-      params[spring_constant]
-             * Helpers::pow( r_ij.distance - params[equilibrium_distance], 2 );
-  return std::vector<result_t>(1, result);
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    ASSERT( arguments.size() == 1,
+        "PairPotential_Harmonic::operator() - requires exactly one argument.");
+    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+        arguments, getParticleTypes()),
+        "PairPotential_Harmonic::operator() - types don't match with ones in arguments.");
+    const argument_t &r_ij = arguments[0];
+    result +=
+        params[spring_constant]
+               * Helpers::pow( r_ij.distance - params[equilibrium_distance], 2 );
+  }
+  return results_t(1, result);
 }
 
 PairPotential_Harmonic::derivative_components_t
 PairPotential_Harmonic::derivative(
-    const arguments_t &arguments
+    const list_of_arguments_t &listarguments
     ) const
 {
-  ASSERT( arguments.size() == 1,
-      "PairPotential_Harmonic::operator() - requires exactly one argument.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "PairPotential_Harmonic::operator() - types don't match with ones in arguments.");
-  derivative_components_t result;
-  const argument_t &r_ij = arguments[0];
-  result.push_back( 2. * params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]) );
-  ASSERT( result.size() == 1,
-      "PairPotential_Harmonic::operator() - we did not create exactly one component.");
-  return result;
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+      ASSERT( arguments.size() == 1,
+          "PairPotential_Harmonic::operator() - requires exactly one argument.");
+      ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+          arguments, getParticleTypes()),
+          "PairPotential_Harmonic::operator() - types don't match with ones in arguments.");
+      const argument_t &r_ij = arguments[0];
+      result +=
+          2. * params[spring_constant] *
+          ( r_ij.distance - params[equilibrium_distance]);
+  }
+  return derivative_components_t(1, result);
 }
 
 PairPotential_Harmonic::results_t
 PairPotential_Harmonic::parameter_derivative(
-    const arguments_t &arguments,
+    const list_of_arguments_t &listarguments,
     const size_t index
     ) const
 {
-  ASSERT( arguments.size() == 1,
-      "PairPotential_Harmonic::parameter_derivative() - requires exactly one argument.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "PairPotential_Harmonic::operator() - types don't match with ones in arguments.");
-  const argument_t &r_ij = arguments[0];
-  switch (index) {
-    case spring_constant:
-    {
-      const result_t result =
-                 Helpers::pow( r_ij.distance - params[equilibrium_distance], 2 );
-      return std::vector<result_t>(1, result);
-      break;
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    ASSERT( arguments.size() == 1,
+        "PairPotential_Harmonic::parameter_derivative() - requires exactly one argument.");
+    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+        arguments, getParticleTypes()),
+        "PairPotential_Harmonic::operator() - types don't match with ones in arguments.");
+    const argument_t &r_ij = arguments[0];
+    switch (index) {
+      case spring_constant:
+      {
+        result +=
+            Helpers::pow( r_ij.distance - params[equilibrium_distance], 2 );
+        break;
+      }
+      case equilibrium_distance:
+      {
+        result +=
+            -2. * params[spring_constant]
+                   * ( r_ij.distance - params[equilibrium_distance]);
+        break;
+      }
+      default:
+        ASSERT(0, "PairPotential_Harmonic::parameter_derivative() - derivative to unknown parameter desired.");
+        break;
     }
-    case equilibrium_distance:
-    {
-      const result_t result =
-          -2. * params[spring_constant]
-                 * ( r_ij.distance - params[equilibrium_distance]);
-      return std::vector<result_t>(1, result);
-      break;
-    }
-    default:
-      ASSERT(0, "PairPotential_Harmonic::parameter_derivative() - derivative to unknown parameter desired.");
-      break;
   }
-
-  return PairPotential_Harmonic::results_t(1, 0.);
+  return results_t(1, result);
 }
 
Index: src/Potentials/Specifics/PairPotential_Harmonic.hpp
===================================================================
--- src/Potentials/Specifics/PairPotential_Harmonic.hpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/PairPotential_Harmonic.hpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -35,4 +35,5 @@
   friend class PotentialFactory;
   // some repeated typedefs to avoid ambiguities
+  typedef FunctionModel::list_of_arguments_t list_of_arguments_t;
   typedef FunctionModel::arguments_t arguments_t;
   typedef FunctionModel::result_t result_t;
@@ -90,25 +91,25 @@
   /** Evaluates the harmonic potential function for the given arguments.
    *
-   * @param arguments single distance
+   * @param listarguments list of single distances
    * @return value of the potential function
    */
-  results_t operator()(const arguments_t &arguments) const;
+  results_t operator()(const list_of_arguments_t &listarguments) const;
 
   /** Evaluates the derivative of the potential function.
    *
-   * @param arguments single distance
+   * @param listarguments list of single distances
    * @return vector with derivative with respect to the input degrees of freedom
    */
-  derivative_components_t derivative(const arguments_t &arguments) const;
+  derivative_components_t derivative(const list_of_arguments_t &listarguments) const;
 
   /** Evaluates the derivative of the function with the given \a arguments
    * with respect to a specific parameter indicated by \a index.
    *
-   * \param arguments set of arguments as input variables to the function
+   * \param listarguments list of single distances
    * \param index derivative of which parameter
    * \return result vector containing the derivative with respect to the given
    *         input
    */
-  results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;
+  results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const;
 
   /** Returns the functor that converts argument_s into the
Index: src/Potentials/Specifics/PairPotential_LennardJones.cpp
===================================================================
--- src/Potentials/Specifics/PairPotential_LennardJones.cpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/PairPotential_LennardJones.cpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -123,74 +123,86 @@
 PairPotential_LennardJones::results_t
 PairPotential_LennardJones::operator()(
-    const arguments_t &arguments
+    const list_of_arguments_t &listarguments
     ) const
 {
-  ASSERT( arguments.size() == 1,
-      "PairPotential_LennardJones::operator() - requires one argument.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "PairPotential_LennardJones::operator() - types don't match with ones in arguments.");
-  const double &r = arguments[0].distance;
-  const double temp = Helpers::pow(params[sigma]/r, 6);
-  const result_t result = 4.*params[epsilon] * (temp*temp - temp);
-  return std::vector<result_t>(1, result);
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    ASSERT( arguments.size() == 1,
+        "PairPotential_LennardJones::operator() - requires one argument.");
+    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+        arguments, getParticleTypes()),
+        "PairPotential_LennardJones::operator() - types don't match with ones in arguments.");
+    const double &r = arguments[0].distance;
+    const double temp = Helpers::pow(params[sigma]/r, 6);
+    result += 4.*params[epsilon] * (temp*temp - temp);
+  }
+  return results_t(1, result);
 }
 
 PairPotential_LennardJones::derivative_components_t
 PairPotential_LennardJones::derivative(
-    const arguments_t &arguments
+    const list_of_arguments_t &listarguments
     ) const
 {
-  ASSERT( arguments.size() == 1,
-      "PairPotential_LennardJones::operator() - requires no argument.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "PairPotential_LennardJones::operator() - types don't match with ones in arguments.");
-  const double &r = arguments[0].distance;
   const double sigma6 = Helpers::pow(params[sigma], 6);
-  const result_t result =
-      4.*params[epsilon] * (
-          sigma6*sigma6*(-12.) / Helpers::pow(r,13)
-          - sigma6*(-6.) /Helpers::pow(r,7)
-      );
-  derivative_components_t results(1, result);
-  return results;
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    ASSERT( arguments.size() == 1,
+        "PairPotential_LennardJones::operator() - requires no argument.");
+    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+        arguments, getParticleTypes()),
+        "PairPotential_LennardJones::operator() - types don't match with ones in arguments.");
+    const double &r = arguments[0].distance;
+    result +=
+        4.*params[epsilon] * (
+            sigma6*sigma6*(-12.) / Helpers::pow(r,13)
+            - sigma6*(-6.) /Helpers::pow(r,7)
+        );
+  }
+  return derivative_components_t(1, result);
 }
 
 PairPotential_LennardJones::results_t
 PairPotential_LennardJones::parameter_derivative(
-    const arguments_t &arguments,
+    const list_of_arguments_t &listarguments,
     const size_t index
     ) const
 {
-  ASSERT( arguments.size() == 1,
-      "PairPotential_LennardJones::parameter_derivative() - requires no argument.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "PairPotential_LennardJones::operator() - types don't match with ones in arguments.");
-  const double &r = arguments[0].distance;
-  switch (index) {
-    case epsilon:
-    {
-      const double temp = Helpers::pow(params[sigma]/r, 6);
-      const result_t result = 4. * (temp*temp - temp);
-      return std::vector<result_t>(1, result);
-      break;
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    ASSERT( arguments.size() == 1,
+        "PairPotential_LennardJones::parameter_derivative() - requires no argument.");
+    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+        arguments, getParticleTypes()),
+        "PairPotential_LennardJones::operator() - types don't match with ones in arguments.");
+    const double &r = arguments[0].distance;
+    switch (index) {
+      case epsilon:
+      {
+        const double temp = Helpers::pow(params[sigma]/r, 6);
+        result += 4. * (temp*temp - temp);
+        break;
+      }
+      case sigma:
+      {
+        const double r6 = Helpers::pow(r, 6);
+        result +=
+            4.*params[epsilon] * (
+                12. * Helpers::pow(params[sigma],11)/(r6*r6)
+                - 6. * Helpers::pow(params[sigma],5)/r6
+            );
+        break;
+      }
+      default:
+        break;
     }
-    case sigma:
-    {
-      const double r6 = Helpers::pow(r, 6);
-      const result_t result =
-          4.*params[epsilon] * (
-              12. * Helpers::pow(params[sigma],11)/(r6*r6)
-              - 6. * Helpers::pow(params[sigma],5)/r6
-          );
-      return std::vector<result_t>(1, result);
-      break;
-    }
-    default:
-      break;
   }
-  return std::vector<result_t>(1, 0.);
+  return results_t(1, result);
 }
 
Index: src/Potentials/Specifics/PairPotential_LennardJones.hpp
===================================================================
--- src/Potentials/Specifics/PairPotential_LennardJones.hpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/PairPotential_LennardJones.hpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -36,4 +36,5 @@
   friend class PotentialFactory;
   // some repeated typedefs to avoid ambiguities
+  typedef FunctionModel::list_of_arguments_t list_of_arguments_t;
   typedef FunctionModel::arguments_t arguments_t;
   typedef FunctionModel::result_t result_t;
@@ -92,25 +93,25 @@
   /** Evaluates the harmonic potential function for the given arguments.
    *
-   * @param arguments single distance
+   * @param listarguments list of single distances
    * @return value of the potential function
    */
-  results_t operator()(const arguments_t &arguments) const;
+  results_t operator()(const list_of_arguments_t &listarguments) const;
 
   /** Evaluates the derivative of the potential function.
    *
-   * @param arguments single distance
+   * @param listarguments list of single distances
    * @return vector with derivative with respect to the input degrees of freedom
    */
-  derivative_components_t derivative(const arguments_t &arguments) const;
+  derivative_components_t derivative(const list_of_arguments_t &listarguments) const;
 
   /** Evaluates the derivative of the function with the given \a arguments
    * with respect to a specific parameter indicated by \a index.
    *
-   * \param arguments set of arguments as input variables to the function
+   * \param listarguments list of single distances
    * \param index derivative of which parameter
    * \return result vector containing the derivative with respect to the given
    *         input
    */
-  results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;
+  results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const;
 
   /** Returns the functor that converts argument_s into the
Index: src/Potentials/Specifics/PairPotential_Morse.cpp
===================================================================
--- src/Potentials/Specifics/PairPotential_Morse.cpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/PairPotential_Morse.cpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -124,17 +124,22 @@
 PairPotential_Morse::results_t
 PairPotential_Morse::operator()(
-    const arguments_t &arguments
+    const list_of_arguments_t &listarguments
     ) const
 {
-  ASSERT( arguments.size() == 1,
-      "PairPotential_Morse::operator() - requires exactly one argument.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "PairPotential_Morse::operator() - types don't match with ones in arguments.");
-  const argument_t &r_ij = arguments[0];
-  // Maple: f(r,D,k,R,c) := D * (1 - exp(-k*(r-R)))^(2)+c
-  const result_t result =
-      params[dissociation_energy] * Helpers::pow( 1.
-          - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])),2);
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    ASSERT( arguments.size() == 1,
+        "PairPotential_Morse::operator() - requires exactly one argument.");
+    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+        arguments, getParticleTypes()),
+        "PairPotential_Morse::operator() - types don't match with ones in arguments.");
+    const argument_t &r_ij = arguments[0];
+    // Maple: f(r,D,k,R,c) := D * (1 - exp(-k*(r-R)))^(2)+c
+    result +=
+        params[dissociation_energy] * Helpers::pow( 1.
+            - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])),2);
+  }
   return std::vector<result_t>(1, result);
 }
@@ -142,75 +147,78 @@
 PairPotential_Morse::derivative_components_t
 PairPotential_Morse::derivative(
-    const arguments_t &arguments
+    const list_of_arguments_t &listarguments
     ) const
 {
-  ASSERT( arguments.size() == 1,
-      "PairPotential_Morse::operator() - requires exactly one argument.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "PairPotential_Morse::operator() - types don't match with ones in arguments.");
-  derivative_components_t result;
-  const argument_t &r_ij = arguments[0];
-  // Maple result: 2*D*(1-exp(-k*(r-R)))*k*exp(-k*(r-R))
-  result.push_back(
-      2. * params[dissociation_energy]
-      * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])))
-      * (- params[spring_constant]) * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))
-  );
-  ASSERT( result.size() == 1,
-      "PairPotential_Morse::operator() - we did not create exactly one component.");
-  return result;
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    ASSERT( arguments.size() == 1,
+        "PairPotential_Morse::operator() - requires exactly one argument.");
+    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+        arguments, getParticleTypes()),
+        "PairPotential_Morse::operator() - types don't match with ones in arguments.");
+    const argument_t &r_ij = arguments[0];
+    // Maple result: 2*D*(1-exp(-k*(r-R)))*k*exp(-k*(r-R))
+    result +=
+        2. * params[dissociation_energy]
+        * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])))
+        * (- params[spring_constant]) * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]));
+  }
+  return derivative_components_t(1, result);
 }
 
 PairPotential_Morse::results_t
 PairPotential_Morse::parameter_derivative(
-    const arguments_t &arguments,
+    const list_of_arguments_t &listarguments,
     const size_t index
     ) const
 {
-  ASSERT( arguments.size() == 1,
-      "PairPotential_Morse::parameter_derivative() - requires exactly one argument.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "PairPotential_Morse::operator() - types don't match with ones in arguments.");
-  const argument_t &r_ij = arguments[0];
-  switch (index) {
-    case spring_constant:
-    {
-      // Maple result: -2*D*(1-exp(-k*(r-R)))*(-r+R)*exp(-k*(r-R))
-      const result_t result =
-          - 2. * params[dissociation_energy]
-          * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])))
-          * (- r_ij.distance + params[equilibrium_distance])
-          * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))
-          ;
-      return std::vector<result_t>(1, result);
-      break;
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    ASSERT( arguments.size() == 1,
+        "PairPotential_Morse::parameter_derivative() - requires exactly one argument.");
+    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+        arguments, getParticleTypes()),
+        "PairPotential_Morse::operator() - types don't match with ones in arguments.");
+    const argument_t &r_ij = arguments[0];
+    switch (index) {
+      case spring_constant:
+      {
+        // Maple result: -2*D*(1-exp(-k*(r-R)))*(-r+R)*exp(-k*(r-R))
+        result +=
+            - 2. * params[dissociation_energy]
+            * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])))
+            * (- r_ij.distance + params[equilibrium_distance])
+            * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))
+            ;
+        break;
+      }
+      case equilibrium_distance:
+      {
+        // Maple result: -2*D*(1-exp(-k*(r-R)))*k*exp(-k*(r-R))
+        result +=
+            - 2. * params[dissociation_energy]
+            * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])))
+            * params[spring_constant] * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))
+            ;
+        break;
+      }
+      case dissociation_energy:
+      {
+        // Maple result: (1-exp(-k*(r-R)))^2
+        result +=
+            Helpers::pow(1.
+                - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])),2);
+        break;
+      }
+      default:
+        ASSERT(0, "PairPotential_Morse::parameter_derivative() - derivative to unknown parameter desired.");
+        break;
     }
-    case equilibrium_distance:
-    {
-      // Maple result: -2*D*(1-exp(-k*(r-R)))*k*exp(-k*(r-R))
-      const result_t result =
-          - 2. * params[dissociation_energy]
-          * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])))
-          * params[spring_constant] * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))
-          ;
-      return std::vector<result_t>(1, result);
-      break;
-    }
-    case dissociation_energy:
-    {
-      // Maple result: (1-exp(-k*(r-R)))^2
-      const result_t result =
-          Helpers::pow(1.
-              - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])),2);
-      return std::vector<result_t>(1, result);
-      break;
-    }
-    default:
-      ASSERT(0, "PairPotential_Morse::parameter_derivative() - derivative to unknown parameter desired.");
-      break;
   }
-  return std::vector<result_t>(1, 0.);
+  return results_t(1, result);
 }
 
Index: src/Potentials/Specifics/PairPotential_Morse.hpp
===================================================================
--- src/Potentials/Specifics/PairPotential_Morse.hpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/PairPotential_Morse.hpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -35,4 +35,5 @@
   friend class PotentialFactory;
   // some repeated typedefs to avoid ambiguities
+  typedef FunctionModel::list_of_arguments_t list_of_arguments_t;
   typedef FunctionModel::arguments_t arguments_t;
   typedef FunctionModel::result_t result_t;
@@ -91,25 +92,25 @@
   /** Evaluates the harmonic potential function for the given arguments.
    *
-   * @param arguments single distance
+   * @param listarguments list of single distances
    * @return value of the potential function
    */
-  results_t operator()(const arguments_t &arguments) const;
+  results_t operator()(const list_of_arguments_t &listarguments) const;
 
   /** Evaluates the derivative of the potential function.
    *
-   * @param arguments single distance
+   * @param listarguments list of single distances
    * @return vector with derivative with respect to the input degrees of freedom
    */
-  derivative_components_t derivative(const arguments_t &arguments) const;
+  derivative_components_t derivative(const list_of_arguments_t &listarguments) const;
 
   /** Evaluates the derivative of the function with the given \a arguments
    * with respect to a specific parameter indicated by \a index.
    *
-   * \param arguments set of arguments as input variables to the function
+   * \param listarguments list of single distances
    * \param index derivative of which parameter
    * \return result vector containing the derivative with respect to the given
    *         input
    */
-  results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;
+  results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const;
 
   /** Returns the functor that converts argument_s into the
Index: src/Potentials/Specifics/ThreeBodyPotential_Angle.cpp
===================================================================
--- src/Potentials/Specifics/ThreeBodyPotential_Angle.cpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/ThreeBodyPotential_Angle.cpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -136,76 +136,90 @@
 ThreeBodyPotential_Angle::results_t
 ThreeBodyPotential_Angle::operator()(
-    const arguments_t &arguments
+    const list_of_arguments_t &listarguments
     ) const
 {
-  ASSERT( arguments.size() == 3,
-      "ThreeBodyPotential_Angle::operator() - requires exactly three arguments.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments.");
-  const argument_t &r_ij = arguments[0]; // 01
-  const argument_t &r_jk = arguments[2]; // 12
-  const argument_t &r_ik = arguments[1]; // 02
-  const result_t result =
-      params[spring_constant]
-             * Helpers::pow( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance], 2 );
-  return std::vector<result_t>(1, result);
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    ASSERT( arguments.size() == 3,
+        "ThreeBodyPotential_Angle::operator() - requires exactly three arguments.");
+    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+        arguments, getParticleTypes()),
+        "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments.");
+    const argument_t &r_ij = arguments[0]; // 01
+    const argument_t &r_jk = arguments[2]; // 12
+    const argument_t &r_ik = arguments[1]; // 02
+    result +=
+        params[spring_constant]
+               * Helpers::pow( function_theta(r_ij.distance, r_jk.distance, r_ik.distance)
+                   - params[equilibrium_distance], 2 );
+  }
+  return results_t(1, result);
 }
 
 ThreeBodyPotential_Angle::derivative_components_t
 ThreeBodyPotential_Angle::derivative(
-    const arguments_t &arguments
+    const list_of_arguments_t &listarguments
     ) const
 {
-  ASSERT( arguments.size() == 3,
-      "ThreeBodyPotential_Angle::operator() - requires exactly three arguments.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments.");
-  derivative_components_t result;
-  const argument_t &r_ij = arguments[0]; //01
-  const argument_t &r_jk = arguments[2]; //12
-  const argument_t &r_ik = arguments[1]; //02
-  result.push_back( 2. * params[spring_constant] * ( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance]) );
-  ASSERT( result.size() == 1,
-      "ThreeBodyPotential_Angle::operator() - we did not create exactly one component.");
-  return result;
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    ASSERT( arguments.size() == 3,
+        "ThreeBodyPotential_Angle::operator() - requires exactly three arguments.");
+    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+        arguments, getParticleTypes()),
+        "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments.");
+    const argument_t &r_ij = arguments[0]; //01
+    const argument_t &r_jk = arguments[2]; //12
+    const argument_t &r_ik = arguments[1]; //02
+    result +=
+        2. * params[spring_constant] *
+        ( function_theta(r_ij.distance, r_jk.distance, r_ik.distance)
+            - params[equilibrium_distance]);
+  }
+  return derivative_components_t(1, result);
 }
 
 ThreeBodyPotential_Angle::results_t
 ThreeBodyPotential_Angle::parameter_derivative(
-    const arguments_t &arguments,
+    const list_of_arguments_t &listarguments,
     const size_t index
     ) const
 {
-  ASSERT( arguments.size() == 3,
-      "ThreeBodyPotential_Angle::parameter_derivative() - requires exactly three arguments.");
-  ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
-      arguments, getParticleTypes()),
-      "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments.");
-  const argument_t &r_ij = arguments[0]; //01
-  const argument_t &r_jk = arguments[2]; //12
-  const argument_t &r_ik = arguments[1]; //02
-  switch (index) {
-    case spring_constant:
-    {
-      const result_t result =
-                 Helpers::pow( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance], 2 );
-      return std::vector<result_t>(1, result);
-      break;
+  result_t result = 0.;
+  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
+      iter != listarguments.end(); ++iter) {
+    const arguments_t &arguments = *iter;
+    ASSERT( arguments.size() == 3,
+        "ThreeBodyPotential_Angle::parameter_derivative() - requires exactly three arguments.");
+    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
+        arguments, getParticleTypes()),
+        "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments.");
+    const argument_t &r_ij = arguments[0]; //01
+    const argument_t &r_jk = arguments[2]; //12
+    const argument_t &r_ik = arguments[1]; //02
+    switch (index) {
+      case spring_constant:
+      {
+        result +=
+            Helpers::pow( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance], 2 );
+        break;
+      }
+      case equilibrium_distance:
+      {
+        result +=
+            -2. * params[spring_constant]
+                   * ( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance]);
+        break;
+      }
+      default:
+        ASSERT(0, "ThreeBodyPotential_Angle::parameter_derivative() - derivative to unknown parameter desired.");
+        break;
     }
-    case equilibrium_distance:
-    {
-      const result_t result =
-          -2. * params[spring_constant]
-                 * ( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance]);
-      return std::vector<result_t>(1, result);
-      break;
-    }
-    default:
-      ASSERT(0, "ThreeBodyPotential_Angle::parameter_derivative() - derivative to unknown parameter desired.");
-      break;
   }
-  return std::vector<result_t>(1);
+  return results_t(1, result);
 }
 
Index: src/Potentials/Specifics/ThreeBodyPotential_Angle.hpp
===================================================================
--- src/Potentials/Specifics/ThreeBodyPotential_Angle.hpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/ThreeBodyPotential_Angle.hpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -35,4 +35,5 @@
   friend class PotentialFactory;
   // some repeated typedefs to avoid ambiguities
+  typedef FunctionModel::list_of_arguments_t list_of_arguments_t;
   typedef FunctionModel::arguments_t arguments_t;
   typedef FunctionModel::result_t result_t;
@@ -90,25 +91,25 @@
   /** Evaluates the harmonic potential function for the given arguments.
    *
-   * @param arguments single distance
+   * @param listarguments list of three distances
    * @return value of the potential function
    */
-  results_t operator()(const arguments_t &arguments) const;
+  results_t operator()(const list_of_arguments_t &listarguments) const;
 
   /** Evaluates the derivative of the potential function.
    *
-   * @param arguments single distance
+   * @param listarguments list of three distances
    * @return vector with derivative with respect to the input degrees of freedom
    */
-  derivative_components_t derivative(const arguments_t &arguments) const;
+  derivative_components_t derivative(const list_of_arguments_t &listarguments) const;
 
   /** Evaluates the derivative of the function with the given \a arguments
    * with respect to a specific parameter indicated by \a index.
    *
-   * \param arguments set of arguments as input variables to the function
+   * \param listarguments list of three distances
    * \param index derivative of which parameter
    * \return result vector containing the derivative with respect to the given
    *         input
    */
-  results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;
+  results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const;
 
   /** Returns the functor that converts argument_s into the
Index: src/Potentials/Specifics/unittests/ConstantPotentialUnitTest.cpp
===================================================================
--- src/Potentials/Specifics/unittests/ConstantPotentialUnitTest.cpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/unittests/ConstantPotentialUnitTest.cpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -82,5 +82,5 @@
       Helpers::isEqual(
           offset,
-          constant( FunctionModel::arguments_t() )[0]
+          constant( FunctionModel::list_of_arguments_t() )[0]
       )
   );
@@ -97,5 +97,5 @@
           0.,
           constant.derivative(
-              FunctionModel::arguments_t()
+              FunctionModel::list_of_arguments_t()
           )[0]
       )
@@ -114,5 +114,5 @@
           1.,
           constant.parameter_derivative(
-              FunctionModel::arguments_t(),
+              FunctionModel::list_of_arguments_t(),
               0
           )[0]
Index: src/Potentials/Specifics/unittests/FourBodyPotential_ImproperUnitTest.cpp
===================================================================
--- src/Potentials/Specifics/unittests/FourBodyPotential_ImproperUnitTest.cpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/unittests/FourBodyPotential_ImproperUnitTest.cpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -121,5 +121,5 @@
   for (size_t index = 0; index < input.size(); ++index) {
     const FourBodyPotential_Improper::results_t result =
-        angle( input[index] );
+        angle( FunctionModel::list_of_arguments_t(1, input[index]) );
     CPPUNIT_ASSERT(
         Helpers::isEqual(
@@ -144,5 +144,5 @@
           0.,
           angle.derivative(
-              input[5]
+              FunctionModel::list_of_arguments_t(1, input[5])
           )[0],
           10.
@@ -165,5 +165,5 @@
           0.,
           angle.parameter_derivative(
-              input[5],
+              FunctionModel::list_of_arguments_t(1, input[5]),
               0
           )[0],
@@ -175,5 +175,5 @@
           0.,
           angle.parameter_derivative(
-              input[5],
+              FunctionModel::list_of_arguments_t(1, input[5]),
               1
           )[0],
Index: src/Potentials/Specifics/unittests/FourBodyPotential_TorsionUnitTest.cpp
===================================================================
--- src/Potentials/Specifics/unittests/FourBodyPotential_TorsionUnitTest.cpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/unittests/FourBodyPotential_TorsionUnitTest.cpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -121,5 +121,5 @@
   for (size_t index = 0; index < input.size(); ++index) {
     const FourBodyPotential_Torsion::results_t result =
-        angle( input[index] );
+        angle( FunctionModel::list_of_arguments_t(1, input[index]) );
     CPPUNIT_ASSERT(
         Helpers::isEqual(
@@ -144,5 +144,5 @@
           0.,
           angle.derivative(
-              input[5]
+              FunctionModel::list_of_arguments_t(1, input[5])
           )[0],
           10.
@@ -165,5 +165,5 @@
           0.,
           angle.parameter_derivative(
-              input[5],
+              FunctionModel::list_of_arguments_t(1, input[5]),
               0
           )[0],
@@ -175,5 +175,5 @@
           0.,
           angle.parameter_derivative(
-              input[5],
+              FunctionModel::list_of_arguments_t(1, input[5]),
               1
           )[0],
Index: src/Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.cpp
===================================================================
--- src/Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.cpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.cpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -299,5 +299,6 @@
         arg.globalid = index; // this is needed for the triplefunction to the configuration
         FunctionModel::arguments_t args(1,arg);
-        const ManyBodyPotential_Tersoff::results_t res = tersoff(args);
+        const ManyBodyPotential_Tersoff::results_t res =
+            tersoff( FunctionModel::list_of_arguments_t(1, args) );
         temp += res[0];
       }
@@ -336,5 +337,5 @@
 //          0.,
 //          tersoff.derivative(
-//              FunctionModel::arguments_t(1,argument_t(1.))
+//              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,argument_t(1.)))
 //          )[0]
 //      )
@@ -361,5 +362,5 @@
 //          0.,
 //          tersoff.parameter_derivative(
-//              FunctionModel::arguments_t(1,argument_t(1.)),
+//              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,argument_t(1.))),
 //              0
 //          )[0]
@@ -370,5 +371,5 @@
 //          0.,
 //          tersoff.parameter_derivative(
-//              FunctionModel::arguments_t(1,argument_t(1.)),
+//              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,argument_t(1.))),
 //              1
 //          )[0]
@@ -379,5 +380,5 @@
 //          1.,
 //          tersoff.parameter_derivative(
-//              FunctionModel::arguments_t(1,argument_t(1.)),
+//              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,argument_t(1.))),
 //              2
 //          )[0]
Index: src/Potentials/Specifics/unittests/PairPotential_HarmonicUnitTest.cpp
===================================================================
--- src/Potentials/Specifics/unittests/PairPotential_HarmonicUnitTest.cpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/unittests/PairPotential_HarmonicUnitTest.cpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -102,5 +102,5 @@
         Helpers::isEqual(
             output[index],
-            harmonic( FunctionModel::arguments_t(1,arg) )[0]
+            harmonic( FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)) )[0]
         )
     );
@@ -122,5 +122,5 @@
           0.,
           harmonic.derivative(
-              FunctionModel::arguments_t(1,arg)
+              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg))
           )[0]
       )
@@ -143,5 +143,5 @@
           0.,
           harmonic.parameter_derivative(
-              FunctionModel::arguments_t(1,arg),
+              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)),
               0
           )[0]
@@ -152,5 +152,5 @@
           0.,
           harmonic.parameter_derivative(
-              FunctionModel::arguments_t(1,arg),
+              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)),
               1
           )[0]
Index: src/Potentials/Specifics/unittests/PairPotential_LennardJonesUnitTest.cpp
===================================================================
--- src/Potentials/Specifics/unittests/PairPotential_LennardJonesUnitTest.cpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/unittests/PairPotential_LennardJonesUnitTest.cpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -118,5 +118,5 @@
         Helpers::isEqual(
             output[index],
-            lj( FunctionModel::arguments_t(1,arg) )[0],
+            lj( FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)) )[0],
             1.e-4/std::numeric_limits<double>::epsilon() // only compare four digits
         )
@@ -140,5 +140,5 @@
           0.,
           lj.derivative(
-              FunctionModel::arguments_t(1,arg)
+              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg))
           )[0],
           1.e+6
@@ -162,5 +162,5 @@
           -1.,
           lj.parameter_derivative(
-              FunctionModel::arguments_t(1,arg),
+              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)),
               0
           )[0],
@@ -172,5 +172,5 @@
           0.,
           lj.parameter_derivative(
-              FunctionModel::arguments_t(1,arg),
+              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)),
               1
           )[0],
Index: src/Potentials/Specifics/unittests/PairPotential_MorseUnitTest.cpp
===================================================================
--- src/Potentials/Specifics/unittests/PairPotential_MorseUnitTest.cpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/unittests/PairPotential_MorseUnitTest.cpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -120,5 +120,5 @@
         Helpers::isEqual(
             output[index],
-            Morse( FunctionModel::arguments_t(1,arg) )[0],
+            Morse( FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)) )[0],
             1.e-4/std::numeric_limits<double>::epsilon() // only compare four digits
         )
@@ -141,5 +141,5 @@
           0.,
           Morse.derivative(
-              FunctionModel::arguments_t(1,arg)
+              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg))
           )[0],
           1.e+6
@@ -163,5 +163,5 @@
           0.,
           Morse.parameter_derivative(
-              FunctionModel::arguments_t(1,arg),
+              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)),
               0
           )[0],
@@ -173,5 +173,5 @@
           0.,
           Morse.parameter_derivative(
-              FunctionModel::arguments_t(1,arg),
+              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)),
               1
           )[0],
@@ -183,5 +183,5 @@
           0.,
           Morse.parameter_derivative(
-              FunctionModel::arguments_t(1,arg),
+              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)),
               2
           )[0],
Index: src/Potentials/Specifics/unittests/ThreeBodyPotential_AngleUnitTest.cpp
===================================================================
--- src/Potentials/Specifics/unittests/ThreeBodyPotential_AngleUnitTest.cpp	(revision 16227a32089e461901d97290d3e2fb626919cf26)
+++ src/Potentials/Specifics/unittests/ThreeBodyPotential_AngleUnitTest.cpp	(revision e1fe7e9e5ab6e5acab3403d31d4139a02a1c1754)
@@ -112,5 +112,5 @@
   for (size_t index = 0; index < input.size(); ++index) {
     const ThreeBodyPotential_Angle::results_t result =
-        angle( input[index] );
+        angle( FunctionModel::list_of_arguments_t(1, input[index]) );
     CPPUNIT_ASSERT(
         Helpers::isEqual(
@@ -135,5 +135,5 @@
           0.,
           angle.derivative(
-              input[5]
+              FunctionModel::list_of_arguments_t(1, input[5])
           )[0],
           10.
@@ -156,5 +156,5 @@
           0.,
           angle.parameter_derivative(
-              input[5],
+              FunctionModel::list_of_arguments_t(1, input[5]),
               0
           )[0],
@@ -166,5 +166,5 @@
           0.,
           angle.parameter_derivative(
-              input[5],
+              FunctionModel::list_of_arguments_t(1, input[5]),
               1
           )[0],
