Index: src/Actions/FragmentationAction/AnalyseFragmentationResultsAction.cpp
===================================================================
--- src/Actions/FragmentationAction/AnalyseFragmentationResultsAction.cpp	(revision 3fb9ab59e65e8b93b5eeb0d1bf0ccdedd5cc8e06)
+++ src/Actions/FragmentationAction/AnalyseFragmentationResultsAction.cpp	(revision a8961da7eefe1d31d2558ca8bd178aaa32621cec)
@@ -392,15 +392,17 @@
       value.potential_distribution =
           boost::fusion::at_key<VMGDataFused::both_sampled_potential>(potentialiter->second.second); // contributions
-//      // and re-average to zero (integral is times volume_element which we don't need here)
-//      const double sum =
-//          std::accumulate(
-//              value.potential_distribution.sampled_grid.begin(),
-//              value.potential_distribution.sampled_grid.end(),
-//              0.);
-//      const double offset = sum/(double)value.potential_distribution.sampled_grid.size();
-//      for (SamplingGrid::sampledvalues_t::iterator iter = value.potential_distribution.sampled_grid.begin();
-//          iter != value.potential_distribution.sampled_grid.end();
-//          ++iter)
-//        *iter -= offset;
+      value.potential_distribution +=
+          boost::fusion::at_key<VMGDataFused::sampled_potential>(potentialiter->second.second); // contributions
+      // and re-average to zero (integral is times volume_element which we don't need here)
+      const double sum =
+          std::accumulate(
+              value.potential_distribution.sampled_grid.begin(),
+              value.potential_distribution.sampled_grid.end(),
+              0.);
+      const double offset = sum/(double)value.potential_distribution.sampled_grid.size();
+      for (SamplingGrid::sampledvalues_t::iterator iter = value.potential_distribution.sampled_grid.begin();
+          iter != value.potential_distribution.sampled_grid.end();
+          ++iter)
+        *iter -= offset;
 #else
       ELOG(2, "Long-range information in homology desired but long-range analysis capability not compiled in.");
Index: src/Actions/FragmentationAction/FragmentationAutomationAction.cpp
===================================================================
--- src/Actions/FragmentationAction/FragmentationAutomationAction.cpp	(revision 3fb9ab59e65e8b93b5eeb0d1bf0ccdedd5cc8e06)
+++ src/Actions/FragmentationAction/FragmentationAutomationAction.cpp	(revision a8961da7eefe1d31d2558ca8bd178aaa32621cec)
@@ -170,6 +170,7 @@
     )
 {
-  ASSERT( _shortrangedata.size() == _longrangedata.size(),
-      "calculateNucleiNucleiLongRangeContribution() - shortrange and longrange data have inequal size.");
+  // we assume here that the fragmentjobs in short- and longrange case are in the
+  // same order and that the full jobs per level are appended to the end of
+  // longrangedata
   std::map<JobId_t, VMGData>::iterator longrangeiter = _longrangedata.begin();
   std::map<JobId_t, MPQCData>::const_iterator iter = _shortrangedata.begin();
@@ -178,7 +179,8 @@
     VMGData &longrange_data = longrangeiter->second;
     // set long-range contributions to zero
+    const double old_energy = longrange_data.nuclei_long;
     longrange_data.nuclei_long = 0.;
     longrange_data.forces.clear();
-    longrange_data.forces.resize(data.positions.size(), FragmentForces::force_t(3,0.));
+    longrange_data.forces.resize(data.positions.size(), FragmentForces::force_t(NDIM,0.));
     longrange_data.hasForces = true;
     // go through positions and evaluate sum naively
@@ -206,12 +208,14 @@
         const Vector otherposition((*otherpositer)[0], (*otherpositer)[1], (*otherpositer)[2]);
         const Vector distance = position - otherposition;
-        const double invsqrdist = 1./distance.Norm();
-        const double factor = (*otherchargeiter)*invsqrdist*invsqrdist*invsqrdist;
+        const double invdist = 1./distance.Norm();
+        const double factor = (*otherchargeiter)*invdist*invdist*invdist;
         (*forceiter)[0] = factor*distance[0];
         (*forceiter)[1] = factor*distance[1];
         (*forceiter)[2] = factor*distance[2];
-        longrange_data.nuclei_long += 0.5*(*chargeiter)*(*otherchargeiter)*invsqrdist;
-      }
-    }
+        longrange_data.nuclei_long += 0.5*(*chargeiter)*(*otherchargeiter)*invdist;
+      }
+    }
+    LOG(2, "DEBUG: fragment #" << longrangeiter->first << ": old is " << old_energy
+        << ", new is " << longrange_data.nuclei_long);
   }
 }
