     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>());

  // TODO: fix momentum at step 0
  inc = blockDim.x * gridDim.x;
  for (i = threadIdx.x + blockIdx.x * blockDim.x; i < N; i += inc) {
      tmpx = 4.0f * (attr_forces[i] - (rep_forces[i] / norm));
      tmpy = 4.0f * (attr_forces[i + N] - (rep_forces[nnodes + 1 + i] / norm));
      tmpx = momentum * tmpx + (1 - momentum) * old_forces[i];
      tmpy = momentum * tmpy + (1 - momentum) * old_forces[i + N];
      pts[i] -= eta * tmpx;
