Index: src/Potentials/CompoundPotential.cpp
===================================================================
--- src/Potentials/CompoundPotential.cpp	(revision c5e63a4efcf5739f8379980891648d24490fe8c4)
+++ src/Potentials/CompoundPotential.cpp	(revision 7c10913cd7643f0bb610d97273b66628c8faa10d)
@@ -84,5 +84,5 @@
     Extractors::elementcounts_t::const_iterator countiter = counts_per_charge.begin();
     for (; countiter != counts_per_charge.end(); ++countiter)
-      if (!graph.hasTimesAtomicNumber(
+      if (!graph.hasGreaterEqualTimesAtomicNumber(
           static_cast<size_t>(countiter->first),
           static_cast<size_t>(countiter->second))
@@ -115,6 +115,5 @@
   parameters_t::const_iterator iter = _params.begin();
   BOOST_FOREACH( FunctionModel* model, models) {
-    const parameters_t &model_params = model->getParameters();
-    const size_t model_dim = model_params.size();
+    const size_t model_dim = model->getParameterDimension();
     if (dim > 0) {
       parameters_t subparams;
@@ -131,4 +130,12 @@
     }
   }
+
+#ifndef NDEBUG
+  parameters_t check_params(getParameters());
+  check_params.resize(_params.size()); // truncate to same size
+  ASSERT( check_params == _params,
+      "CompoundPotential::setParameters() - failed, mismatch in to be set "
+      +toString(_params)+" and set "+toString(check_params)+" params.");
+#endif
 }
 
@@ -140,8 +147,10 @@
   BOOST_FOREACH( const FunctionModel* model, models) {
     const CompoundPotential::parameters_t &params = model->getParameters();
-    std::copy(params.begin(), params.end(), iter);
     ASSERT( iter != parameters.end(),
         "CompoundPotential::getParameters() - iter already at end.");
-  }
+    iter = std::copy(params.begin(), params.end(), iter);
+  }
+  ASSERT( iter == parameters.end(),
+      "CompoundPotential::getParameters() - iter not at end.");
   return parameters;
 }
@@ -169,20 +178,122 @@
 }
 
