 */
 extern "C" __global__ void gpu_rigid_force_sliding_kernel(float4* rdata_force, 
                                                  float4* rdata_torque,
+                                                 unsigned int *d_rigid_group,
+                                                 float4* d_rigid_orientation,
+                                                 unsigned int* d_rigid_particle_idx,
+                                                 float4* d_rigid_particle_dis,
+                                                 float4* d_net_force,
                                                  unsigned int n_group_bodies, 
                                                  unsigned int n_bodies, 
                                                  unsigned int nmax,

    {
    int group_idx = blockIdx.x;
    
    float4 *body_force = sum;
    float4 *body_torque = &sum[window_size];
        
    __shared__ float4 ex_space, ey_space, ez_space;
    float4 force2 = make_float4(0.0f, 0.0f, 0.0f, 0.0f); 
    float4 torque2 = make_float4(0.0f, 0.0f, 0.0f, 0.0f); 
    
    int idx_body = -1;
            
    if (group_idx < n_group_bodies)
        {
        idx_body = tex1Dfetch(rigid_data_body_indices_tex, group_idx);
                  
        if (idx_body < n_bodies && threadIdx.x == 0)
            {
            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);
            }
        }
                
    __syncthreads();
        
    unsigned int n_windows = nmax / window_size;
        
    // slide the window throughout the block
    for (unsigned int start = 0; start < n_windows; start++)
        {
        body_force[threadIdx.x] = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
        body_torque[threadIdx.x] = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
        float4 fi = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
        float4 torquei = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
        
        if (idx_body >= 0 && idx_body < n_bodies)
            { 
            int idx_particle = idx_body * nmax + start * window_size + threadIdx.x;
            unsigned int idx_particle_index = tex1Dfetch(rigid_data_particle_indices_tex, idx_particle);
                     
            if (idx_particle_index != INVALID_INDEX)
                {
                // calculate body force and torques
                float4 particle_pos = tex1Dfetch(rigid_data_particle_pos_tex, idx_particle);
                fi = tex1Dfetch(net_force_tex, idx_particle_index);
                                
                body_force[threadIdx.x].x = fi.x;
                body_force[threadIdx.x].y = fi.y;
                body_force[threadIdx.x].z = fi.z;
                body_force[threadIdx.x].w = fi.w;
                
                // This might require more calculations but more stable
                // particularly when rigid bodies are bigger than half the box
                float4 ri;
                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[threadIdx.x].x = torquei.x;
                body_torque[threadIdx.x].y = torquei.y;
                body_torque[threadIdx.x].z = torquei.z;
                body_torque[threadIdx.x].w = torquei.w;
                }
            }
            
        __syncthreads();

        // reduction within the current window
        unsigned int offset = window_size >> 1;
        while (offset > 0)
            {
            if (threadIdx.x < offset)
                {
                body_force[threadIdx.x].x += body_force[threadIdx.x + offset].x;
                body_force[threadIdx.x].y += body_force[threadIdx.x + offset].y;
                body_force[threadIdx.x].z += body_force[threadIdx.x + offset].z;
                body_force[threadIdx.x].w += body_force[threadIdx.x + offset].w;
                
                body_torque[threadIdx.x].x += body_torque[threadIdx.x + offset].x;
                body_torque[threadIdx.x].y += body_torque[threadIdx.x + offset].y;
                body_torque[threadIdx.x].z += body_torque[threadIdx.x + offset].z;
                body_torque[threadIdx.x].w += body_torque[threadIdx.x + 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 >= 0 && idx_body < n_bodies)
        {
        if (threadIdx.x == 0)
            {
            rdata_force[idx_body] = force2;
            rdata_torque[idx_body] = torque2;
            }
        }
    
			
    }

/*! \param pdata Particle data to step forward 1/2 step
    \param rigid_data Rigid body data to step forward 1/2 step

#pragma mark RIGID_FORCE_KERNEL

//! Shared memory for body force and torque reduction, required allocation when the kernel is called
extern __shared__ float3 sum[];

//! Calculates the body forces and torques by summing the constituent particle forces using a fixed sliding window size
/*! \param rdata_force Body forces
