 	float4 vel2 = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
 	float4 particle_accel = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
 	
-	// ri is used throughout the kernel
+	// ri, body_mass and moment_inertia is used throughout the kernel
 	if (idx_body < rigid_data.n_bodies)
 		{
 		ex_space = tex1Dfetch(rigid_data_exspace_tex, idx_body);

		advanceQuaternion(angmom2, moment_inertia, angvel, ex_space, ey_space, ez_space, orientation, deltaT); 
			
		// write out the results (MEM_TRANSFER: ? bytes)
		rigid_data.com[idx_body] = pos2;
		rigid_data.vel[idx_body] = vel2;
		rigid_data.angmom[idx_body] = angmom2;
		rigid_data.angvel[idx_body] = angvel;
		rigid_data.orientation[idx_body] = orientation;
		rigid_data.ex_space[idx_body] = ex_space;
		rigid_data.ey_space[idx_body] = ey_space;
		rigid_data.ez_space[idx_body] = ez_space;
		rigid_data.body_imagex[idx_body] = body_imagex;
		rigid_data.body_imagey[idx_body] = body_imagey;
		rigid_data.body_imagez[idx_body] = body_imagez;
		}
	}

#pragma mark PARTICLE_PRE_STEP_KERNEL

extern "C" __global__ void gpu_nve_rigid_particle_pre_step_kernel(gpu_pdata_arrays pdata, gpu_rigid_data_arrays rigid_data, gpu_boxsize box, float deltaT, bool limit, float limit_val)
	{
	unsigned int idx_particle = blockIdx.x * blockDim.x + threadIdx.x;	
	unsigned int idx_body = blockIdx.x + rigid_data.local_beg; 	

	unsigned int idx_particle_index = tex1Dfetch(rigid_data_particle_indices_tex, idx_particle);
	// Since we use nmax for all rigid bodies, there might be some empty slot for particles in a rigid body
	// the particle index of these empty slots is set to be INVALID_INDEX.
	if (idx_body < rigid_data.n_bodies && idx_particle_index != INVALID_INDEX) 
		{
		float4 com = tex1Dfetch(rigid_data_com_tex, idx_body);
		float4 vel = tex1Dfetch(rigid_data_vel_tex, idx_body);
		float4 angvel = tex1Dfetch(rigid_data_angvel_tex, idx_body);
		float4 ex_space = tex1Dfetch(rigid_data_exspace_tex, idx_body);
		float4 ey_space = tex1Dfetch(rigid_data_eyspace_tex, idx_body);
		float4 ez_space = tex1Dfetch(rigid_data_ezspace_tex, idx_body);
		float4 particle_pos = tex1Dfetch(rigid_data_particle_pos_tex, idx_particle);
		
		// project the position in the body frame to the space frame: ri = rotation_matrix * particle_pos
		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;
		
		// time to fix the periodic boundary conditions (FLOPS: 15)
		int4 image = tex1Dfetch(pdata_image_tex, idx_particle_index);
		float4 pos = tex1Dfetch(pdata_pos_tex, idx_particle_index);
		
		// x_particle = com + ri
		float4 pos2;
		pos2.x = com.x + ri.x;
		pos2.y = com.y + ri.y;
		pos2.z = com.z + ri.z;
		pos2.w = pos.w;
		
		float x_shift = rintf(pos2.x * box.Lxinv);
		pos2.x -= box.Lx * x_shift;
		image.x += (int)x_shift;
		
		float y_shift = rintf(pos2.y * box.Lyinv);
		pos2.y -= box.Ly * y_shift;
		image.y += (int)y_shift;
		
		float z_shift = rintf(pos2.z * box.Lzinv);
		pos2.z -= box.Lz * z_shift;
		image.z += (int)z_shift;
		
		// v_particle = vel + angvel x ri
		float4 particle_vel;
		particle_vel.x = vel.x + angvel.y * ri.z - angvel.z * ri.y;
		particle_vel.y = vel.y + angvel.z * ri.x - angvel.x * ri.z;
		particle_vel.z = vel.z + angvel.x * ri.y - angvel.y * ri.x;
		
		// write out the results (MEM_TRANSFER: ? bytes)
		pdata.pos[idx_particle_index] = pos2;
		pdata.vel[idx_particle_index] = particle_vel;
		pdata.image[idx_particle_index] = image;
		}

	}	


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

	if (error != cudaSuccess)
		return error;

	// run the kernel for bodies
    gpu_nve_rigid_body_pre_step_kernel<<< grid, threads >>>(pdata, rigid_data, box, deltaT, limit, limit_val);

		
	if (!g_gpu_error_checking)
		{
