Index: src/FunctionApproximation/Extractors.cpp
===================================================================
--- src/FunctionApproximation/Extractors.cpp	(revision f5dbea01dfa29040cfe86148e9597b8da56187fd)
+++ src/FunctionApproximation/Extractors.cpp	(revision 1155ba032fdb209f9d0641f2d81b1d0f593c0c05)
@@ -531,6 +531,6 @@
 }
 
-FunctionModel::list_of_arguments_t Extractors::reorderArgumentsByParticleTypes(
-    const FunctionModel::list_of_arguments_t &listargs,
+FunctionModel::list_of_arguments_t Extractors::filterArgumentsByBindingModel(
+    const FunctionModel::arguments_t &args,
     const HomologyGraph &_graph,
     const ParticleTypes_t &_types,
@@ -539,221 +539,217 @@
 {
   FunctionModel::list_of_arguments_t returnargs;
-  for (FunctionModel::list_of_arguments_t::const_iterator listiter = listargs.begin();
-      listiter != listargs.end(); ++listiter) {
-    const FunctionModel::arguments_t &args = *listiter;
-
-    // deal with the case when there are no distances (ConstantPotential)
-    if (_bindingmodel.getNodes().size() < 2) {
-      LOG(3, "DEBUG: Potential requires only one or no particle types, needs no distances.");
-      continue;
-    }
-    if (_bindingmodel.getGraph().getEdges().empty()) {
-      LOG(3, "DEBUG: Potential represents non-bonded interactions, gets all distances.");
-      // TODO: Here we need to constrain to all distances matching the types?
-      returnargs.push_back(args);
-      continue;
-    }
-
-    /// 0. place all nodes in a lookup map for their type
-    type_index_lookup_t type_index_lookup;
-    for(FunctionModel::arguments_t::const_iterator iter = args.begin();
-        iter != args.end(); ++iter) {
-      if (type_index_lookup.left.find(iter->indices.first) == type_index_lookup.left.end())
-        type_index_lookup.left.insert( std::make_pair(iter->indices.first, iter->types.first) );
-      else {
-        ASSERT(type_index_lookup.left.at(iter->indices.first) == iter->types.first,
-            "Extractors::reorderArgumentsByParticleTypes() - entry " +toString(iter->indices.first)+
-            " is already present with different type "+toString(iter->types.first)
-            +" than "+toString(type_index_lookup.left.at(iter->indices.first)));
-      }
-      if (type_index_lookup.left.find(iter->indices.second) == type_index_lookup.left.end())
-        type_index_lookup.left.insert( std::make_pair(iter->indices.second, iter->types.second) );
-      else {
-        ASSERT(type_index_lookup.left.at(iter->indices.second) == iter->types.second,
-            "Extractors::reorderArgumentsByParticleTypes() - entry " +toString(iter->indices.second)+
-            " is already present with different type "+toString(iter->types.second)
-            +" than "+toString(type_index_lookup.left.at(iter->indices.second)));
-      }
-    }
-    if (DoLog(3)) {
-      std::stringstream output;
-      for (type_index_lookup_t::left_const_iterator indexiter = type_index_lookup.left.begin();
-          indexiter != type_index_lookup.left.end(); ++indexiter) {
-        output << " [" << indexiter->first << "," << indexiter->second << "]";
-      }
-      LOG(3, "DEBUG: index to type map:" << output.str());
-    }
-    if (DoLog(3)) {
-      std::stringstream output;
-      for (type_index_lookup_t::right_const_iterator typeiter = type_index_lookup.right.begin();
-          typeiter != type_index_lookup.right.end(); ++typeiter) {
-        output << " [" << typeiter->first << "," << typeiter->second << "]";
-      }
-      LOG(3, "DEBUG: type to index map " << output.str());
-    }
-
-    /// 1. construct boost::graph from arguments_t (iter)
-    const size_t N = type_index_lookup.left.size();
-    UndirectedGraph fragmentgraph(N);
-    for(FunctionModel::arguments_t::const_iterator iter = args.begin();
-        iter != args.end(); ++iter) {
-      if (iter->bonded)
-        boost::add_edge(iter->indices.first, iter->indices.second, fragmentgraph);
-    }
-    index_map_t index_map = boost::get(boost::vertex_index, fragmentgraph);
-    LOG(2, "DEBUG: We have " << boost::num_vertices(fragmentgraph) << " nodes in the fragment graph.");
-
-    /// 2. grab first node of the binding model (gives a type)
-    const BindingModel::vector_nodes_t::const_iterator firstiter = _bindingmodel.getNodes().begin();
-    const FragmentNode &firstnode = *firstiter;
-    const Extractors::ParticleType_t &firsttype = firstnode.getAtomicNumber();
-
-    /// 3. grab all candidate nodes contained in arguments_t
-    std::pair<
-        type_index_lookup_t::right_const_iterator,
-        type_index_lookup_t::right_const_iterator> range = type_index_lookup.right.equal_range(firsttype);
-
-    /// 4. go over all candidate nodes (gives index)
-    const size_t nodes_in_bindingmodel = _bindingmodel.getNodes().size();
-    LOG(2, "DEBUG: There are " << nodes_in_bindingmodel << " nodes in the binding model.");
-    set_of_nodes_t set_of_nodes;
-    for (type_index_lookup_t::right_const_iterator rangeiter = range.first;
-        rangeiter != range.second; ++rangeiter) {
-      const size_t &rootindex = rangeiter->second;
-      LOG(2, "DEBUG: Current root index is " << rootindex);
-
-      /// 5a. from node in graph (with this index) perform BFS till n-1 (#nodes in binding model)
-      std::vector<size_t> distances(boost::num_vertices(fragmentgraph));
-      boost::breadth_first_search(
-          fragmentgraph,
-          boost::vertex(rootindex, fragmentgraph),
-          boost::visitor(record_distance(&distances[0])));
-      LOG(3, "DEBUG: BFS discovered the following distances " << distances);
-
-      /// 5b. and store all node in map with distance to root as key
-      nodes_per_level_t nodes_per_level;
-      for (size_t i=0;i<distances.size();++i) {
-        nodes_per_level.insert( std::make_pair(level_t(distances[i]), node_t(i)) );
-      }
-      LOG(3, "DEBUG: Nodes per level are " << nodes_per_level);
-
-      /// 6. construct all possible induced connected subgraphs with this map (see fragmentation)
-      nodes_t nodes;
-      // rootindex is always contained
-      nodes.insert( rootindex );
-      level_t level = 0;
-
-      generateAllInducedConnectedSubgraphs(
-          nodes_in_bindingmodel-1, level, nodes, set_of_nodes, nodes_per_level, fragmentgraph, distances, index_map);
-      LOG(2, "DEBUG: We have found " << set_of_nodes.size() << " unique induced, connected subgraphs.");
-    }
-
-    /// 8. go through each induced, connected subgraph
-    for (set_of_nodes_t::const_iterator nodesiter = set_of_nodes.begin();
-        nodesiter != set_of_nodes.end(); ++nodesiter) {
-      /// 9. convert its nodes into a HomologyGraph using the subgraph and arguments_t (types for index)
-      const nodes_t &nodes = *nodesiter;
-      const HomologyGraph nodes_graph =
-          createHomologyGraphFromNodes(nodes, type_index_lookup, fragmentgraph, index_map);
-
-      /// 10. compare each converted HomologyGraph with _bindingmodel: Accept or Reject
-      if (nodes_graph == _bindingmodel.getGraph()) {
-        LOG(2, "ACCEPT: " << nodes_graph << " is identical to " << _bindingmodel.getGraph());
-        /// 11. for each accepted keyset, pick _all_ symmetric distances from arguments_t
-        FunctionModel::arguments_t argumentbunch;
-        argumentbunch.reserve(args.size());
-        for(FunctionModel::arguments_t::const_iterator iter = args.begin();
-            iter != args.end(); ++iter) {
-          if ((nodes.count(iter->indices.first) != 0) && (nodes.count(iter->indices.second) != 0)) {
-            argumentbunch.push_back(*iter);
+
+  // deal with the case when there are no distances (ConstantPotential)
+  if (_bindingmodel.getNodes().size() < 2) {
+    LOG(3, "DEBUG: Potential requires only one or no particle types, needs no distances.");
+    return returnargs;
+  }
+  if (_bindingmodel.getGraph().getEdges().empty()) {
+    LOG(3, "DEBUG: Potential represents non-bonded interactions, gets all distances.");
+    // TODO: Here we need to constrain to all distances matching the types?
+    returnargs.push_back(args);
+    return returnargs;
+  }
+
+  /// 0. place all nodes in a lookup map for their type
+  type_index_lookup_t type_index_lookup;
+  for(FunctionModel::arguments_t::const_iterator iter = args.begin();
+      iter != args.end(); ++iter) {
+    if (type_index_lookup.left.find(iter->indices.first) == type_index_lookup.left.end())
+      type_index_lookup.left.insert( std::make_pair(iter->indices.first, iter->types.first) );
+    else {
+      ASSERT(type_index_lookup.left.at(iter->indices.first) == iter->types.first,
+          "Extractors::reorderArgumentsByParticleTypes() - entry " +toString(iter->indices.first)+
+          " is already present with different type "+toString(iter->types.first)
+          +" than "+toString(type_index_lookup.left.at(iter->indices.first)));
+    }
+    if (type_index_lookup.left.find(iter->indices.second) == type_index_lookup.left.end())
+      type_index_lookup.left.insert( std::make_pair(iter->indices.second, iter->types.second) );
+    else {
+      ASSERT(type_index_lookup.left.at(iter->indices.second) == iter->types.second,
+          "Extractors::reorderArgumentsByParticleTypes() - entry " +toString(iter->indices.second)+
+          " is already present with different type "+toString(iter->types.second)
+          +" than "+toString(type_index_lookup.left.at(iter->indices.second)));
+    }
+  }
+  if (DoLog(3)) {
+    std::stringstream output;
+    for (type_index_lookup_t::left_const_iterator indexiter = type_index_lookup.left.begin();
+        indexiter != type_index_lookup.left.end(); ++indexiter) {
+      output << " [" << indexiter->first << "," << indexiter->second << "]";
+    }
+    LOG(3, "DEBUG: index to type map:" << output.str());
+  }
+  if (DoLog(3)) {
+    std::stringstream output;
+    for (type_index_lookup_t::right_const_iterator typeiter = type_index_lookup.right.begin();
+        typeiter != type_index_lookup.right.end(); ++typeiter) {
+      output << " [" << typeiter->first << "," << typeiter->second << "]";
+    }
+    LOG(3, "DEBUG: type to index map " << output.str());
+  }
+
+  /// 1. construct boost::graph from arguments_t (iter)
+  const size_t N = type_index_lookup.left.size();
+  UndirectedGraph fragmentgraph(N);
+  for(FunctionModel::arguments_t::const_iterator iter = args.begin();
+      iter != args.end(); ++iter) {
+    if (iter->bonded)
+      boost::add_edge(iter->indices.first, iter->indices.second, fragmentgraph);
+  }
+  index_map_t index_map = boost::get(boost::vertex_index, fragmentgraph);
+  LOG(2, "DEBUG: We have " << boost::num_vertices(fragmentgraph) << " nodes in the fragment graph.");
+
+  /// 2. grab first node of the binding model (gives a type)
+  const BindingModel::vector_nodes_t::const_iterator firstiter = _bindingmodel.getNodes().begin();
+  const FragmentNode &firstnode = *firstiter;
+  const Extractors::ParticleType_t &firsttype = firstnode.getAtomicNumber();
+
+  /// 3. grab all candidate nodes contained in arguments_t
+  std::pair<
+      type_index_lookup_t::right_const_iterator,
+      type_index_lookup_t::right_const_iterator> range = type_index_lookup.right.equal_range(firsttype);
+
+  /// 4. go over all candidate nodes (gives index)
+  const size_t nodes_in_bindingmodel = _bindingmodel.getNodes().size();
+  LOG(2, "DEBUG: There are " << nodes_in_bindingmodel << " nodes in the binding model.");
+  set_of_nodes_t set_of_nodes;
+  for (type_index_lookup_t::right_const_iterator rangeiter = range.first;
+      rangeiter != range.second; ++rangeiter) {
+    const size_t &rootindex = rangeiter->second;
+    LOG(2, "DEBUG: Current root index is " << rootindex);
+
+    /// 5a. from node in graph (with this index) perform BFS till n-1 (#nodes in binding model)
+    std::vector<size_t> distances(boost::num_vertices(fragmentgraph));
+    boost::breadth_first_search(
+        fragmentgraph,
+        boost::vertex(rootindex, fragmentgraph),
+        boost::visitor(record_distance(&distances[0])));
+    LOG(3, "DEBUG: BFS discovered the following distances " << distances);
+
+    /// 5b. and store all node in map with distance to root as key
+    nodes_per_level_t nodes_per_level;
+    for (size_t i=0;i<distances.size();++i) {
+      nodes_per_level.insert( std::make_pair(level_t(distances[i]), node_t(i)) );
+    }
+    LOG(3, "DEBUG: Nodes per level are " << nodes_per_level);
+
+    /// 6. construct all possible induced connected subgraphs with this map (see fragmentation)
+    nodes_t nodes;
+    // rootindex is always contained
+    nodes.insert( rootindex );
+    level_t level = 0;
+
+    generateAllInducedConnectedSubgraphs(
+        nodes_in_bindingmodel-1, level, nodes, set_of_nodes, nodes_per_level, fragmentgraph, distances, index_map);
+    LOG(2, "DEBUG: We have found " << set_of_nodes.size() << " unique induced, connected subgraphs.");
+  }
+
+  /// 8. go through each induced, connected subgraph
+  for (set_of_nodes_t::const_iterator nodesiter = set_of_nodes.begin();
+      nodesiter != set_of_nodes.end(); ++nodesiter) {
+    /// 9. convert its nodes into a HomologyGraph using the subgraph and arguments_t (types for index)
+    const nodes_t &nodes = *nodesiter;
+    const HomologyGraph nodes_graph =
+        createHomologyGraphFromNodes(nodes, type_index_lookup, fragmentgraph, index_map);
+
+    /// 10. compare each converted HomologyGraph with _bindingmodel: Accept or Reject
+    if (nodes_graph == _bindingmodel.getGraph()) {
+      LOG(2, "ACCEPT: " << nodes_graph << " is identical to " << _bindingmodel.getGraph());
+      /// 11. for each accepted keyset, pick _all_ symmetric distances from arguments_t
+      FunctionModel::arguments_t argumentbunch;
+      argumentbunch.reserve(args.size());
+      for(FunctionModel::arguments_t::const_iterator iter = args.begin();
+          iter != args.end(); ++iter) {
+        if ((nodes.count(iter->indices.first) != 0) && (nodes.count(iter->indices.second) != 0)) {
+          argumentbunch.push_back(*iter);
+        }
+      }
+      const size_t num_symmetric_distances = nodes.size()*(nodes.size()-1)/2;
+      ASSERT( argumentbunch.size() == num_symmetric_distances,
+          "Extractors::reorderArgumentsByParticleTypes() - only "+toString(argumentbunch.size())+
+          " found instead of "+toString(num_symmetric_distances));
+      LOG(3, "DEBUG: Final bunch of arguments is " << argumentbunch << ".");
+
+      /**
+       * We still need to bring the arguments in the correct order for the potential.
+       *
+       * The potential gives the desired order via the nodes in the HomologyGraph! Their
+       * sequence matches with the user-defined particle types and because of the properties
+       * of the homology graph (FragmentNode stores type and number of edges) it also
+       * matches with the desired bonding graph.
+       */
+
+      /// So, we need to extract from the argumentbunch the information contained in each
+      /// FragmentNode, namely its type and the number of connected edges
+      node_FragmentNode_map_t node_FragmentNode_map = fillNodeFragmentMap(argumentbunch);
+
+      /// Then, we step through the nodes of the bindingmodel ...
+      /// ... and find a suitable mapping from indices in the arguments to
+      /// the index in the order of the HomologyGraph's nodes
+      const BindingModel::vector_nodes_t bindingmodel_nodes = _bindingmodel.getNodes();
+      argindex_to_nodeindex_t argindex_to_nodeindex;
+      size_t nodeindex = 0;
+      for (BindingModel::vector_nodes_t::const_iterator nodeiter = bindingmodel_nodes.begin();
+          nodeiter != bindingmodel_nodes.end(); ++nodeiter) {
+        const FragmentNode &node = *nodeiter;
+        LOG(3, "DEBUG: Binding model's node #" << node << ".");
+        /// ... and pick for each the first (and unique) from these stored nodes.
+        node_FragmentNode_map_t::iterator mapiter = node_FragmentNode_map.begin();
+        for (;mapiter != node_FragmentNode_map.end(); ++mapiter) {
+          if ((mapiter->second.first == node.getAtomicNumber())
+              && (mapiter->second.second == node.getConnectedEdges())) {
+            LOG(3, "DEBUG: #" << mapiter->first << " with type " << mapiter->second.first
+                << " and " << mapiter->second.second << " connected edges matches.");
+            break;
           }
         }
-        const size_t num_symmetric_distances = nodes.size()*(nodes.size()-1)/2;
-        ASSERT( argumentbunch.size() == num_symmetric_distances,
-            "Extractors::reorderArgumentsByParticleTypes() - only "+toString(argumentbunch.size())+
-            " found instead of "+toString(num_symmetric_distances));
-        LOG(3, "DEBUG: Final bunch of arguments is " << argumentbunch << ".");
-
-        /**
-         * We still need to bring the arguments in the correct order for the potential.
-         *
-         * The potential gives the desired order via the nodes in the HomologyGraph! Their
-         * sequence matches with the user-defined particle types and because of the properties
-         * of the homology graph (FragmentNode stores type and number of edges) it also
-         * matches with the desired bonding graph.
-         */
-
-        /// So, we need to extract from the argumentbunch the information contained in each
-        /// FragmentNode, namely its type and the number of connected edges
-        node_FragmentNode_map_t node_FragmentNode_map = fillNodeFragmentMap(argumentbunch);
-
-        /// Then, we step through the nodes of the bindingmodel ...
-        /// ... and find a suitable mapping from indices in the arguments to
-        /// the index in the order of the HomologyGraph's nodes
-        const BindingModel::vector_nodes_t bindingmodel_nodes = _bindingmodel.getNodes();
-        argindex_to_nodeindex_t argindex_to_nodeindex;
-        size_t nodeindex = 0;
-        for (BindingModel::vector_nodes_t::const_iterator nodeiter = bindingmodel_nodes.begin();
-            nodeiter != bindingmodel_nodes.end(); ++nodeiter) {
-          const FragmentNode &node = *nodeiter;
-          LOG(3, "DEBUG: Binding model's node #" << node << ".");
-          /// ... and pick for each the first (and unique) from these stored nodes.
-          node_FragmentNode_map_t::iterator mapiter = node_FragmentNode_map.begin();
-          for (;mapiter != node_FragmentNode_map.end(); ++mapiter) {
-            if ((mapiter->second.first == node.getAtomicNumber())
-                && (mapiter->second.second == node.getConnectedEdges())) {
-              LOG(3, "DEBUG: #" << mapiter->first << " with type " << mapiter->second.first
-                  << " and " << mapiter->second.second << " connected edges matches.");
-              break;
-            }
-          }
-          ASSERT( mapiter != node_FragmentNode_map.end(),
-              "Extractors::reorderArgumentsByParticleTypes() - could not find a suitable node for #"+
-              toString(mapiter->first)+" with type "+toString(mapiter->second.first)+" and "+
-              toString(mapiter->second.second)+" connected edges");
-          std::pair<argindex_to_nodeindex_t::iterator, bool> inserter =
-              argindex_to_nodeindex.insert( std::make_pair(mapiter->first, nodeindex++) );
-          ASSERT( inserter.second,
-              "Extractors::reorderArgumentsByParticleTypes() - node #"+toString(mapiter->first)+
-              " is already present?");
-          // remove to ensure uniqueness of choice
-          node_FragmentNode_map.erase(mapiter);
-        }
-        LOG(4, "DEBUG: argument's indices to node index map is " << argindex_to_nodeindex);
-        // i.e. this is not the arg's index in argumentbunch, but the index of the position
-        // contained in the argument_t
-
-        /// Finally, we only need to bring the arguments in the typical order:
-        /// 01 02 03 04 ... 0n, 12 13 14 ... 1n, 23 24 25 ... 2n, ...,  (n-1)n
-        /// These ordering we store in a map for each argument's indices
-        const size_t num_args = argindex_to_nodeindex.size();
-        argument_placement_map_t argument_placement_map = fillArgumentsPlacementMap(num_args);
-        LOG(4, "DEBUG: Placement map is " << argument_placement_map);
-        ASSERT( argument_placement_map.size() == argumentbunch.size(),
-            "Extractors::reorderArgumentsByParticleTypes() - placement map has size "+
-            toString(argument_placement_map.size())+" and we expected "+toString(argumentbunch.size()));
-
-        // and finally resort the arguments with the constructed placement map
-        FunctionModel::arguments_t sortedargs(argumentbunch.size());
-        for (FunctionModel::arguments_t::const_iterator argiter = argumentbunch.begin();
-            argiter != argumentbunch.end(); ++argiter) {
-          const argument_t &arg = *argiter;
-          const argument_t::indices_t translatedIndices =
-              translateIndices(argindex_to_nodeindex, arg.indices);
-          const argument_placement_map_t::const_iterator indexiter =
-              argument_placement_map.find( translatedIndices );
-          ASSERT( indexiter != argument_placement_map.end(),
-              "Extractors::reorderArgumentsByParticleTypes() - could not find place for edge "+
-              toString(arg.indices));
-          sortedargs[indexiter->second] = arg;
-          LOG(3, "DEBUG: Placed argument " << arg << " at #" << indexiter->second
-              << " as translated indices are " << translatedIndices);
-        }
-        LOG(2, "DEBUG: Sorted arguments are " << sortedargs << ".");
-
-        returnargs.push_back(sortedargs);
-      } else {
-        LOG(2, "REJECT: " << nodes_graph << " is not identical to " << _bindingmodel.getGraph());
-      }
+        ASSERT( mapiter != node_FragmentNode_map.end(),
+            "Extractors::reorderArgumentsByParticleTypes() - could not find a suitable node for #"+
+            toString(mapiter->first)+" with type "+toString(mapiter->second.first)+" and "+
+            toString(mapiter->second.second)+" connected edges");
+        std::pair<argindex_to_nodeindex_t::iterator, bool> inserter =
+            argindex_to_nodeindex.insert( std::make_pair(mapiter->first, nodeindex++) );
+        ASSERT( inserter.second,
+            "Extractors::reorderArgumentsByParticleTypes() - node #"+toString(mapiter->first)+
+            " is already present?");
+        // remove to ensure uniqueness of choice
+        node_FragmentNode_map.erase(mapiter);
+      }
+      LOG(4, "DEBUG: argument's indices to node index map is " << argindex_to_nodeindex);
+      // i.e. this is not the arg's index in argumentbunch, but the index of the position
+      // contained in the argument_t
+
+      /// Finally, we only need to bring the arguments in the typical order:
+      /// 01 02 03 04 ... 0n, 12 13 14 ... 1n, 23 24 25 ... 2n, ...,  (n-1)n
+      /// These ordering we store in a map for each argument's indices
+      const size_t num_args = argindex_to_nodeindex.size();
+      argument_placement_map_t argument_placement_map = fillArgumentsPlacementMap(num_args);
+      LOG(4, "DEBUG: Placement map is " << argument_placement_map);
+      ASSERT( argument_placement_map.size() == argumentbunch.size(),
+          "Extractors::reorderArgumentsByParticleTypes() - placement map has size "+
+          toString(argument_placement_map.size())+" and we expected "+toString(argumentbunch.size()));
+
+      // and finally resort the arguments with the constructed placement map
+      FunctionModel::arguments_t sortedargs(argumentbunch.size());
+      for (FunctionModel::arguments_t::const_iterator argiter = argumentbunch.begin();
+          argiter != argumentbunch.end(); ++argiter) {
+        const argument_t &arg = *argiter;
+        const argument_t::indices_t translatedIndices =
+            translateIndices(argindex_to_nodeindex, arg.indices);
+        const argument_placement_map_t::const_iterator indexiter =
+            argument_placement_map.find( translatedIndices );
+        ASSERT( indexiter != argument_placement_map.end(),
+            "Extractors::reorderArgumentsByParticleTypes() - could not find place for edge "+
+            toString(arg.indices));
+        sortedargs[indexiter->second] = arg;
+        LOG(3, "DEBUG: Placed argument " << arg << " at #" << indexiter->second
+            << " as translated indices are " << translatedIndices);
+      }
+      LOG(2, "DEBUG: Sorted arguments are " << sortedargs << ".");
+
+      returnargs.push_back(sortedargs);
+    } else {
+      LOG(2, "REJECT: " << nodes_graph << " is not identical to " << _bindingmodel.getGraph());
     }
   }
@@ -769,8 +765,8 @@
     )
 {
+  // list allows for quicker removal than vector
   typedef std::list< argument_t > ListArguments_t;
   ListArguments_t availableList(args.begin(), args.end());
   LOG(2, "DEBUG: Initial list of args is " << args << ".");
-
 
   // TODO: fill a lookup map such that we don't have O(M^3) scaling, if M is number
@@ -779,12 +775,15 @@
 //  typedef ParticleTypes_t firsttype;
 //  typedef ParticleTypes_t secondtype;
-//  typedef std::map< firsttype, std::map< secondtype, boost::ref(args) > > ArgsLookup_t;
+//  typedef std::map< firsttype, std::map< secondtype, FunctionModel::arguments_t > > ArgsLookup_t;
 //  ArgsLookup_t ArgsLookup;
+
+  ASSERT( _types.size() <= 2,
+      "Extractors::filterArgumentsByParticleTypes() - this only filters and does not reorder."
+      +std::string(" Hence, it is not useful for multiple arguments per model."));
 
   // basically, we have two choose any two pairs out of types but only those
   // 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.
-  FunctionModel::arguments_t allargs;
-  allargs.reserve(args.size());
+  FunctionModel::list_of_arguments_t returnargs;
   for (ParticleTypes_t::const_iterator firstiter = _types.begin();
       firstiter != _types.end();
@@ -804,10 +803,10 @@
         if ((iter->types.first == *firstiter)
               && (iter->types.second == *seconditer)) {
-          allargs.push_back( *iter );
+          returnargs.push_back( FunctionModel::arguments_t(1, *iter) );
           iter = availableList.erase(iter);
           LOG(4, "DEBUG: Accepted argument.");
         } else if ((iter->types.first == *seconditer)
               && (iter->types.second == *firstiter)) {
-          allargs.push_back( *iter );
+          returnargs.push_back( FunctionModel::arguments_t(1, *iter) );
           iter = availableList.erase(iter);
           LOG(4, "DEBUG: Accepted (flipped) argument.");
@@ -819,30 +818,4 @@
     }
   }
-  LOG(2, "DEBUG: Final list of args is " << allargs << ".");
-
-  // first, we bring together tuples of distances that belong together
-  FunctionModel::list_of_arguments_t singlelist_allargs;
-  singlelist_allargs.push_back(allargs);
-  FunctionModel::list_of_arguments_t returnargs =
-      reorderArgumentsByParticleTypes(singlelist_allargs, _graph, _types, _bindingmodel);
-//  ASSERT( sortedargs.size() == (size_t)1,
-//      "Extractors::filterArgumentsByParticleTypes() - reordering did not generate a single list.");
-//  // then we split up the tuples of arguments and place each into single list
-//  FunctionModel::list_of_arguments_t returnargs;
-//  FunctionModel::arguments_t::const_iterator argiter = sortedargs.begin()->begin();
-//  const size_t num_types = _types.size();
-//  const size_t args_per_tuple = num_types * (num_types-1) / 2;
-//  while (argiter != sortedargs.begin()->end()) {
-//    FunctionModel::arguments_t currenttuple(args_per_tuple);
-//    const FunctionModel::arguments_t::const_iterator startiter = argiter;
-//    std::advance(argiter, args_per_tuple);
-//#ifndef NDEBUG
-//    FunctionModel::arguments_t::const_iterator endoutiter =
-//#endif
-//        std::copy(startiter, argiter, currenttuple.begin());
-//    ASSERT( endoutiter == currenttuple.end(),
-//        "Extractors::filterArgumentsByParticleTypes() - currenttuple not initialized to right size.");
-//    returnargs.push_back(currenttuple);
-//  }
 
   LOG(2, "DEBUG: We have generated " << returnargs.size() << " tuples of distances.");
Index: src/FunctionApproximation/Extractors.hpp
===================================================================
--- src/FunctionApproximation/Extractors.hpp	(revision f5dbea01dfa29040cfe86148e9597b8da56187fd)
+++ src/FunctionApproximation/Extractors.hpp	(revision 1155ba032fdb209f9d0641f2d81b1d0f593c0c05)
@@ -103,10 +103,8 @@
   }
 
-  /** Reorder the arguments to bring adjacent ones together.
+  /** Filter the arguments to select only these required by the model.
    *
-   * 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()
+   * \warning this is meant as a faster way of getting the arguments for simple
+   * pair potentials. In any other case, one should use filterArgumentsByBindingModel()
    *
    * \param listargs list of arguments to reorder each
@@ -115,6 +113,6 @@
    * \return reordered args
    */
-  FunctionModel::list_of_arguments_t reorderArgumentsByParticleTypes(
-      const FunctionModel::list_of_arguments_t &eachargs,
+  FunctionModel::list_of_arguments_t filterArgumentsByParticleTypes(
+      const FunctionModel::arguments_t &args,
       const HomologyGraph &_graph,
       const ParticleTypes_t &_types,
@@ -122,15 +120,17 @@
       );
 
-  /** Filter arguments according to types, allowing multiples.
+  /** Filter and reorder the arguments to bring adjacent ones together.
    *
-   * 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).
+   * We need to find all matching subgraphs (given by \a _bindingmodel) in the
+   * given homology graph (given by \a _graph) of the fragment molecule.
+   * This means filtering down to the desired particle types and then find
+   * all possible matching subgraphs in each of argument lists, \a eachargs.
    *
-   * \param args arguments to reorder
+   * \param listargs list of arguments to filter and order appropriately
    * \param _graph contains binding model of graph
    * \param _types particle type vector
-   * \return filtered list of args
+   * \return reordered args
    */
-  FunctionModel::list_of_arguments_t filterArgumentsByParticleTypes(
+  FunctionModel::list_of_arguments_t filterArgumentsByBindingModel(
       const FunctionModel::arguments_t &args,
       const HomologyGraph &_graph,
Index: src/Potentials/Specifics/FourBodyPotential_Torsion.cpp
===================================================================
--- src/Potentials/Specifics/FourBodyPotential_Torsion.cpp	(revision f5dbea01dfa29040cfe86148e9597b8da56187fd)
+++ src/Potentials/Specifics/FourBodyPotential_Torsion.cpp	(revision 1155ba032fdb209f9d0641f2d81b1d0f593c0c05)
@@ -275,9 +275,6 @@
 {
   FunctionModel::filter_t returnfunction =
-      boost::bind(&Extractors::reorderArgumentsByParticleTypes,
-        boost::bind(&Extractors::filterArgumentsByParticleTypes,
-            _2, _1,
-            boost::cref(getParticleTypes()), boost::cref(getBindingModel())),
-        _1,
+      boost::bind(&Extractors::filterArgumentsByBindingModel,
+        _2, _1,
         boost::cref(getParticleTypes()), boost::cref(getBindingModel())
       );
Index: src/Potentials/Specifics/ThreeBodyPotential_Angle.cpp
===================================================================
--- src/Potentials/Specifics/ThreeBodyPotential_Angle.cpp	(revision f5dbea01dfa29040cfe86148e9597b8da56187fd)
+++ src/Potentials/Specifics/ThreeBodyPotential_Angle.cpp	(revision 1155ba032fdb209f9d0641f2d81b1d0f593c0c05)
@@ -257,9 +257,6 @@
 {
   FunctionModel::filter_t returnfunction =
-      boost::bind(&Extractors::reorderArgumentsByParticleTypes,
-        boost::bind(&Extractors::filterArgumentsByParticleTypes,
-            _2, _1,
-          boost::cref(getParticleTypes()), boost::cref(getBindingModel())),
-        _1,
+      boost::bind(&Extractors::filterArgumentsByBindingModel,
+        _2, _1,
         boost::cref(getParticleTypes()), boost::cref(getBindingModel())
       );
Index: tests/regression/Potential/FitPotential/testsuite-potential-fit-potential.at
===================================================================
--- tests/regression/Potential/FitPotential/testsuite-potential-fit-potential.at	(revision f5dbea01dfa29040cfe86148e9597b8da56187fd)
+++ tests/regression/Potential/FitPotential/testsuite-potential-fit-potential.at	(revision 1155ba032fdb209f9d0641f2d81b1d0f593c0c05)
@@ -22,5 +22,4 @@
 AT_KEYWORDS([potential parse-homologies fit-potential morse])
 AT_SKIP_IF([../../molecuilder --help fit-potential; if test $? -eq 5; then /bin/true; else /bin/false; fi])
-AT_XFAIL_IF([/bin/true])
 # homology files probably originate from Labspace/PotentialFitting/WaterAngle
 # recreate homology file: PYTHONPATH=... /usr/bin/python .../tests/regression/Potential/FitPotential/pre/measurement_function.py length_common.pdb 2 0 0 5 4 127.0.0.1 2050 6-311G* 1 2. 10 .../share/molecuilder/data/bondtables/bondtable.dat
@@ -51,5 +50,4 @@
 AT_KEYWORDS([potential parse-homologies fit-potential harmonic])
 AT_SKIP_IF([../../molecuilder --help fit-potential; if test $? -eq 5; then /bin/true; else /bin/false; fi])
-AT_XFAIL_IF([/bin/true])
 # homology files probably originate from Labspace/PotentialFitting/Water
 # see morse fitting
