diff --git a/cpp/particles/interpolation/particles_field_computer.hpp b/cpp/particles/interpolation/particles_field_computer.hpp
index 1a7b3d1816189aeefa018d378e5362bf9f3d1240..172ca37a0b3a303d9a25a966b682920f40c5822b 100644
--- a/cpp/particles/interpolation/particles_field_computer.hpp
+++ b/cpp/particles/interpolation/particles_field_computer.hpp
@@ -89,15 +89,6 @@ public:
     /// Computation related
     ////////////////////////////////////////////////////////////////////////
 
-    template <int size_particle_rhs>
-    void init_result_array(real_number particles_current_rhs[],
-                                   const partsize_t nb_particles) const {
-        // Set values to zero initialy
-        std::fill(particles_current_rhs,
-                  particles_current_rhs+nb_particles*size_particle_rhs,
-                  0);
-    }
-
     const std::array<real_number, 3> getSpatialBoxWidth() const
     {
         return this->spatial_box_width;
diff --git a/cpp/particles/p2p/p2p_computer.hpp b/cpp/particles/p2p/p2p_computer.hpp
index 49115bd2f5cc626fbae7758e45cf2ec6ef9d8344..5a399e515da2bf59bdfcb3874c7591daf0af1b4c 100644
--- a/cpp/particles/p2p/p2p_computer.hpp
+++ b/cpp/particles/p2p/p2p_computer.hpp
@@ -53,11 +53,6 @@ class p2p_computer{
 public:
     p2p_computer() : isActive(true){}
 
-    template <int size_particle_rhs>
-    void init_result_array(real_number rhs[], const partsize_t nbParticles) const{
-        memset(rhs, 0, sizeof(real_number)*nbParticles*size_particle_rhs);
-    }
-
     template <int size_particle_rhs>
     void reduce_particles_rhs(real_number rhs_dst[], const real_number rhs_src[], const partsize_t nbParticles) const{
         static_assert(size_particle_rhs == 6, "This kernel works only with 6 values per particle's rhs");
diff --git a/cpp/particles/p2p/p2p_computer_empty.hpp b/cpp/particles/p2p/p2p_computer_empty.hpp
index 5c43ef3c8a13d3095e3b172787bca8a8bdbe8c3b..b13a95e1f2bb8fd4f51f51c4cb3b6f74e00b079b 100644
--- a/cpp/particles/p2p/p2p_computer_empty.hpp
+++ b/cpp/particles/p2p/p2p_computer_empty.hpp
@@ -31,9 +31,6 @@
 template <class real_number, class partsize_t>
 class p2p_computer_empty{
 public:
-    template <int size_particle_rhs>
-    void init_result_array(real_number /*rhs*/[], const partsize_t /*nbParticles*/) const{
-    }
 
     template <int size_particle_rhs>
     void reduce_particles_rhs(real_number /*rhs_dst*/[], const real_number /*rhs_src*/[], const partsize_t /*nbParticles*/) const{
diff --git a/cpp/particles/p2p/p2p_distr_mpi.hpp b/cpp/particles/p2p/p2p_distr_mpi.hpp
index 694cae9af622b17f4d0151fb84ab936c82aee980..6639fb70e5d3df4f1540b8025c00020c7849de07 100644
--- a/cpp/particles/p2p/p2p_distr_mpi.hpp
+++ b/cpp/particles/p2p/p2p_distr_mpi.hpp
@@ -102,9 +102,11 @@ protected:
     std::array<long int,3> nb_cell_levels;
 
     template <class DataType, int sizeElement>
-    static void permute_copy(const partsize_t offsetIdx, const partsize_t nbElements,
-                             const std::pair<long int,partsize_t> permutation[],
-                             DataType data[], std::vector<unsigned char>* buffer){
+    static void permute_copy(const partsize_t offsetIdx,
+                             const partsize_t nbElements,
+                             const std::pair<long int, partsize_t> permutation[],
+                             DataType data[],
+                             std::vector<unsigned char>* buffer){
         buffer->resize(nbElements*sizeof(DataType)*sizeElement);
         DataType* dataBuffer = reinterpret_cast<DataType*>(buffer->data());
 
@@ -127,6 +129,8 @@ protected:
                         = dataBuffer[srcData*sizeElement + idxVal];
             }
         }
+
+        buffer->resize(0);
     }
 
     static int foundGridFactor(const real_number in_cutoff_radius, const std::array<real_number,3>& in_spatial_box_width){
@@ -340,9 +344,10 @@ public:
                                 current_offset_particles_for_partition[idxPartition],
                                 current_my_nb_particles_per_partition[idxPartition],
                                 part_to_sort.data(),
-                                       particles_positions,
-                                       &buffer);
-                for(int idx_rhs = 0 ; idx_rhs < nb_rhs ; ++idx_rhs){
+                                particles_positions,
+                                &buffer);
+                for(int idx_rhs = 0 ; idx_rhs < nb_rhs ; ++idx_rhs)
+                {
                     permute_copy<real_number, size_particle_rhs>(
                                     current_offset_particles_for_partition[idxPartition],
                                     current_my_nb_particles_per_partition[idxPartition],
@@ -354,14 +359,14 @@ public:
                                 current_offset_particles_for_partition[idxPartition],
                                 current_my_nb_particles_per_partition[idxPartition],
                                 part_to_sort.data(),
-                                       inout_index_particles,
-                                       &buffer);
+                                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);
+                                particles_coord.get(),
+                                &buffer);
             }
         }
 
@@ -495,27 +500,40 @@ public:
                     whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
                     mpiRequests.emplace_back();
                     assert(descriptor.nbParticlesToExchange*size_particle_positions < std::numeric_limits<int>::max());
-                    AssertMpi(MPI_Isend(const_cast<real_number*>(&particles_positions[particles_offset_layers[my_nb_cell_levels-descriptor.nbLevelsToExchange]*size_particle_positions]),
-                              int(descriptor.nbParticlesToExchange*size_particle_positions), particles_utils::GetMpiType(real_number()),
-                              descriptor.destProc, TAG_POSITION_PARTICLES,
-                              current_com, &mpiRequests.back()));
+                    AssertMpi(MPI_Isend(
+                                const_cast<real_number*>(&particles_positions[particles_offset_layers[my_nb_cell_levels-descriptor.nbLevelsToExchange]*size_particle_positions]),
+                                int(descriptor.nbParticlesToExchange*size_particle_positions),
+                                particles_utils::GetMpiType(real_number()),
+                                descriptor.destProc,
+                                TAG_POSITION_PARTICLES,
+                                current_com,
+                                &mpiRequests.back()));
 
                     whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
                     mpiRequests.emplace_back();
                     assert(descriptor.nbParticlesToExchange*size_particle_positions < std::numeric_limits<int>::max());
-                    AssertMpi(MPI_Isend(const_cast<partsize_t*>(&inout_index_particles[particles_offset_layers[my_nb_cell_levels-descriptor.nbLevelsToExchange]]),
-                              int(descriptor.nbParticlesToExchange), particles_utils::GetMpiType(partsize_t()),
-                              descriptor.destProc, TAG_INDEX_PARTICLES,
-                              current_com, &mpiRequests.back()));
+                    AssertMpi(MPI_Isend(
+                                const_cast<partsize_t*>(&inout_index_particles[particles_offset_layers[my_nb_cell_levels-descriptor.nbLevelsToExchange]]),
+                                int(descriptor.nbParticlesToExchange),
+                                particles_utils::GetMpiType(partsize_t()),
+                                descriptor.destProc,
+                                TAG_INDEX_PARTICLES,
+                                current_com,
+                                &mpiRequests.back()));
 
                     assert(descriptor.toRecvAndMerge == nullptr);
                     descriptor.toRecvAndMerge.reset(new real_number[descriptor.nbParticlesToExchange*size_particle_rhs]);
                     whatNext.emplace_back(std::pair<Action,int>{MERGE_PARTICLES, idxDescr});
                     mpiRequests.emplace_back();
                     assert(descriptor.nbParticlesToExchange*size_particle_rhs < std::numeric_limits<int>::max());
-                    AssertMpi(MPI_Irecv(descriptor.toRecvAndMerge.get(), int(descriptor.nbParticlesToExchange*size_particle_rhs),
-                                        particles_utils::GetMpiType(real_number()), descriptor.destProc, TAG_RESULT_PARTICLES,
-                                        current_com, &mpiRequests.back()));
+                    AssertMpi(MPI_Irecv(
+                                descriptor.toRecvAndMerge.get(),
+                                int(descriptor.nbParticlesToExchange*size_particle_rhs),
+                                particles_utils::GetMpiType(real_number()),
+                                descriptor.destProc,
+                                TAG_RESULT_PARTICLES,
+                                current_com,
+                                &mpiRequests.back()));
                 }
             }
             else{
@@ -524,9 +542,14 @@ public:
 #endif
                 whatNext.emplace_back(std::pair<Action,int>{RECV_PARTICLES, idxDescr});
                 mpiRequests.emplace_back();
-                AssertMpi(MPI_Irecv(&descriptor.nbParticlesToExchange,
-                      1, particles_utils::GetMpiType(partsize_t()), descriptor.destProc, TAG_NB_PARTICLES,
-                      current_com, &mpiRequests.back()));
+                AssertMpi(MPI_Irecv(
+                            &descriptor.nbParticlesToExchange,
+                            1,
+                            particles_utils::GetMpiType(partsize_t()),
+                            descriptor.destProc,
+                            TAG_NB_PARTICLES,
+                            current_com,
+                            &mpiRequests.back()));
             }
         }
 
