                                                                                             n_bodies,
                                                                                             nmax,
                                                                                             window_size,
+                                                                                            thread_mask,
+                                                                                            n_bodies_per_block,
                                                                                             box);
 
     return cudaSuccess;

                                                 unsigned int n_bodies, 
                                                 unsigned int nmax,
                                                 unsigned int window_size,
                                                 gpu_boxsize box)
    {
    int group_idx = blockIdx.x;
    
    float3 *body_force = sum;
    float3 *body_torque = &sum[window_size];
        
    __shared__ float4 ex_space, ey_space, ez_space;
    float3 sum_force = make_float3(0.0f, 0.0f, 0.0f);
    float3 sum_torque = make_float3(0.0f, 0.0f, 0.0f);
    
    int idx_body = d_rigid_group[group_idx];

    if (threadIdx.x == 0)
        {
        float4 orientation = d_rigid_orientation[idx_body];
        exyzFromQuaternion(orientation, ex_space, ey_space, ez_space);
        }
    
    __syncthreads();

                // 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;
