     
 */
 
-extern "C" __global__ void gpu_nph_rigid_step_one_body_kernel(float4* rdata_com, 
-                                                            float4* rdata_vel, 
-                                                            float4* rdata_angmom, 
-                                                            float4* rdata_angvel,
-                                                            float4* rdata_orientation, 
+extern "C" __global__ void gpu_nph_rigid_step_one_body_kernel(Scalar4* rdata_com, 
+                                                            Scalar4* rdata_vel, 
+                                                            Scalar4* rdata_angmom, 
+                                                            Scalar4* rdata_angvel,
+                                                            Scalar4* rdata_orientation, 
                                                             int3* rdata_body_image, 
-                                                            float4* rdata_conjqm,
-                                                            float *d_rigid_mass,
-                                                            float4 *d_rigid_mi,
-                                                            float4 *d_rigid_force,
-                                                            float4 *d_rigid_torque,
+                                                            Scalar4* rdata_conjqm,
+                                                            Scalar *d_rigid_mass,
+                                                            Scalar4 *d_rigid_mi,
+                                                            Scalar4 *d_rigid_force,
+                                                            Scalar4 *d_rigid_torque,
                                                             unsigned int *d_rigid_group, 
                                                             unsigned int n_group_bodies,  
                                                             unsigned int n_bodies, 
-                                                            float npt_rdata_epsilon_dot, 
-                                                            float* npt_rdata_partial_Ksum_t, 
-                                                            float* npt_rdata_partial_Ksum_r,
+                                                            Scalar npt_rdata_epsilon_dot, 
+                                                            Scalar* npt_rdata_partial_Ksum_t, 
+                                                            Scalar* npt_rdata_partial_Ksum_r,
                                                             unsigned int npt_rdata_nf_t,
                                                             unsigned int npt_rdata_nf_r,
                                                             unsigned int npt_rdata_dimension, 
                                                             BoxDim box, 
-                                                            float deltaT)
+                                                            Scalar deltaT)
     {
     unsigned int group_idx = blockIdx.x * blockDim.x + threadIdx.x;
 

    \param x Point to take the expansion

*/
__device__ float nph_maclaurin_series(float x)
    {
    float x2, x4;
    x2 = x * x;
    x4 = x2 * x2;
    return (1.0f + (1.0f/6.0f) * x2 + (1.0f/120.0f) * x4 + (1.0f/5040.0f) * x2 * x4 + (1.0f/362880.0f) * x4 * x4);
    }

/*! Kernel to zero virial contribution from particles from rigid bodies
    \param d_virial_rigid Virial contribution from particles in rigid bodies
    \param local_num Number of particles in this card
*/
extern "C" __global__ void gpu_nph_rigid_zero_virial_rigid_kernel(float *d_virial_rigid, 
                                                                 unsigned int local_num)
    {
    int idx = blockIdx.x * blockDim.x + threadIdx.x; // particle's index

    if (idx < local_num)
        {
        d_virial_rigid[idx] = 0.0f;
        }

    }

    dim3 threads(block_size, 1, 1);
    
    // run the kernel: double the block size to accomodate Ksum_t and Ksum_r
    gpu_nph_rigid_reduce_ksum_kernel<<< grid, threads, 2 * block_size * sizeof(Scalar) >>>(npt_rdata.partial_Ksum_t,
                                                                                          npt_rdata.partial_Ksum_r,
                                                                                          npt_rdata.Ksum_t,
                                                                                          npt_rdata.Ksum_r,
