Index: src/FunctionApproximation/Extractors.cpp
===================================================================
--- src/FunctionApproximation/Extractors.cpp	(revision 0f5d38e3ab5cb4ea3ef27e9a7fab1d2a07dc2437)
+++ src/FunctionApproximation/Extractors.cpp	(revision 51e0e3079be933cbf6a3e3dc4a527a0575e2d1fc)
@@ -45,5 +45,7 @@
 
 #include "CodePatterns/Assert.hpp"
+#include "CodePatterns/IteratorAdaptors.hpp"
 #include "CodePatterns/Log.hpp"
+#include "CodePatterns/toString.hpp"
 
 #include "LinearAlgebra/Vector.hpp"
@@ -538,5 +540,188 @@
 }
 
+
+struct ParticleTypesComparator {
+  bool operator()(const argument_t::types_t &a, const argument_t::types_t &b)
+  {
+    if (a.first < a.second) {
+      if (b.first < b.second) {
+        if (a.first < b.first)
+          return true;
+        else if (a.first > b.first)
+          return false;
+        else
+          return (a.second < b.second);
+      } else {
+        if (a.first < b.second)
+          return true;
+        else if (a.first > b.second)
+          return false;
+        else
+          return (a.second < b.first);
+      }
+    } else {
+      if (b.first < b.second) {
+        if (a.second < b.first)
+          return true;
+        else if (a.second > b.first)
+          return false;
+        else
+          return (a.first < b.second);
+      } else {
+        if (a.second < b.second)
+          return true;
+        else if (a.second > b.second)
+          return false;
+        else
+          return (a.first < b.first);
+      }
+    }
+  }
+};
+
+std::ostream& operator<<(std::ostream &out, const argument_t::types_t &a)
+{
+  out << "[" << a.first << "," << a.second << "]";
+  return out;
+}
+
 FunctionModel::arguments_t Extractors::reorderArgumentsByParticleTypes(
+    const FunctionModel::arguments_t &args,
+    const ParticleTypes_t &_types
+    )
+{
+  /// We  place all arguments into multimap according to particle type pair.
+  // here, we need a special comparator such that types in key pair are always
+  // properly ordered.
+  typedef std::multimap<
+      argument_t::types_t,
+      argument_t,
+      ParticleTypesComparator> TypePair_Argument_Map_t;
+  TypePair_Argument_Map_t argument_map;
+  for(FunctionModel::arguments_t::const_iterator iter = args.begin();
+      iter != args.end(); ++iter) {
+    argument_map.insert( std::make_pair(iter->types, *iter) );
+  }
+  LOG(4, "DEBUG: particle_type map is " << argument_map << ".");
+
+  /// Then, we create the desired unique keys
+  typedef std::vector<argument_t::types_t> UniqueTypes_t;
+  UniqueTypes_t UniqueTypes;
+  for (ParticleTypes_t::const_iterator firstiter = _types.begin();
+      firstiter != _types.end();
+      ++firstiter) {
+    for (ParticleTypes_t::const_iterator seconditer = firstiter;
+        seconditer != _types.end();
+        ++seconditer) {
+      if (seconditer == firstiter)
+        continue;
+      UniqueTypes.push_back( std::make_pair(*firstiter, *seconditer) );
+    }
+  }
+  LOG(4, "DEBUG: Created unique types as keys " << UniqueTypes << ".");
+
+  /// Finally, we use the unique key list to pick corresponding arguments from the map
+  FunctionModel::arguments_t returnargs;
+  returnargs.reserve(args.size());
+  while (!argument_map.empty()) {
+    // note that particle_types_t may be flipped, i.e. 1,8 is equal to 8,1, but we
+    // must maintain the correct order in indices in accordance with the order
+    // in _types, i.e. 1,8,1 must match with e.g. ids 1,0,2 where 1 has type 1,
+    // 0 has type 8, and 2 has type 2.
+    // In other words: We do not want to flip/modify arguments such that they match
+    // with the specific type pair we seek but then this comes at the price that we
+    // have flip indices when the types in a pair are flipped.
+
+    typedef std::vector<size_t> indices_t;
+    //!> here, we gather the indices as we discover them
+    indices_t indices;
+    indices.resize(_types.size(), (size_t)-1);
+
+    // these are two iterators that create index pairs in the same way as we have
+    // created type pairs. If a -1 is still present in indices, then the index is
+    // still arbitrary but is then set by the next found index
+    indices_t::iterator firstindex = indices.begin();
+    indices_t::iterator secondindex = firstindex+1;
+
+    //!> here, we gather the current bunch of arguments as we find them
+    FunctionModel::arguments_t argumentbunch;
+    argumentbunch.reserve(UniqueTypes.size());
+
+    for (UniqueTypes_t::const_iterator typeiter = UniqueTypes.begin();
+        typeiter != UniqueTypes.end(); ++typeiter) {
+      // have all arguments to same type pair as list within the found range
+      std::pair<
+          TypePair_Argument_Map_t::iterator,
+          TypePair_Argument_Map_t::iterator> range_t =
+              argument_map.equal_range(*typeiter);
+      LOG(4, "DEBUG: Set of arguments to current key [" << typeiter->first << ","
+          << typeiter->second << "] is " << std::list<argument_t>(
+              MapValueIterator<TypePair_Argument_Map_t::iterator>(range_t.first),
+              MapValueIterator<TypePair_Argument_Map_t::iterator>(range_t.second)
+              ) << ".");
+      // the first key is always easy and is pivot which the rest has to be associated to
+      if (typeiter == UniqueTypes.begin()) {
+        const argument_t & arg = range_t.first->second;
+        if ((typeiter->first == arg.types.first) && (typeiter->second == arg.types.second)) {
+          // store in correct order
+          *firstindex = arg.indices.first;
+          *secondindex = arg.indices.second;
+        } else {
+          // store in flipped order
+          *firstindex = arg.indices.second;
+          *secondindex = arg.indices.first;
+        }
+        argumentbunch.push_back(arg);
+        argument_map.erase(range_t.first);
+        LOG(4, "DEBUG: Gathered first argument " << arg << ".");
+      } else {
+        // go through the range and pick the first argument matching the index constraints
+        for (TypePair_Argument_Map_t::iterator argiter = range_t.first;
+            argiter != range_t.second; ++argiter) {
+          // seconditer may be -1 still
+          const argument_t &arg = argiter->second;
+          if (arg.indices.first == *firstindex) {
+            if ((arg.indices.second == *secondindex) || (*secondindex == (size_t)-1)) {
+              if (*secondindex == -1)
+                *secondindex = arg.indices.second;
+              argumentbunch.push_back(arg);
+              argument_map.erase(argiter);
+              LOG(4, "DEBUG: Gathered another argument " << arg << ".");
+              break;
+            }
+          } else if ((arg.indices.first == *secondindex) || (*secondindex == (size_t)-1)) {
+            if (arg.indices.second == *firstindex) {
+              if (*secondindex == (size_t)-1)
+                *secondindex = arg.indices.first;
+              argumentbunch.push_back(arg);
+              argument_map.erase(argiter);
+              LOG(4, "DEBUG: Gathered another (flipped) argument " << arg << ".");
+              break;
+            }
+          }
+        }
+      }
+      // move along in indices and check bounds
+      ++secondindex;
+      if (secondindex == indices.end()) {
+        ++firstindex;
+        if (firstindex != indices.end()-1)
+          secondindex = firstindex+1;
+      }
+    }
+    ASSERT( (firstindex == indices.end()-1) && (secondindex == indices.end()),
+        "Extractors::reorderArgumentsByParticleTypes() - we have not gathered enough arguments.");
+    ASSERT( argumentbunch.size() == UniqueTypes.size(),
+        "Extractors::reorderArgumentsByParticleTypes() - we have not gathered enough arguments.");
+    // place bunch of arguments in return args
+    LOG(3, "DEBUG: Given types " << _types << " and found indices " << indices << ".");
+    LOG(3, "DEBUG: Final bunch of arguments is " << argumentbunch << ".");
+    returnargs.insert(returnargs.end(), argumentbunch.begin(), argumentbunch.end());
+  }
+
+  return returnargs;
+}
+
+FunctionModel::arguments_t Extractors::filterArgumentsByParticleTypes(
     const FunctionModel::arguments_t &args,
     const ParticleTypes_t &_types
@@ -557,5 +742,5 @@
 
   // basically, we have two choose any two pairs out of types but only those
-  // where the first is less than the letter. Hence, we start the second
+  // where the first is less than the latter. Hence, we start the second
   // iterator at the current position of the first one and skip the equal case.
   for (ParticleTypes_t::const_iterator firstiter = _types.begin();
@@ -580,9 +765,13 @@
         else if ((iter->types.first == *seconditer)
               && (iter->types.second == *firstiter)) {
-          argument_t flippedtypes(*iter);
-          std::swap( flippedtypes.indices.first, flippedtypes.indices.second );
-          std::swap( flippedtypes.types.first, flippedtypes.types.second );
-          returnargs.push_back( flippedtypes );
-          break;
+//          argument_t flippedtypes(*iter);
+//          std::swap( flippedtypes.indices.first, flippedtypes.indices.second );
+//          std::swap( flippedtypes.types.first, flippedtypes.types.second );
+          returnargs.push_back( *iter );
+          iter = availableList.erase(iter);
+          LOG(4, "DEBUG: Accepted (flipped) argument.");
+        } else {
+          ++iter;
+          LOG(4, "DEBUG: Rejected argument.");
         }
       }
Index: src/FunctionApproximation/Extractors.hpp
===================================================================
--- src/FunctionApproximation/Extractors.hpp	(revision 0f5d38e3ab5cb4ea3ef27e9a7fab1d2a07dc2437)
+++ src/FunctionApproximation/Extractors.hpp	(revision 51e0e3079be933cbf6a3e3dc4a527a0575e2d1fc)
@@ -301,8 +301,10 @@
       );
 
