Index: src/Fragmentation/Summation/SetValues/Fragment.cpp
===================================================================
--- src/Fragmentation/Summation/SetValues/Fragment.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/Fragment.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,198 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2012 University of Bonn. All rights reserved.
+ * 
+ *
+ *   This file is part of MoleCuilder.
+ *
+ *    MoleCuilder is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 2 of the License, or
+ *    (at your option) any later version.
+ *
+ *    MoleCuilder is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * Fragment.cpp
+ *
+ *  Created on: Aug 8, 2012
+ *      Author: heber
+ */
+
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+// include headers that implement a archive in simple text format
+// otherwise BOOST_CLASS_EXPORT_IMPLEMENT has no effect
+#include <boost/archive/text_oarchive.hpp>
+#include <boost/archive/text_iarchive.hpp>
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include "Fragment.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+#include <limits>
+
+#include "CodePatterns/Assert.hpp"
+
+Fragment::Fragment()
+{}
+
+Fragment::Fragment(const positions_t &_positions, const charges_t &_charges)
+{
+  ASSERT( _positions.size() == _charges.size(),
+      "Fragment::Fragment() - given charges and positions don't match in size.");
+  positions_t::const_iterator piter = _positions.begin();
+  charges_t::const_iterator citer = _charges.begin();
+  for (;piter != _positions.end(); ++piter, ++citer)
+    nuclei.push_back( make_pair( *piter, *citer) );
+}
+
+
+Fragment& Fragment::operator+=(const Fragment &other)
+{
+  for (nuclei_t::const_iterator iter = other.nuclei.begin();
+      iter != other.nuclei.end(); ++iter) {
+    if (!containsNuclei(*iter)) {
+      nuclei.push_back(*iter);
+    }
+  }
+  return *this;
+}
+
+Fragment& Fragment::operator=(const Fragment &other)
+{
+  if (this != &other) {
+    nuclei.clear();
+    nuclei = other.nuclei;
+  }
+  return *this;
+}
+
+Fragment& Fragment::operator-=(const Fragment &other)
+{
+  for (nuclei_t::const_iterator iter = other.nuclei.begin();
+      iter != other.nuclei.end(); ++iter)
+    removeNuclei(*iter);
+  return *this;
+}
+
+bool Fragment::containsNuclei(const nucleus_t &n) const
+{
+  for (nuclei_t::const_iterator iter = nuclei.begin();
+      iter != nuclei.end(); ++iter)
+    if (Fragment::isPositionEqual(iter->first, n.first))
+      return true;
+  return false;
+}
+
+void Fragment::removeNuclei(const nucleus_t &n)
+{
+  for (nuclei_t::iterator iter = nuclei.begin();
+      iter != nuclei.end(); ++iter)
+    if (Fragment::isPositionEqual(iter->first, n.first)) {
+      nuclei.erase(iter);
+      break;
+    }
+}
+
+Fragment::positions_t Fragment::getPositions() const
+{
+  positions_t positions;
+  for (nuclei_t::const_iterator iter = nuclei.begin();
+      iter != nuclei.end(); ++iter)
+    positions.push_back(iter->first);
+  return positions;
+}
+
+Fragment::charges_t Fragment::getCharges() const
+{
+  charges_t charges;
+  for (nuclei_t::const_iterator iter = nuclei.begin();
+      iter != nuclei.end(); ++iter)
+    charges.push_back(iter->second);
+  return charges;
+}
+
+Fragment::nucleus_t Fragment::createNucleus(const position_t &position, const double charge)
+{
+  return nucleus_t( make_pair( position, charge ) );
+}
+
+bool Fragment::isPositionEqual(const position_t &a, const position_t &b)
+{
+  bool status = true;
+  for (size_t i=0;i<3;++i)
+    status &= fabs(a[i] - b[i]) < std::numeric_limits<double>::epsilon();
+  return status;
+}
+
+bool Fragment::operator==(const Fragment& other) const
+{
+  bool status = true;
+  // compare sizes
+  if (nuclei.size() != other.nuclei.size()) {
+    return false;
+  } else {
+    // if equal, sort, and compare each nuclei
+    Fragment::nuclei_t thisnuclei(nuclei);
+    Fragment::nuclei_t othernuclei(other.nuclei);
+    std::sort(thisnuclei.begin(), thisnuclei.end());
+    std::sort(othernuclei.begin(), othernuclei.end());
+    Fragment::nuclei_t::const_iterator iter = thisnuclei.begin();
+    Fragment::nuclei_t::const_iterator oiter = othernuclei.begin();
+    for (; iter != thisnuclei.end(); ++iter,++oiter)
+      status &= (*iter == *oiter);
+  }
+  return status;
+}
+
+std::ostream & operator<<(std::ostream &ost, const Fragment &f)
+{
+  size_t NucleiCounter = 1;
+  for (Fragment::nuclei_t::const_iterator iter = f.nuclei.begin();
+      iter != f.nuclei.end(); ++iter)
+    ost << "[" << NucleiCounter++ << "]:" << *iter << ", ";
+  return ost;
+}
+
+std::ostream & operator<<(std::ostream &ost, const Fragment::nucleus_t &n)
+{
+  ost << n.first << " with charge " << n.second;
+  return ost;
+}
+
+bool operator==(const Fragment::nucleus_t &a, const Fragment::nucleus_t &b)
+{
+  bool status = true;
+  // check positions
+  status &= Fragment::isPositionEqual(a.first, b.first);
+  status &= fabs(a.second - b.second) < std::numeric_limits<double>::epsilon();
+  return status;
+}
+
+template<> Fragment ZeroInstance<Fragment>()
+{
+  Fragment::positions_t positions;
+  Fragment::charges_t charges;
+  Fragment returnvalue(positions, charges);
+  return returnvalue;
+}
+
+
+// we need to explicitly instantiate the serialization functions
+BOOST_CLASS_EXPORT_IMPLEMENT(Fragment)
Index: src/Fragmentation/Summation/SetValues/Fragment.hpp
===================================================================
--- src/Fragmentation/Summation/SetValues/Fragment.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/Fragment.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,169 @@
+/*
+ * Fragment.hpp
+ *
+ *  Created on: Aug 8, 2012
+ *      Author: heber
+ */
+
+#ifndef FRAGMENT_HPP_
+#define FRAGMENT_HPP_
+
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <boost/serialization/export.hpp>
+#include <boost/serialization/vector.hpp>
+#include <boost/serialization/utility.hpp>
+
+#include <iosfwd>
+#include <vector>
+
+class FragmentTest;
+
+class Fragment {
+  //!> grant ostream operator access
+  friend std::ostream & operator<<(std::ostream &ost, const Fragment &f);
+  //!> grant unit test access
+  friend class FragmentTest;
+public:
+  typedef std::vector<double> position_t;
+  typedef std::vector< position_t > positions_t;
+  typedef double charge_t;
+  typedef std::vector< charge_t > charges_t;
+  typedef std::pair< position_t, charge_t> nucleus_t;
+  typedef std::vector< nucleus_t > nuclei_t;
+
+  /** Default constructor of class Fragment.
+   *
+   */
+  Fragment();
+
+  /** Default constructor of class Fragment.
+   *
+   */
+  Fragment(const nuclei_t &_nuclei) :
+    nuclei(_nuclei)
+  {}
+
+  /** Constructor of class Fragment.
+   *
+   * @param _positions given positions
+   * @param _charges given charges
+   */
+  Fragment(const positions_t &_positions, const charges_t &_charges);
+
+  /** Adding another fragment onto this one.
+   *
+   * \note The operation is area-conserving, i.e. the new area is the sum of
+   * both areas.
+   *
+   * @param other other fragment
+   * @return ref to this instance
+   */
+  Fragment& operator+=(const Fragment &other);
+
+  /** Assignment operator.
+   *
+   * @param other other fragment to make ourselves equal to
+   * @return ref to this instance
+   */
+  Fragment& operator=(const Fragment &other);
+
+  /** Subtracting another fragment from this one.
+   *
+   * @param other other fragment
+   * @return ref to this instance
+   */
+  Fragment& operator-=(const Fragment &other);
+
+  /** Getter for all stored positions.
+   *
+   * @return vector of positions
+   */
+  positions_t getPositions() const;
+
+  /** Getter for all stored charges.
+   *
+   * @return vector of charges
+   */
+  charges_t getCharges() const;
+
+  /** Equality operator.
+   *
+   * @param other other instance to check against
+   * @return true - both are equal, false - some nucleus_t differ
+   */
+  bool operator==(const Fragment& other) const;
+
+  bool operator!=(const Fragment& other) const
+  {
+    return (!(*this == other));
+  }
+
+  /** Creates type nucleus_t from given \a position and \a charge.
+   *
+   * @param position position of nucleus to create
+   * @param charge charge of nucleus to create
+   * @return nucleus with given \a position and \a charge
+   */
+  static nucleus_t createNucleus(const position_t &position, const double charge);
+
+  /** Helper function to check whether two positions are equal.
+   *
+   * @param a first position
+   * @param b second position
+   * @return a equals b within numerical precision
+   */
+  static bool isPositionEqual(const position_t &a, const position_t &b);
+
+private:
+  /** Helper function that checks whether this nuclei \b position is present.
+   *
+   * This operation is \f${\cal O}(n)\f$
+   *
+   * @param n nuclei to check
+   * @return true - is contained, false - is not contained
+   */
+  bool containsNuclei(const nucleus_t &n) const;
+
+  /** Seeks through all nuclei and removes one with matching \b position if found.
+   *
+   * @param n nuclei to remove
+   */
+  void removeNuclei(const nucleus_t &n);
+
+private:
+  nuclei_t nuclei;
+
+private:
+  friend class boost::serialization::access;
+  // serialization
+  template <typename Archive>
+  void serialize(Archive& ar, const unsigned int version)
+  {
+    ar & nuclei;
+  }
+};
+
+// we need to give this class a unique key for serialization
+BOOST_CLASS_EXPORT_KEY(Fragment)
+
+/** Equality operator for two nuclei.
+ *
+ * @param a first nuclei
+ * @param b second nuclei
+ * @return true - both have same position and charge, false - either charge or position is different
+ */
+bool operator==(const Fragment::nucleus_t &a, const Fragment::nucleus_t &b);
+
+std::ostream & operator<<(std::ostream &ost, const Fragment::nucleus_t &n);
+
+std::ostream & operator<<(std::ostream &ost, const Fragment &f);
+
+template<typename T> T ZeroInstance();
+template<> Fragment ZeroInstance<Fragment>();
+
+#endif /* FRAGMENT_HPP_ */
Index: src/Fragmentation/Summation/SetValues/Histogram.cpp
===================================================================
--- src/Fragmentation/Summation/SetValues/Histogram.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/Histogram.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,354 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2012 University of Bonn. All rights reserved.
+ * 
+ *
+ *   This file is part of MoleCuilder.
+ *
+ *    MoleCuilder is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 2 of the License, or
+ *    (at your option) any later version.
+ *
+ *    MoleCuilder is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * Histogram.cpp
+ *
+ *  Created on: Jul 26, 2012
+ *      Author: heber
+ */
+
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include "Histogram.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+#include <sstream>
+#include <string>
+
+#include <boost/bind.hpp>
+#include <boost/foreach.hpp>
+
+#include "CodePatterns/Assert.hpp"
+#include "CodePatterns/Log.hpp"
+
+/** Tiny helper struct to create equally spaced bins with count of zero.
+ *
+ */
+struct BinCreator_t {
+  BinCreator_t( const double lowerend, const double _width ) :
+    currentstart(lowerend),
+    width(_width)
+  {}
+  Histogram::Bin_t operator()() {
+    Histogram::Bin_t bin( make_pair( currentstart, 0.) );
+    currentstart+=width;
+    return bin;
+  }
+private:
+  double currentstart;
+  const double width;
+};
+
+
+// see http://stackoverflow.com/questions/634087/stdcopy-to-stdcout-for-stdpair
+// printing a pair is not easy, especially so with ostream_iterator as it cannot find
+// overload operator<<() as it is not in namespace std.
+template <class T>
+struct PrintPair : public std::unary_function<T, void>
+{
+    std::ostream& os;
+    PrintPair(std::ostream& strm) : os(strm) {}
+
+    void operator()(const T& elem) const
+    {
+      os << "(" << elem.first << "," << elem.second << ") ";
+    }
+};
+
+std::ostream & operator<<(std::ostream &ost, const Histogram::Bin_t &elem)
+{
+  ost << "(" << elem.first << "," << elem.second << ") ";
+  return ost;
+}
+
+std::ostream & operator<<(std::ostream &ost, const Histogram &histogram)
+{
+  for (Histogram::Bins_t::const_iterator iter = histogram.bins.begin();
+      iter != histogram.bins.end(); ++iter)
+    ost << *iter;
+  return ost;
+}
+
+Histogram::Histogram(const samples_t &samples) :
+  binwidth(0.5),
+  offset(0.)
+{
+  if (!samples.empty()) {
+    // set offset to first value
+    const_cast<BinLowerEnd &>(offset) = *samples.begin();
+    // set binwidth to statistical sensible value
+    MinMax_t MinMax = getMinMaxFromSamples(samples);
+    const_cast<double &>(binwidth) = (*(MinMax.second) - *(MinMax.first))/pow(samples.size(), 1./3.);
+    // and add all samples to these histogram bins
+    addSamples(samples);
+  }
+}
+
+Histogram::Histogram(const samples_t &samples, const BinLowerEnd _offset, const double _binwidth) :
+    binwidth(_binwidth),
+    offset(_offset)
+{
+  addSamples(samples);
+}
+
+Histogram::MinMax_t Histogram::getMinMaxFromSamples(const samples_t &samples) const
+{
+  std::pair<samples_t::const_iterator, samples_t::const_iterator> returnpair;
+  returnpair.first = min_element(samples.begin(), samples.end());
+  returnpair.second = max_element(samples.begin(), samples.end());
+  ASSERT((returnpair.second != samples.end()) || (returnpair.first != samples.end()),
+      "Histogram::Histogram() - cannot find min/max despite non-empty range.");
+  return returnpair;
+}
+
+void Histogram::addSamples(const samples_t &samples)
+{
+  // build histogram from samples
+  if (!samples.empty()) {
+    // 1. get min and max and determine width
+    MinMax_t MinMax = getMinMaxFromSamples(samples);
+//    LOG(1, "DEBUG: min is " << *(MinMax.first) << " and max is " << *(MinMax.second) << ".");
+
+    // 2. create each bin
+    {
+      std::vector<Bin_t> vectorbins;
+      const BinLowerEnd HistogramStart = getLowerEnd(*(MinMax.first));
+      BinCreator_t BinCreator( HistogramStart, binwidth );
+      // we need one extra bin for get...Bin()'s find to work properly
+      const int CountBins = ceil((*(MinMax.second) - HistogramStart)/binwidth)+1;
+      vectorbins.resize(CountBins+1, Bin_t( make_pair(0., 0.) ) );
+      std::generate( vectorbins.begin(), vectorbins.end(), BinCreator );
+      bins.insert(vectorbins.begin(), vectorbins.end());
+    }
+
+    // 3. place each sample into bin
+    BOOST_FOREACH( double value, samples) {
+      const Bins_t::iterator biniter = getLowerEndBin(value);
+      ASSERT( biniter != bins.end(),
+          "Histogram::Histogram() - cannot find bin for value from given samples.");
+      // (make bin count normalized, i.e. always equal to surface area of 1)
+      biniter->second += 1./binwidth;
+    }
+    LOG(2, "DEBUG: Bins are " << printBins() << ".");
+  } else {
+//    LOG(1, "INFO: Given vector of samples is empty.");
+  }
+}
+
+std::string Histogram::printBins() const
+{
+  std::stringstream output;
+  std::for_each( bins.begin(), bins.end(), PrintPair<Bin_t>(output));
+  return output.str();
+}
+
+void Histogram::extendMissingBins(const BinLowerEnd LowerEnd, const BinLowerEnd NextLowerEnd)
+{
+  Bins_t::const_iterator loweriter = getLowerEndBin(LowerEnd);
+  Bins_t::const_iterator upperiter = getHigherEndBin(NextLowerEnd);
+  if (loweriter == bins.end()) {
+    // we need bins beneath our first, add them
+    for(BinLowerEnd offset = getLowerEnd(LowerEnd);
+        offset < getLowerEnd(NextLowerEnd);
+        offset += binwidth) {
+      LOG(4, "DEBUG: Extending below at offset " << offset << ".");
+      bins.insert( std::make_pair (offset, 0. ) );
+    }
+    LOG(3, "DEBUG: Bins after adding empties before start are " << printBins() << ".");
+    loweriter = getLowerEndBin(LowerEnd);
+    ASSERT( loweriter != bins.end(),
+        "Histogram::extendMissingBins() - still no lower bound after adding empties.");
+  }
+  if (upperiter == bins.end()) {
+    // we need bins after our last, add them
+    for(BinLowerEnd offset = getLowerEnd(LowerEnd);
+        offset <= getLowerEnd(NextLowerEnd)+(1.5*binwidth); /* for safety we go a little further */
+        offset += binwidth) {
+      LOG(4, "DEBUG: Extending above at offset " << offset << ".");
+      bins.insert( std::make_pair (offset, 0. ) );
+    }
+    LOG(3, "DEBUG: Bins after adding empties after end are " << printBins() << ".");
+    upperiter = getHigherEndBin(NextLowerEnd);
+    ASSERT( upperiter != bins.end(),
+        "Histogram::extendMissingBins() - still no upper bound after adding empties.");
+  }
+}
+
+void Histogram::superposeOtherHistogram(const Histogram &other, const double prefactor)
+{
+  // go through each of the other histogram's bins
+  Bins_t::const_iterator enditer = --other.bins.end(); // (except internal last one)
+  for (Bins_t::const_iterator biniter = other.bins.begin();
+      biniter != enditer; /* we advance ourselves in loop */) {
+    const Bin_t &bin = *biniter;
+    ++biniter;
+    const Bin_t &nextbin = *biniter;
+
+    LOG(4, "DEBUG: Current bin is " << bin << ", next bin is " << nextbin << ".");
+
+    // Check first whether start or end actually fit into our histogram, if not extend.
+    extendMissingBins(bin.first, nextbin.first);
+
+    // The bin will in general not fit into one bin in this histogram, but overlap.
+    // Hence, we determine the contribution of the bin to each bin in this histogram
+    // its overlaps into and add this weight to the bin.
+    Bins_t::const_iterator loweriter = getLowerEndBin(bin.first);
+    Bins_t::const_iterator upperiter = getHigherEndBin(nextbin.first);
+
+    ASSERT( loweriter->first < upperiter->first,
+        "Histogram::superposeOtherHistogram() - the bin range is invalid.");
+    LOG(4, "DEBUG: bin range here is ["
+        << loweriter->first << ","  << upperiter->first << ").");
+
+    // Next, we create a vector of offsets
+    typedef std::vector< BinLowerEnd > offsets_t;
+    offsets_t offsets;
+    {
+      offsets.push_back(bin.first);
+      Bins_t::const_iterator iter = loweriter;
+      for (++iter; iter != upperiter; ++iter)
+        if (offsets.back() != iter->first)
+          offsets.push_back(iter->first);
+      if (offsets.back() != nextbin.first)
+        offsets.push_back(nextbin.first);
+      LOG(3, "DEBUG: Offsets are " << offsets << ".");
+    }
+
+    // then, we go through the offsets but the last one and add the respective area
+    {
+      offsets_t::const_iterator iter = offsets.begin();
+      offsets_t::const_iterator nextiter = ++offsets.begin();
+      for (; iter != --offsets.end(); ++iter, ++nextiter) {
+        const double length = *nextiter - *iter;
+        const double weight = bin.second * (length/binwidth);
+        Bins_t::iterator filliter = getLowerEndBin(*iter);
+        filliter->second += prefactor * weight;
+      }
+    }
+
+    {
+      std::stringstream output;
+      std::for_each( bins.begin(), bins.end(), PrintPair<Bin_t>(output));
+      LOG(2, "DEBUG: Bins are after summation " << output.str() << ".");
+    }
+  }
+}
+
+Histogram& Histogram::operator=(const Histogram &other)
+{
+  // check for self-assigment
+  if (&other != this) {
+    bins.clear();
+    const_cast<BinLowerEnd &>(offset) = other.offset;
+    const_cast<double &>(binwidth) = other.binwidth;
+    bins.insert(other.bins.begin(), other.bins.end());
+  }
+  return *this;
+}
+
+Histogram& Histogram::operator+=(const Histogram &other)
+{
+  superposeOtherHistogram(other, +1.);
+  LOG(2, "DEBUG: Bins after addition are " << printBins() << ".");
+  return *this;
+}
+
+Histogram& Histogram::operator-=(const Histogram &other)
+{
+  superposeOtherHistogram(other, -1.);
+  LOG(2, "DEBUG: Bins after subtraction are " << printBins() << ".");
+  return *this;
+}
+
+bool Histogram::isEmpty() const
+{
+  bool status = true;
+  for (Bins_t::const_iterator iter = bins.begin(); iter != bins.end(); ++iter)
+    status &= iter->second == 0;
+  return status;
+}
+
+Histogram::Bins_t::iterator Histogram::getLowerEndBin(const double _value)
+{
+  // lower bound returns key that is equal or greater
+  Bins_t::iterator iter = bins.lower_bound(_value);
+  if (iter != bins.end()) {
+    // if we are not on the boundary we always have to step back by one
+    if (_value != iter->first) {
+      if (iter != bins.begin()) {
+          --iter;
+      } else {
+        iter = bins.end();
+      }
+    } else if (iter == --bins.end()) {
+      // if we actually are on boundary of "last bin", set to end
+      iter = bins.end();
+    }
+  }
+  return iter;
+}
+
+Histogram::Bins_t::iterator Histogram::getHigherEndBin(const double _value)
+{
+  // upper bound return key that is strictly greater
+  Bins_t::iterator iter = bins.upper_bound(_value);
+  // if we are on the boundary we have to step back by one
+  if (iter != bins.end())
+    if (_value == iter->first)
+      if (iter != bins.begin())
+          --iter;
+
+  return iter;
+}
+
+Histogram::BinLowerEnd Histogram::getLowerEnd(const double _value) const
+{
+  BinLowerEnd start = _value - offset;
+  // then divide by binwidth
+  const int integral = floor(start/binwidth);
+  //const double fraction = fmod(start,binwidth);
+  start = offset + binwidth*integral;
+  return start;
+}
+
+double Histogram::area() const
+{
+  double area = 0.;
+  for(Bins_t::const_iterator iter = bins.begin(); iter != bins.end(); ++iter)
+    area += iter->second*binwidth;
+  return area;
+}
+
+template<> Histogram ZeroInstance<Histogram>()
+{
+  Histogram returnvalue;
+  return returnvalue;
+}
Index: src/Fragmentation/Summation/SetValues/Histogram.hpp
===================================================================
--- src/Fragmentation/Summation/SetValues/Histogram.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/Histogram.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,216 @@
+/*
+ * Histogram.hpp
+ *
+ *  Created on: Jul 26, 2012
+ *      Author: heber
+ */
+
+#ifndef HISTOGRAM_HPP_
+#define HISTOGRAM_HPP_
+
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <iosfwd>
+#include <map>
+#include <vector>
+
+class HistogramTest;
+
+/** This class generates a histogram from a given vector of sampled values.
+ *
+ * Most importantly, it also contains operator+=() and operator-=() to perform
+ * sum on the \b histograms.
+ *
+ * This is to be used with the OrthogonalSummation which requires these specific
+ * operator implementations.
+ *
+ */
+class Histogram
+{
+  //!> grant unit test access to private parts
+  friend class HistogramTest;
+  //!> grant output operator access
+  friend std::ostream & operator<<(std::ostream &ost, const Histogram &histogram);
+public:
+  //!> named type for a vector of sample values
+  typedef std::vector<double> samples_t;
+  //!> named type for the start of a bin
+  typedef double BinLowerEnd;
+  //!> named type for the count of a bin, may be fractional
+  typedef double BinWeight;
+  //!> named type for the pair identifying a bin by its lower end and count
+  typedef std::pair< BinLowerEnd, BinWeight> Bin_t;
+private:
+  //!> named type for a vector of bins
+  typedef std::map< BinLowerEnd, BinWeight > Bins_t;
+public:
+  /** Cpnstructor from a given set of sampled data.
+   *
+   * \param samples sampled data to construct histogram from
+   */
+  Histogram(const samples_t &samples);
+
+  /** Default constructor.
+   *
+   */
+  Histogram() :
+    binwidth(0.5),
+    offset(0.)
+  {}
+
+  /** Constructor for class Histogram.
+   *
+   * @param offset where the bins should start
+   * @param binwidth width of the bins
+   */
+  Histogram(const BinLowerEnd _offset, const double _binwidth) :
+    binwidth(_binwidth),
+    offset(_offset)
+  {}
+
+  /** Constructor for class Histogram.
+   *
+   * @param samples samples to put in the histogram.
+   * @param offset where the bins should start
+   * @param binwidth width of the bins
+   */
+  Histogram(const samples_t &samples, const BinLowerEnd _offset, const double _binwidth);
+
+  /** Adding another histogram onto this one.
+   *
+   * \note The operation is area-conserving, i.e. the new area is the sum of
+   * both areas.
+   *
+   * @param other other histogram
+   * @return ref to this instance
+   */
+  Histogram& operator+=(const Histogram &other);
+
+  /** Assignment operator.
+   *
+   * @param other other histogram to make ourselves equal to
+   * @return ref to this instance
+   */
+  Histogram& operator=(const Histogram &other);
+
+  /** Subtracting another histogram from this one.
+   *
+   * \note The operation is area-conserving, i.e. the new area is the
+   * difference of both areas.
+   *
+   * @param other other histogram
+   * @return ref to this instance
+   */
+  Histogram& operator-=(const Histogram &other);
+
+  /** States whether each bin of this histogram has count of zero.
+   *
+   * @return true - all bins are zero, false - else
+   */
+  bool isEmpty() const;
+
+  /** Sums all found weights times the Histogram::binwidth.
+   *
+   * @return sum over all weights times width.
+   */
+  double area() const;
+
+  /** Adds given vector of samples to the histogram.
+   *
+   * @param samples vector of samples to add
+   */
+  void addSamples(const samples_t &samples);
+
+private:
+  /** Returns the iterator to the bin representing the lower end into which the \a value fits.
+   *
+   * @param _value value to fit.
+   * @return iterator to bin or to bins.end() if not matching
+   */
+  Bins_t::iterator getLowerEndBin(const double _value);
+
+  /** Returns the iterator to the bin representing the upper end into which the \a value fits.
+   *
+   * @param _value value to fit.
+   * @return iterator to bin or to bins.end() if not matching
+   */
+  Bins_t::iterator getHigherEndBin(const double _value);
+
+  /** Returns the lower end regardless of whether such a bin exists.
+   *
+   * @param _value value to place into a bin
+   * @return start of would-be bin to contain this \a _value
+   */
+  BinLowerEnd getLowerEnd(const double _value) const;
+
+  /** Helper function that contains all the logic of how to superpose two
+   * histograms.
+   *
+   * Is called by Histogram::operator+=() and Histogram::operator-=()
+   *
+   * @param other other histogram
+   * @param prefactor +1. is then addition, -1. is subtraction.
+   */
+  void superposeOtherHistogram(const Histogram &other, const double prefactor);
+
+  /** Helper function to add missing bins in superposition operation.
+   *
+   *  We add here enough bins, initialized to weight zero such that the bin
+   *  of [LowerEnd, NextLowerEnd) fully fits. Does nothing if the bins are
+   *  already present.
+   *
+   * @param LowerEnd lowerend of the other bin (to superpose)
+   * @param NextLowerEnd lower end of bin (to superpose) next adjacent to other
+   */
+  void extendMissingBins(const BinLowerEnd LowerEnd, const BinLowerEnd NextLowerEnd);
+
+  /** Helper function to print BinLowerEnd and BinWeight to a string for every bin.
+   *
+   * @return string containing information on each consecutive bin
+   */
+  std::string printBins() const;
+
+  //!> pair for min and max of given samples
+  typedef std::pair<samples_t::const_iterator, samples_t::const_iterator> MinMax_t;
+
+  /** Helper function that returns minimum and maximum in given \a samples
+   *
+   * @param samples sampled values
+   * @return pair of min and max value
+   */
+  MinMax_t getMinMaxFromSamples(const samples_t &samples) const;
+
+private:
+  //!> vector of bins containing the histogram
+  Bins_t bins;
+  //!> width of bin
+  const double binwidth;
+  //!> offset for the start of the bins
+  const BinLowerEnd offset;
+};
+
+/** Function to print an arbitrary pair to ostream.
+ *
+ * @param ost output stream
+ * @param elem element to print
+ * @return ref to given ostream for concatenation
+ */
+std::ostream & operator<<(std::ostream &ost, const Histogram::Bin_t &elem);
+
+/** Function to print histogram to ostream.
+ *
+ * @param ost output stream
+ * @param histogram histogram to print
+ * @return ref to given ostream for concatenation
+ */
+std::ostream & operator<<(std::ostream &ost, const Histogram &histogram);
+
+template<typename T> T ZeroInstance();
+template<> Histogram ZeroInstance<Histogram>();
+
+
+#endif /* HISTOGRAM_HPP_ */
Index: src/Fragmentation/Summation/SetValues/IndexedVectors.cpp
===================================================================
--- src/Fragmentation/Summation/SetValues/IndexedVectors.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/IndexedVectors.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,116 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2012 University of Bonn. All rights reserved.
+ * 
+ *
+ *   This file is part of MoleCuilder.
+ *
+ *    MoleCuilder is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 2 of the License, or
+ *    (at your option) any later version.
+ *
+ *    MoleCuilder is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * IndexedVectors.cpp
+ *
+ *  Created on: 29.07.2012
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include "Fragmentation/Summation/SetValues/IndexedVectors.hpp"
+
+#include <iostream>
+
+#include "CodePatterns/Assert.hpp"
+
+// static member variables
+const size_t IndexedVectors::FixedSize = 3;
+const IndexedVectors::vector_t IndexedVectors::nullvector( FixedSize, 0.);
+
+IndexedVectors::IndexedVectors(const indices_t &_indices, const vectors_t &_vectors)
+{
+  ASSERT(_indices.size() == _vectors.size(),
+      "IndexedVectors::IndexedVectors() - vector of indices and vectors don't match in size.");
+  indices_t::const_iterator indexiter = _indices.begin();
+  vectors_t::const_iterator vectoriter = _vectors.begin();
+  for (; vectoriter != _vectors.end(); ++indexiter, ++vectoriter) {
+    ASSERT( vectoriter->size() == FixedSize,
+        "IndexedVectors::IndexedVectors() - vector to instance "
+        +toString(*indexiter)+" has size "
+        +toString(vectoriter->size())+" different to FixedSize "
+        +toString(FixedSize)+".");
+    if (*indexiter != (size_t)DropIndex) { // skip all force vectors associated to -1
+#ifndef NDEBUG
+      std::pair<indexedvectors_t::iterator, bool> inserter =
+#endif
+      vectors.insert( std::make_pair( *indexiter, *vectoriter) );
+      ASSERT( inserter.second,
+          "IndexedVectors::IndexedVectors() - index "
+          +toString(inserter.first->first)+" already present with vector "
+          +toString(inserter.first->second)+".");
+    }
+  }
+}
+
+IndexedVectors& IndexedVectors::operator=(const IndexedVectors &other)
+{
+  // check for self-assignment
+  if (this != &other) {
+    vectors.clear();
+    vectors = other.vectors;
+  }
+  return *this;
+}
+
+void IndexedVectors::superposeOtherIndexedVectors(const IndexedVectors &other, const double prefactor)
+{
+  for (indexedvectors_t::const_iterator otheriter = other.vectors.begin();
+      otheriter != other.vectors.end(); ++otheriter) {
+    indexedvectors_t::iterator iter = vectors.find(otheriter->first);
+    if (iter == vectors.end()) {
+      // index does not exist
+      std::pair<indexedvectors_t::iterator, bool> inserter =
+          vectors.insert( std::make_pair( otheriter->first, nullvector) );
+      ASSERT( inserter.second,
+          "IndexedVectors::superposeOtherIndexedVectors() - index is present though unfound before?");
+      iter = inserter.first;
+    }
+    // now simply vector to index from other instance add with prefactor
+    vector_t::iterator vectoriter = iter->second.begin();
+    vector_t::const_iterator othervectoriter = otheriter->second.begin();
+    for (;vectoriter != iter->second.end(); ++vectoriter, ++othervectoriter) {
+      *vectoriter += prefactor * (*othervectoriter);
+    }
+  }
+}
+
+std::ostream & operator<<(std::ostream &ost, const IndexedVectors &other)
+{
+  for (IndexedVectors::indexedvectors_t::const_iterator iter = other.vectors.begin();
+      iter != other.vectors.end(); ++iter)
+    ost << "(" << iter->first << "," << iter->second << ") ";
+  return ost;
+}
+
+template<> IndexedVectors ZeroInstance<IndexedVectors>()
+{
+  IndexedVectors returnvalue;
+  return returnvalue;
+}
Index: src/Fragmentation/Summation/SetValues/IndexedVectors.hpp
===================================================================
--- src/Fragmentation/Summation/SetValues/IndexedVectors.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/IndexedVectors.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,134 @@
+/*
+ * IndexedVectors.hpp
+ *
+ *  Created on: 29.07.2012
+ *      Author: heber
+ */
+
+#ifndef INDEXEDVECTORS_HPP_
+#define INDEXEDVECTORS_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <iosfwd>
+#include <map>
+#include <vector>
+
+class IndexedVectorsTest;
+
+/** IndexedVectors represents a class that contains a a set of vectors,
+ * each associated to a specific index. When adding or subtracting only
+ * the ones are combined that have matching indices.
+ *
+ * This is needed for summing up force vectors per nuclei obtained from
+ * fragment calculations.
+ *
+ */
+class IndexedVectors
+{
+  //!> grant unit test access to private parts
+  friend class IndexedVectorsTest;
+public:
+  //!> typedef for a single vector
+  typedef std::vector<double> vector_t;
+  //!> typedef for the index type
+  typedef size_t index_t;
+  //!> typedef for the indices matching the bunch of vectors
+  typedef std::vector<vector_t> vectors_t;
+  //!> typedef for the ordered indices matching the bunch of vectors
+  typedef std::vector<index_t> indices_t;
+  //!> typedef for a bunch of indexed vectors
+  typedef std::map<index_t, vector_t> indexedvectors_t;
+
+  enum SpecificIndices_t {
+    DropIndex = -1
+  };
+
+  /** Default constructor for class IndexedVectors.
+   *
+   */
+  IndexedVectors() {}
+
+  /** Constructor for class IndexedVectors.
+   *
+   * We construct the internal map from \a _indices and \a _vectors. For
+   * every index -1 contained in \a _indices the respective vector in
+   * \a _vectors is \b not added but silently dropped.
+   *
+   * \param _indices index to each vector
+   * \param _vectors vectors
+   */
+  IndexedVectors(const indices_t &_indices, const vectors_t &_vectors);
+
+  /** Assignment operator.
+   *
+   * \note This is required to place IndexedVectors in STL containers.
+   *
+   * \param other other instance to assign this one to
+   * \return ref to this instance
+   */
+  IndexedVectors& operator=(const IndexedVectors &other);
+
+  /** Addition operator with another IndexedVector instance \a other.
+   *
+   * \param other other instance to sum onto this one.
+   * \return ref to this instance
+   */
+  IndexedVectors& operator+=(const IndexedVectors &other)
+  {
+    superposeOtherIndexedVectors(other, +1.);
+    return *this;
+  }
+
+  /** Subtraction operator with another IndexedVector instance \a other.
+   *
+   * \param other other instance to subtract from this one.
+   * \return ref to this instance
+   */
+  IndexedVectors& operator-=(const IndexedVectors &other)
+  {
+    superposeOtherIndexedVectors(other, -1.);
+    return *this;
+  }
+
+private:
+  /** Helper function that contains all the logic of how to superpose two
+   * indexed vectors.
+   *
+   * Is called by IndexedVectors::operator+=() and IndexedVectors::operator-=()
+   *
+   * @param other other histogram
+   * @param prefactor +1. is then addition, -1. is subtraction.
+   */
+  void superposeOtherIndexedVectors(const IndexedVectors &other, const double prefactor);
+
+private:
+  //!> internal map with all indexed vectors
+  indexedvectors_t vectors;
+  //!> fixed size of all vector_t
+  static const size_t FixedSize;
+  //!> static instance representing a null vector
+  static const vector_t nullvector;
+
+  //!> grant access to output operator
+  friend std::ostream & operator<<(std::ostream &ost, const IndexedVectors &other);
+};
+
+/** Output operator for IndexedVector.
+ *
+ *  Prints a space-separated list of all members as "(index, vector)".
+ *
+ * \param ost output stream to print to
+ * \param other instance to print
+ * \return ref to ost for concatenation
+ */
+std::ostream & operator<<(std::ostream &ost, const IndexedVectors &other);
+
+template<typename T> T ZeroInstance();
+template<> IndexedVectors ZeroInstance<IndexedVectors>();
+
+
+#endif /* INDEXEDVECTORS_HPP_ */
Index: src/Fragmentation/Summation/SetValues/Makefile.am
===================================================================
--- src/Fragmentation/Summation/SetValues/Makefile.am	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/Makefile.am	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,62 @@
+# PLEASE adhere to the alphabetical ordering in this Makefile!
+# Also indentation by a single tab
+
+FRAGMENTATIONSETVALUESSOURCE = \
+	Fragmentation/Summation/SetValues/Fragment.cpp \
+	Fragmentation/Summation/SetValues/Histogram.cpp \
+	Fragmentation/Summation/SetValues/IndexedVectors.cpp \
+	Fragmentation/Summation/SetValues/SamplingGrid.cpp \
+	Fragmentation/Summation/SetValues/SamplingGridProperties.cpp
+
+FRAGMENTATIONSETVALUESHEADER = \
+	Fragmentation/Summation/SetValues/Fragment.hpp \
+	Fragmentation/Summation/SetValues/Histogram.hpp \
+	Fragmentation/Summation/SetValues/IndexedVectors.hpp \
+	Fragmentation/Summation/SetValues/SamplingGrid.hpp \
+	Fragmentation/Summation/SetValues/SamplingGrid_inline.hpp \
+	Fragmentation/Summation/SetValues/SamplingGridProperties.hpp
+
+noinst_LTLIBRARIES += libMolecuilderFragmentationSetValues.la
+libMolecuilderFragmentationSetValues_la_includedir = $(includedir)/MoleCuilder/Fragmentation/Summation/SetValues/
+libMolecuilderFragmentationSetValues_la_CPPFLAGS = ${BOOST_CPPFLAGS} ${CodePatterns_CFLAGS}
+libMolecuilderFragmentationSetValues_la_LDFLAGS = -ldl ${BOOST_LDFLAGS}
+libMolecuilderFragmentationSetValues_la_LIBADD = \
+	${CodePatterns_LIBS}
+
+nobase_libMolecuilderFragmentationSetValues_la_include_HEADERS = ${FRAGMENTATIONSETVALUESHEADER}
+
+## Define the source file list for the "libexample-@MOLECUILDER_API_VERSION@.la"
+## target.  Note that @MOLECUILDER_API_VERSION@ is not interpreted by Automake and
+## will therefore be treated as if it were literally part of the target name,
+## and the variable name derived from that.
+## The file extension .cc is recognized by Automake, and makes it produce
+## rules which invoke the C++ compiler to produce a libtool object file (.lo)
+## from each source file.  Note that it is not necessary to list header files
+## which are already listed elsewhere in a _HEADERS variable assignment.
+libMolecuilderFragmentationSetValues_la_SOURCES = ${FRAGMENTATIONSETVALUESSOURCE}
+
+## Instruct libtool to include ABI version information in the generated shared
+## library file (.so).  The library ABI version is defined in configure.ac, so
+## that all version information is kept in one place.
+#libMolecuilderFragmentationSetValues_la_LDFLAGS += -version-info $(MOLECUILDER_SO_VERSION)
+
+## The generated configuration header is installed in its own subdirectory of
+## $(libdir).  The reason for this is that the configuration information put
+## into this header file describes the target platform the installed library
+## has been built for.  Thus the file must not be installed into a location
+## intended for architecture-independent files, as defined by the Filesystem
+## Hierarchy Standard (FHS).
+## The nodist_ prefix instructs Automake to not generate rules for including
+## the listed files in the distribution on 'make dist'.  Files that are listed
+## in _HEADERS variables are normally included in the distribution, but the
+## configuration header file is generated at configure time and should not be
+## shipped with the source tarball.
+#libMolecuilderFragmentationSetValues_libincludedir = $(libdir)/MoleCuilder/include
+#nodist_libMolecuilderFragmentationSetValues_libinclude_HEADERS = $(top_builddir)/libmolecuilder_config.h
+
+## Install the generated pkg-config file (.pc) into the expected location for
+## architecture-dependent package configuration information.  Occasionally,
+## pkg-config files are also used for architecture-independent data packages,
+## in which case the correct install location would be $(datadir)/pkgconfig.
+#pkgconfigdir = $(libdir)/pkgconfig
+#pkgconfig_DATA = $(top_builddir)/MoleCuilder.pc
Index: src/Fragmentation/Summation/SetValues/SamplingGrid.cpp
===================================================================
--- src/Fragmentation/Summation/SetValues/SamplingGrid.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/SamplingGrid.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,599 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2012 University of Bonn. All rights reserved.
+ * 
+ *
+ *   This file is part of MoleCuilder.
+ *
+ *    MoleCuilder is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 2 of the License, or
+ *    (at your option) any later version.
+ *
+ *    MoleCuilder is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * SamplingGrid.cpp
+ *
+ *  Created on: 25.07.2012
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+// include headers that implement a archive in simple text format
+// otherwise BOOST_CLASS_EXPORT_IMPLEMENT has no effect
+#include <boost/archive/text_oarchive.hpp>
+#include <boost/archive/text_iarchive.hpp>
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include "Fragmentation/Summation/SetValues/SamplingGrid.hpp"
+
+#include <boost/bind.hpp>
+#include <limits>
+
+#include "CodePatterns/Assert.hpp"
+#include "CodePatterns/Log.hpp"
+
+// static instances
+const double SamplingGrid::zeroOffset[3] = { 0., 0., 0. };
+
+SamplingGrid::SamplingGrid() :
+    SamplingGridProperties()
+{
+  setWindowSize(zeroOffset, zeroOffset);
+  ASSERT( getWindowGridPoints() == (size_t)0,
+      "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
+}
+
+SamplingGrid::SamplingGrid(const double _begin[3],
+    const double _end[3],
+    const int _level) :
+  SamplingGridProperties(_begin, _end, _level)
+{
+  setWindowSize(zeroOffset, zeroOffset);
+  ASSERT( getWindowGridPoints() == (size_t)0,
+      "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
+}
+
+SamplingGrid::SamplingGrid(const double _begin[3],
+    const double _end[3],
+    const int _level,
+    const sampledvalues_t &_sampled_grid) :
+  SamplingGridProperties(_begin, _end, _level),
+  sampled_grid(_sampled_grid)
+{
+  setWindowSize(_begin, _end);
+  ASSERT( getWindowGridPoints() == (size_t)_sampled_grid.size(),
+      "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
+}
+
+SamplingGrid::SamplingGrid(const SamplingGrid &_grid) :
+  SamplingGridProperties(_grid),
+  sampled_grid(_grid.sampled_grid)
+{
+  setWindowSize(_grid.begin_window, _grid.end_window);
+  ASSERT( getWindowGridPoints() == _grid.getWindowGridPoints(),
+      "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
+}
+
+SamplingGrid::SamplingGrid(const SamplingGridProperties &_props) :
+  SamplingGridProperties(_props)
+{
+  setWindowSize(zeroOffset, zeroOffset);
+  ASSERT( getWindowGridPoints() == (size_t)0,
+      "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
+}
+
+SamplingGrid::SamplingGrid(
+    const SamplingGridProperties &_props,
+    const sampledvalues_t &_sampled_grid) :
+  SamplingGridProperties(_props),
+  sampled_grid(_sampled_grid)
+{
+  setWindowSize(_props.begin, _props.end);
+  ASSERT( getWindowGridPoints() == (size_t)_sampled_grid.size(),
+      "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
+}
+
+SamplingGrid::~SamplingGrid()
+{}
+
+bool SamplingGrid::isCongruent(const SamplingGrid &_props) const
+{
+  bool status = true;
+  status &= (static_cast<const SamplingGridProperties &>(*this) ==
+      static_cast<const SamplingGridProperties &>(_props));
+  for(size_t i = 0; i<3; ++i) {
+    status &= begin_window[i] == _props.begin_window[i];
+    status &= end_window[i] == _props.end_window[i];
+  }
+  return status;
+}
+
+SamplingGrid& SamplingGrid::operator=(const SamplingGrid& other)
+{
+  // check for self-assignment
+  if (this != &other) {
+    static_cast<SamplingGridProperties &>(*this) = other;
+    setWindowSize(other.begin_window, other.end_window);
+    sampled_grid = other.sampled_grid;
+  }
+  return *this;
+}
+
+static void multiplyElements(
+    double &dest,
+    const double &source,
+    const double prefactor)
+{
+  dest *= prefactor*(source);
+}
+
+
+SamplingGrid& SamplingGrid::operator*=(const SamplingGrid& other)
+{
+  // check that grids are compatible
+  if (isCompatible(other)) {
+    /// get minimum of window
+    double min_begin_window[3];
+    double min_end_window[3];
+    bool doShrink = false;
+    for (size_t index=0; index<3;++index) {
+      if (begin_window[index] <= other.begin_window[index]) {
+        min_begin_window[index] = other.begin_window[index];
+        doShrink = true;
+      } else {
+        min_begin_window[index] = begin_window[index];
+      }
+      if (end_window[index] <= other.end_window[index]) {
+        min_end_window[index] = end_window[index];
+      } else {
+        min_end_window[index] = other.end_window[index];
+        doShrink = true;
+      }
+    }
+    LOG(2, "DEBUG: min begin is " << min_begin_window[0] << "," << min_begin_window[1] << "," << min_begin_window[2] << ".");
+    LOG(2, "DEBUG: min end is " << min_end_window[0] << "," << min_end_window[1] << "," << min_end_window[2] << ".");
+    if (doShrink)
+      shrinkWindow(min_begin_window, min_end_window);
+    addWindowOntoWindow(
+        other.begin_window,
+        other.end_window,
+        begin_window,
+        end_window,
+        sampled_grid,
+        other.sampled_grid,
+        boost::bind(multiplyElements, _1, _2, 1.),
+        sourcewindow);
+  } else {
+    ASSERT(0, "SamplingGrid::operator*=() - multiplying incongruent grids is so far not in the cards.");
+  }
+  return *this;
+}
+
+void SamplingGrid::superposeOtherGrids(const SamplingGrid &other, const double prefactor)
+{
+  /// check that grids are compatible
+  if (isCompatible(other)) {
+    /// get maximum of window
+    double max_begin_window[3];
+    double max_end_window[3];
+    bool doExtend = false;
+    for (size_t index=0; index<3;++index) {
+      if (begin_window[index] >= other.begin_window[index]) {
+        max_begin_window[index] = other.begin_window[index];
+        doExtend = true;
+      } else {
+        max_begin_window[index] = begin_window[index];
+      }
+      if (end_window[index] >= other.end_window[index]) {
+        max_end_window[index] = end_window[index];
+      } else {
+        max_end_window[index] = other.end_window[index];
+        doExtend = true;
+      }
+    }
+    LOG(2, "DEBUG: max begin is " << max_begin_window[0] << "," << max_begin_window[1] << "," << max_begin_window[2] << ".");
+    LOG(2, "DEBUG: max end is " << max_end_window[0] << "," << max_end_window[1] << "," << max_end_window[2] << ".");
+    if (doExtend)
+      extendWindow(max_begin_window, max_end_window);
+    /// and copy other into larger window, too
+    addOntoWindow(other.begin_window, other.end_window, other.sampled_grid, prefactor);
+  } else {
+    ASSERT(0, "SamplingGrid::superposeOtherGrids() - superposing incompatible grids is so far not in the cards.");
+  }
+}
+
+const size_t SamplingGrid::getWindowGridPointsPerAxis(const size_t axis) const
+{
+  static const double round_offset(
+      (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
+          0.5 : 0.); // need offset to get to round_toward_nearest behavior
+//  const double total = getTotalLengthPerAxis(axis);
+  const double delta = getDeltaPerAxis(axis);
+  if (delta == 0)
+    return 0;
+  const double length = getWindowLengthPerAxis(axis);
+  if (length == 0)
+    return 0;
+  return (size_t)(length/delta+round_offset);
+}
+
+double SamplingGrid::integral() const
+{
+  const double volume_element = getVolume()/(double)getTotalGridPoints();
+  double int_value = 0.;
+  for (sampledvalues_t::const_iterator iter = sampled_grid.begin();
+      iter != sampled_grid.end();
+      ++iter)
+    int_value += *iter;
+  int_value *= volume_element;
+  LOG(2, "DEBUG: SamplingGrid::integral() is " << scientific << setprecision(13) << int_value << ".");
+  return int_value;
+}
+
+double SamplingGrid::integral(const SamplingGrid &weight) const
+{
+  if (isCompatible(weight)) {
+    const double volume_element = getVolume()/(double)getTotalGridPoints();
+    double int_value = 0.;
+    sampledvalues_t::const_iterator iter = sampled_grid.begin();
+    sampledvalues_t::const_iterator weightiter = weight.sampled_grid.begin();
+    for (;iter != sampled_grid.end();++iter,++weightiter)
+      int_value += (*weightiter) * (*iter);
+    int_value *= volume_element;
+    //LOG(2, "DEBUG: SamplingGrid::integral() is " << scientific << setprecision(13) << int_value << ".");
+    return int_value;
+  } else
+    return 0.;
+}
+
+double SamplingGrid::getNearestLowerGridPoint(
+    const double value, const size_t axis) const
+{
+  const double length = end[axis] - begin[axis];
+  if ((fabs(length) < std::numeric_limits<double>::epsilon()) || (getGridPointsPerAxis() == 0))
+    return begin[axis];
+  if (value > end[axis])
+    return end[axis];
+  const double offset = value - begin[axis];
+  if (offset < 0)
+    return begin[axis];
+  // modify value a little if value actually has been on a grid point but
+  // perturbed by numerical rounding: here up, as we always go lower
+  const double factor =
+      floor((double)getGridPointsPerAxis()*(offset/length)
+          +std::numeric_limits<double>::epsilon()*1.e4);
+  return factor*(length/(double)getGridPointsPerAxis());
+}
+
+double SamplingGrid::getNearestHigherGridPoint(
+    const double value, const size_t axis) const
+{
+  const double length = end[axis] - begin[axis];
+  if ((fabs(length) < std::numeric_limits<double>::epsilon()) || (getGridPointsPerAxis() == 0))
+    return end[axis];
+  if (value > end[axis])
+    return end[axis];
+  const double offset = value - begin[axis];
+  if (offset < 0)
+    return begin[axis];
+  // modify value a little if value actually has been on a grid point but
+  // perturbed by numerical rounding: here down, as we always go higher
+  const double factor =
+      ceil((double)getGridPointsPerAxis()*(offset/length)
+          -std::numeric_limits<double>::epsilon()*1.e4);
+  return factor*(length/(double)getGridPointsPerAxis());
+}
+
+void SamplingGrid::setWindowSize(
+    const double _begin_window[3],
+    const double _end_window[3])
+{
+  for (size_t index=0;index<3;++index) {
+    begin_window[index] = getNearestLowerGridPoint(_begin_window[index], index);
+    ASSERT( begin_window[index] >= begin[index],
+        "SamplingGrid::setWindowSize() - window starts earlier than domain on "
+        +toString(index)+"th component.");
+    end_window[index] = getNearestHigherGridPoint(_end_window[index], index);
+    ASSERT( end_window[index] <= end[index],
+        "SamplingGrid::setWindowSize() - window ends later than domain on "
+        +toString(index)+"th component.");
+  }
+}
+
+void SamplingGrid::setWindow(
+    const double _begin_window[3],
+    const double _end_window[3])
+{
+  setWindowSize(_begin_window, _end_window);
+  const size_t gridpoints_window = getWindowGridPoints();
+  sampled_grid.clear();
+  sampled_grid.resize(gridpoints_window, 0.);
+}
+
+void SamplingGrid::setDomain(
+    const double _begin[3],
+    const double _end[3])
+{
+  setDomainSize(_begin, _end);
+  setWindowSize(_begin, _end);
+  const size_t gridpoints = getTotalGridPoints();
+  sampled_grid.resize(gridpoints, 0.);
+}
+
+void SamplingGrid::extendWindow(
+    const double _begin_window[3],
+    const double _end_window[3])
+{
+#ifndef NDEBUG
+  for(size_t index=0;index < 3; ++index) {
+    // check that we truly have to extend the window
+    ASSERT ( begin_window[index] >= _begin_window[index],
+        "SamplingGrid::extendWindow() - component "+toString(index)+
+        " of window start is greater than old value.");
+    ASSERT ( end_window[index] <= _end_window[index],
+        "SamplingGrid::extendWindow() - component "+toString(index)+
+        " of window end is less than old value.");
+
+    // check that we are still less than domain
+    ASSERT ( _begin_window[index] >= begin[index],
+        "SamplingGrid::extendWindow() - component "+toString(index)+
+        " of window start is less than domain start.");
+    ASSERT ( _end_window[index] <= end[index],
+        "SamplingGrid::extendWindow() - component "+toString(index)+
+        " of window end is greater than domain end.");
+  }
+#endif
+  // copy old window size and values
+  double old_begin_window[3];
+  double old_end_window[3];
+  for(size_t index=0;index<3;++index) {
+    old_begin_window[index] = begin_window[index];
+    old_end_window[index] = end_window[index];
+  }
+  sampledvalues_t old_values(sampled_grid);
+  // set new window
+  setWindow(_begin_window,_end_window);
+  // now extend it ...
+  addOntoWindow(old_begin_window, old_end_window, old_values, +1.);
+  LOG(2, "DEBUG: Grid after extension is " << sampled_grid << ".");
+}
+
+void SamplingGrid::shrinkWindow(
+    const double _begin_window[3],
+    const double _end_window[3])
+{
+#ifndef NDEBUG
+  for(size_t index=0;index < 3; ++index) {
+    // check that we truly have to shrink the window
+    ASSERT ( begin_window[index] <= _begin_window[index],
+        "SamplingGrid::shrinkWindow() - component "+toString(index)+
+        " of window start is less than old value.");
+    ASSERT ( end_window[index] >= _end_window[index],
+        "SamplingGrid::shrinkWindow() - component "+toString(index)+
+        " of window end is greater than old value.");
+
+    // check that we are still less than domain
+    ASSERT ( _begin_window[index] >= begin[index],
+        "SamplingGrid::shrinkWindow() - component "+toString(index)+
+        " of window start is less than domain start.");
+    ASSERT ( _end_window[index] <= end[index],
+        "SamplingGrid::shrinkWindow() - component "+toString(index)+
+        " of window end is greater than domain end.");
+  }
+#endif
+  // copy old window size and values
+  double old_begin_window[3];
+  double old_end_window[3];
+  for(size_t index=0;index<3;++index) {
+    old_begin_window[index] = begin_window[index];
+    old_end_window[index] = end_window[index];
+  }
+  sampledvalues_t old_values(sampled_grid);
+  // set new window
+  setWindow(_begin_window,_end_window);
+  // now extend it ...
+  addIntoWindow(old_begin_window, old_end_window, old_values, +1.);
+  LOG(2, "DEBUG: Grid after extension is " << sampled_grid << ".");
+}
+
+static void addElements(
+    double &dest,
+    const double &source,
+    const double prefactor)
+{
+  dest += prefactor*(source);
+}
+
+void SamplingGrid::addOntoWindow(
+    const double _begin_window[3],
+    const double _end_window[3],
+    const sampledvalues_t &_sampled_grid,
+    const double prefactor)
+{
+  addWindowOntoWindow(
+      begin_window,
+      end_window,
+      _begin_window,
+      _end_window,
+      sampled_grid,
+      _sampled_grid,
+      boost::bind(addElements, _1, _2, boost::cref(prefactor)),
+      destwindow);
+}
+
+void SamplingGrid::addIntoWindow(
+    const double _begin_window[3],
+    const double _end_window[3],
+    const sampledvalues_t &_sampled_grid,
+    const double prefactor)
+{
+  addWindowOntoWindow(
+      _begin_window,
+      _end_window,
+      begin_window,
+      end_window,
+      sampled_grid,
+      _sampled_grid,
+      boost::bind(addElements, _1, _2, boost::cref(prefactor)),
+      sourcewindow);
+}
+
+void SamplingGrid::addWindowOntoWindow(
+    const double larger_wbegin[3],
+    const double larger_wend[3],
+    const double smaller_wbegin[3],
+    const double smaller_wend[3],
+    sampledvalues_t &dest_sampled_grid,
+    const sampledvalues_t &source_sampled_grid,
+    boost::function<void (double &, const double &)> op,
+    enum eLargerWindow larger_window)
+{
+#ifndef NDEBUG
+  for(size_t index=0;index<3;++index) {
+    ASSERT( smaller_wbegin[index] >= larger_wbegin[index],
+        "SamplingGrid::addWindowOntoWindow() - given smaller window starts earlier than larger window in component "
+        +toString(index)+".");
+    ASSERT( smaller_wend[index] <= larger_wend[index],
+        "SamplingGrid::addWindowOntoWindow() - given smaller window ends later than larger window in component "
+        +toString(index)+".");
+  }
+#endif
+  // the only issue are indices
+  const size_t gridpoints_axis = getGridPointsPerAxis();
+  size_t pre_offset[3];
+  size_t post_offset[3];
+  size_t length[3];
+  size_t total[3];
+  const double round_offset =
+      (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
+          0.5 : 0.; // need offset to get to round_toward_nearest behavior
+  for(size_t index=0;index<3;++index) {
+    // we refrain from using floor/ceil as the window's starts and ends,
+    // the grids have to be compatible (equal level), should always be on
+    // discrete grid point locations.
+    const double delta = (double)gridpoints_axis/(end[index] - begin[index]);
+    // delta is conversion factor from box length to discrete length, i.e. number of points
+    pre_offset[index] = delta*(smaller_wbegin[index] - larger_wbegin[index])+round_offset;
+    length[index] = delta*(smaller_wend[index] - smaller_wbegin[index])+round_offset;
+    post_offset[index] = delta*(larger_wend[index] - smaller_wend[index])+round_offset;
+    total[index] = delta*(larger_wend[index] - larger_wbegin[index])+round_offset;
+    // total is used as safe-guard against loss due to discrete conversion
+    ASSERT( pre_offset[index]+post_offset[index]+length[index] == total[index],
+        "SamplingGrid::addWindowOntoWindow() - pre, post, and length don't sum up to total for "
+        +toString(index)+"th component.");
+  }
+  // assert that calculated lengths match with given vector sizes
+#ifndef NDEBUG
+  const size_t calculated_size = length[0]*length[1]*length[2];
+  if (larger_window == destwindow) {
+    ASSERT( calculated_size == source_sampled_grid.size(),
+        "SamplingGrid::addWindowOntoWindow() - not enough source sampled values given: "
+        +toString(calculated_size)+" != "+toString(source_sampled_grid.size())+".");
+    ASSERT( calculated_size <= dest_sampled_grid.size(),
+        "SamplingGrid::addWindowOntoWindow() - not enough sampled values available: "
+        +toString(calculated_size)+" <= "+toString(dest_sampled_grid.size())+".");
+  } else {
+    ASSERT( calculated_size == dest_sampled_grid.size(),
+        "SamplingGrid::addWindowOntoWindow() - not enough dest sampled values given: "
+        +toString(calculated_size)+" != "+toString(dest_sampled_grid.size())+".");
+    ASSERT( calculated_size <= source_sampled_grid.size(),
+        "SamplingGrid::addWindowOntoWindow() - not enough source sampled values available: "
+        +toString(calculated_size)+" <= "+toString(source_sampled_grid.size())+".");
+  }
+  const size_t total_size = total[0]*total[1]*total[2];
+  if (larger_window == destwindow) {
+    ASSERT( total_size == dest_sampled_grid.size(),
+        "SamplingGrid::addWindowOntoWindow() - total size is not equal to number of present dest points: "
+        +toString(total_size)+" != "+toString(dest_sampled_grid.size())+".");
+  } else {
+    ASSERT( total_size == source_sampled_grid.size(),
+        "SamplingGrid::addWindowOntoWindow() - total size is not equal to number of present source points: "
+        +toString(total_size)+" != "+toString(source_sampled_grid.size())+".");
+  }
+#endif
+  size_t N[3];
+//  size_t counter = 0;
+  sampledvalues_t::iterator destiter = dest_sampled_grid.begin();
+  sampledvalues_t::const_iterator sourceiter = source_sampled_grid.begin();
+  if (larger_window == destwindow)
+    std::advance(destiter, pre_offset[0]*total[1]*total[2]);
+  else
+    std::advance(sourceiter, pre_offset[0]*total[1]*total[2]);
+  for(N[0]=0; N[0] < length[0]; ++N[0]) {
+    if (larger_window == destwindow)
+      std::advance(destiter, pre_offset[1]*total[2]);
+    else
+      std::advance(sourceiter, pre_offset[1]*total[2]);
+    for(N[1]=0; N[1] < length[1]; ++N[1]) {
+      if (larger_window == destwindow)
+        std::advance(destiter, pre_offset[2]);
+      else
+        std::advance(sourceiter, pre_offset[2]);
+      for(N[2]=0; N[2] < length[2]; ++N[2]) {
+        ASSERT( destiter != dest_sampled_grid.end(),
+            "SamplingGrid::addWindowOntoWindow() - destiter is already at end of window.");
+        ASSERT( sourceiter != source_sampled_grid.end(),
+            "SamplingGrid::addWindowOntoWindow() - destiter is already at end of window.");
+        op(*destiter, *sourceiter);
+        ++destiter;
+        ++sourceiter;
+      }
+      if (larger_window == destwindow)
+        std::advance(destiter, post_offset[2]);
+      else
+        std::advance(sourceiter, post_offset[2]);
+    }
+    if (larger_window == destwindow)
+      std::advance(destiter, post_offset[1]*total[2]);
+    else
+      std::advance(sourceiter, post_offset[1]*total[2]);
+  }
+#ifndef NDEBUG
+  if (larger_window == destwindow)
+    std::advance(destiter, post_offset[0]*total[1]*total[2]);
+  else
+    std::advance(sourceiter, post_offset[0]*total[1]*total[2]);
+  ASSERT( destiter == dest_sampled_grid.end(),
+      "SamplingGrid::addWindowOntoWindow() - destiter is not at end of window.");
+  ASSERT( sourceiter == source_sampled_grid.end(),
+      "SamplingGrid::addWindowOntoWindow() - sourceiter is not at end of window.");
+#endif
+  LOG(2, "DEBUG: Grid after adding other is " << dest_sampled_grid << ".");
+}
+
+std::ostream & operator<<(std::ostream &ost, const SamplingGrid& other)
+{
+  ost << "SamplingGrid";
+  ost << " starting at " << other.begin[0] << "," << other.begin[1] << "," << other.begin[2];
+  ost << " ending at " << other.end[0] << "," << other.end[1] << "," << other.end[2];
+  ost << ", window starting at " << other.begin_window[0] << "," << other.begin_window[1] << "," << other.begin_window[2];
+  ost << ", window ending at " << other.end_window[0] << "," << other.end_window[1] << "," << other.end_window[2];
+  ost << ", level of " << other.level;
+  ost << " and integrated value of " << other.integral();
+  return ost;
+}
+
+template<> SamplingGrid ZeroInstance<SamplingGrid>()
+{
+  SamplingGrid returnvalue;
+  return returnvalue;
+}
+
+// we need to explicitly instantiate the serialization functions as
+// its is only serialized through its base class FragmentJob
+BOOST_CLASS_EXPORT_IMPLEMENT(SamplingGrid)
Index: src/Fragmentation/Summation/SetValues/SamplingGrid.hpp
===================================================================
--- src/Fragmentation/Summation/SetValues/SamplingGrid.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/SamplingGrid.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,416 @@
+/*
+ * SamplingGrid.hpp
+ *
+ *  Created on: 25.07.2012
+ *      Author: heber
+ */
+
+#ifndef SAMPLINGGRID_HPP_
+#define SAMPLINGGRID_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <boost/function.hpp>
+#include <iosfwd>
+#include <vector>
+
+#include "boost/serialization/export.hpp"
+#include "boost/serialization/vector.hpp"
+
+#include "Fragmentation/Summation/SetValues/SamplingGridProperties.hpp"
+
+class MPQCData;
+class SamplingGridTest;
+
+/** This class stores a sample function on a three-dimensional grid.
+ *
+ *  \note We do not use boost::multi_array because it is not trivial to serialize.
+ *
+ */
+class SamplingGrid : public SamplingGridProperties
+{
+  //!> grant unit test access to private parts
+  friend class SamplingGridTest;
+  //!> grant output operator access
+  friend std::ostream & operator<<(std::ostream &ost, const SamplingGrid& other);
+public:
+  //!> typedef for sampled values
+  typedef std::vector< double > sampledvalues_t;
+
+  /** Constructor for class SamplingGrid for full window.
+   *
+   * Here, the window of sampled values spans the given domain.
+   *
+   * \param _begin offset for grid per axis
+   * \param _end edge length of grid per axis
+   * \param _level number of grid points in \f$2^{\mathrm{level}}\f$
+   * \param _sampled_grid sample points
+   */
+  SamplingGrid(const double _begin[3],
+      const double _end[3],
+      const int _level,
+      const sampledvalues_t &_sampled_grid);
+
+  /** Constructor for class SamplingGrid for empty window.
+   *
+   * Here, the window is initially of size zero.
+   *
+   * \param _begin offset for grid per axis
+   * \param _end edge length of grid per axis
+   * \param _level number of grid points in \f$2^{\mathrm{level}}\f$
+   */
+  SamplingGrid(const double _begin[3],
+      const double _end[3],
+      const int _level);
+
+  /** Copy constructor for class SamplingGrid with full window from SamplingGridProperties.
+   *
+   * Here, the window is initially empty.
+   *
+   * \param _props properties to copy
+   */
+  SamplingGrid(const SamplingGridProperties &_props);
+
+  /** Copy constructor for class SamplingGrid with empty window from SamplingGridProperties.
+   *
+   * Here, the window must span the whole domain
+   *
+   * \param _props properties to copy
+   * \param _sampled_grid sample points
+   */
+  SamplingGrid(
+      const SamplingGridProperties &_props,
+      const sampledvalues_t &_sampled_grid);
+
+  /** Copy constructor for class SamplingGrid.
+   *
+   * The window of sampled values corresponds to the one on \a _grid.
+   *
+   * \param _grid grid to copy
+   */
+  SamplingGrid(const SamplingGrid &_grid);
+
+  /** default cstor.
+   */
+  SamplingGrid();
+
+  virtual ~SamplingGrid();
+
+  /** Checks whether another instance is consistent with this one.
+   *
+   *  \note Conistency is stronger as grids must have the same window.
+   *
+   * \param _props other properties to check against
+   * \return true - are consistent, false - else
+   */
+  bool isCongruent(const SamplingGrid &_props) const;
+
+  /** Assignment operator.
+   *
+   * \param other other instance to assign ourselves to
+   */
+  SamplingGrid& operator=(const SamplingGrid& other);
+
+  /** Addition operator with another SamplingGrid instance \a other.
+   *
+   * \param other other instance to sum onto this one.
+   * \return ref to this instance
+   */
+  SamplingGrid& operator+=(const SamplingGrid& other)
+  {
+    superposeOtherGrids(other, +1.);
+    return *this;
+  }
+
+  /** Element-wise multiplication operator with another SamplingGrid instance \a other.
+   *
+   * With non-zero windows we have to pay some more attention here.
+   * Now, the windows may not be congruent but we have to find the intersection
+   * of the two windows and then construct the new window only of this size and
+   * multiply. The trick then is not to copy&change the other grid but to
+   * access it properly.
+   *
+   * \param other other instance to sum onto this one.
+   * \return ref to this instance
+   */
+  SamplingGrid& operator*=(const SamplingGrid& other);
+
+  /** Subtraction operator with another SamplingGrid instance \a other.
+   *
+   * \param other other instance to subtract from this one.
+   * \return ref to this instance
+   */
+  SamplingGrid& operator-=(const SamplingGrid& other)
+  {
+    superposeOtherGrids(other, -1.);
+    return *this;
+  }
+
+  /** Returns the numeric integral over the grid.
+   *
+   * @return sum of grid values times volume element
+   */
+  double integral() const;
+
+  /** Returns the numeric integral over the grid where the grid is element-wise multiplied with \a weight.
+   *
+   * @param weight grid of weights
+   * @return sum of grid values weighted by respective element in weight times volume element
+   */
+  double integral(const SamplingGrid &weight) const;
+
+  /** Returns the volume of the domain for this sampled function.
+   *
+   * @return volume
+   */
+  const double getVolume() const;
+
+  /** Returns the total number of gridpoints of the discrete mesh covering the (window) volume.
+   *
+   * @return number of gridpoints sampled_values should have
+   */
+  const size_t getWindowGridPoints() const;
+
+  /** Returns the number of gridpoints of the discrete mesh for the current
+   * window size for given axis \axis.
+   *
+   * \param axis axis to calculate number of gridpoints for
+   * \return number of gridpoints along this axis
+   */
+  const size_t getWindowGridPointsPerAxis(const size_t axis) const;
+
+  /** Returns the total number of gridpoints of the discrete mesh covering the volume.
+   *
+   * @return number of gridpoints sampled_values should have
+   */
+  const size_t getTotalGridPoints() const;
+
+  /** Returns the number of grid points per axis.
+   *
+   * @return number of grid points per unit length
+   */
+  const size_t getGridPointsPerAxis() const;
+
+  /** Returns the length of the domain for the given \a axis.
+   *
+   * \param axis axis for which to get step length
+   * \return domain length for the given axis, i.e. end - begin
+   */
+  const double getTotalLengthPerAxis(const size_t axis) const;
+
+  /** Returns the length of the window for the given \a axis.
+   *
+   * \param axis axis for which to get step length
+   * \return window length for the given axis, i.e. end - begin
+   */
+  const double getWindowLengthPerAxis(const size_t axis) const;
+
+  /** Returns the real step length from one discrete grid point to the next.
+   *
+   * \param axis axis for which to get step length
+   * \return step length for the given axis, as domain length over getGridPointsPerAxis()
+   */
+  const double getDeltaPerAxis(const size_t axis) const;
+
+  /** Returns the volume of the domain covered by the current window.
+   *
+   * @return volume of window
+   */
+  const double getWindowVolume() const;
+
+  /** Sets the size of the window.
+   *
+   * \note also resets the sampled points so far.
+   *
+   * \param _begin_window start of new window
+   * \param _end_window end of window
+   */
+  void setWindow(const double _begin_window[3], const double _end_window[3]);
+
+private:
+  /** Sets the size of the domain.
+   *
+   * \note also resets the sampled points so far and the window.
+   *
+   * \param _begin start of new window
+   * \param _end end of window
+   */
+  void setDomain(const double _begin[3], const double _end[3]);
+
+  /** Sets the size of the domain.
+   *
+   * \note this is just internally used for easing the array setting.
+   *
+   * \param _begin start of domain
+   * \param _end end of domain
+   */
+  void setDomainSize(const double _begin[3], const double _end[3]);
+
+  /** Extends the window while keeping the values.
+   *
+   * \param _begin_window new start of window
+   * \param _end_window new end of window
+   */
+  void extendWindow(const double _begin_window[3], const double _end_window[3]);
+
+  /** Shrinks the window while keeping the values.
+   *
+   * \param _begin_window new start of window
+   * \param _end_window new end of window
+   */
+  void shrinkWindow(const double _begin_window[3], const double _end_window[3]);
+
+  /** Adds another (smaller) window onto the one in this instance.
+   *
+   * \note We assume here that the given window fits on the this one.
+   *
+   * \param _begin_window start of other window
+   * \param _end_window end of other window
+   * \param _sampled_grid other set of sampled values
+   * @param prefactor +1. is then addition, -1. is subtraction.
+   */
+  void addOntoWindow(
+      const double _begin_window[3],
+      const double _end_window[3],
+      const sampledvalues_t &_sampled_grid,
+      const double prefactor);
+
+  /** Adds another (larger) window into the one in this instance.
+   *
+   * \note We assume here that the given window is larger than this one.
+   *
+   * \param _begin_window start of other window
+   * \param _end_window end of other window
+   * \param _sampled_grid other set of sampled values
+   * @param prefactor +1. is then addition, -1. is subtraction.
+   */
+  void addIntoWindow(
+      const double _begin_window[3],
+      const double _end_window[3],
+      const sampledvalues_t &_sampled_grid,
+      const double prefactor);
+
+  /** Enum to help in addWindowOntoWindow() decide which iterator needs to be
+   * advanced.
+   */
+  enum eLargerWindow {
+    destwindow,
+    sourcewindow
+  };
+
+  /** Helper function to copy one (larger) window into a (smaller) window.
+   *
+   * \note Why do we need the extra \a choice? We need to know which window
+   * tuples is associated with which sampled values that are constrained by
+   * one of them being constant, hence the source values
+   *
+   * \param larger_wbegin start of larger window
+   * \param larger_wend end of larger window
+   * \param smaller_wbegin start of smaller window
+   * \param smaller_wend end of smaller window
+   * \param dest_sampled_grid larger set of sampled values
+   * \param source_sampled_grid smaller set of sampled values
+   * \param op operation to perform with the two elements
+   * \param larger_window indicates which is the larger window
+   */
+  void addWindowOntoWindow(
+      const double larger_wbegin[3],
+      const double larger_wend[3],
+      const double smaller_wbegin[3],
+      const double smaller_wend[3],
+      sampledvalues_t &dest_sampled_grid,
+      const sampledvalues_t &source_sampled_grid,
+      boost::function<void (double &, const double &)> op,
+      enum eLargerWindow larger_window);
+
+  /** Helper function that contains all the logic of how to superpose two
+   * grids.
+   *
+   * Is called by SamplingGrid::operator+=() and SamplingGrid::operator-=()
+   *
+   * @param other other histogram
+   * @param prefactor +1. is then addition, -1. is subtraction.
+   */
+  void superposeOtherGrids(const SamplingGrid &other, const double prefactor);
+
+  /** Sets the size of the window.
+   *
+   * \note also resets the sampled points so far.
+   *
+   * \param _begin_window start of new window
+   * \param _end_window end of window
+   */
+  void setWindowSize(const double _begin_window[3], const double _end_window[3]);
+
+  /** Helper function to get point at grid point for given \a axis and less than value.
+   *
+   * @param value value to find nearest grid point to
+   * @param axis axis of the value
+   * @return nearest lower grid point
+   */
+  double getNearestLowerGridPoint(
+      const double value, const size_t axis) const;
+
+  /** Helper function to get point at grid point for given \a axis and greater than value.
+   *
+   * @param value value to find nearest grid point to
+   * @param axis axis of the value
+   * @return nearest lower grid point
+   */
+  double getNearestHigherGridPoint(
+      const double value, const size_t axis) const;
+
+public:
+  /// We do not store the whole grid if many entries are actually zero
+  /// but only a window wherein the sampled function is non-zero.
+
+  //!> sample points of the window
+  sampledvalues_t sampled_grid;
+
+  //!> start of the window relative to SamplingGridProperties::begin and SamplingGridProperties::size
+  double begin_window[3];
+  //!> end of the window relative to SamplingGridProperties::begin and SamplingGridProperties::size
+  double end_window[3];
+
+private:
+  friend class MPQCData;
+
+  friend class boost::serialization::access;
+  // serialization
+  template <typename Archive>
+  void serialize(Archive& ar, const unsigned int version)
+  {
+    ar & boost::serialization::base_object<SamplingGridProperties>(*this);
+    ar & const_cast< sampledvalues_t &>(sampled_grid);
+    for(size_t i=0;i<3;++i) {
+      ar & begin_window[i];
+      ar & end_window[i];
+    }
+  }
+
+  //!> static typedef to use in cstor when no initial values are given
+  static const double zeroOffset[3];
+};
+
+/** Output operator for class SamplingGrid.
+ *
+ * \param ost output stream to print to
+ * \param other instance to print
+ * \return ref to stream for concatenation
+ */
+std::ostream & operator<<(std::ostream &ost, const SamplingGrid& other);
+
+template<typename T> T ZeroInstance();
+template<> SamplingGrid ZeroInstance<SamplingGrid>();
+
+// we need to give this class a unique key for serialization
+// its is only serialized through its base class FragmentJob
+BOOST_CLASS_EXPORT_KEY(SamplingGrid)
+
+// define inline functions
+#include "SamplingGrid_inline.hpp"
+
+#endif /* SAMPLINGGRID_HPP_ */
Index: src/Fragmentation/Summation/SetValues/SamplingGridProperties.cpp
===================================================================
--- src/Fragmentation/Summation/SetValues/SamplingGridProperties.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/SamplingGridProperties.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,103 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2012 University of Bonn. All rights reserved.
+ * 
+ *
+ *   This file is part of MoleCuilder.
+ *
+ *    MoleCuilder is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 2 of the License, or
+ *    (at your option) any later version.
+ *
+ *    MoleCuilder is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * SamplingGrid.cpp
+ *
+ *  Created on: 25.07.2012
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+// include headers that implement a archive in simple text format
+// otherwise BOOST_CLASS_EXPORT_IMPLEMENT has no effect
+#include <boost/archive/text_oarchive.hpp>
+#include <boost/archive/text_iarchive.hpp>
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include "Fragmentation/Summation/SetValues/SamplingGridProperties.hpp"
+
+SamplingGridProperties::SamplingGridProperties(
+    const double _begin[3],
+    const double _end[3],
+    const int _level) :
+  level(_level)
+{
+  for(size_t i=0; i<3; ++i) {
+    begin[i] = _begin[i];
+    end[i] = _end[i];
+  }
+}
+
+SamplingGridProperties::SamplingGridProperties(const SamplingGridProperties &_props) :
+  level(_props.level)
+{
+  for(size_t i=0; i<3; ++i) {
+    begin[i] = _props.begin[i];
+    end[i] = _props.end[i];
+  }
+}
+
+SamplingGridProperties::SamplingGridProperties() :
+  level(0)
+{
+  for(size_t i=0; i<3; ++i) {
+    begin[i] = 0.;
+    end[i] = 0.;
+  }
+}
+
+SamplingGridProperties::~SamplingGridProperties()
+{}
+
+SamplingGridProperties& SamplingGridProperties::operator=(const SamplingGridProperties& other)
+{
+  // check for self-assignment
+  if (this != &other) {
+    for(size_t index=0; index<3; ++index) {
+      begin[index] = other.begin[index];
+      end[index] = other.end[index];
+    }
+    level = other.level;
+  }
+  return *this;
+}
+
+bool SamplingGridProperties::operator==(const SamplingGridProperties &_props) const
+{
+  bool status = true;
+  for (size_t i=0; i<3; ++i) {
+    status &= begin[i] == _props.begin[i];
+    status &= end[i] == _props.end[i];
+  }
+  status &= level == _props.level;
+  return status;
+}
+
+// we need to explicitly instantiate the serialization functions as
+// its is only serialized through its base class FragmentJob
+BOOST_CLASS_EXPORT_IMPLEMENT(SamplingGridProperties)
Index: src/Fragmentation/Summation/SetValues/SamplingGridProperties.hpp
===================================================================
--- src/Fragmentation/Summation/SetValues/SamplingGridProperties.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/SamplingGridProperties.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,109 @@
+/*
+ * SamplingGridProperties.hpp
+ *
+ *  Created on: 25.07.2012
+ *      Author: heber
+ */
+
+#ifndef SAMPLINGGRIDPROPERTIES_HPP_
+#define SAMPLINGGRIDPROPERTIES_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "boost/serialization/export.hpp"
+#include "boost/serialization/array.hpp"
+
+/** This class stores a sample function on a three-dimensional grid.
+ *
+ */
+class SamplingGridProperties {
+public:
+  /** Constructor for class SamplingGridProperties.
+   *
+   * \param _begin offset of grid
+   * \param _end edge length per axis
+   * \param _level number of gridpoints as \f$2^{\mathrm{level}}\f$ per unit(!) length
+   */
+  SamplingGridProperties(
+      const double _begin[3],
+      const double _end[3],
+      const int _level);
+
+  /** Copy constructor for class SamplingGridProperties.
+   *
+   * \param _props instance to copy from
+   */
+  SamplingGridProperties(const SamplingGridProperties &_props);
+
+  /** Default constructor.
+   */
+  SamplingGridProperties();
+
+  virtual ~SamplingGridProperties();
+
+  /** Checks whether another instance is compatible with this one.
+   *
+   *  \note Compatibility implies only that both grids have same grid spacing.
+   *
+   * \param _props other properties to check against
+   * \return true - are compatible, false - else
+   */
+  bool isCompatible(const SamplingGridProperties &_props) const
+  {
+    return level == _props.level;
+  }
+
+  /** Assignment operator.
+   *
+   * \param other other instance to assign ourselves to
+   */
+  SamplingGridProperties& operator=(const SamplingGridProperties& other);
+
+  /** Equality operator for class SamplingGridProperties.
+   *
+   * \param _props other object to compare to
+   */
+  bool operator==(const SamplingGridProperties &_props) const;
+
+  /** Inequality operator for class SamplingGridProperties.
+   *
+   * \param _props other object to compare to
+   */
+  bool operator!=(const SamplingGridProperties &_props) const
+  {
+    return (!(*this == _props));
+  }
+
+public:
+  //!> offset of grid
+  double begin[3];
+  //!> size of grid, i.e. edge length per axis of domain
+  double end[3];
+  //!> level, i.e. \f$2^{\mathrm{level}}\f$ grid points per axis per unit(!) length
+  int level;
+
+private:
+
+  friend class boost::serialization::access;
+  // serialization
+  template <typename Archive>
+  void serialize(Archive& ar, const unsigned int version)
+  {
+    int i;
+    for (i=0; i<3; ++i)
+      ar & begin[i];
+    for (i=0; i<3; ++i)
+      ar & end[i];
+    ar & level;
+  }
+
+};
+
+// we need to give this class a unique key for serialization
+// its is only serialized through its base class FragmentJob
+BOOST_CLASS_EXPORT_KEY(SamplingGridProperties)
+
+#endif /* SAMPLINGGRIDPROPERTIES_HPP_ */
Index: src/Fragmentation/Summation/SetValues/SamplingGrid_inline.hpp
===================================================================
--- src/Fragmentation/Summation/SetValues/SamplingGrid_inline.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/SamplingGrid_inline.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,106 @@
+/*
+ * SamplingGrid_inline.hpp
+ *
+ *  Created on: Jan 27, 2013
+ *      Author: heber
+ */
+
+#ifndef SAMPLINGGRID_INLINE_HPP_
+#define SAMPLINGGRID_INLINE_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "SamplingGrid.hpp"
+
+#ifdef HAVE_INLINE
+inline
+#endif
+const double SamplingGrid::getVolume() const
+{
+  double volume = 1.;
+  for (size_t i=0;i<3;++i)
+    volume *= end[i]-begin[i];
+  return volume;
+}
+
+#ifdef HAVE_INLINE
+inline
+#endif
+const double SamplingGrid::getWindowVolume() const
+{
+  double volume = 1.;
+  for (size_t i=0;i<3;++i)
+    volume *= end_window[i]-begin_window[i];
+  return volume;
+}
+
+
+#ifdef HAVE_INLINE
+inline
+#endif
+const size_t SamplingGrid::getTotalGridPoints() const
+{
+  return pow(getGridPointsPerAxis(),3);
+}
+
+#ifdef HAVE_INLINE
+inline
+#endif
+const size_t SamplingGrid::getGridPointsPerAxis() const
+{
+  return pow(2, level);
+}
+
+#ifdef HAVE_INLINE
+inline
+#endif
+const double SamplingGrid::getDeltaPerAxis(const size_t axis) const
+{
+  return getTotalLengthPerAxis(axis)/(double)getGridPointsPerAxis();
+}
+
+#ifdef HAVE_INLINE
+inline
+#endif
+const double SamplingGrid::getTotalLengthPerAxis(const size_t axis) const
+{
+  return (end[axis] - begin[axis]);
+}
+
+#ifdef HAVE_INLINE
+inline
+#endif
+const double SamplingGrid::getWindowLengthPerAxis(const size_t axis) const
+{
+  return (end_window[axis] - begin_window[axis]);
+}
+
+#ifdef HAVE_INLINE
+inline
+#endif
+const size_t SamplingGrid::getWindowGridPoints() const
+{
+  return getWindowGridPointsPerAxis(0)
+      *getWindowGridPointsPerAxis(1)
+      *getWindowGridPointsPerAxis(2);
+}
+
+
+#ifdef HAVE_INLINE
+inline
+#endif
+void SamplingGrid::setDomainSize(
+    const double _begin[3],
+    const double _end[3])
+{
+  for (size_t index=0;index<3;++index) {
+    begin[index] = _begin[index];
+    end[index] = _end[index];
+  }
+}
+
+
+#endif /* SAMPLINGGRID_INLINE_HPP_ */
Index: src/Fragmentation/Summation/SetValues/unittests/FragmentUnitTest.cpp
===================================================================
--- src/Fragmentation/Summation/SetValues/unittests/FragmentUnitTest.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/unittests/FragmentUnitTest.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,509 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2012 University of Bonn. All rights reserved.
+ * 
+ *
+ *   This file is part of MoleCuilder.
+ *
+ *    MoleCuilder is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 2 of the License, or
+ *    (at your option) any later version.
+ *
+ *    MoleCuilder is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * FragmentUnitTest.cpp
+ *
+ *  Created on: Aug 09, 2012
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+// include headers that implement a archive in simple text format
+#include <boost/archive/text_oarchive.hpp>
+#include <boost/archive/text_iarchive.hpp>
+
+#include "FragmentUnitTest.hpp"
+
+#include <algorithm>
+#include <limits>
+
+#include <boost/assign.hpp>
+#include <boost/foreach.hpp>
+
+#include <sstream>
+
+#include "CodePatterns/Assert.hpp"
+
+#ifdef HAVE_TESTRUNNER
+#include "UnitTestMain.hpp"
+#endif /*HAVE_TESTRUNNER*/
+
+using namespace boost::assign;
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( FragmentTest );
+
+
+void FragmentTest::setUp()
+{
+  // failing asserts should be thrown
+  ASSERT_DO(Assert::Throw);
+
+  Fragment::position_t pos(3,0.);
+  positions += pos;
+  pos[0] = 1.;
+  positions += pos;
+  pos[1] = 1.;
+  positions += pos;
+  pos[2] = 1.;
+  positions += pos;
+  charges += 1., 2., 3., 4.;
+  CPPUNIT_ASSERT_EQUAL( (size_t)3, pos.size() );
+  CPPUNIT_ASSERT( positions.size() == charges.size() );
+
+  fragment = new Fragment(positions, charges);
+}
+
+
+void FragmentTest::tearDown()
+{
+  delete fragment;
+}
+
+/** UnitTest for isPositionEqual()
+ */
+void FragmentTest::isPositionEqual_Test()
+{
+  // same position
+  for (Fragment::positions_t::const_iterator iter = positions.begin();
+      iter != positions.end(); ++iter)
+    CPPUNIT_ASSERT( Fragment::isPositionEqual(*iter, *iter) );
+
+  // other position
+  Fragment::position_t unequalpos(3,2.);
+  for (Fragment::positions_t::const_iterator iter = positions.begin();
+      iter != positions.end(); ++iter)
+    CPPUNIT_ASSERT( !Fragment::isPositionEqual(*iter, unequalpos) );
+}
+
+
+/** UnitTest for containsNuclei()
+ */
+void FragmentTest::containsNuclei_Test()
+{
+  {
+    // tests present ones for containment
+    Fragment::nuclei_t validnuclei(positions.size());
+    std::transform(positions.begin(), positions.end(),
+        charges.begin(), validnuclei.begin(), Fragment::createNucleus);
+    BOOST_FOREACH( Fragment::nucleus_t nucleus, validnuclei) {
+      CPPUNIT_ASSERT( fragment->containsNuclei(nucleus) );
+    }
+  }
+
+  {
+    // test some others
+    Fragment::nuclei_t invalidnuclei;
+    Fragment::position_t pos(3, -1.);
+    invalidnuclei += Fragment::createNucleus(pos, 1.);
+    pos[0] = 0.;
+    invalidnuclei += Fragment::createNucleus(pos, 1.);
+    pos[1] = 0.;
+    invalidnuclei += Fragment::createNucleus(pos, 1.);
+    pos[2] = -std::numeric_limits<double>::epsilon()*1e+4;
+    invalidnuclei += Fragment::createNucleus(pos, 1.);
+    BOOST_FOREACH( Fragment::nucleus_t nucleus, invalidnuclei) {
+      CPPUNIT_ASSERT( !fragment->containsNuclei(nucleus) );
+    }
+  }
+}
+
+/** UnitTest for removeNuclei()
+ */
+void FragmentTest::removeNuclei_Test()
+{
+  {
+    // tests present ones for removal
+    size_t size = fragment->nuclei.size();
+    Fragment::nuclei_t validnuclei(positions.size());
+    std::transform(positions.begin(), positions.end(),
+        charges.begin(), validnuclei.begin(), Fragment::createNucleus);
+    BOOST_FOREACH( Fragment::nucleus_t nucleus, validnuclei) {
+      CPPUNIT_ASSERT_NO_THROW( fragment->removeNuclei(nucleus) );
+      CPPUNIT_ASSERT_EQUAL( --size, fragment->nuclei.size() );
+    }
+  }
+  {
+    // test some others
+    Fragment::nuclei_t invalidnuclei;
+    Fragment::position_t pos(3, -1.);
+    invalidnuclei += Fragment::createNucleus(pos, 1.);
+    pos[0] = 0;
+    invalidnuclei += Fragment::createNucleus(pos, 1.);
+    pos[1] = 0;
+    invalidnuclei += Fragment::createNucleus(pos, 1.);
+    pos[2] = -std::numeric_limits<double>::epsilon()*1e+4;
+    invalidnuclei += Fragment::createNucleus(pos, 1.);
+    BOOST_FOREACH( Fragment::nucleus_t nucleus, invalidnuclei) {
+      CPPUNIT_ASSERT_NO_THROW( fragment->removeNuclei(nucleus) );
+    }
+  }
+}
+
+/** UnitTest for operator==() for Fragment::nucleus_t
+ */
+void FragmentTest::equalityNucleus_Test()
+{
+  Fragment::nuclei_t validnuclei(positions.size());
+  std::transform(positions.begin(), positions.end(),
+      charges.begin(), validnuclei.begin(), Fragment::createNucleus);
+  {
+    // create some nuclei
+    BOOST_FOREACH( Fragment::nucleus_t nucleus, validnuclei) {
+      CPPUNIT_ASSERT( nucleus == nucleus );
+    }
+  }
+
+  {
+    // create nucleus at other position
+    Fragment::position_t pos(3, 2.);
+    Fragment::nucleus_t unequalposnucleus(pos, 1.);
+    BOOST_FOREACH( Fragment::nucleus_t nucleus, validnuclei) {
+      CPPUNIT_ASSERT( nucleus != unequalposnucleus );
+    }
+  }
+
+  {
+    // create nucleus with different charge
+    Fragment::position_t pos(3, 1.);
+    Fragment::nucleus_t unequalchargenucleus(pos, 5.);
+    BOOST_FOREACH( Fragment::nucleus_t nucleus, validnuclei) {
+      CPPUNIT_ASSERT( nucleus != unequalchargenucleus );
+    }
+  }
+}
+
+/** UnitTest for operator==()
+ */
+void FragmentTest::equality_Test()
+{
+  // create different fragment
+  Fragment::positions_t otherpositions;
+  Fragment::position_t otherpos(3, 2.);
+  otherpositions += otherpos;
+  otherpos[0] = 0.;
+  otherpositions += otherpos;
+  otherpos[1] = 0.;
+  otherpositions += otherpos;
+  Fragment::charges_t othercharges;
+  othercharges += 1., 2., 3.;
+  Fragment otherfragment(otherpositions, othercharges);
+
+  CPPUNIT_ASSERT( !(*fragment == otherfragment) );
+  CPPUNIT_ASSERT( *fragment != otherfragment );
+
+  // test against empty fragment
+  Fragment emptyfragment;
+  CPPUNIT_ASSERT( !(*fragment == emptyfragment) );
+  CPPUNIT_ASSERT( *fragment != emptyfragment );
+
+  // tests against themselves
+  CPPUNIT_ASSERT( *fragment == *fragment );
+  CPPUNIT_ASSERT( otherfragment == otherfragment );
+  CPPUNIT_ASSERT( emptyfragment == emptyfragment );
+}
+
+/** UnitTest for operator+=()
+ */
+void FragmentTest::assignment_Test()
+{
+  // create different fragment
+  Fragment::positions_t otherpositions;
+  Fragment::position_t otherpos(3, 2.);
+  otherpositions += otherpos;
+  otherpos[0] = 0.;
+  otherpositions += otherpos;
+  otherpos[1] = 0.;
+  otherpositions += otherpos;
+  Fragment::charges_t othercharges;
+  othercharges += 1., 2., 3.;
+  Fragment otherfragment(otherpositions, othercharges);
+
+  // check for inequality
+  CPPUNIT_ASSERT( otherfragment.nuclei.size() != fragment->nuclei.size() );
+  CPPUNIT_ASSERT( otherfragment != *fragment );
+
+  //assign
+  otherfragment = *fragment;
+
+  // check for equality
+  CPPUNIT_ASSERT( otherfragment.nuclei.size() == fragment->nuclei.size() );
+  CPPUNIT_ASSERT( otherfragment == *fragment );
+}
+
+/** UnitTest for operator+=()
+ */
+void FragmentTest::operatorPlusEqual_NonOverlapping_Test()
+{
+  {
+    // create non-overlapping set
+    Fragment::positions_t otherpositions;
+    Fragment::position_t otherpos(3, 2.);
+    otherpositions += otherpos;
+    otherpos[0] = 0.;
+    otherpositions += otherpos;
+    otherpos[1] = 0.;
+    otherpositions += otherpos;
+    Fragment::charges_t othercharges;
+    othercharges += 1., 2., 3.;
+    Fragment otherfragment(otherpositions, othercharges);
+    const size_t othersize = otherfragment.nuclei.size();
+    const size_t size = fragment->nuclei.size();
+    *fragment += otherfragment;
+    CPPUNIT_ASSERT_EQUAL( othersize, otherfragment.nuclei.size() );
+    CPPUNIT_ASSERT_EQUAL( size+othersize, fragment->nuclei.size() );
+    {
+      // tests all for containment
+      Fragment::nuclei_t validnuclei(positions.size()+otherpositions.size());
+      Fragment::nuclei_t::iterator iter =
+          std::transform(positions.begin(), positions.end(),
+          charges.begin(), validnuclei.begin(), Fragment::createNucleus);
+      std::transform(otherpositions.begin(), otherpositions.end(),
+          othercharges.begin(), iter, Fragment::createNucleus);
+      BOOST_FOREACH( Fragment::nucleus_t nucleus, validnuclei) {
+        CPPUNIT_ASSERT( fragment->containsNuclei(nucleus) );
+      }
+    }
+    {
+      // tests positions for no containment in otherfragment
+      Fragment::nuclei_t invalidnuclei(positions.size());
+      std::transform(positions.begin(), positions.end(),
+          charges.begin(), invalidnuclei.begin(), Fragment::createNucleus);
+      BOOST_FOREACH( Fragment::nucleus_t nucleus, invalidnuclei) {
+        CPPUNIT_ASSERT( !otherfragment.containsNuclei(nucleus) );
+      }
+    }
+    {
+      // tests otherpositions for containment in otherfragment
+      Fragment::nuclei_t validnuclei(otherpositions.size());
+      std::transform(otherpositions.begin(), otherpositions.end(),
+          othercharges.begin(), validnuclei.begin(), Fragment::createNucleus);
+      BOOST_FOREACH( Fragment::nucleus_t nucleus, validnuclei) {
+        CPPUNIT_ASSERT( otherfragment.containsNuclei(nucleus) );
+      }
+    }
+  }
+}
+
+/** UnitTest for operator+=()
+ */
+void FragmentTest::operatorPlusEqual_Test()
+{
+  {
+    // create overlapping set (first overlaps)
+    Fragment::positions_t otherpositions;
+    Fragment::position_t otherpos(3, 1.);
+    otherpositions += otherpos;
+    otherpos[0] = 2.;
+    otherpositions += otherpos;
+    otherpos[1] = 2.;
+    otherpositions += otherpos;
+    Fragment::charges_t othercharges;
+    othercharges += 1., 2., 3.;
+    Fragment otherfragment(otherpositions, othercharges);
+    const size_t othersize = otherfragment.nuclei.size();
+    const size_t size = fragment->nuclei.size();
+    *fragment += otherfragment;
+    CPPUNIT_ASSERT_EQUAL( othersize, otherfragment.nuclei.size() );
+    CPPUNIT_ASSERT_EQUAL( size+othersize-1, fragment->nuclei.size() ); // one for already present
+    {
+      // tests all for containment
+      Fragment::nuclei_t validnuclei(positions.size()+otherpositions.size());
+      Fragment::nuclei_t::iterator iter =
+          std::transform(positions.begin(), positions.end(),
+          charges.begin(), validnuclei.begin(), Fragment::createNucleus);
+      std::transform(otherpositions.begin(), otherpositions.end(),
+          othercharges.begin(), iter, Fragment::createNucleus);
+      BOOST_FOREACH( Fragment::nucleus_t nucleus, validnuclei) {
+        CPPUNIT_ASSERT( fragment->containsNuclei(nucleus) );
+      }
+    }
+    {
+      // tests positions for no containment in otherfragment (but last)
+      Fragment::nuclei_t invalidnuclei(positions.size()-1);
+      std::transform(positions.begin(), --positions.end(),
+          charges.begin(), invalidnuclei.begin(), Fragment::createNucleus);
+      BOOST_FOREACH( Fragment::nucleus_t nucleus, invalidnuclei) {
+        CPPUNIT_ASSERT( !otherfragment.containsNuclei(nucleus) );
+      }
+    }
+    {
+      // tests otherpositions for containment in otherfragment
+      Fragment::nuclei_t validnuclei(otherpositions.size());
+      std::transform(otherpositions.begin(), otherpositions.end(),
+          othercharges.begin(), validnuclei.begin(), Fragment::createNucleus);
+      BOOST_FOREACH( Fragment::nucleus_t nucleus, validnuclei) {
+        CPPUNIT_ASSERT( otherfragment.containsNuclei(nucleus) );
+      }
+    }
+  }
+}
+
+/** UnitTest for operator-=()
+ */
+void FragmentTest::operatorMinusEqual_NonOverlapping_Test()
+{
+  {
+    // create non-overlapping set
+    Fragment::positions_t otherpositions;
+    Fragment::position_t otherpos(3, 2.);
+    otherpositions += otherpos;
+    otherpos[0] = 0.;
+    otherpositions += otherpos;
+    otherpos[1] = 0.;
+    otherpositions += otherpos;
+    Fragment::charges_t othercharges;
+    othercharges += 1., 2., 3.;
+    Fragment otherfragment(otherpositions, othercharges);
+    const size_t othersize = otherfragment.nuclei.size();
+    const size_t size = fragment->nuclei.size();
+    *fragment -= otherfragment;
+    CPPUNIT_ASSERT_EQUAL( othersize, otherfragment.nuclei.size() );
+    CPPUNIT_ASSERT_EQUAL( size, fragment->nuclei.size() );
+    {
+      // tests positions for containment
+      Fragment::nuclei_t validnuclei(positions.size());
+      std::transform(positions.begin(), positions.end(),
+          charges.begin(), validnuclei.begin(), Fragment::createNucleus);
+      BOOST_FOREACH( Fragment::nucleus_t nucleus, validnuclei) {
+        CPPUNIT_ASSERT( fragment->containsNuclei(nucleus) );
+      }
+    }
+    {
+      // tests otherpositions for no containment
+      Fragment::nuclei_t invalidnuclei(otherpositions.size());
+      std::transform(otherpositions.begin(), otherpositions.end(),
+          othercharges.begin(), invalidnuclei.begin(), Fragment::createNucleus);
+      BOOST_FOREACH( Fragment::nucleus_t nucleus, invalidnuclei) {
+        CPPUNIT_ASSERT( !fragment->containsNuclei(nucleus) );
+      }
+    }
+    {
+      // tests positions for no containment in otherfragment
+      Fragment::nuclei_t invalidnuclei(positions.size());
+      std::transform(positions.begin(), positions.end(),
+          charges.begin(), invalidnuclei.begin(), Fragment::createNucleus);
+      BOOST_FOREACH( Fragment::nucleus_t nucleus, invalidnuclei) {
+        CPPUNIT_ASSERT( !otherfragment.containsNuclei(nucleus) );
+      }
+    }
+    {
+      // tests otherpositions for containment in otherfragment
+      Fragment::nuclei_t validnuclei(otherpositions.size());
+      std::transform(otherpositions.begin(), otherpositions.end(),
+          othercharges.begin(), validnuclei.begin(), Fragment::createNucleus);
+      BOOST_FOREACH( Fragment::nucleus_t nucleus, validnuclei) {
+        CPPUNIT_ASSERT( otherfragment.containsNuclei(nucleus) );
+      }
+    }
+  }
+}
+
+/** UnitTest for operator-=()
+ */
+void FragmentTest::operatorMinusEqual_Test()
+{
+  {
+    // create overlapping set (first overlaps although with different charge)
+    Fragment::positions_t otherpositions;
+    Fragment::position_t otherpos(3, 1.);
+    otherpositions += otherpos;
+    otherpos[0] = 2.;
+    otherpositions += otherpos;
+    otherpos[1] = 2.;
+    otherpositions += otherpos;
+    Fragment::charges_t othercharges;
+    othercharges += 1., 2., 3.;
+    Fragment otherfragment(otherpositions, othercharges);
+    const size_t othersize = otherfragment.nuclei.size();
+    const size_t size = fragment->nuclei.size();
+    CPPUNIT_ASSERT( Fragment::isPositionEqual(otherpositions[0],positions[3]) );
+    *fragment -= otherfragment;
+    CPPUNIT_ASSERT_EQUAL( othersize, otherfragment.nuclei.size() );
+    CPPUNIT_ASSERT_EQUAL( size-1, fragment->nuclei.size() ); // just one overlaps
+    {
+      // tests all but last for containment
+      Fragment::nuclei_t validnuclei(positions.size());
+      std::transform(positions.begin(), positions.end(),
+          charges.begin(), validnuclei.begin(), Fragment::createNucleus);
+      BOOST_FOREACH( Fragment::nucleus_t nucleus, validnuclei) {
+        if (Fragment::isPositionEqual(nucleus.first, otherpositions[0])) // only test position
+          CPPUNIT_ASSERT( !fragment->containsNuclei(nucleus) );
+        else
+          CPPUNIT_ASSERT( fragment->containsNuclei(nucleus) );
+      }
+    }
+    {
+      // tests positions for no containment in otherfragment
+      Fragment::nuclei_t invalidnuclei(positions.size()-1);
+      std::transform(positions.begin(), --positions.end(),
+          charges.begin(), invalidnuclei.begin(), Fragment::createNucleus);
+      BOOST_FOREACH( Fragment::nucleus_t nucleus, invalidnuclei) {
+        CPPUNIT_ASSERT( !otherfragment.containsNuclei(nucleus) );
+      }
+    }
+    {
+      // tests otherpositions for containment in otherfragment
+      Fragment::nuclei_t validnuclei(otherpositions.size());
+      std::transform(otherpositions.begin(), otherpositions.end(),
+          othercharges.begin(), validnuclei.begin(), Fragment::createNucleus);
+      BOOST_FOREACH( Fragment::nucleus_t nucleus, validnuclei) {
+        CPPUNIT_ASSERT( otherfragment.containsNuclei(nucleus) );
+      }
+    }
+  }
+}
+
+
+
+/** UnitTest for serialization
+ */
+void FragmentTest::serializeTest()
+{
+  // serialize
+  std::stringstream outputstream;
+  boost::archive::text_oarchive oa(outputstream);
+  oa << fragment;
+
+  // deserialize
+  Fragment *samefragment = NULL;
+  std::stringstream returnstream(outputstream.str());
+  boost::archive::text_iarchive ia(returnstream);
+  ia >> samefragment;
+
+  CPPUNIT_ASSERT( *fragment == *samefragment );
+
+  delete samefragment;
+}
Index: src/Fragmentation/Summation/SetValues/unittests/FragmentUnitTest.hpp
===================================================================
--- src/Fragmentation/Summation/SetValues/unittests/FragmentUnitTest.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/unittests/FragmentUnitTest.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,60 @@
+/*
+ * FragmentUnitTest.hpp
+ *
+ *  Created on: Aug 09, 2012
+ *      Author: heber
+ */
+
+#ifndef FRAGMENTUNITTEST_HPP_
+#define FRAGMENTUNITTEST_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "Fragmentation/Summation/SetValues/Fragment.hpp"
+
+/********************************************** Test classes **************************************/
+
+class FragmentTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( FragmentTest) ;
+    CPPUNIT_TEST ( isPositionEqual_Test );
+    CPPUNIT_TEST ( containsNuclei_Test );
+    CPPUNIT_TEST ( removeNuclei_Test );
+    CPPUNIT_TEST ( equalityNucleus_Test );
+    CPPUNIT_TEST ( equality_Test );
+    CPPUNIT_TEST ( assignment_Test );
+    CPPUNIT_TEST ( operatorPlusEqual_NonOverlapping_Test );
+    CPPUNIT_TEST ( operatorPlusEqual_Test );
+    CPPUNIT_TEST ( operatorMinusEqual_NonOverlapping_Test );
+    CPPUNIT_TEST ( operatorMinusEqual_Test );
+    CPPUNIT_TEST ( serializeTest );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+      void setUp();
+      void tearDown();
+      void isPositionEqual_Test();
+      void containsNuclei_Test();
+      void removeNuclei_Test();
+      void equalityNucleus_Test();
+      void equality_Test();
+      void assignment_Test();
+      void operatorPlusEqual_NonOverlapping_Test();
+      void operatorPlusEqual_Test();
+      void operatorMinusEqual_NonOverlapping_Test();
+      void operatorMinusEqual_Test();
+      void serializeTest();
+
+private:
+      Fragment *fragment;
+      Fragment::positions_t positions;
+      Fragment::charges_t charges;
+};
+
+#endif /* FRAGMENTUNITTEST_HPP_ */
Index: src/Fragmentation/Summation/SetValues/unittests/HistogramUnitTest.cpp
===================================================================
--- src/Fragmentation/Summation/SetValues/unittests/HistogramUnitTest.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/unittests/HistogramUnitTest.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,318 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2012 University of Bonn. All rights reserved.
+ * 
+ *
+ *   This file is part of MoleCuilder.
+ *
+ *    MoleCuilder is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 2 of the License, or
+ *    (at your option) any later version.
+ *
+ *    MoleCuilder is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * HistogramUnitTest.cpp
+ *
+ *  Created on: Jul 26, 2012
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include "Fragmentation/Summation/SetValues/Histogram.hpp"
+
+#include "HistogramUnitTest.hpp"
+
+#include <cmath>
+#include <limits>
+
+#include <boost/assign.hpp>
+
+#include "CodePatterns/Assert.hpp"
+
+#ifdef HAVE_TESTRUNNER
+#include "UnitTestMain.hpp"
+#endif /*HAVE_TESTRUNNER*/
+
+using namespace boost::assign;
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( HistogramTest );
+
+
+void HistogramTest::setUp()
+{
+  // failing asserts should be thrown
+  ASSERT_DO(Assert::Throw);
+
+  histogram = NULL;
+}
+
+
+void HistogramTest::tearDown()
+{
+  delete histogram;
+}
+
+/** UnitTest for checking internal state after given samples
+ */
+void HistogramTest::internalState_Test()
+{
+  // generate non-empty histogramgram
+  Histogram::samples_t sampled_values;
+  sampled_values += 1.,2.,3.,4.;
+  CPPUNIT_ASSERT_EQUAL( (size_t)4, sampled_values.size() );
+  histogram = new Histogram(sampled_values, 1., 1.);
+  CPPUNIT_ASSERT_EQUAL( (size_t)4+1, histogram->bins.size() );
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->bins.begin() );
+  CPPUNIT_ASSERT_EQUAL( 1., histogram->bins.begin()->second );
+  CPPUNIT_ASSERT_EQUAL( 0., (--histogram->bins.end())->second );
+}
+
+/** UnitTest for getLowerEndBin()
+ */
+void HistogramTest::getLowerEndBin_Test()
+{
+  // generate non-empty histogramgram
+  Histogram::samples_t sampled_values;
+  sampled_values += 1.,2.,3.,4.;
+  CPPUNIT_ASSERT_EQUAL( (size_t)4, sampled_values.size() );
+  histogram = new Histogram(sampled_values, 1., 1.);
+
+  // check values inside
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getLowerEndBin(1.) );
+  CPPUNIT_ASSERT_EQUAL( 1., histogram->getLowerEndBin(1.)->first );
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getLowerEndBin(2.) );
+  CPPUNIT_ASSERT_EQUAL( 2., histogram->getLowerEndBin(2.)->first );
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getLowerEndBin(3.) );
+  CPPUNIT_ASSERT_EQUAL( 3., histogram->getLowerEndBin(3.)->first );
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getLowerEndBin(4.) );
+  CPPUNIT_ASSERT_EQUAL( 4., histogram->getLowerEndBin(4.)->first );
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getLowerEndBin(1.5) );
+  CPPUNIT_ASSERT_EQUAL( 1., histogram->getLowerEndBin(1.5)->first );
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getLowerEndBin(2.5) );
+  CPPUNIT_ASSERT_EQUAL( 2., histogram->getLowerEndBin(2.5)->first );
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getLowerEndBin(3.5) );
+  CPPUNIT_ASSERT_EQUAL( 3., histogram->getLowerEndBin(3.5)->first );
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getLowerEndBin(4.5) );
+  CPPUNIT_ASSERT_EQUAL( 4., histogram->getLowerEndBin(4.5)->first );
+
+  // check values outside
+  CPPUNIT_ASSERT( histogram->bins.end() == histogram->getLowerEndBin(.9) );
+  CPPUNIT_ASSERT( histogram->bins.end() == histogram->getLowerEndBin(5.01) );
+  CPPUNIT_ASSERT( histogram->bins.end() == histogram->getLowerEndBin(5.) );
+}
+
+/** UnitTest for getHigherEndBin()
+ */
+void HistogramTest::getHigherEndBin_Test()
+{
+  // generate non-empty histogramgram
+  Histogram::samples_t sampled_values;
+  sampled_values += 1.,2.,3.,4.;
+  CPPUNIT_ASSERT_EQUAL( (size_t)4, sampled_values.size() );
+  histogram = new Histogram(sampled_values, 1., 1.);
+
+  // check values inside
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getHigherEndBin(.5) );
+  CPPUNIT_ASSERT_EQUAL( 1., histogram->getHigherEndBin(.5)->first );
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getHigherEndBin(1.) );
+  CPPUNIT_ASSERT_EQUAL( 2., histogram->getHigherEndBin(1.)->first );
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getHigherEndBin(2.) );
+  CPPUNIT_ASSERT_EQUAL( 3., histogram->getHigherEndBin(2.)->first );
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getHigherEndBin(3.) );
+  CPPUNIT_ASSERT_EQUAL( 4., histogram->getHigherEndBin(3.)->first );
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getHigherEndBin(4.) );
+  CPPUNIT_ASSERT_EQUAL( 5., histogram->getHigherEndBin(4.)->first );
+  CPPUNIT_ASSERT( histogram->bins.end() == histogram->getHigherEndBin(5.) );
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getHigherEndBin(1.5) );
+  CPPUNIT_ASSERT_EQUAL( 2., histogram->getHigherEndBin(1.5)->first );
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getHigherEndBin(2.5) );
+  CPPUNIT_ASSERT_EQUAL( 3., histogram->getHigherEndBin(2.5)->first );
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getHigherEndBin(3.5) );
+  CPPUNIT_ASSERT_EQUAL( 4., histogram->getHigherEndBin(3.5)->first );
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getHigherEndBin(4.5) );
+  CPPUNIT_ASSERT_EQUAL( 5., histogram->getHigherEndBin(4.5)->first );
+  CPPUNIT_ASSERT( histogram->bins.end() == histogram->getHigherEndBin(5.5) );
+
+  // check values outside
+  CPPUNIT_ASSERT( histogram->bins.end() != histogram->getHigherEndBin(-.1) );
+  CPPUNIT_ASSERT_EQUAL( 1., histogram->getHigherEndBin(-.1)->first );
+  CPPUNIT_ASSERT( histogram->bins.end() == histogram->getHigherEndBin(5.01) );
+  CPPUNIT_ASSERT( histogram->bins.end() == histogram->getHigherEndBin(5.) );
+}
+
+
+/** UnitTest for getLowerEnd()
+ */
+void HistogramTest::getLowerEnd_Test()
+{
+  // create non-empty histogram and test
+  Histogram::samples_t sampled_values;
+  sampled_values += 1.,2.,3.,4.;
+  histogram = new Histogram(sampled_values, 1., 1.);
+
+  // go through each bin and check against getLowerEnd()
+  for (Histogram::Bins_t::const_iterator iter = histogram->bins.begin();
+      iter != histogram->bins.end(); ++iter)
+    CPPUNIT_ASSERT_EQUAL( iter->first, histogram->getLowerEnd(iter->first) );
+
+  CPPUNIT_ASSERT_EQUAL( 7., histogram->getLowerEnd(7.2) );
+  CPPUNIT_ASSERT_EQUAL( -5., histogram->getLowerEnd(-4.5) );
+  CPPUNIT_ASSERT_EQUAL( 4000., histogram->getLowerEnd(4000.1) );
+  CPPUNIT_ASSERT_EQUAL( -4001., histogram->getLowerEnd(-4000.9) );
+}
+
+/** UnitTest for isEmpty()
+ */
+void HistogramTest::isEmpty_Test()
+{
+  histogram = new Histogram(1., 1.);
+  // test on empty histogramgram
+  CPPUNIT_ASSERT( histogram->isEmpty() );
+
+  // create non-empty histogram and test
+  Histogram::samples_t sampled_values;
+  sampled_values += 1.,2.,3.,4.;
+  delete histogram;
+  histogram = new Histogram(sampled_values, 1., 1.);
+  CPPUNIT_ASSERT( !histogram->isEmpty() );
+
+}
+
+/** UnitTest for area()
+ */
+void HistogramTest::areaTest()
+{
+  histogram = new Histogram(1., 1.);
+  // test on empty histogramgram
+  CPPUNIT_ASSERT_EQUAL( 0., histogram->area() );
+
+  // create non-empty histogram and sum up
+  Histogram::samples_t sampled_values;
+  sampled_values += 1.,2.,3.,4., 4.;
+  delete histogram;
+  histogram = new Histogram(sampled_values, 1., 1.);
+  CPPUNIT_ASSERT_EQUAL( 5., histogram->area() );
+  Histogram::samples_t more_values;
+  more_values += 1.75,2.5,3.;
+  Histogram otherhistogram(more_values, 1.75, .625);
+  CPPUNIT_ASSERT_EQUAL( 3., otherhistogram.area() );
+
+}
+
+/** UnitTest for superposeOtherHistogram()
+ */
+void HistogramTest::superposeOtherHistogram_Test()
+{
+  // create two histograms, one is larger
+  Histogram::samples_t sampled_values;
+  sampled_values += 1.,2.,3.,4.;
+  histogram = new Histogram(sampled_values, 1., 1.);
+  Histogram::samples_t more_values;
+  more_values += 1.75,2.5,3.;
+  Histogram otherhistogram(more_values, 1.75, .625);
+
+  // check that size increases
+  CPPUNIT_ASSERT_EQUAL( (size_t)3+1, otherhistogram.bins.size() );
+  otherhistogram += *histogram;
+  const size_t number = (size_t)((
+      otherhistogram.getLowerEnd(4.+histogram->binwidth)
+      + otherhistogram.binwidth
+      - otherhistogram.getLowerEnd(1.))/otherhistogram.binwidth)+1;
+  CPPUNIT_ASSERT_EQUAL( number, otherhistogram.bins.size() );
+}
+
+/** UnitTest for operator+=()
+ */
+void HistogramTest::operatorPlusEqual_Test()
+{
+  // create non-empty histograms and sum them
+  Histogram::samples_t sampled_values;
+  sampled_values += 1.,2.,3.,4.;
+  histogram = new Histogram(sampled_values, 1., 1.);
+
+  // add upon larger
+  {
+    Histogram::samples_t more_values;
+    more_values += 0.5, 1.25, 2.5,3., 4.2;
+    Histogram otherhistogram(more_values, 0.5, .625);
+    const double areas = histogram->area() + otherhistogram.area();
+
+    // check that sum is now sum of both
+    otherhistogram += *histogram;
+//    CPPUNIT_ASSERT_EQUAL( areas, otherhistogram.area() );
+    CPPUNIT_ASSERT( fabs(areas - otherhistogram.area()) < std::numeric_limits<double>::epsilon()*1e+1 );
+  }
+
+  // add upon smaller
+  {
+    Histogram::samples_t more_values;
+    more_values += 1.75,2.5,3.;
+    Histogram otherhistogram(more_values, 1.75, .625);
+    const double areas = histogram->area() + otherhistogram.area();
+
+    // check that sum is now sum of both
+    otherhistogram += *histogram;
+//    CPPUNIT_ASSERT_EQUAL( areas, otherhistogram.area() );
+    CPPUNIT_ASSERT( fabs(areas - otherhistogram.area()) < std::numeric_limits<double>::epsilon()*1e+1 );
+  }
+}
+
+/** UnitTest for operator-=()
+ */
+void HistogramTest::operatorMinusEqual_Test()
+{
+  // create non-empty histograms and sum them
+  Histogram::samples_t sampled_values;
+  sampled_values += 1.,2.,3.,4.;
+  histogram = new Histogram(sampled_values, 1., 1.);
+
+  // subtract smaller
+  {
+    Histogram::samples_t more_values;
+    more_values += 0.5, 1.25, 2.5,3., 4.2;
+    Histogram otherhistogram(more_values, 0.5, .625);
+    const double difference = otherhistogram.area() - histogram->area();
+
+    // check that sum is now difference of both
+    otherhistogram -= *histogram;
+//    CPPUNIT_ASSERT_EQUAL( difference, otherhistogram.area() );
+    CPPUNIT_ASSERT( fabs(difference - otherhistogram.area()) < std::numeric_limits<double>::epsilon()*1e+1 );
+  }
+
+  // subtract larger
+  {
+    Histogram::samples_t more_values;
+    more_values += 1.75,2.5,3.;
+    Histogram otherhistogram(more_values, 1.75, .625);
+    const double difference = otherhistogram.area() - histogram->area();
+
+    // check that sum is now difference of both
+    otherhistogram -= *histogram;
+//    CPPUNIT_ASSERT_EQUAL( difference, otherhistogram.area() );
+    CPPUNIT_ASSERT( fabs(difference - otherhistogram.area()) < std::numeric_limits<double>::epsilon()*1e+1 );
+  }
+}
+
Index: src/Fragmentation/Summation/SetValues/unittests/HistogramUnitTest.hpp
===================================================================
--- src/Fragmentation/Summation/SetValues/unittests/HistogramUnitTest.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/unittests/HistogramUnitTest.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,55 @@
+/*
+ * HistogramUnitTest.hpp
+ *
+ *  Created on: Jul 26, 2012
+ *      Author: heber
+ */
+
+#ifndef HISTOGRAMUNITTEST_HPP_
+#define HISTOGRAMUNITTEST_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <cppunit/extensions/HelperMacros.h>
+
+class Histogram;
+
+/********************************************** Test classes **************************************/
+
+class HistogramTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( HistogramTest) ;
+    CPPUNIT_TEST ( internalState_Test );
+    CPPUNIT_TEST ( getLowerEndBin_Test );
+    CPPUNIT_TEST ( getHigherEndBin_Test );
+    CPPUNIT_TEST ( getLowerEnd_Test );
+    CPPUNIT_TEST ( isEmpty_Test );
+    CPPUNIT_TEST ( areaTest );
+    CPPUNIT_TEST ( superposeOtherHistogram_Test );
+    CPPUNIT_TEST ( operatorPlusEqual_Test );
+    CPPUNIT_TEST ( operatorMinusEqual_Test );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+      void setUp();
+      void tearDown();
+      void internalState_Test();
+      void getLowerEndBin_Test();
+      void getHigherEndBin_Test();
+      void getLowerEnd_Test();
+      void isEmpty_Test();
+      void areaTest();
+      void superposeOtherHistogram_Test();
+      void operatorPlusEqual_Test();
+      void operatorMinusEqual_Test();
+
+private:
+      Histogram *histogram;
+
+};
+
+#endif /* HISTOGRAMUNITTEST_HPP_ */
Index: src/Fragmentation/Summation/SetValues/unittests/IndexedVectorsUnitTest.cpp
===================================================================
--- src/Fragmentation/Summation/SetValues/unittests/IndexedVectorsUnitTest.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/unittests/IndexedVectorsUnitTest.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,183 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2012 University of Bonn. All rights reserved.
+ * 
+ *
+ *   This file is part of MoleCuilder.
+ *
+ *    MoleCuilder is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 2 of the License, or
+ *    (at your option) any later version.
+ *
+ *    MoleCuilder is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * IndexedVectorsUnitTest.cpp
+ *
+ *  Created on: Jul 29, 2012
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include "IndexedVectorsUnitTest.hpp"
+
+#include <cmath>
+#include <limits>
+
+#include <boost/assign.hpp>
+
+
+#include "CodePatterns/Assert.hpp"
+
+#ifdef HAVE_TESTRUNNER
+#include "UnitTestMain.hpp"
+#endif /*HAVE_TESTRUNNER*/
+
+using namespace boost::assign;
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( IndexedVectorsTest );
+
+
+void IndexedVectorsTest::setUp()
+{
+  // failing asserts should be thrown
+  ASSERT_DO(Assert::Throw);
+
+  // create two vector_t
+  ones += 1.,1.,1.;
+  twos += 2.,2.,2.;
+  threes += 3.,3.,3.;
+  CPPUNIT_ASSERT_EQUAL( IndexedVectors::FixedSize, ones.size() );
+  CPPUNIT_ASSERT_EQUAL( IndexedVectors::FixedSize, twos.size() );
+  CPPUNIT_ASSERT_EQUAL( IndexedVectors::FixedSize, threes.size() );
+
+  // create vectors
+
+  IndexedVectors::vectors_t vectors;
+  vectors.push_back(ones);
+  vectors.push_back(twos);
+  IndexedVectors::vectors_t othervectors;
+  othervectors.push_back(threes);
+  othervectors.push_back(threes);
+
+  // create two indices
+  IndexedVectors::indices_t indices;
+  indices += 1,2;
+  IndexedVectors::indices_t otherindices;
+  otherindices += 1,3;
+
+  // create indexed vectors
+  ivectors = new IndexedVectors(indices, vectors);
+  otherivectors = new IndexedVectors(otherindices, othervectors);
+}
+
+
+void IndexedVectorsTest::tearDown()
+{
+  delete ivectors;
+  delete otherivectors;
+}
+
+
+/** UnitTest for cstor's
+ */
+void IndexedVectorsTest::Constructor_Test()
+{
+  // check whether -1 is dropped
+  IndexedVectors::indices_t indices;
+  indices += 1,-1,3;
+  IndexedVectors::vectors_t vectors;
+  vectors.push_back(ones);
+  vectors.push_back(twos);
+  vectors.push_back(threes);
+  IndexedVectors testivectors(indices, vectors);
+
+  CPPUNIT_ASSERT_EQUAL( (size_t)2, testivectors.vectors.size() );
+  CPPUNIT_ASSERT_EQUAL( ones, testivectors.vectors[1] );
+  CPPUNIT_ASSERT_EQUAL( threes, testivectors.vectors[3] );
+  CPPUNIT_ASSERT( testivectors.vectors.find(-1) == testivectors.vectors.end() );
+}
+
+/** UnitTest for operator+=()
+ */
+void IndexedVectorsTest::operatorPlusEqual_Test()
+{
+  // safeguard initial sizes
+  CPPUNIT_ASSERT_EQUAL( (size_t)2, ivectors->vectors.size() );
+  CPPUNIT_ASSERT_EQUAL( (size_t)2, otherivectors->vectors.size() );
+
+  // perform operation
+  *ivectors += *otherivectors;
+
+  // check new and ole sizes
+  CPPUNIT_ASSERT_EQUAL( (size_t)3, ivectors->vectors.size() );
+  CPPUNIT_ASSERT_EQUAL( (size_t)2, otherivectors->vectors.size() );
+
+  // then check result
+  IndexedVectors::vector_t result( IndexedVectors::nullvector );
+  CPPUNIT_ASSERT_EQUAL( IndexedVectors::FixedSize, result.size() );
+  for (size_t i=0; i<IndexedVectors::FixedSize; ++i)
+    result[i] = ones[i] + threes[i];
+  for (IndexedVectors::indexedvectors_t::const_iterator iter = ivectors->vectors.begin();
+      iter != ivectors->vectors.end(); ++iter) {
+    CPPUNIT_ASSERT_EQUAL( IndexedVectors::FixedSize, iter->second.size() );
+  }
+  CPPUNIT_ASSERT_EQUAL( result, ivectors->vectors[1] );
+  CPPUNIT_ASSERT_EQUAL( twos, ivectors->vectors[2] );
+  CPPUNIT_ASSERT_EQUAL( threes, ivectors->vectors[3] );
+}
+
+/** UnitTest for operator-=()
+ */
+void IndexedVectorsTest::operatorMinusEqual_Test()
+{
+  // safeguard initial sizes
+  CPPUNIT_ASSERT_EQUAL( (size_t)2, ivectors->vectors.size() );
+  CPPUNIT_ASSERT_EQUAL( (size_t)2, otherivectors->vectors.size() );
+
+  // perform operation
+  *ivectors -= *otherivectors;
+
+  // check new and ole sizes
+  CPPUNIT_ASSERT_EQUAL( (size_t)3, ivectors->vectors.size() );
+  CPPUNIT_ASSERT_EQUAL( (size_t)2, otherivectors->vectors.size() );
+
+  // then check result
+  IndexedVectors::vector_t result( IndexedVectors::nullvector );
+  IndexedVectors::vector_t thirdresult( IndexedVectors::nullvector );
+  CPPUNIT_ASSERT_EQUAL( IndexedVectors::FixedSize, result.size() );
+  for (size_t i=0; i<IndexedVectors::FixedSize; ++i) {
+    result[i] = ones[i] - threes[i];
+    thirdresult[i] = -threes[i];
+  }
+  for (IndexedVectors::indexedvectors_t::const_iterator iter = ivectors->vectors.begin();
+      iter != ivectors->vectors.end(); ++iter) {
+    CPPUNIT_ASSERT_EQUAL( IndexedVectors::FixedSize, iter->second.size() );
+  }
+  CPPUNIT_ASSERT_EQUAL( result, ivectors->vectors[1] );
+  CPPUNIT_ASSERT_EQUAL( twos, ivectors->vectors[2] );
+  CPPUNIT_ASSERT_EQUAL( thirdresult, ivectors->vectors[3] );
+}
+
Index: src/Fragmentation/Summation/SetValues/unittests/IndexedVectorsUnitTest.hpp
===================================================================
--- src/Fragmentation/Summation/SetValues/unittests/IndexedVectorsUnitTest.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/unittests/IndexedVectorsUnitTest.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,47 @@
+/*
+ * IndexedVectorsUnitTest.hpp
+ *
+ *  Created on: Jul 29, 2012
+ *      Author: heber
+ */
+
+#ifndef INDEXEDVECTORSUNITTEST_HPP_
+#define INDEXEDVECTORSUNITTEST_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "Fragmentation/Summation/SetValues/IndexedVectors.hpp"
+
+/********************************************** Test classes **************************************/
+
+class IndexedVectorsTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( IndexedVectorsTest) ;
+    CPPUNIT_TEST ( Constructor_Test );
+    CPPUNIT_TEST ( operatorPlusEqual_Test );
+    CPPUNIT_TEST ( operatorMinusEqual_Test );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+      void setUp();
+      void tearDown();
+      void Constructor_Test();
+      void operatorPlusEqual_Test();
+      void operatorMinusEqual_Test();
+
+private:
+      IndexedVectors *ivectors;
+      IndexedVectors *otherivectors;
+      IndexedVectors::vector_t ones;
+      IndexedVectors::vector_t twos;
+      IndexedVectors::vector_t threes;
+
+};
+
+#endif /* INDEXEDVECTORSUNITTEST_HPP_ */
Index: src/Fragmentation/Summation/SetValues/unittests/Makefile.am
===================================================================
--- src/Fragmentation/Summation/SetValues/unittests/Makefile.am	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/unittests/Makefile.am	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,61 @@
+# PLEASE adhere to the alphabetical ordering in this Makefile!
+# Also indentation by a single tab
+
+FRAGMENTATIONSETVALUESTESTSSOURCES = \
+	../Fragmentation/Summation/SetValues/unittests/FragmentUnitTest.cpp \
+	../Fragmentation/Summation/SetValues/unittests/HistogramUnitTest.cpp \
+	../Fragmentation/Summation/SetValues/unittests/IndexedVectorsUnitTest.cpp \
+	../Fragmentation/Summation/SetValues/unittests/SamplingGridPropertiesUnitTest.cpp \
+	../Fragmentation/Summation/SetValues/unittests/SamplingGridUnitTest.cpp
+
+FRAGMENTATIONSETVALUESTESTSHEADERS = \
+	../Fragmentation/Summation/SetValues/unittests/FragmentUnitTest.hpp \
+	../Fragmentation/Summation/SetValues/unittests/HistogramUnitTest.hpp \
+	../Fragmentation/Summation/SetValues/unittests/IndexedVectorsUnitTest.hpp \
+	../Fragmentation/Summation/SetValues/unittests/SamplingGridPropertiesUnitTest.hpp \
+	../Fragmentation/Summation/SetValues/unittests/SamplingGridUnitTest.hpp
+
+FRAGMENTATIONSETVALUESTESTS = \
+	FragmentUnitTest \
+	HistogramUnitTest \
+	IndexedVectorsUnitTest \
+	SamplingGridPropertiesUnitTest \
+	SamplingGridUnitTest
+
+TESTS += $(FRAGMENTATIONSETVALUESTESTS)
+check_PROGRAMS += $(FRAGMENTATIONSETVALUESTESTS) 
+noinst_PROGRAMS += $(FRAGMENTATIONSETVALUESTESTS)
+
+FRAGMENTATIONSETVALUESLIBS = \
+	../libMolecuilderFragmentationSetValues.la \
+	${CodePatterns_LIBS} \
+	$(BOOST_LIB)
+
+FragmentUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
+	../Fragmentation/Summation/SetValues/unittests/FragmentUnitTest.cpp \
+	../Fragmentation/Summation/SetValues/unittests/FragmentUnitTest.hpp
+FragmentUnitTest_LDADD = \
+	${FRAGMENTATIONSETVALUESLIBS} \
+	$(BOOST_SERIALIZATION_LDFLAGS) $(BOOST_SERIALIZATION_LIBS)
+
+HistogramUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
+	../Fragmentation/Summation/SetValues/unittests/HistogramUnitTest.cpp \
+	../Fragmentation/Summation/SetValues/unittests/HistogramUnitTest.hpp
+HistogramUnitTest_LDADD = ${FRAGMENTATIONSETVALUESLIBS}
+
+IndexedVectorsUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
+	../Fragmentation/Summation/SetValues/unittests/IndexedVectorsUnitTest.cpp \
+	../Fragmentation/Summation/SetValues/unittests/IndexedVectorsUnitTest.hpp
+IndexedVectorsUnitTest_LDADD = ${FRAGMENTATIONSETVALUESLIBS}
+
+SamplingGridPropertiesUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
+	../Fragmentation/Summation/SetValues/unittests/SamplingGridPropertiesUnitTest.cpp \
+	../Fragmentation/Summation/SetValues/unittests/SamplingGridPropertiesUnitTest.hpp
+SamplingGridPropertiesUnitTest_LDADD = ${FRAGMENTATIONSETVALUESLIBS}
+
+SamplingGridUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
+	../Fragmentation/Summation/SetValues/unittests/SamplingGridUnitTest.cpp \
+	../Fragmentation/Summation/SetValues/unittests/SamplingGridUnitTest.hpp
+SamplingGridUnitTest_LDADD = ${FRAGMENTATIONSETVALUESLIBS}
+
+#AUTOMAKE_OPTIONS = parallel-tests
Index: src/Fragmentation/Summation/SetValues/unittests/SamplingGridPropertiesUnitTest.cpp
===================================================================
--- src/Fragmentation/Summation/SetValues/unittests/SamplingGridPropertiesUnitTest.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/unittests/SamplingGridPropertiesUnitTest.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,125 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2012 University of Bonn. All rights reserved.
+ * 
+ *
+ *   This file is part of MoleCuilder.
+ *
+ *    MoleCuilder is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 2 of the License, or
+ *    (at your option) any later version.
+ *
+ *    MoleCuilder is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * SamplingGridPropertiesUnitTest.cpp
+ *
+ *  Created on: Jul 29, 2012
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include "SamplingGridPropertiesUnitTest.hpp"
+
+#include "CodePatterns/Assert.hpp"
+
+#ifdef HAVE_TESTRUNNER
+#include "UnitTestMain.hpp"
+#endif /*HAVE_TESTRUNNER*/
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( SamplingGridPropertiesTest );
+
+
+void SamplingGridPropertiesTest::setUp()
+{
+  // failing asserts should be thrown
+  ASSERT_DO(Assert::Throw);
+
+  // create the grid
+  const double begin[3] = { 0., 0., 0. };
+  const double end[3] = { 1., 1., 1. };
+  props = new SamplingGridProperties(begin, end, 2);
+}
+
+
+void SamplingGridPropertiesTest::tearDown()
+{
+  delete props;
+}
+
+/** UnitTest for isCompatible()
+ */
+void SamplingGridPropertiesTest::isCompatible_Test()
+{
+  const double begin[3] = { 0., 0., 0. };
+  const double otherbegin[3] = { 1., 0.1, -0.5 };
+  const double end[3] = { 1., 1., 1. };
+  const double otherend[3] = { 2., 2., 2. };
+
+  // create other props
+  SamplingGridProperties sameprops(begin, end, 2);
+  SamplingGridProperties otherprops(otherbegin, end, 2);
+  SamplingGridProperties anotherprops(begin, otherend, 2);
+  SamplingGridProperties moreotherprops(begin, end, 4);
+
+  // only checks for same level
+  CPPUNIT_ASSERT( props->isCompatible(*props) );
+  CPPUNIT_ASSERT( props->isCompatible(sameprops) );
+  CPPUNIT_ASSERT( sameprops.isCompatible(*props) );
+  CPPUNIT_ASSERT( props->isCompatible(otherprops) );
+  CPPUNIT_ASSERT( otherprops.isCompatible(*props) );
+  CPPUNIT_ASSERT( props->isCompatible(anotherprops) );
+  CPPUNIT_ASSERT( anotherprops.isCompatible(*props) );
+  CPPUNIT_ASSERT( !props->isCompatible(moreotherprops) );
+  CPPUNIT_ASSERT( !moreotherprops.isCompatible(*props) );
+}
+
+/** UnitTest for operator==()
+ */
+void SamplingGridPropertiesTest::equality_Test()
+{
+  const double begin[3] = { 0., 0., 0. };
+  const double otherbegin[3] = { 1., 0.1, -0.5 };
+  const double end[3] = { 1., 1., 1. };
+  const double otherend[3] = { 2., 2., 2. };
+
+  // create other props
+  SamplingGridProperties sameprops(begin, end, 2);
+  SamplingGridProperties otherprops(otherbegin, end, 2);
+  SamplingGridProperties anotherprops(begin, otherend, 2);
+  SamplingGridProperties moreotherprops(begin, end, 4);
+
+  // check for same level, begin, and end
+  CPPUNIT_ASSERT( *props == *props );
+  CPPUNIT_ASSERT( *props == sameprops );
+  CPPUNIT_ASSERT( sameprops == *props);
+  CPPUNIT_ASSERT( *props != otherprops);
+  CPPUNIT_ASSERT( otherprops != *props);
+  CPPUNIT_ASSERT( *props != anotherprops );
+  CPPUNIT_ASSERT( anotherprops != *props );
+  CPPUNIT_ASSERT( *props != moreotherprops );
+  CPPUNIT_ASSERT( moreotherprops != *props );
+}
+
Index: src/Fragmentation/Summation/SetValues/unittests/SamplingGridPropertiesUnitTest.hpp
===================================================================
--- src/Fragmentation/Summation/SetValues/unittests/SamplingGridPropertiesUnitTest.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/unittests/SamplingGridPropertiesUnitTest.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,41 @@
+/*
+ * SamplingGridPropertiesUnitTest.hpp
+ *
+ *  Created on: Jul 29, 2012
+ *      Author: heber
+ */
+
+#ifndef SAMPLINGGRIDPROPERTIESUNITTEST_HPP_
+#define SAMPLINGGRIDPROPERTIESUNITTEST_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "Fragmentation/Summation/SetValues/SamplingGridProperties.hpp"
+
+/********************************************** Test classes **************************************/
+
+class SamplingGridPropertiesTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( SamplingGridPropertiesTest) ;
+    CPPUNIT_TEST ( isCompatible_Test );
+    CPPUNIT_TEST ( equality_Test );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+      void setUp();
+      void tearDown();
+      void isCompatible_Test();
+      void equality_Test();
+
+private:
+      SamplingGridProperties *props;
+
+};
+
+#endif /* SAMPLINGGRIDPROPERTIESUNITTEST_HPP_ */
Index: src/Fragmentation/Summation/SetValues/unittests/SamplingGridUnitTest.cpp
===================================================================
--- src/Fragmentation/Summation/SetValues/unittests/SamplingGridUnitTest.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/unittests/SamplingGridUnitTest.cpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,422 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2012 University of Bonn. All rights reserved.
+ * 
+ *
+ *   This file is part of MoleCuilder.
+ *
+ *    MoleCuilder is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 2 of the License, or
+ *    (at your option) any later version.
+ *
+ *    MoleCuilder is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * SamplingGridUnitTest.cpp
+ *
+ *  Created on: Jul 29, 2012
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include "SamplingGridUnitTest.hpp"
+
+#include "CodePatterns/Assert.hpp"
+
+#include <boost/assign.hpp>
+#include <cmath>
+#include <numeric>
+
+#ifdef HAVE_TESTRUNNER
+#include "UnitTestMain.hpp"
+#endif /*HAVE_TESTRUNNER*/
+
+using namespace boost::assign;
+
+/********************************************** Test classes **************************************/
+
+const double grid_value=1.;
+
+#define NUMBEROFSAMPLES(n) (size_t)(pow(pow(2,n),3))
+#define DOMAINVOLUME(l) (size_t)pow(l,3)
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( SamplingGridTest );
+
+
+void SamplingGridTest::setUp()
+{
+  // failing asserts should be thrown
+  ASSERT_DO(Assert::Throw);
+
+  // create the grid
+  const double begin[3] = { 0., 0., 0. };
+  const double end[3] = { 1., 1., 1. };
+  for (size_t i=0; i< DOMAINVOLUME(1)*NUMBEROFSAMPLES(2); ++i)
+    values += grid_value;
+  grid = new SamplingGrid(begin, end, 2, values);
+  CPPUNIT_ASSERT_EQUAL( grid_value, *(grid->sampled_grid.begin()) );
+}
+
+
+void SamplingGridTest::tearDown()
+{
+  delete grid;
+}
+
+/** UnitTest on compatible combination of props and values
+ */
+void SamplingGridTest::compatibleGrids_Test()
+{
+  // check illegal grid
+  const double begin[3] = { 0., 0., 0. };
+  const double end[3] = { 2., 2., 2. };
+  SamplingGridProperties illegal_props(begin, end, 5);
+  SamplingGridProperties legal_props(begin, end, 2);
+  CPPUNIT_ASSERT( !grid->isCompatible(illegal_props) );
+  CPPUNIT_ASSERT( grid->isCompatible(legal_props) );
+  SamplingGrid::sampledvalues_t illegal_values;
+  for (size_t i=0; i< NUMBEROFSAMPLES(1); ++i)
+    illegal_values += 1.5;
+  SamplingGrid::sampledvalues_t legal_values;
+  for (size_t i=0; i< NUMBEROFSAMPLES(2); ++i)
+    legal_values += 1.5;
+#ifndef NDEBUG
+  // throws because props and size of illegal_values don't match
+  std::cout << "The following assertion is intended and does not indicate a failure of the test." << std::endl;
+  CPPUNIT_ASSERT_THROW( SamplingGrid illegal_grid(illegal_props, illegal_values), Assert::AssertionFailure );
+  std::cout << "The following assertion is intended and does not indicate a failure of the test." << std::endl;
+  CPPUNIT_ASSERT_THROW( SamplingGrid illegal_grid(legal_props, illegal_values), Assert::AssertionFailure );
+  std::cout << "The following assertion is intended and does not indicate a failure of the test." << std::endl;
+  CPPUNIT_ASSERT_THROW( SamplingGrid illegal_grid(illegal_props, legal_values), Assert::AssertionFailure );
+#endif
+  // check that grid is still the same
+  for (SamplingGrid::sampledvalues_t::const_iterator iter = grid->sampled_grid.begin();
+      iter != grid->sampled_grid.end(); ++iter)
+    CPPUNIT_ASSERT_EQUAL( grid_value, *iter );
+}
+
+/** UnitTest for isCongruent()
+ */
+void SamplingGridTest::isCongruent_Test()
+{
+  const double begin[3] = { 0., 0., 0. };
+  const double end[3] = { 2., 2., 2. };
+  const double otherbegin[3] = { 0.1, 0.1, 0.1 };
+  const double otherend[3] = { 1., 1., 1. };
+  SamplingGridProperties illegal_begin_props(otherbegin, end, 5);
+  SamplingGridProperties illegal_end_props(begin, otherend, 5);
+  SamplingGridProperties illegal_level_props(begin, end, 5);
+  SamplingGridProperties legal_props(begin, end, 2);
+
+  // differing windows
+  const double begin_window[3] = { 0.5, 0.5, 0.5 };
+  const double end_window[3] = { 1., 1., 1. };
+  const double otherbegin_window[3] = { 0.45, 0.45, 0.45 };
+  const double otherend_window[3] = { 1.05, 1.05, 1.05 };
+
+  // check that incompatible grid are also incongruent
+  SamplingGrid default_grid(legal_props);
+  // note that we always construct a temporary SamplingGrid from given props
+  CPPUNIT_ASSERT( default_grid.isCompatible(illegal_begin_props) == default_grid.isCongruent(illegal_begin_props));
+  CPPUNIT_ASSERT( default_grid.isCompatible(illegal_end_props) == default_grid.isCongruent(illegal_end_props));
+  CPPUNIT_ASSERT( default_grid.isCompatible(illegal_level_props) == default_grid.isCongruent(illegal_level_props));
+  CPPUNIT_ASSERT( default_grid.isCompatible(legal_props) == default_grid.isCongruent(legal_props) );
+
+  default_grid.setWindowSize(begin, end);
+  // same window
+  {
+    SamplingGrid illegal_begin_grid(illegal_begin_props);
+    SamplingGrid illegal_end_grid(illegal_end_props);
+    SamplingGrid illegal_level_grid(illegal_level_props);
+    SamplingGrid legal_grid(legal_props);
+    illegal_begin_grid.setWindowSize(begin, end);
+    illegal_end_grid.setWindowSize(begin, end);
+    illegal_level_grid.setWindowSize(begin, end);
+    legal_grid.setWindowSize(begin, end);
+    CPPUNIT_ASSERT( !illegal_begin_grid.isCongruent(default_grid) );
+    CPPUNIT_ASSERT( !illegal_end_grid.isCongruent(default_grid) );
+    CPPUNIT_ASSERT( !illegal_level_grid.isCongruent(default_grid) );
+    CPPUNIT_ASSERT( legal_grid.isCongruent(default_grid) );
+  }
+
+  // different begin
+  {
+    SamplingGrid illegal_begin_grid(illegal_begin_props);
+    SamplingGrid illegal_end_grid(illegal_end_props);
+    SamplingGrid illegal_level_grid(illegal_level_props);
+    SamplingGrid legal_grid(legal_props);
+    illegal_begin_grid.setWindowSize(otherbegin_window, end);
+    illegal_end_grid.setWindowSize(otherbegin_window, end);
+    illegal_level_grid.setWindowSize(otherbegin_window, end);
+    legal_grid.setWindowSize(begin, end);
+    CPPUNIT_ASSERT( !illegal_begin_grid.isCongruent(legal_grid) );
+    CPPUNIT_ASSERT( !illegal_end_grid.isCongruent(legal_grid) );
+    CPPUNIT_ASSERT( !illegal_level_grid.isCongruent(legal_grid) );
+    CPPUNIT_ASSERT( legal_grid.isCongruent(default_grid) );
+  }
+
+  // different end
+  {
+    SamplingGrid illegal_begin_grid(illegal_begin_props);
+    SamplingGrid illegal_end_grid(illegal_end_props);
+    SamplingGrid illegal_level_grid(illegal_level_props);
+    SamplingGrid legal_grid(legal_props);
+    illegal_begin_grid.setWindowSize(begin, otherend_window);
+    illegal_end_grid.setWindowSize(begin, otherend_window);
+    illegal_level_grid.setWindowSize(begin, otherend_window);
+    legal_grid.setWindowSize(begin, end);
+    CPPUNIT_ASSERT( !illegal_begin_grid.isCongruent(legal_grid) );
+    CPPUNIT_ASSERT( !illegal_end_grid.isCongruent(legal_grid) );
+    CPPUNIT_ASSERT( !illegal_level_grid.isCongruent(legal_grid) );
+    CPPUNIT_ASSERT( legal_grid.isCongruent(default_grid) );
+  }
+
+  // different begin and end
+  {
+    SamplingGrid illegal_begin_grid(illegal_begin_props);
+    SamplingGrid illegal_end_grid(illegal_end_props);
+    SamplingGrid illegal_level_grid(illegal_level_props);
+    SamplingGrid legal_grid(legal_props);
+    illegal_begin_grid.setWindowSize(otherbegin_window, otherend_window);
+    illegal_end_grid.setWindowSize(otherbegin_window, otherend_window);
+    illegal_level_grid.setWindowSize(otherbegin_window, otherend_window);
+    legal_grid.setWindowSize(begin, end);
+    CPPUNIT_ASSERT( !illegal_begin_grid.isCongruent(legal_grid) );
+    CPPUNIT_ASSERT( !illegal_end_grid.isCongruent(legal_grid) );
+    CPPUNIT_ASSERT( !illegal_level_grid.isCongruent(legal_grid) );
+    CPPUNIT_ASSERT( legal_grid.isCongruent(default_grid) );
+  }
+}
+
+/** UnitTest for integral()
+ */
+void SamplingGridTest::integral_Test()
+{
+  double sum = 0.;
+  sum = std::accumulate( grid->sampled_grid.begin(), grid->sampled_grid.end(), sum );
+  CPPUNIT_ASSERT_EQUAL( sum*grid->getVolume()/grid->getWindowGridPoints(), grid->integral() );
+}
+
+/** UnitTest for getVolume()
+ */
+void SamplingGridTest::getVolume_Test()
+{
+  CPPUNIT_ASSERT_EQUAL( 1., grid->getVolume() );
+}
+
+/** UnitTest for getWindowSize()
+ */
+void SamplingGridTest::getWindowSize_Test()
+{
+  // check size of default grid
+  CPPUNIT_ASSERT_EQUAL( (size_t)(NUMBEROFSAMPLES(grid->level)), grid->getWindowGridPoints() );
+
+  // create another one and check its size, too
+  const double begin[3] = { 0., 0., 0. };
+  const double end[3] = { 1., 1., 1. };
+  for (size_t level = 3; level<=6; ++level) {
+    values.clear();
+    for (size_t i=0; i< NUMBEROFSAMPLES(level); ++i)
+      values += grid_value;
+    delete grid;
+    // use other pointer in case something fails
+    SamplingGrid *tmpgrid = new SamplingGrid(begin, end, level, values);
+    grid = tmpgrid;
+    CPPUNIT_ASSERT_EQUAL( (size_t)NUMBEROFSAMPLES(level), grid->getWindowGridPoints() );
+  }
+}
+
+/** UnitTest for extendWindow()
+ */
+void SamplingGridTest::extendWindow_Test()
+{
+  // we have a grid with size of one, extend to twice the size and check
+  const double begin[3] = { 0., 0., 0. };
+  const double size = 2.;
+  const double end[3] = { size, size, size };
+  double offset[3];
+  for (offset[0] = 0.; offset[0] <= 1.; offset[0] += .5)
+    for (offset[1] = 0.; offset[1] <= 1.; offset[1] += .5)
+      for (offset[2] = 0.; offset[2] <= 1.; offset[2] += .5) {
+        const double window_begin[3] = { 0.+offset[0], 0.+offset[1], 0.+offset[2]};
+        const double window_end[3] = { 1.+offset[0], 1.+offset[1], 1.+offset[2]};
+        SamplingGrid newgrid(begin, end, 2);
+        newgrid.setWindowSize(window_begin, window_end);
+        // resize values by hand to new window size. Otherwise they get zero'd.
+        newgrid.sampled_grid = values;
+        newgrid.sampled_grid.resize(NUMBEROFSAMPLES(1));
+        newgrid.extendWindow(begin, end);
+
+        // check integral
+        CPPUNIT_ASSERT_EQUAL(
+            grid_value/(NUMBEROFSAMPLES(grid->level)/NUMBEROFSAMPLES(grid->level)),
+            grid->integral()
+            );
+
+        // check number of points
+        CPPUNIT_ASSERT_EQUAL(
+            (size_t)NUMBEROFSAMPLES(grid->level),
+            grid->getWindowGridPoints()
+            );
+      }
+}
+
+/** UnitTest for extendWindow() with asymmetric values
+ */
+void SamplingGridTest::extendWindow_asymmetric_Test()
+{
+  std::cout << "SamplingGridTest::extendWindow_asymmetric_Test()" << std::endl;
+  const double begin[3] = { 0., 0., 0. };
+  const double end[3] = { 2., 2., 2. };
+  double offset[3];
+  for (offset[0] = 0.; offset[0] <= 1.; offset[0] += .5)
+    for (offset[1] = 0.; offset[1] <= 1.; offset[1] += .5)
+      for (offset[2] = 0.; offset[2] <= 1.; offset[2] += .5) {
+        const double window_begin[3] = { 0.+offset[0], 0.+offset[1], 0.+offset[2]};
+        const double window_end[3] = { 1.+offset[0], 1.+offset[1], 1.+offset[2]};
+        SamplingGrid newgrid(begin, end, 2);
+        CPPUNIT_ASSERT_EQUAL( (size_t)0, newgrid.getWindowGridPoints() );
+        newgrid.setWindowSize(window_begin, window_end);
+        // window size is only half of domain size
+        const size_t max_samples = NUMBEROFSAMPLES(newgrid.level)*pow(0.5,3);
+        for (size_t i=0; i< max_samples; ++i)
+          newgrid.sampled_grid += grid_value*i;
+        const size_t sum_weight = (max_samples)*(max_samples-1)/2;
+        const double integral = newgrid.integral();
+        newgrid.extendWindow(begin, end);
+
+        // check that integral has remained the same
+        CPPUNIT_ASSERT_EQUAL( integral, newgrid.integral() );
+        CPPUNIT_ASSERT_EQUAL( grid_value*sum_weight/DOMAINVOLUME(2), newgrid.integral() );
+      }
+}
+
+/** UnitTest for addOntoWindow()
+ */
+void SamplingGridTest::addOntoWindow_Test()
+{
+  // first window is from (0,0,0) to (1,1,1)
+  CPPUNIT_ASSERT_EQUAL( 1.*grid_value, grid->integral() );
+
+  // create values for half-sized window
+  values.clear();
+  for (size_t i=0; i< (size_t)pow(.5*pow(2,2),3); ++i)
+    values += grid_value;
+
+  // check that too large a window throws
+#ifndef NDEBUG
+  const double begin[3] = { .5, .5, .5 };
+  const double wrongend[3] = { 1.5, 1.5, 1.5 };
+  std::cout << "The following assertion is intended and does not indicate a failure of the test." << std::endl;
+  CPPUNIT_ASSERT_THROW( grid->addOntoWindow(begin, wrongend, values, +1.), Assert::AssertionFailure );
+#endif
+
+  // create another window from (.5,.5,.5) to (1., 1., 1.)
+  double offset[3];
+  for (offset[0] = 0.; offset[0] <= .5; offset[0] += .5)
+    for (offset[1] = 0.; offset[1] <= .5; offset[1] += .5)
+      for (offset[2] = 0.; offset[2] <= .5; offset[2] += .5) {
+        const double window_begin[3] = { 0.+offset[0], 0.+offset[1], 0.+offset[2]};
+        const double window_end[3] = { .5+offset[0], .5+offset[1], .5+offset[2]};
+
+        SamplingGrid newgrid(*grid);
+        // now perform working operation
+        newgrid.addOntoWindow(window_begin, window_end, values, +1.);
+
+        // check integral to be one and one eighth times the old value
+        CPPUNIT_ASSERT_EQUAL( (1.+pow(.5,3))*grid_value, newgrid.integral() );
+      }
+}
+
+/** UnitTest for addOntoWindow() with asymmetric values
+ */
+void SamplingGridTest::addOntoWindow_asymmetric_Test()
+{
+  const size_t size = grid->end[0]-grid->begin[0];
+  // check with asymmetric values
+  grid->sampled_grid.clear();
+  grid->sampled_grid.resize(DOMAINVOLUME(size)*NUMBEROFSAMPLES(grid->level), 0.);
+
+  for (size_t i=0;i<grid->level*(grid->end[0]-grid->begin[0]);++i)
+    grid->sampled_grid[(i*2+0)*2+0] += .5*grid_value;
+  for (size_t i=0;i<grid->level*(grid->end[1]-grid->begin[1]);++i)
+    grid->sampled_grid[(0*2+i)*2+0] += 1.*grid_value;
+  for (size_t i=0;i<grid->level*(grid->end[2]-grid->begin[2]);++i)
+    grid->sampled_grid[(0*2+0)*2+i] += 1.5*grid_value;
+
+  const double integral = grid->integral();
+
+  // now perform working operation
+  grid->addOntoWindow(grid->begin, grid->end, values, +1.);
+  // values is equal to integral of 1.
+  CPPUNIT_ASSERT_EQUAL( 1.+integral, grid->integral() );
+}
+
+/** UnitTest for operator+=()
+ */
+void SamplingGridTest::operatorPlusEqual_Test()
+{
+  // create other grid
+  const double begin[3] = { 0., 0., 0. };
+  const double end[3] = { 1., 1., 1. };
+  SamplingGrid::sampledvalues_t othervalues;
+  const double othergrid_value = 1.5;
+  for (size_t i=0; i< NUMBEROFSAMPLES(2); ++i)
+    othervalues += othergrid_value;
+  SamplingGrid othergrid(begin, end, 2, othervalues);
+  CPPUNIT_ASSERT_EQUAL( othergrid_value, *(othergrid.sampled_grid.begin()) );
+
+  // perform operation
+  CPPUNIT_ASSERT_NO_THROW( *grid += othergrid );
+
+  // check the contents of the grid
+  const double sum = grid_value+othergrid_value;
+  for (SamplingGrid::sampledvalues_t::const_iterator iter = grid->sampled_grid.begin();
+      iter != grid->sampled_grid.end(); ++iter)
+    CPPUNIT_ASSERT_EQUAL( sum, *iter );
+}
+
+/** UnitTest for operator-=()
+ */
+void SamplingGridTest::operatorMinusEqual_Test()
+{
+  // create other grid
+  const double begin[3] = { 0., 0., 0. };
+  const double end[3] = { 1., 1., 1. };
+  SamplingGrid::sampledvalues_t othervalues;
+  const double othergrid_value = 1.5;
+  for (size_t i=0; i< NUMBEROFSAMPLES(2); ++i)
+    othervalues += othergrid_value;
+  SamplingGrid othergrid(begin, end, 2, othervalues);
+  CPPUNIT_ASSERT_EQUAL( othergrid_value, *(othergrid.sampled_grid.begin()) );
+
+  // perform operation
+  CPPUNIT_ASSERT_NO_THROW( *grid -= othergrid );
+
+  // check the contents of the grid
+  const double difference = grid_value-othergrid_value;
+  for (SamplingGrid::sampledvalues_t::const_iterator iter = grid->sampled_grid.begin();
+      iter != grid->sampled_grid.end(); ++iter)
+    CPPUNIT_ASSERT_EQUAL( difference, *iter );
+}
+
Index: src/Fragmentation/Summation/SetValues/unittests/SamplingGridUnitTest.hpp
===================================================================
--- src/Fragmentation/Summation/SetValues/unittests/SamplingGridUnitTest.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
+++ src/Fragmentation/Summation/SetValues/unittests/SamplingGridUnitTest.hpp	(revision fbf143b9fb3d409937862a091c997ab0428c8b3e)
@@ -0,0 +1,59 @@
+/*
+ * SamplingGridUnitTest.hpp
+ *
+ *  Created on: Jul 29, 2012
+ *      Author: heber
+ */
+
+#ifndef SAMPLINGGRIDUNITTEST_HPP_
+#define SAMPLINGGRIDUNITTEST_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "Fragmentation/Summation/SetValues/SamplingGrid.hpp"
+
+/********************************************** Test classes **************************************/
+
+class SamplingGridTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( SamplingGridTest) ;
+    CPPUNIT_TEST ( compatibleGrids_Test );
+    CPPUNIT_TEST ( operatorPlusEqual_Test );
+    CPPUNIT_TEST ( operatorMinusEqual_Test );
+    CPPUNIT_TEST ( isCongruent_Test );
+    CPPUNIT_TEST ( integral_Test );
+    CPPUNIT_TEST ( getVolume_Test );
+    CPPUNIT_TEST ( getWindowSize_Test );
+    CPPUNIT_TEST ( extendWindow_Test );
+    CPPUNIT_TEST ( extendWindow_asymmetric_Test );
+    CPPUNIT_TEST ( addOntoWindow_Test );
+    CPPUNIT_TEST ( addOntoWindow_asymmetric_Test );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+      void setUp();
+      void tearDown();
+      void compatibleGrids_Test();
+      void isCongruent_Test();
+      void integral_Test();
+      void getVolume_Test();
+      void getWindowSize_Test();
+      void extendWindow_Test();
+      void extendWindow_asymmetric_Test();
+      void addOntoWindow_Test();
+      void addOntoWindow_asymmetric_Test();
+      void operatorPlusEqual_Test();
+      void operatorMinusEqual_Test();
+
+private:
+      SamplingGrid *grid;
+      SamplingGrid::sampledvalues_t values;
+};
+
+#endif /* SAMPLINGGRIDUNITTEST_HPP_ */
