                                                  unsigned int window_size,
                                                  gpu_boxsize box)
     {
-/*	
-    int idx_particle = blockIdx.x * blockDim.x + threadIdx.x;
-    int idx_body = blockIdx.x + local_beg;
-    
-    float4 *body_force = sum;
-    float4 *body_torque = &sum[window_size];
-    
-    unsigned int idx_particle_index = tex1Dfetch(rigid_data_particle_indices_tex, idx_particle);
-//  float4 com, pos;
-    float4 ri, ex_space, ey_space, ez_space, particle_pos;
-    float4 fi = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
-    float4 torquei = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
-    float4 force2 = make_float4(0.0f, 0.0f, 0.0f, 0.0f); 
-    float4 torque2 = make_float4(0.0f, 0.0f, 0.0f, 0.0f); 
-    
-    unsigned int n_windows = nmax / window_size + 1;
-    
-    // slide the window throughout the block
-    for (unsigned int start = 0; start < n_windows; start++)
-    {
-        unsigned int iwindow = threadIdx.x / window_size;
-        unsigned int ilocal = threadIdx.x % window_size;
-        
-        body_force[ilocal] = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
-        body_torque[ilocal] = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
-            
-        if (idx_body < n_bodies && idx_particle_index != INVALID_INDEX && iwindow == start)
-            {
-            // calculate body force and torques
-       //   com = tex1Dfetch(rigid_data_com_tex, idx_body);
-       //   pos = tex1Dfetch(pdata_pos_tex, idx_particle_index);
-            ex_space = tex1Dfetch(rigid_data_exspace_tex, idx_body);
-            ey_space = tex1Dfetch(rigid_data_eyspace_tex, idx_body);
-            ez_space = tex1Dfetch(rigid_data_ezspace_tex, idx_body);
-            particle_pos = tex1Dfetch(rigid_data_particle_pos_tex, idx_particle);
-        
-            fi = tex1Dfetch(net_force_tex, idx_particle_index);
-            
-            body_force[ilocal].x = fi.x;
-            body_force[ilocal].y = fi.y;
-            body_force[ilocal].z = fi.z;
-            body_force[ilocal].w = fi.w;
-            
-                     
-            // This might require more calculations but more stable
-            ri.x = ex_space.x * particle_pos.x + ey_space.x * particle_pos.y 
-                    + ez_space.x * particle_pos.z;
-            ri.y = ex_space.y * particle_pos.x + ey_space.y * particle_pos.y 
-                    + ez_space.y * particle_pos.z;
-            ri.z = ex_space.z * particle_pos.x + ey_space.z * particle_pos.y 
-                    + ez_space.z * particle_pos.z;
-
-            torquei.x = ri.y * fi.z - ri.z * fi.y;
-            torquei.y = ri.z * fi.x - ri.x * fi.z;
-            torquei.z = ri.x * fi.y - ri.y * fi.x;
-            torquei.w = 0.0;
-            
-            body_torque[ilocal].x = torquei.x;
-            body_torque[ilocal].y = torquei.y;
-            body_torque[ilocal].z = torquei.z;
-            body_torque[ilocal].w = torquei.w;
-            }
-            
-        __syncthreads();
-
-        // reduction within the current window
-        unsigned int offset = window_size >> 1;
-        while (offset > 0)
-            {
-            if (threadIdx.x < offset)
-                {
-                body_force[ilocal].x += body_force[ilocal + offset].x;
-                body_force[ilocal].y += body_force[ilocal + offset].y;
-                body_force[ilocal].z += body_force[ilocal + offset].z;
-                body_force[ilocal].w += body_force[ilocal + offset].w;
-                
-                body_torque[ilocal].x += body_torque[ilocal + offset].x;
-                body_torque[ilocal].y += body_torque[ilocal + offset].y;
-                body_torque[ilocal].z += body_torque[ilocal + offset].z;
-                body_torque[ilocal].w += body_torque[ilocal + offset].w;
-                }
-                
-            offset >>= 1;
-            
-            __syncthreads();
-            }
-        
-        // accumulate the body force into the thread-local variables
-        force2.x += body_force[0].x;
-        force2.y += body_force[0].y;
-        force2.z += body_force[0].z;
-        force2.w += body_force[0].w;
-        
-        torque2.x += body_torque[0].x;
-        torque2.y += body_torque[0].y;
-        torque2.z += body_torque[0].z;
-        torque2.w += body_torque[0].w;
-        
-        __syncthreads();
-    }
-        
-    if (idx_body < n_bodies)
-        {
-        if (threadIdx.x == 0)
-            {
-            rdata_force[idx_body] = force2;
-            rdata_torque[idx_body] = torque2;
-            }
-        }
-*/
     int idx_body = blockIdx.x + local_beg;
     
     float4 *body_force = sum;

    \param n_bodies Number of rigid bodies
    \param local_beg Starting body index in this card
    \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
*/
extern "C" __global__ void gpu_rigid_force_kernel(float4* rdata_force, 

                                                                                             n_bodies, 
                                                                                             local_beg,
                                                                                             nmax,
                                                                                             box);
		}
	else	// large rigid bodies
