 }
 
 
-extern "C" __global__ void
+  extern "C" __global__ void
 __launch_bounds__(NTHREAD)
-dev_approximate_gravity(const int n_active_groups,
-                        int    n_bodies,
-                        float eps2,
-                        uint2 node_begend,
-                        int    *active_groups,
-                        real4  *body_pos,
-                        real4  *multipole_data,
-                        float4 *acc_out,
-                        int    *ngb_out,
-                        int    *active_inout,
-                        int2   *interactions,
-                        float4  *boxSizeInfo,
-                        float4  *groupSizeInfo,
-                        float4  *boxCenterInfo,
-                        float4  *groupCenterInfo,
-                        real4   *body_vel,
-                        int     *MEM_BUF) {
-//                                                    int     grpOffset){
-
-
-  const int blockDim2 = NTHREAD2;
-   __shared__ int shmem[15*(1 << blockDim2)];
-//    __shared__ int shmem[24*(1 << blockDim2)]; is possible on FERMI
-//    int             lmem[LMEM_STACK_SIZE];
-
-
-
-  /*********** check if this block is linked to a leaf **********/
-
-  int bid = gridDim.x * blockIdx.y + blockIdx.x;
-
-  while(true)
-  {
-
-    if(threadIdx.x == 0)
+  dev_approximate_gravity(const int n_active_groups,
+      int    n_bodies,
+      float eps2,
+      uint2 node_begend,
+      int    *active_groups,
+      real4  *body_pos,
+      real4  *multipole_data,
+      float4 *acc_out,
+      int    *ngb_out,
+      int    *active_inout,
+      int2   *interactions,
+      float4  *boxSizeInfo,
+      float4  *groupSizeInfo,
+      float4  *boxCenterInfo,
+      float4  *groupCenterInfo,
+      real4   *body_vel,
+      int     *MEM_BUF) {
+    //                                                    int     grpOffset){
+
+
+    const int blockDim2 = NTHREAD2;
+    __shared__ int shmem[15*(1 << blockDim2)];
+    //    __shared__ int shmem[24*(1 << blockDim2)]; is possible on FERMI
+    //    int             lmem[LMEM_STACK_SIZE];
+
+
+
+    /*********** check if this block is linked to a leaf **********/
+
+    int bid = gridDim.x * blockIdx.y + blockIdx.x;
+
+    while(true)
     {
-      bid         = atomicAdd(&active_inout[n_bodies], 1);
-      shmem[0]    = bid;
-    }
-  __syncthreads();
-
-  bid   = shmem[0];
-
-  if (bid >= n_active_groups) return;
-
-
-  int tid = threadIdx.y * blockDim.x + threadIdx.x;
-
-  int grpOffset = 0;
-
-//   volatile int *lmem = &MEM_BUF[blockIdx.x*LMEM_STACK_SIZE*blockDim.x + threadIdx.x*LMEM_STACK_SIZE];
-//   int *lmem = &MEM_BUF[blockIdx.x*LMEM_STACK_SIZE*blockDim.x + threadIdx.x*LMEM_STACK_SIZE];
-   int *lmem = &MEM_BUF[blockIdx.x*LMEM_STACK_SIZE*blockDim.x];
-
-
-  /*********** set necessary thread constants **********/
-
-  real4 curGroupSize    = groupSizeInfo[active_groups[bid + grpOffset]];
-  int   groupData       = __float_as_int(curGroupSize.w);
-  uint body_i           =   groupData & CRITMASK;
-  uint nb_i             = ((groupData & INVCMASK) >> CRITBIT) + 1;
-
-  real4 group_pos       = groupCenterInfo[active_groups[bid + grpOffset]];
-
-//   if(tid == 0)
-//   printf("[%f %f %f %f ] \n [%f %f %f %f ] %d %d \n",
-//           curGroupSize.x, curGroupSize.y, curGroupSize.z, curGroupSize.w,
-//           group_pos.x, group_pos.y, group_pos.z, group_pos.w, body_i, nb_i);
-
-
-  int DIM2x = 0;
-  while (((nb_i - 1) >> DIM2x) > 0) DIM2x++;
-
-  DIM2x     = max(DIM2x,4);
-  int DIM2y = blockDim2 - DIM2x;
-
-  int tx = tid & ((1 << DIM2x) - 1);
-  int ty = tid >> DIM2x;
-
-  body_i += tx%nb_i;
-
-  //float4 pos_i = tex1Dfetch(bodies_pos_ref, body_i);   // texture read: 4 floats
-
-
-  float4 pos_i = body_pos[body_i];
-
-
-  int ngb_i;
-
-  float4 acc_i = {0.0f, 0.0f, 0.0f, 0.0f};
-
-  #ifdef INDSOFT
-    eps2 = body_vel[body_i].w;
-    float group_eps = eps2;
-
-    volatile float *reduc = (float*) &shmem[0];
-    reduc[threadIdx.x] = eps2;
-
-    //Find the maximum softening value for the particles in this group
-    __syncthreads();
-    // do reduction in shared mem
-    if(blockDim.x >= 512) if (tid < 256) {reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 256]);} __syncthreads();
-    if(blockDim.x >= 256) if (tid < 128) {reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 128]);} __syncthreads();
-    if(blockDim.x >= 128) if (tid < 64)  {reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 64]);} __syncthreads();
-    if(blockDim.x >= 64) if (tid < 32) { reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 32]);}
-    if(blockDim.x >= 32) if (tid < 16) { reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 16]);}
-
-    if(tid < 8)
-    {
-      reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 8]);
-      reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 4]);
-      reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 2]);
-      reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 1]);
-    }
-    __syncthreads();
-
-    group_eps = reduc[0];
-  #else
-    float group_eps  = 0;
-  #endif
-
-  int apprCount = 0;
-  int direCount = 0;
-
-
-  acc_i = approximate_gravity<blockDim2, 0>( DIM2x, DIM2y, tid, tx, ty,
-                                          body_i, pos_i, group_pos,
-                                          eps2, node_begend,
-                                          multipole_data, body_pos,
-                                          shmem, lmem, ngb_i, apprCount, direCount, boxSizeInfo, curGroupSize, boxCenterInfo,
-                                         group_eps, body_vel);
-  if(apprCount < 0)
-  {
-
-    //Try to get access to the big stack, only one block per time is allowed
-    if(threadIdx.x == 0)
-    {
-      int res = atomicExch(&active_inout[n_bodies+1], 1); //If the old value (res) is 0 we can go otherwise sleep
-      int waitCounter  = 0;
-      while(res != 0)
+
+      if(threadIdx.x == 0)
       {
-          //Sleep
-          for(int i=0; i < (1024); i++)
+        bid         = atomicAdd(&active_inout[n_bodies], 1);
+        shmem[0]    = bid;
+      }
+      __syncthreads();
+
+      bid   = shmem[0];
+
+      if (bid >= n_active_groups) return;
+
+
+      int tid = threadIdx.y * blockDim.x + threadIdx.x;
+
+      int grpOffset = 0;
+
+      //   volatile int *lmem = &MEM_BUF[blockIdx.x*LMEM_STACK_SIZE*blockDim.x + threadIdx.x*LMEM_STACK_SIZE];
+      //   int *lmem = &MEM_BUF[blockIdx.x*LMEM_STACK_SIZE*blockDim.x + threadIdx.x*LMEM_STACK_SIZE];
+      int *lmem = &MEM_BUF[blockIdx.x*LMEM_STACK_SIZE*blockDim.x];
+
+
+      /*********** set necessary thread constants **********/
+
+      real4 curGroupSize    = groupSizeInfo[active_groups[bid + grpOffset]];
+      int   groupData       = __float_as_int(curGroupSize.w);
+      uint body_i           =   groupData & CRITMASK;
+      uint nb_i             = ((groupData & INVCMASK) >> CRITBIT) + 1;
+
+      real4 group_pos       = groupCenterInfo[active_groups[bid + grpOffset]];
+
+      //   if(tid == 0)
+      //   printf("[%f %f %f %f ] \n [%f %f %f %f ] %d %d \n",
+      //           curGroupSize.x, curGroupSize.y, curGroupSize.z, curGroupSize.w,
+      //           group_pos.x, group_pos.y, group_pos.z, group_pos.w, body_i, nb_i);
+
+
+      int DIM2x = 0;
+      while (((nb_i - 1) >> DIM2x) > 0) DIM2x++;
+
+      DIM2x     = max(DIM2x,4);
+      int DIM2y = blockDim2 - DIM2x;
+
+      int tx = tid & ((1 << DIM2x) - 1);
+      int ty = tid >> DIM2x;
+
+      body_i += tx%nb_i;
+
+      //float4 pos_i = tex1Dfetch(bodies_pos_ref, body_i);   // texture read: 4 floats
+
+
+      float4 pos_i = body_pos[body_i];
+
+
+      int ngb_i;
+
+      float4 acc_i = {0.0f, 0.0f, 0.0f, 0.0f};
+
+#ifdef INDSOFT
+      eps2 = body_vel[body_i].w;
+      float group_eps = eps2;
+
+      volatile float *reduc = (float*) &shmem[0];
+      reduc[threadIdx.x] = eps2;
+
+      //Find the maximum softening value for the particles in this group
+      __syncthreads();
+      // do reduction in shared mem
+      if(blockDim.x >= 512) if (tid < 256) {reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 256]);} __syncthreads();
+      if(blockDim.x >= 256) if (tid < 128) {reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 128]);} __syncthreads();
+      if(blockDim.x >= 128) if (tid < 64)  {reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 64]);} __syncthreads();
+      if(blockDim.x >= 64) if (tid < 32) { reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 32]);}
+      if(blockDim.x >= 32) if (tid < 16) { reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 16]);}
+
+      if(tid < 8)
+      {
+        reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 8]);
+        reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 4]);
+        reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 2]);
+        reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 1]);
+      }
+      __syncthreads();
+
+      group_eps = reduc[0];
+#else
+      float group_eps  = 0;
+#endif
+
+      int apprCount = 0;
+      int direCount = 0;
+
+
+      acc_i = approximate_gravity<blockDim2, 0>( DIM2x, DIM2y, tid, tx, ty,
+          body_i, pos_i, group_pos,
+          eps2, node_begend,
+          multipole_data, body_pos,
+          shmem, lmem, ngb_i, apprCount, direCount, boxSizeInfo, curGroupSize, boxCenterInfo,
+          group_eps, body_vel);
+      if(apprCount < 0)
+      {
+
+        //Try to get access to the big stack, only one block per time is allowed
+        if(threadIdx.x == 0)
+        {
+          int res = atomicExch(&active_inout[n_bodies+1], 1); //If the old value (res) is 0 we can go otherwise sleep
+          int waitCounter  = 0;
+          while(res != 0)
           {
-                  waitCounter += 1;
+            //Sleep
+            for(int i=0; i < (1024); i++)
+            {
+              waitCounter += 1;
+            }
+            //Test again
+            shmem[0] = waitCounter;
+            res = atomicExch(&active_inout[n_bodies+1], 1); 
           }
-          //Test again
-          shmem[0] = waitCounter;
-          res = atomicExch(&active_inout[n_bodies+1], 1); 
+        }
+
+        __syncthreads();
+
+        lmem = &MEM_BUF[gridDim.x*LMEM_STACK_SIZE*blockDim.x];    //Use the extra large buffer
+        apprCount = direCount = 0;
+        acc_i = approximate_gravity<blockDim2, 8>( DIM2x, DIM2y, tid, tx, ty,
+            body_i, pos_i, group_pos,
+            eps2, node_begend,
+            multipole_data, body_pos,
+            shmem, lmem, ngb_i, apprCount, direCount, boxSizeInfo, curGroupSize, boxCenterInfo,
+            group_eps, body_vel);
+
+        lmem = &MEM_BUF[blockIdx.x*LMEM_STACK_SIZE*blockDim.x]; //Back to normal location
+
+        if(threadIdx.x == 0)
+        {
+          atomicExch(&active_inout[n_bodies+1], 0); //Release the lock
+        }
+      }//end if apprCount < 0
+
+      if (tid < nb_i) {
+        acc_out     [body_i] = acc_i;
+        ngb_out     [body_i] = ngb_i;
+        active_inout[body_i] = 1;
+        interactions[body_i].x = apprCount;
+        interactions[body_i].y = direCount ;
       }