@@ -423,5 +427,5 @@
           interpolation_degree,
           VMGFragmentController::DoSampleParticles,
-          VMGFragmentController::DoTreatGrid,
+          VMGFragmentController::DontTreatGrid,
           params.DoValenceOnly.get() ? MPQCData::DoSampleValenceOnly : MPQCData::DontSampleValenceOnly,
           params.DoPrintDebug.get(),
@@ -468,6 +472,6 @@
         destiter->second.both_sampled_potential = srciter->second.sampled_potential;
         destiter->second.nuclei_long = srciter->second.nuclei_long;
-        destiter->second.forces = srciter->second.forces;
-        destiter->second.hasForces = srciter->second.hasForces;
+        destiter->second.forces += srciter->second.forces;
+        destiter->second.hasForces &= srciter->second.hasForces;
       }
     }
Index: src/Fragmentation/Summation/Containers/FragmentationLongRangeResults.cpp
===================================================================
--- src/Fragmentation/Summation/Containers/FragmentationLongRangeResults.cpp	(revision 3fb9ab59e65e8b93b5eeb0d1bf0ccdedd5cc8e06)
+++ src/Fragmentation/Summation/Containers/FragmentationLongRangeResults.cpp	(revision a8961da7eefe1d31d2558ca8bd178aaa32621cec)
@@ -121,5 +121,5 @@
         Result_perIndexSet_Grid);
 
-    // multiply each short-range potential with the respective charge
+    // multiply each short-range e-e potential with the respective charge
     std::map<JobId_t,MPQCData>::const_iterator mpqciter = fragmentData.begin();
     std::map<JobId_t,VMGData>::iterator vmgiter = longrangeData.begin();
@@ -214,5 +214,5 @@
 #endif
 
-      // then, we obtain the e-n+n-n full solution in the same way
+      // then, we obtain the n-n full solution directly
       double nuclei_solution_energy = fullsolutionData[level-1].nuclei_long;
       double nuclei_short_range_energy =
@@ -229,9 +229,4 @@
       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;
Index: src/Jobs/InterfaceVMGJob.cpp
===================================================================
--- src/Jobs/InterfaceVMGJob.cpp	(revision 3fb9ab59e65e8b93b5eeb0d1bf0ccdedd5cc8e06)
+++ src/Jobs/InterfaceVMGJob.cpp	(revision a8961da7eefe1d31d2558ca8bd178aaa32621cec)
@@ -116,6 +116,6 @@
 
   Grid& grid = multigrid(multigrid.MaxLevel());
-  grid.Clear();
-  //grid.ClearBoundary(); // we don't have a boundary under periodic boundary conditions
+//  grid.Clear();
+//  grid.ClearBoundary(); // we don't have a boundary under periodic boundary conditions
 
   // print debugging info on grid size
@@ -393,6 +393,6 @@
   comm.PrintOnce(Debug, "E_total*:       %e", e_long + e_short_peak + e_short_spline - e_self);
 
-  returndata.nuclei_long = e_long;
-  returndata.electron_long = e_long;
+  returndata.nuclei_long = e;   // nuclei needs self-energy subtracted, ...
+  returndata.electron_long = e_long;  // charge grid gets only long-range interaction
 
   // calculate residual
