/*****************************************************************************
File Name      : "phi-GRAPE/GPU.c"              // BH (1 || 2) + ACC + EJECT
                :
Contents       : N-body code with integration by individual block time step
                : together with the parallel using of GRAPE6a board's.
                :
                : Added the GPU support via SAPPORO library.
                :
                : Normalization to the physical units!!!
                :
                : External Potential added
                : Plummer-Kuzmin: Bulge, Disk, Halo
                : Kharchenko+Andreas...
                :
                : SC extra POT for Bek SC test runs...
                :
                : Rebuced to the Single BH -> Plummer
                : Andreas+Fazeel...
                :
                : Stellar evolution added
                : Stellar lifetimes: Raiteri, Villata & Navarro (1996)
                : IMS mass loss: van den Hoeg & Groenewegen (1997)
                :
                : STARDESTR_EXT: Tidal disruption of stars by external BH...
                : Chingis, Denis & Maxim...
                :
                : STARDESTR: Tidal disruption of stars by BH...
                : Jose, Li Shuo & Shiyan Zhong
                :
                : STARDISK: Drag force...
                : Chingis, Denis & Maxim...
                :
                : STARDISK: variable hz = HZ*(R/R_crit) up to R_crit...
                : Taras, Andreas...
                :
                : Live BH (1 || 2) + ACC + EJECT...
                : Li Shuo & Shiyan Zhong
                :
                : dt_min for BH (1 || 2)...
                :
                : added the PN calculus for the BBH
                : PN0, PN1, PN2, PN2.5 (coded on the base of
                : Gabor Kupi original routine)
                :
                : added the "name" array...
                :
                : added the GMC's calculus (GMC on CPU; GMC2 on GPU)
                : for Alexey SC runs... and also for Fazeel Zurich runs...
                :
                : CPU_TIMELIMIT added for the Julich MW cluster runs...
                :
Coded by       : Peter Berczik
Version number : 19.04
Last redaction : 2019.04.16 12:55
*****************************************************************************/
#include "config.h"
Config *config;

#define NORM            // Physical normalization

//#define EXTPOT_BH        // BH - usually NB units


//#define EXTPOT_GAL_DEH      // Dehnen Galactic - usually physical units

#ifdef ETICS
#include "grapite.h"
// why do we need CEP as a compilaion flag... just have it always on when ETICS is on. IF there is no CEP, there should be a graceful skipping of those operations.
//#define ETICS_CEP
#ifndef ETICS_DTSCF
#error "ETICS_DTSCF must be defined"
#endif
#endif

#define TIMING

#define ETA_S_CORR   4.0
#define ETA_BH_CORR  4.0

#define DTMAXPOWER  -3.0
#define DTMINPOWER -36.0

/*
-3.0     0.125
-4.0   0.0625
-5.0   0.03125
-7.0     ~1e-2
-10.0    ~1e-3
.............
-20.0     ~1e-6
-23.0    ~1e-7
-26.0    ~1e-8
-30.0    ~1e-9
*/

//#define ACT_DEF_LL

#if defined(ACT_DEF_LL) && defined(ACT_DEF_GRAPITE)
#error "Contradicting preprocessor flags!"
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <sys/time.h>
#include <sys/resource.h>

#include <algorithm>
#include <numeric>
#include <stdexcept>

#define G6_NPIPE 2048
#include "grape6.h"

#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 Pi          3.14159265358979323846
#define TWOPi       6.283185307179
#define sqrt_TWOPi  2.506628274631

#ifdef NORM
//http://pdg.lbl.gov/2015/reviews/rpp2015-rev-astrophysical-constants.pdf

#define G      6.67388E-11        // (m/s^2) * (m^2/kg)
#define Msol   1.988489E+30        // kg
#define Rsol   6.957E+08        // m
#define AU     149597870700.0        // m
#define pc     3.08567758149E+16    // m
#define Year   31556925.2        // s
#define c_feny 299792458.0        // m/s

#define kpc    (1.0E+03*pc)        // m
#define km     1.0E+03            // km   -> m
#define cm3    1.0E-06            // cm^3 -> m^3
#define Myr    (1.0E+06*Year)        // s
#define Gyr    (1.0E+09*Year)        // s
#define R_gas  8.31447215        // J/(K*mol)
#define k_gas  1.380650424E-23        // J/K
#define N_A    6.022141510E+23      // 1/mol
#define mu     1.6605388628E-27        // kg
#define mp     1.67262163783E-27     // kg
#define me     9.1093821545E-31     // kg

#define pc2    (pc*pc)
#define pc3    (pc*pc*pc)
#define kpc2   (kpc*kpc)
#define kpc3   (kpc*kpc*kpc)
#endif

/*
    1KB =    1024
    2KB =    2048
    4KB =    4096
    8KB =    8192
16KB =   16384
32KB =   32768
64KB =   65536
128KB =  131072
256KB =  262144
512KB =  524288
1024KB = 1048576   ->   1MB
*/

#define KB         1024
#define MB         (KB*KB)

#define N_MAX      (6*MB)
#define N_MAX_loc  (2*MB)

struct double3 {
    double data[3];
    double3() {}
    double3(const double x, const double y, const double z)
    {
        data[0] = x;
        data[1] = y;
        data[2] = z;
    }
    double& operator[](int i) {return data[i];}
    double3& operator=(const double3& a)
    {
        data[0] = a.data[0];
        data[1] = a.data[1];
        data[2] = a.data[2];
        return *this;
    }
    double3& operator+=(const double3& a)
    {
        data[0] += a.data[0];
        data[1] += a.data[1];
        data[2] += a.data[2];
        return *this;
    }
    double3& operator-=(const double3& a)
    {
        data[0] -= a.data[0];
        data[1] -= a.data[1];
        data[2] -= a.data[2];
        return *this;
    }
    double3& operator/=(const double& c)
    {
        data[0] /= c;
        data[1] /= c;
        data[2] /= c;
        return *this;
    }
    double norm2() const
    {
        return data[0]*data[0]+data[1]*data[1]+data[2]*data[2];
    }
    double norm() const
    {
        return sqrt(data[0]*data[0]+data[1]*data[1]+data[2]*data[2]);
    }
    operator double*() {return data;}
};

double3 operator*(const double& c, const double3& a)
{
    return double3(a.data[0]*c, a.data[1]*c, a.data[2]*c);
}

double3 operator*(const double3& a, const double& c)
{
    return double3(a.data[0]*c, a.data[1]*c, a.data[2]*c);
}

double operator*(const double3& a, const double3& b)
{
    return a.data[0]*b.data[0]+a.data[1]*b.data[1]+a.data[2]*b.data[2];
}

double3 operator/(const double3& a, const double& c)
{
    return double3(a.data[0]/c, a.data[1]/c, a.data[2]/c);
}

double3 operator+(const double3& a, const double3& b)
{
    return double3(a.data[0]+b.data[0], a.data[1]+b.data[1], a.data[2]+b.data[2]);
}

double3 operator-(const double3& a, const double3& b)
{
    return double3(a.data[0]-b.data[0], a.data[1]-b.data[1], a.data[2]-b.data[2]);
}



double L[3]; // needed in pn_bh_spin.c
// Needed for things related to BHs
#include "debug.h"
// int    ind_sort[N_MAX];
// double var_sort[N_MAX];


double CPU_time_real0, CPU_time_user0, CPU_time_syst0;
double CPU_time_real,  CPU_time_user,  CPU_time_syst;

#ifdef TIMING

double CPU_tmp_real0, CPU_tmp_user0, CPU_tmp_syst0;
double CPU_tmp_real,  CPU_tmp_user,  CPU_tmp_syst;

double DT_TOT,
    DT_ACT_DEF1, DT_ACT_DEF2, DT_ACT_DEF3, DT_ACT_PRED,
    DT_ACT_GRAV, DT_EXT_GRAV,
    DT_GMC_GRAV, DT_GMC_GMC_GRAV, DT_EXT_GMC_GRAV,
    DT_ACT_CORR, DT_ACT_LOAD,
    DT_STEVOL, DT_STARDISK, DT_STARDESTR;

double DT_ACT_REDUCE;

#endif

/* some local settings for G6a board's */

int clusterid, ii, nn, numGPU;

int    npipe=G6_NPIPE, index_i[G6_NPIPE];
double h2_i[G6_NPIPE], p_i[G6_NPIPE];
double3 x_i[G6_NPIPE], v_i[G6_NPIPE],
        a_i[G6_NPIPE], jerk_i[G6_NPIPE];
int    new_tunit=51, new_xunit=51;

double ti=0.0;
double3 a2by18, a1by6, aby2;

/* normalization... */

double m_norm, r_norm, v_norm, t_norm;

double eps_BH=0.0;

/* external potential... */


#ifdef EXTPOT_GAL_DEH
double m_ext, r_ext, g_ext,
    tmp_r2, tmp_r3, dum, dum2, dum3, dum_g, tmp_g;
#endif

#ifdef EXTPOT_BH
double  r2, rv_ij,
        m_bh, b_bh, eps_bh; // NOTE different from eps_BH
        // NOTE there are other m_bh and b_bh defined outside this ifdef block
        // NOTE this is a mess. Looks like eps_bh is never used, and the m_bh from outside the block is never used.
#endif

// double x_bh1[3], x_bh2[3], v_bh1[3], v_bh2[3];
double3 x_bh1, x_bh2, v_bh1, v_bh2;

// double pot_bh1, a_bh1[3], adot_bh1[3],
//     pot_bh2, a_bh2[3], adot_bh2[3];
double pot_bh1, pot_bh2;
double3 a_bh1, adot_bh1, a_bh2, adot_bh2;


//double eps_BH = 0.0;

#include "n_bh.c"

double3 a_pn1[7], adot_pn1[7], a_pn2[7], adot_pn2[7];


#include "pn_bh_spin.c"

