                                                         , nh
                                                         , free_energy_threshold
                                                         , n_rows
-                                                        , {});
-//TODO: debug
-                                              //          , initial_clusters);
+                                                        , initial_clusters);
     // write log
 //    screening_log(sigma2
 //                , first_frame_above_threshold

      cudaFree(d_coords_sorted[i_gpu]);
      cudaFree(d_clustering[i_gpu]);
    }

//TODO
    std::cout << "#  " << clustering_sorted[0] << std::endl;

//    return clustering;
    return normalized_cluster_names(first_frame_above_threshold
                                  , clustering
                                  , fe_sorted);

        if (i_gpu == 0) {
          clustering_sorted[i] = tmp_clust[i];
        } else {
          clustering_sorted[i] = std::max(clustering_sorted[i]
                                        , tmp_clust[i]);
        }
      }
    }
    // reduce clustering to min number of ids
    clustering_sorted = sanitize_state_names(clustering_sorted);
    //TODO: refined clustering (by rerunning with alternative kernel)




    // convert state trajectory from
    // FE-sorted order to original order
    std::vector<std::size_t> clustering(n_rows, 0);
    for (unsigned int i=0; i < n_rows; ++i) {
      clustering[fe_sorted[i].first] = clustering_sorted[i];
    }
