 
       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);
-
-#ifdef INDSOFT
-      node_eps[tid]  = octopole0.w;
-#endif
 
     } else {
+
       //Set non-active memory locations to zero
       node_mon0[tid] = 0.0f;
       node_mon1[tid] = make_float3(1.0e10f, 1.0e10f, 1.0e10f);
-      node_oct0[tid] = make_float3(0.0f, 0.0f, 0.0f);
-      node_oct1[tid] = make_float3(0.0f, 0.0f, 0.0f);
-
-#ifdef INDSOFT
-      node_eps[tid]  = 0.01f;
-#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
-    }
+      acc_i = add_acc(acc_i, pos_i, node_mon0[offs + i], node_mon1[offs+i],eps2);
+    apprCount += DIMx;
 
     __syncthreads();
-#ifdef INDSOFT
-    node_eps[tid]  = temp;
-    __syncthreads();
-#endif
   } //if n_approx > 0
 
   if(n_direct > 0)
   {
-    if (tid < n_direct) {
-
-      float4 posj = body_pos[direct[tid]];
-
-      //     float4 posj  = tex1Dfetch(texBody, direct[tid]);
+    if (tid < n_direct) 
+    {
+      const float4 posj = body_pos[direct[tid]];
+#if 0
+      const float4 posj  = tex1Dfetch(texBody, direct[tid]);
+#endif
       sh_mass[tid] = posj.w;
       sh_pos [tid] = make_float3(posj.x, posj.y, posj.z);
       sh_jid [tid] = direct[tid];
-
-#ifdef INDSOFT
-      sh_eps [tid] = body_vel[direct[tid]].w;
-#endif
     } else {
       sh_mass[tid] = 0.0f;
       sh_pos [tid] = make_float3(1.0e10f, 1.0e10f, 1.0e10f);
       sh_jid [tid] = -1;
-
-#ifdef INDSOFT
-      sh_eps [tid] = 0.01f;
-#endif
     }
 
     __syncthreads();
 #pragma unroll
-    for (int j = 0; j < DIMx; j++) {
-      if ((sh_jid[offs + j] >= 0)) {
-        int selfGrav = (body_i != sh_jid[offs + j]);
-        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;
-        }
-      }
-    }
+    for (int j = 0; j < DIMx; j++) 
+      if ((sh_jid[offs + j] >= 0)) 
+        acc_i = add_acc(acc_i, pos_i, sh_mass[offs + j], sh_pos[offs + j], eps2);
+#if 0
+    direCount += DIMx;
+#endif
     __syncthreads();
   }
 


  if(n_approx > 0)
  {
#ifdef INDSOFT
    float temp = node_eps[tid];
#endif

    if (tid < n_approx) {
      int address      = (approx[tid] << 1) + approx[tid];
#ifndef TEXTURES
      float4 monopole  = multipole_data[address    ];
      float4 octopole0 = multipole_data[address + 1];

      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
