diff --git a/phigrape.cpp b/phigrape.cpp
index 1d18693..2976c59 100644
--- a/phigrape.cpp
+++ b/phigrape.cpp
@@ -163,7 +163,7 @@ public:
     {
         ind_act_loc.resize(n_loc);
     }
-    double get_minimum_time(const std::vector<double> &t, std::vector<double> &dt)
+    double get_minimum_time(const std::vector<double> &t, const std::vector<double> &dt)
     {
         double min_t_loc, min_t;
 #ifdef ETICS
@@ -249,6 +249,29 @@ inline double aarseth_step(const double eta, const double dt, const double3 &a,
     return sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
 }
 
+inline double blockize_step(double dt, double dt_prev, double min_t, double dt_min, double dt_max)
+{
+    double dt_new = dt_prev;
+    if (dt < dt_min) dt_prev = dt_min;
+    if ((dt < dt_prev) && (dt > dt_min)) {
+        int power = log(dt)/M_LN2 - 1;
+        dt_new = pow(2.0, power);
+    }
+    if ((dt > 2*dt_new) && (fmod(min_t, 2*dt_new) == 0) && (2*dt_new <= dt_max)) dt_new *= 2;
+    return dt_new;
+}
+
+inline void predictor(double min_t, const int n_act, const std::vector<int> &ind_act, const std::vector<double> &t, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double3> &a, const std::vector<double3> &adot, std::vector<double3> &x_act_new, std::vector<double3> &v_act_new)
+{
+    for (int i=0; i<n_act; i++) {
+        int j_act = ind_act[i];
+        double dt = min_t - t[j_act];
+        double dt2half = 0.5*dt*dt;
+        double dt3over6 = (1./3.)*dt*dt2half;
+        x_act_new[i] = x[j_act] + v[j_act]*dt + a[j_act]*dt2half + adot[j_act]*dt3over6;
+        v_act_new[i] = v[j_act] + a[j_act]*dt + adot[j_act]*dt2half;
+    } /* i */
+}
 
 int main(int argc, char *argv[])
 {
@@ -521,6 +544,8 @@ int main(int argc, char *argv[])
         /* Find the BH(s) indices in the active list */
         int i_bh1=0, i_bh2=1;
 #ifdef ETICS
+        /* Unlike with the simple active search, with GPU accelerated GRAPite
+        active search, the list of active indices is not sorted. */
         int n_bh = config.live_smbh_count;
         if (config.grapite_active_search && (n_bh>0)) {
             int act_def_grapite_bh_count = 0;
@@ -537,14 +562,7 @@ int main(int argc, char *argv[])
 #endif
 
         /* predict the active particles positions etc... on all the nodes */
-        for (int i=0; i<n_act; i++) {
-            int j_act = ind_act[i];
-            double dt = min_t - t[j_act];
-            double dt2half = 0.5*dt*dt;
-            double dt3over6 = (1./3.)*dt*dt2half;
-            x_act_new[i] = x[j_act] + v[j_act]*dt + a[j_act]*dt2half + adot[j_act]*dt3over6;
-            v_act_new[i] = v[j_act] + a[j_act]*dt + adot[j_act]*dt2half;
-        } /* i */
+        predictor(min_t, n_act, ind_act, t, x, v, a, adot, x_act_new, v_act_new);
 
         /* Calculate gravity on active particles */
         calc_self_grav(min_t, n_act, ind_act, x_act_new, v_act_new, pot_act_new, a_act_new, adot_act_new);
@@ -562,46 +580,37 @@ int main(int argc, char *argv[])
         }
 
         /* correct the active particles positions etc... on all the nodes */
-        double min_dt = dt_max;
+        double min_dt = dt_max; // notice that min_dt is not the same as dt_min; this one is to store the minimum timestep among currently active particles
         for (int i=0; i<n_act; i++) {
             int j_act = ind_act[i];
-            double dt_tmp = min_t - t[j_act];
+            double dt_cur = dt[j_act];
 
             double3 a2, a3;
-            calc_high_derivatives(dt_tmp, a[j_act], a_act_new[i], adot[j_act], adot_act_new[i], a2, a3);
+            calc_high_derivatives(dt_cur, a[j_act], a_act_new[i], adot[j_act], adot_act_new[i], a2, a3);
 
-            corrector(dt_tmp, a2, a3, x_act_new[i], v_act_new[i]);
+            corrector(dt_cur, a2, a3, x_act_new[i], v_act_new[i]);
 
             //TODO make beautiful
             double eta_curr;
             if ((config.live_smbh_count > 0) && (ind_act[i] < config.live_smbh_count)) eta_curr = config.eta/config.eta_bh_corr;
             else eta_curr = config.eta;
 
-            double dt_new = aarseth_step(eta_curr, dt_tmp, a_act_new[i], adot_act_new[i], a2, a3);
+            double dt_new = aarseth_step(eta_curr, dt_cur, a_act_new[i], adot_act_new[i], a2, a3);
 
-            //TODO the below should be moved to a function
-            if (dt_new < dt_min) dt_tmp = dt_min;
-            if ((dt_new < dt_tmp) && (dt_new > dt_min)) {
-                int power = log(dt_new)/log(2.0) - 1;
-                dt_tmp = pow(2.0, (double)power); // TODO why is this casting needed here?
-            }
-
-            if ((dt_new > 2*dt_tmp) && (fmod(min_t, 2*dt_tmp) == 0) && (2*dt_tmp <= dt_max)) {
-                dt_tmp *= 2;
-            }
+            dt_new = blockize_step(dt_new, dt_cur, min_t, dt_min, dt_max);
 
             if (config.dt_min_warning && (myRank == 0)) {
-                if (dt_tmp == dt_min) {
-                    printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d time_cur=%.16E\n", dt_tmp, ind_act[i], time_cur);
+                if (dt_new == dt_min) {
+                    printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d time_cur=%.16E\n", dt_cur, ind_act[i], time_cur);
                     fflush(stdout);
                 }
             }
-            if (dt_tmp < min_dt) min_dt = dt_tmp;
+            if (dt_new < min_dt) min_dt = dt_new;
 
             x[j_act] = x_act_new[i];
             v[j_act] = v_act_new[i];
             t[j_act] = min_t;
-            dt[j_act] = dt_tmp;
+            dt[j_act] = dt_new;
             pot[j_act] = pot_act_new[i];
             pot_ext[j_act] = pot_act_ext[i];
             a[j_act] = a_act_new[i];
@@ -614,10 +623,8 @@ int main(int argc, char *argv[])
             if (config.live_smbh_count==2) dt[1] = min_dt;
         }
 
-        if (config.binary_smbh_influence_sphere_output && (myRank == rootRank)) {
-            //TODO clean up this mass. We don't want to have all these _act arrays; they are only needed for this lame function.
+        if (config.binary_smbh_influence_sphere_output && (myRank == rootRank))
             binary_smbh_influence_sphere_output(ind_act, n_act, timesteps, time_cur);
-        }
 
         /* load the new values for active particles to the local GRAPE's */
         for (int i=0; i<n_act; i++) {
@@ -629,7 +636,6 @@ int main(int argc, char *argv[])
                 int j_act = ind_act[i];
                 int address = ind_act[i] - myRank*n_loc;
                 g6_set_j_particle(clusterid, address, ind_act[i], t[j_act], dt[j_act], m[j_act], zeros, adot[j_act]*(1./6.), a[j_act]*0.5, v[j_act], x[j_act]);
-
             } /* if (myRank == cur_rank) */
         } /* i */