     float norm;
     
     // These variables currently govern the tolerance (whether it recurses on a cell)
-    float theta = 0.5f;
+    float theta = 0.0f;
     float epssq = 0.05 * 0.05;
     // float itolsq = 1.0f / (0.5 * 0.5);
 

        gpuErrchk(cudaDeviceSynchronize());

        // compute attractive forces
        computeAttrForce(N_POINTS, sparsePij.size(), nnodes, sparse_handle, descr, sparsePij, pijRowPtr, pijColInd, forceProd, pts, attr_forces, ones);
        gpuErrchk(cudaDeviceSynchronize());
        
        norm = thrust::reduce(norml.begin(), norml.begin() + N_POINTS, 0.0f, thrust::plus<float>());

              // from sptree.cpp
              tmp = 1 / (1 + tmp);
              mult = massd[n] * tmp;
              // printf("%0.2f\n", massd[n]);
              normsum += mult;
              mult *= tmp;
              vx += dx * mult;
