Index: src/LevMartester.cpp
===================================================================
--- src/LevMartester.cpp	(revision 775dd1aa628462c5a79c38c46bfea677334fad29)
+++ src/LevMartester.cpp	(revision 0dad5f9b6ef661b1712806b83f1f1f3da56f8bb3)
@@ -369,5 +369,6 @@
 
       SaturationPotential saturation(types);
-      saturation.setTriplefunction(2.5, triplefunction);
+      // saturation has internal triplefunction now
+//      saturation.setTriplefunction(2.5, triplefunction);
       FunctionModel &model = saturation;
       FunctionApproximation approximator(TersoffData, model); // CH4 has 5 atoms, hence 5*4/2 distances
Index: src/Potentials/Specifics/SaturationPotential.cpp
===================================================================
--- src/Potentials/Specifics/SaturationPotential.cpp	(revision 775dd1aa628462c5a79c38c46bfea677334fad29)
+++ src/Potentials/Specifics/SaturationPotential.cpp	(revision 0dad5f9b6ef661b1712806b83f1f1f3da56f8bb3)
@@ -176,4 +176,75 @@
 }
 
+std::vector<FunctionModel::arguments_t>
+triplefunction(
+    const argument_t &argument,
+    const FunctionModel::arguments_t& args)
+{
+  const size_t firstindex = argument.indices.first;
+  const size_t secondindex = argument.indices.second;
+//  LOG(2, "DEBUG: first index is " << firstindex << ", second index is " << secondindex << ".");
+
+  // place all arguments that share either index into a lookup map
+  typedef std::map< size_t, FunctionModel::arguments_t::const_iterator > IndexLookup_t;
+  IndexLookup_t LookuptoFirst;
+  IndexLookup_t LookuptoSecond;
+  for (FunctionModel::arguments_t::const_iterator iter = args.begin();
+      iter != args.end();
+      ++iter) {
+    if (((*iter).indices.first == argument.indices.first)
+        && ((*iter).indices.second == argument.indices.second))
+      continue;
+    if (firstindex == (*iter).indices.first) {
+      LookuptoFirst.insert( std::make_pair( (*iter).indices.second, iter) );
+    }
+    else if (firstindex == (*iter).indices.second) {
+      LookuptoFirst.insert( std::make_pair( (*iter).indices.first, iter) );
+    }
+    if (secondindex == (*iter).indices.first) {
+      LookuptoSecond.insert( std::make_pair( (*iter).indices.second, iter) );
+    }
+    else if (secondindex == (*iter).indices.second) {
+      LookuptoSecond.insert( std::make_pair((*iter).indices.first, iter) );
+    }
+  }
+//  {
+//    std::stringstream lookupstream;
+//    for (IndexLookup_t::const_iterator iter = LookuptoFirst.begin();
+//        iter != LookuptoFirst.end();
+//        ++iter) {
+//      lookupstream << "(" << iter->first << "," << *(iter->second) << ") ";
+//    }
+//    LOG(2, "DEBUG: LookupToFirst is " << lookupstream.str() << ".");
+//  }
+//  {
+//    std::stringstream lookupstream;
+//    for (IndexLookup_t::const_iterator iter = LookuptoSecond.begin();
+//        iter != LookuptoSecond.end();
+//        ++iter) {
+//      lookupstream << "(" << iter->first << "," << *(iter->second) << ") ";
+//    }
+//    LOG(2, "DEBUG: LookuptoSecond is " << lookupstream.str() << ".");
+//  }
+
+  // now go through the first lookup as the second argument and pick the
+  // corresponding third argument by the matching index
+  std::vector<FunctionModel::arguments_t> results;
+  for (IndexLookup_t::const_iterator iter = LookuptoFirst.begin();
+      iter != LookuptoFirst.end();
+      ++iter) {
+    IndexLookup_t::const_iterator otheriter = LookuptoSecond.find(iter->first);
+    ASSERT( otheriter != LookuptoSecond.end(),
+        "triplefunction() - cannot find index "+toString(iter->first)
+        +" in LookupToSecond");
+    FunctionModel::arguments_t result(1, argument);
+    result.reserve(3);
+    result.push_back(*(iter->second));
+    result.push_back(*(otheriter->second));
+    results.push_back(result);
+  }
+
+  return results;
+}
+
 SaturationPotential::results_t
 SaturationPotential::operator()(
@@ -197,13 +268,23 @@
       if (result != result)
         ELOG(1, "result is NAN.");
-    }
-  }
-  {
-    // Angle contribution
-    const double tmp = angle(arguments)[0];  // as we have all distances we get both jk and kj
-//    LOG(2, "DEBUG: angle yields " << tmp << " for << " << arguments << ".");
-    result += tmp;
-    if (result != result)
-      ELOG(1, "result is NAN.");
+
+      // Angle contribution
+      {
+        typedef std::vector<FunctionModel::arguments_t> tripleargs_t;
+        tripleargs_t tripleargs =
+            triplefunction(r_ij, arguments);
+        for (tripleargs_t::const_iterator iter = tripleargs.begin();
+            iter != tripleargs.end();
+            ++iter) {
+          FunctionModel::arguments_t tempargs =
+              Extractors::reorderArgumentsByParticleTypes(*iter, angle.getParticleTypes());
+          const double tmp = angle(tempargs)[0];  // as we have all distances we get both jk and kj
+//          LOG(2, "DEBUG: angle yields " << tmp << " for << " << arguments << ".");
+          result += tmp;
+          if (result != result)
+            ELOG(1, "result is NAN.");
+        }
+      }
+    }
   }
   return std::vector<result_t>(1, energy_offset + result);
