diff --git a/bfps/cpp/rFFTW_interpolator.cpp b/bfps/cpp/rFFTW_interpolator.cpp
index 5e95d742c414c6ee694cad3839a8b4aa4125bd97..46673eda1e4dff0b63e2bd91488429078ba5d633 100644
--- a/bfps/cpp/rFFTW_interpolator.cpp
+++ b/bfps/cpp/rFFTW_interpolator.cpp
@@ -121,12 +121,35 @@ int rFFTW_interpolator<rnumber, interp_neighbours>::read_rFFTW(void *void_src)
 }
 
 template <class rnumber, int interp_neighbours>
-void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
-        double t,
+void rFFTW_interpolator<rnumber, interp_neighbours>::get_grid_coordinates(
+        const int nparticles,
+        const int pdimension,
+        const double *x,
         int *xg,
-        double *xx,
+        double *xx)
+{
+    static double grid_size[] = {this->dx, this->dy, this->dz};
+    double tval;
+    std::fill_n(xg, nparticles*3, 0);
+    std::fill_n(xx, nparticles*3, 0.0);
+    for (int p=0; p<nparticles; p++)
+    {
+        for (int c=0; c<3; c++)
+        {
+            tval = floor(x[p*pdimension+c]/grid_size[c]);
+            xg[p*3+c] = MOD(int(tval), this->descriptor->sizes[2-c]);
+            xx[p*3+c] = (x[p*pdimension+c] - tval*grid_size[c]) / grid_size[c];
+        }
+    }
+}
+
+template <class rnumber, int interp_neighbours>
+void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
+        const double t,
+        const int *xg,
+        const double *xx,
         double *dest,
-        int *deriv)
+        const int *deriv)
 {
     double bx[interp_neighbours*2+2], by[interp_neighbours*2+2], bz[interp_neighbours*2+2];
     if (deriv == NULL)
diff --git a/bfps/cpp/rFFTW_interpolator.hpp b/bfps/cpp/rFFTW_interpolator.hpp
index 0e9be89607fcb3546e7189519772cf0755a194b3..c7efb9d8f60f616c131509a260c76ba3d61507a4 100644
--- a/bfps/cpp/rFFTW_interpolator.hpp
+++ b/bfps/cpp/rFFTW_interpolator.hpp
@@ -71,22 +71,24 @@ class rFFTW_interpolator
 
         /* map real locations to grid coordinates */
         void get_grid_coordinates(
+                const int nparticles,
+                const int pdimension,
                 const double *__restrict__ x,
-                const int *__restrict__ xg,
-                const double *__restrict__ xx);
+                int *__restrict__ xg,
+                double *__restrict__ xx);
         /* interpolate field at an array of locations */
         void sample(
                 const int nparticles,
                 const double t,
                 const double *__restrict__ x,
-                const double *__restrict__ y,
+                double *__restrict__ y,
                 const int *deriv = NULL);
         /* interpolate 1 point */
         void operator()(
                 const double t,
                 const int *__restrict__ xg,
                 const double *__restrict__ xx,
-                const double *__restrict__ dest,
+                double *__restrict__ dest,
                 const int *deriv = NULL);
         int read_rFFTW(void *src);
 };
diff --git a/bfps/cpp/rFFTW_particles.cpp b/bfps/cpp/rFFTW_particles.cpp
index 42a308f2d791a83df5208965e18a57898fbe406a..4d0ab85b1a8cc0c16ac2424bdc7eafb823867a2a 100644
--- a/bfps/cpp/rFFTW_particles.cpp
+++ b/bfps/cpp/rFFTW_particles.cpp
@@ -132,7 +132,7 @@ void rFFTW_particles<particle_type, rnumber, interp_neighbours>::sample_vec_fiel
     double *xx = new double[3*this->nparticles];
     double *yy =  new double[3*this->nparticles];
     std::fill_n(yy, 3*this->nparticles, 0.0);
-    this->get_grid_coordinates(x, xg, xx);
+    vec->get_grid_coordinates(this->nparticles, this->ncomponents, x, xg, xx);
     /* perform interpolation */
     for (int p=0; p<this->nparticles; p++)
         (*vec)(t, xg + p*3, xx + p*3, yy + p*3, deriv);
diff --git a/done.txt b/done.txt
index 787392010cce3f4cb0f5921eba4a7fa8ba31e89a..8dd974c1f1353734f7a4f17dc79c6ca548c88cba 100644
--- a/done.txt
+++ b/done.txt
@@ -1,2 +1,3 @@
 x 2015-12-04 make code py3 compatible                                @python3
 x 2015-12-23 decide on versioning system                                           +merge0
+x 2015-12-24 move get_grid coords to interpolator                    @optimization +v1.0
diff --git a/todo.txt b/todo.txt
index 71c2677e5fa6af740a768f8fbee676c3357e1a75..dcc583669300cd001fe198ec6ac15659bcdcdfa3 100644
--- a/todo.txt
+++ b/todo.txt
@@ -1,7 +1,6 @@
 (B) FFTW interpolator doesn't need its own field            @optimization +v1.0
 (B) compute z polynomials only when needed                  @optimization +v1.0
 (B) get rid of temporal interpolation                       @optimization +v1.0
-(B) move get_grid coords to interpolator                    @optimization +v1.0
 (B) use less memory                                         @optimization
 (C) clean up machine_settings mess                          @design @documentation +v2.0
 (C) code overview                                           @documentation