     gpu_nve_rigid_step_two_body_kernel<<< body_grid, body_threads >>>(rigid_data.vel, 
                                                                       rigid_data.angmom, 
                                                                       rigid_data.angvel,
+                                                                      rigid_data.orientation,
                                                                       rigid_data.conjqm,
+                                                                      rigid_data.body_mass,
+                                                                      rigid_data.moment_inertia,
+                                                                      rigid_data.force,
+                                                                      rigid_data.torque,
+                                                                      rigid_data.body_indices,
                                                                       n_group_bodies,
                                                                       n_bodies, 
                                                                       nmax, 
                                                                       box, 
                                                                       deltaT);
+    
     block_size = 192;
     dim3 particle_grid(group_size/block_size+1, 1, 1);
     dim3 particle_threads(block_size, 1, 1);

    {
    unsigned int group_idx = blockIdx.x * blockDim.x + threadIdx.x;
    
    if (group_idx < n_group_bodies)
        {
        float body_mass;
        float4 moment_inertia, vel, angmom, orientation, ex_space, ey_space, ez_space, force, torque;
        float dt_half = 0.5 * deltaT;
        
        unsigned int idx_body = tex1Dfetch(rigid_data_body_indices_tex, group_idx);
            
        // Update body velocity and angmom
        if (idx_body < n_bodies)
            {        
            // update the velocity
            body_mass = tex1Dfetch(rigid_data_body_mass_tex, idx_body);
            vel = tex1Dfetch(rigid_data_vel_tex, idx_body);
            angmom = tex1Dfetch(rigid_data_angmom_tex, idx_body);
            force = tex1Dfetch(rigid_data_force_tex, idx_body);
            torque = tex1Dfetch(rigid_data_torque_tex, idx_body);
            moment_inertia = tex1Dfetch(rigid_data_moment_inertia_tex, idx_body);
            orientation = tex1Dfetch(rigid_data_orientation_tex, idx_body);
            
            exyzFromQuaternion(orientation, ex_space, ey_space, ez_space);
              
            float dtfm = dt_half / body_mass;
            float4 vel2;
            vel2.x = vel.x + dtfm * force.x;
            vel2.y = vel.y + dtfm * force.y;
            vel2.z = vel.z + dtfm * force.z;
            vel2.w = 0.0;

            // update angular momentum
            float4 angmom2;
            angmom2.x = angmom.x + dt_half * torque.x;
            angmom2.y = angmom.y + dt_half * torque.y;
            angmom2.z = angmom.z + dt_half * torque.z;
            angmom2.w = 0.0;
            
            // update angular velocity        
            float4 angvel2;
            computeAngularVelocity(angmom2, moment_inertia, ex_space, ey_space, ez_space, angvel2);
            
            // write out results
            rdata_vel[idx_body] = vel2;
            rdata_angmom[idx_body] = angmom2;
            rdata_angvel[idx_body] = angvel2;
            }
        }
    }

// Take the second 1/2 step forward in the NVE integration step

extern "C" __global__ void gpu_nve_rigid_step_two_body_kernel(float4* rdata_vel, 
                                                         float4* rdata_angmom, 
                                                         float4* rdata_angvel,
                                                         float4* rdata_orientation,
                                                         float4* rdata_conjqm,
                                                         float *d_rigid_mass,
                                                         float4 *d_rigid_mi,
                                                         float4 *d_rigid_force,
                                                         float4 *d_rigid_torque,
                                                         unsigned int *d_rigid_group,
                                                         unsigned int n_group_bodies,
                                                         unsigned int n_bodies, 
                                                         unsigned int nmax,
