                                                                       thrust::raw_pointer_cast(pts.data()),
                                                                       thrust::raw_pointer_cast(attr_forces.data()),
                                                                       thrust::raw_pointer_cast(rep_forces.data()),
+                                                                      thrust::raw_pointer_cast(gains.data()),
                                                                       thrust::raw_pointer_cast(old_forces.data()));
           gpuErrchk(cudaDeviceSynchronize());
           thrust::transform(pts.begin(), pts.end(), random_vec.begin(), pts.begin(), thrust::plus<float>());

        thrust::device_vector<float> forceProd(sparsePij.size());
        thrust::device_vector<float> rep_forces((nnodes + 1) * 2, 0);
        thrust::device_vector<float> attr_forces(N_POINTS * 2, 0);
        thrust::device_vector<float> old_forces(N_POINTS * 2, 0); // for momentum
        thrust::device_vector<int> errl(1);
        thrust::device_vector<int> startl(nnodes + 1);

            y += dy;
            ch = childd[n*4+j];
            // repeat until the two bodies are different children
            // if (depth >= 100)
            //     printf("parent: %d, child: %d, depth: %d, pos_parent: (%0.5f, %0.5f), pos_child: (%0.5f, %0.5f), dir: %d\n", n, ch, cell, depth, x, y, px, py, j);
          } while (ch >= 0 && r > 1e-10);
          childd[n*4+j] = i;

