diff --git a/bfps/cpp/interpolator.cpp b/bfps/cpp/interpolator.cpp
index 28312c3141350764facfbbe86074092496c273ae..cf04af668c0cee55b1e4181ebc2654d1d4612965 100644
--- a/bfps/cpp/interpolator.cpp
+++ b/bfps/cpp/interpolator.cpp
@@ -163,9 +163,11 @@ void interpolator<rnumber, interp_neighbours>::operator()(
     }
     std::fill_n(dest, 3, 0);
     ptrdiff_t bigiz, bigiy, bigix;
+    double tval[3];
     for (int iz = -interp_neighbours; iz <= interp_neighbours+1; iz++)
     {
         bigiz = ptrdiff_t(xg[2]+iz);
+        std::fill_n(tval, 3, 0);
         for (int iy = -interp_neighbours; iy <= interp_neighbours+1; iy++)
         {
             bigiy = ptrdiff_t(MOD(xg[1]+iy, this->descriptor->sizes[1]));
@@ -180,9 +182,13 @@ void interpolator<rnumber, interp_neighbours>::operator()(
                     dest[c] += (this->f0[tindex]*(1-t) + t*this->f1[tindex])*(bz[iz+interp_neighbours]*
                                                                               by[iy+interp_neighbours]*
                                                                               bx[ix+interp_neighbours]);
+                    tval[c] += (this->f0[tindex]*(1-t) + t*this->f1[tindex])*(bz[iz+interp_neighbours]*
+                                                                              by[iy+interp_neighbours]*
+                                                                              bx[ix+interp_neighbours]);
                 }
             }
         }
+        DEBUG_MSG("%ld %d %d %g %g %g\n", bigiz, xg[1], xg[0], tval[0], tval[1], tval[2]);
     }
 }
 
diff --git a/bfps/cpp/rFFTW_interpolator.cpp b/bfps/cpp/rFFTW_interpolator.cpp
index 5e5f00eea504ded77d6dd0718443f9cfe80b5455..e4bd10a831311e6d38cd75aaa7c71db7d9af2dd9 100644
--- a/bfps/cpp/rFFTW_interpolator.cpp
+++ b/bfps/cpp/rFFTW_interpolator.cpp
@@ -103,10 +103,13 @@ void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
     }
     std::fill_n(dest, 3, 0);
     ptrdiff_t bigiz, bigiy, bigix;
+    double tval[3];
     for (int iz = -interp_neighbours; iz <= interp_neighbours+1; iz++)
     {
         bigiz = ptrdiff_t(((xg[2]+iz) + this->descriptor->sizes[0]) % this->descriptor->sizes[0]);
         if (this->descriptor->myrank == this->descriptor->rank[bigiz])
+        {
+            std::fill_n(tval, 3, 0);
             for (int iy = -interp_neighbours; iy <= interp_neighbours+1; iy++)
             {
                 bigiy = ptrdiff_t(MOD(xg[1]+iy, this->descriptor->sizes[1]));
@@ -116,14 +119,19 @@ void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
                     for (int c=0; c<3; c++)
                     {
                         ptrdiff_t tindex = ((bigiz *this->descriptor->sizes[1] +
-                                             bigiy)*this->descriptor->sizes[2] +
+                                             bigiy)*(this->descriptor->sizes[2]+2) +
                                              bigix)*3+c;
                         dest[c] += (this->f0[tindex]*(1-t) + t*this->f1[tindex])*(bz[iz+interp_neighbours]*
                                                                                   by[iy+interp_neighbours]*
                                                                                   bx[ix+interp_neighbours]);
+                        tval[c] += (this->f0[tindex]*(1-t) + t*this->f1[tindex])*(bz[iz+interp_neighbours]*
+                                                                                  by[iy+interp_neighbours]*
+                                                                                  bx[ix+interp_neighbours]);
                     }
                 }
             }
+            DEBUG_MSG("%ld %d %d %g %g %g\n", bigiz, xg[1], xg[0], tval[0], tval[1], tval[2]);
+        }
     }
 }
 
diff --git a/bfps/cpp/rFFTW_particles.cpp b/bfps/cpp/rFFTW_particles.cpp
index 401aa72c2bba17431ec40edf4ab6d4d02db77d5f..af225f9541d62adf93f1ff89aeadfb9a7cbe0cf1 100644
--- a/bfps/cpp/rFFTW_particles.cpp
+++ b/bfps/cpp/rFFTW_particles.cpp
@@ -141,13 +141,7 @@ void rFFTW_particles<particle_type, rnumber, interp_neighbours>::sample_vec_fiel
     this->get_grid_coordinates(x, xg, xx);
     /* perform interpolation */
     for (int p=0; p<this->nparticles; p++)
-    {
-        if (this->myrank == this->get_rank(x[3*p+2]))
-        {
-            DEBUG_MSG("particle %d\n", p);
-            (*vec)(t, xg + p*3, xx + p*3, yy + p*3, deriv);
-        }
-    }
+        (*vec)(t, xg + p*3, xx + p*3, yy + p*3, deriv);
     MPI_Allreduce(
             yy,
             y,
@@ -292,10 +286,10 @@ void rFFTW_particles<particle_type, rnumber, interp_neighbours>::get_grid_coordi
             xg[p*3+c] = MOD(int(tval), this->fs->rd->sizes[2-c]);
             xx[p*3+c] = (x[p*this->ncomponents+c] - tval*grid_size[c]) / grid_size[c];
         }
-        xg[p*3+2] -= this->fs->rd->starts[0];
+        /*xg[p*3+2] -= this->fs->rd->starts[0];
         if (this->myrank == this->fs->rd->rank[0] &&
             xg[p*3+2] > this->fs->rd->subsizes[0])
-            xg[p*3+2] -= this->fs->rd->sizes[0];
+            xg[p*3+2] -= this->fs->rd->sizes[0];*/
     }
 }