       sh_pos [laneId] = make_float3(1.0e10f, 1.0e10f, 1.0e10f);
 
     }
-#pragma unroll 16
+#ifndef _QUADRUPOLE_
     for (int i = 0; i < WARP_SIZE; i++)
       for (int k = 0; k < NI; k++)
-        acc_i[k] = add_acc(acc_i[k], pos_i[k], sh_mass[offs + i], sh_pos[offs+i],eps2);
+        acc_i[k] = add_acc(acc_i[k], pos_i[k], sh_mass[i], sh_pos[i],eps2);
+#else
+    for (int i = 0; i < WARP_SIZE; i++)
+    {
+      float4 Q0, Q1;
+      Q0 = Q1 = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
+      if (i < n_approx)
+      {
+        const int address = approxM[i] * 3;
+        Q0 = tex1Dfetch(texMultipole, address + 1);
+        Q1 = tex1Dfetch(texMultipole, address + 2);
+      }
+      for (int k = 0; k < NI; k++)
+        acc_i[k] = add_acc(acc_i[k], pos_i[k], sh_mass[i], sh_pos[i], Q0, Q1, eps2);
+    }
+#endif
     apprCount += WARP_SIZE*NI;
 
   } //if n_approx > 0


  /*********** set necessary thread constants **********/

  const int offs = 0;
  const int DIM2 = WARP_SIZE2;
  const int DIM  = WARP_SIZE;



/*********** Forces *************/

__device__ __forceinline__ float4 add_acc(
    float4 acc,  const float4 pos,
    const float massj, const float3 posj,
    const float eps2)