#ifdef ETICS
double t_exp, dt_exp=ETICS_DTSCF; // t_exp is just the initial value
#ifdef ETICS_CEP
int grapite_cep_index;
#endif
#endif

void get_CPU_time(double *time_real, double *time_user, double *time_syst)
{
    struct rusage xxx;
    double sec_u, microsec_u, sec_s, microsec_s;
    struct timeval tv;

    getrusage(RUSAGE_SELF,&xxx);

    sec_u = xxx.ru_utime.tv_sec;
    sec_s = xxx.ru_stime.tv_sec;

    microsec_u = xxx.ru_utime.tv_usec;
    microsec_s = xxx.ru_stime.tv_usec;

    *time_user = sec_u + microsec_u * 1.0E-06;
    *time_syst = sec_s + microsec_s * 1.0E-06;

    gettimeofday(&tv, NULL);
    *time_real = tv.tv_sec + 1.0E-06 * tv.tv_usec;

    *time_user = *time_real;
}

void read_data(char inp_fname[], int *diskstep, int *N, double *time_cur, int ind[], double m[], double3 x[], double3 v[])
{
    auto inp = fopen(inp_fname, "r");
    fscanf(inp, "%d \n", diskstep);
    fscanf(inp, "%d \n", N);
    fscanf(inp, "%lE \n", time_cur);
    for (int i=0; i<*N; i++) fscanf(inp,"%d %lE %lE %lE %lE %lE %lE %lE \n", &ind[i], &m[i], &x[i][0], &x[i][1], &x[i][2], &v[i][0], &v[i][1], &v[i][2]);
    fclose(inp);
}

void write_snap_data(char out_fname[], int diskstep, int N, double time_cur, int ind[], double m[], double3 x[], double3 v[])
{
    auto out = fopen(out_fname, "w");
    fprintf(out,"%06d \n", diskstep);
    fprintf(out,"%07d \n", N);
    fprintf(out,"%.10E \n", time_cur);
    for (int i=0; i<N; i++) {
        fprintf(out,"%07d \t %.10E \t % .10E % .10E % .10E \t % .10E % .10E % .10E \n",
            ind[i],
            m[i],
            x[i][0], x[i][1], x[i][2],
            v[i][0], v[i][1], v[i][2]);
        }
    fclose(out);
}

void write_bh_data(double time_cur, double m[], double3 x[], double3 v[], double pot[], double3 a[], double3 adot[], double dt[])
{
    if (config->live_smbh_count == 2) {

        auto out = fopen("bh.dat", "a");

        for (int i=0; i < 2; i++) {
            double3 *a_pn, *adot_pn, a_bh, adot_bh;
            double pot_bh;
            if (i==0) {
                a_pn    = a_pn1;
                adot_pn = adot_pn1;
                pot_bh  = pot_bh1;
                a_bh    = a_bh1;
                adot_bh = adot_bh1;
            } else {
                a_pn    = a_pn2;
                adot_pn = adot_pn2;
                pot_bh  = pot_bh2;
                a_bh    = a_bh2;
                adot_bh = adot_bh2;
            }
            if (config->live_smbh_custom_eps >= 0) {
                fprintf(out, "%.16E \t %.8E \t % .16E % .16E % .16E \t %.16E \t % .16E % .16E % .16E \t %.16E \t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t %.8E \t\t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E ",
                        time_cur, m[i],
                        x[i][0], x[i][1], x[i][2], x[i].norm(),
                        v[i][0], v[i][1], v[i][2], v[i].norm(),
                        pot[i],
                        a[i][0],       a[i][1],       a[i][2],       a[i].norm(),
                        adot[i][0],    adot[i][1],    adot[i][2],    adot[i].norm(),
                        dt[i],
                        pot_bh,
                        a_bh[0],       a_bh[1],       a_bh[2],       a_bh.norm(),
                        adot_bh[0],    adot_bh[1],    adot_bh[2],    adot_bh.norm());
                if (config->binary_smbh_pn) {
                    fprintf(out, "\t");
                    for (int pn_idx=0; pn_idx < 7; pn_idx++) {
                        fprintf(out, "\t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E ", a_pn[pn_idx][0], a_pn[pn_idx][1], a_pn[pn_idx][2], a_pn[pn_idx].norm(), adot_pn[pn_idx][0], adot_pn[pn_idx][1], adot_pn[pn_idx][2], adot_pn[pn_idx].norm());
                    }
                }
                fprintf(out, "\n");
            } else {
                fprintf(out,"%.16E \t %.8E \t % .16E % .16E % .16E \t %.16E \t % .16E % .16E % .16E \t %.16E \t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t %.8E \n",
                        time_cur, m[i],
                        x[i][0], x[i][1], x[i][2], x[i].norm(),
                        v[i][0], v[i][1], v[i][2], v[i].norm(),
                        pot[i],
                        a[i][0], a[i][1], a[i][2], a[i].norm(),
                        adot[i][0], adot[i][1], adot[i][2], adot[i].norm(),
                        dt[i]);
            }
        }
        fprintf(out, "\n");
        fclose(out);
    } else if (config->live_smbh_count == 1) {
        auto out = fopen("bh.dat", "a");
        double tmp_r = sqrt( SQR(x[0][0]) + SQR(x[0][1]) + SQR(x[0][2]) );
        double tmp_v = sqrt( SQR(v[0][0]) + SQR(v[0][1]) + SQR(v[0][2]) );
        double tmp_a = sqrt( SQR(a[0][0]) + SQR(a[0][1]) + SQR(a[0][2]) );
        double tmp_adot = sqrt( SQR(adot[0][0]) + SQR(adot[0][1]) + SQR(adot[0][2]) );
        fprintf(out,"%.16E \t %.8E \t % .16E % .16E % .16E \t %.16E \t % .16E % .16E % .16E \t %.16E \t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t %.8E \n",
                time_cur, m[0],
                x[0][0], x[0][1], x[0][2], tmp_r,
                v[0][0], v[0][1], v[0][2], tmp_v,
                pot[0],
                a[0][0], a[0][1], a[0][2], tmp_a,
                adot[0][0], adot[0][1], adot[0][2], tmp_adot,
                dt[0]);
        fprintf(out,"\n");
        fclose(out);
    }
}

void write_bh_nb_data(double time_cur, int N, double m[], double3 x[], double3 v[])
{
    int nb = config->live_smbh_neighbor_number;
    int    ind_sort[N_MAX];
    double var_sort[N_MAX];

    auto out = fopen("bh_neighbors.dat", "a");

    /* 1st BH */

    for (int i_bh=0; i_bh < config->live_smbh_count; i_bh++) {
        for (int i=0; i<N; i++) var_sort[i] = (x[i]-x[i_bh]).norm();
        std::iota(ind_sort, ind_sort+N, 0);
        std::partial_sort(ind_sort, ind_sort + nb, ind_sort + N, [&](int i, int j) {return var_sort[i] < var_sort[j];});

        fprintf(out,"%.16E \t %07d \t %.8E \t % .8E % .8E % .8E \t % .8E % .8E % .8E \t",
                time_cur,
                i_bh,
                m[i_bh],
                x[i_bh][0], x[i_bh][1], x[i_bh][2],
                v[i_bh][0], v[i_bh][1], v[i_bh][2]);

        for (int j=1; j < nb; j++) {
            int i = ind_sort[j];
            fprintf(out,"%02d %07d   %.8E   % .8E % .8E % .8E   %.8E   % .8E % .8E % .8E   %.8E \t",
                    j, i,
                    m[i],
                    x[i][0], x[i][1], x[i][2], (x[i]-x[i_bh]).norm(),
                    v[i][0], v[i][1], v[i][2], (v[i]-v[i_bh]).norm());
        }

        fprintf(out,"\n");
    }
    fprintf(out, "\n"); // this is redundant
    fclose(out);
}

class External_gravity {
public:
    void apply(const int n_act, const double3 x[], const double3 v[], double pot[], double3 a[], double3 adot[])
    {
        for (int i=0; i<n_act; i++) {
            this->set_coordinates(x[i], v[i]);
            this->calc_gravity();
            pot[i]  += potential;
            a[i]    += acceleration;
            adot[i] += jerk;
        }
    }
    virtual void calc_gravity() = 0;
    bool is_active;
protected:
    double potential;
    double3 acceleration, jerk;
    double3 x, v;
    void set_coordinates(double3 x, double3 v)
    {
        this->x = x;
        this->v = v;
    }
};

class Plummer : public External_gravity {
public:
    Plummer(double m, double b) : m(m), b(b) {is_active=(m>0);}
private:
    void calc_gravity()
    {
        double r2 = SQR(b);
        r2 += x.norm2();
        double r = sqrt(r2);
        double rv_ij = v*x;
        double tmp = m / r;
        potential = -tmp;
        tmp /= r2;
        acceleration = - tmp * x;
        jerk = - tmp * (v - 3*rv_ij * x/r2);
    }
    double m, b;
};

