 
     int block_count_nodes = (num_movable_nodes + thread_count - 1) / thread_count;
     computeElectricForce<<<block_count_nodes, thread_count, 0, stream_movable>>>(
+        num_bins_x, num_bins_y,
+        num_movable_impacted_bins_x, num_movable_impacted_bins_y,
+        field_map_x_tensor, field_map_y_tensor,
+        x_tensor, y_tensor,
+        node_size_x_tensor, node_size_y_tensor,
+        bin_center_x_tensor, bin_center_y_tensor,
+        xl, yl, xh, yh,
+        bin_size_x, bin_size_y,
+        num_movable_nodes,
+        grad_x_tensor, grad_y_tensor);
+
+    if (num_filler_nodes)
+    {
+        status = cudaStreamCreate(&stream_filler);
+        if (status != cudaSuccess)
+        {
+            printf("cudaStreamCreate failed for stream_filler\n");
+            fflush(stdout);
+            return 1;
+        }
+
+        block_count_nodes = (num_filler_nodes + thread_count - 1) / thread_count;
+        int offset = num_nodes - num_filler_nodes;
+        computeElectricForce<<<block_count_nodes, thread_count, 0, stream_filler>>>(
             num_bins_x, num_bins_y,
-            num_movable_impacted_bins_x, num_movable_impacted_bins_y,
+            num_filler_impacted_bins_x, num_filler_impacted_bins_y,
             field_map_x_tensor, field_map_y_tensor,
-            x_tensor, y_tensor,
-            node_size_x_tensor, node_size_y_tensor,
+            x_tensor + offset, y_tensor + offset,
+            node_size_x_tensor + offset, node_size_y_tensor + offset,
             bin_center_x_tensor, bin_center_y_tensor,
             xl, yl, xh, yh,
             bin_size_x, bin_size_y,
-            num_movable_nodes,
-            grad_x_tensor, grad_y_tensor
-            );
-
-    if (num_filler_nodes)
-    {
-        status = cudaStreamCreate(&stream_filler);
-        if (status != cudaSuccess)
-        {
-            printf("cudaStreamCreate failed for stream_filler\n");
-            fflush(stdout);
-            return 1;
-        }
-
-        block_count_nodes = (num_filler_nodes + thread_count - 1) / thread_count;
-        computeElectricForce<<<block_count_nodes, thread_count, 0, stream_filler>>>(
-                num_bins_x, num_bins_y,
-                num_filler_impacted_bins_x, num_filler_impacted_bins_y,
-                field_map_x_tensor, field_map_y_tensor,
-                x_tensor+num_nodes-num_filler_nodes, y_tensor+num_nodes-num_filler_nodes,
-                node_size_x_tensor+num_nodes-num_filler_nodes, node_size_y_tensor+num_nodes-num_filler_nodes,
-                bin_center_x_tensor, bin_center_y_tensor,
-                xl, yl, xh, yh,
-                bin_size_x, bin_size_y,
-                num_filler_nodes,
-                grad_x_tensor+num_nodes-num_filler_nodes, grad_y_tensor+num_nodes-num_filler_nodes
-                );
+            num_filler_nodes,
+            grad_x_tensor + offset, grad_y_tensor + offset);
 
         status = cudaStreamDestroy(stream_filler);
         stream_filler = 0;

                            x_tensor[i], y_tensor[i], x_tensor[i]+node_size_x_tensor[i], y_tensor[i]+node_size_y_tensor[i]
                          );
                }
#endif

            }
        }
    }
}

template <typename T>
int computeElectricForceCudaLauncher(
        int num_bins_x, int num_bins_y,
        int num_movable_impacted_bins_x, int num_movable_impacted_bins_y,
        int num_filler_impacted_bins_x, int num_filler_impacted_bins_y,
        const T* field_map_x_tensor, const T* field_map_y_tensor,
        const T* x_tensor, const T* y_tensor,
        const T* node_size_x_tensor, const T* node_size_y_tensor,
        const T* bin_center_x_tensor, const T* bin_center_y_tensor,
        T xl, T yl, T xh, T yh,
        T bin_size_x, T bin_size_y,
        int num_nodes, int num_movable_nodes, int num_filler_nodes,
        T* grad_x_tensor, T* grad_y_tensor
        )
{
    int thread_count = 512;


    cudaError_t status;
    cudaStream_t stream_movable;


template <typename T>
__global__ void computeElectricForceAtomic(
    int num_bins_x, int num_bins_y,
    int num_impacted_bins_x, int num_impacted_bins_y,
    const T *field_map_x_tensor, const T *field_map_y_tensor,
    const T *x_tensor, const T *y_tensor,
    const T *node_size_x_tensor, const T *node_size_y_tensor,
    const T *bin_center_x_tensor, const T *bin_center_y_tensor,
    T xl, T yl, T xh, T yh,
    T bin_size_x, T bin_size_y,
    int num_nodes,
    T *grad_x_tensor, T *grad_y_tensor)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < num_nodes * num_impacted_bins_x * num_impacted_bins_y)
    {
        // density overflow function
        auto computeDensityFunc = [](T x, T node_size, T bin_center, T bin_size) {
            //return max(T(0.0), min(x+node_size, bin_center+bin_size/2) - max(x, bin_center-bin_size/2));
            // Yibo: cannot understand why negative overlap is allowed in RePlAce
            return min(x + node_size, bin_center + bin_size / 2) - max(x, bin_center - bin_size / 2);
        };
        int node_id = i / (num_impacted_bins_x * num_impacted_bins_y);
        int residual_index = i - node_id * num_impacted_bins_x * num_impacted_bins_y;

        // stretch node size to bin size
        T node_size_x = max(bin_size_x * SQRT2, node_size_x_tensor[node_id]);
        T node_size_y = max(bin_size_y * SQRT2, node_size_y_tensor[node_id]);
        T node_x = x_tensor[node_id] + node_size_x_tensor[node_id] / 2 - node_size_x / 2;
        T node_y = y_tensor[node_id] + node_size_y_tensor[node_id] / 2 - node_size_y / 2;
        // boundary condition
        //node_x = max(node_x, xl);
        //node_x = min(node_x, xh-node_size_x);
