diff --git a/cpp/particles/p2p_distr_mpi.hpp b/cpp/particles/p2p_distr_mpi.hpp
index 274dc8a64137427aeb3ce76aa0efe801a30e586e..8d2bb40919022c1a113bd519a6d3f8507d156df0 100644
--- a/cpp/particles/p2p_distr_mpi.hpp
+++ b/cpp/particles/p2p_distr_mpi.hpp
@@ -333,19 +333,32 @@ public:
                 });
 
                 // Permute array using buffer
+		// permute 4th function parameter using buffer, based on information in first 3 parameters
                 std::vector<unsigned char> buffer;
-                permute_copy<real_number, size_particle_positions>(current_offset_particles_for_partition[idxPartition],
-                                                                   current_my_nb_particles_per_partition[idxPartition],
-                                                                   part_to_sort.data(), particles_positions, &buffer);
-                permute_copy<real_number, size_particle_rhs>(current_offset_particles_for_partition[idxPartition],
-                                                             current_my_nb_particles_per_partition[idxPartition],
-                                                             part_to_sort.data(), particles_current_rhs, &buffer);
-                permute_copy<partsize_t, 1>(current_offset_particles_for_partition[idxPartition],
-                                            current_my_nb_particles_per_partition[idxPartition],
-                                            part_to_sort.data(), inout_index_particles, &buffer);
-                permute_copy<long int, 1>(current_offset_particles_for_partition[idxPartition],
-                                            current_my_nb_particles_per_partition[idxPartition],
-                                            part_to_sort.data(), particles_coord.get(), &buffer);
+                permute_copy<real_number, size_particle_positions>(
+				current_offset_particles_for_partition[idxPartition],
+                                current_my_nb_particles_per_partition[idxPartition],
+                                part_to_sort.data(),
+			       	particles_positions,
+			       	&buffer);
+                permute_copy<real_number, size_particle_rhs>(
+				current_offset_particles_for_partition[idxPartition],
+                                current_my_nb_particles_per_partition[idxPartition],
+                                part_to_sort.data(),
+			       	particles_current_rhs,
+			       	&buffer);
+                permute_copy<partsize_t, 1>(
+				current_offset_particles_for_partition[idxPartition],
+                                current_my_nb_particles_per_partition[idxPartition],
+                                part_to_sort.data(),
+			       	inout_index_particles,
+			       	&buffer);
+                permute_copy<long int, 1>(
+				current_offset_particles_for_partition[idxPartition],
+                                current_my_nb_particles_per_partition[idxPartition],
+                                part_to_sort.data(),
+			       	particles_coord.get(),
+			       	&buffer);
             }
         }
 
@@ -639,7 +652,12 @@ public:
                                     const std::vector<std::pair<partsize_t,partsize_t>>* neighbors[27];
                                     long int neighbors_indexes[27];
                                     std::array<real_number,3> shift[27];