-    }
-
-    __syncthreads();
-
-    lmem = &MEM_BUF[gridDim.x*LMEM_STACK_SIZE*blockDim.x];    //Use the extra large buffer
-    apprCount = direCount = 0;
-    acc_i = approximate_gravity<blockDim2, 8>( DIM2x, DIM2y, tid, tx, ty,
-                                            body_i, pos_i, group_pos,
-                                            eps2, node_begend,
-                                            multipole_data, body_pos,
-                                            shmem, lmem, ngb_i, apprCount, direCount, boxSizeInfo, curGroupSize, boxCenterInfo,
-                                          group_eps, body_vel);
-
-    lmem = &MEM_BUF[blockIdx.x*LMEM_STACK_SIZE*blockDim.x]; //Back to normal location
-
-    if(threadIdx.x == 0)
-    {
-            atomicExch(&active_inout[n_bodies+1], 0); //Release the lock
-    }
-  }//end if apprCount < 0
-
-  if (tid < nb_i) {
-      acc_out     [body_i] = acc_i;
-      ngb_out     [body_i] = ngb_i;
-      active_inout[body_i] = 1;
-      interactions[body_i].x = apprCount;
-      interactions[body_i].y = direCount ;
-    }
-
-
-  }     //end while
-}
-
-
+
+
+    }     //end while
+  }
+
+

      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
    }

    __syncthreads();
    #ifdef INDSOFT
      node_eps[tid]  = temp;
      __syncthreads();
    #endif
  } //if n_approx > 0

  if(n_direct > 0)

#include <stdio.h>

#include "node_specs.h"

__forceinline__ __device__ float Wkernel(const float q)
{
  const float sigma = 8.0f/M_PI;

  const float qm = 1.0f - q;
  const float f1 = sigma * (1.0f + (-6.0f)*q*q*qm);
  const float f2 = sigma * 2.0f*qm*qm*qm;

  return fmaxf(0.0f, fminf(f1, f2));
}

__forceinline__ __device__ float interact(
    const float3 ipos,
    const float  h,
    const float  hinv,
    const float3 jpos,
    const float  jmass)
{
  const float3 dr = make_float3(jpos.x - ipos.x, jpos.y - ipos.y, jpos.z - ipos.z);
  const float  r2 = dr.x*dr.x + dr.y*dr.y + dr.z*dr.z;
  if (r2 >= h*h) return 0.0f;
  const float q  = sqrtf(r2) * hinv;
  const float hinv3 = hinv*hinv*hinv;

  return jmass * Wkernel(q) * hinv3;
}


/***
