                                 , unsigned int n_cols
                                 , float max_dist2
                                 , unsigned int* clustering
+                                , unsigned int* ref_states
                                 , unsigned int i_from
                                 , unsigned int i_to) {
-    extern __shared__ float smem[];
+    // dynamic shared mem for ref coords
+    extern __shared__ float smem_coords[];
+    // static shared mem for cluster ids
+    __shared__ unsigned int smem_ref_states[BSIZE_SCR];
+    // thread dimensions
     unsigned int bid = blockIdx.x;
     unsigned int tid = threadIdx.x;
     unsigned int bsize = blockDim.x;

                                                        , n_rows
                                                        , initial_clusters);
    // write log
//    screening_log(sigma2
//                , first_frame_above_threshold
//                , fe_sorted);
    float max_dist2 = 4*sigma2;
    // prepare CUDA environment
    int n_gpus = get_num_gpus();
    std::vector<float*> d_coords_sorted(n_gpus);
    std::vector<unsigned int*> d_clustering(n_gpus);
    // sort coords (and previous clustering results)
    // according to free energies
    std::vector<float> tmp_coords_sorted(n_rows * n_cols);
    std::vector<unsigned int> prev_clustering_sorted(n_rows);
    for (unsigned int i=0; i < n_rows; ++i) {
      for (unsigned int j=0; j < n_cols; ++j) {
        tmp_coords_sorted[i*n_cols+j] = coords[fe_sorted[i].first*n_cols+j];
      }
      prev_clustering_sorted[i] = prev_clustering[fe_sorted[i].first];
    }
    unsigned int block_rng;
//TODO: check: what happens if no new frames are added compared to prev run?
    unsigned int i_from;
    unsigned int i_to;
    unsigned int i;
    // re-use initial clustering results
    unsigned int prev_last_frame = 0;
    if (prev_clustering_sorted[0] != 0) {
      for (i=1; i < n_rows; ++i) {
        if (prev_clustering_sorted[i] == 0) {
          prev_last_frame = i;
          break;
        }
      }
    }
    int i_gpu;
    unsigned int gpu_rng =
      min_multiplicator(first_frame_above_threshold - prev_last_frame
                      , n_gpus);

                         , 0);
    check_error("getting max shared mem size");
    unsigned int shared_mem = 2 * BSIZE_SCR * n_cols * sizeof(float);
    unsigned int block_rng, i_from, i_to, i, i_gpu;
    // lump microstates until nothing changes
    bool microstates_lumped = true;
    unsigned int i_loop = 0;
    // initialize GPUs
    for (i_gpu=0; i_gpu < n_gpus; ++i_gpu) {
      cudaSetDevice(i_gpu);
      // allocate memory on GPUs
      cudaMalloc((void**) &d_coords_sorted[i_gpu]
               , sizeof(float) * n_rows * n_cols);
      cudaMalloc((void**) &d_clustering[i_gpu]
               , sizeof(unsigned int) * n_rows);
      cudaMalloc((void**) &d_ref_states[i_gpu]
               , sizeof(unsigned int) * n_rows);
      // copy sorted coords and previous clustering results to GPUs
      cudaMemcpy(d_coords_sorted[i_gpu]
               , tmp_coords_sorted.data()
               , sizeof(float) * n_rows * n_cols
               , cudaMemcpyHostToDevice);
    }
    while (microstates_lumped) {
      std::vector<unsigned int> clustering_sorted_orig = clustering_sorted;
      std::cerr << "microstate lumping iteration " << ++i_loop << std::endl;
      #pragma omp parallel for\
        default(none)\
        private(i,i_gpu,block_rng,i_from,i_to)\
        firstprivate(n_gpus,n_rows,n_cols,gpu_rng,max_dist2,\
                     prev_last_frame,\
                     shared_mem,first_frame_above_threshold)\
        shared(d_coords_sorted,d_clustering,d_ref_states,\
               tmp_coords_sorted,clustering_sorted)\
        num_threads(n_gpus)
      for (i_gpu=0; i_gpu < n_gpus; ++i_gpu) {
        cudaSetDevice(i_gpu);
        cudaMemcpy(d_clustering[i_gpu]
                 , clustering_sorted.data()
                 , sizeof(unsigned int) * n_rows
                 , cudaMemcpyHostToDevice);
        cudaMemcpy(d_ref_states[i_gpu]
                 , clustering_sorted.data()
                 , sizeof(unsigned int) * n_rows
                 , cudaMemcpyHostToDevice);
        i_from = prev_last_frame + i_gpu * gpu_rng;
        i_to = (i_gpu == (n_gpus-1))
             ? first_frame_above_threshold
             : prev_last_frame + (i_gpu+1) * gpu_rng;
        block_rng = min_multiplicator(i_to-i_from
                                    , BSIZE_SCR);
        for (i=0; i*BSIZE_SCR < first_frame_above_threshold; ++i) {
          initial_density_clustering_krnl
            <<< block_rng
              , BSIZE_SCR
              , shared_mem >>>
            (i*BSIZE_SCR
           , d_coords_sorted[i_gpu]
           , n_rows
           , n_cols
           , max_dist2
           , d_clustering[i_gpu]
           , d_ref_states[i_gpu]
           , i_from
           , i_to);
        }
        cudaDeviceSynchronize();
        check_error("after kernel loop");
      }
      // collect & merge clustering results from GPUs
      for (int i_gpu=0; i_gpu < n_gpus; ++i_gpu) {
        std::vector<unsigned int> tmp_clust(n_rows, 0);
        cudaMemcpy(tmp_clust.data()
                 , d_clustering[i_gpu]
                 , sizeof(unsigned int) * n_rows
                 , cudaMemcpyDeviceToHost);
        for (i=0; i < first_frame_above_threshold; ++i) {
          if (i_gpu == 0) {
            clustering_sorted[i] = tmp_clust[i];
          } else {
            clustering_sorted[i] = std::min(clustering_sorted[i]
                                          , tmp_clust[i]);
          }
        }
      }
      std::cout << "lumping CUDA-results ..." << std::endl;
      // lump microstates
      clustering_sorted = lumped_clusters(clustering_sorted
                                        , clustering_sorted_orig
                                        , first_frame_above_threshold);
      std::cout << "     ... finished" << std::endl;
      // compare prev_clustering to clustering_sorted:
      // if equal, end while
      microstates_lumped = (clustering_sorted_orig != clustering_sorted);
    } // end while (if microstates_lumped == false)
    // cleanup CUDA environment
    for (int i_gpu=0; i_gpu < n_gpus; ++i_gpu) {
      cudaFree(d_coords_sorted[i_gpu]);
      cudaFree(d_clustering[i_gpu]);
      cudaFree(d_ref_states[i_gpu]);
    }
    // convert state trajectory from
    // FE-sorted order to original order
    std::vector<std::size_t> clustering(n_rows);
    for (unsigned int i=0; i < n_rows; ++i) {
      clustering[fe_sorted[i].first] = clustering_sorted[i];
    }
    return normalized_cluster_names(first_frame_above_threshold
                                  , clustering
