Index: src/Jobs/InterfaceVMGJob.cpp
===================================================================
--- src/Jobs/InterfaceVMGJob.cpp	(revision 5a5196ce72216409382904043b36c0b0e8c3d422)
+++ src/Jobs/InterfaceVMGJob.cpp	(revision a82602d89e2fa235117faebd8927979ae430d56d)
@@ -41,4 +41,7 @@
 #include "grid/grid.hpp"
 #include "grid/multigrid.hpp"
+#include "units/particle/comm_mpi_particle.hpp"
+#include "units/particle/interpolation.hpp"
+#include "units/particle/linked_cell_list.hpp"
 #include "mg.hpp"
 
@@ -98,4 +101,5 @@
 
   Grid& grid = multigrid(multigrid.MaxLevel());
+  grid.Clear();
   //grid.ClearBoundary(); // we don't have a boundary under periodic boundary conditions
 
@@ -114,8 +118,11 @@
    * Charge assignment on the grid
    */
-  Comm& comm = *MG::GetComm();
+  Particle::CommMPI& comm = *dynamic_cast<Particle::CommMPI*>(MG::GetComm());
   Grid& particle_grid = comm.GetParticleGrid();
-
   particle_grid.Clear();
+
+  // distribute particles
+  particles.clear();
+  comm.CommParticles(grid, particles);
 
   assert(particle_grid.Global().LocalSize().IsComponentwiseGreater(
@@ -180,6 +187,9 @@
 void InterfaceVMGJob::ExportSolution(Grid& grid)
 {
+  /// sample the obtained potential to evaluate with the electron charge density
+
   // grid now contains the sough-for potential
-  Comm& comm = *MG::GetComm();
+  //Comm& comm = *MG::GetComm();
+  Particle::CommMPI& comm = *dynamic_cast<Particle::CommMPI*>(MG::GetComm());
 
   const Index begin_local = grid.Global().LocalBegin() - grid.Local().HaloSize1();
@@ -207,5 +217,131 @@
   comm.PrintStringOnce("Grid potential sum: %e", potential_sum);
 
-  //Particle::CommMPI& comm = *dynamic_cast<Particle::CommMPI*>(MG::GetComm());
-  returndata.e_long = potential_sum;
+  {
+    Grid::iterator grid_iter = grid.Iterators().Local().Begin();
+    comm.PrintStringOnce("Grid potential at (0,0,0): %e", grid.GetVal(*grid_iter));
+  }
+
+  //Particle::CommMPI& comm = *dynamic_cast<Particle::CommMPI*>(MG::GetComm());  returndata.e_long = potential_sum;
+
+  /// Calculate potential energy of nuclei
+
+  vmg_float e = 0.0;
+  vmg_float e_long = 0.0;
+  vmg_float e_self = 0.0;
+  vmg_float e_short_peak = 0.0;
+  vmg_float e_short_spline = 0.0;
+
+  Factory& factory = MG::GetFactory();
+
+  /*
+   * Get parameters and arrays
+   */
+  const vmg_int& near_field_cells = factory.GetObjectStorageVal<int>("PARTICLE_NEAR_FIELD_CELLS");
+  const vmg_int& interpolation_degree = factory.GetObjectStorageVal<int>("PARTICLE_INTERPOLATION_DEGREE");
+
+  Particle::Interpolation ip(interpolation_degree);
+
+  const vmg_float r_cut = near_field_cells * grid.Extent().MeshWidth().Max();
+
+  /*
+   * Copy potential values to a grid with sufficiently large halo size.
+   * This may be optimized in future.
+   * The parameters of this grid have been set in the import step.
+   */
+  Grid& particle_grid = comm.GetParticleGrid();
+
+  for (i.X()=0; i.X()<grid.Local().Size().X(); ++i.X())
+    for (i.Y()=0; i.Y()<grid.Local().Size().Y(); ++i.Y())
+      for (i.Z()=0; i.Z()<grid.Local().Size().Z(); ++i.Z())
+        particle_grid(i + particle_grid.Local().Begin()) = grid.GetVal(i + grid.Local().Begin());
+
+  comm.CommToGhosts(particle_grid);
+
+  /*
+   * Compute potentials
+   */
+  Particle::LinkedCellList lc(particles, near_field_cells, grid);
+  Particle::LinkedCellList::iterator p1, p2;
+  Grid::iterator iter;
+
+  comm.CommLCListToGhosts(lc);
+
+  for (int i=lc.Local().Begin().X(); i<lc.Local().End().X(); ++i)
+    for (int j=lc.Local().Begin().Y(); j<lc.Local().End().Y(); ++j)
+      for (int k=lc.Local().Begin().Z(); k<lc.Local().End().Z(); ++k) {
+
+  if (lc(i,j,k).size() > 0)
+    ip.ComputeCoefficients(particle_grid, Index(i,j,k) - lc.Local().Begin() + particle_grid.Local().Begin());
+
+  for (p1=lc(i,j,k).begin(); p1!=lc(i,j,k).end(); ++p1) {
+
+    // Interpolate long-range part of potential and electric field
+    ip.Evaluate(**p1);
+
+    // Subtract self-induced potential
+    (*p1)->Pot() -= (*p1)->Charge() * spl.GetAntiDerivativeAtZero();
+
+    e_long += 0.5 * (*p1)->Charge() * ip.EvaluatePotentialLR(**p1);
+    e_self += 0.5 * (*p1)->Charge() * (*p1)->Charge() * spl.GetAntiDerivativeAtZero();
+
+    for (int dx=-1*near_field_cells; dx<=near_field_cells; ++dx)
+      for (int dy=-1*near_field_cells; dy<=near_field_cells; ++dy)
+        for (int dz=-1*near_field_cells; dz<=near_field_cells; ++dz) {
+
+    for (p2=lc(i+dx,j+dy,k+dz).begin(); p2!=lc(i+dx,j+dy,k+dz).end(); ++p2)
+
+      if (*p1 != *p2) {
+
+        const Vector dir = (*p1)->Pos() - (*p2)->Pos();
+        const vmg_float length = dir.Length();
+
+        if (length < r_cut) {
+
+          (*p1)->Pot() += (*p2)->Charge() / length * (1.0 + spl.EvaluatePotential(length));
+          (*p1)->Field() += (*p2)->Charge() * dir * spl.EvaluateField(length);
+
+          e_short_peak += 0.5 * (*p1)->Charge() * (*p2)->Charge() / length;
+          e_short_spline += 0.5 * (*p1)->Charge() * (*p2)->Charge() / length * spl.EvaluatePotential(length);
+        }
+      }
+        }
+  }
+      }
+
+  /* Remove average force term */
+  Vector average_force = 0.0;
+  for (std::list<Particle::Particle>::const_iterator iter=particles.begin(); iter!=particles.end(); ++iter)
+    average_force += iter->Charge() * iter->Field();
+  const vmg_int& npl = MG::GetFactory().GetObjectStorageVal<vmg_int>("PARTICLE_NUM_LOCAL");
+  const vmg_int num_particles_global = comm.GlobalSum(npl);
+  average_force /= num_particles_global;
+  comm.GlobalSumArray(average_force.vec(), 3);
+  for (std::list<Particle::Particle>::iterator iter=particles.begin(); iter!=particles.end(); ++iter)
+    iter->Field() -= average_force / iter->Charge();
+
+  comm.CommParticlesBack(particles);
+
+  vmg_float* q = factory.GetObjectStorageArray<vmg_float>("PARTICLE_CHARGE_ARRAY");
+  const vmg_int& num_particles_local = factory.GetObjectStorageVal<vmg_int>("PARTICLE_NUM_LOCAL");
+  const vmg_float* p = factory.GetObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY");
+//  const vmg_float* f = factory.GetObjectStorageArray<vmg_float>("PARTICLE_FIELD_ARRAY");
+
+
+  e_long = comm.GlobalSumRoot(e_long);
+  e_short_peak = comm.GlobalSumRoot(e_short_peak);
+  e_short_spline = comm.GlobalSumRoot(e_short_spline);
+  e_self = comm.GlobalSumRoot(e_self);
+
+  for (int j=0; j<num_particles_local; ++j)
+    e += 0.5 * p[j] * q[j];
+  e = comm.GlobalSumRoot(e);
+
+  comm.PrintStringOnce("E_long:         %e", e_long);
+  comm.PrintStringOnce("E_short_peak:   %e", e_short_peak);
+  comm.PrintStringOnce("E_short_spline: %e", e_short_spline);
+  comm.PrintStringOnce("E_self:         %e", e_self);
+  comm.PrintStringOnce("E_total:        %e", e);
+  comm.PrintStringOnce("E_total*:       %e", e_long + e_short_peak + e_short_spline - e_self);
+
+  returndata.e_long = e;
 }
Index: src/Jobs/VMGJob.cpp
===================================================================
--- src/Jobs/VMGJob.cpp	(revision 5a5196ce72216409382904043b36c0b0e8c3d422)
+++ src/Jobs/VMGJob.cpp	(revision a82602d89e2fa235117faebd8927979ae430d56d)
@@ -63,5 +63,5 @@
 //#include "solver/solver_regular.hpp"
 #include "solver/solver_singular.hpp"
-//#include "units/particle/comm_mpi_particle.hpp"
+#include "units/particle/comm_mpi_particle.hpp"
 
 #include "CodePatterns/MemDebug.hpp"
@@ -71,4 +71,6 @@
 #include "LinearAlgebra/defs.hpp"
 #include "Jobs/InterfaceVMGJob.hpp"
+
+#include "CodePatterns/Assert.hpp"
 
 using namespace VMG;
@@ -85,14 +87,63 @@
   particle_charges(_particle_charges),
   near_field_cells(_near_field_cells),
-  returndata(static_cast<const SamplingGridProperties &>(_density_grid))
+  returndata(static_cast<const SamplingGridProperties &>(_density_grid)),
+  num_particles(0),
+  f(NULL),
+  x(NULL),
+  p(NULL),
+  q(NULL)
 {}
 
 VMGJob::VMGJob() :
   FragmentJob(JobId::IllegalJob),
-  near_field_cells(5)
+  near_field_cells(5),
+  num_particles(0),
+  f(NULL),
+  x(NULL),
+  p(NULL),
+  q(NULL)
 {}
 
 VMGJob::~VMGJob()
-{}
+{
+  delete[] f;
+  delete[] x;
+  delete[] p;
+  delete[] q;
+}
+
+void VMGJob::InitVMGArrays()
+{
+  num_particles = particle_charges.size();
+  f = new double[num_particles*3];
+  x = new double[num_particles*3];
+  p = new double[num_particles];
+  q = new double[num_particles];
+
+  {
+    size_t index=0;
+    for (std::vector< std::vector< double > >::const_iterator iter = particle_positions.begin();
+        iter != particle_positions.end(); ++iter) {
+      for (std::vector< double >::const_iterator positer = (*iter).begin();
+          positer != (*iter).end(); ++positer) {
+        f[index] = 0.;
+        x[index++] = *positer;
+      }
+    }
+    ASSERT( index == num_particles*3,
+        "VMGJob::VMGJob() - too many particles or components in particle_positions.");
+  }
+  {
+    size_t index = 0;
+    for (std::vector< double >::const_iterator iter = particle_charges.begin();
+        iter != particle_charges.end(); ++iter) {
+      p[index] = 0.;
+      q[index++] = *iter;
+    }
+    ASSERT( index == num_particles,
+        "VMGJob::VMGJob() - too many charges in particle_charges.");
+  }
+}
+
 
 FragmentResult::ptr VMGJob::Work()
@@ -137,5 +188,5 @@
    */
 #ifdef HAVE_MPI
-  new CommMPI(boundary, new DomainDecompositionMPI());
+  new Particle::CommMPI(boundary, new DomainDecompositionMPI());
 #else
   new CommSerial(boundary);
@@ -169,4 +220,13 @@
   new ObjectStorage<int>("PARTICLE_INTERPOLATION_DEGREE", 3);
 
+  // first initialize f,x,p,q,... from STL vectors
+  InitVMGArrays();
+  // then initialize as objects with VMG's storage
+  new ObjectStorage<vmg_float*>("PARTICLE_POS_ARRAY", x);
+  new ObjectStorage<vmg_float*>("PARTICLE_CHARGE_ARRAY", q);
+  new ObjectStorage<vmg_float*>("PARTICLE_POTENTIAL_ARRAY", p);
+  new ObjectStorage<vmg_float*>("PARTICLE_FIELD_ARRAY", f);
+  new ObjectStorage<vmg_int>("PARTICLE_NUM_LOCAL", num_particles);
+
   /*
    * Post init
Index: src/Jobs/VMGJob.hpp
===================================================================
--- src/Jobs/VMGJob.hpp	(revision 5a5196ce72216409382904043b36c0b0e8c3d422)
+++ src/Jobs/VMGJob.hpp	(revision a82602d89e2fa235117faebd8927979ae430d56d)
@@ -53,4 +53,6 @@
   void InitVMG();
 
+  void InitVMGArrays();
+
 private:
   //!> sampled density required as input
@@ -66,4 +68,15 @@
   //!> temporary instance to hold return data
   VMGData returndata;
+
+  //!> number of particles
+  size_t num_particles;
+  //!> forces array
+  double *f;
+  //!> position array
+  double *x;
+  //!> potential array
+  double *p;
+  //!> charge array
+  double *q;
 
 private:
