 
         /*** in the following 8 lines, we calculate indexes of all the children that have to be walked from the index of the first child***/
         if (flag && nodes[prefix[tid] + 1] == 0) nodes[prefix[tid] + 1] = child + 1; __syncthreads();
-	if (flag && nodes[prefix[tid] + 2] == 0) nodes[prefix[tid] + 2] = child + 2; __syncthreads();
-	if (flag && nodes[prefix[tid] + 3] == 0) nodes[prefix[tid] + 3] = child + 3; __syncthreads();
-	if (flag && nodes[prefix[tid] + 4] == 0) nodes[prefix[tid] + 4] = child + 4; __syncthreads();
-	if (flag && nodes[prefix[tid] + 5] == 0) nodes[prefix[tid] + 5] = child + 5; __syncthreads();
-	if (flag && nodes[prefix[tid] + 6] == 0) nodes[prefix[tid] + 6] = child + 6; __syncthreads();
-	if (flag && nodes[prefix[tid] + 7] == 0) nodes[prefix[tid] + 7] = child + 7; __syncthreads();
-
-	n_offset += n_total;    //Increase the offset in the array by the number of newly added nodes
-
-
-	/***
-	**** --> save list of nodes to LMEM
-	***/
+        if (flag && nodes[prefix[tid] + 2] == 0) nodes[prefix[tid] + 2] = child + 2; __syncthreads();
+        if (flag && nodes[prefix[tid] + 3] == 0) nodes[prefix[tid] + 3] = child + 3; __syncthreads();
+        if (flag && nodes[prefix[tid] + 4] == 0) nodes[prefix[tid] + 4] = child + 4; __syncthreads();
+        if (flag && nodes[prefix[tid] + 5] == 0) nodes[prefix[tid] + 5] = child + 5; __syncthreads();
+        if (flag && nodes[prefix[tid] + 6] == 0) nodes[prefix[tid] + 6] = child + 6; __syncthreads();
+        if (flag && nodes[prefix[tid] + 7] == 0) nodes[prefix[tid] + 7] = child + 7; __syncthreads();
+
+        n_offset += n_total;    //Increase the offset in the array by the number of newly added nodes
+
+
+        /***
+         **** --> save list of nodes to LMEM
+         ***/
 
         /*** if half of shared memory or more is filled with the the nodes, dump these into slowmem stack ***/
