diff --git a/phigrape.cpp b/phigrape.cpp
index 70804df..a919f70 100644
--- a/phigrape.cpp
+++ b/phigrape.cpp
@@ -145,6 +145,8 @@ Last redaction : 2019.04.16 12:55
 #include <sys/time.h>
 #include <sys/resource.h>
 
+#include <algorithm>
+
 /*
 double aaa;
 double aaapars[5];
@@ -158,12 +160,12 @@ extern void qwerty_(double *aaa, double *aaapars);
 #include <mpi.h>
 
 /* Some "good" functions and constants... */
-#define SIG(x)   ( ((x)<0) ? (-1):(1) )
-#define ABS(x)   ( ((x)<0) ? (-x):(x) )
-#define MAX(a,b) ( ((a)>(b)) ? (a):(b) )
-#define MIN(a,b) ( ((a)<(b)) ? (a):(b) )
-#define SQR(x)   ( (x)*(x) )
-#define POW3(x)  ( (x)*SQR(x) )
+#define SIG(x)   (((x)<0) ? (-1):(1) )
+#define ABS(x)   (((x)<0) ? (-x):(x) )
+#define MAX(a,b) (((a)>(b)) ? (a):(b) )
+#define MIN(a,b) (((a)<(b)) ? (a):(b) )
+#define SQR(x)   ((x)*(x) )
+#define POW3(x)  ((x)*SQR(x) )
 
 #define Pi          3.14159265358979323846
 #define TWOPi       6.283185307179
@@ -477,7 +479,7 @@ return -   0 if everything OK
 
 double my_rand(void)
 {
-return( (double)(rand()/(double)RAND_MAX) );
+return((double)(rand()/(double)RAND_MAX) );
 }
 
 double my_rand2(void)
@@ -952,101 +954,6 @@ fclose(out);
 
 #endif
 
-
-
-
-void calc_self_grav_zero()
-{
-
-/* calc the grav for the active particles */
-
-g6_set_ti(clusterid, time_cur);
-
-ni = n_act;
-
-/* define the local phi, a, adot for these active particles */
-
-for (i=0; i<ni; i+=npipe)
-    {
-
-    nn = npipe;
-    if (ni-i < npipe) nn = ni - i;
-
-    for (ii=0; ii<nn; ii++)
-    {
-
-    index_i[ii] = ind_act[i+ii];
-    h2_i[ii]    = eps2;
-
-    for (k=0;k<3;k++)
-        {
-        x_i[ii][k] = x_act[i+ii][k]; v_i[ii][k] = v_act[i+ii][k];
-        } /* k */
-
-    p_i[ii] = -1.0;
-
-    for (k=0;k<3;k++)
-        {
-        a_i[ii][k]    = 10.0;
-        jerk_i[ii][k] = 100.0;
-        } /* k */
-
-    } /* ii */
-
-    g6calc_firsthalf(clusterid, n_loc, nn, index_i, x_i, v_i , a_i, jerk_i, p_i, eps2, h2_i);
-    g6calc_lasthalf(clusterid, n_loc, nn, index_i, x_i, v_i, eps2, h2_i, a_i, jerk_i, p_i);
-
-    for (ii=0; ii<nn; ii++)
-        {
-        if (ABS(jerk_i[ii][0]) < 1.0E-05) jerk_i[ii][0] = 1.0E-05;
-        if (ABS(jerk_i[ii][1]) < 1.0E-05) jerk_i[ii][1] = 1.0E-05;
-        if (ABS(jerk_i[ii][2]) < 1.0E-05) jerk_i[ii][2] = 1.0E-05;
-
-        if (ABS(jerk_i[ii][0]) > 100.0) jerk_i[ii][0] = 100.0;
-        if (ABS(jerk_i[ii][1]) > 100.0) jerk_i[ii][1] = 100.0;
-        if (ABS(jerk_i[ii][2]) > 100.0) jerk_i[ii][2] = 100.0;
-        }
-
-    g6calc_firsthalf(clusterid, n_loc, nn, index_i, x_i, v_i , a_i, jerk_i, p_i, eps2, h2_i);
-    g6calc_lasthalf(clusterid, n_loc, nn, index_i, x_i, v_i, eps2, h2_i, a_i, jerk_i, p_i);
-
-    for (ii=0; ii<nn; ii++)
-        {
-
-        pot_act_tmp[i+ii] = p_i[ii];
-
-        for (k=0;k<3;k++)
-        {
-        a_act_tmp[i+ii][k] = a_i[ii][k];
-        adot_act_tmp[i+ii][k] = jerk_i[ii][k];
-        } /* k */
-
-        } /* ii */
-
-    } /* i */
-
-/* Store the value of the local partial force etc... */
-
-for (i=0; i<n_act; i++)
-    {
-
-    pot_act_tmp_loc[i] = pot_act_tmp[i];
-
-    for (k=0;k<3;k++)
-    {
-    a_act_tmp_loc[i][k] = a_act_tmp[i][k];
-    adot_act_tmp_loc[i][k] = adot_act_tmp[i][k];
-    } /* k */
-
-    } /* i */
-
-}
-
-
-
-
-
-
 void calc_ext_grav_zero()
 {
 
@@ -1194,7 +1101,7 @@ for (i=0; i<ni; i++)
     dum3  = POW3(dum);
     dum_g = pow( dum , -g_ext );
 
-    pot_ext[i] -= ( (m_ext/r_ext) / (2.0-g_ext) ) * ( 1.0 - dum2 * dum_g );
+    pot_ext[i] -= ((m_ext/r_ext) / (2.0-g_ext) ) * ( 1.0 - dum2 * dum_g );
 
     tmp = (m_ext/tmp_r3) * dum3 * dum_g;
 
@@ -1205,7 +1112,7 @@ for (i=0; i<ni; i++)
     rv_ij = ( vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij );
 
     tmp = ( m_ext/POW3(tmp_r+r_ext) ) * dum_g;
-    tmp_g = ( (r_ext*g_ext + 3.0*tmp_r)/(tmp_r2*(tmp_r+r_ext)) ) * rv_ij;
+    tmp_g = ((r_ext*g_ext + 3.0*tmp_r)/(tmp_r2*(tmp_r+r_ext)) ) * rv_ij;
 
     adot[i][0] += tmp * ( tmp_g * x_ij - vx_ij );
     adot[i][1] += tmp * ( tmp_g * y_ij - vy_ij );
@@ -1296,17 +1203,15 @@ for (i=0; i<ni; i++)
 
 }
 
