     #pragma omp parallel for default(none)\
       private(i)\
       firstprivate(n_gpus,n_rows,n_cols,gpu_range)\
-      shared(partial_pops,radii,coords,sorted_coords,blimits)\
+      shared(partial_pops,radii,coords)\
       num_threads(n_gpus)\
       schedule(dynamic,1)
     for (i=0; i < n_gpus; ++i) {


// for pops
#define BSIZE_POPS 64
//#define N_STREAMS 8

// for neighborhood search
#define BSIZE_NH 128

    using Clustering::Tools::min_multiplicator;
    ASSUME_ALIGNED(coords);
    unsigned int n_radii = radii.size();
    std::vector<float> rad2(n_radii);
    for (std::size_t i=0; i < n_radii; ++i) {
      rad2[i] = radii[i]*radii[i];
    }
    // GPU setup
    cudaSetDevice(i_gpu);
    float* d_coords;
    float* d_rad2;
    unsigned int* d_pops;
    cudaMalloc((void**) &d_coords
             , sizeof(float) * n_rows * n_cols);
    cudaMalloc((void**) &d_pops
             , sizeof(unsigned int) * n_rows * n_radii);
    cudaMalloc((void**) &d_rad2
             , sizeof(float) * n_radii);
    cudaMemset(d_pops
             , 0