-	while(n_offset >= DIM) {
-	  n_offset -= DIM;
-	  const int offs1 = ACCS<SHIFT>(n_stack1);
-	  nstack[offs1] = nodes[n_offset + tid];   n_stack1++;
-	  n_nodes1 += DIM;
+        while(n_offset >= DIM) {
+          n_offset -= DIM;
+          const int offs1 = ACCS<SHIFT>(n_stack1);
+          nstack[offs1] = nodes[n_offset + tid];   n_stack1++;
+          n_nodes1 += DIM;
 
           if((n_stack1 - c_stack0) >= (LMEM_STACK_SIZE << SHIFT))
           {
             //We overwrote our current stack
-	    apprCount = -1; return acc_i;	 
-          }
-	}
-
-	__syncthreads();
-
-
-
-	/******************************/
-	/******************************/
-	/*****     EVALUATION     *****/
-	/******************************/
-	/******************************/
-#if 1
-	/***********************************/
-	/******       APPROX          ******/
-	/***********************************/
-
-        #ifdef OLDPREFIX
-          n_total = calc_prefix<DIM2>(prefix, tid,  1 - (split || !use_node));
-        #else
-          inclusive_scan_block<ADDOP<int>, int>(prefix, 1 - (split || !use_node), tid);
-          n_total = prefix[blockDim.x - 1];
-        #endif
-
-
-// 	n_total = calc_prefix<DIM2>(prefix, tid,  !split && use_node);         // for some unkown reason this does not work right on the GPU
-	if (!split && use_node) approx[n_approx + prefix[tid] - 1] = node;
- 	__syncthreads();
-	n_approx += n_total;
-
- 	while (n_approx >= DIM) {
- 	  n_approx -= DIM;
-	  int address      = (approx[n_approx + tid] << 1) + approx[n_approx + tid];
-          #ifndef TEXTURES
-            float4 monopole  = multipole_data[address    ];
-            float4 octopole0 = multipole_data[address + 1];
-            float4 octopole1 = multipole_data[address + 2];
-          #else
-            float4 monopole  = tex1Dfetch(texMultipole, address);
-            float4 octopole0 = tex1Dfetch(texMultipole, address + 1);
-            float4 octopole1 = tex1Dfetch(texMultipole, address + 2);
-          #endif
-
- 	  node_mon0[tid] = monopole.w;
- 	  node_mon1[tid] = make_float3(monopole.x,  monopole.y,  monopole.z);
- 	  node_oct0[tid] = make_float3(octopole0.x, octopole0.y, octopole0.z);
- 	  node_oct1[tid] = make_float3(octopole1.x, octopole1.y, octopole1.z);
-
-#if 0
-    const float f_dm   = 0.0f;
-    const float f_star = 1.0f
-    const float darkMatterMass = f_dm   * octopole1.w;
-    /* eg: we need to be careful with the line below to avoid truncation error due to 
-       subtraction of two large numbers, monopole.w and darkMatterMass both could be
-       very large.
-       Instead, we can use octopole1.w to be stellar mass, and DM mass to be 
-       monopole.w, then we add the two together to get total mass, but this will
-       require more changes to the kernel */
-    const float    stellarMass = f_star * (monopole.w - darkMatterMass);
-    const float hinv = 1.0f/hi;   /* eg: this can be precomputing to avoid division */
-    density += interact(
-        make_float3(pos_i.x, pos_i.y, pos_i.z), h, hinv,
-        make_float3(monopole.x, monople.y, monopole.z), darkMatterMass + stellarMass);
-    /* eg: the interact function still calls sqrtf(f), which invloves 1 div and 1 rsqrtf,
-       so ideally we would like to take advantage of rsqrtf in add_acc, and then we only
-       do 1 div */
-#endif
-
-#ifdef INDSOFT
-    float temp          = node_eps[tid]; //Backup value in the shmem into register
-    node_eps[tid]       = octopole0.w;
-#endif
-
-    __syncthreads();
-#pragma unroll
-    for (int i = 0; i < DIMx; i++)
-    {
-      apprCount++;
-#ifdef DOGRAV
-#ifdef INDSOFT
-      acc_i = add_acc(acc_i, pos_i,
-          node_mon0[offs + i], node_mon1[offs + i],
-          node_oct0[offs + i], node_oct1[offs + i], node_eps[offs + i], eps2);
-
-#else
-      acc_i = add_acc(acc_i, pos_i,
-          node_mon0[offs + i], node_mon1[offs + i],
-          node_oct0[offs + i], node_oct1[offs + i], eps2);
-#endif
-
-#endif
-    }
-    __syncthreads();
-#ifdef INDSOFT
-    node_eps[tid] = temp; //Restore original value in shmem
-    __syncthreads();
-#endif
-  }
-  __syncthreads();
-#endif
-
-#if 1
-  /***********************************/
-  /******       DIRECT          ******/
-  /***********************************/
-
-  int *sh_body = &approx[DIM];
-
-  flag         = split && leaf && use_node;                                //flag = split + leaf + use_node
-  int  jbody   = node_data & BODYMASK;                                     //the first body in the leaf
-  int  nbody   = (((node_data & INVBMASK) >> LEAFBIT)+1) & BTEST(flag);    //number of bodies in the leaf masked with the flag
-
-  body_list[tid] = direct[tid];                                            //copy list of bodies from previous pass to body_list
-  sh_body  [tid] = jbody;                                                  //store the leafs first body id into shared memory
-
-  // step 1
-#ifdef OLDPREFIX
-  calc_prefix<DIM2>(prefix0, tid, flag);
-#else
-  inclusive_scan_block<ADDOP<int>, int>(prefix0, (int)flag, tid);       // inclusive scan on flags to construct array
-#endif
-
-  if (flag) prefix1[prefix0[tid] - 1] = tid;                             //with tidś whose leaves have to be opened
-  __syncthreads();                                                      //thread barrier, make sure all warps completed the job
-
-  // step 2
-#ifdef OLDPREFIX
-  int n_bodies  = calc_prefix<DIM2>(prefix0, tid, nbody);
-#else
-  inclusive_scan_block<ADDOP<int>, int>(prefix0, nbody, tid);        // inclusive scan to compute memory offset for each body
-  int n_bodies = prefix0[blockDim.x - 1];                            //Total number of bides extract from the leaves
-  __syncthreads();                                                   // thread barrier to make sure that warps completed their jobs
-#endif
-
-  direct [tid]  = prefix0[tid];                                       //Store a copy of inclusive scan in direct
-  prefix0[tid] -= nbody;                                              //convert inclusive int oexclusive scan
-  prefix0[tid] += 1;                                                  //add unity, since later prefix0[tid] == 0 used to check barrier
-
-  int nl_pre = 0;                                                     //Number of leaves that have already been processed
-
-#define NJMAX (DIM*4)
-  while (n_bodies > 0) {
-    int nb    = min(n_bodies, NJMAX - n_direct);                    //Make sure number of bides to be extracted does not exceed
-    //the amount of allocated shared memory
-
-    // step 0                                                      //nullify part of the body_list that will be filled with bodies
-    for (int i = n_direct; i < n_direct + nb; i += DIM){           //from the leaves that are being processed
-      body_list[i + tid] = 0;
-    }
-    __syncthreads();
-
-    //step 1:
-    if (flag && (direct[tid] <= nb) && (prefix0[tid] > 0))        //make sure that the thread indeed carries a leaf
-      body_list[n_direct + prefix0[tid] - 1] = 1;                 //whose bodies will be extracted
-    __syncthreads();
-
-    //step 2:
-#ifdef OLDPREFIX
-    int nl = calc_prefix<DIM2>(nb, &body_list[n_direct], tid);
-#else
-    int nl = inclusive_scan_array<ADDOP<int>, int>              // inclusive scan to compute number of leaves to process
-      (&body_list[n_direct], nb, tid);            // to make sure that there is enough shared memory for bodies
-#endif
-    nb = direct[prefix1[nl_pre + nl - 1]];                        // number of bodies stored in these leaves
-
-    // step 3:
-    for (int i = n_direct; i < n_direct + nb; i += DIM) {          //segmented fill of the body_list
-      int j = prefix1[nl_pre + body_list[i + tid] - 1];            // compute the first body in shared j-body array
-      body_list[i + tid] = (i + tid - n_direct) -                 //add to the index of the first j-body in a child
-        (prefix0[j] - 1) + sh_body[j];         //the index of the first child in body_list array
-    }
-    __syncthreads();
-
-
-    /**************************************************
-     *  example of what is accomplished in steps 0-4   *
-     *       ---------------------------               *
-     * step 0: body_list = 000000000000000000000       *
-     * step 1: body_list = 100010001000000100100       *
-     * step 2: body_list = 111122223333333444555       *
-     * step 3: body_list = 012301230123456012012       *
-     *         assuming that sh_body[j] = 0            *
-     ***************************************************/
-
-    n_bodies     -= nb;                                   //subtract from n_bodies number of bodies that have been extracted
-    nl_pre       += nl;                                   //increase the number of leaves that where processed
-    direct [tid] -= nb;                                   //subtract the number of extracted bodies in this pass
-    prefix0[tid] = max(prefix0[tid] - nb, 0);             //same here, but do not let the number be negative (GT200 bug!?)
-    n_direct     += nb;                                  //increase the number of bodies to be procssed
-
-    while(n_direct >= DIM) {
-      n_direct -= DIM;
-
-
-      float4 posj  = body_pos[body_list[n_direct + tid]];
-      //             float4 posj  = tex1Dfetch(texBody, body_list[n_direct + tid]);
-      sh_mass[tid] = posj.w;
-      sh_pos [tid] = make_float3(posj.x, posj.y, posj.z);
-      sh_jid [tid] = body_list[n_direct + tid];
-
-#ifdef INDSOFT
-      float   temp = sh_eps [tid];        //Store the value from Shmem into a register
-      sh_eps [tid] = body_vel[body_list[n_direct + tid]].w;  //Load the softening
-#endif
-      __syncthreads();
-#pragma unroll
-      for (int j = 0; j < DIMx; j++)
-      {
-        int selfGrav = (body_i != sh_jid[offs + j]);
-        // 	      if (body_i != sh_jid[offs + j]) //If statement replaced by multiplication
-        {
-          float ds2 =  1.0e10f;
-          direCount++;
-#ifdef DOGRAV
-#ifdef INDSOFT
-          acc_i = add_acc(acc_i, pos_i, sh_mass[offs + j], sh_pos[offs + j], sh_eps[offs + j], ds2, eps2, selfGrav);
-#else
-          acc_i = add_acc(acc_i, pos_i, sh_mass[offs + j], sh_pos[offs + j], ds2, eps2, selfGrav);
-#endif
-#endif
-          ds2 += selfGrav*1.0e10f;
-          if (ds2 < ds2_min) {
-            ngb     = sh_jid[offs + j];
-            ds2_min = ds2;
+            apprCount = -1; return acc_i;	 
           }
         }
-      }
-#ifdef INDSOFT
-      //Restore the shmem value after all threads have used it
-      __syncthreads();
-      sh_eps[tid] = temp;
-#endif
-      __syncthreads();
-    }
-
-  }
-  direct[tid] = body_list[tid];
-  __syncthreads();
+
+        __syncthreads();
+
+
+
+        /******************************/
+        /******************************/
+        /*****     EVALUATION     *****/
+        /******************************/
+        /******************************/
+#if 1
+        /***********************************/
+        /******       APPROX          ******/
+        /***********************************/
+
+#ifdef OLDPREFIX
+        n_total = calc_prefix<DIM2>(prefix, tid,  1 - (split || !use_node));
+#else
+        inclusive_scan_block<ADDOP<int>, int>(prefix, 1 - (split || !use_node), tid);
+        n_total = prefix[blockDim.x - 1];
+#endif
+
+
+        // 	n_total = calc_prefix<DIM2>(prefix, tid,  !split && use_node);         // for some unkown reason this does not work right on the GPU
+        if (!split && use_node) approx[n_approx + prefix[tid] - 1] = node;
+        __syncthreads();
+        n_approx += n_total;
+
+        while (n_approx >= DIM) 
+        {
+          n_approx -= DIM;
+          const int address      = (approx[n_approx + tid] << 1) + approx[n_approx + tid];
+#ifndef TEXTURES
+          const float4 monopole  = multipole_data[address    ];
+#if 0
+          float4 octopole0 = multipole_data[address + 1];
+          float4 octopole1 = multipole_data[address + 2];
+#endif
+#else
+          const float4 monopole  = tex1Dfetch(texMultipole, address);
+#if 0
+          float4 octopole0 = tex1Dfetch(texMultipole, address + 1);
+          float4 octopole1 = tex1Dfetch(texMultipole, address + 2);
+#endif
+#endif
+
+          node_mon0[tid] = monopole.w;
+          node_mon1[tid] = make_float3(monopole.x,  monopole.y,  monopole.z);
+          __syncthreads();
+
+#if 0
+          const float f_dm   = 0.0f;
+          const float f_star = 1.0f
+            const float darkMatterMass = f_dm   * octopole1.w;
+          /* eg: we need to be careful with the line below to avoid truncation error due to 
+             subtraction of two large numbers, monopole.w and darkMatterMass both could be
+             very large.
+             Instead, we can use octopole1.w to be stellar mass, and DM mass to be 
+             monopole.w, then we add the two together to get total mass, but this will
+             require more changes to the kernel */
+          const float    stellarMass = f_star * (monopole.w - darkMatterMass);
+          const float hinv = 1.0f/hi;   /* eg: this can be precomputing to avoid division */
+          density += interact(
+              make_float3(pos_i.x, pos_i.y, pos_i.z), h, hinv,
+              make_float3(monopole.x, monople.y, monopole.z), darkMatterMass + stellarMass);
+          /* eg: the interact function still calls sqrtf(f), which invloves 1 div and 1 rsqrtf,
+             so ideally we would like to take advantage of rsqrtf in add_acc, and then we only
+             do 1 div */
+#endif
+
+
+#pragma unroll
+          for (int i = 0; i < DIMx; i++)
+            acc_i = add_acc(acc_i, pos_i, node_mon0[offs + i], node_mon1[offs+i], eps2);
+          apprCount += DIMx;
+          __syncthreads();
+        }
+        __syncthreads();
+#endif
+
+#if 1
+        /***********************************/
+        /******       DIRECT          ******/
+        /***********************************/
+
+        int *sh_body = &approx[DIM];
+
+        flag         = split && leaf && use_node;                                //flag = split + leaf + use_node
+        int  jbody   = node_data & BODYMASK;                                     //the first body in the leaf
+        int  nbody   = (((node_data & INVBMASK) >> LEAFBIT)+1) & BTEST(flag);    //number of bodies in the leaf masked with the flag
+
+        body_list[tid] = direct[tid];                                            //copy list of bodies from previous pass to body_list
+        sh_body  [tid] = jbody;                                                  //store the leafs first body id into shared memory
+
+        // step 1
+#ifdef OLDPREFIX
+        calc_prefix<DIM2>(prefix0, tid, flag);
+#else
+        inclusive_scan_block<ADDOP<int>, int>(prefix0, (int)flag, tid);       // inclusive scan on flags to construct array
+#endif
+
+        if (flag) prefix1[prefix0[tid] - 1] = tid;                             //with tidś whose leaves have to be opened
+        __syncthreads();                                                      //thread barrier, make sure all warps completed the job
+
+        // step 2
+#ifdef OLDPREFIX
+        int n_bodies  = calc_prefix<DIM2>(prefix0, tid, nbody);
+#else
+        inclusive_scan_block<ADDOP<int>, int>(prefix0, nbody, tid);        // inclusive scan to compute memory offset for each body
+        int n_bodies = prefix0[blockDim.x - 1];                            //Total number of bides extract from the leaves
+        __syncthreads();                                                   // thread barrier to make sure that warps completed their jobs
+#endif
+
+        direct [tid]  = prefix0[tid];                                       //Store a copy of inclusive scan in direct
+        prefix0[tid] -= nbody;                                              //convert inclusive int oexclusive scan
+        prefix0[tid] += 1;                                                  //add unity, since later prefix0[tid] == 0 used to check barrier
+
+        int nl_pre = 0;                                                     //Number of leaves that have already been processed
+
+#define NJMAX (DIM*4)
+        while (n_bodies > 0) 
+        {
+          int nb    = min(n_bodies, NJMAX - n_direct);                    //Make sure number of bides to be extracted does not exceed
+          //the amount of allocated shared memory
+
+          // step 0                                                      //nullify part of the body_list that will be filled with bodies
+          for (int i = n_direct; i < n_direct + nb; i += DIM){           //from the leaves that are being processed
+            body_list[i + tid] = 0;
+          }
+          __syncthreads();
+
+          //step 1:
+          if (flag && (direct[tid] <= nb) && (prefix0[tid] > 0))        //make sure that the thread indeed carries a leaf
+            body_list[n_direct + prefix0[tid] - 1] = 1;                 //whose bodies will be extracted
+          __syncthreads();
+
+          //step 2:
+#ifdef OLDPREFIX
+          int nl = calc_prefix<DIM2>(nb, &body_list[n_direct], tid);
+#else
+          int nl = inclusive_scan_array<ADDOP<int>, int>              // inclusive scan to compute number of leaves to process
+            (&body_list[n_direct], nb, tid);            // to make sure that there is enough shared memory for bodies
+#endif
+          nb = direct[prefix1[nl_pre + nl - 1]];                        // number of bodies stored in these leaves
+
+          // step 3:
+          for (int i = n_direct; i < n_direct + nb; i += DIM) {          //segmented fill of the body_list
+            int j = prefix1[nl_pre + body_list[i + tid] - 1];            // compute the first body in shared j-body array
+            body_list[i + tid] = (i + tid - n_direct) -                 //add to the index of the first j-body in a child
+              (prefix0[j] - 1) + sh_body[j];         //the index of the first child in body_list array
+          }
+          __syncthreads();
+
+
+          /**************************************************
+           *  example of what is accomplished in steps 0-4   *
+           *       ---------------------------               *
+           * step 0: body_list = 000000000000000000000       *
+           * step 1: body_list = 100010001000000100100       *
+           * step 2: body_list = 111122223333333444555       *
+           * step 3: body_list = 012301230123456012012       *
+           *         assuming that sh_body[j] = 0            *
+           ***************************************************/
+
+          n_bodies     -= nb;                                   //subtract from n_bodies number of bodies that have been extracted
+          nl_pre       += nl;                                   //increase the number of leaves that where processed
+          direct [tid] -= nb;                                   //subtract the number of extracted bodies in this pass
+          prefix0[tid] = max(prefix0[tid] - nb, 0);             //same here, but do not let the number be negative (GT200 bug!?)
+          n_direct     += nb;                                  //increase the number of bodies to be procssed
+
+          while(n_direct >= DIM) 
+          {
+            n_direct -= DIM;
+
+
+            const float4 posj  = body_pos[body_list[n_direct + tid]];
+#if 0
+            const float4 posj  = tex1Dfetch(texBody, body_list[n_direct + tid]);
+#endif
+            sh_mass[tid] = posj.w;
+            sh_pos [tid] = make_float3(posj.x, posj.y, posj.z);
+            sh_jid [tid] = body_list[n_direct + tid];  /* we need this to distinghuis between DM and *-particles */
+
+            __syncthreads();
+#pragma unroll
+            for (int j = 0; j < DIMx; j++)
+              acc_i = add_acc(acc_i, pos_i, sh_mass[offs + j], sh_pos[offs + j], eps2);
+#if 0
+            direCount += DIMx;
+#endif
+            __syncthreads();
+          }
+
+        }
+        direct[tid] = body_list[tid];
+        __syncthreads();
 #endif
       } //end lvl
 


}