class Miyamoto_Nagai : public External_gravity {
public:
    Miyamoto_Nagai(double m, double a, double b) : m(m), a(a), b(b) {is_active=(m>0);}
private:
    void calc_gravity()
    {
        double  x_ij=x[0],  y_ij=x[1],  z_ij=x[2];
        double vx_ij=v[0], vy_ij=v[1], vz_ij=v[2];
        auto z2_tmp = z_ij*z_ij + SQR(b);
        auto z_tmp  = sqrt(z2_tmp);
        auto r2_tmp = x_ij*x_ij + y_ij*y_ij + SQR(z_tmp + a);
        auto r_tmp  = sqrt(r2_tmp);
        potential = - m / r_tmp;
        auto tmp = m / (r2_tmp*r_tmp);
        acceleration[0] = - tmp * x_ij;
        acceleration[1] = - tmp * y_ij;
        acceleration[2] = - tmp * z_ij * (z_tmp + a)/z_tmp;
        tmp = m / (z_tmp*r2_tmp*r2_tmp*r_tmp);
        jerk[0] = tmp * (- vx_ij*z_tmp*r2_tmp
                + 3.0*x_ij*vx_ij*x_ij*z_tmp
                + 3.0*x_ij*vy_ij*y_ij*z_tmp
                + 3.0*x_ij*vz_ij*z_ij*SQR(z_tmp + a));
        jerk[1] = tmp * (- vy_ij*z_tmp*r2_tmp
                + 3.0*y_ij*vx_ij*x_ij*z_tmp
                + 3.0*y_ij*vy_ij*y_ij*z_tmp
                + 3.0*y_ij*vz_ij*z_ij*SQR(z_tmp + a));
        jerk[2] = tmp * (- vz_ij*(z_tmp + a)*(x_ij*x_ij*(z2_tmp*z_tmp + a*SQR(b)) +
                y_ij*y_ij*(z2_tmp*z_tmp + a*SQR(b)) -
                (2.0*a*(SQR(z_ij) - SQR(b))*z_tmp +
                2.0*SQR(z_ij*z_ij) +
                SQR(b)*SQR(z_ij) -
                SQR(b)*(SQR(a) + SQR(b))))
                + 3.0*vx_ij*x_ij*z_ij*z2_tmp*(z_tmp + a)
                + 3.0*vy_ij*y_ij*z_ij*z2_tmp*(z_tmp + a)) / z2_tmp;
    }
    double m, a, b;
};

class Logarithmic_halo : public External_gravity {
public:
    Logarithmic_halo(double v_halo, double r_halo) : v2_halo(v_halo*v_halo), r2_halo(r_halo*r_halo) {is_active=(r_halo>0);}
private:
    void calc_gravity()
    {
        auto r2 = x.norm2();
        auto rv_ij = 2.0*(v*x);
        auto r2_r2_halo = (r2 + r2_halo);
        potential = - 0.5*v2_halo * log(1.0 + r2/r2_halo);
        auto tmp = v2_halo/r2_r2_halo;
        acceleration = - tmp * x;
        tmp /= (r2 + r2_halo);
        jerk = tmp * (rv_ij * x - r2_r2_halo * v);
    }
    double v2_halo, r2_halo;
};

void calc_self_grav(double t, double eps2, double &g6_calls, int n_loc,
                    int n_act, int ind_act[],
                    double3 x_act_new[], double3 v_act_new[],
                    double pot_act_tmp[],
                    double3 a_act_tmp[],
                    double3 adot_act_tmp[],
                    double h2_i[])
{
    /* calc the new grav for the active particles */

#ifdef TIMING
    get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif

    g6_set_ti(clusterid, t);

    int ni = n_act; // TODO why is this needed?

    /* define the local phi, a, adot for these active particles */

    for (int i=0; i<ni; i+=npipe) {
        nn = npipe;
        if (ni-i < npipe) nn = ni - i;
        for (ii=0; ii<nn; ii++) {
            h2_i[ii]    = eps2; // TODO This should be a global or something
        } /* ii */
        //TODO any way we can clean up this ugly casting?
        g6calc_firsthalf(clusterid, n_loc, nn, ind_act+i, (double(*)[3])x_act_new+i, (double(*)[3])v_act_new+i, (double(*)[3])a_act_tmp+i, (double(*)[3])adot_act_tmp+i, pot_act_tmp+i, eps2, h2_i);
        g6calc_lasthalf( clusterid, n_loc, nn, ind_act+i, (double(*)[3])x_act_new+i, (double(*)[3])v_act_new+i, eps2, h2_i, (double(*)[3])a_act_tmp+i, (double(*)[3])adot_act_tmp+i, pot_act_tmp+i);
        g6_calls++;
    } /* i */

#ifdef TIMING
    get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
    DT_ACT_GRAV += (CPU_tmp_user - CPU_tmp_user0);
#endif
}

void calc_ext_grav(std::vector<External_gravity*> &external_gravity_components, int n_act, double3 *x_act_new, double3 *v_act_new, double *pot_act_ext, double3 *a_act_new, double3* adot_act_new)
{

/* Define the external potential for all active particles on all nodes */

int ni = n_act; // TODO redundant?

    std::fill(pot_act_ext, pot_act_ext+n_act, 0.);

    for (auto component : external_gravity_components) {
        if (component->is_active)
            component->apply(n_act, x_act_new, v_act_new, pot_act_ext, a_act_new, adot_act_new);
    }

#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif

#ifdef EXTPOT_GAL_DEH

for (i=0; i<ni; i++)
    {

    x_ij = x_act_new[i][0];
    y_ij = x_act_new[i][1];
    z_ij = x_act_new[i][2];

    vx_ij = v_act_new[i][0];
    vy_ij = v_act_new[i][1];
    vz_ij = v_act_new[i][2];

    tmp_r = sqrt( SQR(x_ij) + SQR(y_ij) + SQR(z_ij) );
    tmp_r2 = SQR(tmp_r);
    tmp_r3 = POW3(tmp_r);

    dum   = tmp_r/(tmp_r + r_ext);
    dum2  = SQR(dum);
    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 );

    tmp = (m_ext/tmp_r3) * dum3 * dum_g;

    a_act_new[i][0] -= tmp * x_ij;
    a_act_new[i][1] -= tmp * y_ij;
    a_act_new[i][2] -= tmp * z_ij;

    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;

    adot_act_new[i][0] += tmp * ( tmp_g * x_ij - vx_ij );
    adot_act_new[i][1] += tmp * ( tmp_g * y_ij - vy_ij );
    adot_act_new[i][2] += tmp * ( tmp_g * z_ij - vz_ij );

    } /* i */

#endif

    /* For simple Plummer potential... */

#ifdef EXTPOT_BH

for (i=0; i<ni; i++)
    {

    x_ij = x_act_new[i][0];
    y_ij = x_act_new[i][1];
    z_ij = x_act_new[i][2];

    vx_ij = v_act_new[i][0];
    vy_ij = v_act_new[i][1];
    vz_ij = v_act_new[i][2];

    r2 = SQR(b_bh);

    r2 += x_ij*x_ij + y_ij*y_ij + z_ij*z_ij;   r = sqrt(r2);

    rv_ij = vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij;

    tmp = m_bh / r;

    pot_act_ext[i] -= tmp;

    tmp /= r2;

    a_act_new[i][0] -= tmp * x_ij;
    a_act_new[i][1] -= tmp * y_ij;
    a_act_new[i][2] -= tmp * z_ij;

    adot_act_new[i][0] -= tmp * (vx_ij - 3.0*rv_ij * x_ij/r2);
    adot_act_new[i][1] -= tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
    adot_act_new[i][2] -= tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);

    } /* i */

#endif

#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_EXT_GRAV += (CPU_tmp_user - CPU_tmp_user0);
#endif

}

