Index: src/Actions/FragmentationAction/FragmentationAutomationAction.cpp
===================================================================
--- src/Actions/FragmentationAction/FragmentationAutomationAction.cpp	(revision 2388688c3e24cb54fc978d6c8efd1ddead979d8a)
+++ src/Actions/FragmentationAction/FragmentationAutomationAction.cpp	(revision a2295a3134eea526357331db860e7c4fa6cbf28c)
@@ -247,5 +247,5 @@
     const std::map<JobId_t,MPQCData> &fragmentData,
     const std::string &KeySetFilename,
-    SamplingGrid &full_sample,
+    std::vector<SamplingGrid> &full_sample,
     Fragment &full_fragment)
 {
@@ -286,5 +286,10 @@
         fragmentData, MatrixNrLookup, container, subsetmap));
   // obtain full grid
-  full_sample = boost::fusion::at_key<MPQCDataFused::sampled_grid>(Result_Grid_fused.back());
+  full_sample.clear();
+  full_sample.reserve(Result_Grid_fused.size());
+  for (std::vector<MPQCDataGridMap_t>::const_iterator iter = ++Result_Grid_fused.begin();
+      iter !=Result_Grid_fused.end();
+      ++iter)
+    full_sample.push_back(boost::fusion::at_key<MPQCDataFused::sampled_grid>((*iter)));
   full_fragment = boost::fusion::at_key<MPQCDataFused::fragment>(Result_Fragment_fused.back());
 
@@ -459,5 +464,5 @@
  * @param fragmentData MPQCData resulting from the jobs
  * @param longrangeData VMGData resulting from long-range jobs
- * @param fullsolutionData VMGData resulting from long-range of full problem
+ * @param fullsolutionData VMGData resulting from long-range of full problem from level 2 onward
  * @param KeySetFilename filename with keysets to associate forces correctly
  * @param NoAtoms total number of atoms
@@ -467,8 +472,8 @@
     const std::map<JobId_t,MPQCData> &fragmentData,
     const std::map<JobId_t,VMGData> &longrangeData,
-    const VMGData &fullsolutionData,
+    const std::vector<VMGData> &fullsolutionData,
     const std::string &KeySetFilename,
     size_t NoAtoms,
-    SamplingGrid &full_sample)
+    std::vector<SamplingGrid> &full_sample)
 {
   // create lookup from job nr to fragment number
@@ -550,11 +555,46 @@
                 VMGMatrixNrLookup);
     epotentialSummer(subsetmap->getMaximumSubsetLevel());
