     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,

        __syncthreads();
        }
    
    if (threadIdx.x == 0)
        {
        rdata_force[idx_body] = make_float4(body_force[0].x, body_force[0].y, body_force[0].z, 0.0f);
        rdata_torque[idx_body] = make_float4(body_torque[0].x, body_torque[0].y, body_torque[0].z, 0.0f);
        }
    }


        {
        float4 fi = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
        
        unsigned int k = start * window_size + (threadIdx.x & thread_mask);
        
        if (k < nmax && idx_body[m] != -1)
            {
            int localidx = idx_body[m] * nmax + k;
            unsigned int pidx = d_rigid_particle_idx[localidx];
                        
            if (pidx != INVALID_INDEX)