-std::vector<CompoundPotential::arguments_t> CompoundPotential::splitUpArgumentsByModels(
+bool CompoundPotential::areValidArguments(
+    const SerializablePotential::ParticleTypes_t &_types,
+    const arguments_t &args) const
+{
+  // /this function does much the same as Extractors::reorderArgumentsByParticleTypes()
+  typedef std::list< argument_t > ListArguments_t;
+  ListArguments_t availableList(args.begin(), args.end());
+
+  /// 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
+  /// iterator at the current position of the first one and skip the equal case.
+  for (SerializablePotential::ParticleTypes_t::const_iterator firstiter = _types.begin();
+      firstiter != _types.end();
+      ++firstiter) {
+    for (SerializablePotential::ParticleTypes_t::const_iterator seconditer = firstiter;
+        seconditer != _types.end();
+        ++seconditer) {
+      if (seconditer == firstiter)
+        continue;
+
+      // search the right one in _args (we might allow switching places of
+      // firstiter and seconditer, as distance is symmetric).
+      // we remove the matching argument to make sure we don't pick it twice
+      ListArguments_t::iterator iter = availableList.begin();
+      for (;iter != availableList.end(); ++iter) {
+        LOG(3, "DEBUG: Current args is " << *iter << ".");
+        if ((iter->types.first == *firstiter)
+              && (iter->types.second == *seconditer)) {
+          availableList.erase(iter);
+          break;
+        }
+        else if ((iter->types.first == *seconditer)
+              && (iter->types.second == *firstiter)) {
+          availableList.erase(iter);
+          break;
+        }
+      }
+      if ( iter == availableList.end())
+        return false;
+    }
+  }
+
+  return true;
+}
+
+CompoundPotential::arguments_by_model_t CompoundPotential::splitUpArgumentsByModels(
     const arguments_t &arguments) const
 {
-  std::vector<arguments_t> partial_args(models.size());
-  std::vector<arguments_t>::iterator partialiter = partial_args.begin();
+  arguments_by_model_t partial_args;
   arguments_t::const_iterator argiter = arguments.begin();
-  BOOST_FOREACH( const SerializablePotential::ParticleTypes_t &types, particletypes_per_model) {
-    // we always expect N(N-1)/2 distances for N particle types
-    arguments_t::const_iterator enditer = argiter+(types.size()*(types.size()-1)/2);
-    ASSERT( argiter != arguments.end(),
-        "CompoundPotential::splitUpArgumentsByModels() - incorrect number of arguments.");
-    std::copy(argiter, enditer, std::back_inserter(*partialiter++));
-    argiter = enditer;
-  }
-  ASSERT( argiter == arguments.end(),
-      "CompoundPotential::splitUpArgumentsByModels() - incorrect number of arguments.");
+  particletypes_per_model_t::const_iterator typesiter = particletypes_per_model.begin();
+  models_t::const_iterator modeliter = models.begin();
+
+  // add constant model (which is always first model) with empty args if present
+  if (typesiter->empty()) {
+    partial_args.push_back(
+        std::pair<FunctionModel *, arguments_t>(*modeliter, arguments_t())
+        );
+    ++modeliter;
+    ++typesiter;
+  }
+  // then check other models
+  while (argiter != arguments.end()) {
+    if (typesiter+1 != particletypes_per_model.end()) {
+      // check whether next argument bunch is for same model or different one
+      // we extract both partial_arguments, if the latter fits, we take the latter.
+      const SerializablePotential::ParticleTypes_t &types = *typesiter;
+      const SerializablePotential::ParticleTypes_t &nexttypes = *(typesiter+1);
+
+      // we always expect N(N-1)/2 distances for N particle types
+      arguments_t::const_iterator enditer = argiter+(types.size()*(types.size()-1)/2);
+      arguments_t::const_iterator nextenditer = argiter+(nexttypes.size()*(nexttypes.size()-1)/2);
+      arguments_t args(argiter, enditer);
+      arguments_t nextargs(argiter, nextenditer);
+      if (types.size() < nexttypes.size()) {
+        if (areValidArguments(nexttypes, nextargs)) {
+          // latter matches, increment
+          ++typesiter;
+          partial_args.push_back(
+              std::make_pair(*(++modeliter), arguments_t(argiter, nextenditer))
+              );
+          argiter = nextenditer;
+        } else if (areValidArguments(types, args)) {
+          // only former matches, don't increment
+          partial_args.push_back(
+              std::make_pair(*modeliter, arguments_t(argiter, enditer))
+              );
+          argiter = enditer;
+        } else
+          ASSERT(0,
+              "CompoundPotential::splitUpArgumentsByModels() - neither type matches its argument bunch.");
+      } else {
+        if (areValidArguments(types, args)) {
+          // only former matches, don't increment
+          partial_args.push_back(
+              std::make_pair(*modeliter, arguments_t(argiter, enditer))
+              );
+          argiter = enditer;
+        } else if (areValidArguments(nexttypes, nextargs)) {
+          // latter matches, increment
+          ++typesiter;
+          partial_args.push_back(
+              std::make_pair(*(++modeliter), arguments_t(argiter, nextenditer))
+              );
+          argiter = nextenditer;
+        }
+      }
+    } else {
+      const SerializablePotential::ParticleTypes_t &types = *typesiter;
+      // we always expect N(N-1)/2 distances for N particle types
+      arguments_t::const_iterator enditer = argiter+(types.size()*(types.size()-1)/2);
+      partial_args.push_back(
+          std::make_pair(*modeliter, arguments_t(argiter, enditer))
+          );
+      argiter = enditer;
+    }
+  }
+
   return partial_args;
 }
@@ -190,18 +301,37 @@
 CompoundPotential::results_t CompoundPotential::operator()(const arguments_t &arguments) const
 {
-  // first, we have to split up the given arguments
-  std::vector<arguments_t> partial_args =
+  /// first, we have to split up the given arguments
+  arguments_by_model_t partial_args =
       splitUpArgumentsByModels(arguments);
-  // then, with each bunch of arguments, we call the specific model
+  // print split up argument list for debugging
+  if (DoLog(4)) {
+    LOG(4, "Arguments by model are: ");
+    for(arguments_by_model_t::const_iterator iter = partial_args.begin();
+        iter != partial_args.end(); ++iter) {
+      LOG(4, "\tModel with " << iter->first->getParameterDimension()
+          << " parameters " << iter->first->getParameters()
+          << " and arguments: " << iter->second);
+    }
+  }
+
+  /// then, with each bunch of arguments, we call the specific model
   results_t results(1,0.);
-  std::vector<results_t> partial_results(models.size());
-  std::transform(
-      models.begin(), models.end(),
-      partial_args.begin(),
-      partial_results.begin(),
-      boost::bind(&FunctionModel::operator(), _1, _2)
-  );
-  std::for_each(partial_results.begin(), partial_results.end(),
-      std::cout << (boost::lambda::_1)[0] << "\t");
+  std::vector<results_t> partial_results;
+  for(arguments_by_model_t::const_iterator iter = partial_args.begin();
+      iter != partial_args.end(); ++iter) {
+    partial_results.push_back(
+        (*iter->first)(iter->second)
+    );
+  }
+  // print partial results for debugging
+  if (DoLog(4)) {
+    std::stringstream output;
+    output << "Partial results are: ";
+    std::for_each(partial_results.begin(), partial_results.end(),
+        output << (boost::lambda::_1)[0] << "\t");
+    LOG(4, output.str());
+  }
+
+  /// Finally, sum up all results and return
   std::for_each(partial_results.begin(), partial_results.end(),
       results[0] += (boost::lambda::_1)[0]);
@@ -212,5 +342,5 @@
 {
   // first, we have to split up the given arguments
-  std::vector<arguments_t> partial_args =
+  arguments_by_model_t partial_args =
       splitUpArgumentsByModels(arguments);
   // then, with each bunch of arguments, we call the specific model
@@ -223,7 +353,7 @@
   std::partial_sum(dimensions.begin(), dimensions.end(), dimensions.begin());
 
-  // look for value not less than index
+  // look for first value greater than index
   std::vector<size_t>::const_iterator iter =
-      std::lower_bound(dimensions.begin(), dimensions.end(), index);
+      std::upper_bound(dimensions.begin(), dimensions.end(), index);
 
   // step forward to same model
@@ -231,13 +361,30 @@
   std::advance(modeliter,
       std::distance(const_cast<const std::vector<size_t> &>(dimensions).begin(), iter) );
-  std::vector<arguments_t>::const_iterator argiter = partial_args.begin();
-  std::advance(argiter,
-      std::distance(const_cast<const std::vector<size_t> &>(dimensions).begin(), iter) );
-
-  // evaluate with correct relative index and return
-  const size_t indexbase = (iter == dimensions.begin()) ? 0 : *(--iter);
-  CompoundPotential::results_t results =
-      (*modeliter)->parameter_derivative(*argiter, index-indexbase);
-  return results;
+
+  CompoundPotential::results_t returnresults;
+  for(arguments_by_model_t::const_iterator argiter = partial_args.begin();
+      argiter != partial_args.end(); ++argiter) {
+    const FunctionModel *model = argiter->first;
+
+    // for every matching model evaluate
+    if (model == *modeliter) {
+      // evaluate with correct relative index and return
+      const size_t indexbase = (iter == dimensions.begin()) ? 0 : *(iter-1);
+      CompoundPotential::results_t results =
+          model->parameter_derivative(argiter->second, index-indexbase);
+
+      // either set results or add
+      if (returnresults.empty())
+        returnresults = results;
+      else
+        std::transform(
+            results.begin(), results.end(),
+            returnresults.begin(),
+            returnresults.begin(),
+            std::plus<FunctionModel::result_t>());
+    }
+  }
+
+  return returnresults;
 }
 
@@ -258,8 +405,10 @@
   BOOST_FOREACH( FunctionModel* model, models) {
     const CompoundPotential::parameters_t params = model->getLowerBoxConstraints();
-    std::copy(params.begin(), params.end(), iter);
     ASSERT( iter != constraints.end(),
-        "CompoundPotential::getParameters() - iter already at end.");
-  }
+        "CompoundPotential::getLowerBoxConstraints() - iter already at end.");
+    iter = std::copy(params.begin(), params.end(), iter);
+  }
+  ASSERT( iter == constraints.end(),
+      "CompoundPotential::getLowerBoxConstraints() - iter not at end.");
   return constraints;
 }
