Index: src/Actions/FragmentationAction/AnalyseFragmentationResultsAction.cpp
===================================================================
--- src/Actions/FragmentationAction/AnalyseFragmentationResultsAction.cpp	(revision b6b21a5af3ef8085052c11ae4c00c67da135353a)
+++ src/Actions/FragmentationAction/AnalyseFragmentationResultsAction.cpp	(revision e2925fd7f1c9a7ab18fe71b5a40598b65ddb05dc)
@@ -388,4 +388,7 @@
     }
   }
+#else
+  if (DoLongrange)
+    ELOG(2, "File contains long-range information but long-range analysis capability not compiled in.");
 #endif
 
Index: src/Actions/FragmentationAction/FragmentationAutomationAction.cpp
===================================================================
--- src/Actions/FragmentationAction/FragmentationAutomationAction.cpp	(revision b6b21a5af3ef8085052c11ae4c00c67da135353a)
+++ src/Actions/FragmentationAction/FragmentationAutomationAction.cpp	(revision e2925fd7f1c9a7ab18fe71b5a40598b65ddb05dc)
@@ -153,7 +153,7 @@
     vmgcontroller.setPort(params.port.get());
     const size_t NoJobs = shortrangedata.size()+full_sample.size();
-    vmgcontroller.requestIds(NoJobs);
-
-    // Phase Five: create VMGJobs
+    vmgcontroller.requestIds(2*NoJobs);
+
+    // Phase Five a: create VMGJobs for electronic charge distribution
     const size_t near_field_cells = params.near_field_cells.get();
     const size_t interpolation_degree = params.interpolation_degree.get();
@@ -163,11 +163,11 @@
         near_field_cells,
         interpolation_degree,
+        VMGFragmentController::DontSampleParticles,
+        VMGFragmentController::DoTreatGrid,
         params.DoValenceOnly.get() ? MPQCData::DoSampleValenceOnly : MPQCData::DontSampleValenceOnly,
         params.DoPrintDebug.get()))
-      // can't use the summed up fragment due to saturation hydrogens
-//      summedChargeDensity.getFragment(),
       return Action::failure;
 
-    // Phase Six: calculate result
+    // Phase Six a: calculate result
     vmgcontroller.waitforResults(NoJobs);
     vmgcontroller.getResults(longrangedata);
@@ -175,6 +175,40 @@
         "FragmentationFragmentationAutomationAction::performCall() - number of MPQCresults+"
         +toString(full_sample.size())+"="+toString(NoJobs)
-        +" and VMGresults "+toString(longrangedata.size())+" don't match.");
+        +" and first VMGresults "+toString(longrangedata.size())+" don't match.");
     Exitflag += vmgcontroller.getExitflag();
+
+    {
+      std::map<JobId_t, VMGData> longrangedata_both;
+      // Phase Five b: create VMGJobs for nuclei charge distributions
+      const size_t near_field_cells = params.near_field_cells.get();
+      const size_t interpolation_degree = params.interpolation_degree.get();
+      if (!vmgcontroller.createLongRangeJobs(
+          shortrangedata,
+          full_sample,
+          near_field_cells,
+          interpolation_degree,
+          VMGFragmentController::DoSampleParticles,
+          VMGFragmentController::DoTreatGrid,
+          params.DoValenceOnly.get() ? MPQCData::DoSampleValenceOnly : MPQCData::DontSampleValenceOnly,
+          params.DoPrintDebug.get()))
+        return Action::failure;
+
+      // Phase Six b: calculate result
+      vmgcontroller.waitforResults(NoJobs);
+      vmgcontroller.getResults(longrangedata_both);
+      ASSERT( NoJobs == longrangedata_both.size(),
+          "FragmentationFragmentationAutomationAction::performCall() - number of MPQCresults+"
+          +toString(full_sample.size())+"="+toString(NoJobs)
+          +" and second VMGresults "+toString(longrangedata_both.size())+" don't match.");
+      Exitflag += vmgcontroller.getExitflag();
+
+      // go through either data and replace nuclei_long with contribution from both
+      ASSERT( longrangedata.size() == longrangedata_both.size(),
+          "FragmentationFragmentationAutomationAction::performCall() - longrange results have different sizes.");
+      std::map<JobId_t, VMGData>::iterator destiter = longrangedata.begin();
+      std::map<JobId_t, VMGData>::iterator srciter = longrangedata_both.begin();
+      for (;destiter != longrangedata.end(); ++srciter, ++destiter)
+        destiter->second.nuclei_long = srciter->second.nuclei_long;
+    }
   }
 
