Index: ThirdParty/vmg/src/units/particle/bspline.cpp
===================================================================
--- ThirdParty/vmg/src/units/particle/bspline.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/units/particle/bspline.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,247 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   bspline.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Nov 21 13:27:22 2011
+ *
+ * @brief  B-Splines for molecular dynamics.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#ifndef BSPLINE_DEGREE
+#error BSPLINE_DEGREE not defined.
+#endif
+
+#if BSPLINE_DEGREE < 3 || BSPLINE_DEGREE > 6
+#error Please choose a B-Spline degree between 3 and 6
+#endif
+
+#include "base/helper.hpp"
+#include "base/math.hpp"
+#include "units/particle/bspline.hpp"
+
+#define POW(x,y) Helper::pow(x,y)
+
+using namespace VMG;
+
+Particle::BSpline::BSpline(const int& near_field_cells, const vmg_float& h) :
+  spline_nom((BSPLINE_DEGREE+1)/2),
+  spline_denom((BSPLINE_DEGREE+1)/2),
+  potential_nom((BSPLINE_DEGREE+1)/2+1),
+  potential_denom((BSPLINE_DEGREE+1)/2+1),
+  field_nom((BSPLINE_DEGREE+1)/2),
+  field_denom((BSPLINE_DEGREE+1)/2),
+  intervals((BSPLINE_DEGREE+1)/2),
+  R(near_field_cells*h),
+  near_field_cells(near_field_cells)
+{
+  for (unsigned int i=0; i<intervals.size(); ++i)
+    intervals[i] = R * ( -1.0 + 2.0 / static_cast<vmg_float>(BSPLINE_DEGREE) * (i + BSPLINE_DEGREE/2 + 1));
+
+#if BSPLINE_DEGREE == 3
+
+  spline_nom[0] = Polynomial(2, 81.0*POW(R,2), 0.0, -243.0);
+  spline_nom[1] = Polynomial(2, 243.0*POW(R,2), -486.0*R, 243.0);
+
+  spline_denom[0] = Polynomial(0, 16.0 * Math::pi * POW(R,5));
+  spline_denom[1] = Polynomial(0, 32.0 * Math::pi * POW(R,5));
+
+  potential_nom[0] = Polynomial(5,
+				0.0,
+				-195.0*POW(R,4),
+				0.0,
+				270.0*POW(R,2),
+				0.0,
+				-243.0);
+
+  potential_nom[1] = Polynomial(5,
+				2.0     * POW(R,5),
+				-405.0  * POW(R,4),
+				0.0,
+				810.0   * POW(R,2),
+				-810.0  * R,
+				243.0);
+
+  potential_nom[2] = Polynomial(0, -1.0);
+
+  potential_denom[0] = Polynomial(0, 80.0*POW(R,5));
+  potential_denom[1] = Polynomial(0, 160.0*POW(R,5));
+  potential_denom[2] = Polynomial(0, 1.0);
+
+  field_nom[0] = Polynomial(5,
+			    20.0 * POW(R,5),
+			    0.0,
+			    0.0,
+			    -135.0 * POW(R,2),
+			    0.0,
+			    243.0);
+
+  field_nom[1] = Polynomial(5,
+			    81.0 * POW(R,5),
+			    0.0,
+			    0.0,
+			    -810.0 * POW(R,2),
+			    1215.0 * R,
+			    -486.0);
+
+  field_denom[0] = Polynomial(3,
+			      0.0,
+			      0.0,
+			      0.0,
+			      20.0 * POW(R,5));
+
+  field_denom[1] = Polynomial(3,
+			      0.0,
+			      0.0,
+			      0.0,
+			      80.0 * POW(R,5));
+
+  antid = 39.0 / (16.0 * R);
+
+#elif BSPLINE_DEGREE == 6
+
+  spline_nom[0] = Polynomial(5,
+                             297.0     * POW(R,5),
+                             0.0,
+                             -2430.0   * POW(R,3),
+                             0.0,
+                             10935.0   * R,
+                             -10935.0);
+
+  spline_nom[1] = Polynomial(5,
+                             459.0    * POW(R,5),
+                             2025.0   * POW(R,4),
+                             -17010.0 * POW(R,3),
+                             36450.0  * POW(R,2),
+                             -32805.0 * R,
+                             10935.0);
+
+  spline_nom[2] = Polynomial(5,
+                             2187.0   * POW(R,5),
+                             -10935.0 * POW(R,4),
+                             21870.0  * POW(R,3),
+                             -21870.0 * POW(R,2),
+                             10935.0  * R,
+                             -2187.0);
+
+  spline_denom[0] = Polynomial(0, 20.0 * Math::pi * POW(R,8));
+  spline_denom[1] = Polynomial(0, 40.0 * Math::pi * POW(R,8));
+  spline_denom[2] = Polynomial(0, 40.0 * Math::pi * POW(R,8));
+
+  potential_nom[0] = Polynomial(8,
+                                0.0,
+                                -956.0   * POW(R,7),
+                                0.0,
+                                2772.0   * POW(R,5),
+                                0.0,
+                                -6804.0  * POW(R,3),
+                                0.0,
+                                14580.0  * R,
+                                -10935.0);
+
+  potential_nom[1] = Polynomial(8,
+                                -5.0      * POW(R,8),
+                                -5676.0   * POW(R,7),
+                                0.0,
+                                12852.0   * POW(R,5),
+                                28350.0   * POW(R,4),
+                                -142884.0 * POW(R,3),
+                                204120.0  * POW(R,2),
+                                -131220.0 * R,
+                                32805.0);
+
+  potential_nom[2] = Polynomial(8,
+                                169.0    * POW(R,8),
+                                -2916.0  * POW(R,7),
+                                0.0,
+                                20412.0  * POW(R,5),
+                                -51030.0 * POW(R,4),
+                                61236.0  * POW(R,3),
+                                -40824.0 * POW(R,2),
+                                14580.0  * R,
+                                -2187.0);
+
+  potential_nom[3]  = Polynomial(0, -1.0);
+
+  potential_denom[0] = Polynomial(0, 280.0 * POW(R,8));
+  potential_denom[1] = Polynomial(0, 1680.0 * POW(R,8));
+  potential_denom[2] = Polynomial(0, 560.0 * POW(R,8));
+  potential_denom[3] = Polynomial(0, 0.0);
+
+  field_nom[0] = Polynomial(8,
+			    280.0    * POW(R,8),
+			    0.0,
+			    0.0,
+			    -5544.0  * POW(R,5),
+			    0.0,
+			    27216.0  * POW(R,3),
+			    0.0,
+			    -87480.0 * R,
+			    76545.0);
+
+  field_nom[1] = Polynomial(8,
+			    1675.0 * POW(R,8),
+			    0.0,
+			    0.0,
+			    -25704.0 * POW(R,5),
+			    -85050.0 * POW(R,4),
+			    571536.0 * POW(R,3),
+			    -1020600.0 * POW(R,2),
+			    787320.0 * R,
+			    -229635.0);
+
+  field_nom[2] = Polynomial(8,
+			    729.0 * POW(R,8),
+			    0.0,
+			    0.0,
+			    -40824.0 * POW(R,5),
+			    153090.0 * POW(R,4),
+			    -244944.0 * POW(R,3),
+			    204120.0 * POW(R,2),
+			    -87480.0 * R,
+			    15309.0);
+
+  field_denom[0] = Polynomial(3,
+			      0.0,
+			      0.0,
+			      0.0,
+			      280.0 * POW(R,8));
+  field_denom[1] = Polynomial(3,
+			      0.0,
+			      0.0,
+			      0.0,
+			      1680.0 * POW(R,8));
+
+  field_denom[2] = Polynomial(3,
+			      0.0,
+			      0.0,
+			      0.0,
+			      560.0 * POW(R,8));
+
+  antid = 239.0 / (70.0 * R);
+
+#else
+#error B-Spline degree not supported. Choose 3 or 6.
+#endif /* BSPLINE_DEGREE */
+}
Index: ThirdParty/vmg/src/units/particle/bspline.hpp
===================================================================
--- ThirdParty/vmg/src/units/particle/bspline.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/units/particle/bspline.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,313 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   bspline.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Nov 21 13:27:22 2011
+ *
+ * @brief  B-Splines for molecular dynamics.
+ *
+ */
+
+#ifndef BSPLINE_HPP_
+#define BSPLINE_HPP_
+
+#include <cmath>
+
+#include "base/helper.hpp"
+#include "base/index.hpp"
+#include "base/polynomial.hpp"
+#include "base/vector.hpp"
+#include "grid/grid.hpp"
+#include "units/particle/particle.hpp"
+
+namespace VMG
+{
+
+namespace Particle
+{
+
+class BSpline
+{
+public:
+  BSpline(const int& near_field_cells, const vmg_float& h);
+
+  vmg_float EvaluateSpline(const vmg_float& val) const
+  {
+    for (unsigned int i=0; i<intervals.size(); ++i)
+      if (val < intervals[i])
+        return spline_nom[i](val) / spline_denom[i](val);
+    return 0.0;
+  }
+
+  void SetSpline(Grid& grid, const Particle& p) const
+  {
+    assert(p.Pos().X() >= grid.Extent().Begin().X() && p.Pos().X() < grid.Extent().End().X());
+    assert(p.Pos().Y() >= grid.Extent().Begin().Y() && p.Pos().Y() < grid.Extent().End().Y());
+    assert(p.Pos().Z() >= grid.Extent().Begin().Z() && p.Pos().Z() < grid.Extent().End().Z());
+
+    vmg_float *vals = new vmg_float[Helper::intpow(2*near_field_cells+1,3)];
+
+    vmg_float temp_val;
+    vmg_float int_val = 0.0;
+    int c = 0;
+
+    const int index_global_x = grid.Global().GlobalBegin().X() + std::floor((p.Pos().X() - grid.Extent().Begin().X()) / grid.Extent().MeshWidth().X());
+    const int index_global_y = grid.Global().GlobalBegin().Y() + std::floor((p.Pos().Y() - grid.Extent().Begin().Y()) / grid.Extent().MeshWidth().Y());
+    const int index_global_z = grid.Global().GlobalBegin().Z() + std::floor((p.Pos().Z() - grid.Extent().Begin().Z()) / grid.Extent().MeshWidth().Z());
+
+    assert(index_global_x >= grid.Global().LocalBegin().X() && index_global_x < grid.Global().LocalEnd().X());
+    assert(index_global_y >= grid.Global().LocalBegin().Y() && index_global_y < grid.Global().LocalEnd().Y());
+    assert(index_global_z >= grid.Global().LocalBegin().Z() && index_global_z < grid.Global().LocalEnd().Z());
+
+    const int index_local_x = index_global_x - grid.Global().LocalBegin().X() + grid.Local().HaloSize1().X();
+    const int index_local_y = index_global_y - grid.Global().LocalBegin().Y() + grid.Local().HaloSize1().Y();
+    const int index_local_z = index_global_z - grid.Global().LocalBegin().Z() + grid.Local().HaloSize1().Z();
+
+    assert(index_local_x >= grid.Local().Begin().X() && index_local_x < grid.Local().End().X());
+    assert(index_local_y >= grid.Local().Begin().Y() && index_local_y < grid.Local().End().Y());
+    assert(index_local_z >= grid.Local().Begin().Z() && index_local_z < grid.Local().End().Z());
+
+    const vmg_float pos_beg_x =  grid.Extent().Begin().X() + grid.Extent().MeshWidth().X() * (index_global_x - grid.Global().GlobalBegin().X() - near_field_cells) - p.Pos().X();
+    const vmg_float pos_beg_y =  grid.Extent().Begin().Y() + grid.Extent().MeshWidth().Y() * (index_global_y - grid.Global().GlobalBegin().Y() - near_field_cells) - p.Pos().Y();
+    const vmg_float pos_beg_z =  grid.Extent().Begin().Z() + grid.Extent().MeshWidth().Z() * (index_global_z - grid.Global().GlobalBegin().Z() - near_field_cells) - p.Pos().Z();
+
+    const vmg_float& h_x = grid.Extent().MeshWidth().X();
+    const vmg_float& h_y = grid.Extent().MeshWidth().Y();
+    const vmg_float& h_z = grid.Extent().MeshWidth().Z();
+
+    // Iterate over all grid points which lie in the support of the interpolating B-Spline
+    vmg_float dir_x = pos_beg_x;
+    for (int i=-1*near_field_cells; i<=near_field_cells; ++i) {
+      vmg_float dir_y = pos_beg_y;
+      for (int j=-1*near_field_cells; j<=near_field_cells; ++j) {
+        vmg_float dir_z = pos_beg_z;
+        for (int k=-1*near_field_cells; k<=near_field_cells; ++k) {
+
+          // Compute distance from grid point to particle
+          temp_val = EvaluateSpline(std::sqrt(dir_x*dir_x+dir_y*dir_y+dir_z*dir_z));
+          vals[c++] = temp_val * p.Charge();
+          int_val += temp_val;
+
+          dir_z += h_z;
+        }
+        dir_y += h_y;
+      }
+      dir_x += h_x;
+    }
+
+    // Reciprocal value of the numerically integrated spline
+    int_val = 1.0 / (int_val * h_x * h_y * h_z);
+
+	c = 0;
+    for (int i=-1*near_field_cells; i<=near_field_cells; ++i)
+      for (int j=-1*near_field_cells; j<=near_field_cells; ++j)
+        for (int k=-1*near_field_cells; k<=near_field_cells; ++k)
+          grid(index_local_x + i,
+              index_local_y + j,
+              index_local_z + k) += vals[c++] * int_val;
+    assert( c == Helper::intpow(2*near_field_cells+1,3) );
+
+    delete [] vals;
+  }
+
+  void changeGridBySelfInducedPotential(Grid& grid, Particle& p, const vmg_float &sign) const
+  {
+	assert(p.Pos().X() >= grid.Extent().Begin().X() && p.Pos().X() < grid.Extent().End().X());
+	assert(p.Pos().Y() >= grid.Extent().Begin().Y() && p.Pos().Y() < grid.Extent().End().Y());
+	assert(p.Pos().Z() >= grid.Extent().Begin().Z() && p.Pos().Z() < grid.Extent().End().Z());
+
+	const int index_global_x = grid.Global().GlobalBegin().X() + std::floor((p.Pos().X() - grid.Extent().Begin().X()) / grid.Extent().MeshWidth().X());
+	const int index_global_y = grid.Global().GlobalBegin().Y() + std::floor((p.Pos().Y() - grid.Extent().Begin().Y()) / grid.Extent().MeshWidth().Y());
+	const int index_global_z = grid.Global().GlobalBegin().Z() + std::floor((p.Pos().Z() - grid.Extent().Begin().Z()) / grid.Extent().MeshWidth().Z());
+
+	assert(index_global_x >= grid.Global().LocalBegin().X() && index_global_x < grid.Global().LocalEnd().X());
+	assert(index_global_y >= grid.Global().LocalBegin().Y() && index_global_y < grid.Global().LocalEnd().Y());
+	assert(index_global_z >= grid.Global().LocalBegin().Z() && index_global_z < grid.Global().LocalEnd().Z());
+
+	const int index_local_x = index_global_x - grid.Global().LocalBegin().X() + grid.Local().HaloSize1().X();
+	const int index_local_y = index_global_y - grid.Global().LocalBegin().Y() + grid.Local().HaloSize1().Y();
+	const int index_local_z = index_global_z - grid.Global().LocalBegin().Z() + grid.Local().HaloSize1().Z();
+
+	assert(index_local_x >= grid.Local().Begin().X() && index_local_x < grid.Local().End().X());
+	assert(index_local_y >= grid.Local().Begin().Y() && index_local_y < grid.Local().End().Y());
+	assert(index_local_z >= grid.Local().Begin().Z() && index_local_z < grid.Local().End().Z());
+
+	const vmg_float pos_beg_x =  grid.Extent().Begin().X() + grid.Extent().MeshWidth().X() * (index_global_x - grid.Global().GlobalBegin().X() - near_field_cells) - p.Pos().X();
+	const vmg_float pos_beg_y =  grid.Extent().Begin().Y() + grid.Extent().MeshWidth().Y() * (index_global_y - grid.Global().GlobalBegin().Y() - near_field_cells) - p.Pos().Y();
+	const vmg_float pos_beg_z =  grid.Extent().Begin().Z() + grid.Extent().MeshWidth().Z() * (index_global_z - grid.Global().GlobalBegin().Z() - near_field_cells) - p.Pos().Z();
+
+	const vmg_float& h_x = grid.Extent().MeshWidth().X();
+	const vmg_float& h_y = grid.Extent().MeshWidth().Y();
+	const vmg_float& h_z = grid.Extent().MeshWidth().Z();
+
+	vmg_float dir_x = pos_beg_x;
+    for (int i=-1*near_field_cells; i<=near_field_cells; ++i) {
+  	  vmg_float dir_y = pos_beg_y;
+      for (int j=-1*near_field_cells; j<=near_field_cells; ++j) {
+  		vmg_float dir_z = pos_beg_z;
+        for (int k=-1*near_field_cells; k<=near_field_cells; ++k) {
+  		  const double length_sq = dir_x*dir_x+dir_y*dir_y+dir_z*dir_z;
+		  const double length = std::sqrt(length_sq);
+          grid(index_local_x + i,
+              index_local_y + j,
+              index_local_z + k) += sign * p.Charge() * EvaluatePotential(length);
+
+          dir_z += h_z;
+		}
+	    dir_y += h_y;
+	  }
+	  dir_x += h_x;
+	}
+  }
+
+  void SubtractSelfInducedForces(const Grid& grid, Particle& p) const
+  {
+	assert(p.Pos().X() >= grid.Extent().Begin().X() && p.Pos().X() < grid.Extent().End().X());
+	assert(p.Pos().Y() >= grid.Extent().Begin().Y() && p.Pos().Y() < grid.Extent().End().Y());
+	assert(p.Pos().Z() >= grid.Extent().Begin().Z() && p.Pos().Z() < grid.Extent().End().Z());
+
+	vmg_float temp_val = 0.;
+
+	const int index_global_x = grid.Global().GlobalBegin().X() + std::floor((p.Pos().X() - grid.Extent().Begin().X()) / grid.Extent().MeshWidth().X());
+	const int index_global_y = grid.Global().GlobalBegin().Y() + std::floor((p.Pos().Y() - grid.Extent().Begin().Y()) / grid.Extent().MeshWidth().Y());
+	const int index_global_z = grid.Global().GlobalBegin().Z() + std::floor((p.Pos().Z() - grid.Extent().Begin().Z()) / grid.Extent().MeshWidth().Z());
+
+	assert(index_global_x >= grid.Global().LocalBegin().X() && index_global_x < grid.Global().LocalEnd().X());
+	assert(index_global_y >= grid.Global().LocalBegin().Y() && index_global_y < grid.Global().LocalEnd().Y());
+	assert(index_global_z >= grid.Global().LocalBegin().Z() && index_global_z < grid.Global().LocalEnd().Z());
+
+	const int index_local_x = index_global_x - grid.Global().LocalBegin().X() + grid.Local().HaloSize1().X();
+	const int index_local_y = index_global_y - grid.Global().LocalBegin().Y() + grid.Local().HaloSize1().Y();
+	const int index_local_z = index_global_z - grid.Global().LocalBegin().Z() + grid.Local().HaloSize1().Z();
+
+	assert(index_local_x >= grid.Local().Begin().X() && index_local_x < grid.Local().End().X());
+	assert(index_local_y >= grid.Local().Begin().Y() && index_local_y < grid.Local().End().Y());
+	assert(index_local_z >= grid.Local().Begin().Z() && index_local_z < grid.Local().End().Z());
+
+	const vmg_float pos_beg_x =  grid.Extent().Begin().X() + grid.Extent().MeshWidth().X() * (index_global_x - grid.Global().GlobalBegin().X() - near_field_cells) - p.Pos().X();
+	const vmg_float pos_beg_y =  grid.Extent().Begin().Y() + grid.Extent().MeshWidth().Y() * (index_global_y - grid.Global().GlobalBegin().Y() - near_field_cells) - p.Pos().Y();
+	const vmg_float pos_beg_z =  grid.Extent().Begin().Z() + grid.Extent().MeshWidth().Z() * (index_global_z - grid.Global().GlobalBegin().Z() - near_field_cells) - p.Pos().Z();
+
+	const vmg_float& h_x = grid.Extent().MeshWidth().X();
+	const vmg_float& h_y = grid.Extent().MeshWidth().Y();
+	const vmg_float& h_z = grid.Extent().MeshWidth().Z();
+
+	vmg_float int_val = 0.;
+    vmg_float dir_x = pos_beg_x;
+    for (int i=-1*near_field_cells; i<=near_field_cells; ++i) {
+      vmg_float dir_y = pos_beg_y;
+      for (int j=-1*near_field_cells; j<=near_field_cells; ++j) {
+        vmg_float dir_z = pos_beg_z;
+        for (int k=-1*near_field_cells; k<=near_field_cells; ++k) {
+
+  		  const double length_sq = dir_x*dir_x+dir_y*dir_y+dir_z*dir_z;
+  		  if (fabs(length_sq) > std::numeric_limits<double>::epsilon()) {
+			  // Compute distance from grid point to particle
+			  int_val += EvaluateSpline(std::sqrt(length_sq));
+  		  }
+
+          dir_z += h_z;
+        }
+        dir_y += h_y;
+      }
+      dir_x += h_x;
+    }
+    // Reciprocal value of the numerically integrated spline
+    if (fabs(int_val) > std::numeric_limits<double>::epsilon())
+    	int_val = 1. / (int_val * h_x * h_y * h_z);
+    else
+    	int_val = 1. / (h_x * h_y * h_z);
+
+    // Iterate over all grid points which lie in the support of the interpolating B-Spline
+    vmg_float test_int_val = 0.;
+	dir_x = pos_beg_x;
+	for (int i=-1*near_field_cells; i<=near_field_cells; ++i) {
+	  vmg_float dir_y = pos_beg_y;
+	  for (int j=-1*near_field_cells; j<=near_field_cells; ++j) {
+		vmg_float dir_z = pos_beg_z;
+		for (int k=-1*near_field_cells; k<=near_field_cells; ++k) {
+
+		  // Compute distance from grid point to particle
+		  const double length_sq = dir_x*dir_x+dir_y*dir_y+dir_z*dir_z;
+		  if (fabs(length_sq) > std::numeric_limits<double>::epsilon()) {
+		    const double length = std::sqrt(length_sq);
+
+		    temp_val = h_x * h_y * h_z * int_val * EvaluateSpline(length);
+		    test_int_val += temp_val;
+		    p.Field()[0] -= p.Charge()* temp_val * dir_x / (length_sq*length);
+		    p.Field()[1] -= p.Charge()* temp_val * dir_y / (length_sq*length);
+		    p.Field()[2] -= p.Charge()* temp_val * dir_z / (length_sq*length);
+		  } else {
+			  std::cerr << "Value very small " << length_sq << "=("
+					  << dir_x << ")^2+(" << dir_y << ")^2+(" << dir_z << ")^2 for particle at ("
+					  << p.Pos().X() << ","
+					  << p.Pos().Y() << ","
+					  << p.Pos().Z() << ")"
+					  << "\n";
+		  }
+		  dir_z += h_z;
+		}
+		dir_y += h_y;
+	  }
+	  dir_x += h_x;
+	}
+	if ( fabs(test_int_val -1.) > std::numeric_limits<double>::epsilon()*1e4 ) {
+		std::cerr << "Integrated spline value should be 1 but is " << test_int_val << std::endl;
+		assert(0);
+	}
+  }
+
+  vmg_float EvaluatePotential(const vmg_float& val) const
+  {
+    for (unsigned int i=0; i<intervals.size(); ++i)
+      if (val < intervals[i])
+        return potential_nom[i](val) / potential_denom[i](val);
+    return potential_nom.back()(val) / potential_denom.back()(val);
+  }
+
+  vmg_float EvaluateField(const vmg_float& val) const
+  {
+    for (unsigned int i=0; i<intervals.size(); ++i)
+      if (val < intervals[i])
+        return field_nom[i](val) / field_denom[i](val);
+    return 0.0;
+  }
+
+  const vmg_float& GetAntiDerivativeAtZero() const
+  {
+    return antid;
+  }
+
+private:
+  std::vector<Polynomial> spline_nom, spline_denom;
+  std::vector<Polynomial> potential_nom, potential_denom;
+  std::vector<Polynomial> field_nom, field_denom;
+  vmg_float antid;
+  std::vector<vmg_float> intervals;
+
+  const vmg_float R;
+  const int near_field_cells;
+};
+
+}
+
+}
+
+#endif /* BSPLINE_HPP_ */
Index: ThirdParty/vmg/src/units/particle/comm_mpi_particle.cpp
===================================================================
--- ThirdParty/vmg/src/units/particle/comm_mpi_particle.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/units/particle/comm_mpi_particle.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,291 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   comm_mpi_particle.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Tue May 8 15:27:06 2012
+ *
+ * @brief  Class for MPI-based particle-related communication.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#ifdef HAVE_MPI
+
+#include <mpi.h>
+#ifdef HAVE_MARMOT
+#include <enhancempicalls.h>
+#include <sourceinfompicalls.h>
+#endif
+
+#include "units/particle/comm_mpi_particle.hpp"
+#include "units/particle/linked_cell_list.hpp"
+#include "units/particle/particle.hpp"
+
+using namespace VMG;
+
+void Particle::CommMPI::CommParticles(const Grid& grid, std::list<Particle>& particles)
+{
+  Factory& factory = MG::GetFactory();
+
+  const vmg_int& num_particles_local = factory.GetObjectStorageVal<vmg_int>("PARTICLE_NUM_LOCAL");
+  vmg_float* x = factory.GetObjectStorageArray<vmg_float>("PARTICLE_POS_ARRAY");
+  vmg_float* q = factory.GetObjectStorageArray<vmg_float>("PARTICLE_CHARGE_ARRAY");
+
+  int rank, size;
+  MPI_Comm_rank(comm_global, &rank);
+  MPI_Comm_size(comm_global, &size);
+
+#ifndef VMG_ONE_SIDED
+  vmg_int* receiver;
+  if (!factory.TestObject("PARTICLE_RECEIVER_ARRAY"))
+    new ObjectStorageArray<vmg_int>("PARTICLE_RECEIVER_ARRAY", size);
+  receiver = factory.GetObjectStorageArray<vmg_int>("PARTICLE_RECEIVER_ARRAY");
+#endif
+
+  Index index;
+  std::vector<int> global_extent(6*size);
+  std::vector<int> send_sizes(size);
+  std::vector<int> recv_sizes(size);
+  std::vector<Index> begin_remote(size);
+  std::vector<Index> end_remote(size);
+  std::vector< std::vector<vmg_float> > send_buffer_x(size);
+  std::vector< std::vector<vmg_float> > send_buffer_q(size);
+  std::vector< std::vector<vmg_int> > send_buffer_ind(size);
+  std::vector< std::vector<vmg_float> > recv_buffer_x(size);
+  std::vector< std::vector<vmg_float> > recv_buffer_q(size);
+  std::vector< std::vector<vmg_int> > recv_buffer_ind(size);
+
+  std::memcpy(&global_extent[6*rank], grid.Global().LocalBegin().vec(), 3*sizeof(int));
+  std::memcpy(&global_extent[6*rank+3], grid.Global().LocalEnd().vec(), 3*sizeof(int));
+
+  MPI_Allgather(MPI_IN_PLACE, 6, MPI_INT, &global_extent.front(), 6, MPI_INT, comm_global);
+
+  for (int i=0; i<size; ++i) {
+    begin_remote[i] = static_cast<Index>(&global_extent[6*i]);
+    end_remote[i] = static_cast<Index>(&global_extent[6*i+3]);
+  }
+
+  for (int i=0; i<num_particles_local; ++i) {
+    index = grid.Global().GlobalBegin() + static_cast<Index>(((Vector(&x[3*i]) - grid.Extent().Begin()) / grid.Extent().MeshWidth()).Floor());
+    for (int j=0; j<size; ++j)
+      if (index.IsInBounds(begin_remote[j], end_remote[j])) {
+        send_buffer_x[j].push_back(x[3*i+0]);
+        send_buffer_x[j].push_back(x[3*i+1]);
+        send_buffer_x[j].push_back(x[3*i+2]);
+        send_buffer_q[j].push_back(q[i]);
+        send_buffer_ind[j].push_back(i);
+        break;
+      }
+  }
+
+  /*
+   * Communicate which process gets how many particles
+   */
+  for (int i=0; i<size; ++i)
+    send_sizes[i] = send_buffer_q[i].size();
+
+  MPI_Alltoall(&send_sizes.front(), 1, MPI_INT, &recv_sizes.front(), 1, MPI_INT, comm_global);
+
+  assert(RequestsPending() == 0);
+
+  /*
+   * Send particles
+   */
+  for (int i=0; i<size; ++i) {
+
+    if (!send_buffer_q[i].empty()) {
+      MPI_Isend(&send_buffer_x[i].front(), send_buffer_x[i].size(), MPI_DOUBLE, i, 0, comm_global, &Request());
+      MPI_Isend(&send_buffer_q[i].front(), send_buffer_q[i].size(), MPI_DOUBLE, i, 1, comm_global, &Request());
+      MPI_Isend(&send_buffer_ind[i].front(), send_buffer_ind[i].size(), MPI_INT, i, 2, comm_global, &Request());
+    }
+
+#ifndef VMG_ONE_SIDED
+    receiver[i] = send_buffer_q[i].size();
+#endif
+  }
+
+  /*
+   * Receive particles
+   */
+  for (int i=0; i<size; ++i) {
+
+    if (recv_sizes[i] > 0) {
+
+      recv_buffer_x[i].resize(3*recv_sizes[i]);
+      recv_buffer_q[i].resize(recv_sizes[i]);
+      recv_buffer_ind[i].resize(recv_sizes[i]);
+
+      MPI_Irecv(&recv_buffer_x[i].front(), 3*recv_sizes[i], MPI_DOUBLE, i, 0, comm_global, &Request());
+      MPI_Irecv(&recv_buffer_q[i].front(), recv_sizes[i], MPI_DOUBLE, i, 1, comm_global, &Request());
+      MPI_Irecv(&recv_buffer_ind[i].front(), recv_sizes[i], MPI_INT, i, 2, comm_global, &Request());
+
+    }
+
+  }
+
+  WaitAll();
+
+  particles.clear();
+
+  for (int i=0; i<size; ++i)
+    for (int j=0; j<recv_sizes[i]; ++j)
+      particles.push_back(Particle(&recv_buffer_x[i][3*j], recv_buffer_q[i][j], 0.0, 0.0, i, recv_buffer_ind[i][j]));
+}
+
+void Particle::CommMPI::CommParticlesBack(std::list<Particle>& particles)
+{
+  std::list<Particle>::iterator iter;
+
+#ifdef VMG_ONE_SIDED
+  if (!win_created) {
+    vmg_float* p = MG::GetFactory().GetObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY");
+    const vmg_int& num_particles_local = MG::GetFactory().GetObjectStorageVal<vmg_int>("PARTICLE_NUM_LOCAL");
+    MPI_Win_create(p, num_particles_local*sizeof(vmg_float), sizeof(vmg_float), info, comm_global, &win);
+    win_created = true;
+  }
+
+  MPI_Win_fence(MPI_MODE_NOPRECEDE, win);
+
+  for (iter=particles.begin(); iter!=particles.end(); ++iter)
+    MPI_Put(&iter->Pot(), 1, MPI_DOUBLE, iter->Rank(), iter->Index(), 1, MPI_DOUBLE, win);
+
+  MPI_Win_fence(MPI_MODE_NOSTORE | MPI_MODE_NOSUCCEED, win);
+#else
+  int rank, size;
+  MPI_Comm_rank(comm_global, &rank);
+  MPI_Comm_size(comm_global, &size);
+
+  std::vector< std::vector<vmg_float> > send_buffer_float(size);
+  std::vector< std::vector<vmg_float> > recv_buffer_float(size);
+  std::vector< std::vector<vmg_int> > send_buffer_index(size);
+  std::vector< std::vector<vmg_int> > recv_buffer_index(size);
+
+  vmg_int* size_receive = MG::GetFactory().GetObjectStorageArray<vmg_int>("PARTICLE_RECEIVER_ARRAY");
+  vmg_float* p = MG::GetFactory().GetObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY");
+  vmg_float* f = MG::GetFactory().GetObjectStorageArray<vmg_float>("PARTICLE_FIELD_ARRAY");
+
+  // Build send buffer
+  for (iter=particles.begin(); iter!=particles.end(); ++iter) {
+    send_buffer_float[iter->Rank()].push_back(iter->Pot());
+    send_buffer_float[iter->Rank()].push_back(iter->Field()[0]);
+    send_buffer_float[iter->Rank()].push_back(iter->Field()[1]);
+    send_buffer_float[iter->Rank()].push_back(iter->Field()[2]);
+    send_buffer_index[iter->Rank()].push_back(iter->Index());
+  }
+
+  // Send potentials
+  for (int i=0; i<size; ++i) {
+    if (!send_buffer_float[i].empty()) {
+      MPI_Isend(&send_buffer_float[i].front(), send_buffer_float[i].size(), MPI_DOUBLE, i, 699+rank, comm_global, &Request());
+      MPI_Isend(&send_buffer_index[i].front(), send_buffer_index[i].size(), MPI_INT, i, 32111+rank, comm_global, &Request());
+    }
+  }
+
+  //Receive potentials
+  for (int i=0; i<size; ++i) {
+    if (size_receive[i] > 0) {
+      recv_buffer_float[i].resize(4*size_receive[i]);
+      recv_buffer_index[i].resize(size_receive[i]);
+      MPI_Irecv(&recv_buffer_float[i].front(), 4*size_receive[i], MPI_DOUBLE, i, 699+i, comm_global, &Request());
+      MPI_Irecv(&recv_buffer_index[i].front(), size_receive[i], MPI_INT, i, 32111+i, comm_global, &Request());
+    }
+  }
+
+  WaitAll();
+
+  // Add potential values
+  for (int i=0; i<size; ++i)
+    for (vmg_int j=0; j<size_receive[i]; ++j) {
+      p[recv_buffer_index[i][j]] = recv_buffer_float[i][4*j];
+      std::memcpy(&f[3*recv_buffer_index[i][j]], &recv_buffer_float[i][4*j+1], 3*sizeof(vmg_float));
+    }
+#endif
+
+}
+
+void Particle::CommMPI::CommLCListToGhosts(LinkedCellList& lc)
+{
+  VMG::MPI::DatatypesLocal types(lc, comm_global, false);
+  std::vector<int> send_size(types.NB().size());
+  vmg_int recv_size;
+  std::list<Particle*>::iterator iter;
+  Index ind;
+  Vector offset;
+
+  const Vector halo_length = lc.Local().HaloSize1() * lc.Extent().MeshWidth();
+
+  lc.ClearHalo();
+
+  for (unsigned int i=0; i<types.NB().size(); ++i)
+    if (types.NB()[i].Feasible()) {
+
+      for (int j=0; j<3; ++j)
+        if ((types.Offset()[i][j] < 0 && lc.Global().LocalBegin()[j] == lc.Global().GlobalBegin()[j]) ||
+            (types.Offset()[i][j] > 0 && lc.Global().LocalEnd()[j] == lc.Global().GlobalEnd()[j]))
+          offset[j] = -1.0 * types.Offset()[i][j] * lc.Extent().Size()[j];
+        else
+          offset[j] = 0.0;
+
+      for (ind.X() = types.NB()[i].Starts().X(); ind.X() < types.NB()[i].Starts().X()+types.NB()[i].Subsizes().X(); ++ind.X())
+        for (ind.Y() = types.NB()[i].Starts().Y(); ind.Y() < types.NB()[i].Starts().Y()+types.NB()[i].Subsizes().Y(); ++ind.Y())
+          for (ind.Z() = types.NB()[i].Starts().Z(); ind.Z() < types.NB()[i].Starts().Z()+types.NB()[i].Subsizes().Z(); ++ind.Z())
+            for (iter=lc(ind).begin(); iter!=lc(ind).end(); ++iter) {
+
+              for (int j=0; j<3; ++j)
+                types.NB()[i].Buffer().push_back((*iter)->Pos()[j] + offset[j]);
+              types.NB()[i].Buffer().push_back((*iter)->Charge());
+
+              assert(lc.Extent().Begin().IsComponentwiseLessOrEqual((*iter)->Pos()));
+              assert(lc.Extent().End().IsComponentwiseGreaterOrEqual((*iter)->Pos()));
+              assert(lc.Extent().Begin().IsComponentwiseLessOrEqual((*iter)->Pos() + offset + halo_length));
+              assert(lc.Extent().End().IsComponentwiseGreaterOrEqual((*iter)->Pos() + offset - halo_length));
+            }
+
+      send_size[i] = types.NB()[i].Buffer().size();
+      MPI_Isend(&send_size[i], 1, MPI_INT, types.NB()[i].Rank(), 2048+types.NB()[i].TagSend(), comm_global, &Request());
+
+      if (send_size[i] > 0)
+        MPI_Isend(&types.NB()[i].Buffer().front(), send_size[i], MPI_DOUBLE,
+            types.NB()[i].Rank(), 4096+types.NB()[i].TagSend(),
+            comm_global, &Request());
+    }
+
+  for (unsigned int i=0; i<types.Halo().size(); ++i)
+    if (types.Halo()[i].Feasible()) {
+      MPI_Recv(&recv_size, 1, MPI_INT, types.Halo()[i].Rank(), 2048+types.Halo()[i].TagReceive(), comm_global, MPI_STATUS_IGNORE);
+      if (recv_size > 0) {
+        types.Halo()[i].Buffer().resize(recv_size);
+        MPI_Irecv(&types.Halo()[i].Buffer().front(), recv_size, MPI_DOUBLE,
+            types.Halo()[i].Rank(), 4096+types.Halo()[i].TagReceive(),
+            comm_global, &Request());
+      }
+    }
+
+  WaitAll();
+
+  for (unsigned int i=0; i<types.Halo().size(); ++i)
+    for (unsigned int j=0; j<types.Halo()[i].Buffer().size(); j+=4)
+      lc.AddParticleToHalo(&types.Halo()[i].Buffer()[j], types.Halo()[i].Buffer()[j+3]);
+}
+
+#endif /* HAVE_MPI */
Index: ThirdParty/vmg/src/units/particle/comm_mpi_particle.hpp
===================================================================
--- ThirdParty/vmg/src/units/particle/comm_mpi_particle.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/units/particle/comm_mpi_particle.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,57 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef COMM_MPI_PARTICLE_HPP_
+#define COMM_MPI_PARTICLE_HPP_
+
+#include <list>
+
+#include "comm/comm_mpi.hpp"
+#include "units/particle/particle.hpp"
+
+namespace VMG
+{
+
+namespace Particle
+{
+
+class LinkedCellList;
+
+class CommMPI : public VMG::CommMPI
+{
+public:
+  CommMPI(const Boundary& boundary, DomainDecomposition* domain_dec, const MPI_Comm& mpi_comm) :
+    VMG::CommMPI(boundary, domain_dec, mpi_comm)
+  {}
+
+  CommMPI(const Boundary& boundary, DomainDecomposition* domain_dec) :
+    VMG::CommMPI(boundary, domain_dec)
+  {}
+
+  virtual ~CommMPI() {}
+
+  void CommParticles(const Grid& grid, std::list<Particle>& particles);
+  void CommParticlesBack(std::list<Particle>& particles);
+  void CommLCListToGhosts(LinkedCellList& lc);
+};
+
+}
+
+}
+
+#endif /* COMM_MPI_PARTICLE_HPP_ */
Index: ThirdParty/vmg/src/units/particle/interface_fcs.cpp
===================================================================
--- ThirdParty/vmg/src/units/particle/interface_fcs.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/units/particle/interface_fcs.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,247 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   interface_fcs.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Apr 18 12:56:20 2011
+ *
+ * @brief  VMG::InterfaceFCS
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#ifndef HAVE_MPI
+#error MPI is needed to use the Scafacos interface.
+#endif
+
+#include <mpi.h>
+#ifdef HAVE_MARMOT
+#include <enhancempicalls.h>
+#include <sourceinfompicalls.h>
+#endif
+
+#include "base/object.hpp"
+#include "base/timer.hpp"
+#include "comm/domain_decomposition_mpi.hpp"
+#ifdef DEBUG
+#include "comm/mpi/error_handler.hpp"
+#endif
+#include "cycles/cycle_cs_periodic.hpp"
+#include "cycles/cycle_fas_dirichlet.hpp"
+#include "discretization/discretization_poisson_fd.hpp"
+#include "discretization/discretization_poisson_fv.hpp"
+#include "level/level_operator_cs.hpp"
+#include "level/level_operator_fas.hpp"
+#include "level/stencils.hpp"
+#include "smoother/gsrb_poisson_2.hpp"
+#include "smoother/gsrb_poisson_4.hpp"
+#include "solver/givens.hpp"
+#include "solver/solver_regular.hpp"
+#include "solver/solver_singular.hpp"
+#include "units/particle/comm_mpi_particle.hpp"
+#include "units/particle/interface_fcs.h"
+#include "units/particle/interface_particles.hpp"
+
+#include "cycles/cycle_cs_periodic_debug.hpp"
+//TODO: Remove Debug
+
+using namespace VMG;
+
+namespace VMGBackupSettings
+{
+  static vmg_int level = -1;
+  static vmg_int periodic[3] = {-1, -1, -1};
+  static vmg_int max_iter = -1;
+  static vmg_int smoothing_steps = -1;
+  static vmg_int cycle_type = -1;
+  static vmg_float precision = -1;
+  static vmg_float box_offset[3];
+  static vmg_float box_size = -1.0;
+  static vmg_int near_field_cells = -1;
+  static vmg_int interpolation_degree = -1;
+  static vmg_int discretization_order = -1;
+  static MPI_Comm mpi_comm;
+}
+
+static void VMG_fcs_init(vmg_int level, vmg_int* periodic,vmg_int max_iter,
+			 vmg_int smoothing_steps, vmg_int cycle_type, vmg_float precision,
+			 vmg_float* box_offset, vmg_float box_size,
+			 vmg_int near_field_cells, vmg_int interpolation_degree,
+                         vmg_int discretization_order, MPI_Comm mpi_comm)
+{
+  VMGBackupSettings::level = level;
+  std::memcpy(VMGBackupSettings::periodic, periodic, 3*sizeof(vmg_int));
+  VMGBackupSettings::max_iter = max_iter;
+  VMGBackupSettings::smoothing_steps = smoothing_steps;
+  VMGBackupSettings::cycle_type = cycle_type;
+  VMGBackupSettings::precision = precision;
+  std::memcpy(VMGBackupSettings::box_offset, box_offset, 3*sizeof(vmg_float));
+  VMGBackupSettings::box_size = box_size;
+  VMGBackupSettings::near_field_cells = near_field_cells;
+  VMGBackupSettings::interpolation_degree = interpolation_degree;
+  VMGBackupSettings::discretization_order = discretization_order;
+  VMGBackupSettings::mpi_comm = mpi_comm;
+
+#ifdef DEBUG
+  MPI_Errhandler mpiErrorHandler;
+  MPI_Comm_create_errhandler(VMG::MPI::ConvertToException, &mpiErrorHandler);
+  MPI_Comm_set_errhandler(mpi_comm, mpiErrorHandler);
+#endif
+
+  const Boundary boundary(periodic[0] ? Periodic : Open,
+			  periodic[1] ? Periodic : Open,
+			  periodic[2] ? Periodic : Open);
+
+  const bool singular = boundary[0] == Periodic &&
+                        boundary[1] == Periodic &&
+                        boundary[2] == Periodic;
+
+  /*
+   * Choose multigrid components
+   */
+  if (singular) {
+
+    new Particle::CommMPI(boundary, new DomainDecompositionMPI(), mpi_comm);
+    new DiscretizationPoissonFD(discretization_order);
+    new InterfaceParticles(boundary, 2, level, Vector(box_offset), box_size, near_field_cells);
+    new LevelOperatorCS(Stencils::RestrictionFullWeight, Stencils::InterpolationTrilinear);
+    new Givens<SolverSingular>();
+    new CycleCSPeriodic(cycle_type);
+
+  }else {
+
+    new Particle::CommMPI(boundary, new DomainDecompositionMPI(), mpi_comm);
+    new DiscretizationPoissonFV(discretization_order);
+    new InterfaceParticles(boundary, 2, level, Vector(box_offset), box_size, near_field_cells, 9, 1.6);
+    new LevelOperatorFAS(Stencils::RestrictionFullWeight, Stencils::Injection, Stencils::InterpolationTrilinear);
+    new Givens<SolverRegular>();
+    new CycleFASDirichlet(cycle_type);
+
+  }
+
+  /*
+   * Use Gauss-Seidel Red-Black ordering
+   */
+  if (discretization_order == 2)
+    new GaussSeidelRBPoisson2();
+  else
+    new GaussSeidelRBPoisson4();
+
+  /*
+   * Register required parameters
+   */
+  new ObjectStorage<int>("PRESMOOTHSTEPS", smoothing_steps);
+  new ObjectStorage<int>("POSTSMOOTHSTEPS", smoothing_steps);
+  new ObjectStorage<vmg_float>("PRECISION", precision);
+  new ObjectStorage<int>("MAX_ITERATION", max_iter);
+  new ObjectStorage<int>("PARTICLE_NEAR_FIELD_CELLS", near_field_cells);
+  new ObjectStorage<int>("PARTICLE_INTERPOLATION_DEGREE", interpolation_degree);
+
+  /*
+   * Post init
+   */
+  MG::PostInit();
+
+  /*
+   * Check whether the library is correctly initialized now.
+   */
+  MG::IsInitialized();
+}
+
+void VMG_fcs_setup(vmg_int level, vmg_int* periodic, vmg_int max_iter,
+		   vmg_int smoothing_steps, vmg_int cycle_type, vmg_float precision,
+		   vmg_float* box_offset, vmg_float box_size,
+		   vmg_int near_field_cells, vmg_int interpolation_degree,
+                   vmg_int discretization_order, MPI_Comm mpi_comm)
+{
+  if (VMGBackupSettings::level != level ||
+      VMGBackupSettings::periodic[0] != periodic[0] ||
+      VMGBackupSettings::periodic[1] != periodic[1] ||
+      VMGBackupSettings::periodic[2] != periodic[2] ||
+      VMGBackupSettings::max_iter != max_iter ||
+      VMGBackupSettings::smoothing_steps != smoothing_steps ||
+      VMGBackupSettings::cycle_type != cycle_type ||
+      VMGBackupSettings::precision != precision ||
+      VMGBackupSettings::box_offset[0] != box_offset[0] ||
+      VMGBackupSettings::box_offset[1] != box_offset[1] ||
+      VMGBackupSettings::box_offset[2] != box_offset[2] ||
+      VMGBackupSettings::box_size != box_size ||
+      VMGBackupSettings::near_field_cells != near_field_cells ||
+      VMGBackupSettings::interpolation_degree != interpolation_degree ||
+      VMGBackupSettings::discretization_order != discretization_order ||
+      VMGBackupSettings::mpi_comm != mpi_comm) {
+
+    VMG_fcs_destroy();
+    VMG_fcs_init(level, periodic, max_iter,
+		 smoothing_steps, cycle_type, precision,
+		 box_offset, box_size, near_field_cells,
+                 interpolation_degree, discretization_order,
+		 mpi_comm);
+
+  }
+}
+
+int VMG_fcs_check()
+{
+  const int& near_field_cells = MG::GetFactory().GetObjectStorageVal<int>("PARTICLE_NEAR_FIELD_CELLS");
+  const Multigrid& multigrid = *MG::GetRhs();
+  const Grid& grid = multigrid(multigrid.MaxLevel());
+
+  vmg_int error_code = 0;
+
+  if (!grid.Global().LocalSize().IsComponentwiseGreater(near_field_cells))
+    error_code = 1;
+
+  return MG::GetComm()->GlobalMax(error_code);
+}
+
+void VMG_fcs_run(vmg_float* x, vmg_float* q, vmg_float* p, vmg_float* f, vmg_int num_particles_local)
+{
+  /*
+   * Register parameters for later use.
+   */
+  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_local);
+
+  /*
+   * Start the multigrid solver
+   */
+  MG::Solve();
+
+  Timer::Print();
+}
+
+void VMG_fcs_print_timer()
+{
+  Timer::Print();
+}
+
+void VMG_fcs_destroy(void)
+{
+  /*
+   * Delete all data.
+   */
+  MG::Destroy();
+}
Index: ThirdParty/vmg/src/units/particle/interface_fcs.h
===================================================================
--- ThirdParty/vmg/src/units/particle/interface_fcs.h	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/units/particle/interface_fcs.h	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,49 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   interface_fcs.h
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Tue Apr 12 17:39:27 2011
+ *
+ * @brief  Scafacos C interface.
+ *
+ *
+ */
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void VMG_fcs_setup(fcs_int level, fcs_int* periodic, fcs_int max_iter,
+		   fcs_int smoothing_steps, fcs_int cycle_type, fcs_float precision,
+		   fcs_float* box_offset, fcs_float box_size,
+		   fcs_int near_field_cells, fcs_int interpolation_degree,
+                   fcs_int discretization_order, MPI_Comm mpi_comm);
+
+int VMG_fcs_check();
+
+void VMG_fcs_run(fcs_float* pos, fcs_float* charge, fcs_float* potential, fcs_float* f, fcs_int num_particles_local);
+
+void VMG_fcs_print_timer(void);
+
+void VMG_fcs_destroy(void);
+
+#ifdef __cplusplus
+} /* extern "C" */
+#endif
Index: ThirdParty/vmg/src/units/particle/interface_particles.cpp
===================================================================
--- ThirdParty/vmg/src/units/particle/interface_particles.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/units/particle/interface_particles.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,261 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   interface_particles.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Apr 18 12:56:48 2011
+ *
+ * @brief  VMG::InterfaceParticles
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#ifdef HAVE_MPI
+#include <mpi.h>
+#ifdef HAVE_MARMOT
+#include <enhancempicalls.h>
+#include <sourceinfompicalls.h>
+#endif
+#endif
+
+#include <algorithm>
+#include <cmath>
+#include <cstring>
+
+#include "base/helper.hpp"
+#include "base/index.hpp"
+#include "base/math.hpp"
+#include "base/vector.hpp"
+#include "comm/comm.hpp"
+#include "grid/grid.hpp"
+#include "grid/multigrid.hpp"
+#include "grid/tempgrid.hpp"
+#include "units/particle/comm_mpi_particle.hpp"
+#include "units/particle/interface_particles.hpp"
+#include "units/particle/interpolation.hpp"
+#include "units/particle/linked_cell_list.hpp"
+#include "mg.hpp"
+
+using namespace VMG;
+
+void InterfaceParticles::ImportRightHandSide(Multigrid& multigrid)
+{
+  Index index_global, index_local, index;
+  Vector pos_rel, pos_abs, grid_val;
+
+  Factory& factory = MG::GetFactory();
+  Particle::CommMPI& comm = *dynamic_cast<Particle::CommMPI*>(MG::GetComm());
+
+  const int& near_field_cells = factory.GetObjectStorageVal<int>("PARTICLE_NEAR_FIELD_CELLS");
+
+  Grid& grid = multigrid(multigrid.MaxLevel());
+  Grid& particle_grid = comm.GetParticleGrid();
+
+//  grid.Clear();
+  particle_grid.Clear();
+
+  assert(particle_grid.Global().LocalSize().IsComponentwiseGreater(near_field_cells));
+
+  /*
+   * Distribute particles to their processes
+   */
+  particles.clear();
+  comm.CommParticles(grid, particles);
+
+  /*
+   * Charge assignment on the grid
+   */
+  std::list<Particle::Particle>::iterator iter;
+
+#ifdef OUTPUT_DEBUG
+  vmg_float particle_charges = 0.0;
+  for (iter=particles.begin(); iter!=particles.end(); ++iter)
+    particle_charges += iter->Charge();
+  particle_charges = MG::GetComm()->GlobalSumRoot(particle_charges);
+  comm.PrintOnce(Debug, "Particle list charge sum: %e", particle_charges);
+  comm.Print(Debug, "Local number of particles: %d", particles.size());
+#endif
+
+  for (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 * Math::pi *
+	  particle_grid.GetVal(particle_grid.Local().Begin().X() + i,
+			       particle_grid.Local().Begin().Y() + j,
+			       particle_grid.Local().Begin().Z() + k);
+
+#ifdef OUTPUT_DEBUG
+  Grid::iterator grid_iter;
+  vmg_float charge_sum = 0.0;
+  for (grid_iter=grid.Iterators().Local().Begin(); grid_iter!=grid.Iterators().Local().End(); ++grid_iter)
+    charge_sum += grid.GetVal(*grid_iter);
+  charge_sum = MG::GetComm()->GlobalSum(charge_sum);
+  comm.PrintOnce(Debug, "Grid charge sum: %e", charge_sum);
+#endif
+}
+
+void InterfaceParticles::ExportSolution(Grid& grid)
+{
+  Index i;
+
+#ifdef OUTPUT_DEBUG
+  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;
+#endif
+
+  Factory& factory = MG::GetFactory();
+  Particle::CommMPI& comm = *dynamic_cast<Particle::CommMPI*>(MG::GetComm());
+
+  /*
+   * 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();
+//      spl.SubtractSelfInducedForces(particle_grid, **p1);
+
+#ifdef OUTPUT_DEBUG
+	  e_long += 0.5 * (*p1)->Charge() * ip.EvaluatePotentialLR(**p1);
+	  e_self += 0.5 * (*p1)->Charge() * (*p1)->Charge() * spl.GetAntiDerivativeAtZero();
+#endif
+
+	  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);
+
+#ifdef OUTPUT_DEBUG
+		      e_short_peak += 0.5 * (*p1)->Charge() * (*p2)->Charge() / length;
+		      e_short_spline += 0.5 * (*p1)->Charge() * (*p2)->Charge() / length * spl.EvaluatePotential(length);
+#endif
+		    }
+		  }
+	      }
+	}
+      }
+
+  /* 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);
+
+#ifdef OUTPUT_DEBUG
+  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");
+
+  // extract forces
+  if (num_particles_local != 0) {
+    size_t index = 0;
+    comm.PrintOnce(Debug, "%d force vector:        %e %e %e", (index/3)+1, f[index++], f[index++], f[index++]);
+  }
+
+  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.PrintOnce(Debug, "E_long:         %e", e_long);
+  comm.PrintOnce(Debug, "E_short_peak:   %e", e_short_peak);
+  comm.PrintOnce(Debug, "E_short_spline: %e", e_short_spline);
+  comm.PrintOnce(Debug, "E_self:         %e", e_self);
+  comm.PrintOnce(Debug, "E_total:        %e", e);
+  comm.PrintOnce(Debug, "E_total*:       %e", e_long + e_short_peak + e_short_spline - e_self);
+#endif
+}
Index: ThirdParty/vmg/src/units/particle/interface_particles.hpp
===================================================================
--- ThirdParty/vmg/src/units/particle/interface_particles.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/units/particle/interface_particles.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,68 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   interface_particles.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Tue Apr 12 17:40:36 2011
+ *
+ * @brief  Interface for computing forces and energies in
+ *         particle systems.
+ *
+ */
+
+#ifndef INTERFACE_PARTICLES_HPP
+#define INTERFACE_PARTIVLES_HPP
+
+#include <list>
+
+#include "base/defs.hpp"
+#include "base/interface.hpp"
+#include "units/particle/bspline.hpp"
+#include "units/particle/particle.hpp"
+
+namespace VMG
+{
+
+class Grid;
+class Multigrid;
+
+class InterfaceParticles : public Interface
+{
+public:
+  InterfaceParticles(const Boundary& boundary, const int& levelMin, const int& levelMax,
+		     const Vector& box_offset, const vmg_float& box_size,
+		     const int& near_field_cells,
+		     const int& max_boundary_nodes = 9, const vmg_float& alpha = 1.6) :
+    Interface(boundary, levelMin, levelMax, box_offset, box_size, max_boundary_nodes, alpha),
+    spl(near_field_cells, Extent(MaxLevel()).MeshWidth().Max())
+  {}
+
+  void ImportRightHandSide(Multigrid& multigrid);
+  void ExportSolution(Grid& grid);
+
+protected:
+  Particle::BSpline spl;
+
+private:
+  std::list<Particle::Particle> particles;
+};
+
+}
+
+#endif /* INTERFACE_PARTICLES_HPP */
Index: ThirdParty/vmg/src/units/particle/interpolation.cpp
===================================================================
--- ThirdParty/vmg/src/units/particle/interpolation.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/units/particle/interpolation.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,174 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include <vector>
+
+#include "base/helper.hpp"
+#include "base/index.hpp"
+#include "base/vector.hpp"
+#include "grid/grid.hpp"
+#include "units/particle/interpolation.hpp"
+#include "units/particle/particle.hpp"
+
+#include "mg.hpp"
+#include "comm/comm.hpp"
+
+using namespace VMG;
+
+Particle::Interpolation::Interpolation(const int& degree) :
+  deg(degree),
+  deg_1(degree+1),
+  buffer(degree+1),
+  buffer_diff(degree)
+{
+  coeff = new vmg_float[Helper::intpow(deg_1, 3)];
+  coeff_buffer = new vmg_float[deg_1];
+  for (int i=0; i<degree; ++i)
+    buffer_diff[i].resize(i+1);
+}
+
+Particle::Interpolation::~Interpolation()
+{
+  delete [] coeff;
+  delete [] coeff_buffer;
+}
+
+void Particle::Interpolation::ComputeCoefficients(const Grid& grid, const Index& index)
+{
+  Index i;
+
+  const Index begin = index - deg/2;
+
+  h = grid.Extent().MeshWidth();
+
+  for (i[0]=0; i[0]<deg_1; ++i[0])
+    for (i[1]=0; i[1]<deg_1; ++i[1])
+      for (i[2]=0; i[2]<deg_1; ++i[2])
+        _access_coeff(i) = grid.GetVal(i+begin);
+
+  pos_begin = grid.Extent().Begin()
+    + (begin - grid.Local().Begin() + grid.Global().LocalBegin() - grid.Global().GlobalBegin()) * grid.Extent().MeshWidth();
+
+  // compute coefficients x-direction
+  for (i=0; i[1]<deg_1; ++i[1])
+    for (i[2]=0; i[2]<deg_1; ++i[2])
+      _compute_coefficients_1d(i, 0);
+
+  // compute coefficients y-direction
+  for (i=0; i[0]<deg_1; ++i[0])
+    for (i[2]=0; i[2]<deg_1; ++i[2])
+      _compute_coefficients_1d(i, 1);
+
+  // compute coefficients z-direction
+  for (i=0; i[0]<deg_1; ++i[0])
+    for (i[1]=0; i[1]<deg_1; ++i[1])
+      _compute_coefficients_1d(i, 2);
+}
+
+void Particle::Interpolation::_compute_coefficients_1d(const Index& index, const unsigned int& direction)
+{
+  vmg_float power = 1.0;
+  unsigned long faculty = 1;
+  int c;
+  Index i;
+
+  for (i=index, c=0; c<deg_1; ++i[direction], ++c)
+    coeff_buffer[c] = _access_coeff(i);
+
+  i=index;
+  ++i[direction];
+  for (c=1; c<deg_1; ++i[direction], ++c) {
+    for (int j=0; j<deg_1-c; ++j)
+      coeff_buffer[j] = coeff_buffer[j+1] - coeff_buffer[j];
+    faculty *= c;
+    power *= h[direction];
+    _access_coeff(i) = coeff_buffer[0] / (faculty*power);
+  }
+}
+
+void Particle::Interpolation::Evaluate(Particle& p)
+{
+  const Vector& pos = p.Pos();
+  vmg_float& pot = p.Pot();
+  Vector& field = p.Field();
+
+  pot = 0.0;
+  field = 0.0;
+
+  Vector offset = pos - pos_begin;
+  buffer[0] = 1.0;
+  for (int i=0; i<deg; ++i) {
+    buffer[i+1] = buffer[i] * offset;
+    for (int j=0; j<i; ++j)
+      buffer_diff[i][j] = buffer_diff[i-1][j] * offset;
+    buffer_diff[i][i] = buffer[i];
+    offset -= h;
+  }
+
+  for (int i=1; i<deg; ++i)
+    for (int j=1; j<=i; ++j)
+      buffer_diff[i][0] += buffer_diff[i][j];
+
+  for (int i=0; i<deg_1; ++i)
+    for (int j=0; j<deg_1; ++j)
+      for (int k=0; k<deg_1; ++k)
+ 	pot += _access_coeff(i,j,k) * buffer[i][0] * buffer[j][1] * buffer[k][2];
+
+  for (int i=0; i<deg_1; ++i)
+    for (int j=0; j<deg_1; ++j)
+      for (int k=0; k<deg; ++k) {
+	field[0] -= _access_coeff(k+1, i, j) * buffer[i][1] * buffer[j][2] * buffer_diff[k][0][0];
+	field[1] -= _access_coeff(i, k+1, j) * buffer[i][0] * buffer[j][2] * buffer_diff[k][0][1];
+	field[2] -= _access_coeff(i, j, k+1) * buffer[i][0] * buffer[j][1] * buffer_diff[k][0][2];
+      }
+}
+
+vmg_float Particle::Interpolation::EvaluatePotentialLR(const Particle& p)
+{
+  vmg_float result = 0.0;
+  Vector prod, offset;
+  Index i;
+
+  const Vector& pos = p.Pos();
+
+  prod[0] = 1.0;
+  offset[0] = pos[0] - pos_begin[0];
+  for (i[0]=0; i[0]<deg_1; ++i[0]) {
+    prod[1] = prod[0];
+    offset[1] = pos[1] - pos_begin[1];
+    for (i[1]=0; i[1]<deg_1; ++i[1]) {
+      prod[2] = prod[1];
+      offset[2] = pos[2] - pos_begin[2];
+      for (i[2]=0; i[2]<deg_1; ++i[2]) {
+        result += _access_coeff(i) * prod[2];
+        prod[2] *= offset[2];
+        offset[2] -= h[2];
+      }
+      prod[1] *= offset[1];
+      offset[1] -= h[1];
+    }
+    prod[0] *= offset[0];
+    offset[0] -= h[0];
+  }
+
+  return result;
+}
Index: ThirdParty/vmg/src/units/particle/interpolation.hpp
===================================================================
--- ThirdParty/vmg/src/units/particle/interpolation.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/units/particle/interpolation.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,74 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef INTERPOLATION_HPP_
+#define INTERPOLATION_HPP_
+
+#include <vector>
+
+#include "base/index.hpp"
+#include "base/vector.hpp"
+
+namespace VMG
+{
+
+class Grid;
+
+namespace Particle
+{
+
+class Particle;
+
+class Interpolation
+{
+public:
+  Interpolation(const int& degree);
+  ~Interpolation();
+
+  void ComputeCoefficients(const Grid& grid, const Index& index);
+  void Evaluate(Particle& p);
+
+  vmg_float EvaluatePotentialLR(const Particle& p);
+
+private:
+  vmg_float& _access_coeff(const Index& index)
+  {
+    return coeff[index.Z() + deg_1 * (index.Y() + deg_1 * index.X())];
+  }
+
+  vmg_float& _access_coeff(const int& i, const int& j, const int& k)
+  {
+    return coeff[k + deg_1 * (j + deg_1 * i)];
+  }
+
+  void _compute_coefficients_1d(const Index& index, const unsigned int& direction);
+
+  vmg_float* coeff;
+  vmg_float* coeff_buffer;
+  int deg, deg_1;
+  Vector pos_begin;
+  Vector h;
+  std::vector<Vector> buffer;
+  std::vector< std::vector<Vector> > buffer_diff;
+};
+
+}
+
+}
+
+#endif /* INTERPOLATION_HPP_ */
Index: ThirdParty/vmg/src/units/particle/linked_cell_list.cpp
===================================================================
--- ThirdParty/vmg/src/units/particle/linked_cell_list.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/units/particle/linked_cell_list.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,131 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   linked_cell_list.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Nov 21 13:27:22 2011
+ *
+ * @brief  A linked cell list implementation.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include "comm/comm.hpp"
+#include "units/particle/linked_cell_list.hpp"
+#include "mg.hpp"
+
+using namespace VMG;
+
+Particle::LinkedCellList::LinkedCellList(std::list<Particle>& particles,
+					 const int& near_field_cells, const Grid& grid)
+{
+  std::list<Particle>::iterator iter;
+  Index index_global, index_local;
+  LocalIndices local = grid.Local();
+
+  local.BoundaryBegin1() = 0;
+  local.BoundaryEnd1() = 0;
+  local.BoundarySize1() = 0;
+
+  local.BoundaryBegin2() = 0;
+  local.BoundaryEnd2() = 0;
+  local.BoundarySize2() = 0;
+
+  for (int i=0; i<3; ++i) {
+
+    if (local.HaloSize1()[i] > 0) {
+      local.HaloBegin1()[i] = 0;
+      local.HaloEnd1()[i] = near_field_cells+1;
+      local.HaloSize1()[i] = near_field_cells+1;
+      local.Begin()[i] = near_field_cells+1;
+      local.End()[i] = local.Begin()[i] + local.Size()[i];
+      local.SizeTotal()[i] = local.Size()[i] +
+	                     local.HaloEnd1()[i] - local.HaloBegin1()[i] +
+	                     local.HaloEnd2()[i] - local.HaloBegin2()[i];
+    }else {
+      local.Begin()[i] = 0;
+      local.End()[i] = local.Size()[i];
+    }
+
+    if (local.HaloSize2()[i] > 0) {
+      local.HaloBegin2()[i] = local.End()[i];
+      local.HaloEnd2()[i] = local.HaloBegin2()[i] + near_field_cells+1;
+      local.HaloSize2()[i] = near_field_cells+1;
+      local.SizeTotal()[i] = local.Size()[i] +
+	                     local.HaloEnd1()[i] - local.HaloBegin1()[i] +
+	                     local.HaloEnd2()[i] - local.HaloBegin2()[i];
+    }
+  }
+
+  SetGridSize(grid.Global(), local, grid.Extent());
+
+  for (iter=particles.begin(); iter!=particles.end(); ++iter)
+    AddParticle(&(*iter));
+}
+
+Particle::LinkedCellList::~LinkedCellList()
+{
+  ClearHalo();
+}
+
+void Particle::LinkedCellList::AddParticle(Particle* p)
+{
+  const Index global_index = global.GlobalBegin() + ((p->Pos() - Extent().Begin()) / Extent().MeshWidth()).Floor();
+  const Index local_index = global_index - Global().LocalBegin() + Local().Begin();
+
+  assert(local_index.IsInBounds(Local().Begin(), Local().End()));
+  (*this)(local_index).push_back(p);
+}
+
+void Particle::LinkedCellList::AddParticleToHalo(const vmg_float* x, const vmg_float& q)
+{
+  const Index global_index = Global().GlobalBegin() + ((Vector(x) - Extent().Begin()) / Extent().MeshWidth()).Floor();
+  const Index local_index = global_index - Global().LocalBegin() + Local().Begin();
+
+  assert(local_index.IsInBounds(0, Local().SizeTotal()));
+  assert((local_index[0] >= Local().HaloBegin1()[0] && local_index[0] < Local().HaloEnd1()[0]) ||
+	 (local_index[1] >= Local().HaloBegin1()[1] && local_index[1] < Local().HaloEnd1()[1]) ||
+	 (local_index[2] >= Local().HaloBegin1()[2] && local_index[2] < Local().HaloEnd1()[2]) ||
+	 (local_index[0] >= Local().HaloBegin2()[0] && local_index[0] < Local().HaloEnd2()[0]) ||
+	 (local_index[1] >= Local().HaloBegin2()[1] && local_index[1] < Local().HaloEnd2()[1]) ||
+	 (local_index[2] >= Local().HaloBegin2()[2] && local_index[2] < Local().HaloEnd2()[2]));
+
+  (*this)(local_index).push_back(new Particle(x, q));
+}
+
+void Particle::LinkedCellList::ClearHalo()
+{
+  Grid::iterator g_iter;
+  std::list<Particle*>::iterator p_iter;
+
+  for (int i=0; i<3; ++i) {
+
+    for (g_iter = Iterators().Halo1()[i].Begin(); g_iter != Iterators().Halo1()[i].End(); ++g_iter)
+      for (p_iter = (*this)(*g_iter).begin(); p_iter!=(*this)(*g_iter).end(); ++p_iter)
+	delete *p_iter;
+
+    for (g_iter = Iterators().Halo2()[i].Begin(); g_iter != Iterators().Halo2()[i].End(); ++g_iter)
+      for (p_iter = (*this)(*g_iter).begin(); p_iter!=(*this)(*g_iter).end(); ++p_iter)
+	delete *p_iter;
+
+  }
+}
Index: ThirdParty/vmg/src/units/particle/linked_cell_list.hpp
===================================================================
--- ThirdParty/vmg/src/units/particle/linked_cell_list.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/units/particle/linked_cell_list.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,71 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   linked_cell_list.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Nov 21 13:27:22 2011
+ *
+ * @brief  A linked cell list implementation.
+ *
+ */
+
+#ifndef LINKED_CELL_LIST_HPP_
+#define LINKED_CELL_LIST_HPP_
+
+#include <list>
+#include <vector>
+
+#include "base/index.hpp"
+#include "grid/is_grid.hpp"
+#include "units/particle/particle.hpp"
+
+namespace VMG
+{
+
+namespace Particle
+{
+
+class LinkedCellList : public IsGrid< std::list<Particle*> >
+{
+public:
+  typedef std::list<Particle*>::iterator iterator;
+  typedef std::list<Particle*>::const_iterator const_iterator;
+  typedef std::list<Particle*>::reverse_iterator reverse_iterator;
+  typedef std::list<Particle*>::const_reverse_iterator const_reverse_iterator;
+
+  LinkedCellList(std::list<Particle>& particles,
+		 const int& near_field_cells, const Grid& grid);
+
+  ~LinkedCellList();
+
+  void AddParticle(Particle* particle);
+  void AddParticleToHalo(const vmg_float* x, const vmg_float& q);
+
+  void ClearHalo();
+
+  const Index& NearFieldCells() const {return near_field_cells_;}
+private:
+  Index near_field_cells_;
+};
+
+}
+
+}
+
+#endif /* LINKED_CELL_LIST_HPP_ */
Index: ThirdParty/vmg/src/units/particle/particle.cpp
===================================================================
--- ThirdParty/vmg/src/units/particle/particle.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/units/particle/particle.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,37 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   particle.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Sat Sep 17 14:45:56 2011
+ *
+ * @brief  Class to represent a particle
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include "libvmg_config.h"
+#endif
+
+#include "units/particle/particle.hpp"
+
+std::ostream& VMG::Particle::operator<<(std::ostream& out, const VMG::Particle::Particle& particle)
+{
+  return out << particle;
+}
Index: ThirdParty/vmg/src/units/particle/particle.hpp
===================================================================
--- ThirdParty/vmg/src/units/particle/particle.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/units/particle/particle.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,145 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   particle.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Sat Sep 17 14:45:56 2011
+ *
+ * @brief  Class to represent a particle
+ *
+ */
+
+#ifndef PARTICLE_HPP_
+#define PARTICLE_HPP_
+
+#include <sstream>
+#include <string>
+
+#include "base/vector.hpp"
+
+namespace VMG
+{
+
+namespace Particle
+{
+
+class Particle
+{
+public:
+  Particle() :
+    x_(0.0),
+    f_(0.0),
+    q_(0.0),
+    p_(0.0),
+    rank_(-1),
+    index_(-1)
+  {}
+
+  Particle(const vmg_float* x, const vmg_float& q) :
+    x_(x),
+    f_(0.0),
+    q_(q),
+    p_(0.0),
+    rank_(-1),
+    index_(-1)
+  {}
+
+  Particle(const vmg_float* x, const vmg_float& q, const vmg_float& p, const vmg_float* f, const int& rank, const vmg_int& index) :
+    x_(x),
+    f_(f),
+    q_(q),
+    p_(p),
+    rank_(rank),
+    index_(index)
+  {}
+
+  Particle(const Vector& x, const vmg_float& q, const vmg_float& p, const Vector& f, const int& rank, const vmg_int& index) :
+    x_(x),
+    f_(f),
+    q_(q),
+    p_(p),
+    rank_(rank),
+    index_(index)
+  {}
+
+  Particle(const Particle& other) :
+    x_(other.x_),
+    f_(other.f_),
+    q_(other.q_),
+    p_(other.p_),
+    rank_(other.rank_),
+    index_(other.index_)
+  {}
+
+  Vector& Pos() {return x_;}
+  const Vector& Pos() const {return x_;}
+
+  vmg_float& Charge() {return q_;}
+  const vmg_float& Charge() const {return q_;}
+
+  vmg_float& Pot() {return p_;}
+  const vmg_float& Pot() const {return p_;}
+
+  Vector& Field() {return f_;}
+  const Vector& Field() const {return f_;}
+
+  int& Rank() {return rank_;}
+  const int& Rank() const {return rank_;}
+
+  vmg_int& Index() {return index_;}
+  const vmg_int& Index() const {return index_;}
+
+  bool operator==(const Particle& rhs)
+  {
+    return (this->rank_ == rhs.rank_) && (this->index_ == rhs.index_);
+  }
+
+  bool operator!=(const Particle& rhs)
+  {
+    return (this->rank_ != rhs.rank_) || (this->index_ != rhs.index_);
+  }
+
+  std::string ToString() const
+  {
+    std::stringstream str;
+
+    str << "VMG::Particle::Particle {" << std::endl
+        << "  Pos:       " << x_ << std::endl
+        << "  Charge:    " << q_ << std::endl
+        << "  Potential: " << p_ << std::endl
+        << "  Field:     " << f_ << std::endl
+        << "}";
+
+    return str.str();
+  }
+
+private:
+  Vector x_, f_;
+  vmg_float q_, p_;
+  int rank_;
+  vmg_int index_;
+};
+
+std::ostream& operator<<(std::ostream& out, const VMG::Particle::Particle& particle);
+
+}
+
+}
+
+#endif /* PARTICLE_HPP_ */