@@ -536,8 +559,14 @@ public:
                 std::vector<int> willsendall(nb_processes_involved*nb_processes_involved, 0);// TODO debug
                 std::vector<int> willrecvall(nb_processes_involved*nb_processes_involved, 0);// TODO debug
 
-                MPI_Gather(willrecv.data(), nb_processes_involved, MPI_INT, willrecvall.data(),
-                            nb_processes_involved, MPI_INT, 0, MPI_COMM_WORLD);
+                MPI_Gather(willrecv.data(),
+                           nb_processes_involved,
+                           MPI_INT,
+                           willrecvall.data(),
+                           nb_processes_involved,
+                           MPI_INT,
+                           0,
+                           MPI_COMM_WORLD);
                 MPI_Gather(willsend.data(), nb_processes_involved, MPI_INT, willsendall.data(),
                             nb_processes_involved, MPI_INT, 0, MPI_COMM_WORLD);
 
@@ -603,17 +632,27 @@ public:
                             whatNext.emplace_back(std::pair<Action,int>{COMPUTE_PARTICLES, releasedAction.second});
                             mpiRequests.emplace_back();
                             assert(NbParticlesToReceive*size_particle_positions < std::numeric_limits<int>::max());
-                            AssertMpi(MPI_Irecv(descriptor.toCompute.get(), int(NbParticlesToReceive*size_particle_positions),
-                                                particles_utils::GetMpiType(real_number()), destProc, TAG_POSITION_PARTICLES,
-                                                current_com, &mpiRequests.back()));
+                            AssertMpi(MPI_Irecv(
+                                        descriptor.toCompute.get(),
+                                        int(NbParticlesToReceive*size_particle_positions),
+                                        particles_utils::GetMpiType(real_number()),
+                                        destProc,
+                                        TAG_POSITION_PARTICLES,
+                                        current_com,
+                                        &mpiRequests.back()));
 
                             descriptor.indexes.reset(new partsize_t[NbParticlesToReceive]);
                             whatNext.emplace_back(std::pair<Action,int>{COMPUTE_PARTICLES, releasedAction.second});
                             mpiRequests.emplace_back();
                             assert(NbParticlesToReceive*size_particle_positions < std::numeric_limits<int>::max());
