                          , 0);
     check_error("getting max shared mem size");
     unsigned int shared_mem = 2 * BSIZE_SCR * n_cols * sizeof(float);
-    #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,prev_max_state,\
-                   shared_mem,first_frame_above_threshold)\
-      shared(d_coords_sorted,d_clustering,\
-             tmp_coords_sorted,prev_clustering_sorted)\
-      num_threads(n_gpus)
+    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) {
-      // allocate memory and copy sorted coords to GPUs
       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);
-
-
-      //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
-               , cudaMemcpyHostToDevice);
-      // perform initial clustering on yet unclustered frames
-      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]
-         , i_from
-         , i_to);
-      }
-      cudaDeviceSynchronize();
-      check_error("after kernel loop");
-    }
-    // collect & merge clustering results from GPUs
-    std::vector<unsigned int> clustering_sorted = prev_clustering_sorted;
-    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) * first_frame_above_threshold
-               , cudaMemcpyDeviceToHost);
-      for (i=prev_last_frame; i < first_frame_above_threshold; ++i) {
-        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];
-    }
+    }
+    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

#include "logger.hpp"

#include <algorithm>

#include <cuda.h>
#include <omp.h>

    int comp_size = min(bsize, n_rows - offset);
    if (tid < comp_size) {
      for (unsigned int j=0; j < n_cols; ++j) {
        smem_coords[tid*n_cols+j] = sorted_coords[(tid+offset)*n_cols+j];
      }
      smem_ref_states[tid] = ref_states[(tid+offset)];
    }
    __syncthreads();
    if (gid < i_to) {
      unsigned int min_id = clustering[gid];
      // load reference coordinates for re-use into shared memory
      unsigned int ref_id = tid+bsize;
      for (unsigned int j=0; j < n_cols; ++j) {
        smem_coords[ref_id*n_cols+j] = sorted_coords[gid*n_cols+j];
      }
      // check against reference structures: if close enough, they'll be
      // considered as the same state (choosing min. state id).
      // state lumping will be done on host.
      for (unsigned int i=0; i < comp_size; ++i) {
        if (min_id != smem_ref_states[i]
         && smem_ref_states[i] != 0) {
          float dist2 = 0.0f;
          for (unsigned int j=0; j < n_cols; ++j) {
            float c = smem_coords[ref_id*n_cols+j] - smem_coords[i*n_cols+j];
            dist2 = fma(c, c, dist2);
          }
          if (dist2 < max_dist2
           && smem_ref_states[i] < min_id) {
            min_id = smem_ref_states[i];
          }
        }
      }
      clustering[gid] = min_id;
    }
  }

  std::vector<unsigned int>
  lumped_clusters(std::vector<unsigned int> clustering
                , std::vector<unsigned int> ref_states
                , unsigned int first_frame_above_threshold) {
    std::unordered_set<unsigned int> sstates;
    sstates.insert(clustering.begin()
                 , clustering.begin() + first_frame_above_threshold);
    sstates.insert(ref_states.begin()
                 , ref_states.begin() + first_frame_above_threshold);
    sstates.erase(0);
    std::vector<unsigned int> states(sstates.begin(), sstates.end());
    unsigned int n_states = states.size();
    // identify lumps
    unsigned int i, j, s, s1, s2;
    std::map<unsigned int, unsigned int> state_mapping;
    for (i=0; i < n_states; ++i) {
      state_mapping[states[i]] = states[i];
    }
    #pragma omp parallel for\
      default(none)\
      private(i,j,s1,s2,s)\
      firstprivate(n_states,first_frame_above_threshold)\
      shared(state_mapping,states,clustering,ref_states)
    for (i=0; i < n_states; ++i) {
      s = states[i];
      for (j=0; j < first_frame_above_threshold; ++j) {
        s1 = clustering[j];
        s2 = ref_states[j];
        //TODO: no loop over states, just over frames,
        //      then select 's' from {s1,s2}, compute min and update
        //      state_mapping in critical section
        if (s1 == s || s2 == s) {
          state_mapping[s] = std::min(state_mapping[s]
                                    , std::min(s1, s2));
        }
      }
    }
    // assign new (ordered) names to lumps
    std::set<unsigned int> distinct_states;
    for (auto sm: state_mapping) {
      distinct_states.insert(sm.second);
    }
    std::map<unsigned int, unsigned int> remapped_distinct;
    unsigned int i_state = 0;
    for (unsigned int s: distinct_states) {
      remapped_distinct[s] = ++i_state;
    }
    for (auto& sm: state_mapping) {
      sm.second = remapped_distinct[sm.second];
    }
    // remap states to lumps
    #pragma omp parallel for\
      default(none)\
      private(i)\
      firstprivate(first_frame_above_threshold)\
      shared(clustering,state_mapping)
    for (i=0; i < first_frame_above_threshold; ++i) {
      clustering[i] = state_mapping[clustering[i]];
    }
    return clustering;
  }

  std::vector<unsigned int>
  lumped_clusters__old(std::vector<unsigned int> clustering
                , std::vector<unsigned int> ref_states
                , unsigned int first_frame_above_threshold) {
    unsigned int i_state = 0;
    std::map<unsigned int, std::unordered_set<unsigned int>> lumps;
    std::set<std::pair<unsigned int, unsigned int>> pairings;
    for (unsigned int i=0; i < first_frame_above_threshold; ++i) {
      unsigned int s1 = clustering[i];
      unsigned int s2 = ref_states[i];
      pairings.emplace(s1, s2);
    }
    std::unordered_set<unsigned int> seen_states;
    for (auto p: pairings) {
      unsigned int s1 = p.first;
      unsigned int s2 = p.second;
      if (s1 == s2 && seen_states.count(s1) > 0) {
        continue;
      } else {
        seen_states.insert(s1);
        seen_states.insert(s2);
      }
      unsigned int l1=0, l2=0;
      // find lumps of states
      for (auto l: lumps) {
        if (l.second.count(s1) > 0) {
          l1 = l.first;
        }
        if (l.second.count(s2) > 0) {
          l2 = l.first;
        }
        if (l1 != 0 && l2 != 0) {
          break;
        }
      }
      if (l1 == 0 && l2 == 0) {
        // no lump found: make new one
        lumps[++i_state] = {s1, s2};
      } else if (l1 == 0 && l2 != 0) {
        // s2 has lump: add s1
        lumps[l2].insert(s1);
      } else if (l1 != 0 && l2 == 0) {
        // s1 has lump: add s2
        lumps[l1].insert(s2);
      } else if (l1 != l2) {
        // s1 and s2 have different lumps: join them
        auto mm = std::minmax(l1,l2);
        lumps[mm.first].insert(lumps[mm.second].begin()
                             , lumps[mm.second].end());
        lumps.erase(mm.second);
      }
    }
    // recluster trajectory
    i_state = 0;
    std::unordered_map<unsigned int, unsigned int> dict;
    for (auto l: lumps) {
      ++i_state;
      for (auto s: l.second) {
        dict[s] = i_state;
      }
    }
    for (unsigned int i=0; i < clustering.size(); ++i) {
      if (clustering[i] > 0) {
        clustering[i] = dict[clustering[i]];
      }
    }
    return clustering;
  }

  std::vector<std::size_t>
  initial_density_clustering(const std::vector<float>& free_energy
