Index: src/Jobs/Grid/SamplingGrid.hpp
===================================================================
--- src/Jobs/Grid/SamplingGrid.hpp	(revision 64bafe0d0cd15e6cd2e50cb7a27b730cfec01778)
+++ src/Jobs/Grid/SamplingGrid.hpp	(revision 8f3cdda52c1833250b3652985b1c78fc60b21a05)
@@ -177,5 +177,4 @@
   const double getWindowVolume() const;
 
-private:
   /** Sets the size of the window.
    *
@@ -187,4 +186,5 @@
   void setWindow(const double _begin_window[3], const double _end_window[3]);
 
+private:
   /** Sets the size of the domain.
    *
Index: src/Jobs/InterfaceVMGJob.cpp
===================================================================
--- src/Jobs/InterfaceVMGJob.cpp	(revision 64bafe0d0cd15e6cd2e50cb7a27b730cfec01778)
+++ src/Jobs/InterfaceVMGJob.cpp	(revision 8f3cdda52c1833250b3652985b1c78fc60b21a05)
@@ -52,12 +52,14 @@
 #include <cmath>
 #include <iostream>
+#include <limits>
 
 #include "CodePatterns/Log.hpp"
 
+#include "Jobs/Grid/SamplingGrid.hpp"
 
 using namespace VMG;
 using VMGInterfaces::InterfaceVMGJob;
 
-InterfaceVMGJob::InterfaceVMGJob(const std::vector< double > &_sampled_input,
+InterfaceVMGJob::InterfaceVMGJob(const SamplingGrid &_sampled_input,
     VMGData &_returndata,
     const std::vector< std::vector<double> > &_particle_positions,
@@ -66,11 +68,11 @@
     int levelMin,
     int levelMax,
-    const VMG::Vector &box_begin,
-    vmg_float box_end,
+    const VMG::Vector &_box_begin,
+    vmg_float _box_end,
     const int& near_field_cells,
     int coarseningSteps,
     double alpha) :
   VMG::Interface(boundary, levelMin, levelMax,
-      box_begin, box_end, coarseningSteps, alpha),
+      _box_begin, _box_end, coarseningSteps, alpha),
   spl(near_field_cells, Extent(MaxLevel()).MeshWidth().Max()),
   sampled_input(_sampled_input),
@@ -78,4 +80,8 @@
   level(levelMax)
 {
+  for (size_t i=0;i<3;++i) {
+    box_begin[i] = _box_begin[i];
+    box_end[i] = _box_end;
+  }
   std::vector< std::vector<double> >::const_iterator positer = _particle_positions.begin();
   std::vector<double>::const_iterator chargeiter = _particle_charges.begin();
@@ -94,4 +100,160 @@
 {}
 
+void InterfaceVMGJob::addGridOntoWindow(
+    VMG::Grid &grid,
+    SamplingGrid &window,
+    const double prefactor) const
+{
+#ifndef NDEBUG
+  for(size_t index=0;index<3;++index) {
+    ASSERT( box_begin[index] >= window.begin_window[index],
+        "InterfaceVMGJob::addGridOntoWindow() - given grid starts earlier than window in component "
+        +toString(index)+".");
+    ASSERT( box_end[index] <= window.end_window[index],
+        "InterfaceVMGJob::addGridOntoWindow() - given grid ends later than window in component "
+        +toString(index)+".");
+  }
+#endif
+  // the only issue are indices
+  const size_t gridpoints_axis = window.getGridPointsPerAxis();
+  size_t pre_offset[3];
+  size_t post_offset[3];
+  size_t length[3];
+  size_t total[3];
+  const double round_offset =
+      (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
+          0.5 : 0.; // need offset to get to round_toward_nearest behavior
+  for(size_t index=0;index<3;++index) {
+    const double delta = (double)gridpoints_axis/(window.end[index] - window.begin[index]);
+    // delta is conversion factor from box length to discrete length, i.e. number of points
+    pre_offset[index] = delta*(box_begin[index] - window.begin_window[index])+round_offset;
+    length[index] = delta*(box_end[index] - box_begin[index])+round_offset;
+    post_offset[index] = delta*(window.end_window[index] - box_end[index])+round_offset;
+    total[index] = delta*(window.end_window[index] - window.begin_window[index])+round_offset;
+    // total is used as safe-guard against loss due to discrete conversion
+    ASSERT( pre_offset[index]+post_offset[index]+length[index] == total[index],
+        "InterfaceVMGJob::addGridOntoWindow() - pre, post, and length don't sum up to total for "
+        +toString(index)+"th component.");
+  }
+#ifndef NDEBUG
+  const size_t calculated_size = length[0]*length[1]*length[2];
+//  const size_t given_size = std::distance(grid_begin, grid_end);
+//  ASSERT( calculated_size == given_size,
+//      "InterfaceVMGJob::addGridOntoWindow() - not enough sampled values given: "
+//      +toString(calculated_size)+" != "+toString(given_size)+".");
+  ASSERT( calculated_size <=  window.sampled_grid.size(),
+      "InterfaceVMGJob::addGridOntoWindow() - not enough sampled values available: "
+      +toString(calculated_size)+" <= "+toString(window.sampled_grid.size())+".");
+  const size_t total_size = total[0]*total[1]*total[2];
+  ASSERT( total_size == window.sampled_grid.size(),
+      "InterfaceVMGJob::addGridOntoWindow() - total size is not equal to number of present points: "
+      +toString(total_size)+" != "+toString(window.sampled_grid.size())+".");
+#endif
+  size_t N[3];
+  SamplingGrid::sampledvalues_t::iterator griditer = window.sampled_grid.begin();
+  std::advance(griditer, pre_offset[0]*total[1]*total[2]);
+  Grid::iterator copyiter = grid.Iterators().Local().Begin();
+  for(N[0]=0; N[0] < length[0]; ++N[0]) {
+    std::advance(griditer, pre_offset[1]*total[2]);
+    for(N[1]=0; N[1] < length[1]; ++N[1]) {
+      std::advance(griditer, pre_offset[2]);
+      for(N[2]=0; N[2] < length[2]; ++N[2]) {
+        ASSERT( griditer != window.sampled_grid.end(),
+            "InterfaceVMGJob::addGridOntoWindow() - griditer is already at end of window.");
+        ASSERT( copyiter != grid.Iterators().Local().End(),
+            "InterfaceVMGJob::addGridOntoWindow() - griditer is already at end of window.");
+        *griditer++ += prefactor*grid(*copyiter++);
+      }
+      std::advance(griditer, post_offset[2]);
+    }
+    std::advance(griditer, post_offset[1]*total[2]);
+  }
+#ifndef NDEBUG
+  std::advance(griditer, post_offset[0]*total[1]*total[2]);
+  ASSERT( griditer == window.sampled_grid.end(),
+      "InterfaceVMGJob::addGridOntoWindow() - griditer is not at end of window.");
+  ASSERT( copyiter == grid.Iterators().Local().End(),
+      "InterfaceVMGJob::addGridOntoWindow() - copyiter is not at end of window.");
+#endif
+}
+
+void InterfaceVMGJob::addWindowOntoGrid(
+    VMG::Grid& window,
+    const SamplingGrid &grid,
+    const double prefactor) const
+{
+#ifndef NDEBUG
+  for(size_t index=0;index<3;++index) {
+    ASSERT( grid.begin_window[index] >= box_begin[index],
+        "InterfaceVMGJob::addWindowOntoGrid() - given window starts earlier than grid in component "
+        +toString(index)+".");
+    ASSERT( grid.end_window[index] <= box_end[index],
+        "InterfaceVMGJob::addWindowOntoGrid() - given window ends later than grid in component "
+        +toString(index)+".");
+  }
+#endif
+  // the only issue are indices
+  const size_t gridpoints_axis = grid.getGridPointsPerAxis();
+  size_t pre_offset[3];
+  size_t post_offset[3];
+  size_t length[3];
+  size_t total[3];
+  const double round_offset =
+      (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
+          0.5 : 0.; // need offset to get to round_toward_nearest behavior
+  for(size_t index=0;index<3;++index) {
+    const double delta = (double)gridpoints_axis/(grid.end[index] - grid.begin[index]);
+    // delta is conversion factor from box length to discrete length, i.e. number of points
+    pre_offset[index] = delta*(grid.begin_window[index] - box_begin[index])+round_offset;
+    length[index] = delta*(grid.end_window[index] - grid.begin_window[index])+round_offset;
+    post_offset[index] = delta*(box_end[index] - grid.end_window[index])+round_offset;
+    total[index] = delta*(box_end[index] - box_begin[index])+round_offset;
+    // total is used as safe-guard against loss due to discrete conversion
+    ASSERT( pre_offset[index]+post_offset[index]+length[index] == total[index],
+        "InterfaceVMGJob::addWindowOntoGrid() - pre, post, and length don't sum up to total for "
+        +toString(index)+"th component.");
+  }
+#ifndef NDEBUG
+  const size_t calculated_size = length[0]*length[1]*length[2];
+  ASSERT( calculated_size == grid.sampled_grid.size(),
+      "InterfaceVMGJob::addWindowOntoGrid() - not enough sampled values given: "
+      +toString(calculated_size)+" != "+toString(grid.sampled_grid.size())+".");
+//  ASSERT( calculated_size <=  given_size,
+//      "InterfaceVMGJob::addWindowOntoGrid() - not enough sampled values available: "
+//      +toString(calculated_size)+" <= "+toString(given_size)+".");
+//  const size_t total_size = total[0]*total[1]*total[2];
+//  ASSERT( total_size == given_size,
+//      "InterfaceVMGJob::addWindowOntoGrid() - total size is not equal to number of present points: "
+//      +toString(total_size)+" != "+toString(given_size)+".");
+#endif
+  size_t N[3];
+  Grid::iterator griditer = window.Iterators().Local().Begin();
+  griditer += pre_offset[0]*total[1]*total[2];
+  SamplingGrid::sampledvalues_t::const_iterator copyiter = grid.sampled_grid.begin();
+  for(N[0]=0; N[0] < length[0]; ++N[0]) {
+    griditer += pre_offset[1]*total[2];
+    for(N[1]=0; N[1] < length[1]; ++N[1]) {
+      griditer += pre_offset[2];
+      for(N[2]=0; N[2] < length[2]; ++N[2]) {
+        ASSERT( griditer != window.Iterators().Local().End(),
+            "InterfaceVMGJob::addWindowOntoGrid() - griditer is already at end of window.");
+        ASSERT( copyiter != grid.sampled_grid.end(),
+            "InterfaceVMGJob::addWindowOntoGrid() - griditer is already at end of window.");
+        window(*griditer++) += prefactor*(*copyiter++);
+      }
+      griditer += post_offset[2];
+    }
+    griditer += post_offset[1]*total[2];
+  }
+#ifndef NDEBUG
+  griditer += post_offset[0]*total[1]*total[2];
+  ASSERT( griditer == window.Iterators().Local().End(),
+      "InterfaceVMGJob::addWindowOntoGrid() - griditer is not at end of window.");
+  ASSERT( copyiter == grid.sampled_grid.end(),
+      "InterfaceVMGJob::addWindowOntoGrid() - copyiter is not at end of window.");
+#endif
+//  LOG(2, "DEBUG: Grid after adding other is " << grid << ".");
+}
+
 void InterfaceVMGJob::ImportRightHandSide(Multigrid& multigrid)
 {
@@ -142,16 +304,14 @@
 
   // print nuclei grid to vtk
-//  comm.PrintGrid(particle_grid, "Sampled Nuclei Density");
+  comm.PrintGrid(particle_grid, "Sampled Nuclei Density");
 
   // add sampled electron charge density onto grid
-  std::vector<double>::const_iterator sample_iter = sampled_input.begin();
-  for (Grid::iterator iter = grid.Iterators().Local().Begin();
-      iter != grid.Iterators().Local().End();
-      ++iter)
-    grid(*iter) = -1. * (*sample_iter++);
-  assert( sample_iter == sampled_input.end() );
+  addWindowOntoGrid(
+      grid,
+      sampled_input,
+      -1.);
 
   // print electron grid to vtk
-//  comm.PrintGrid(grid, "Sampled Electron Density");
+  comm.PrintGrid(grid, "Sampled Electron Density");
 
   // add particle_grid onto grid
@@ -180,5 +340,5 @@
 
   // print total grid to vtk
-//  comm.PrintGrid(grid, "Total Charge Density");
+  comm.PrintGrid(grid, "Total Charge Density");
 
 //  delete temp_grid;
@@ -200,10 +360,13 @@
 
   // obtain sampled potential from grid
-  returndata.sampled_potential.sampled_grid.clear();
-  for (i.X()=grid.Local().Begin().X(); i.X()<grid.Local().End().X(); ++i.X())
-    for (i.Y()=grid.Local().Begin().Y(); i.Y()<grid.Local().End().Y(); ++i.Y())
-      for (i.Z()=grid.Local().Begin().Z(); i.Z()<grid.Local().End().Z(); ++i.Z()) {
-        returndata.sampled_potential.sampled_grid.push_back(grid(i));
-      }
+  returndata.sampled_potential.setWindow(
+      box_begin,
+      box_end
+      );
+  addGridOntoWindow(
+      grid,
+      returndata.sampled_potential,
+      +1.
+      );
 
   // calculate integral over potential as long-range energy contribution
Index: src/Jobs/InterfaceVMGJob.hpp
===================================================================
--- src/Jobs/InterfaceVMGJob.hpp	(revision 64bafe0d0cd15e6cd2e50cb7a27b730cfec01778)
+++ src/Jobs/InterfaceVMGJob.hpp	(revision 8f3cdda52c1833250b3652985b1c78fc60b21a05)
@@ -22,6 +22,9 @@
 #include "Jobs/VMGData.hpp"
 
+class SamplingGrid;
+
 namespace VMG
 {
+class Grid;
 class MGGrid;
 class MGMultigrid;
@@ -34,5 +37,5 @@
 {
 public:
-  InterfaceVMGJob(const std::vector< double > &_sampled_input,
+  InterfaceVMGJob(const SamplingGrid &_sampled_input,
       VMGData &returndata,
       const std::vector< std::vector< double > > &_particle_positions,
@@ -41,6 +44,6 @@
       int levelMin,
       int levelMax,
-      const VMG::Vector &box_begin,
-      vmg_float box_end,
+      const VMG::Vector &_box_begin,
+      vmg_float _box_end,
       const int& near_field_cells,
       int coarseningSteps=0,
@@ -56,6 +59,30 @@
 
 private:
+  /** Helper function to copy VMG's Grid onto the non-zero window of a SamplingGrid
+   *
+   * The code is very similar to SamplingGrid::addOntoWindow().
+   *
+   * @param grid grid with non-zero window
+   * @param prefactor factor for each added value, e.g. -1. makes a subtraction
+   */
+  void addGridOntoWindow(
+      VMG::Grid& grid,
+      SamplingGrid &window,
+      const double prefactor) const;
+
+  /** Helper function to copy non-zero window in SamplingGrid onto VMG's Grid.
+   *
+   * The code is very similar to SamplingGrid::addOntoWindow().
+   *
+   * @param grid grid with non-zero window
+   * @param prefactor factor for each added value, e.g. -1. makes a subtraction
+   */
+  void addWindowOntoGrid(
+      VMG::Grid& window,
+      const SamplingGrid &grid,
+      const double prefactor) const;
+private:
   //!> sampled density on the grid as input
-  const std::vector< double > sampled_input;
+  const SamplingGrid sampled_input;
   //!> sampled potential on the grid as output
   VMGData &returndata;
@@ -64,4 +91,9 @@
   //!> nuclei charges
   std::list<VMG::Particle::Particle> particles;
+
+  //!> contains start of grid intervals
+  double box_begin[3];
+  //!> contains end of grid intervals
+  double box_end[3];
 };
 
Index: src/Jobs/VMGJob.cpp
===================================================================
--- src/Jobs/VMGJob.cpp	(revision 64bafe0d0cd15e6cd2e50cb7a27b730cfec01778)
+++ src/Jobs/VMGJob.cpp	(revision 8f3cdda52c1833250b3652985b1c78fc60b21a05)
@@ -207,5 +207,5 @@
   new DiscretizationPoissonFD(4);
   new VMGInterfaces::InterfaceVMGJob(
-      density_grid.sampled_grid,
+      density_grid,
       returndata,
       particle_positions,