-                            AssertMpi(MPI_Irecv(descriptor.indexes.get(), int(NbParticlesToReceive),
-                                                particles_utils::GetMpiType(partsize_t()), destProc, TAG_INDEX_PARTICLES,
-                                                current_com, &mpiRequests.back()));
+                            AssertMpi(MPI_Irecv(
+                                        descriptor.indexes.get(),
+                                        int(NbParticlesToReceive),
+                                        particles_utils::GetMpiType(partsize_t()),
+                                        destProc,
+                                        TAG_INDEX_PARTICLES,
+                                        current_com,
+                                        &mpiRequests.back()));
                         }
                     }
 
@@ -633,8 +672,12 @@ public:
 
                             assert(descriptor.toCompute != nullptr);
                             assert(descriptor.indexes != nullptr);
+                            // allocate rhs buffer
                             descriptor.results.reset(new real_number[NbParticlesToReceive*size_particle_rhs]);
-                            computer_thread.template init_result_array<size_particle_rhs>(descriptor.results.get(), NbParticlesToReceive);
+                            // clean up rhs buffer
+                            set_particle_data_to_zero<partsize_t, real_number, size_particle_rhs>(
+                                    descriptor.results.get(),
+                                    NbParticlesToReceive);
 
                             // Compute
                             partsize_t idxPart = 0;