void energy_contr(const double time_cur, const double timesteps, const double n_act_sum, const double g6_calls, double& rcm_sum, double& vcm_sum, double& E_tot_0, double& E_tot_corr_0, double& E_tot_corr_sd_0, int &skip_con, int N, double m[], double3 x[], double3 v[], double pot[], double pot_ext[])
{
    // TODO maybe use static variables here for the previous step energy?
    double E_pot = 0.0;
    for (int i=0; i<N; i++) E_pot += m[i]*pot[i];
    E_pot *= 0.5;

    double E_kin = 0.0;
    for (int i=0; i<N; i++) E_kin += m[i]*( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
    E_kin *= 0.5;

    double E_pot_ext = 0.0;

    for (int i=0; i<N; i++) E_pot_ext += m[i]*pot_ext[i];

    if (timesteps == 0.0) { //TODO what is this thing?
        rcm_sum = 0.0;
        vcm_sum = 0.0;
    }

    double mcm     = 0.0;
    double rcm_mod = 0.0;
    double vcm_mod = 0.0;

    double3 xcm = {0, 0, 0};
    double3 vcm = {0, 0, 0};

    for (int i=0; i<N; i++) {
        mcm += m[i];
        xcm += m[i] * x[i];
        vcm += m[i] * v[i];
    } /* i */

    xcm /= mcm;
    vcm /= mcm;

    rcm_mod = xcm.norm();
    vcm_mod = vcm.norm();

    rcm_sum += rcm_mod;
    vcm_sum += vcm_mod;

    double3 mom = {0, 0, 0};

    for (int i=0; i<N; i++) {
        mom[0] += m[i] * ( x[i][1] * v[i][2] - x[i][2] * v[i][1] );
        mom[1] += m[i] * ( x[i][2] * v[i][0] - x[i][0] * v[i][2] );
        mom[2] += m[i] * ( x[i][0] * v[i][1] - x[i][1] * v[i][0] );
    } /* i */

    get_CPU_time(&CPU_time_real, &CPU_time_user, &CPU_time_syst);

    double E_tot = E_pot + E_kin + E_pot_ext;
    double E_tot_corr    = E_tot /* + E_corr */; // TODO what is this?

    if (timesteps == 0.0) {
        /* initialize the E_tot_0 + etc... */

        E_tot_0 = E_tot;
        E_tot_corr_0 = E_tot;
        E_tot_corr_sd_0 = E_tot;

        skip_con = 1;
    } else {
        /* read the E_tot_0 + etc... */
        if (skip_con == 0) {
            auto inp = fopen("contr.con","r");
            double tmp;
            fscanf(inp,"%lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE",
                   &tmp, &tmp, &tmp, &tmp,
                   &tmp, &tmp, &tmp, &tmp, &tmp, &tmp, &tmp,
                   &E_tot_0, &tmp,
                   &E_tot_corr_0, &tmp,
                   &E_tot_corr_sd_0, &tmp,
                   &tmp, &tmp, &tmp, &tmp, &tmp, &tmp, &tmp, &tmp);
            fclose(inp);
            skip_con = 1;
        }
    }
    double DE_tot         = E_tot         - E_tot_0;

    /* This is the only output to screen */
    printf("%.3E %.3E  % .4E %.4E % .4E  % .4E % .4E  %.2E\n",
            time_cur, timesteps,
            E_pot, E_kin, E_pot_ext, E_tot, DE_tot,
            CPU_time_user-CPU_time_user0);

    fflush(stdout);

    auto out = fopen("contr.dat", "a");
    fprintf(out,"%.8E \t %.8E %.8E %.8E \t % .8E % .8E % .8E % .8E % .8E \t % .8E % .8E \t % .8E % .8E % .8E \t %.8E %.8E %.8E \n",
            time_cur, timesteps, n_act_sum, g6_calls,
            E_pot, E_kin, E_pot_ext,
            E_tot, DE_tot,
            rcm_mod, vcm_mod,
            mom[0], mom[1], mom[2],
            CPU_time_real-CPU_time_real0, CPU_time_user-CPU_time_user0, CPU_time_syst-CPU_time_syst0);
    fclose(out);

    out = fopen("contr.con", "w");
    fprintf(out,"%.8E \t %.8E %.8E %.8E \t % .8E % .8E % .8E % .8E % .8E \t % .8E % .8E \t % .8E % .8E % .8E \t %.8E %.8E %.8E \n",
            time_cur, timesteps, n_act_sum, g6_calls,
            E_pot, E_kin, E_pot_ext,
            E_tot, DE_tot,
            rcm_mod, vcm_mod,
            mom[0], mom[1], mom[2],
            CPU_time_real-CPU_time_real0, CPU_time_user-CPU_time_user0, CPU_time_syst-CPU_time_syst0);
    fclose(out);

    E_tot_0 = E_tot;
    E_tot_corr_0 = E_tot_corr;
    E_tot_corr_sd_0 = 0;
}

int main(int argc, char *argv[])
{
    int name_proc, n_proc=1, myRank=0, rootRank=0, cur_rank,
        nj, diskstep=0, power, jjj, iii,
        skip_con=0, tmp_i;

    double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
           dt_bh, t_bh=0.0, dt_bh_tmp,
           t_end, time_cur, dt_min, dt_max, min_t, min_t_loc, dt_new,
           eta_s, eta, eta_bh,
           E_tot_0, E_tot_corr_0, E_tot_corr_sd_0,
           rcm_sum=0.0, vcm_sum=0.0,
           eps=0.0, eps2,
           a2_mod, adot2_mod,
           dt_tmp, dt2half, dt3over6, dt4over24, dt5over120,
           dtinv, dt2inv, dt3inv,
           a1abs, adot1abs, a2dot1abs, a3dot1abs,
           timesteps=0.0, n_act_sum=0.0, n_act_distr[N_MAX], g6_calls=0.0, g6_calls_sum=0.0,
           tmp;

    double3 xcm, vcm, mom,
            xdc, vdc,
            a2, a3, a2dot1;

    char processor_name[MPI_MAX_PROCESSOR_NAME],
         inp_fname[30];

    /* global variables */

    int    N, ind[N_MAX];
    double m[N_MAX], pot[N_MAX], t[N_MAX], dt[N_MAX];
    double3 x[N_MAX], v[N_MAX], a[N_MAX], adot[N_MAX];

    /* local variables */

    int    n_loc, ind_loc[N_MAX_loc];
    double m_loc[N_MAX_loc], pot_loc[N_MAX_loc], t_loc[N_MAX_loc], dt_loc[N_MAX_loc];
    double3 x_loc[N_MAX_loc], v_loc[N_MAX_loc],
            a_loc[N_MAX_loc], adot_loc[N_MAX_loc];

    /* data for active particles */

    int    n_act, ind_act[N_MAX];
    double m_act[N_MAX],
           pot_act[N_MAX], t_act[N_MAX], dt_act[N_MAX],
           pot_act_new[N_MAX],
           pot_act_tmp[N_MAX];

    double3 x_act[N_MAX], v_act[N_MAX],
            a_act[N_MAX], adot_act[N_MAX],
            x_act_new[N_MAX], v_act_new[N_MAX],
            a_act_new[N_MAX], adot_act_new[N_MAX],
            a_act_tmp[N_MAX], adot_act_tmp[N_MAX];;

    FILE   *out;

    double C_NB;
    int    usedOrNot[7];

    double pot_ext[N_MAX], pot_act_ext[N_MAX];    // for all EXTPOT

    int    i_bh1, i_bh2;

    double m_bh1, m_bh2;

    int    inf_event[N_MAX];
    double DR2, tmp_r2;
    double3 x_bbhc, v_bbhc;
    double DV2, EB, SEMI_a, SEMI_a2;

    double s_bh1[3] = {0.0, 0.0, 1.0};
    double s_bh2[3] = {0.0, 0.0, 1.0};


    /* INIT the rand() !!! */
    srand(19640916);                 /* it is just my birthday :-) */

    /* Init MPI */
    MPI_Init(&argc, &argv);

    /* Define the total number of processors */
    MPI_Comm_size(MPI_COMM_WORLD, &n_proc);

    /* Define of the Rank of each processors */
    MPI_Comm_rank(MPI_COMM_WORLD, &myRank);

    /* Define the processors names */
    MPI_Get_processor_name(processor_name, &name_proc);

    /* Print the Rank and the names of processors */
    printf("Rank of the processor %03d on %s \n", myRank, processor_name);

    /* Wait to all processors to finish his works... */
    MPI_Barrier(MPI_COMM_WORLD);

    /* read the input parameters to the rootRank */

    config = new Config("phigrape.conf");
    C_NB = config->pn_c;
    std::copy(config->pn_usage.begin(), config->pn_usage.end(), usedOrNot);

    if (myRank == rootRank) {

        //TODO move it out of (myRank == rootRank) so you don't need to communicate them.
        eps = config->eps;
        t_end = config->t_end;
        dt_disk = config->dt_disk;
        dt_contr = config->dt_contr;
        dt_bh = config->dt_bh;
        eta = config->eta;
        strcpy(inp_fname, config->input_file_name.c_str());

        if (config->binary_smbh_influence_sphere_output) for (int i=0; i<N; i++) inf_event[i] = 0; // WARNING N wasn't set yet!

        /*
        eps      : Plummer softening parameter (can be even 0)
        t_end    : end time of calculation
        dt_disk  : interval of snapshot files output (0xxx.dat)
        dt_contr : interval for the energy control output (contr.dat)
        dt_bh    : interval for BH output (bh.dat & bh_neighbors.dat)
        eta      : parameter for timestep determination
        inp_data : name of the input file (data.inp)
        */

        /* Read the data for EXTernal POTential */

#ifdef EXTPOT_BH
        fscanf(inp,"%lE %lE", &m_bh, &b_bh);
#endif

#ifdef EXTPOT_GAL_DEH
        fscanf(inp,"%lE %lE %lE", &m_ext, &r_ext, &g_ext);
#endif

        /* read the global data for particles to the rootRank */

        read_data(inp_fname, &diskstep, &N, &time_cur, ind, m, x, v);

        /* possible coordinate & velocity limits for ALL particles !!! */

        printf("\n");
        printf("Begin the calculation of phi-GRAPE program on %03d processors\n", n_proc);
        printf("\n");
        printf("N       = %07d \t eps      = %.6E \n", N, eps);
        printf("t_beg   = %.6E \t t_end    = %.6E \n", time_cur, t_end);
        printf("dt_disk = %.6E \t dt_contr = %.6E \n", dt_disk, dt_contr);
        printf("dt_bh   = %.6E \n", dt_bh);
        printf("eta     = %.6E \n", eta);
        printf("\n");

#ifdef NORM
        printf("Normalization: \n");
        printf("\n");
        printf("m_norm = %.4E [Msol]  r_norm = %.4E [pc] \n", m_norm/Msol, r_norm/pc);
        printf("v_morm = %.4E [km/s]  t_morm = %.4E [Myr] \n", v_norm/km, t_norm/Myr);
        printf("\n");
#endif

        printf("External Potential: \n");
        printf("\n");

#ifdef EXTPOT_BH
        printf("m_bh = %.4E  b_bh = %.4E \n", m_bh, b_bh);
#endif

#ifdef EXTPOT_GAL_DEH
        printf("m_ext  = %.6E r_ext  = %.6E \t g_ext = %.3E \n", m_ext, r_ext, g_ext);
#endif

        if ((diskstep == 0) && (time_cur == 0.0)) {
            out = fopen("contr.dat", "w");
            fclose(out);

#ifdef TIMING
            out = fopen("timing.dat", "w");
            fclose(out);
#endif

            if (config->live_smbh_output && (config->live_smbh_count > 0)) {
                out = fopen("bh.dat", "w");
                fclose(out);
            }
            if ((config->live_smbh_neighbor_output) && (config->live_smbh_count > 0)) {
                out = fopen("bh_neighbors.dat", "w");
                fclose(out);
            }
            if (config->binary_smbh_influence_sphere_output) {
                out = fopen("bbh_inf.dat", "w");
                fclose(out);
            }

        } else {   // if (diskstep == 0)

        }        // if (diskstep == 0)

        get_CPU_time(&CPU_time_real0, &CPU_time_user0, &CPU_time_syst0);

    } /* if (myRank == rootRank) */

    /* Wait to all processors to finish his works... */
    MPI_Barrier(MPI_COMM_WORLD);

    /* Broadcast all useful values to all processors... */
    MPI_Bcast(&N,        1, MPI_INT,    rootRank, MPI_COMM_WORLD);
    MPI_Bcast(&eps,      1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Bcast(&eta,      1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Bcast(&t_end,    1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Bcast(&dt_disk,  1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Bcast(&dt_contr, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Bcast(&dt_bh,    1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Bcast(&time_cur, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);

#ifdef NORM
    MPI_Bcast(&m_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Bcast(&r_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Bcast(&v_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Bcast(&t_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
#endif        // NORM

#ifdef EXTPOT_BH
    MPI_Bcast(&m_bh, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Bcast(&b_bh, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
#endif

#ifdef EXTPOT_GAL_DEH
    MPI_Bcast(&m_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Bcast(&r_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Bcast(&g_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
#endif

    /* Wait to all processors to finish his works... */
    MPI_Barrier(MPI_COMM_WORLD);

    double normalization_mass=1, normalization_length=1, normalization_velocity=1;
    if (config->ext_units_physical) {
        normalization_mass = 1/config->unit_mass;
        normalization_length = 1000/config->unit_length;
        normalization_velocity = 1.52484071426404437233e+01*sqrt(config->unit_length/config->unit_mass);
    }
    Plummer ext_bulge(config->ext_m_bulge*normalization_mass, config->ext_b_bulge*normalization_length);
    Miyamoto_Nagai ext_disk(config->ext_m_disk*normalization_mass, config->ext_a_disk*normalization_length, config->ext_b_disk*normalization_length);
    Plummer ext_halo_plummer(config->ext_m_halo_plummer*normalization_mass, config->ext_b_halo_plummer*normalization_length);
    Logarithmic_halo ext_log_halo(config->ext_log_halo_v*normalization_velocity, config->ext_log_halo_r*normalization_length);

    std::vector<External_gravity*> external_gravity_components;
    external_gravity_components.push_back(&ext_bulge);
    external_gravity_components.push_back(&ext_disk);
    external_gravity_components.push_back(&ext_halo_plummer);
    external_gravity_components.push_back(&ext_log_halo);

    if (ext_bulge.is_active)        printf("m_bulge = %.4E                        b_bulge = %.4E\n", config->ext_m_bulge*normalization_mass, config->ext_b_bulge*normalization_length);
    if (ext_disk.is_active)         printf("m_disk  = %.4E  a_disk  = %.4E  b_disk  = %.4E\n", config->ext_m_disk*normalization_mass, config->ext_a_disk*normalization_length, config->ext_b_disk*normalization_length);
    if (ext_halo_plummer.is_active) printf("m_halo  = %.4E                        b_halo = %.4E\n", config->ext_m_halo_plummer*normalization_mass, config->ext_b_halo_plummer*normalization_length);
    if (ext_log_halo.is_active)     printf("v_halo = %.6E r_halo = %.6E \n", config->ext_log_halo_v*normalization_velocity, config->ext_log_halo_r*normalization_length);
    printf("\n");
    fflush(stdout);

    eta_s  = eta/ETA_S_CORR;
    eta_bh = eta/ETA_BH_CORR;

    eps2 = SQR(eps);

    dt_min = 1.0*pow(2.0, DTMINPOWER);
    dt_max = 1.0*pow(2.0, DTMAXPOWER);

    if (dt_disk == dt_contr)
        dt_max = dt_contr;
    else
        dt_max = MIN(dt_disk, dt_contr);

    if (dt_max > 1.0) dt_max = 1.0;

    t_disk  = dt_disk*(1.0+floor(time_cur/dt_disk));
    t_contr = dt_contr*(1.0+floor(time_cur/dt_contr));
    t_bh    = dt_bh*(1.0+floor(time_cur/dt_bh));

    if (myRank == rootRank) {
        printf("t_disk = %.6E   t_contr = %.6E   t_bh = %.6E \n", t_disk, t_contr, t_bh);
        printf("\n");
        fflush(stdout);
    } /* if (myRank == rootRank) */

    std::fill(t,  t+N,  time_cur);
    std::fill(dt, dt+N, dt_min);

    n_loc = N/n_proc;

    /* Wait to all processors to finish his works... */
    MPI_Barrier(MPI_COMM_WORLD);

    /* Broadcast the values of all particles to all processors... */
    MPI_Bcast(ind, N, MPI_INT,    rootRank, MPI_COMM_WORLD);
    MPI_Bcast(m,   N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Bcast(x, 3*N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Bcast(v, 3*N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);

    /* Wait to all processors to finish his works... */
    MPI_Barrier(MPI_COMM_WORLD);

    /* Scatter the "local" vectors from "global" */
    MPI_Scatter(ind, n_loc, MPI_INT,    ind_loc, n_loc, MPI_INT,    rootRank, MPI_COMM_WORLD);
    MPI_Scatter(m,   n_loc, MPI_DOUBLE, m_loc,   n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Scatter(t,   n_loc, MPI_DOUBLE, t_loc,   n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Scatter(dt,  n_loc, MPI_DOUBLE, dt_loc,  n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Scatter(x, 3*n_loc, MPI_DOUBLE, x_loc, 3*n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Scatter(v, 3*n_loc, MPI_DOUBLE, v_loc, 3*n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);

    /* Wait to all processors to finish his works... */
    MPI_Barrier(MPI_COMM_WORLD);

    /* init the local GRAPE's */

#ifdef MPI_OVERRIDE
    numGPU = 1; // TODO get this from config file
    clusterid = myRank % numGPU;
#else
    MPI_Comm shmcomm;
    MPI_Comm_split_type(MPI_COMM_WORLD, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL, &shmcomm);
    MPI_Comm_size(shmcomm, &numGPU);
    MPI_Comm_rank(shmcomm, &clusterid);
#endif
    printf("Rank of the processor %03d : Number of GPUs %01d : Cluster ID %01d \n", myRank, numGPU, clusterid);
    fflush(stdout);

    g6_open(clusterid);
    npipe = g6_npipes();
    g6_set_tunit(new_tunit);
    g6_set_xunit(new_xunit);

#ifdef ETICS
    grapite_read_particle_tags(N, "grapite.mask", myRank, n_loc);
    grapite_set_dt_exp(dt_exp);
    dt_exp = time_cur; // if we don't have a binary restart then we don't remember the coefficients, and we need to calculate them now.
    grapite_set_t_exp(t_exp);
#endif

    /* initial load the particles to the local GRAPE's */

    nj=n_loc;

    /* load the nj particles to the G6 */

    for (int j=0; j<nj; j++) {
        a2by18 = {0, 0, 0};
        a1by6  = {0, 0, 0};
        aby2   = {0, 0, 0};
        g6_set_j_particle(clusterid, j, ind_loc[j], t_loc[j], dt_loc[j], m_loc[j], a2by18, a1by6, aby2, v_loc[j], x_loc[j]);
    } /* j */

#ifdef ETICS
    double etics_length_scale;
    if (myRank == rootRank) etics_length_scale = grapite_get_length_scale(N, m, x, v); // We don't want all ranks to do it, because they need to write a file and might confuse each other
    MPI_Bcast(&etics_length_scale, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    grapite_set_length_scale(etics_length_scale);
#endif

#ifdef ETICS_CEP
    // First time only: get the CEP index
    grapite_cep_index = grapite_get_cep_index();

    // First calculate the DC
    grapite_calc_center(N, m, x, v, xcm, vcm, xdc, vdc);

    // Now copy it to the global particle list
    memcpy(x[grapite_cep_index], xdc, 3*sizeof(double));
    memcpy(v[grapite_cep_index], vdc, 3*sizeof(double));

    // Now copy it to the local particle list for tha appropriate rank
    if (myRank == grapite_cep_index/n_loc) {
        memcpy(x_loc[grapite_cep_index - myRank*n_loc], xdc, 3*sizeof(double));
        memcpy(v_loc[grapite_cep_index - myRank*n_loc], vdc, 3*sizeof(double));
    }

    double zeros[3] = {0., 0., 0.}; // We haven't calculated the force yet!
    grapite_update_cep(time_cur, xdc, vdc, zeros, zeros);
#endif

    /* define the all particles as a active on all the processors for the first time grav calc. */

    n_act = N;

    for (int i=0; i<n_act; i++) {
        ind_act[i]  = ind[i];
        iii         = ind_act[i];
        m_act[i]    = m[iii];

        x_act[i] = x[iii];
        v_act[i] = v[iii];

        t_act[i]  = t[iii];
        dt_act[i] = dt[iii];

    } /* i */

    // NOTE this is where calc_self_grav_zero() used to be.
    calc_self_grav(time_cur, eps2, g6_calls, n_loc,
                   n_act, ind_act,
                   x_act, v_act,
                   pot_act_tmp,
                   a_act_tmp,
                   adot_act_tmp,
                   h2_i);

    /* Wait to all processors to finish his works... */
    MPI_Barrier(MPI_COMM_WORLD);

    /* Reduce the "global" vectors from "local" on all processors) */
    MPI_Allreduce(pot_act_tmp,  pot,    n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(a_act_tmp,    a,    3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(adot_act_tmp, adot, 3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

    /* Wait to all processors to finish his works... */
    MPI_Barrier(MPI_COMM_WORLD);

    if (config->live_smbh_count == 2) {
        i_bh1 = 0;
        i_bh2 = 1;

        if (config->live_smbh_custom_eps >= 0) {

            m_bh1 = m[i_bh1];
            m_bh2 = m[i_bh2];

            x_bh1 = x[i_bh1];
            v_bh1 = v[i_bh1];

            x_bh2 = x[i_bh2];
            v_bh2 = v[i_bh2];

            // calculate and "minus" the BH <-> BH _softened_ pot, acc & jerk

            tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1,
                                    m_bh2, x_bh2, v_bh2,
                                    eps,
                                    &pot_bh1, a_bh1, adot_bh1,
                                    &pot_bh2, a_bh2, adot_bh2);

            pot[i_bh1] -= pot_bh1;
            pot[i_bh2] -= pot_bh2;

            a[i_bh1] -= a_bh1;
            a[i_bh2] -= a_bh2;

            adot[i_bh1] -= adot_bh1;
            adot[i_bh2] -= adot_bh2;

            // calculate and "plus" the new BH <-> BH _unsoftened_ pot, acc, jerk

            tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1,
                                    m_bh2, x_bh2, v_bh2,
                                    config->live_smbh_custom_eps,
                                    &pot_bh1, a_bh1, adot_bh1,
                                    &pot_bh2, a_bh2, adot_bh2);

            pot[i_bh1] += pot_bh1;
            pot[i_bh2] += pot_bh2;

            a[i_bh1] += a_bh1;
            a[i_bh2] += a_bh2;

            adot[i_bh1] += adot_bh1;
            adot[i_bh2] += adot_bh2;
        }

        if (config->binary_smbh_pn) {

            // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk

            dt_bh_tmp = dt[0];

            tmp_i = calc_force_pn_BH(m_bh1, x_bh1, v_bh1, s_bh1,
                                    m_bh2, x_bh2, v_bh2, s_bh2,
                                    C_NB,  dt_bh_tmp, usedOrNot,
                                    (double(*)[3])a_pn1, (double(*)[3])adot_pn1,
                                    (double(*)[3])a_pn2, (double(*)[3])adot_pn2, myRank, rootRank);

            a[i_bh1] += a_pn1[1] + a_pn1[2] + a_pn1[3] + a_pn1[4] + a_pn1[5] + a_pn1[6];
            a[i_bh2] += a_pn2[1] + a_pn2[2] + a_pn2[3] + a_pn2[4] + a_pn2[5] + a_pn2[6];

            adot[i_bh1] += adot_pn1[1] + adot_pn1[2] + adot_pn1[3] + adot_pn1[4] + adot_pn1[5] + adot_pn1[6];
            adot[i_bh2] += adot_pn2[1] + adot_pn2[2] + adot_pn2[3] + adot_pn2[4] + adot_pn2[5] + adot_pn2[6];

            if (myRank == rootRank) {
                if (tmp_i == 505) {
                    printf("PN RSDIST: %.8E \t %.8E \n", timesteps, time_cur);
                    fflush(stdout);
                    exit(-1);
                }
            }
        }
    }

    calc_ext_grav(external_gravity_components, n_act, x_act, v_act, pot_ext, a, adot);

    /* Wait to all processors to finish his works... */
    MPI_Barrier(MPI_COMM_WORLD);

    /* Wait to all processors to finish his works... */
    MPI_Barrier(MPI_COMM_WORLD);

    /* Energy control... */
    if (myRank == rootRank) {
        energy_contr(time_cur, timesteps, n_act_sum, g6_calls, rcm_sum, vcm_sum, E_tot_0, E_tot_corr_0, E_tot_corr_sd_0, skip_con, N, m, x, v, pot, pot_ext);
    } /* if (myRank == rootRank) */

#ifdef ETICS_DUMP
    if (diskstep==0) {
        sprintf(out_fname, "coeffs.%06d.%02d.dat", 0, myRank);
        grapite_dump(out_fname, 2);
    }
#endif

    /* Wait to all processors to finish his works... */
    MPI_Barrier(MPI_COMM_WORLD);

    /* Scatter the "local" vectors from "global" */
    MPI_Scatter(pot,    n_loc, MPI_DOUBLE, pot_loc,    n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Scatter(a,    3*n_loc, MPI_DOUBLE, a_loc,    3*n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
    MPI_Scatter(adot, 3*n_loc, MPI_DOUBLE, adot_loc, 3*n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);

#ifdef ETICS_CEP
    // First calculate the DC
    grapite_calc_center(N, m, x, v, xcm, vcm, xdc, vdc);

    // Now copy it to the global particle list
    memcpy(x[grapite_cep_index], xdc, 3*sizeof(double));
    memcpy(v[grapite_cep_index], vdc, 3*sizeof(double));

    // Now copy it to the local particle list for tha appropriate rank
    if (myRank == grapite_cep_index/n_loc) {
        memcpy(x_loc[grapite_cep_index - myRank*n_loc], xdc, 3*sizeof(double));
        memcpy(v_loc[grapite_cep_index - myRank*n_loc], vdc, 3*sizeof(double));
    }

    grapite_update_cep(time_cur, xdc, vdc, a[grapite_cep_index], adot[grapite_cep_index]);
#endif

    /* Wait to all processors to finish his works... */
    MPI_Barrier(MPI_COMM_WORLD);

    /* Define initial timestep for all particles on all nodes */

    for (int i=0; i<N; i++) {
        a2_mod = a[i][0]*a[i][0]+a[i][1]*a[i][1]+a[i][2]*a[i][2];
        adot2_mod = adot[i][0]*adot[i][0]+adot[i][1]*adot[i][1]+adot[i][2]*adot[i][2];

        if (adot2_mod == 0)
            dt_tmp = eta_s;
        else
            dt_tmp = eta_s*sqrt(a2_mod/adot2_mod);

        power = log(dt_tmp)/log(2.0) - 1;

        dt_tmp = pow(2.0, (double)power);

        if (dt_tmp < dt_min) dt_tmp = dt_min;
        if (dt_tmp > dt_max) dt_tmp = dt_max;

        dt[i] = dt_tmp;

        if (config->dt_min_warning && (myRank == 0)) {
            if (dt[i] == dt_min) {
                printf("!!! Warning0: dt = dt_min = %.6E \t ind = %07d \n", dt[i], ind[i]);
                fflush(stdout);
            }
        }
    } /* i */

    if (config->live_smbh_count > 0) {
        double min_dt = *std::min_element(dt, dt+N);
        for (int i=0; i<config->live_smbh_count; i++) dt[i] = min_dt;
    }

    /* Wait to all processors to finish his works... */
    MPI_Barrier(MPI_COMM_WORLD);

    /* Scatter the "local" vectors from "global" */
    MPI_Scatter(dt,  n_loc, MPI_DOUBLE, dt_loc,  n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);

    /* Wait to all processors to finish his works... */
    MPI_Barrier(MPI_COMM_WORLD);

    /* load the new values for particles to the local GRAPE's */

    nj=n_loc;

    /* load the nj particles to the G6 */

    for (int j=0; j<nj; j++) {
        a2by18 = {0, 0, 0};
        a1by6 = adot_loc[j]*(1./6.);
        aby2 = a_loc[j]*0.5;
        g6_set_j_particle(clusterid, j, ind_loc[j], t_loc[j], dt_loc[j], m_loc[j], a2by18, a1by6, aby2, v_loc[j], x_loc[j]);
    } /* j */

    if (myRank == rootRank) {

        /* Write BH data... */
        if (config->live_smbh_output) write_bh_data(time_cur, m, x, v, pot, a, adot, dt);

        /* Write BH NB data... */
        if (config->live_smbh_neighbor_output) write_bh_nb_data(time_cur, N, m, x, v);

    } /* if (myRank == rootRank) */

    /* Get the Starting time on rootRank */

    if (myRank == rootRank) {
        get_CPU_time(&CPU_time_real0, &CPU_time_user0, &CPU_time_syst0);
        get_CPU_time(&CPU_time_real,  &CPU_time_user,  &CPU_time_syst);
    } /* if (myRank == rootRank) */

    timesteps = 0.0;
    n_act_sum = 0.0;

    for (int i=1; i<N+1; i++) { // TODO why like this?
        n_act_distr[i-1] = 0.0;
    }

    g6_calls = 0.0;

#ifdef TIMING
    DT_TOT      = 0.0;

    DT_ACT_DEF1 = 0.0;
    DT_ACT_DEF2 = 0.0;
    DT_ACT_DEF3 = 0.0;
    DT_ACT_PRED = 0.0;
    DT_ACT_GRAV = 0.0;
    DT_EXT_GRAV = 0.0;
    DT_EXT_GMC_GRAV = 0.0;
    DT_GMC_GMC_GRAV = 0.0;
    DT_ACT_CORR = 0.0;
    DT_ACT_LOAD = 0.0;

    DT_STEVOL    = 0.0;
    DT_STARDISK  = 0.0;
    DT_STARDESTR = 0.0;

    DT_ACT_REDUCE = 0.0;
#endif

    /* The main integration loop */

    while (time_cur <= t_end) {

        /* Define the minimal time and the active particles on all the nodes (exclude the ZERO masses!!!) */

#ifdef TIMING
        get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif

#ifdef ACT_DEF_GRAPITE
        min_t_loc = grapite_get_minimum_time();
#else
        min_t_loc = t[0]+dt[0];

        for (int j=0; j<n_loc; j++) {
            jjj = j + myRank*n_loc;
            tmp = t[jjj] + dt[jjj];
            if (tmp < min_t_loc) min_t_loc = tmp;
        }
#endif

        /* Wait to all processors to finish his works... */
        MPI_Barrier(MPI_COMM_WORLD);

        /* Reduce the "global" min_t from min_t_loc "local" on all processors) */
        MPI_Allreduce(&min_t_loc, &min_t, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);

        /* Wait to all processors to finish his works... */
        MPI_Barrier(MPI_COMM_WORLD);

#ifdef TIMING
        get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
        DT_ACT_DEF1 += (CPU_tmp_user - CPU_tmp_user0);
#endif

#ifdef TIMING
        get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif

#ifdef ACT_DEF_GRAPITE
        int ind_act_loc[N_MAX], n_act_loc;
        grapite_active_search(min_t, ind_act_loc, &n_act_loc);
        if (myRank > 0)
            for (int i=0; i<n_act_loc; i++)
                ind_act_loc[i] += myRank*n_loc;
        int n_act_arr[256], displs[256]; // Assuming maximum of 256 processes... seems safe.
        MPI_Allgather(&n_act_loc, 1, MPI_INT, n_act_arr, 1, MPI_INT, MPI_COMM_WORLD);
        n_act = n_act_arr[0];
        for (int i=1; i<n_proc; i++)
            n_act += n_act_arr[i];
        displs[0] = 0;
        for (int i=1; i<n_proc; i++)
            displs[i]=displs[i-1]+n_act_arr[i-1];
        MPI_Allgatherv(ind_act_loc, n_act_loc, MPI_INT, ind_act, n_act_arr, displs, MPI_INT, MPI_COMM_WORLD);
#else
        n_act = 0;

        for (int i=0; i<N; i++) {
            if (t[i]+dt[i] == min_t ) {
                ind_act[n_act] = i;
                n_act++;
            }
        } /* i */
#endif      // ACT_DEF_GRAPITE

// TODO deal with it below
#ifdef ACT_DEF_GRAPITE
#error please fix here
#endif
#if defined(ACT_DEF_GRAPITE) && (defined(ADD_BH1) || defined(ADD_BH2))
#ifdef ADD_BH1
#define ACT_DEF_GRAPITE_NUMBH 1
#else
#define ACT_DEF_GRAPITE_NUMBH 2
#endif
        int acf_def_grapite_bh_count = 0;
        for (i=0; i<n_act; i++) {
            if (ind_act[i]==0) {
                i_bh1 = i;
                acf_def_grapite_bh_count++;
        }
#ifdef ADD_BH2
            else if (ind_act[i]==1) {
                i_bh2 = i;
                acf_def_grapite_bh_count++;
            }
#endif
            if (acf_def_grapite_bh_count==ACT_DEF_GRAPITE_NUMBH) break;
        }
        if (i==n_act) {
            fprintf(stderr, "ERROR: black holes were not found in the active particle list");
            exit(1);
        }
#else
        if (config->live_smbh_count > 0) {
            i_bh1 = 0;
            i_bh2 = 1;
        }
#endif

#ifdef TIMING
        get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
        DT_ACT_DEF2 += (CPU_tmp_user - CPU_tmp_user0);
#endif

#ifdef TIMING
        get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
        for (int i=0; i<n_act; i++) {
            iii = ind_act[i];

            m_act[i] = m[iii];

            x_act[i] = x[iii];
            v_act[i] = v[iii];

            t_act[i]  = t[iii];
            dt_act[i] = dt[iii];

            pot_act[i] = pot[iii];

            pot_act_ext[i] = pot_ext[iii];

            a_act[i] = a[iii];
            adot_act[i] = adot[iii];

        } /* i */

#ifdef TIMING
        get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
        DT_ACT_DEF3 += (CPU_tmp_user - CPU_tmp_user0);
#endif

        /* predict the active particles positions etc... on all the nodes */

#ifdef TIMING
        get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif

        for (int i=0; i<n_act; i++) {
            dt_tmp = min_t - t_act[i];
            dt2half = 0.5*dt_tmp*dt_tmp;
            dt3over6 = (1./3.)*dt_tmp*dt2half;
            x_act_new[i] = x_act[i] + v_act[i]*dt_tmp + a_act[i]*dt2half + adot_act[i]*dt3over6;
            v_act_new[i] = v_act[i] + a_act[i]*dt_tmp + adot_act[i]*dt2half;
        } /* i */

#ifdef TIMING
        get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
        DT_ACT_PRED += (CPU_tmp_user - CPU_tmp_user0);
#endif

        calc_self_grav(min_t, eps2, g6_calls, n_loc,
                n_act, ind_act,
                x_act_new, v_act_new,
                pot_act_tmp,
                a_act_tmp,
                adot_act_tmp,
                h2_i);

        /* Reduce the "global" vectors from "local" on all the nodes */

#ifdef TIMING
        get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif

        MPI_Allreduce(pot_act_tmp,  pot_act_new,    n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

        MPI_Allreduce(a_act_tmp,    a_act_new,    3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
        MPI_Allreduce(adot_act_tmp, adot_act_new, 3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

        /* Wait to all processors to finish his works... */
        MPI_Barrier(MPI_COMM_WORLD);

#ifdef TIMING
        get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
        DT_ACT_REDUCE += (CPU_tmp_user - CPU_tmp_user0);
#endif

        if (config->live_smbh_count == 2) {
            if (config->live_smbh_custom_eps >= 0) {
                m_bh1 = m_act[i_bh1];
                m_bh2 = m_act[i_bh2];

                x_bh1 = x_act_new[i_bh1];
                v_bh1 = v_act_new[i_bh1];

                x_bh2 = x_act_new[i_bh2];
                v_bh2 = v_act_new[i_bh2];

                // calculate and "minus" the BH <-> BH softened pot, acc & jerk

                tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1,
                                        m_bh2, x_bh2, v_bh2,
                                        eps,
                                        &pot_bh1, a_bh1, adot_bh1,
                                        &pot_bh2, a_bh2, adot_bh2);

                pot_act_new[i_bh1] -= pot_bh1;
                pot_act_new[i_bh2] -= pot_bh2;

                a_act_new[i_bh1] -= a_bh1;
                a_act_new[i_bh2] -= a_bh2;

                adot_act_new[i_bh1] -= adot_bh1;
                adot_act_new[i_bh2] -= adot_bh2;

                // calculate and "plus" the new BH <-> BH unsoftened pot, acc, jerk

                tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1,
                                        m_bh2, x_bh2, v_bh2,
                                        config->live_smbh_custom_eps,
                                        &pot_bh1, a_bh1, adot_bh1,
                                        &pot_bh2, a_bh2, adot_bh2);

                pot_act_new[i_bh1] += pot_bh1;
                pot_act_new[i_bh2] += pot_bh2;

                a_act_new[i_bh1] += a_bh1;
                a_act_new[i_bh2] += a_bh2;

                adot_act_new[i_bh1] += adot_bh1;
                adot_act_new[i_bh2] += adot_bh2;
            }

            if (config->binary_smbh_pn) {
                // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk

                dt_bh_tmp = dt[0];

                tmp_i = calc_force_pn_BH(m_bh1, x_bh1, v_bh1, s_bh1,
                                        m_bh2, x_bh2, v_bh2, s_bh2,
                                        C_NB,  dt_bh_tmp, usedOrNot,
                                        (double(*)[3])a_pn1, (double(*)[3])adot_pn1,
                                        (double(*)[3])a_pn2, (double(*)[3])adot_pn2, myRank, rootRank);

                a_act_new[i_bh1] += a_pn1[1] + a_pn1[2] + a_pn1[3] + a_pn1[4] + a_pn1[5] + a_pn1[6];
                a_act_new[i_bh2] += a_pn2[1] + a_pn2[2] + a_pn2[3] + a_pn2[4] + a_pn2[5] + a_pn2[6];

                adot_act_new[i_bh1] += adot_pn1[1] + adot_pn1[2] + adot_pn1[3] + adot_pn1[4] + adot_pn1[5] + adot_pn1[6];
                adot_act_new[i_bh2] += adot_pn2[1] + adot_pn2[2] + adot_pn2[3] + adot_pn2[4] + adot_pn2[5] + adot_pn2[6];

                if (myRank == rootRank) {
                    if (tmp_i == 505) {
                        printf("PN RSDIST: TS = %.8E \t t  = %.8E \n", timesteps, time_cur);
                        fflush(stdout);
                        exit(-1);
                    }
                }
            }
        }

        calc_ext_grav(external_gravity_components, n_act, x_act_new, v_act_new, pot_act_ext, a_act_new, adot_act_new);

        /* correct the active particles positions etc... on all the nodes */

#ifdef TIMING
        get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif

        for (int i=0; i<n_act; i++) {
            dt_tmp = min_t - t_act[i];

            dt3over6 = dt_tmp*dt_tmp*dt_tmp/6.0;
            dt4over24 = dt3over6*dt_tmp/4.0;
            dt5over120 = dt4over24*dt_tmp/5.0;

            dtinv = 1.0/dt_tmp;
            dt2inv = dtinv*dtinv;
            dt3inv = dt2inv*dtinv;

            double3 a0mia1 = a_act[i]-a_act_new[i];
            double3 ad04plad12 = 4.0*adot_act[i] + 2.0*adot_act_new[i];
            double3 ad0plad1 = adot_act[i] + adot_act_new[i];

            a2 = -6.0*a0mia1*dt2inv - ad04plad12*dtinv;
            a3 = 12.0*a0mia1*dt3inv + 6.0*ad0plad1*dt2inv;

            x_act_new[i] += dt4over24*a2 + dt5over120*a3;
            v_act_new[i] += dt3over6*a2 + dt4over24*a3;

            a1abs = sqrt(a_act_new[i][0]*a_act_new[i][0]+a_act_new[i][1]*a_act_new[i][1]+a_act_new[i][2]*a_act_new[i][2]);
            adot1abs = sqrt(adot_act_new[i][0]*adot_act_new[i][0]+adot_act_new[i][1]*adot_act_new[i][1]+adot_act_new[i][2]*adot_act_new[i][2]);

            a2dot1 = a2 + dt_tmp*a3;

            a2dot1abs = a2dot1.norm();
            a3dot1abs = a3.norm();

            if ((config->live_smbh_count > 0) && (ind_act[i] < config->live_smbh_count))
                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));

            if (dt_new < dt_min) dt_tmp = 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);
            }

            if ((dt_new > 2.0*dt_tmp) && (fmod(min_t, 2.0*dt_tmp) == 0.0) && (2.0*dt_tmp <= dt_max)) {
                dt_tmp *= 2.0;
            }

            dt_act[i] = dt_tmp;

            t_act[i] = min_t;

            pot_act[i] = pot_act_new[i];

            x_act[i]    =    x_act_new[i];
            v_act[i]    =    v_act_new[i];
            a_act[i]    =    a_act_new[i];
            adot_act[i] = adot_act_new[i];

            if (config->dt_min_warning && (myRank == 0)) {
                if (dt_act[i] == dt_min) {
                    printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d \n", dt[i], ind_act[i]);
                    fflush(stdout);
                }
            }

        } /* i */

        /* define the min. dt over all the act. part. and set it also for the BH... */

        if (config->live_smbh_count > 0) {
            double min_dt = *std::min_element(dt_act, dt_act+n_act);
            if (config->live_smbh_count>=1) dt_act[i_bh1] = min_dt;
            if (config->live_smbh_count==2) dt_act[i_bh2] = min_dt;
        }

        if (config->binary_smbh_influence_sphere_output) {
            if (myRank == rootRank) {

                out = fopen("bbh_inf.dat", "a");

                m_bh1 = m_act[i_bh1];
                m_bh2 = m_act[i_bh2];

                x_bh1 = x_act[i_bh1];
                x_bh2 = x_act[i_bh2];

                v_bh1 = v_act[i_bh1];
                v_bh2 = v_act[i_bh2];

                x_bbhc = (m_bh1*x_bh1 + m_bh2*x_bh2)/(m_bh1 + m_bh2);
                v_bbhc = (m_bh1*v_bh1 + m_bh2*v_bh2)/(m_bh1 + m_bh2);

                DR2 = SQR(x_bh1[0] - x_bh2[0]) + SQR(x_bh1[1] - x_bh2[1]) + SQR(x_bh1[2] - x_bh2[2]);
                DV2 = SQR(v_bh1[0] - v_bh2[0]) + SQR(v_bh1[1] - v_bh2[1]) + SQR(v_bh1[2] - v_bh2[2]);

                EB = -(m_bh1 + m_bh2) / sqrt(DR2) + 0.5 * DV2;

                SEMI_a = -0.5 * (m_bh1 + m_bh2)/EB;
                SEMI_a2 = SQR(SEMI_a);

                for (int i=2; i<n_act; i++) {

                    tmp_r2 = SQR(x_act[i][0] - x_bbhc[0]) + SQR(x_act[i][1] - x_bbhc[1]) + SQR(x_act[i][2] - x_bbhc[2]);

                    iii = ind_act[i];

                    if (tmp_r2 < SEMI_a2*SQR(config->binary_smbh_influence_radius_factor)) {

                        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],
                                sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2],
                                m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_act[0],
                                m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_act[1],
                                sqrt(tmp_r2),
                                m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i],
                                dt_act[i]);

                        inf_event[iii] = 1;
                        }

                    } else {
                        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],
                                m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_act[0],
                                m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_act[1],
                                sqrt(tmp_r2),
                                m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i],
                                dt_act[i]);
                        }

                        inf_event[iii] = 0;

                    } /* if (tmp_r2 < DR2*R_INF2) */

                } /* i */

                fclose(out);

            } /* if (myRank == rootRank) */
        }

        /* Return back the new coordinates + etc... of active part. to the global data... */

        for (int i=0; i<n_act; i++) {

            iii = ind_act[i];

            m[iii] = m_act[i];

            x[iii] = x_act[i];
            v[iii] = v_act[i];

            t[iii] = t_act[i];
            dt[iii] = dt_act[i];

            pot[iii] = pot_act[i];

            pot_ext[iii] = pot_act_ext[i];

            a[iii] = a_act[i];
            adot[iii] = adot_act[i];
        } /* i */

#ifdef TIMING
        get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
        DT_ACT_CORR += (CPU_tmp_user - CPU_tmp_user0);
#endif

        /* load the new values for active particles to the local GRAPE's */

#ifdef TIMING
        get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif

        for (int j=0; j<n_act; j++) { // TODO would be nicer to use i instead of j here
#ifdef ETICS_CEP
            if (ind_act[j] == grapite_cep_index) grapite_update_cep(t_act[j], x_act[j], v_act[j], a_act[j], adot_act[j]); // All ranks should do it.
#endif
            cur_rank = ind_act[j]/n_loc;

            if (myRank == cur_rank) {

                jjj = ind_act[j] - myRank*n_loc;

                a2by18 = {0, 0, 0};
                a1by6 = adot_act[j]*(1./6.);
                aby2 = a_act[j]*0.5;

                g6_set_j_particle(clusterid, jjj, ind_act[j], t_act[j], dt_act[j], m_act[j], a2by18, a1by6, aby2, v_act[j], x_act[j]);

            } /* if (myRank == cur_rank) */

        } /* j */

#ifdef TIMING
        get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
        DT_ACT_LOAD += (CPU_tmp_user - CPU_tmp_user0);
#endif

        /* Current time set to min_t */

        time_cur = min_t;

        timesteps += 1.0;
        n_act_sum += n_act;
        n_act_distr[n_act-1]++;

        if (time_cur >= t_bh) {
            if (myRank == rootRank) {
                /* Write BH data... */
                if (config->live_smbh_output) write_bh_data(time_cur, m, x, v, pot, a, adot, dt);

                /* Write BH NB data... */
                if (config->live_smbh_neighbor_output) write_bh_nb_data(time_cur, N, m, x, v);

            } /* if (myRank == rootRank) */

            t_bh += dt_bh;
        } /* if (time_cur >= t_bh) */

        if (time_cur >= t_contr) {
            if (myRank == rootRank) {

                energy_contr(time_cur, timesteps, n_act_sum, g6_calls, rcm_sum, vcm_sum, E_tot_0, E_tot_corr_0, E_tot_corr_sd_0, skip_con, N, m, x, v, pot, pot_ext);

                /* write cont data */

                write_snap_data((char*)"data.con", diskstep, N, time_cur, ind, m, x, v);

                /* possible OUT for timing !!! */

#ifdef TIMING
                out = fopen("timing.dat", "a");

                DT_TOT = DT_ACT_DEF1 + DT_ACT_DEF2 + DT_ACT_DEF3 + DT_ACT_PRED +
                         DT_ACT_GRAV + DT_EXT_GRAV + DT_GMC_GRAV +
                         DT_GMC_GMC_GRAV + DT_EXT_GMC_GRAV +
                         DT_ACT_CORR + DT_ACT_LOAD +
                         DT_STEVOL + DT_STARDISK + DT_STARDESTR +
                         DT_ACT_REDUCE;

                fprintf(out,"%.8E \t %.6E \t %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f \t %.3f \t %.8E %.8E %.8E \t %.8E %.8E %.8E \n",
                        time_cur, DT_TOT,
                        100.0*DT_ACT_DEF1/DT_TOT, 100.0*DT_ACT_DEF2/DT_TOT, 100.0*DT_ACT_DEF3/DT_TOT, 100.0*DT_ACT_PRED/DT_TOT,
                        100.0*DT_ACT_GRAV/DT_TOT, 100.0*DT_EXT_GRAV/DT_TOT, 100.0*DT_GMC_GRAV/DT_TOT,
                        100.0*DT_GMC_GMC_GRAV/DT_TOT, 100.0*DT_EXT_GMC_GRAV/DT_TOT,
                        100.0*DT_ACT_CORR/DT_TOT, 100.0*DT_ACT_LOAD/DT_TOT,
                        100.0*DT_STEVOL/DT_TOT, 100.0*DT_STARDISK/DT_TOT, 100.0*DT_STARDESTR/DT_TOT,
                        100.0*DT_ACT_REDUCE/DT_TOT,
                        CPU_time_real-CPU_time_real0, CPU_time_user-CPU_time_user0, CPU_time_syst-CPU_time_syst0,
                        timesteps, n_act_sum, 57.0*N*n_act_sum/(CPU_time_user-CPU_time_user0)/1.0E+09);

                fclose(out);
#endif

            } /* if (myRank == rootRank) */

#ifdef ETICS_CEP
            // We are /inside/ a control step, so all particles must be synchronized; we can safely calculate their density centre. The acceleration and jerk currently in the memory are for the predicted position of the CEP, by calling grapite_calc_center we "correct" the position and velocity, but not the gravity at that point.

            // First calculate the DC
            grapite_calc_center(N, m, x, v, xcm, vcm, xdc, vdc);

            // Now copy it to the global particle list
            memcpy(x[grapite_cep_index], xdc, 3*sizeof(double));
            memcpy(v[grapite_cep_index], vdc, 3*sizeof(double));

            // Now copy it to the local particle list for tha appropriate rank
            if (myRank == grapite_cep_index/n_loc) {
                memcpy(x_loc[grapite_cep_index - myRank*n_loc], xdc, 3*sizeof(double));
                memcpy(v_loc[grapite_cep_index - myRank*n_loc], vdc, 3*sizeof(double));
            }
            grapite_update_cep(time_cur, xdc, vdc, a[grapite_cep_index], adot[grapite_cep_index]);
#endif

            t_contr += dt_contr;
        } /* if (time_cur >= t_contr) */

        if (time_cur >= t_disk) {
            if (myRank == rootRank) {
                diskstep++;
                char out_fname[256];
                sprintf(out_fname, "%06d.dat", diskstep);
                write_snap_data(out_fname, diskstep, N, time_cur, ind, m, x, v);
            } /* if (myRank == rootRank) */

#ifdef ETICS_DUMP
            sprintf(out_fname, "coeffs.%06d.%02d.dat", diskstep, myRank);
            grapite_dump(out_fname, 2);
#endif

            t_disk += dt_disk;
        } /* if (time_cur >= t_disk) */

    } /* while (time_cur < t_end) */

    /* close the local GRAPEs */

    g6_close(clusterid);

    /* Wait to all processors to finish his works... */
    MPI_Barrier(MPI_COMM_WORLD);

    MPI_Reduce(&g6_calls, &g6_calls_sum, 1, MPI_DOUBLE, MPI_SUM, rootRank, MPI_COMM_WORLD);

    /* Wait to all processors to finish his works... */
    MPI_Barrier(MPI_COMM_WORLD);

    if (myRank == rootRank) {

        /* Write some output for the timestep annalize... */

        printf("\n");
        printf("timesteps = %.0f   Total sum of integrated part. = %.0f   g6_calls on all nodes = %.0f \n", timesteps, n_act_sum, g6_calls);
        printf("\n");
        printf("Real Speed = %.3f GFlops \n", 57.0*N*n_act_sum/(CPU_time_user-CPU_time_user0)/1.0E+09);
        fflush(stdout);

    } /* if (myRank == rootRank) */

    /* Wait to all processors to finish his works... */
    MPI_Barrier(MPI_COMM_WORLD);

    /* Finalize the MPI work */
    MPI_Finalize();

    return 0;
}