diff --git a/modelling/doc/algorithms.rst b/modelling/doc/algorithms.rst
index 439fd2d9009650a24933a6ccaebb0a6d567c8528..7ff5ae3a33e08ccccdf6d98f87abe2748e1d0d19 100644
--- a/modelling/doc/algorithms.rst
+++ b/modelling/doc/algorithms.rst
@@ -24,11 +24,9 @@ iterative superposition multiple times by using a sliding window to select the
 initial subset and gathers all unique results. These results can be very 
 similar and only differ by single positions. The algorithm therefore reduces
 the amount of solutions by merging them based on a threshold of similarity.
-If the sum of matching positions within the distance threshold divided by
-the maximum length of the two solutions is above a cluster thresh, the two
-solutions get merged by producing a common solution containing the shared
-positions. As a final result, the algorithm therefore detects common rigid
-subsets of positions.
+The similarity is defined by the fraction of positions in solution A that are
+also present in solution B. As a final result, the algorithm therefore detects 
+common rigid subsets of positions.
 
 .. method:: RigidBlocks(bb_list_one, bb_list_two, [window_length = 12, max_iterations=20, distance_thresh=3.0, cluster_thresh=0.9])
 
diff --git a/modelling/src/rigid_blocks.cc b/modelling/src/rigid_blocks.cc
index ef24abeb7a846f4039e1fd5331a5fc13b954691b..df46c90fc4195cf3c47b6f4bb09a2b9ec5301e7f 100644
--- a/modelling/src/rigid_blocks.cc
+++ b/modelling/src/rigid_blocks.cc
@@ -30,7 +30,7 @@ inline Real IndexListSimilarity(const std::vector<uint>& v1,
     sum += complete_v1[i] * complete_v2[i]; 
   }
 
-  return static_cast<Real>(sum) / std::max(v1.size(),v2.size());
+  return static_cast<Real>(sum) / std::min(v1.size(),v2.size());
 }
 
 void DoIt(promod3::core::EMatX3& pos_mat_one, 
@@ -76,18 +76,15 @@ void DoIt(promod3::core::EMatX3& pos_mat_one,
 
     for(uint j = 1; j < clusters[i].size(); ++j){
 
-      uint size = std::max(current_indices.size(),
-                           initial_indices[clusters[i][j]].size());
+      std::vector<uint> union_result;
 
-      std::vector<uint> intersection_result(size, 0);
+      std::set_union(current_indices.begin(),
+                     current_indices.end(),
+                     initial_indices[clusters[i][j]].begin(),
+                     initial_indices[clusters[i][j]].end(),
+                     std::back_inserter(union_result));
 
-      it = std::set_intersection(current_indices.begin(),
-                                 current_indices.end(),
-                                 initial_indices[clusters[i][j]].begin(),
-                                 initial_indices[clusters[i][j]].end(),
-                                 intersection_result.begin());
-
-      current_indices.assign(intersection_result.begin(),it);
+      current_indices = union_result;
     }
     indices.push_back(current_indices);
   }