         ke += (double)(0.5) * mass_i * (vel_i.x * vel_i.x + vel_i.y * vel_i.y + vel_i.z * vel_i.z);
         }
 
-    d_cell_vel[idx] = momentum;
-    d_cell_energy[idx] = make_scalar3(ke, 0.0, __int_as_scalar(np));
+    if (tpp > 1)
+        {
+        momentum.x = warp_reduce<tpp>(momentum.x);
+        momentum.y = warp_reduce<tpp>(momentum.y);
+        momentum.z = warp_reduce<tpp>(momentum.z);
+        momentum.w = warp_reduce<tpp>(momentum.w);
+        ke = warp_reduce<tpp>(ke);
+        }
+
+    if (tpp == 1 || (idx % tpp == 0))
+        {
+        d_cell_vel[cell_id] = momentum;
+        d_cell_energy[cell_id] = make_scalar3(ke, 0.0, __int_as_scalar(np));
+        }
     }
 
 //! Finalizes the cell thermo compute by properly averaging cell quantities

    if (max_block_size == UINT_MAX)
        {
        cudaFuncAttributes attr;
        cudaFuncGetAttributes(&attr, (const void*)mpcd::gpu::kernel::begin_cell_thermo);
        max_block_size = attr.maxThreadsPerBlock;
        }

    unsigned int run_block_size = min(block_size, max_block_size);
    dim3 grid(cli.getH() / run_block_size + 1);
    mpcd::gpu::kernel::begin_cell_thermo<<<grid, run_block_size>>>(d_cell_vel,
                                                                   d_cell_energy,
                                                                   d_cell_np,
                                                                   d_cell_list,

 * \todo The cell accumulation could have wider parallelism by using multiple
 *       threads in a strided access pattern to compute the cell properties.
 */
template<unsigned int tpp>
__global__ void begin_cell_thermo(Scalar4 *d_cell_vel,
                                  Scalar3 *d_cell_energy,
                                  const unsigned int *d_cell_np,