-                                    const int nbNeighbors = my_tree.getNeighbors(current_cell_idx, neighbors, neighbors_indexes, shift, true);
+                                    const int nbNeighbors = my_tree.getNeighbors(
+						    current_cell_idx,
+						    neighbors,
+						    neighbors_indexes,
+						    shift,
+						    true);
 
                                     // with other interval
                                     for(int idx_neighbor = 0 ; idx_neighbor < nbNeighbors ; ++idx_neighbor){
@@ -648,13 +666,14 @@ public:
                                         for(size_t idx_2 = 0 ; idx_2 < (*neighbors[idx_neighbor]).size() ; ++idx_2){
                                             for(partsize_t idx_p1 = 0 ; idx_p1 < nb_parts_in_cell ; ++idx_p1){
                                                 for(partsize_t idx_p2 = 0 ; idx_p2 < (*neighbors[idx_neighbor])[idx_2].second ; ++idx_p2){
-                                                    const real_number dist_r2 = compute_distance_r2(descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions + IDXC_X],
-                                                                                                    descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions + IDXC_Y],
-                                                                                                    descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions + IDXC_Z],
-                                                                                                    particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDXC_X],
-                                                                                                    particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDXC_Y],
-                                                                                                    particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDXC_Z],
-                                                                                                    shift[idx_neighbor][IDXC_X], shift[idx_neighbor][IDXC_Y], shift[idx_neighbor][IDXC_Z]);
+                                                    const real_number dist_r2 = compute_distance_r2(
+								    descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions + IDXC_X],
+                                                                    descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions + IDXC_Y],
+                                                                    descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions + IDXC_Z],
+                                                                    particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDXC_X],
+                                                                    particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDXC_Y],
+                                                                    particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDXC_Z],
+                                                                    shift[idx_neighbor][IDXC_X], shift[idx_neighbor][IDXC_Y], shift[idx_neighbor][IDXC_Z]);
                                                     if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
                                                         computer_thread.template compute_interaction<size_particle_positions, size_particle_rhs>(
                                                                             descriptor.indexes[(idxPart+idx_p1)],
@@ -663,7 +682,11 @@ public:
                                                                             inout_index_particles[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)],
                                                                             &particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions],
                                                                             &particles_current_rhs[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_rhs],
-                                                                            dist_r2, cutoff_radius_compute, shift[idx_neighbor][IDXC_X], shift[idx_neighbor][IDXC_Y], shift[idx_neighbor][IDXC_Z]);
+                                                                            dist_r2,
+									    cutoff_radius_compute,
+									    shift[idx_neighbor][IDXC_X],
+									    shift[idx_neighbor][IDXC_Y],
+									    shift[idx_neighbor][IDXC_Z]);
                                                     }
                                                 }
                                             }
diff --git a/cpp/particles/p2p_ghost_collisions.hpp b/cpp/particles/p2p_ghost_collisions.hpp
index cb7270fb76128be90ed1a27255a150cfdec48dd6..e3169bfade1c1bd491a601b7d42dcbe72eb747a8 100644
--- a/cpp/particles/p2p_ghost_collisions.hpp
+++ b/cpp/particles/p2p_ghost_collisions.hpp
@@ -32,8 +32,8 @@ std::vector<partsize_t> pairs2vec(std::set <std::pair <partsize_t,partsize_t>> I
   std::vector<partsize_t> v(2*ID_pairs.size());
   for(unsigned int i=0; i < ID_pairs.size(); i++)
   {
-    v[2*i-2] = std::get<0>(*std::next(ID_pairs.begin(), i));
-    v[2*i-1] = std::get<1>(*std::next(ID_pairs.begin(), i));
+    v[2*i] = (*std::next(ID_pairs.begin(), i)).first;
+    v[2*i+1] = (*std::next(ID_pairs.begin(), i)).second;
   }
   return v;
 }
@@ -43,7 +43,9 @@ std::set <std::pair <partsize_t,partsize_t>> vec2pairs(std::vector<partsize_t> v
   std::set <std::pair <partsize_t,partsize_t>> ID_pairs;
   for(unsigned int i=0; i < v.size()/2; i++)
     {
-      ID_pairs.insert(std::pair <partsize_t,partsize_t> (v[2*i],v[2*i+1]));
+      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;
 }
@@ -69,14 +71,20 @@ public:
 
     template <int size_particle_positions, int size_particle_rhs>
     void compute_interaction(const partsize_t idx_part1,
-                             const real_number /*pos_part1*/[], real_number /*rhs_part1*/[],
+                             const real_number /*pos_part1*/[],
+			     real_number /*rhs_part1*/[],
                              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*/){
+                             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*/){
         collision_counter += 1;
         std::pair <partsize_t, partsize_t> single_collision_pair(idx_part1, idx_part2);
         this->collision_pairs.insert(single_collision_pair);
+        //assert(idx_part1==0 and idx_part2==0);
     }
 
     void merge(const p2p_ghost_collisions& other){