     \param x Point to take the expansion
 
 */
-__device__ float nph_maclaurin_series(float x)
-    {
-    float x2, x4;
+__device__ Scalar nph_maclaurin_series(Scalar x)
+    {
+    Scalar 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);
+    return (Scalar(1.0) + Scalar(1.0/6.0) * x2 + Scalar(1.0/120.0) * x4 + Scalar(1.0/5040.0) * x2 * x4 + Scalar(1.0/362880.0) * 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, 
+extern "C" __global__ void gpu_nph_rigid_zero_virial_rigid_kernel(Scalar *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;
+        d_virial_rigid[idx] = 0;
         }
 
     }

    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(float) >>>(npt_rdata.partial_Ksum_t,
                                                                                          npt_rdata.partial_Ksum_r,
                                                                                          npt_rdata.Ksum_t,
                                                                                          npt_rdata.Ksum_r,

    invquatvec(orientation, conjqm2, mbody);
    transpose_dot(ex_space, ey_space, ez_space, mbody, angmom2);
    
    angmom2.x *= Scalar(0.5);
    angmom2.y *= Scalar(0.5);
    angmom2.z *= Scalar(0.5);
    angmom2.w = 0;
    
    // update angular velocity
    Scalar4 angvel2;
    computeAngularVelocity(angmom2, moment_inertia, ex_space, ey_space, ez_space, angvel2);
    
    akin_r = angmom2.x * angvel2.x + angmom2.y * angvel2.y + angmom2.z * angvel2.z;