-void calc_self_grav()
+void calc_self_grav(double t)
 {
     /* calc the new grav for the active particles */
 
-    // TODO possibly redundant copying in this function
-
 #ifdef TIMING
     get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
 #endif
 
-    g6_set_ti(clusterid, min_t);
+    g6_set_ti(clusterid, t);
 
     ni = n_act;
 
@@ -1316,30 +1221,11 @@ void calc_self_grav()
         nn = npipe;
         if (ni-i < npipe) nn = ni - i;
         for (ii=0; ii<nn; ii++) {
-            iii = ind_act[i+ii];
-            index_i[ii] = iii;
-            h2_i[ii]    = eps2;
-            for (k=0;k<3;k++) {
-                x_i[ii][k] = x_act_new[i+ii][k];
-                v_i[ii][k] = v_act_new[i+ii][k];
-            } /* k */
-            p_i[ii] = pot_act_tmp_loc[iii];
-            for (k=0;k<3;k++) {
-                a_i[ii][k] = a_act_tmp_loc[iii][k];
-                jerk_i[ii][k] = adot_act_tmp_loc[iii][k];
-            } /* k */
+            h2_i[ii]    = eps2; // TODO This should be a global or something
         } /* ii */
-        g6calc_firsthalf(clusterid, n_loc, nn, index_i, x_i, v_i , a_i, jerk_i, p_i, eps2, h2_i);
-        g6calc_lasthalf(clusterid, n_loc, nn, index_i, x_i, v_i, eps2, h2_i, a_i, jerk_i, p_i);
+        g6calc_firsthalf(clusterid, n_loc, nn, ind_act+i, x_act_new+i, v_act_new+i,             a_act_tmp+i, adot_act_tmp+i, pot_act_tmp+i, eps2, h2_i);
+        g6calc_lasthalf( clusterid, n_loc, nn, ind_act+i, x_act_new+i, v_act_new+i, eps2, h2_i, a_act_tmp+i, adot_act_tmp+i, pot_act_tmp+i);
         g6_calls++;
-
-        for (ii=0; ii<nn; ii++) {
-            pot_act_tmp[i+ii]  = p_i[ii];
-            for (k=0;k<3;k++) {
-                a_act_tmp[i+ii][k] = a_i[ii][k];
-                adot_act_tmp[i+ii][k] = jerk_i[ii][k];
-            } /* k */
-        } /* ii */
     } /* i */
 
     /* Store the new value of the local partial force etc... */
@@ -1348,8 +1234,8 @@ void calc_self_grav()
         iii = ind_act[i];
         pot_act_tmp_loc[iii] = pot_act_tmp[i];
         for (k=0;k<3;k++) {
-        a_act_tmp_loc[iii][k] = a_act_tmp[i][k];
-        adot_act_tmp_loc[iii][k] = adot_act_tmp[i][k];
+            a_act_tmp_loc[iii][k] = a_act_tmp[i][k];
+            adot_act_tmp_loc[iii][k] = adot_act_tmp[i][k];
         } /* k */
     } /* i */
 
@@ -1510,7 +1396,7 @@ for (i=0; i<ni; i++)
     dum3  = POW3(dum);
     dum_g = pow( dum , -g_ext );
 
-    pot_act_ext[i] -= ( (m_ext/r_ext) / (2.0-g_ext) ) * ( 1.0 - dum2 * dum_g );
+    pot_act_ext[i] -= ((m_ext/r_ext) / (2.0-g_ext) ) * ( 1.0 - dum2 * dum_g );
 
     tmp = (m_ext/tmp_r3) * dum3 * dum_g;
 
