 
 global void reduce_energy(force_energy_tup *result) {
 	int l = threadIdx.x;
-	shared float energies[warpSize];
-	float e = block_sum<float>(energies, result[l].energy);
+	float e = block_sum<float>(result[l].energy);
 	if ( l == 0 ) {
 		result[0].energy = e;
 	}	

		rec_deriv = diff * dor;
	}
    
    shared float energies[32];
	shared float3 derivs[32];
	float this_e = block_sum<float>(energies, rec_energy);
	float3 deriv = block_sum<float3>(derivs, rec_deriv);
	if (threadIdx.x == 0)
	{
		curl(this_e, (float *) &deriv, v);

    a.energy += b.energy;
}

/* requires blockDim.x <= 1024, blockDim.y == 1 */
template <class T>
device __forceinline__
T block_sum(T mySum) {
	const unsigned int lane = threadIdx.x & 31;
	const unsigned int wid = threadIdx.x>>5;

    __shared__ T scratch[32];

	mySum = warp_sum(mySum);
	if (lane==0)
        scratch[wid] = mySum;
	__syncthreads();

	if (wid == 0) {
		mySum = (threadIdx.x < blockDim.x >> 5) ? scratch[lane] : zero<T>();
		mySum = warp_sum(mySum);
        if (threadIdx.x == 0 && isNotDiv32(blockDim.x))
            mySum += scratch[blockDim.x >> 5];
	}
	return mySum;
}