#ifdef INDSOFT
__device__ float4 get_D04(float ds2, float epsP, float epsQ, int selfGrav = 1) {

//   float eps     = fmaxf(epsP, epsQ);
  float epseff  = epsP + epsQ;
  float ids, ids2, ids3;

  if(ds2 >= (epseff*epseff))
  {
     ids  = rsqrtf(ds2);
     //if(isnan(ids)) ids = 0; not needed if we use non-zero softening
     ids  = ids*selfGrav;       //Prevent selfGravity, instead of using if-statement
     ids3 = ids*ids*ids;
  }
  else
  {
    //these two lines are faster than a real sqrt
    float dist = ds2*selfGrav*rsqrtf(ds2); //Gives NaN is ds is 0
    if(isnan(dist)) dist = 0.0f;
    //float    dist = sqrtf(ds2); //Slower than the two lines above
    //assert(!isnan(dist));

    float epseffi = 1.f/epseff;
    float rhinv   = dist*epseffi;

    if(rhinv < 0.5f)
    {
      ids3  = 4.f/3.f + (rhinv*rhinv)*(4.f*rhinv-4.8f);
      ids   = 1.4f - (rhinv*rhinv)*(8.f/3.f+(rhinv*rhinv)*(3.2f*rhinv-4.8f));
    }
    else
    {
      ids3 = 8.f/3.f-6.f*rhinv+4.8f*(rhinv*rhinv)-4.f/3.f*(rhinv*rhinv*rhinv)-1.f/120.f/(rhinv*rhinv*rhinv);
      ids  = 1.6f-1/(30.f*rhinv)-(rhinv*rhinv)*(16.f/3.f+rhinv*(-8.f+rhinv*(4.8f-rhinv*16.f/15.f)));
    }//end if rhin < 0.5

    ids  = ids*2.f*epseffi;
    ids3 =  ids3*8.f*(epseffi*epseffi*epseffi);

    //Self gravity prevention (and NaN prevention)
    ids  *= selfGrav;
    ids3 *= selfGrav;

  } //end dist >= epseff
  
  ids2 = ids*ids;
  float ids5 = ids3*ids2;
  float ids7 = ids5*ids2;
  return (float4){ids, -ids3, +3.0f*ids5, -15.0f*ids7};
} 