@@ -228,53 +309,54 @@
 {
   double result = 0.;
-  switch (index) {
-    case all_energy_offset:
-    {
-      result = 1.;
-      break;
-    }
-    case morse_spring_constant:
-    case morse_equilibrium_distance:
-    case morse_dissociation_energy:
-    {
-      const ParticleTypes_t &morse_types = morse.getParticleTypes();
-      for(arguments_t::const_iterator argiter = arguments.begin();
-          argiter != arguments.end();
-          ++argiter) {
-        const argument_t &r_ij = *argiter;
-        if (((r_ij.types.first == morse_types[0]) && (r_ij.types.second == morse_types[1]))
-            || ((r_ij.types.first == morse_types[1]) && (r_ij.types.second == morse_types[0]))) {
-          arguments_t args(1, r_ij);
-          switch (index) {
-            case morse_spring_constant:
-              result += morse.parameter_derivative(args, PairPotential_Morse::spring_constant)[0];
-              break;
-            case morse_equilibrium_distance:
-              result += morse.parameter_derivative(args, PairPotential_Morse::equilibrium_distance)[0];
-              break;
-            case morse_dissociation_energy:
-              result += morse.parameter_derivative(args, PairPotential_Morse::dissociation_energy)[0];
-              break;
-            default:
-              ASSERT( 0, "SaturationPotential::parameter_derivative() - impossible to get here.");
-              break;
+  if (index == all_energy_offset) {
+    result = 1.;
+  } else {
+    const ParticleTypes_t &morse_types = morse.getParticleTypes();
+    for(arguments_t::const_iterator argiter = arguments.begin();
+        argiter != arguments.end();
+        ++argiter) {
+      const argument_t &r_ij = *argiter;
+      if (((r_ij.types.first == morse_types[0]) && (r_ij.types.second == morse_types[1]))
+          || ((r_ij.types.first == morse_types[1]) && (r_ij.types.second == morse_types[0]))) {
+        arguments_t args(1, r_ij);
+
+        switch (index) {
+          case morse_spring_constant:
+            result += morse.parameter_derivative(args, PairPotential_Morse::spring_constant)[0];
+            break;
+          case morse_equilibrium_distance:
+            result += morse.parameter_derivative(args, PairPotential_Morse::equilibrium_distance)[0];
+            break;
+          case morse_dissociation_energy:
+            result += morse.parameter_derivative(args, PairPotential_Morse::dissociation_energy)[0];
+            break;
+          case angle_spring_constant:
+          case angle_equilibrium_distance:
+          {
+            typedef std::vector<FunctionModel::arguments_t> tripleargs_t;
+            tripleargs_t tripleargs =
+                triplefunction(r_ij, arguments);
+            for (tripleargs_t::const_iterator iter = tripleargs.begin();
+                iter != tripleargs.end();
+                ++iter) {
+              FunctionModel::arguments_t tempargs =
+                  Extractors::reorderArgumentsByParticleTypes(*iter, angle.getParticleTypes());
+              if (index == angle_spring_constant)
+                result += angle.parameter_derivative(tempargs, PairPotential_Angle::spring_constant)[0];
+              else if (index == angle_equilibrium_distance)
+                result += angle.parameter_derivative(tempargs, PairPotential_Angle::equilibrium_distance)[0];
+//              LOG(2, "DEBUG: angle yields " << tmp << " for << " << arguments << ".");
+//              result += tmp;
+              if (result != result)
+                ELOG(1, "result is NAN.");
+            }
           }
+            break;
+          default:
+            ASSERT( 0, "SaturationPotential::parameter_derivative() - impossible to get here.");
+            break;
         }
       }
-      break;
-    }
-    case angle_spring_constant:
-    {
-      result = angle.parameter_derivative(arguments, PairPotential_Angle::spring_constant)[0];
-      break;
-    }
-    case angle_equilibrium_distance:
-    {
-      result = angle.parameter_derivative(arguments, PairPotential_Angle::equilibrium_distance)[0];
-      break;
-    }
-    default:
-      ELOG(1, "SaturationPotential::parameter_derivative() - index " << index << " invalid.");
-      break;
+    }
   }
   return SaturationPotential::results_t(1, result);
@@ -333,22 +415,24 @@
     // Finally, we also need to arrange them in correct order
     // (for PairPotentiale_Angle).
-    charges_t firstpair(2, boost::cref(charges[0]));
+    //charges_t firstpair(2, boost::cref(charges[0]));
+    // only that saturation potential never has its middle element twice!
+    // hence, we skip the firstpair but keep the code for later generalization
     charges_t secondpair(2, boost::cref(charges[1]));
     const charges_t &thirdpair = charges;
     returnfunction =
-        boost::bind(&Extractors::reorderArgumentsByParticleTypes,
+//        boost::bind(&Extractors::reorderArgumentsByParticleTypes,
           boost::bind(&Extractors::combineArguments,
-            boost::bind(&Extractors::combineArguments,
-              boost::bind(&Extractors::gatherAllDistancesFromFragment,
-                  boost::bind(&Fragment::getPositions, _1),
-                  boost::bind(&Fragment::getCharges, _1),
-                  firstpair,  // no crefs here as are temporaries!
-                  _2),
+//            boost::bind(&Extractors::combineArguments,
+//              boost::bind(&Extractors::gatherAllDistancesFromFragment,
+//                  boost::bind(&Fragment::getPositions, _1),
+//                  boost::bind(&Fragment::getCharges, _1),
+//                  firstpair,  // no crefs here as are temporaries!
+//                  _2),
               boost::bind(&Extractors::gatherAllDistancesFromFragment,
                   boost::bind(&Fragment::getPositions, _1),
                   boost::bind(&Fragment::getCharges, _1),
                   secondpair,  // no crefs here as are temporaries!
-                  _2)
-            ),
+                  _2),
+//            ),
             boost::bind(&Extractors::gatherAllDistancesFromFragment,
                 boost::bind(&Fragment::getPositions, _1),
@@ -356,7 +440,7 @@
                 boost::cref(thirdpair), // only the last one is no temporary
                 _2)
-          ),
-          boost::cref(angle.getParticleTypes())
-        );
+          );
+//          boost::cref(angle.getParticleTypes())
+//        );
 }
   return returnfunction;
Index: src/Potentials/Specifics/SaturationPotential.hpp
===================================================================
--- src/Potentials/Specifics/SaturationPotential.hpp	(revision 775dd1aa628462c5a79c38c46bfea677334fad29)
+++ src/Potentials/Specifics/SaturationPotential.hpp	(revision 0dad5f9b6ef661b1712806b83f1f1f3da56f8bb3)
@@ -90,4 +90,6 @@
   /** Sets the magic triple function that we use for getting angle distances.
    *
+   * \note does nothing as we have an internal triplefunction
+   *
    * \param _saturation_cutoff cutoff for the triplefunction
    * @param _triplefunction function that returns a list of triples (i.e. the
@@ -100,6 +102,4 @@
       )
   {
-    saturation_cutoff = _saturation_cutoff;
-    triplefunction = _triplefunction;
   }
 
@@ -207,8 +207,4 @@
   double energy_offset;
 
-  //!> bound function that obtains the triples for the internal coordinationb summation.
-  boost::function< std::vector< arguments_t >(const argument_t &, const double)> triplefunction;
-  double saturation_cutoff;
-
   //!> static definitions of the parameter name for this potential
   static const ParameterNames_t ParameterNames;
