           if (n >= 0) {
             dx = posxd[n] - px;
             dy = posyd[n] - py;
-            tmp = dx*dx + (dy*dy + epssqd); // distance squared + softening (how does the softening actually interact with things?
+            tmp = dx*dx + dy*dy; // distance squared
             // tmp = dx*dx + (dy*dy + epssqd) (why softening?)
             if ((n < nbodiesd) || __all(tmp >= dq[depth])) {  // check if all threads agree that cell is far enough away (or is a body)
             //   tmp = rsqrtf(tmp);  // compute distance

  __shared__ float dq[MAXDEPTH * THREADS5/WARPSIZE];

  if (0 == threadIdx.x) {
    tmp = radiusd * 2;
    // precompute values that depend only on tree level
    dq[0] = tmp * tmp * itolsqd;
    for (i = 1; i < maxdepthd; i++) {
      dq[i] = dq[i - 1] * 0.25f;
      dq[i - 1] += epssqd;
    }
    dq[i - 1] += epssqd;

    if (maxdepthd > MAXDEPTH) {
      *errd = maxdepthd;

    float norm;
    
    // These variables currently govern the tolerance (whether it recurses on a cell)
    float theta = 0.5f;
    InitializationKernel<<<1, 1>>>(thrust::raw_pointer_cast(errl.data()));
    gpuErrchk(cudaDeviceSynchronize());
    
