     unsigned int n_group_bodies = rigid_data.n_group_bodies;
     unsigned int nmax = rigid_data.nmax;
     
-    // run the kernel: the shared memory size is used for dynamic memory allocation of extern __shared__ sum
-    unsigned int window_size;
-    if (nmax < 64)
-        window_size = 32;          // some fixed value divisble by nmax
+    unsigned int n_bodies_per_block;
+    unsigned int target_num_blocks = 60;
+    if (n_group_bodies / 8 >= target_num_blocks)
+        n_bodies_per_block = 8;
     else
-        window_size = 64;
-    
-    unsigned int block_size = window_size;
-    dim3 force_grid(n_group_bodies, 1, 1);
+    if (n_group_bodies / 4 >= target_num_blocks)
+        n_bodies_per_block = 4;
+    else
+    if (n_group_bodies / 2 >= target_num_blocks)
+        n_bodies_per_block = 2;
+    else
+        n_bodies_per_block = 1;
+
+    unsigned int block_size = 64;
+    unsigned int window_size = block_size / n_bodies_per_block;
+    unsigned int thread_mask = window_size - 1;
+    
+    dim3 force_grid(n_group_bodies / n_bodies_per_block + 1, 1, 1);
     dim3 force_threads(block_size, 1, 1);
 
-    gpu_rigid_force_sliding_kernel<<< force_grid, force_threads, 2 * window_size * sizeof(float3) >>>(rigid_data.force, 
+    gpu_rigid_force_sliding_kernel<<< force_grid, force_threads, 2 * block_size * sizeof(float3) >>>(rigid_data.force, 
                                                                                             rigid_data.torque,
                                                                                             rigid_data.body_indices,
                                                                                             rigid_data.orientation,

                                                                                            n_bodies,
                                                                                            nmax,
                                                                                            window_size,
                                                                                            box);

    return cudaSuccess;

                // This might require more calculations but more stable
                // particularly when rigid bodies are bigger than half the box
                float3 ri;
                ri.x = ex_space[m].x * particle_pos.x + ey_space[m].x * particle_pos.y 
                        + ez_space[m].x * particle_pos.z;
                ri.y = ex_space[m].y * particle_pos.x + ey_space[m].y * particle_pos.y 
                        + ez_space[m].y * particle_pos.z;
                ri.z = ex_space[m].z * particle_pos.x + ey_space[m].z * particle_pos.y 
                        + ez_space[m].z * particle_pos.z;
    
                sum_torque.x += ri.y * fi.z - ri.z * fi.y;
                sum_torque.y += ri.z * fi.x - ri.x * fi.z;