@@ -716,9 +759,14 @@ public:
                             whatNext.emplace_back(std::pair<Action,int>{RELEASE_BUFFER_PARTICLES, releasedAction.second});
                             mpiRequests.emplace_back();
                             assert(NbParticlesToReceive*size_particle_rhs < std::numeric_limits<int>::max());
-                            AssertMpi(MPI_Isend(descriptor.results.get(), int(NbParticlesToReceive*size_particle_rhs),
-                                                particles_utils::GetMpiType(real_number()), destProc, TAG_RESULT_PARTICLES,
-                                                current_com, &mpiRequests.back()));
+                            AssertMpi(MPI_Isend(
+                                        descriptor.results.get(),
+                                        int(NbParticlesToReceive*size_particle_rhs),
+                                        particles_utils::GetMpiType(real_number()),
+                                        destProc,
+                                        TAG_RESULT_PARTICLES,
+                                        current_com,
+                                        &mpiRequests.back()));
                             delete[] descriptor.toCompute.release();
                             delete[] descriptor.indexes.release();
                         }
@@ -740,8 +788,10 @@ public:
                         NeighborDescriptor& descriptor = neigDescriptors[releasedAction.second];
                         assert(descriptor.isRecv == false);
                         assert(descriptor.toRecvAndMerge != nullptr);
-                        computer_thread.template reduce_particles_rhs<size_particle_rhs>(&particles_current_rhs[0][particles_offset_layers[my_nb_cell_levels-descriptor.nbLevelsToExchange]*size_particle_rhs],
-                                descriptor.toRecvAndMerge.get(), descriptor.nbParticlesToExchange);
+                        computer_thread.template reduce_particles_rhs<size_particle_rhs>(
+                                &particles_current_rhs[0][particles_offset_layers[my_nb_cell_levels-descriptor.nbLevelsToExchange]*size_particle_rhs],
+                                descriptor.toRecvAndMerge.get(),
+                                descriptor.nbParticlesToExchange);
                         delete[] descriptor.toRecvAndMerge.release();
                     }
                 }
