Index: src/Fragmentation/Automation/VMGFragmentController.cpp
===================================================================
--- src/Fragmentation/Automation/VMGFragmentController.cpp	(revision 5281ff99623f1f5b7ed8ddd755f2d6daf60cc077)
+++ src/Fragmentation/Automation/VMGFragmentController.cpp	(revision e2925fd7f1c9a7ab18fe71b5a40598b65ddb05dc)
@@ -71,4 +71,6 @@
     const size_t near_field_cells,
     const size_t interpolation_degree,
+    const SampleParticles_t _SampleParticles,
+    const TreatGrid_t _TreatGrid,
     const MPQCData::DoValenceOnly_t _DoValenceOnly,
     const bool _DoPrintDebug)
@@ -86,9 +88,12 @@
         new VMGJob(
             next_id,
-            data.sampled_grid,
+            _TreatGrid == DoTreatGrid ?
+                data.sampled_grid :
+                SamplingGrid(data.sampled_grid.begin, data.sampled_grid.end, data.sampled_grid.level),
             data.positions,
             data.charges,
             near_field_cells,
             interpolation_degree,
+            _SampleParticles == DoSampleParticles,
             _DoPrintDebug) );
     jobs.push_back(testJob);
@@ -122,4 +127,5 @@
       iter != full_sampled_grid.end();
       ++iter) {
+    const SamplingGrid &grid = *iter;
     const JobId_t next_id = getAvailableId();
     LOG(1, "INFO: Creating full VMGJob with " << *iter
@@ -128,9 +134,12 @@
         new VMGJob(
             next_id,
-            *iter,
+            _TreatGrid == DoTreatGrid ?
+                grid :
+                SamplingGrid(grid.begin, grid.end, grid.level),
             positions,
             charges,
             near_field_cells,
             interpolation_degree,
+            _SampleParticles == DoSampleParticles,
             _DoPrintDebug) );
     jobs.push_back(testJob);
Index: src/Fragmentation/Automation/VMGFragmentController.hpp
===================================================================
--- src/Fragmentation/Automation/VMGFragmentController.hpp	(revision 5281ff99623f1f5b7ed8ddd755f2d6daf60cc077)
+++ src/Fragmentation/Automation/VMGFragmentController.hpp	(revision e2925fd7f1c9a7ab18fe71b5a40598b65ddb05dc)
@@ -39,4 +39,16 @@
   {}
 
+  //>! Enumerization of possible treatment of nuclei point charges.
+  enum SampleParticles_t {
+    DontSampleParticles=0,
+    DoSampleParticles=1
+  };
+
+  //>! Enumerization of possible treatment of electronic charge distribution on grid.
+  enum TreatGrid_t {
+    DontTreatGrid=0,
+    DoTreatGrid=1
+  };
+
   /** Creates the long range jobs from finished short-range jobs.
    *
@@ -45,4 +57,6 @@
    * \param near_field_cells radius of b-splines for interpolation charges onto grid
    * \param interpolation_degree interpolation degree for newton interpolation from grid
+   * \param _SampleParticles whether to add nuclei charges to vmg calculation.
+   * \param _TreatGrid whether to add electron charge distribution to vmg calculation.
    * \param _DoValenceOnly whether to neglect core electron and nuclei charges
    * \param _DoPrintDebug whether we do print grid for debug visualization or not
@@ -53,4 +67,6 @@
       const size_t near_field_cells,
       const size_t interpolation_degree,
+      const SampleParticles_t _SampleParticles,
+      const TreatGrid_t _TreatGrid,
       const MPQCData::DoValenceOnly_t _DoValenceOnly,
       const bool _DoPrintDebug);
Index: src/Fragmentation/Summation/Containers/FragmentationLongRangeResults.cpp
===================================================================
--- src/Fragmentation/Summation/Containers/FragmentationLongRangeResults.cpp	(revision 5281ff99623f1f5b7ed8ddd755f2d6daf60cc077)
+++ src/Fragmentation/Summation/Containers/FragmentationLongRangeResults.cpp	(revision e2925fd7f1c9a7ab18fe71b5a40598b65ddb05dc)
@@ -129,4 +129,6 @@
       boost::fusion::at_key<VMGDataFused::electron_longrange>(instance) = 0.;
       boost::fusion::at_key<VMGDataFused::electron_shortrange>(instance) = 0.;
+      boost::fusion::at_key<VMGDataFused::mixed_longrange>(instance) = 0.;
+      boost::fusion::at_key<VMGDataFused::mixed_shortrange>(instance) = 0.;
       boost::fusion::at_key<VMGDataFused::nuclei_longrange>(instance) = 0.;
       boost::fusion::at_key<VMGDataFused::nuclei_shortrange>(instance) = 0.;
@@ -136,36 +138,78 @@
     }
     for (size_t level = 2; level <= MaxLevel; ++level) {
-      // we have to fill in the remainder values in the LongRangeMap by hand
-      // weight times correct charge density of the same level
+      // We have calculated three different contributions: e-e, e-n+n-n, and n-n.
+      // And we want to have e-e+e-n, n-n+n-e (where e-n = n-e).
+      // For each of these three contributions we have a full solution and summed
+      // up short range solutions.
+
+      // first we obtain the full e-e energy as potential times charge on the
+      // respective level.
+      // \note that sampled_potential starts at level 2 because we do not calculate
+      // for level 1 as there saturated hydrogens are still present, leaving the
+      // result to be nonsense.
       const SamplingGrid &charge_weight =
           boost::fusion::at_key<MPQCDataFused::sampled_grid>(Result_Grid_fused[level-1]);
       SamplingGrid full_sample_solution = fullsolutionData[level-2].sampled_potential;
       full_sample_solution *= charge_weight;
+      double electron_solution_energy = full_sample_solution.integral();
+
+      // then we subtract the summed-up short-range e-e interaction energy from
+      // the full solution.
       const SamplingGrid &short_range_correction =
           boost::fusion::at_key<VMGDataFused::sampled_potential>(Result_LongRange_fused[level-1]);
+      double electron_short_range_energy = short_range_correction.integral();
       full_sample_solution -= short_range_correction;
-      double full_solution_energy = fullsolutionData[level-2].e_long;
-      const double short_range_energy =
-          boost::fusion::at_key<VMGDataFused::energy_long>(Result_LongRange_fused[level-1]);
-      full_solution_energy -= short_range_energy;
-
-      // multiply element-wise with charge distribution
+      electron_solution_energy -= electron_short_range_energy;
+      ASSERT( fabs(electron_solution_energy - full_sample_solution.integral()) < 1e-7,
+          "FragmentationLongRangeResults::operator() - integral and energy are not exchangeable.");
+
+      // then, we obtain the e-n+n-n full solution in the same way
+      double nuclei_solution_energy = fullsolutionData[level-2].nuclei_long;
+      double nuclei_short_range_energy =
+          boost::fusion::at_key<VMGDataFused::nuclei_long>(Result_LongRange_fused[level-1]);
+      nuclei_solution_energy -= nuclei_short_range_energy;
+
+      // and also the e-n full solution
+      double both_solution_energy = fullsolutionData[level-2].electron_long;
+      double both_short_range_energy =
+          boost::fusion::at_key<VMGDataFused::electron_long>(Result_LongRange_fused[level-1]);
+      both_solution_energy -= both_short_range_energy;
+
+      // energies from interpolation at nuclei position has factor of 1/2 already
+      electron_solution_energy *= .5;
+      electron_short_range_energy *= .5;
+
+      // At last, we subtract e-n from n-n+e-n for full solution and short-range
+      // correction.
+      nuclei_solution_energy -= both_solution_energy;
+      nuclei_short_range_energy -= both_short_range_energy;
+
       VMGDataLongRangeMap_t instance;
-      boost::fusion::at_key<VMGDataFused::electron_longrange>(instance) = .5*full_sample_solution.integral();
+      boost::fusion::at_key<VMGDataFused::electron_longrange>(instance) = electron_solution_energy;
 //      LOG(0, "Remaining long-range potential integral of level " << level << " is "
 //          << full_sample_solution.integral() << ".");
-      boost::fusion::at_key<VMGDataFused::electron_shortrange>(instance) = .5*short_range_correction.integral();
+      boost::fusion::at_key<VMGDataFused::electron_shortrange>(instance) = electron_short_range_energy;
 //      LOG(0, "Short-range correction potential integral of level " << level << " is "
 //          << short_range_correction.integral() << ".");
-      boost::fusion::at_key<VMGDataFused::nuclei_longrange>(instance) = full_solution_energy;
+      boost::fusion::at_key<VMGDataFused::mixed_longrange>(instance) = both_solution_energy;
 //      LOG(0, "Remaining long-range energy from potential integral of level " << level << " is "
 //          << full_solution_energy << ".");
-      boost::fusion::at_key<VMGDataFused::nuclei_shortrange>(instance) = short_range_energy;
+      boost::fusion::at_key<VMGDataFused::mixed_shortrange>(instance) = both_short_range_energy;
+//      LOG(0, "Short-range correction energy from potential integral of level " << level << " is "
+//          << short_range_energy << ".");
+      boost::fusion::at_key<VMGDataFused::nuclei_longrange>(instance) = nuclei_solution_energy;
+//      LOG(0, "Remaining long-range energy from potential integral of level " << level << " is "
+//          << full_solution_energy << ".");
+      boost::fusion::at_key<VMGDataFused::nuclei_shortrange>(instance) = nuclei_short_range_energy;
 //      LOG(0, "Short-range correction energy from potential integral of level " << level << " is "
 //          << short_range_energy << ".");
       boost::fusion::at_key<VMGDataFused::total_longrange>(instance) =
-          boost::fusion::at_key<VMGDataFused::electron_longrange>(instance) - full_solution_energy;
+          boost::fusion::at_key<VMGDataFused::electron_longrange>(instance)
+          + 2.*boost::fusion::at_key<VMGDataFused::mixed_longrange>(instance)
+          + boost::fusion::at_key<VMGDataFused::nuclei_longrange>(instance);
       boost::fusion::at_key<VMGDataFused::total_shortrange>(instance) =
-          boost::fusion::at_key<VMGDataFused::electron_shortrange>(instance) - short_range_energy;
+          boost::fusion::at_key<VMGDataFused::electron_shortrange>(instance)
+          + 2.*boost::fusion::at_key<VMGDataFused::mixed_shortrange>(instance)
+          + boost::fusion::at_key<VMGDataFused::nuclei_shortrange>(instance);
       Result_LongRangeIntegrated_fused.push_back(instance);
     }
Index: src/Fragmentation/Summation/Containers/VMGData.cpp
===================================================================
--- src/Fragmentation/Summation/Containers/VMGData.cpp	(revision 5281ff99623f1f5b7ed8ddd755f2d6daf60cc077)
+++ src/Fragmentation/Summation/Containers/VMGData.cpp	(revision e2925fd7f1c9a7ab18fe71b5a40598b65ddb05dc)
@@ -42,5 +42,6 @@
 VMGData::VMGData(const SamplingGridProperties &_props) :
   sampled_potential(_props),
-  e_long(0.)
+  nuclei_long(0.),
+  electron_long(0.)
 {}
 
@@ -50,5 +51,6 @@
 std::ostream & operator<<(std::ostream &ost, const VMGData &data)
 {
-  ost << "Long-Range energy: " << data.e_long << std::endl;
+  ost << "Nuclei long-Range energy: " << data.nuclei_long << std::endl;
+  ost << "Electron long-Range energy: " << data.electron_long << std::endl;
   return ost;
 }
Index: src/Fragmentation/Summation/Containers/VMGData.hpp
===================================================================
--- src/Fragmentation/Summation/Containers/VMGData.hpp	(revision 5281ff99623f1f5b7ed8ddd755f2d6daf60cc077)
+++ src/Fragmentation/Summation/Containers/VMGData.hpp	(revision e2925fd7f1c9a7ab18fe71b5a40598b65ddb05dc)
@@ -42,6 +42,8 @@
   //!> sampled potential on the grid as output
   SamplingGrid sampled_potential;
-  //!> long-range contribution to energy
-  double e_long;
+  //!> nuclei long-range contribution to energy
+  double nuclei_long;
+  //!> electron long-range contribution to energy
+  double electron_long;
 
 private:
@@ -52,7 +54,11 @@
   {
     ar & sampled_potential;
-    ar & e_long;
+    if (version > 0)
+      ar & nuclei_long;
+    ar & electron_long;
   }
 };
+
+BOOST_CLASS_VERSION(VMGData, 1)
 
 std::ostream & operator<<(std::ostream &ost, const VMGData &data);
Index: src/Fragmentation/Summation/Containers/VMGDataFused.hpp
===================================================================
--- src/Fragmentation/Summation/Containers/VMGDataFused.hpp	(revision 5281ff99623f1f5b7ed8ddd755f2d6daf60cc077)
+++ src/Fragmentation/Summation/Containers/VMGDataFused.hpp	(revision e2925fd7f1c9a7ab18fe71b5a40598b65ddb05dc)
@@ -27,5 +27,6 @@
   struct sampled_potential {};
   struct energy_potential {};
-  struct energy_long {};
+  struct nuclei_long {};
+  struct electron_long {};
 
   // keys for longrange
@@ -35,4 +36,6 @@
   struct nuclei_longrange {};
   struct nuclei_shortrange {};
+  struct mixed_longrange {};
+  struct mixed_shortrange {};
   struct total_longrange {};
   struct total_shortrange {};
Index: src/Fragmentation/Summation/Containers/VMGDataMap.hpp
===================================================================
--- src/Fragmentation/Summation/Containers/VMGDataMap.hpp	(revision 5281ff99623f1f5b7ed8ddd755f2d6daf60cc077)
+++ src/Fragmentation/Summation/Containers/VMGDataMap.hpp	(revision e2925fd7f1c9a7ab18fe71b5a40598b65ddb05dc)
@@ -32,5 +32,6 @@
     boost::fusion::pair<VMGDataFused::sampled_potential, SamplingGrid >,
     boost::fusion::pair<VMGDataFused::energy_potential, double >,
-    boost::fusion::pair<VMGDataFused::energy_long, double >
+    boost::fusion::pair<VMGDataFused::nuclei_long, double >,
+    boost::fusion::pair<VMGDataFused::electron_long, double >
 > VMGDataMap_t;
 
@@ -38,5 +39,6 @@
     VMGDataFused::sampled_potential,
     VMGDataFused::energy_potential,
-    VMGDataFused::energy_long
+    VMGDataFused::nuclei_long,
+    VMGDataFused::electron_long
 > VMGDataVector_t;
 
@@ -51,4 +53,6 @@
     boost::fusion::pair<VMGDataFused::electron_longrange, double >,
     boost::fusion::pair<VMGDataFused::electron_shortrange, double >,
+    boost::fusion::pair<VMGDataFused::mixed_longrange, double >,
+    boost::fusion::pair<VMGDataFused::mixed_shortrange, double >,
     boost::fusion::pair<VMGDataFused::nuclei_longrange, double >,
     boost::fusion::pair<VMGDataFused::nuclei_shortrange, double >,
@@ -60,4 +64,6 @@
     VMGDataFused::electron_longrange,
     VMGDataFused::electron_shortrange,
+    VMGDataFused::mixed_longrange,
+    VMGDataFused::mixed_shortrange,
     VMGDataFused::nuclei_longrange,
     VMGDataFused::nuclei_shortrange,
Index: src/Fragmentation/Summation/Containers/VMGData_printKeyNames.hpp
===================================================================
--- src/Fragmentation/Summation/Containers/VMGData_printKeyNames.hpp	(revision 5281ff99623f1f5b7ed8ddd755f2d6daf60cc077)
+++ src/Fragmentation/Summation/Containers/VMGData_printKeyNames.hpp	(revision e2925fd7f1c9a7ab18fe71b5a40598b65ddb05dc)
@@ -36,7 +36,10 @@
 	(sampled_potential) \
   (energy_potential) \
-  (energy_long) \
+  (nuclei_long) \
+  (electron_long) \
   (electron_longrange) \
   (electron_shortrange) \
+  (mixed_longrange) \
+  (mixed_shortrange) \
   (nuclei_longrange) \
   (nuclei_shortrange) \
Index: src/Fragmentation/Summation/Converter/DataConverter.hpp
===================================================================
--- src/Fragmentation/Summation/Converter/DataConverter.hpp	(revision 5281ff99623f1f5b7ed8ddd755f2d6daf60cc077)
+++ src/Fragmentation/Summation/Converter/DataConverter.hpp	(revision e2925fd7f1c9a7ab18fe71b5a40598b65ddb05dc)
@@ -90,5 +90,6 @@
     boost::fusion::at_key<VMGDataFused::sampled_potential>(instance) = extractedData.sampled_potential;
     boost::fusion::at_key<VMGDataFused::energy_potential>(instance) = extractedData.sampled_potential.integral();
-    boost::fusion::at_key<VMGDataFused::energy_long>(instance) = extractedData.e_long;
+    boost::fusion::at_key<VMGDataFused::nuclei_long>(instance) = extractedData.nuclei_long;
+    boost::fusion::at_key<VMGDataFused::electron_long>(instance) = extractedData.electron_long;
     VMGData_fused.insert( std::make_pair(dataiter->first, instance) );
   }
