Index: src/Jobs/ChargeSmearer.cpp
===================================================================
--- src/Jobs/ChargeSmearer.cpp	(revision 17e4fde191dda907bf7bd206803886c13e02d440)
+++ src/Jobs/ChargeSmearer.cpp	(revision 62d092e08dbdb67e4f00547587c38322bf8bb828)
@@ -37,4 +37,5 @@
 #include "ChargeSmearer.hpp"
 
+#include "CodePatterns/Log.hpp"
 #include "CodePatterns/Singleton_impl.hpp"
 
@@ -68,13 +69,20 @@
   // reallocate memory for spline values
   delete[] vals;
-  double *vals = new double[VMG::Helper::intpow(2*nfc+1,3)];
+  vals = new double[VMG::Helper::intpow(2*nfc+1,3)];
 
   // recalculate spline values
   unsigned int c=0;
   int_val = 0.;
-  for (double dir_x=-1*nfc*meshwidth; dir_x<=nfc*meshwidth; dir_x+=meshwidth) {
-    for (double dir_y=-1*nfc*meshwidth; dir_y<=nfc*meshwidth; dir_y+=meshwidth) {
-      for (double dir_z=-1*nfc*meshwidth; dir_z<=nfc*meshwidth; dir_z+=meshwidth) {
-
+  std::pair<int, int> bounds[3];
+  for (int dim=0;dim<3;++dim) {
+    bounds[dim].first = -1*nfc;
+    bounds[dim].second = nfc;
+  }
+  for (int i=bounds[0].first; i<=bounds[0].second; ++i)
+    for (int j=bounds[1].first; j<=bounds[1].second; ++j)
+      for (int k=bounds[2].first; k<=bounds[2].second; ++k) {
+        const double dir_x = (double)i*meshwidth;
+        const double dir_y = (double)j*meshwidth;
+        const double dir_z = (double)k*meshwidth;
         // Compute distance from grid point to particle
         const double distance = std::sqrt(dir_x*dir_x+dir_y*dir_y+dir_z*dir_z);
@@ -82,10 +90,47 @@
         vals[c++] = temp_val;
         int_val += temp_val;
+        LOG(4, "DEBUG: Spline evaluated at " << distance << " is " << temp_val);
       }
-    }
-  }
 
   /// then integrate
-  int_val = 1.0 / (int_val * meshwidth * meshwidth * meshwidth);
+  int_val = 1.0 / int_val;
+
+  LOG(3, "DEBUG: Spline for smearing electronic charge distribution has norm factor of " << int_val);
+}
+
+void ChargeSmearer::visitBSplineDomain(
+    const VMG::GridIterator &_iter,
+    const visitor_t &_visitor
+    ) const
+{
+  unsigned int c = 0;
+  VMG::Index index;
+  const VMG::Index end = _iter.GetIndex() + (int)nfc;
+  for (index[0] = _iter.GetIndex()[0] - (int)nfc; index[0]<=end[0]; ++index[0])
+    for (index[1] = _iter.GetIndex()[1] - (int)nfc; index[1]<=end[1]; ++index[1])
+      for (index[2] = _iter.GetIndex()[2] - (int)nfc; index[2]<=end[2]; ++index[2]) {
+        if (index.IsInBounds(_iter.GetBegin(), _iter.GetEnd()))
+          _visitor(index, vals[c++]);
+        else
+          ++c;
+      }
+  assert( c == (unsigned int)VMG::Helper::intpow(2*nfc+1,3) );
+}
+
+static void DomainIntegrator(
+    const VMG::Index & _index,
+    const double _val,
+    double &_int_val)
+{
+  _int_val += _val;
+}
+
+static void SplineSetter(
+    const VMG::Index & _index,
+    const double _val,
+    const double _factor,
+    VMG::Grid& _grid)
+{
+  _grid(_index) += _val * _factor;
 }
 
@@ -95,18 +140,26 @@
     const double _charge) const
 {
-  /// prepare bounds
-  std::pair<int, int> bounds[3];
+  // check whether we go over whole bspline support
+  bool renormalize = false;
   for (int dim=0;dim<3;++dim) {
-    bounds[dim].first = std::max(_iter.GetBegin()[dim], _iter.GetIndex()[dim] - (int)nfc);
-    bounds[dim].second = std::min(_iter.GetEnd()[dim], _iter.GetIndex()[dim] + (int)nfc);
+    renormalize |= (_iter.GetBegin()[dim] > _iter.GetIndex()[dim] - (int)nfc);
+    renormalize |= (_iter.GetEnd()[dim] < _iter.GetIndex()[dim] + (int)nfc);
+  }
+
+  /// renormalize bspline if necessary (bounds don't fit in window)
+  double temp_int_val = 0.;
+  if (renormalize) {
+    const visitor_t integrator = boost::bind(DomainIntegrator, _1, _2, boost::ref(temp_int_val));
+    visitBSplineDomain(_iter, integrator);
+    temp_int_val = 1.0 / temp_int_val;
+  } else {
+    temp_int_val = int_val;
   }
 
   /// then transfer to grid in bounded window
-  unsigned int c = 0;
-  for (int i=bounds[0].first; i<=bounds[0].second; ++i)
-    for (int j=bounds[1].first; j<=bounds[1].second; ++j)
-      for (int k=bounds[2].first; k<=bounds[2].second; ++k)
-        _grid(i,j,k) += vals[c++] * int_val * _charge;
-  assert( c == (unsigned int)VMG::Helper::intpow(2*nfc+1,3) );
+  {
+    const visitor_t setter = boost::bind(SplineSetter, _1, _2, temp_int_val * _charge, boost::ref(_grid));
+    visitBSplineDomain(_iter, setter);
+  }
 }
 
Index: src/Jobs/ChargeSmearer.hpp
===================================================================
--- src/Jobs/ChargeSmearer.hpp	(revision 17e4fde191dda907bf7bd206803886c13e02d440)
+++ src/Jobs/ChargeSmearer.hpp	(revision 62d092e08dbdb67e4f00547587c38322bf8bb828)
@@ -19,4 +19,5 @@
   class Grid;
   class GridIterator;
+  class Index;
 
   namespace Particle {
@@ -109,4 +110,19 @@
 
 private:
+  //!> typedef for the visitor functor
+  typedef boost::function<void (const VMG::Index &, const double)> visitor_t;
+
+  /** Visitor pattern for going over bspline support and executing the given
+   * functor at each valid node.
+   *
+   * \param _iter iterator giving central node for support on grid
+   * \param _visitor visitor functor to execute at each node (that's still inside grid)
+   */
+  void visitBSplineDomain(
+      const VMG::GridIterator &_iter,
+      const visitor_t &_visitor
+      ) const;
+
+private:
   //!> near field radius for bspline
   unsigned int nfc;