diff --git a/cpp/particles/p2p/p2p_ghost_collisions.hpp b/cpp/particles/p2p/p2p_ghost_collisions.hpp
index 166ac5fc1f45d363b6e908d4dd1121b985fd69e0..5dee282ac94c3f8b4a9427477e8a5a14936139b0 100644
--- a/cpp/particles/p2p/p2p_ghost_collisions.hpp
+++ b/cpp/particles/p2p/p2p_ghost_collisions.hpp
@@ -41,12 +41,16 @@ void print_pair_vec(std::vector<partsize_t> vec)
 template <class real_number, class partsize_t>
 class p2p_ghost_collisions
 {
-    private:
+    protected:
         const double pi  = atan(1.0)*4;
         const double pi2 = atan(1.0)*8;
-        // depending on the particle shape, we will be deciding whether particles intersect in different ways
+
         enum particle_shape {CYLINDER, SPHERE, DISK};
-        particle_shape current_particle_shape;
+
+    private:
+        bool isActive;
+
+        bool synchronisation;
 
         // description for cylinders:
         double cylinder_width;
@@ -54,8 +58,11 @@ class p2p_ghost_collisions
 
         // description for disks:
         double disk_width;
-        bool isActive;
-   protected:        
+
+        // depending on the particle shape, we will be deciding whether particles intersect in different ways
+        particle_shape current_particle_shape;
+
+    protected:
         void add_colliding_pair(partsize_t idx_part1, partsize_t idx_part2)
         {
             // store colliding particle ids in order, to be able to identify pairs more easily
@@ -70,15 +77,21 @@ class p2p_ghost_collisions
             this->collision_pairs_local.push_back(idx_part_small);
             this->collision_pairs_local.push_back(idx_part_big);
         }
-        bool synchronisation;
+
         std::vector <partsize_t> collision_pairs_local;
         std::vector <partsize_t> collision_pairs_global;
 
 public:
-    p2p_ghost_collisions(): current_particle_shape(SPHERE), cylinder_width(1.0), cylinder_length(1.0), disk_width(1.0), isActive(true), synchronisation(false) {}
+    p2p_ghost_collisions():
+        isActive(true),
+        synchronisation(false),
+        cylinder_width(1.0),
+        cylinder_length(1.0),
+        disk_width(1.0),
+        current_particle_shape(SPHERE) {}
 
-    // Copy constructor use a counter set to zero
-    p2p_ghost_collisions(const p2p_ghost_collisions& src)
+
+    void copy_from_p2p_ghost_collisions(const p2p_ghost_collisions<real_number, partsize_t>& src)
     {
         this->current_particle_shape = src.current_particle_shape;
         this->cylinder_width = src.cylinder_width;
@@ -89,8 +102,20 @@ public:
         this->collision_pairs_local.reserve(src.collision_pairs_local.capacity());
     }
 
-    template <int size_particle_rhs>
-    void init_result_array(real_number /*rhs*/[], const partsize_t /*nbParticles*/) const{
+    // Copy constructor use a counter set to zero
+    p2p_ghost_collisions(const p2p_ghost_collisions<real_number, partsize_t>& src)
+    {
+        this->copy_from_p2p_ghost_collisions(src);
+    }
+
+    double rod_distance(double px, double py, double pz, double qx, double qy, double qz, double t, double s, double x, double y, double z){
+        double x_dist, y_dist, z_dist;
+        double min_distance;
+        x_dist = this->cylinder_length*0.5*t*px-this->cylinder_length*0.5*s*qx+x;
+        y_dist = this->cylinder_length*0.5*t*py-this->cylinder_length*0.5*s*qy+y;
+        z_dist = this->cylinder_length*0.5*t*pz-this->cylinder_length*0.5*s*qz+z;
+        min_distance = sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist);
+        return min_distance;
     }
 
     template <int size_particle_rhs>
@@ -208,7 +233,8 @@ public:
                 break;
             case CYLINDER:
                 {
-                    double pq, xp, xq, t, s, x_dist, y_dist, z_dist ,min_distance;
+                    double pq, xp, xq, t, s, min_distance, min_distance_current;
+                    double px,py, pz, qx, qy, qz;
 
                     const double x = xseparation;
                     const double y = yseparation;
@@ -217,12 +243,17 @@ public:
 
                     if( r <= this->cylinder_length )
                     {
-
                         /* p and q are the orientation vectors of the first and second particles. */
+                        px = pos_part1[IDXC_X+3];
+                        py = pos_part1[IDXC_Y+3];
+                        pz = pos_part1[IDXC_Z+3];
+                        qx = pos_part2[IDXC_X+3];
+                        qy = pos_part2[IDXC_Y+3];
+                        qz = pos_part2[IDXC_Z+3];
                         /* pq, xp, xq are scalar products of these vectors with x, relative position */
-                        pq = pos_part1[IDXC_X+3] * pos_part2[IDXC_X+3] + pos_part1[IDXC_Y+3] * pos_part2[IDXC_Y+3] + pos_part1[IDXC_Z+3] * pos_part2[IDXC_Z+3];
-                        xp = x * pos_part1[IDXC_X+3] + y * pos_part1[IDXC_Y+3] + z * pos_part1[IDXC_Z+3];
-                        xq = x * pos_part2[IDXC_X+3] + y * pos_part2[IDXC_Y+3] + z * pos_part2[IDXC_Z+3];
+                        pq = px * qx + py * qy + pz * qz;
+                        xp = x * px + y * py + z * pz;
+                        xq = x * qx + y * qy + z * qz;
                         /* t and s parametrize the two rods. Find min distance: */
                         assert(this->cylinder_length > 0);
                         t = 2.0/(this->cylinder_length*(pq*pq-1.0))*(xp-pq*xq);
@@ -231,10 +262,7 @@ public:
                         if( abs(t)<=1.0 and abs(s)<=1.0 )
                         {
                             /* Get minimal distance in case of both t and s in {-1,1}. Else: check edges */
-                            x_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_X+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_X+3]+x;
-                            y_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Y+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Y+3]+y;
-                            z_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Z+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Z+3]+z;
-                            min_distance = sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist);
+                            min_distance = this->rod_distance(px, py, pz, qx, qy, qz, t, s, x, y, z);
                         }
                         else
                         {
@@ -243,34 +271,26 @@ public:
                             t = 1.0;
                             s = t*pq+2.0/this->cylinder_length*xq;
                             if( abs(s)>1.0 ) { s = s / abs(s) ;}
-                            x_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_X+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_X+3]+x;
-                            y_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Y+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Y+3]+y;
-                            z_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Z+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Z+3]+z;
-                            min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
+                                min_distance_current = this->rod_distance(px, py, pz, qx, qy, qz, t, s, x, y, z);
+                                min_distance = fmin( min_distance_current, min_distance );
                             /* t fixed at -1, find min along s */
                             t = -1.0;
                             s = t*pq+2.0/this->cylinder_length*xq;
                             if( abs(s)>1.0 ) { s = s / abs(s) ;}
