         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

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

#include "CellCommunicator.cuh"
#include "ReductionOperators.h"

#include "hoomd/extern/cub/cub/cub.cuh"

namespace mpcd
{
namespace gpu
{

//! CTA reduce, returns result in first thread
template<int LOGICAL_WARP_SIZE, typename T>
__device__ static T warp_reduce(T x)
    {
    static_assert(LOGICAL_WARP_SIZE && !(LOGICAL_WARP_SIZE & (LOGICAL_WARP_SIZE-1)), "Logical warp size must be a power of 2");

    T tmp = x;
    #pragma unroll
    for (int dest_count = LOGICAL_WARP_SIZE/2; dest_count >= 1; dest_count /= 2)
        {
        tmp += cub::ShuffleDown(x, dest_count, LOGICAL_WARP_SIZE);
        }
    return tmp;
    }

namespace kernel
{
//! Begins the cell thermo compute by summing cell quantities