@@ -1521,7 +1407,7 @@ for (i=0; i<ni; i++)
     rv_ij = ( vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij );
 
     tmp = ( m_ext/POW3(tmp_r+r_ext) ) * dum_g;
-    tmp_g = ( (r_ext*g_ext + 3.0*tmp_r)/(tmp_r2*(tmp_r+r_ext)) ) * rv_ij;
+    tmp_g = ((r_ext*g_ext + 3.0*tmp_r)/(tmp_r2*(tmp_r+r_ext)) ) * rv_ij;
 
     adot_act_new[i][0] += tmp * ( tmp_g * x_ij - vx_ij );
     adot_act_new[i][1] += tmp * ( tmp_g * y_ij - vy_ij );
@@ -1735,8 +1621,8 @@ E_tot = E_pot + E_kin + E_pot_ext;
 E_tot_corr    = E_tot + E_corr;
 E_tot_corr_sd = E_tot + E_corr + E_sd;
 
-//  if ( (Timesteps == 0.0) && (time_cur == 0.0) )
-if ( (Timesteps == 0.0) )
+//  if ((Timesteps == 0.0) && (time_cur == 0.0) )
+if (Timesteps == 0.0)
     {
     /* initialize the E_tot_0 + etc... */
 
@@ -1997,7 +1883,7 @@ int main(int argc, char *argv[])
 
         fflush(stdout);
 
-        if ( (diskstep == 0) && (time_cur == 0.0) ) {
+        if ((diskstep == 0) && (time_cur == 0.0)) {
             out = fopen("contr.dat","w");
             fclose(out);
 
@@ -2249,9 +2135,12 @@ int main(int argc, char *argv[])
         t_act[i]  = t[iii];
         dt_act[i] = dt[iii];
 
-        } /* i */
+    } /* i */
 
-    calc_self_grav_zero();
+    // NOTE this is where calc_self_grav_zero() used to be. Hopefully the copying from *_act to *_act_new is just a temporary measure.
+    std::copy(&x_act[0][0], &x_act[0][0]+3*N, &x_act_new[0][0]);
+    std::copy(&v_act[0][0], &v_act[0][0]+3*N, &v_act_new[0][0]);
+    calc_self_grav(time_cur);
 
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
@@ -2604,7 +2493,7 @@ int main(int argc, char *argv[])
         n_act = 0;
 
         for (i=0; i<N; i++) {
-            if ( (t[i]+dt[i]) == min_t ) {
+            if (t[i]+dt[i] == min_t ) {
                 ind_act[n_act] = i;
                 n_act++;
             }
@@ -2702,7 +2591,7 @@ int main(int argc, char *argv[])
         DT_ACT_PRED += (CPU_tmp_user - CPU_tmp_user0);
 #endif
 
-        calc_self_grav();
+        calc_self_grav(min_t);
 
         /* Reduce the "global" vectors from "local" on all the nodes */
 
@@ -2851,7 +2740,7 @@ int main(int argc, char *argv[])
             a3dot1abs = sqrt(a3[0]*a3[0]+a3[1]*a3[1]+a3[2]*a3[2]);
 
 #ifdef ADD_BH2
-            if ( (ind_act[i] == 0) || (ind_act[i] == 1) ) {
+            if ((ind_act[i] == 0) || (ind_act[i] == 1)) {
                 dt_new = sqrt(eta_bh*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
             } else {
                 dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
@@ -2859,7 +2748,7 @@ int main(int argc, char *argv[])
 #else
 
 #ifdef ADD_BH1
-            if ( (ind_act[i] == 0) ) {
+            if (ind_act[i] == 0) {
                 dt_new = sqrt(eta_bh*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
             } else {
                 dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
@@ -2872,7 +2761,7 @@ int main(int argc, char *argv[])
 
             if (dt_new < dt_min) dt_tmp = dt_min;
 
-            if ( (dt_new < dt_tmp) && (dt_new > dt_min) ) {
+            if ((dt_new < dt_tmp) && (dt_new > dt_min)) {
                 power = log(dt_new)/log(2.0) - 1;
 
                 dt_tmp = pow(2.0, (double)power);
@@ -2965,9 +2854,9 @@ int main(int argc, char *argv[])
 
                 iii = ind_act[i];
 
-                if ( (tmp_r2 < SEMI_a2*R_INF2) ) {
+                if (tmp_r2 < SEMI_a2*R_INF2) {
 
-                    if ( (inf_event[iii] == 0) ) {
+                    if (inf_event[iii] == 0) {
 
                         fprintf(out,"INF1 %.6E %.16E %07d %07d   %.6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E \n",
                             Timesteps, time_cur, i, ind_act[i],
@@ -2982,7 +2871,7 @@ int main(int argc, char *argv[])
                     }
 
                 } else {
-                    if ( (inf_event[iii] == 1) ) {
+                    if (inf_event[iii] == 1) {
                         fprintf(out,"INF2 %.6E %.16E %07d %07d   %.6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E \n",
                             Timesteps, time_cur, i, ind_act[i],
                             sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2],