-                            x_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_X+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_X+3]+x;
-                            y_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Y+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Y+3]+y;
-                            z_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Z+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Z+3]+z;
-                            min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
+                                min_distance_current = this->rod_distance(px, py, pz, qx, qy, qz, t, s, x, y, z);
+                                min_distance = fmin( min_distance_current, min_distance );
                             /* s fixed at 1, find min along t */
                             s = 1.0;
                             t = s*pq-2.0/this->cylinder_length*xp;
                             if( abs(t)>1.0 ) { t = t / abs(t) ;}
-                            x_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_X+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_X+3]+x;
-                            y_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Y+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Y+3]+y;
-                            z_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Z+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Z+3]+z;
-                            min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
+                                min_distance_current = this->rod_distance(px, py, pz, qx, qy, qz, t, s, x, y, z);
+                                min_distance = fmin( min_distance_current, min_distance );
                             /* s fixed at -1, find min along t */
                             s = -1.0;
                             t = s*pq-2.0/this->cylinder_length*xp;
                             if( abs(t)>1.0 ) { t = t / abs(t) ;}
-                            x_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_X+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_X+3]+x;
-                            y_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Y+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Y+3]+y;
-                            z_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Z+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Z+3]+z;
-                            min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
+                                min_distance_current = this->rod_distance(px, py, pz, qx, qy, qz, t, s, x, y, z);
+                                min_distance = fmin( min_distance_current, min_distance );
                         }
                         /* If cylinders overlap count it as a collision */
                         if( min_distance<=this->cylinder_width )
@@ -376,12 +396,12 @@ public:
         this->cylinder_length = LENGTH;
     }
 
-    double get_cylinder_width()
+    double get_cylinder_width() const
     {
         return this->cylinder_width;
     }
 
-    double get_cylinder_length()
+    double get_cylinder_length() const
     {
         return this->cylinder_length;
     }
@@ -391,12 +411,17 @@ public:
         this->current_particle_shape = DISK;
     }
 
+    particle_shape get_current_particle_shape() const
+    {
+        return this->current_particle_shape;
+    }
+
     void set_disk_width(const double WIDTH)
     {
         this->disk_width = WIDTH;
     }
 
-    double get_disk_width()
+    double get_disk_width() const
     {
         return this->disk_width;
     }
@@ -405,9 +430,15 @@ public:
         return isActive;
     }
 
-    void setEnable(const bool inIsActive) {
+    void setEnable(const bool inIsActive)
+    {
         isActive = inIsActive;
     }
+
+    bool isSynchronized() const
+    {
+        return this->synchronisation;
+    }
 };
 
 
diff --git a/cpp/particles/p2p/p2p_merge_collisions.hpp b/cpp/particles/p2p/p2p_merge_collisions.hpp
index cf93f3e8001699b3cb8de2db4c043409bae4861b..4b728c8f3c3db1c26e41947fe1a23325fbe688e8 100644
--- a/cpp/particles/p2p/p2p_merge_collisions.hpp
+++ b/cpp/particles/p2p/p2p_merge_collisions.hpp
@@ -37,10 +37,6 @@ public:
     // Copy constructor use a counter set to zero
     p2p_merge_collisions(const p2p_merge_collisions&){}
 