-  /** Reorder arguments according to types.
-   *
-   * If particle types is (0,1,2) and three argumens, each with a pair of types,
-   * are given, then the alignment will be: (0,1), (0,2), and (1,2).
+  /** Reorder the arguments to bring adjacent ones together.
+   *
+   * After filtering via particle types arguments related via same indices
+   * must not necessarily be contained in the same bunch. This reordering
+   * is done here, preserving the alignment given in
+   * \sa filterArgumentsByParticleTypes()
    *
    * \param args arguments to reorder
@@ -315,4 +317,18 @@
       );
 
+  /** Filter arguments according to types.
+   *
+   * If particle types is (0,1,2) and three arguments, each with a pair of types,
+   * are given, then the alignment will be: (0,1), (0,2), and (1,2).
+   *
+   * \param args arguments to reorder
+   * \param _types particle type vector
+   * \return filtered args
+   */
+  FunctionModel::arguments_t filterArgumentsByParticleTypes(
+      const FunctionModel::arguments_t &args,
+      const ParticleTypes_t &_types
+      );
+
   /** Combines two argument lists by sorting and making unique.
    *
Index: src/Potentials/Specifics/ConstantPotential.cpp
===================================================================
--- src/Potentials/Specifics/ConstantPotential.cpp	(revision 0f5d38e3ab5cb4ea3ef27e9a7fab1d2a07dc2437)
+++ src/Potentials/Specifics/ConstantPotential.cpp	(revision 51e0e3079be933cbf6a3e3dc4a527a0575e2d1fc)
@@ -180,6 +180,7 @@
 {
   FunctionModel::filter_t returnfunction =
-      boost::bind(&Helpers::NoOp_Filterfunction,
-          _1);
+      boost::bind(&Extractors::filterArgumentsByParticleTypes,
+          _1,
+          getParticleTypes());
   return returnfunction;
 }
Index: src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp
===================================================================
--- src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp	(revision 0f5d38e3ab5cb4ea3ef27e9a7fab1d2a07dc2437)
+++ src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp	(revision 51e0e3079be933cbf6a3e3dc4a527a0575e2d1fc)
@@ -701,6 +701,7 @@
 {
   FunctionModel::filter_t returnfunction =
-      boost::bind(&Helpers::NoOp_Filterfunction,
-          _1);
+      boost::bind(&Extractors::filterArgumentsByParticleTypes,
+          _1,
+          getParticleTypes());
   return returnfunction;
 }
Index: src/Potentials/Specifics/PairPotential_Angle.cpp
===================================================================
--- src/Potentials/Specifics/PairPotential_Angle.cpp	(revision 0f5d38e3ab5cb4ea3ef27e9a7fab1d2a07dc2437)
+++ src/Potentials/Specifics/PairPotential_Angle.cpp	(revision 51e0e3079be933cbf6a3e3dc4a527a0575e2d1fc)
@@ -225,6 +225,8 @@
   FunctionModel::filter_t returnfunction =
       boost::bind(&Extractors::reorderArgumentsByParticleTypes,
-          _1,
-          getParticleTypes()
+        boost::bind(&Extractors::filterArgumentsByParticleTypes,
+            _1,
+            getParticleTypes()),
+        getParticleTypes()
       );
   return returnfunction;
Index: src/Potentials/Specifics/PairPotential_Harmonic.cpp
===================================================================
--- src/Potentials/Specifics/PairPotential_Harmonic.cpp	(revision 0f5d38e3ab5cb4ea3ef27e9a7fab1d2a07dc2437)
+++ src/Potentials/Specifics/PairPotential_Harmonic.cpp	(revision 51e0e3079be933cbf6a3e3dc4a527a0575e2d1fc)
@@ -201,6 +201,7 @@
 {
   FunctionModel::filter_t returnfunction =
-      boost::bind(&Helpers::NoOp_Filterfunction,
-          _1);
+      boost::bind(&Extractors::filterArgumentsByParticleTypes,
+          _1,
+          getParticleTypes());
   return returnfunction;
 }
Index: src/Potentials/Specifics/PairPotential_Morse.cpp
===================================================================
--- src/Potentials/Specifics/PairPotential_Morse.cpp	(revision 0f5d38e3ab5cb4ea3ef27e9a7fab1d2a07dc2437)
+++ src/Potentials/Specifics/PairPotential_Morse.cpp	(revision 51e0e3079be933cbf6a3e3dc4a527a0575e2d1fc)
@@ -230,6 +230,7 @@
 {
   FunctionModel::filter_t returnfunction =
-      boost::bind(&Helpers::NoOp_Filterfunction,
-          _1);
+      boost::bind(&Extractors::filterArgumentsByParticleTypes,
+          _1,
+          getParticleTypes());
   return returnfunction;
 }