@@ -272,8 +421,10 @@
   BOOST_FOREACH( FunctionModel* model, models) {
     const CompoundPotential::parameters_t params = model->getUpperBoxConstraints();
-    std::copy(params.begin(), params.end(), iter);
     ASSERT( iter != constraints.end(),
-        "CompoundPotential::getParameters() - iter already at end.");
-  }
+        "CompoundPotential::getUpperBoxConstraints() - iter already at end.");
+    iter = std::copy(params.begin(), params.end(), iter);
+  }
+  ASSERT( iter == constraints.end(),
+      "CompoundPotential::getUpperBoxConstraints() - iter not at end.");
   return constraints;
 }
Index: src/Potentials/CompoundPotential.hpp
===================================================================
--- src/Potentials/CompoundPotential.hpp	(revision c5e63a4efcf5739f8379980891648d24490fe8c4)
+++ src/Potentials/CompoundPotential.hpp	(revision 7c10913cd7643f0bb610d97273b66628c8faa10d)
@@ -132,10 +132,22 @@
 
 private:
+  //!> typedef for split up arguments, each associated to a model
+  typedef std::vector< std::pair<FunctionModel *, arguments_t> > arguments_by_model_t;
+
   /** Helper function to split up concatenated arguments for operator() calls.
    *
    * \param arguments arguments to split up
-   * \return vector of partial arguments
+   * \return vector of partial arguments with associated model
    */
-  std::vector<arguments_t> splitUpArgumentsByModels(const arguments_t &arguments) const;
+  arguments_by_model_t splitUpArgumentsByModels(const arguments_t &arguments) const;
+
+  /** Helper function to check whether split up argument bunch matches with types.
+   *
+   * \param types types of potential to check whether args match
+   * \param args vector of argument whose types to check
+   */
+  bool areValidArguments(
+      const SerializablePotential::ParticleTypes_t &types,
+      const arguments_t &args) const;
 
 private:
Index: src/Potentials/unittests/CompoundPotentialUnitTest.cpp
===================================================================
--- src/Potentials/unittests/CompoundPotentialUnitTest.cpp	(revision c5e63a4efcf5739f8379980891648d24490fe8c4)
+++ src/Potentials/unittests/CompoundPotentialUnitTest.cpp	(revision 7c10913cd7643f0bb610d97273b66628c8faa10d)
@@ -92,14 +92,14 @@
   }
 
-  // create graph
+  // create graph (i.e. this simulates a water molecule)
   {
     // add nodes
     nodes +=
-        std::make_pair(FragmentNode(0,1),1),
-        std::make_pair(FragmentNode(1,1),1);
+        std::make_pair(FragmentNode(0,2),1),
+        std::make_pair(FragmentNode(1,1),2);
 
     // add edges
     edges +=
-        std::make_pair(FragmentEdge(0,1),1);
+        std::make_pair(FragmentEdge(0,1),2);
 
     // construct graph
@@ -125,14 +125,14 @@
       4.46595;
   output +=
-      -78.5211,
-      -78.8049,
-      -78.935,
-      -78.9822,
-      -78.9867,
-      -78.971,
-      -78.9475,
-      -78.9225,
-      -78.899,
-      -78.8784;
+      2.*0.467226+d,
+      2.*0.183388+d,
+      2.*0.0532649+d,
+      2.*0.00609808+d,
+      2.*0.00160056+d,
+      2.*0.0172506+d,
+      2.*0.0407952+d,
+      2.*0.0658475+d,
+      2.*0.0893157+d,
+      2.*0.109914+d;
 
   CPPUNIT_ASSERT_EQUAL( input.size(), output.size() );
@@ -157,6 +157,9 @@
   compound.setParameters(params);
   for (size_t index = 0; index < input.size(); ++index) {
-    argument_t arg(argument_t::indices_t(0,1), argument_t::types_t(0,1), input[index]);
-    const double result = compound( FunctionModel::arguments_t(1,arg) )[0];
+    argument_t firstarg(argument_t::indices_t(0,1), argument_t::types_t(0,1), input[index]);
+    argument_t secondarg(argument_t::indices_t(0,2), argument_t::types_t(0,1), input[index]);
+    FunctionModel::arguments_t args;
+    args += firstarg,secondarg;
+    const double result = compound( args )[0];
     CPPUNIT_ASSERT(
         Helpers::isEqual(
@@ -177,6 +180,9 @@
 //  params += d,a,c,D;
 //  compound.setParameters(params);
-//  argument_t arg(argument_t::indices_t(0,1), argument_t::types_t(0,1), c);
-//  const double result = compound.derivative(FunctionModel::arguments_t(1,arg))[0]
+//  argument_t firstarg(argument_t::indices_t(0,1), argument_t::types_t(0,1), input[index]);
+//  argument_t secondarg(argument_t::indices_t(0,2), argument_t::types_t(0,1), input[index]);
+//  FunctionModel::arguments_t args;
+//  args += firstarg,secondarg;
+//  const double result = compound.derivative( args )[0]
 //  CPPUNIT_ASSERT(
 //      Helpers::isEqual(
@@ -197,9 +203,12 @@
   params += d,a,c,D;
   compound.setParameters(params);
-  argument_t arg(argument_t::indices_t(0,1), argument_t::types_t(0,1), c);
-  {
-    const double result =
-        compound.parameter_derivative(
-            FunctionModel::arguments_t(1,arg),
+  argument_t firstarg(argument_t::indices_t(0,1), argument_t::types_t(0,1), c);
+  argument_t secondarg(argument_t::indices_t(0,2), argument_t::types_t(0,1), c);
+  FunctionModel::arguments_t args;
+  args += firstarg,secondarg;
+  {
+    const double result =
+        compound.parameter_derivative(
+            args,
             0)[0];
     CPPUNIT_ASSERT(
@@ -214,5 +223,5 @@
     const double result =
         compound.parameter_derivative(
-            FunctionModel::arguments_t(1,arg),
+            args,
             1)[0];
     CPPUNIT_ASSERT(
@@ -227,5 +236,5 @@
     const double result =
         compound.parameter_derivative(
-            FunctionModel::arguments_t(1,arg),
+            args,
             2)[0];
     CPPUNIT_ASSERT(
@@ -240,5 +249,5 @@
     const double result =
         compound.parameter_derivative(
-            FunctionModel::arguments_t(1,arg),
+            args,
             3)[0];
     CPPUNIT_ASSERT(