-    template <int size_particle_rhs>
-    void init_result_array(real_number /*rhs*/[], const partsize_t /*nbParticles*/) const{
-    }
-
     template <int size_particle_rhs>
     void reduce_particles_rhs(real_number /*rhs_dst*/[], const real_number /*rhs_src*/[], const partsize_t /*nbParticles*/) const{
     }
diff --git a/cpp/particles/particles_distr_mpi.hpp b/cpp/particles/particles_distr_mpi.hpp
index aa6478f0e37640f57b1836f25485ec2c14e2c189..df2dd6b9444d6371ced0f9cfb51925234657b68f 100644
--- a/cpp/particles/particles_distr_mpi.hpp
+++ b/cpp/particles/particles_distr_mpi.hpp
@@ -413,11 +413,19 @@ public:
                         const partsize_t NbParticlesToReceive = descriptor.nbParticlesToRecv;
 
                         assert(descriptor.toCompute != nullptr);
+                        // allocate buffer
                         descriptor.results.reset(new real_number[NbParticlesToReceive*size_particle_rhs]);
-                        in_computer.template init_result_array<size_particle_rhs>(descriptor.results.get(), NbParticlesToReceive);
+                        // clean up buffer
+                        set_particle_data_to_zero<partsize_t, real_number, size_particle_rhs>(
+                                descriptor.results.get(),
+                                NbParticlesToReceive);
 
                         if(more_than_one_thread == false){
-                            in_computer.template apply_computation<field_class, size_particle_positions, size_particle_rhs>(in_field, descriptor.toCompute.get(), descriptor.results.get(), NbParticlesToReceive);
+                            in_computer.template apply_computation<field_class, size_particle_positions, size_particle_rhs>(
+                                    in_field,
+                                    descriptor.toCompute.get(),
+                                    descriptor.results.get(),
+                                    NbParticlesToReceive);
                         }
                         else{
                             TIMEZONE_OMP_INIT_PRETASK(timeZoneTaskKey)
diff --git a/cpp/particles/particles_system.hpp b/cpp/particles/particles_system.hpp
index 9487f870db874731103e9f8cd5609fa334c9d89b..b47763028aa25d040c51ffcde11a3dd3cfe7a87d 100644
--- a/cpp/particles/particles_system.hpp
+++ b/cpp/particles/particles_system.hpp
@@ -298,8 +298,11 @@ public:
         if(computer_p2p.isEnable() == true){
             TIMEZONE("particles_system::compute_p2p");
             distr_p2p.template compute_distr<p2p_computer_class, size_particle_positions, size_particle_rhs>(
-                            computer_p2p, current_my_nb_particles_per_partition.get(),
-                            my_particles_positions.get(), my_particles_rhs.data(), int(my_particles_rhs.size()),
+                            computer_p2p,
+                            current_my_nb_particles_per_partition.get(),
+                            my_particles_positions.get(),
+                            my_particles_rhs.data(),
+                            int(my_particles_rhs.size()),
                             my_particles_positions_indexes.get());
         }
     }
diff --git a/cpp/particles/particles_utils.hpp b/cpp/particles/particles_utils.hpp
index fcb9802231603f7b7f182e679b286857212c4919..4cef3a3b4d4c44d44408a328497ebca94a7df8a9 100644
--- a/cpp/particles/particles_utils.hpp
+++ b/cpp/particles/particles_utils.hpp
@@ -398,4 +398,19 @@ std::vector<real_number> BuildLimitsAllProcesses(
     return spatial_limit_per_proc;
 }
 
-#endif
+template <typename partsize_t, typename rnumber, int size_of_particle>
+int set_particle_data_to_zero(
+        rnumber __restrict__ *data,
+        const partsize_t numberParticles)
+{
+    // TODO: ensure simd.
+    // don't use openmp here, as this function may be called from within openmp parallel regions
+    std::fill_n(
+            data,
+            numberParticles*size_of_particle,
+            0);
+    return EXIT_SUCCESS;
+}
+
+#endif//PARTICLES_UTILS_HPP
+