Index: src/Potentials/Specifics/SaturationPotential.cpp
===================================================================
--- src/Potentials/Specifics/SaturationPotential.cpp	(revision dd966f7cb4448d0fa69df1cba9f9aef8d1e54baa)
+++ src/Potentials/Specifics/SaturationPotential.cpp	(revision edac7ccd053ce7f11d7cd70e65c9058c8ac9a273)
@@ -45,8 +45,9 @@
 
 SaturationPotential::SaturationPotential(
+    const double _saturation_cutoff,
     boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction) :
   energy_offset(0.),
   triplefunction(_triplefunction),
-  saturation_cutoff(2.5)
+  saturation_cutoff(_saturation_cutoff)
 {}
 
@@ -58,8 +59,9 @@
     const double _angle_equilibrium_distance,
     const double _all_energy_offset,
+    const double _saturation_cutoff,
     boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction) :
   energy_offset(_all_energy_offset),
   triplefunction(_triplefunction),
-  saturation_cutoff(2.5)
+  saturation_cutoff(_saturation_cutoff)
 {
   parameters_t morse_params(morse.getParameterDimension());
@@ -148,21 +150,24 @@
       ++argiter) {
     const argument_t &r_ij = *argiter;
-    arguments_t args(1, r_ij);
-    if ((r_ij.indices.first == 0)) { // || (r_ij.indices.second == 0)) {
+    if ((r_ij.indices.first == 0)) { // first item must be the non-hydrogen
+      arguments_t args(1, r_ij);
+
+      // Morse contribution
       result += morse(args)[0];
       if (result != result)
         ELOG(1, "result is NAN.");
+
+      // Angle contribution
       std::vector<arguments_t> triples = triplefunction(r_ij, saturation_cutoff);
+      args.resize(3, r_ij);
       for (std::vector<arguments_t>::const_iterator iter = triples.begin();
           iter != triples.end(); ++iter) {
         ASSERT( iter->size() == 2,
-            "ManyBodyPotential_Tersoff::function_derivative_c() - the triples result must contain exactly two distances.");
+            "SaturationPotential::function_derivative_c() - the triples result must contain exactly two distances.");
         const argument_t &r_ik = (*iter)[0];
         const argument_t &r_jk = (*iter)[1];
-        arguments_t args;
-        args.push_back(r_ij);
-        args.push_back(r_ik);
-        args.push_back(r_jk);
-        result += .5 * angle(args)[0];
+        args[1] = r_ik;
+        args[2] = r_jk;
+        result += .5*angle(args)[0];  // as we have all distances we get both jk and kj
         if (result != result)
           ELOG(1, "result is NAN.");
@@ -191,53 +196,63 @@
 {
   double result = 0.;
-  for(arguments_t::const_iterator argiter = arguments.begin();
-      argiter != arguments.end();
-      ++argiter) {
-    const argument_t &r_ij = *argiter;
-    arguments_t args(1, r_ij);
-    if ((r_ij.indices.first == 0)) { // || (r_ij.indices.second == 0)) {
-      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:
-        {
-          std::vector<arguments_t> triples = triplefunction(r_ij, saturation_cutoff);
-          for (std::vector<arguments_t>::const_iterator iter = triples.begin();
-              iter != triples.end(); ++iter) {
-            ASSERT( iter->size() == 2,
-                "ManyBodyPotential_Tersoff::function_derivative_c() - the triples result must contain exactly two distances.");
-            const argument_t &r_ik = (*iter)[0];
-            const argument_t &r_jk = (*iter)[1];
-            arguments_t args;
-            args.push_back(r_ij);
-            args.push_back(r_ik);
-            args.push_back(r_jk);
-            switch (index) {
-              case angle_spring_constant:
-              {
-                result += .5*angle.parameter_derivative(args, PairPotential_Angle::spring_constant)[0];
-                break;
+  if (index == all_energy_offset) {
+    result = 1.;
+  } else {
+    for(arguments_t::const_iterator argiter = arguments.begin();
+        argiter != arguments.end();
+        ++argiter) {
+      const argument_t &r_ij = *argiter;
+      if ((r_ij.indices.first == 0)) { // first item must be the non-hydrogen
+        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:
+          {
+            args.resize(3, r_ij);
+            std::vector<arguments_t> triples = triplefunction(r_ij, saturation_cutoff);
+            for (std::vector<arguments_t>::const_iterator iter = triples.begin();
+                iter != triples.end(); ++iter) {
+              ASSERT( iter->size() == 2,
+                  "SaturationPotential::parameter_derivative() - the triples result must contain exactly two distances.");
+              const argument_t &r_ik = (*iter)[0];
+              ASSERT( r_ik.indices.first == r_ij.indices.first,
+                  "SaturationPotential::parameter_derivative() - i not same in ij, ik.");
+              const argument_t &r_jk = (*iter)[1];
+              ASSERT( r_jk.indices.first == r_ij.indices.second,
+                  "SaturationPotential::parameter_derivative() - j not same in ij, jk.");
+              ASSERT( r_ik.indices.second == r_jk.indices.second,
+                  "SaturationPotential::parameter_derivative() - k not same in ik, jk.");
+              args[1] = r_ik;
+              args[2] = r_jk;
+              switch (index) {   // .5 due to we have all distances we get both jk and kj
+                case angle_spring_constant:
+                {
+                  result += .5*angle.parameter_derivative(args, PairPotential_Angle::spring_constant)[0];
+                  break;
+                }
+                case angle_equilibrium_distance:
+                {
+                  result += .5*angle.parameter_derivative(args, PairPotential_Angle::equilibrium_distance)[0];
+                  break;
+                }
+                default:
+                  break;
               }
-              case angle_equilibrium_distance:
-              {
-                result += .5*angle.parameter_derivative(args, PairPotential_Angle::equilibrium_distance)[0];
-                break;
-              }
-              default:
-                break;
             }
+            break;
           }
         }
@@ -245,6 +260,4 @@
     }
   }
-  if (index == all_energy_offset)
-    result = 1.;
   return SaturationPotential::results_t(1, result);
 }
Index: src/Potentials/Specifics/SaturationPotential.hpp
===================================================================
--- src/Potentials/Specifics/SaturationPotential.hpp	(revision dd966f7cb4448d0fa69df1cba9f9aef8d1e54baa)
+++ src/Potentials/Specifics/SaturationPotential.hpp	(revision edac7ccd053ce7f11d7cd70e65c9058c8ac9a273)
@@ -38,4 +38,5 @@
 public:
   SaturationPotential(
+      const double _saturation_cutoff,
       boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction
       );
@@ -47,4 +48,5 @@
       const double _angle_equilibrium_distance,
       const double _energy_offset,
+      const double _saturation_cutoff,
       boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction
       );
@@ -127,4 +129,5 @@
 
   enum parameter_enum_t {
+    all_energy_offset,
     morse_spring_constant,
     morse_equilibrium_distance,
@@ -132,5 +135,4 @@
     angle_spring_constant,
     angle_equilibrium_distance,
-    all_energy_offset,
     MAXPARAMS
   };
