     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;
-
+      ptr_list += particle_number;
+      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;
+      if (r2 > CL2) continue;
+      const auto r6 = r2 * r2 * r2;
+      const auto df = ((static_cast<Dtype>(24.0) * r6 - static_cast<Dtype>(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);
-
-      ptr_list += particle_number;
     }
     atomicAdd(&p[tid].x, pfx);
     atomicAdd(&p[tid].y, pfy);

}

template <typename Vec, typename Dtype>
__global__ void force_kernel_memopt(const Vec*   __restrict__ q,
                                    Vec*         __restrict__ p,
                                    const int32_t particle_number,

}

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,