-    SamplingGrid full_sample_solution = fullsolutionData.sampled_potential;
-//    LOG(0, "Remaining long-range energy from energy_potential is " << full_sample_solution.integral()-epotentialSummer.getFullContribution() << ".");
-    full_sample_solution -= potentialSummer.getFullContribution();
-    // multiply element-wise with charge distribution
-    LOG(0, "Remaining long-range potential integral is " << full_sample_solution.integral() << ".");
-    full_sample_solution *= full_sample;
-    LOG(0, "Remaining long-range energy from potential integral is " << full_sample_solution.integral() << ".");
+
+    std::vector<VMGDataLongRangeMap_t> Result_LongRange_fused;
+    Result_LongRange_fused.reserve(subsetmap->getMaximumSubsetLevel());
+    for (size_t level = 1; level <= subsetmap->getMaximumSubsetLevel(); ++level) {
+      // weight times correct charge density of the same level
+      // NOTE: potential for level 1 is not calculated as saturation hydrogen
+      // are not removed on this level yet
+      const size_t potentiallevel = level < 2 ? 0 : (level-2);
+      SamplingGrid charge_weight = boost::fusion::at_key<MPQCDataFused::sampled_grid>(Result_Grid_fused[level-1]);
+      SamplingGrid full_sample_solution = fullsolutionData[potentiallevel].sampled_potential;
+      SamplingGrid short_range_correction = potentialSummer(level);
+  //    LOG(0, "Remaining long-range energy from energy_potential is " << full_sample_solution.integral()-epotentialSummer.getFullContribution() << ".");
+      full_sample_solution -= short_range_correction;
+      // multiply element-wise with charge distribution
+      VMGDataLongRangeMap_t instance;
+      boost::fusion::at_key<VMGDataFused::potential_longrange>(instance) = full_sample_solution.integral();
+      LOG(0, "Remaining long-range potential integral of level " << level << " is "
+          << full_sample_solution.integral() << ".");
+      boost::fusion::at_key<VMGDataFused::potential_shortrange>(instance) = short_range_correction.integral();
+      LOG(0, "Short-range correction potential integral of level " << level << " is "
+          << short_range_correction.integral() << ".");
+      boost::fusion::at_key<VMGDataFused::energy_longrange>(instance) = full_sample_solution.integral(charge_weight);
+      LOG(0, "Remaining long-range energy from potential integral of level " << level << " is "
+          << full_sample_solution.integral(charge_weight) << ".");
+      boost::fusion::at_key<VMGDataFused::energy_shortrange>(instance) = short_range_correction.integral(charge_weight);
+      LOG(0, "Short-range correction energy from potential integral of level " << level << " is "
+          << short_range_correction.integral(charge_weight) << ".");
+      Result_LongRange_fused.push_back(instance);
+    }
+    {
+      //    LOG(0, "Remaining long-range energy from energy_potential is " << full_sample_solution.integral()-epotentialSummer.getFullContribution() << ".");
+      SamplingGrid full_sample_solution = fullsolutionData.back().sampled_potential;
+      SamplingGrid short_range_correction = potentialSummer.getFullContribution();
+      full_sample_solution -= short_range_correction;
+      // multiply element-wise with charge distribution
+      LOG(0, "Remaining long-range potential integral is " << full_sample_solution.integral() << ".");
+      LOG(0, "Short-range correction potential integral of level is " << short_range_correction.integral() << ".");
+      LOG(0, "Remaining long-range energy from potential integral is "
+          << full_sample_solution.integral(full_sample.back()) << ".");
+      LOG(0, "Short-range correction energy from potential integral is "
+          << short_range_correction.integral(full_sample.back()) << ".");
+    }
 
     // TODO: Extract long-range corrections to forces
@@ -567,5 +607,5 @@
                 VMGMatrixNrLookup);
     elongSummer(subsetmap->getMaximumSubsetLevel());
-    double e_long = fullsolutionData.e_long;
+    double e_long = fullsolutionData.back().e_long;
     e_long -= elongSummer.getFullContribution();
     LOG(0, "Remaining long-range energy is " << e_long << ".");
@@ -577,7 +617,13 @@
             Result_Energy_fused, MaxLevel);
     LOG(0, "Energy table is \n" << energyresult);
+
+    const std::string gridresult =
+        writeTable<VMGDataLongRangeMap_t, VMGDataLongRangeVector_t >()(
+            Result_LongRange_fused, MaxLevel);
+    LOG(0, "LongRange table is \n" << gridresult);
+
     const std::string eigenvalueresult;
-
     LOG(0, "Eigenvalue table is \n" << eigenvalueresult);
