Index: src/Jobs/InterfaceVMGJob.cpp
===================================================================
--- src/Jobs/InterfaceVMGJob.cpp	(revision 442ceedac1f17008644f48fc45034d4d28604b23)
+++ src/Jobs/InterfaceVMGJob.cpp	(revision cd77fcf679a54de7deec8445ce3cd9cccd538da9)
@@ -32,4 +32,8 @@
 #endif
 
+#ifdef HAVE_MPI
+#include "mpi.h"
+#endif
+
 #include "CodePatterns/MemDebug.hpp"
 
@@ -40,4 +44,10 @@
 
 #include "base/vector.hpp"
+#include "base/math.hpp"
+#ifdef HAVE_MPI
+#include "comm/comm_mpi.hpp"
+#else
+#include "comm/comm_serial.hpp"
+#endif
 #include "grid/grid.hpp"
 #include "grid/multigrid.hpp"
@@ -48,13 +58,84 @@
 using VMGInterfaces::InterfaceVMGJob;
 
+InterfaceVMGJob::InterfaceVMGJob(const std::vector< double > &_sampled_input,
+    std::vector< double > &_sampled_output,
+    const std::vector< std::vector<double> > &_particle_positions,
+    const std::vector< double > &_particle_charges,
+    VMG::Boundary boundary,
+    int levelMin,
+    int levelMax,
+    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),
+  spl(near_field_cells, Extent(MaxLevel()).MeshWidth().Max()),
+  sampled_input(_sampled_input),
+  sampled_output(_sampled_output),
+  level(levelMax)
+{
+  std::vector< std::vector<double> >::const_iterator positer = _particle_positions.begin();
+  std::vector<double>::const_iterator chargeiter = _particle_charges.begin();
+  double pos[3];
+  for (; positer != _particle_positions.end(); ++positer, ++chargeiter) {
+    ASSERT( (*positer).size() == 3,
+        "InterfaceVMGJob::InterfaceVMGJob() - particle "
+        +toString(distance(_particle_positions.begin(), positer))+" has not exactly 3 coordinates.");
+    for (size_t i=0;i<3;++i)
+      pos[i] = (*positer)[i];
+    particles.push_back(Particle::Particle(pos, *chargeiter));
+  }
+}
+
+InterfaceVMGJob::~InterfaceVMGJob()
+{}
+
 void InterfaceVMGJob::ImportRightHandSide(Multigrid& multigrid)
 {
   Index i;
   Vector pos;
-
-//  VMG::TempGrid *temp_grid = new VMG::TempGrid(129, 0, 0., 1.);
+  //  VMG::TempGrid *temp_grid = new VMG::TempGrid(129, 0, 0., 1.);
 
   Grid& grid = multigrid(multigrid.MaxLevel());
   grid.ClearBoundary();
+
+  /// 1. assign nuclei as smeared-out charges to the grid
+
+  /*
+   * Charge assignment on the grid
+   */
+#ifdef HAVE_MPI
+  CommMPI& comm = *dynamic_cast<CommMPI*>(VMG::MG::GetComm());
+#else
+  CommSerial& comm = *dynamic_cast<CommSerial*>(VMG::MG::GetComm());
+#endif
+  Grid& particle_grid = comm.GetParticleGrid();
+
+  particle_grid.Clear();
+
+  assert(particle_grid.Global().LocalSize().IsComponentwiseGreater(
+      VMG::MG::GetFactory().GetObjectStorageVal<int>("PARTICLE_NEAR_FIELD_CELLS")));
+
+  for (std::list<Particle::Particle>::iterator iter = particles.begin();
+      iter != particles.end(); ++iter)
+    spl.SetSpline(particle_grid, *iter);
+
+  // Communicate charges over halo
+  comm.CommFromGhosts(particle_grid);
+
+  // Assign charge values to the right hand side
+  for (int i=0; i<grid.Local().Size().X(); ++i)
+    for (int j=0; j<grid.Local().Size().Y(); ++j)
+      for (int k=0; k<grid.Local().Size().Z(); ++k)
+  grid(grid.Local().Begin().X() + i,
+       grid.Local().Begin().Y() + j,
+       grid.Local().Begin().Z() + k) = 4.0 * VMG::Math::pi *
+    particle_grid.GetVal(particle_grid.Local().Begin().X() + i,
+             particle_grid.Local().Begin().Y() + j,
+             particle_grid.Local().Begin().Z() + k);
+
+  /// 2. add sampled electron density to the grid
 
   const Index begin_local = grid.Global().LocalBegin() - grid.Local().HaloSize1();
@@ -71,15 +152,11 @@
       << "Z in [" << grid.Local().Begin().Z() << "," << grid.Local().End().Z() << "].");
 
+  const double element_volume =
+      grid.Extent().MeshWidth().X() * grid.Extent().MeshWidth().Y() * grid.Extent().MeshWidth().Z();
   std::vector<double>::const_iterator sample_iter = sampled_input.begin();
   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()) {
-        pos = grid.Extent().MeshWidth() * static_cast<Vector>(begin_local + i);
-//        R.x() = pos.X();
-//        R.y() = pos.Y();
-//        R.z() = pos.Z();
-        grid(i) = *sample_iter; //temp_grid(i);
-        ++sample_iter;
-      }
+      for (i.Z()=grid.Local().Begin().Z(); i.Z()<grid.Local().End().Z(); ++i.Z())
+        grid(i) -= (*sample_iter++) * element_volume; //temp_grid(i);
   assert( sample_iter == sampled_input.end() );
 
