Index: src/Potentials/Specifics/SaturationPotential.cpp
===================================================================
--- src/Potentials/Specifics/SaturationPotential.cpp	(revision 035c0da1f58fb5190abeb335461b2da4ebc8375b)
+++ src/Potentials/Specifics/SaturationPotential.cpp	(revision e9bdeff648884b0f7b2819ea260bf685fa1774fa)
@@ -74,5 +74,5 @@
   SerializablePotential(_ParticleTypes),
   morse(_ParticleTypes),
-  angle(symmetrizeTypes(_ParticleTypes)),
+  angle(addSaturationType(_ParticleTypes)),
   energy_offset(0.)
 {
@@ -81,6 +81,6 @@
   ASSERT( _ParticleTypes.size() == (size_t)2,
       "SaturationPotential::SaturationPotential() - exactly two types must be given.");
-  ASSERT( _ParticleTypes[1] == 1,
-      "SaturationPotential::SaturationPotential() - second type must be hydrogen.");
+//  ASSERT( _ParticleTypes[1] == 1,
+//      "SaturationPotential::SaturationPotential() - second type must be hydrogen.");
 }
 
@@ -95,11 +95,11 @@
   SerializablePotential(_ParticleTypes),
   morse(_ParticleTypes),
-  angle(symmetrizeTypes(_ParticleTypes)),
+  angle(addSaturationType(_ParticleTypes)),
   energy_offset(_all_energy_offset)
 {
   ASSERT( _ParticleTypes.size() == (size_t)2,
       "SaturationPotential::SaturationPotential() - exactly two types must be given.");
-  ASSERT( _ParticleTypes[1] == 1,
-      "SaturationPotential::SaturationPotential() - second type must be hydrogen.");
+//  ASSERT( _ParticleTypes[1] == 1,
+//      "SaturationPotential::SaturationPotential() - second type must be hydrogen.");
   parameters_t morse_params(morse.getParameterDimension());
   morse_params[PairPotential_Morse::spring_constant] = _morse_spring_constant;
@@ -269,4 +269,12 @@
   double result = 0.;
   const ParticleTypes_t &morse_types = morse.getParticleTypes();
+  const ParticleTypes_t &angle_types = angle.getParticleTypes();
+  double multiplicity = 1.;
+  if ((angle_types[0] == angle_types[1]) && (angle_types[1] == angle_types[2]))
+    multiplicity = 1./6.;
+  else if ((angle_types[0] == angle_types[1])  
+        || (angle_types[1] == angle_types[2])
+        || (angle_types[0] == angle_types[2]))
+    multiplicity = .5;
   for(arguments_t::const_iterator argiter = arguments.begin();
       argiter != arguments.end();
@@ -283,5 +291,7 @@
       if (result != result)
         ELOG(1, "result is NAN.");
-
+    }
+    if (((r_ij.types.first == angle_types[0]) && (r_ij.types.second == angle_types[1]))
+        || ((r_ij.types.first == angle_types[1]) && (r_ij.types.second == angle_types[0]))) {
       // Angle contribution
       {
@@ -295,5 +305,5 @@
               Extractors::reorderArgumentsByParticleTypes(*iter, angle.getParticleTypes());
           // We get both angles, e.g. 0-4-1 and 1-4-0, hence multiply with 0.5
-          const double tmp = .5*angle(tempargs)[0];  // as we have all distances we get both jk and kj
+          const double tmp = multiplicity*angle(tempargs)[0];  // as we have all distances we get both jk and kj
 //          LOG(3, "DEBUG: angle yields " << tmp << " for << " << tempargs << ".");
           result += tmp;
@@ -325,55 +335,85 @@
 {
   double result = 0.;
-  if (index == all_energy_offset) {
+  switch (index) {
+  case 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());
-              // We get both angles, e.g. 0-4-1 and 1-4-0, hence multiply with 0.5
-              if (index == angle_spring_constant)
-                result += .5*angle.parameter_derivative(tempargs, PairPotential_Angle::spring_constant)[0];
-              else if (index == angle_equilibrium_distance)
-                result += .5*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;
+  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);
+
+         double tmp = 0.;
+          switch (index) {
+            case morse_spring_constant:
+              tmp += morse.parameter_derivative(args, PairPotential_Morse::spring_constant)[0];
+              break;
+            case morse_equilibrium_distance:
+              tmp += morse.parameter_derivative(args, PairPotential_Morse::equilibrium_distance)[0];
+              break;
+            case morse_dissociation_energy:
+              tmp += morse.parameter_derivative(args, PairPotential_Morse::dissociation_energy)[0];
+              break;
+            default:
+              ASSERT(0, "SaturationPotential::parameter_derivative() - We cannot get here.");
+              break;
           }
-            break;
-          default:
-            ASSERT( 0, "SaturationPotential::parameter_derivative() - impossible to get here.");
-            break;
+//          LOG(2, "DEBUG: morse yields " << tmp << " for << " << args << ".");
+          result += tmp;
         }
       }
     }
+    break;
+  case angle_spring_constant:
+  case angle_equilibrium_distance:
+    {
+      const ParticleTypes_t &angle_types = angle.getParticleTypes();
+      double multiplicity = 1.;
+      if ((angle_types[0] == angle_types[1]) && (angle_types[1] == angle_types[2]))
+        multiplicity = 1./6.;
+      else if ((angle_types[0] == angle_types[1])  
+            || (angle_types[1] == angle_types[2])
+            || (angle_types[0] == angle_types[2]))
+        multiplicity = .5;
+      for(arguments_t::const_iterator argiter = arguments.begin();
+          argiter != arguments.end();
+          ++argiter) {
+        const argument_t &r_ij = *argiter;
+        if (((r_ij.types.first == angle_types[0]) && (r_ij.types.second == angle_types[1]))
+            || ((r_ij.types.first == angle_types[1]) && (r_ij.types.second == angle_types[0]))) {
+          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());
+            // We get both angles, e.g. 0-4-1 and 1-4-0, hence multiply with 0.5
+            double tmp = 0.;
+            if (index == angle_spring_constant)
+              tmp += multiplicity*angle.parameter_derivative(tempargs, PairPotential_Angle::spring_constant)[0];
+            else if (index == angle_equilibrium_distance)
+              tmp += multiplicity*angle.parameter_derivative(tempargs, PairPotential_Angle::equilibrium_distance)[0];
+//            LOG(2, "DEBUG: angle yields " << tmp << " for << " << tempargs << ".");
+            result += tmp;
+            if (result != result)
+              ELOG(1, "result is NAN.");
+          }
+        }
+      }
+    }
+    break;
+  default:
+    ASSERT( 0, "SaturationPotential::parameter_derivative() - impossible to get here.");
+    break;
   }
   return SaturationPotential::results_t(1, result);