+
     const std::string forceresult =
         writeTable<MPQCDataForceMap_t, MPQCDataForceVector_t>()(
@@ -635,5 +681,5 @@
   // obtain combined charge density
   LOG(1, "INFO: Parsing fragment files from " << params.path.get() << ".");
-  SamplingGrid full_sample;
+  std::vector<SamplingGrid> full_sample; // have charges from level 2 onward summed up
   Fragment full_fragment;
   sumUpChargeDensity(
@@ -644,5 +690,5 @@
 
   // Phase Four: obtain more ids
-  vmgcontroller.requestIds(fragmentData.size()+1);
+  vmgcontroller.requestIds(fragmentData.size()+full_sample.size());
 
   // Phase Five: create VMGJobs
@@ -652,14 +698,21 @@
 
   // Phase Six: calculate result
-  vmgcontroller.waitforResults(fragmentData.size()+1);
+  vmgcontroller.waitforResults(fragmentData.size()+full_sample.size());
   std::map<JobId_t, VMGData> longrangeData;
   vmgcontroller.getResults(longrangeData);
-  ASSERT( fragmentData.size()+1 == longrangeData.size(),
-      "FragmentationFragmentationAutomationAction::performCall() - number of MPQCresults+1 "
-      +toString(fragmentData.size()+1)+" and VMGresults "+toString(longrangeData.size())+" don't match.");
-
-  // remove full solution from map (must be highest id), has to be treated extra
-  VMGData fullsolutionData = (--longrangeData.end())->second;
-  longrangeData.erase(--longrangeData.end());
+  ASSERT( fragmentData.size()+full_sample.size() == longrangeData.size(),
+      "FragmentationFragmentationAutomationAction::performCall() - number of MPQCresults+"
+      +toString(full_sample.size())+" "+toString(fragmentData.size()+full_sample.size())
+      +" and VMGresults "+toString(longrangeData.size())+" don't match.");
+
+  // remove full solution corresponding to full_sample from map (must be highest ids), has to be treated extra
+  std::map<JobId_t, VMGData>::iterator iter = longrangeData.end();
+  for (size_t i=0;i<full_sample.size();++i)
+    --iter;
+  std::map<JobId_t, VMGData>::iterator remove_iter = iter;
+  std::vector<VMGData> fullsolutionData;
+  for (; iter != longrangeData.end(); ++iter)
+    fullsolutionData.push_back(iter->second);
+  longrangeData.erase(remove_iter, longrangeData.end());
 
   // Final phase: print result
@@ -674,11 +727,13 @@
         full_sample);
 
-    // create debug jobs to print the summed-up potential to vtk files
-    debugcontroller.requestIds(1);
-    if (!debugcontroller.createDebugJobs(std::vector<SamplingGrid>(1,full_sample)))
-      return Action::failure;
-    debugcontroller.waitforResults(1);
-    std::map<JobId_t, std::string> debugData;
-    debugcontroller.getResults(debugData);
+    if (!full_sample.empty()) {
+      // create debug jobs to print the summed-up potential to vtk files
+      debugcontroller.requestIds(full_sample.size());
+      if (!debugcontroller.createDebugJobs(full_sample))
+        return Action::failure;
+      debugcontroller.waitforResults(full_sample.size());
+      std::map<JobId_t, std::string> debugData;
+      debugcontroller.getResults(debugData);
+    }
   }
   }
Index: src/Fragmentation/Automation/VMGFragmentController.cpp
===================================================================
--- src/Fragmentation/Automation/VMGFragmentController.cpp	(revision 2388688c3e24cb54fc978d6c8efd1ddead979d8a)
+++ src/Fragmentation/Automation/VMGFragmentController.cpp	(revision a2295a3134eea526357331db860e7c4fa6cbf28c)
@@ -50,5 +50,5 @@
 bool VMGFragmentController::createLongRangeJobs(
     const std::map<JobId_t, MPQCData> &fragmentData,
-    const SamplingGrid &full_sampled_grid,
+    const std::vector<SamplingGrid> &full_sampled_grid,
     const Fragment &full_fragment,
     const size_t near_field_cells)
@@ -68,9 +68,10 @@
   }
 
+  // prepare positions and charges of full system
+  std::vector< std::vector<double> > positions;
+  std::vector<double> charges;
   {
     const World::AtomComposite &atoms = World::getInstance().getAllAtoms();
-    std::vector< std::vector<double> > positions;
     positions.reserve(atoms.size());
-    std::vector<double> charges;
     charges.reserve(atoms.size());
     std::vector<double> position(3, 0.);
@@ -83,9 +84,13 @@
       charges.push_back((double)((*iter)->getElement().getAtomicNumber()));
     }
+  }
+  for(std::vector<SamplingGrid>::const_iterator iter = full_sampled_grid.begin();
+      iter != full_sampled_grid.end();
+      ++iter) {
     const JobId_t next_id = getAvailableId();
-    LOG(1, "INFO: Creating full VMGJob with " << full_sampled_grid
+    LOG(1, "INFO: Creating full VMGJob with " << *iter
         << " gridpoints and " << charges.size() << " particle charges.");
     FragmentJob::ptr testJob(
-        new VMGJob(next_id, full_sampled_grid, positions, charges, near_field_cells) );
+        new VMGJob(next_id, *iter, positions, charges, near_field_cells) );
     jobs.push_back(testJob);
   }
