diff --git a/cpp/particles/p2p_cylinder_collisions.hpp b/cpp/particles/p2p_cylinder_collisions.hpp
index 45d781f6f142dea19910a15ff575f50158784dbe..1f4f8c9bf6d0a9ade18108815ebb6fd493431a61 100644
--- a/cpp/particles/p2p_cylinder_collisions.hpp
+++ b/cpp/particles/p2p_cylinder_collisions.hpp
@@ -26,9 +26,40 @@
 #include <cstring>
 #include <math.h>       /* for sqrt, abs */
 
+#include <set>
+#include <utility>
+#include <vector>
+template < class partsize_t>
+std::vector<partsize_t> pairs2vec(std::set <std::pair <partsize_t,partsize_t>> ID_pairs){
+  std::vector<partsize_t> v(2*ID_pairs.size());
+  int i = 0;
+  for(auto p: ID_pairs)
+  {
+    v[2*i] = p.first;
+    v[2*i+1] = p.second;
+    i++;
+  }
+  return v;
+}
+
+template < class partsize_t>
+std::set <std::pair <partsize_t,partsize_t>> vec2pairs(std::vector<partsize_t> v){
+  std::set <std::pair <partsize_t,partsize_t>> ID_pairs;
+  assert(v.size()%2 == 0);
+  for(unsigned int i=0; i < v.size()/2; i++)
+    {
+      //DEBUG_MSG((std::to_string(v[2*i])+" "+std::to_string(v[2*i+1])+"\n").c_str());
+      std::pair <partsize_t, partsize_t> single_collision_pair(v[2*i],v[2*i+1]);
+      ID_pairs.insert(single_collision_pair);
+    }
+  return ID_pairs;
+}
+
+
 template <class real_number, class partsize_t>
 class p2p_cylinder_ghost_collisions{
     long int collision_counter;
+    std::set <std::pair <partsize_t, partsize_t>> collision_pairs;
     double width;
     double length;
     /* Following doubles are needed for the collision computation */
@@ -49,9 +80,9 @@ public:
     }
 
     template <int size_particle_positions, int size_particle_rhs>
-    void compute_interaction(const partsize_t /*idx_part1*/,
+    void compute_interaction(const partsize_t idx_part1,
                              const real_number pos_part1[], real_number /*rhs_part1*/[],
-                             const partsize_t /*idx_part2*/,
+                             const partsize_t idx_part2,
                              const real_number pos_part2[], real_number /*rhs_part2*/[],
                              const real_number /*dist_pow2*/,  const real_number /*cutoff*/,
                              const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/){
@@ -115,10 +146,19 @@ public:
         /* If cylinders overlap count it as a collision */
         if( min_distance<=width )
             collision_counter += 1;
+			std::pair <partsize_t, partsize_t> single_collision_pair(idx_part1, idx_part2);
+			this->collision_pairs.insert(single_collision_pair);
+			//DEBUG_MSG("inside compute interaction idx_part1 = %ld and idx_part2 = %ld\n", idx_part1, idx_part2);
+			assert(idx_part1!=idx_part2);
     }
 
     void merge(const p2p_cylinder_ghost_collisions& other){
         collision_counter += other.collision_counter;
+		std::set <std::pair <partsize_t, partsize_t>> new_collision_pairs;
+		std::set_union(collision_pairs.begin(), collision_pairs.end(),
+			other.collision_pairs.begin(), other.collision_pairs.end(),
+			std::inserter(new_collision_pairs, new_collision_pairs.begin()));
+		collision_pairs = new_collision_pairs;
     }
 
     constexpr static bool isEnable() {
@@ -152,6 +192,15 @@ public:
     {
         return this->length;
     }
+    
+	std::set <std::pair <partsize_t,partsize_t>> get_collision_pairs() const{
+		return collision_pairs;
+	}
+
+	void reset_collision_pairs(){
+		collision_pairs.clear();
+	}
+    
 };