       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;
-
+      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);

}

template <typename Vec, typename Dtype>
__global__ void force_kernel_memopt(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
      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);
    }
    atomicAdd(&p[tid].x, pfx);
    atomicAdd(&p[tid].y, pfy);
