           dist2 = fma(c, c, dist2);
         }
         if (dist2 < max_dist2) {
-        //TODO: this must be the error: always taking the smallest id
-        //      may lead to self-loops. we must keep track of all(!) seen
-        //      neighbors...
-          cluster_id = min(cluster_id, i+offset+1);
+          unsigned int min_id = min(cluster_id, i+offset+1);
+          if (min_id == self_id) {
+            cluster_id = max(cluster_id, i+offset+1);
+          } else {
+            cluster_id = min_id;
+          }
         }
       }
       clustering[gid] = cluster_id;
     }
   }
 
-  unsigned int
-  follow(const std::vector<unsigned int>& clustering
-       , unsigned int id) {
-    if (id == 0
-    ||  clustering[id-1] == id) {
-      return id;
-    } else {
-      return follow(clustering
-                  , clustering[id-1]);
-    }
-  }
+//  unsigned int
+//  follow(const std::vector<unsigned int>& clustering
+//       , unsigned int id
+//       , unsigned int min_id = 0) {
+//    if (id == 0
+//     || id == min_id) {
+//      return id;
+//    } else {
+//      return follow(clustering
+//                  , clustering[id-1]
+//                  , (min_id == 0)
+//                      ? id
+//                      : std::min(id, min_id));
+//    }
+//  }
 
   std::vector<unsigned int>
-  sanitize_state_names(std::vector<unsigned int> clustering
-                     , unsigned int prev_max_state) {
-//    std::vector<unsigned int> unique_names =
-//      Clustering::Tools::unique_elements(clustering);
-//    std::map<unsigned int, unsigned int> end_points;
-//    for (unsigned int u: unique_names) {
-//      unsigned int u_orig = u;
-//      if (u > 0) {
-//        // state ids are strictly ordered,
-//        // thus this will always terminate
-//        while (u > prev_max_state && clustering[u-1] != u) {
-//          u = clustering[u-1];
-//        }
-//      }
-//      end_points[u_orig] = u;
+  sanitize_state_names(std::vector<unsigned int> clustering) {
+
+    //TODO: lump frames into microstates.
+    //      be careful: cluster ids do not necessarily
+    //                  denote frame ids (especially not with previous
+    //                  clustering)
+
+    // map names to ascending index
+    unsigned int id=0;
+    std::map<unsigned int, unsigned int> name_dict;
+    for (unsigned int name: unique_names) {
+      ++id;
+      name_dict[name] = id;
+    }
+    // relabel states
+    for (unsigned int& i: clustering) {
+      i = name_dict[i];
+    }
+    return clustering;
+  }
+
+
+//TODO remove old code
+//    for (unsigned int i=0; i < clustering.size(); ++i) {
+//      clustering[i] = follow(clustering
+//                           , clustering[i]);
 //    }
-//    for (unsigned int& s: clustering) {
-//      s = end_points[s];
-//    }
-
-//TODO test alternative
-    for (unsigned int i=0; i < clustering.size(); ++i) {
-      clustering[i] = follow(clustering
-                           , clustering[i]);
-    }
-    return clustering;
-  }
+//    return clustering;
+//  }
 
 
   std::vector<std::size_t>

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

               , tmp_coords_sorted.data()
               , sizeof(float) * n_rows * n_cols
               , cudaMemcpyHostToDevice);


      //TODO change kernel to use two lists: own state + min next-neighbor state


      // copy prev results to GPU-buffer (and set the rest to zero)
      cudaMemset(d_clustering[i_gpu]
               , 0
               , sizeof(unsigned int) * n_rows);
      cudaMemcpy(d_clustering[i_gpu]
               , prev_clustering_sorted.data()
               , sizeof(unsigned int) * n_rows
