 }
 
 template <typename Vec, typename Dtype>
+__global__ void force_kernel_ifless_with_aar(const Vec* q,
+                                             Vec* p,
+                                             const int32_t particle_number,
+                                             const Dtype dt,
+                                             const Dtype CL2,
+                                             const int32_t* sorted_list,
+                                             const int32_t* number_of_partners,
+                                             const int32_t* pointer) {
+  const auto tid = threadIdx.x + blockIdx.x * blockDim.x;
+  if (tid < particle_number) {
+    const auto qi = q[tid];
+    const auto np = number_of_partners[tid];
+    const auto kp = pointer[tid];
+
+    for (int32_t k = 0; k < np; k++) {
+      const auto j = sorted_list[kp + k];
+      const auto dx = q[j].x - qi.x;
+      const auto dy = q[j].y - qi.y;
+      const auto dz = q[j].z - qi.z;
+      const auto r2 = dx * dx + dy * dy + dz * dz;
+      const auto r6 = r2 * r2 * r2;
+      auto df = ((24.0 * r6 - 48.0) / (r6 * r6 * r2)) * dt;
+      if (r2 > CL2) df = 0.0;
+      atomicAdd(&p[tid].x, df * dx);
+      atomicAdd(&p[tid].y, df * dy);
+      atomicAdd(&p[tid].z, df * dz);
+      atomicAdd(&p[j].x, -df * dx);
+      atomicAdd(&p[j].y, -df * dy);
+      atomicAdd(&p[j].z, -df * dz);
+    }
+  }
+}
+
+template <typename Vec, typename Dtype>
 __global__ void force_kernel_memopt_with_aar(const Vec*   __restrict__ q,
                                              Vec*         __restrict__ p,
                                              const int32_t particle_number,

    Dtype pfx = 0.0, pfy = 0.0, pfz = 0.0;
    for (int32_t k = 0; k < np; k++) {
      const auto j = __ldg(ptr_list); // use ROC
      const auto dx = q[j].x - qi.x;
      const auto dy = q[j].y - qi.y;
      const auto dz = q[j].z - qi.z;
      const auto r2 = dx * dx + dy * dy + dz * dz;
      const auto r6 = r2 * r2 * r2;
      auto df = ((static_cast<Dtype>(24.0) * r6 - static_cast<Dtype>(48.0)) / (r6 * r6 * r2)) * dt;
      if (r2 > CL2) df = 0.0;

      pfx += df * dx;
      pfy += df * dy;
      pfz += df * dz;

      atomicAdd(&p[j].x, -df * dx);
      atomicAdd(&p[j].y, -df * dy);
      atomicAdd(&p[j].z, -df * dz);

      ptr_list += particle_number;
    }
    atomicAdd(&p[tid].x, pfx);
    atomicAdd(&p[tid].y, pfy);

      const auto dy = q[j].y - qi.y;
      const auto dz = q[j].z - qi.z;
      const auto r2 = dx * dx + dy * dy + dz * dz;
      if (r2 > CL2) continue;
      const auto r6 = r2 * r2 * r2;
      const auto df = ((24.0 * r6 - 48.0) / (r6 * r6 * r2)) * dt;
      pfx += df * dx;
      pfy += df * dy;
      pfz += df * dz;
      atomicAdd(&p[j].x, -df * dx);
      atomicAdd(&p[j].y, -df * dy);
      atomicAdd(&p[j].z, -df * dz);