@@ -396,22 +436,30 @@
 }
 
+const SaturationPotential::ParticleTypes_t
+SaturationPotential::addSaturationType(const ParticleTypes_t &_ParticleTypes)
+{
+  ParticleTypes_t types(_ParticleTypes);
+  types.push_back( ParticleType_t(1) );
+  return types;
+}
+
 FunctionModel::extractor_t
 SaturationPotential::getFragmentSpecificExtractor() const
 {
-  Fragment::charges_t charges;
-  charges.resize(getParticleTypes().size());
-  std::transform(getParticleTypes().begin(), getParticleTypes().end(),
-      charges.begin(), boost::lambda::_1);	  
+//  Fragment::charges_t charges;
+//  charges.resize(getParticleTypes().size());
+//  std::transform(getParticleTypes().begin(), getParticleTypes().end(),
+//      charges.begin(), boost::lambda::_1);	  
   FunctionModel::extractor_t returnfunction;
-  if (charges[0] == charges[1]) {
-    // In case both types are equal there is only a single pair of possible
-    // type combinations.
-     returnfunction =
-        boost::bind(&Extractors::gatherAllDistancesFromFragment,
-            boost::bind(&Fragment::getPositions, _1),
-            boost::bind(&Fragment::getCharges, _1),
-            charges, // is only temporarily created, hence copy
-            _2);
-  } else {
+//  if (charges[0] == charges[1]) {
+//    // In case both types are equal there is only a single pair of possible
+//    // type combinations.
+//     returnfunction =
+//        boost::bind(&Extractors::gatherAllDistancesFromFragment,
+//            boost::bind(&Fragment::getPositions, _1),
+//            boost::bind(&Fragment::getCharges, _1),
+//            charges, // is only temporarily created, hence copy
+//            _2);
+//  } else {
     // we have to chain here a rather complex "tree" of functions
     // as we only have a couple of ParticleTypes but need to get
@@ -419,33 +467,45 @@
     // 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
-    Fragment::charges_t secondpair(2, boost::cref(charges[1]));
-    const Fragment::charges_t &thirdpair = charges;
+//    Fragment::charges_t secondpair(2, boost::cref(charges[1]));
+//    const Fragment::charges_t &thirdpair = charges;
+    Fragment::charges_t charges_angle;
+    {
+      charges_angle.resize(angle.getParticleTypes().size());
+      std::transform(angle.getParticleTypes().begin(), angle.getParticleTypes().end(),
+          charges_angle.begin(), boost::lambda::_1);
+    }
+    Fragment::charges_t charges_morse;
+    {
+      charges_morse.resize(morse.getParticleTypes().size());
+      std::transform(morse.getParticleTypes().begin(), morse.getParticleTypes().end(),
+          charges_morse.begin(), boost::lambda::_1);
+    }
     returnfunction =
 //        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::gatherAllDistancesFromFragment,
                   boost::bind(&Fragment::getPositions, _1),
                   boost::bind(&Fragment::getCharges, _1),
-                  secondpair,  // no crefs here as are temporaries!
+                  charges_angle,  // no crefs here as are temporaries!
                   _2),
-//            ),
-            boost::bind(&Extractors::gatherAllDistancesFromFragment,
-                boost::bind(&Fragment::getPositions, _1),
-                boost::bind(&Fragment::getCharges, _1),
-                boost::cref(thirdpair), // only the last one is no temporary
-                _2)
-          );
-//          boost::cref(angle.getParticleTypes())
-//        );
-}
+              boost::bind(&Extractors::gatherAllDistancesFromFragment,
+                  boost::bind(&Fragment::getPositions, _1),
+                  boost::bind(&Fragment::getCharges, _1),
+                  charges_morse,  // no crefs here as are temporaries!
+                  _2)
+//            )
+//            boost::bind(&Extractors::gatherAllDistancesFromFragment,
+//                boost::bind(&Fragment::getPositions, _1),
+//                boost::bind(&Fragment::getCharges, _1),
+//                boost::cref(thirdpair), // only the last one is no temporary
+//                _2)
+//          ),
+//          boost::bind(&PairPotential_Angle::getParticleTypes, boost::cref(angle))
+        );
+//  }
   return returnfunction;
 }
Index: src/Potentials/Specifics/SaturationPotential.hpp
===================================================================
--- src/Potentials/Specifics/SaturationPotential.hpp	(revision 035c0da1f58fb5190abeb335461b2da4ebc8375b)
+++ src/Potentials/Specifics/SaturationPotential.hpp	(revision e9bdeff648884b0f7b2819ea260bf685fa1774fa)
@@ -203,4 +203,11 @@
   static const ParticleTypes_t symmetrizeTypes(const ParticleTypes_t &_ParticleTypes);
 
+  /** Adds saturation type bond partner.
+   *
+   * @param _ParticleTypes types from a pair potential
+   * @return types for a saturation type-enhanced triple potential
+   */
+  static const ParticleTypes_t addSaturationType(const ParticleTypes_t &_ParticleTypes);
+
 private:
   PairPotential_Morse morse;
