Index: src/Actions/FragmentationAction/FragmentationAutomationAction.cpp
===================================================================
--- src/Actions/FragmentationAction/FragmentationAutomationAction.cpp	(revision 964c6386220871ef19613b809fb383601ae9e39c)
+++ src/Actions/FragmentationAction/FragmentationAutomationAction.cpp	(revision 82034f10efe37b9dbb18b67c42337e3e0d789109)
@@ -152,4 +152,69 @@
   return result;
 }
+
+
+#ifdef HAVE_JOBMARKET
+/** This performs a direct sum calculation for the nuclei-nuclei interaction energy
+ * and forces for each fragment contained in \a shortrangedata.
+ *
+ * Results are (re)placed in \a _longrangedata.
+ *
+ * \warning This is only valid for fully open boundary conditions.
+ *
+ * \param _shortrangedata containing information on position and charges per fragment
+ * \param _longrangedata containing energy and forces on return
+ */
+static void calculateNucleiNucleiLongRangeContribution(
+    const std::map<JobId_t, MPQCData> &_shortrangedata,
+    std::map<JobId_t, VMGData>& _longrangedata
+    )
+{
+  ASSERT( _shortrangedata.size() == _longrangedata.size(),
+      "calculateNucleiNucleiLongRangeContribution() - shortrange and longrange data have inequal size.");
+  std::map<JobId_t, VMGData>::iterator longrangeiter = _longrangedata.begin();
+  std::map<JobId_t, MPQCData>::const_iterator iter = _shortrangedata.begin();
+  for (;iter != _shortrangedata.end(); ++iter, ++longrangeiter) {
+    const MPQCData &data = iter->second;
+    VMGData &longrange_data = longrangeiter->second;
+    // set long-range contributions to zero
+    longrange_data.nuclei_long = 0.;
+    longrange_data.forces.clear();
+    longrange_data.forces.resize(data.positions.size(), FragmentForces::force_t(3,0.));
+    longrange_data.hasForces = true;
+    // go through positions and evaluate sum naively
+    ASSERT( data.positions.size() == data.charges.size(),
+        "calculateNucleiNucleiLongRangeContribution() - positions and charges differ in size.");
+    std::vector< std::vector<double> >::const_iterator positer = data.positions.begin();
+    std::vector< double >::const_iterator chargeiter = data.charges.begin();
+    FragmentForces::iterator forceiter = longrange_data.forces.begin();
+    for (; (positer != data.positions.end()) && (chargeiter != data.charges.end());
+        ++positer, ++chargeiter, ++forceiter) {
+      ASSERT( positer->size() == NDIM,
+          "calculateNucleiNucleiLongRangeContribution() - position "
+          +toString(std::distance(data.positions.begin(), positer))+" has not 3 components.");
+      const Vector position((*positer)[0], (*positer)[1], (*positer)[2]);
+
+      std::vector< std::vector<double> >::const_iterator otherpositer = data.positions.begin();
+      std::vector< double >::const_iterator otherchargeiter = data.charges.begin();
+      for (; (otherpositer != data.positions.end()) && (otherchargeiter != data.charges.end());
+          ++otherpositer, ++otherchargeiter) {
+        if ((positer == otherpositer) || (chargeiter == otherchargeiter))
+          continue;
+        ASSERT( otherpositer->size() == NDIM,
+            "calculateNucleiNucleiLongRangeContribution() - other position "
+          +toString(std::distance(data.positions.begin(), otherpositer))+" has not 3 components.");
+        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;
+        (*forceiter)[0] = factor*distance[0];
+        (*forceiter)[1] = factor*distance[1];
+        (*forceiter)[2] = factor*distance[2];
+        longrange_data.nuclei_long += 0.5*(*chargeiter)*(*otherchargeiter)*invsqrdist;
+      }
+    }
+  }
+}
+#endif
 
 ActionState::ptr FragmentationFragmentationAutomationAction::performCall() {
@@ -387,4 +452,12 @@
       Exitflag += vmgcontroller.getExitflag();
 
+      // now replace nuclei-nuclei contribution with values from direct summation (is exact for
+      // open boundary conditions
+      //
+      // NOTE: We need to calculate both_sampled_potential for fitting partial charges,
+      // nowhere else is this quantity needed.
+      //if (OpenBoundaryConditions)
+      //  calculateNucleiNucleiLongRangeContribution(shortrangedata, longrangedata_both);
+
       // go through either data and replace nuclei_long with contribution from both
       ASSERT( longrangedata.size() == longrangedata_both.size(),