Index: src/Fragmentation/Automation/VMGFragmentController.hpp
===================================================================
--- src/Fragmentation/Automation/VMGFragmentController.hpp	(revision 2388688c3e24cb54fc978d6c8efd1ddead979d8a)
+++ src/Fragmentation/Automation/VMGFragmentController.hpp	(revision a2295a3134eea526357331db860e7c4fa6cbf28c)
@@ -16,4 +16,7 @@
 
 #include "Fragmentation/Automation/SpecificFragmentController.hpp"
+
+#include <map>
+#include <vector>
 
 #include "JobMarket/JobId.hpp"
@@ -38,5 +41,5 @@
   bool createLongRangeJobs(
       const std::map<JobId_t, MPQCData> &fragmentData,
-      const SamplingGrid &full_sampled_grid,
+      const std::vector<SamplingGrid> &full_sampled_grid,
       const Fragment &full_fragment,
       const size_t near_field_cells);
Index: src/Fragmentation/Summation/OrthogonalFullSummator.hpp
===================================================================
--- src/Fragmentation/Summation/OrthogonalFullSummator.hpp	(revision 2388688c3e24cb54fc978d6c8efd1ddead979d8a)
+++ src/Fragmentation/Summation/OrthogonalFullSummator.hpp	(revision a2295a3134eea526357331db860e7c4fa6cbf28c)
@@ -84,4 +84,5 @@
   MapValue getFullContribution() const {
     typename SetValueMap<MapValue>::const_iterator iter = --OS.getSetValues().end();
+    LOG(2, "DEBUG: Returning Full contribution of level " << iter->first->size() << ".");
     return iter->second->getContribution();
   }
Index: src/Jobs/VMGDataFused.hpp
===================================================================
--- src/Jobs/VMGDataFused.hpp	(revision 2388688c3e24cb54fc978d6c8efd1ddead979d8a)
+++ src/Jobs/VMGDataFused.hpp	(revision a2295a3134eea526357331db860e7c4fa6cbf28c)
@@ -28,4 +28,11 @@
   struct energy_potential {};
   struct energy_long {};
+
+  // keys for longrange
+
+  struct potential_longrange {};
+  struct potential_shortrange {};
+  struct energy_longrange {};
+  struct energy_shortrange {};
 }
 
Index: src/Jobs/VMGDataMap.hpp
===================================================================
--- src/Jobs/VMGDataMap.hpp	(revision 2388688c3e24cb54fc978d6c8efd1ddead979d8a)
+++ src/Jobs/VMGDataMap.hpp	(revision a2295a3134eea526357331db860e7c4fa6cbf28c)
@@ -41,4 +41,25 @@
 > VMGDataVector_t;
 
+/** This boost::fusion map defines key-value or rather key-type pairs with
+ * which we associate all sampled_grid data members in VMGData and their type.
+ *
+ * This lets us resolves any ambiguitites of types in VMGData, e.g.
+ * to know vector<double> is forces or energy_eigenvalues.
+ *
+ */
+typedef boost::fusion::map<
+    boost::fusion::pair<VMGDataFused::potential_longrange, double >,
+    boost::fusion::pair<VMGDataFused::potential_shortrange, double >,
+    boost::fusion::pair<VMGDataFused::energy_longrange, double >,
+    boost::fusion::pair<VMGDataFused::energy_shortrange, double >
+> VMGDataLongRangeMap_t;
+
+typedef boost::mpl::list<
+    VMGDataFused::potential_longrange,
+    VMGDataFused::potential_shortrange,
+    VMGDataFused::energy_longrange,
+    VMGDataFused::energy_shortrange
+> VMGDataLongRangeVector_t;
+
 
 #endif /* VMGDATAMAP_HPP_ */
Index: src/Jobs/VMGData_printKeyNames.hpp
===================================================================
--- src/Jobs/VMGData_printKeyNames.hpp	(revision 2388688c3e24cb54fc978d6c8efd1ddead979d8a)
+++ src/Jobs/VMGData_printKeyNames.hpp	(revision a2295a3134eea526357331db860e7c4fa6cbf28c)
@@ -31,5 +31,10 @@
 #define tokensequence \
 	(sampled_potential) \
-  (energy_long)
+  (energy_potential) \
+  (energy_long) \
+  (potential_longrange) \
+  (potential_shortrange) \
+  (energy_longrange) \
+  (energy_shortrange)
 
 /// we take note of the number of keys in tokensequence as (local) loop bounds below
