     \param nmax Maximum number of particles in a rigid body
     \param window_size Window size for reduction
     \param box Box dimensions for periodic boundary condition handling
+    
+    Compute the force and torque sum on all bodies in the system from their constituent particles. n_bodies_per_block
+    bodies are handled within each block of execution on the GPU. The reason for this is to decrease
+    over-parallelism and use the GPU cores more effectively when bodies are smaller than the block size. Otherwise,
+    small bodies leave many threads in the block idle with nothing to do.
+    
+    On start, the properties common to each body are read in, computed, and stored in shared memory for all the threads
+    working on that body to access. Then, the threads loop over all particles that are part of the body with
+    a sliding window. Each loop of the window computes the force and torque for block_size/n_bodies_per_block particles
+    in as many threads in parallel. These quantities are summed over enough windows to cover the whole body.
+    
+    The block_size/n_bodies_per_block partial sums are stored in shared memory. Then n_bodies_per_block partial
+    reductions are performed in parallel using all threads to sum the total force and torque on each body. This looks
+    just like a normal reduction, except that it terminates at a certain level in the tree. To make the math
+    for the partial reduction work out, block_size must be a power of 2 as must n_bodies_per_block.
+    
+    Performance testing on GF100 with many different bodies of different sizes ranging from 4-256 particles per body
+    has found that the optimum block size for most bodies is 64 threads. Performance increases for all body sizes
+    as n_bodies_per_block is increased, but only up to 8. n_bodies_per_block=16 slows performance significantly.
+    Based on these performance results, this kernel is hardcoded to handle only 1,2,4,8 n_bodies_per_block
+    with a power of 2 block size (hardcoded to 64 in the kernel launch).
+    
+    However, there is one issue to the n_bodies_per_block parallelism reduction. If the reduction results in too few
+    blocks, performance can actually be reduced. For example, if there are only 64 bodies running at the "most optimal"
+    n_bodies_per_block=8 results in only 8 blocks on the GPU! That isn't even enough to heat up all 15 SMs on GF100.
+    Even though n_bodies_per_block=1 is not fully optimal per block, running 64 slow blocks in parallel is faster than
+    running 8 fast blocks in parallel. Testing on GF100 determines that 60 blocks is the dividing line (makes sense - 
+    that's 4 blocks active on each SM).
 */
 extern "C" __global__ void gpu_rigid_force_sliding_kernel(float4* rdata_force, 
                                                  float4* rdata_torque,

        __syncthreads();
        }
    
    if ((threadIdx.x & thread_mask) == 0 && idx_body[m] != -1)
        {
        rdata_force[idx_body[m]] = make_float4(body_force[threadIdx.x].x, body_force[threadIdx.x].y, body_force[threadIdx.x].z, 0.0f);

    
    __syncthreads();
    
    // compute the number of windows that we need to loop over
    unsigned int n_windows = nmax / window_size;
    if (n_windows == 0)
        n_windows = 1;