#else

__device__ float4 get_D04(float ds2, int selfGrav = 1) {
#if 1
  float ids  = rsqrtf(ds2);  //Does not work with zero-softening
  //   if(isnan(ids)) ids = 0;               //This does work with zero-softening, few percent performance drop
  //float ids  = (1.0f / sqrtf(ds2)) * selfGrav; Slower in Pre CUDA4.1
  ids *= selfGrav;
#else
  const float ids = (ds2 > 0.0f) ? rsqrtf(ds2) : 0.0f;
#endif
  const float ids2 = ids*ids;
  float ids3 = ids *ids2;
  float ids5 = ids3*ids2;
  float ids7 = ids5*ids2;
  return make_float4(ids, -ids3, +3.0f*ids5, -15.0f*ids7);
}  // 9 flops

#endif

#ifdef INDSOFT
__device__ float4 add_acc(float4 acc,  float4 pos,
                          float massj, float3 posj, float epsQ,
                          float &ds2, float eps2P, int selfGrav) {

#else

__device__ float4 add_acc(float4 acc,  float4 pos,
			  float massj, float3 posj,
			  float &ds2, float eps2, int selfGrav) {
#endif

  float3 dr = {pos.x - posj.x,
	       pos.y - posj.y,
	       pos.z - posj.z};

  #ifdef INDSOFT
    ds2 = dr.x*dr.x + dr.y*dr.y + dr.z*dr.z;
    float4 D04 = get_D04(ds2, eps2P, epsQ, selfGrav);
  #else
    ds2 = dr.x*dr.x + dr.y*dr.y + dr.z*dr.z + eps2;
    float4 D04 = get_D04(ds2, selfGrav);
  #endif

  float  D0  = D04.x*massj;
  float  D1  = D04.y*massj;

  acc.w -= D0;
  acc.x += D1*dr.x;
  acc.y += D1*dr.y;
  acc.z += D1*dr.z;

  return acc;
}

#ifdef INDSOFT
__device__ float4 add_acc(float4 acc, float4 pos,
                          float mass, float3 com,
                          float3 Q0,  float3 Q1, float epsNode, float eps2) {
#else

__device__ float4 add_acc(float4 acc, float4 pos,
			  float mass, float3 com,
			  float3 Q0,  float3 Q1, float eps2) {
#endif
  float3 dr = {pos.x - com.x,
	       pos.y - com.y,
	       pos.z - com.z};

  #ifdef INDSOFT
    float  ds2 = dr.x*dr.x + dr.y*dr.y + dr.z*dr.z;
    float4 D04 = get_D04(ds2, epsNode, eps2);
  #else
    float  ds2 = dr.x*dr.x + dr.y*dr.y + dr.z*dr.z + eps2;
    float4 D04 = get_D04(ds2);
  #endif

  float  D0  = D04.x*mass;
  float  D1  = D04.y*mass;
  float  D2  = D04.z*mass;
  float  D3  = D04.w*mass;

  float oct_q11 = Q0.x;
  float oct_q22 = Q0.y;
  float oct_q33 = Q0.z;
  float oct_q12 = Q1.x;
  float oct_q13 = Q1.y;
  float oct_q23 = Q1.z;

  float Qii = oct_q11 + oct_q22 + oct_q33;
  float QijRiRj =
          (oct_q11*dr.x*dr.x + oct_q22*dr.y*dr.y + oct_q33*dr.z*dr.z) +
    2.0f*(oct_q12*dr.y*dr.x + oct_q13*dr.z*dr.x + oct_q23*dr.y*dr.z);
//2.0f was
  //volatile float QijRiRj_1 = oct_q12*dr.y*dr.x + oct_q13*dr.z*dr.x + oct_q23*dr.y*dr.z;
 //QijRiRj = dr.x*dr.x;// + oct_q22*dr.y*dr.y + oct_q33*dr.z*dr.z;  //oct_q11;//*

  acc.w        -= D0 + 0.5f*D1*Qii + 0.5f*D2*QijRiRj;
  float C01a    = D1 + 0.5f*D2*Qii + 0.5f*D3*QijRiRj;
  acc.x         += C01a*dr.x + D2*(oct_q11*dr.x + oct_q12*dr.y + oct_q13*dr.z);
  acc.y         += C01a*dr.y + D2*(oct_q12*dr.x + oct_q22*dr.y + oct_q23*dr.z);
  acc.z         += C01a*dr.z + D2*(oct_q13*dr.x + oct_q23*dr.y + oct_q33*dr.z);

  return acc;
}


//Minimum distance opening criteria
#ifdef INDSOFT
  __device__ bool split_node_grav_md(float4 nodeCenter, float4 nodeSize, float4 groupCenter, float4 groupSize,
                                     float group_eps, float node_eps)
#else
  __device__ bool split_node_grav_md(float4 nodeCenter, float4 nodeSize, float4 groupCenter, float4 groupSize)
#endif
{
  //Compute the distance between the group and the cell
  float3 dr = {fabs(groupCenter.x - nodeCenter.x) - (groupSize.x + nodeSize.x),
               fabs(groupCenter.y - nodeCenter.y) - (groupSize.y + nodeSize.y),
               fabs(groupCenter.z - nodeCenter.z) - (groupSize.z + nodeSize.z)};

  dr.x += fabs(dr.x); dr.x *= 0.5f;
  dr.y += fabs(dr.y); dr.y *= 0.5f;
  dr.z += fabs(dr.z); dr.z *= 0.5f;

  //Distance squared, no need to do sqrt since opening criteria has been squared
  float ds2    = dr.x*dr.x + dr.y*dr.y + dr.z*dr.z;

  #ifdef INDSOFT
    //Naar idee van Inti nu minder overbodige openingen
    if(ds2 <= ((group_eps + node_eps ) * (group_eps + node_eps) )) return true;
  #endif

  return (ds2 <= fabs(nodeCenter.w));
}

//Improved Barnes Hut criterium
#ifdef INDSOFT
__device__ bool split_node_grav_impbh(float4 nodeCOM, float4 groupCenter, float4 groupSize,
                                      float group_eps, float node_eps)
#else
__device__ bool split_node_grav_impbh(float4 nodeCOM, float4 groupCenter, float4 groupSize)
#endif
{
  //Compute the distance between the group and the cell
  float3 dr = {fabs(groupCenter.x - nodeCOM.x) - (groupSize.x),
               fabs(groupCenter.y - nodeCOM.y) - (groupSize.y),
               fabs(groupCenter.z - nodeCOM.z) - (groupSize.z)};

  dr.x += fabs(dr.x); dr.x *= 0.5f;
  dr.y += fabs(dr.y); dr.y *= 0.5f;
  dr.z += fabs(dr.z); dr.z *= 0.5f;

  //Distance squared, no need to do sqrt since opening criteria has been squared
  float ds2    = dr.x*dr.x + dr.y*dr.y + dr.z*dr.z;

  #ifdef INDSOFT
    //Extra test
    if(ds2 <= ((group_eps + node_eps ) * (group_eps + node_eps) )) return true;
  #endif

  return (ds2 <= fabs(nodeCOM.w));
}



      int c_stack0 = n_stack_pre;

      /*********** walk a level **********/
      while(c_stack0 < n_stack0) 
      {

        /***
         **** --> fetch the list of nodes rom LMEM
         ***/
        bool use_node = tid <  n_nodes0;
        { prefix[tid] = nstack[ACCS<SHIFT>(c_stack0)];   c_stack0++; }
        __syncthreads();
        int node  = prefix[min(tid, n_nodes0 - 1)];

        if(n_nodes0 > 0){       //Work around pre 4.1 compiler bug
          n_nodes0 -= DIM;
        }

        /***
         **** --> process each of the nodes in the list in parallel
         ***/

#ifndef TEXTURES
        float4 nodeSize = boxSizeInfo[node];                   //Fetch the size of the box. Size.w = child info
        float4 node_pos = boxCenterInfo[node];                 //Fetch the center of the box. center.w = opening info
#else
        float4 nodeSize =  tex1Dfetch(texNodeSize, node);
        float4 node_pos =  tex1Dfetch(texNodeCenter, node);
#endif

        int node_data = __float_as_int(nodeSize.w);

        //Check if a cell has to be opened
#ifndef TEXTURES
        float4 nodeCOM = multipole_data[node*3];
#else
        float4 nodeCOM = tex1Dfetch(texMultipole,node*3);
#endif

        nodeCOM.w      = node_pos.w;
        bool   split   = split_node_grav_impbh(nodeCOM, group_pos, groupSize);


        bool leaf       = node_pos.w <= 0;  //Small AND equal incase of a 1 particle cell       //Check if it is a leaf
        //         split = true;


        uint mask    = BTEST((split && !leaf) && use_node);               // mask = #FFFFFFFF if use_node+split+not_a_leaf==true, otherwise zero
        int child    =    node_data & 0x0FFFFFFF;                         //Index to the first child of the node
        int nchild   = (((node_data & 0xF0000000) >> 28)) & mask;         //The number of children this node has


        /***
         **** --> calculate prefix
         ***/

        int *prefix0 = &prefix[  0];
        int *prefix1 = &prefix[DIM];

#ifdef OLDPREFIX
        int n_total = calc_prefix<DIM2>(prefix, tid,  nchild);
        prefix[tid] += n_offset - nchild;
        __syncthreads();
#else
        inclusive_scan_block<ADDOP<int>, int>(prefix, nchild, tid);        // inclusive scan to compute memory offset of each child
        int n_total = prefix[blockDim.x - 1];                              // fetch total number of children, i.e. offset of the last child -1
        __syncthreads();                                                   // thread barrier to make sure that warps completed their jobs
        prefix[tid] += n_offset - nchild;                                  // convert inclusive into exclusive scan for referencing purpose
        __syncthreads();                                                   // thread barrier
#endif

        for (int i = n_offset; i < n_offset + n_total; i += DIM)         //nullify part of the array that will be filled with children
          nodes[tid + i] = 0;                                          //but do not touch those parts which has already been filled
        __syncthreads();                                                 //Thread barrier to make sure all warps finished writing data

        bool flag = (split && !leaf) && use_node;                        //Flag = use_node + split + not_a_leaf;Use only non_leaf nodes that are to be split
