diff --git a/cpp/particles/p2p/p2p_computer.hpp b/cpp/particles/p2p/p2p_computer.hpp
index 47d193e095403d329cd704e8623567c0f080a39b..ec8b2bcd35d765939a04d178ed65e3b60888a245 100644
--- a/cpp/particles/p2p/p2p_computer.hpp
+++ b/cpp/particles/p2p/p2p_computer.hpp
@@ -71,11 +71,16 @@ 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{
+                             const real_number pos_part2[],
+                             real_number rhs_part2[],
+                             const real_number dist_pow2,
+                             const real_number cutoff,
+                             const real_number /*xseparation*/,
+                             const real_number /*yseparation*/,
+                             const real_number /*zseparation*/) const{
         static_assert(size_particle_positions == 6, "This kernel works only with 6 values for one particle's position+orientation");
         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 ce4e282c3bdb7ff4bfc479b5a0b16d106f0492c1..03d8cd64995ed6a190f1d444d2846d230985dec8 100644
--- a/cpp/particles/p2p/p2p_computer_empty.hpp
+++ b/cpp/particles/p2p/p2p_computer_empty.hpp
@@ -44,8 +44,11 @@ public:
                              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{
+                             const real_number /*dist_pow2*/,
+                             const real_number /*cutoff*/,
+                             const real_number /*xseparation*/,
+                             const real_number /*yseparation*/,
+                             const real_number /*zseparation*/) const{
     }
 
     void merge(const p2p_computer_empty&){}
diff --git a/cpp/particles/p2p/p2p_distr_mpi.hpp b/cpp/particles/p2p/p2p_distr_mpi.hpp
index a525d154c7fe78d3bf3a4b122d13607ae514737c..34281d6180800fc0ea4bda7d747289947181c1db 100644
--- a/cpp/particles/p2p/p2p_distr_mpi.hpp
+++ b/cpp/particles/p2p/p2p_distr_mpi.hpp
@@ -260,15 +260,16 @@ public:
 
     real_number compute_distance_r2(const real_number x1, const real_number y1, const real_number z1,
                                     const real_number x2, const real_number y2, const real_number z2,
-                                    const real_number xshift_coef, const real_number yshift_coef, const real_number zshift_coef) const {
-        real_number diff_x = std::abs(apply_pbc(x1,IDXC_X)-apply_pbc(x2,IDXC_X)+xshift_coef*spatial_box_width[IDXC_X]);
-        assert(diff_x <= 2*cutoff_radius);
+                                    const real_number xshift_coef, const real_number yshift_coef, const real_number zshift_coef,
+                                    real_number &diff_x, real_number &diff_y, real_number &diff_z) const {
+        diff_x = apply_pbc(x1,IDXC_X)-apply_pbc(x2,IDXC_X)+xshift_coef*spatial_box_width[IDXC_X];
+        assert(std::abs(diff_x) <= 2*cutoff_radius);
 
-        real_number diff_y = std::abs(apply_pbc(y1,IDXC_X)-apply_pbc(y2,IDXC_X)+yshift_coef*spatial_box_width[IDXC_Y]);
-        assert(diff_y <= 2*cutoff_radius);
+        diff_y = apply_pbc(y1,IDXC_X)-apply_pbc(y2,IDXC_X)+yshift_coef*spatial_box_width[IDXC_Y];
+        assert(std::abs(diff_y) <= 2*cutoff_radius);
 
-        real_number diff_z = std::abs(apply_pbc(z1,IDXC_X)-apply_pbc(z2,IDXC_X)+zshift_coef*spatial_box_width[IDXC_Z]);
-        assert(diff_z <= 2*cutoff_radius);
+        diff_z = apply_pbc(z1,IDXC_X)-apply_pbc(z2,IDXC_X)+zshift_coef*spatial_box_width[IDXC_Z];
+        assert(std::abs(diff_z) <= 2*cutoff_radius);
 
         return (diff_x*diff_x) + (diff_y*diff_y) + (diff_z*diff_z);
     }
@@ -333,14 +334,14 @@ public:
                 });
 
                 // Permute array using buffer
-		// permute 4th function parameter using buffer, based on information in first 3 parameters
+                // 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_offset_particles_for_partition[idxPartition],
                                 current_my_nb_particles_per_partition[idxPartition],
                                 part_to_sort.data(),
-			       	particles_positions,
-			       	&buffer);
+                                       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],
@@ -350,17 +351,17 @@ public:
                                     &buffer);
                 }
                 permute_copy<partsize_t, 1>(
-				current_offset_particles_for_partition[idxPartition],
+                                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_offset_particles_for_partition[idxPartition],
                                 current_my_nb_particles_per_partition[idxPartition],
                                 part_to_sort.data(),
-			       	particles_coord.get(),
-			       	&buffer);
+                                       particles_coord.get(),
+                                       &buffer);
             }
         }
 
@@ -654,12 +655,13 @@ public:
                                     const std::vector<std::pair<partsize_t,partsize_t>>* neighbors[27];
                                     long int neighbors_indexes[27];
                                     std::array<real_number,3> shift[27];
+                                    real_number diff_x, diff_y, diff_z;
                                     const int nbNeighbors = my_tree.getNeighbors(
-						    current_cell_idx,
-						    neighbors,
-						    neighbors_indexes,
-						    shift,
-						    true);
+                                                    current_cell_idx,
+                                                    neighbors,
+                                                    neighbors_indexes,
+                                                    shift,
+                                                    true);
 
                                     // with other interval
                                     for(int idx_neighbor = 0 ; idx_neighbor < nbNeighbors ; ++idx_neighbor){
@@ -669,13 +671,18 @@ public:
                                             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_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]);
+                                                                    shift[idx_neighbor][IDXC_X],
+                                                                    shift[idx_neighbor][IDXC_Y],
+                                                                    shift[idx_neighbor][IDXC_Z],
+                                                                    diff_x,
+                                                                    diff_y,
+                                                                    diff_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)],
@@ -685,10 +692,10 @@ public:
                                                                             &particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions],
                                                                             &particles_current_rhs[0][((*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]);
+                                                                            cutoff_radius_compute,
+                                                                            diff_x,
+                                                                            diff_y,
+                                                                            diff_z);
                                                     }
                                                 }
                                             }
@@ -761,13 +768,15 @@ public:
                         // self interval
                         for(partsize_t idx_p1 = 0 ; idx_p1 < intervals[idx_1].second ; ++idx_p1){
                             for(partsize_t idx_p2 = idx_p1+1 ; idx_p2 < intervals[idx_1].second ; ++idx_p2){
+                                real_number diff_x, diff_y, diff_z;
                                 const real_number dist_r2 = compute_distance_r2(particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDXC_X],
                                                                                 particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDXC_Y],
                                                                                 particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDXC_Z],
                                                                                 particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions + IDXC_X],
                                                                                 particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions + IDXC_Y],
                                                                                 particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions + IDXC_Z],
-                                                                                0, 0, 0);
+                                                                                0, 0, 0,
+                                                                                diff_x, diff_y, diff_z);
                                 if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
                                     computer_thread.template compute_interaction<size_particle_positions,size_particle_rhs>(
                                                         inout_index_particles[(intervals[idx_1].first+idx_p1)],
@@ -776,7 +785,7 @@ public:
                                                         inout_index_particles[(intervals[idx_1].first+idx_p2)],
                                                         &particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions],
                                                         &particles_current_rhs[0][(intervals[idx_1].first+idx_p2)*size_particle_rhs],
-                                                        dist_r2, cutoff_radius_compute, 0, 0, 0);
+                                                        dist_r2, cutoff_radius_compute, diff_x, diff_y, diff_z);
                                 }
                             }
                         }
@@ -785,13 +794,15 @@ public:
                         for(size_t idx_2 = idx_1+1 ; idx_2 < intervals.size() ; ++idx_2){
                             for(partsize_t idx_p1 = 0 ; idx_p1 < intervals[idx_1].second ; ++idx_p1){
                                 for(partsize_t idx_p2 = 0 ; idx_p2 < intervals[idx_2].second ; ++idx_p2){
+                                    real_number diff_x, diff_y, diff_z;
                                     const real_number dist_r2 = compute_distance_r2(particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDXC_X],
                                                                                     particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDXC_Y],
                                                                                     particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDXC_Z],
                                                                                     particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions + IDXC_X],
                                                                                     particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions + IDXC_Y],
                                                                                     particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions + IDXC_Z],
-                                                                                    0, 0, 0);
+                                                                                    0, 0, 0,
+                                                                                    diff_x, diff_y, diff_z);
                                     if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
                                         computer_thread.template compute_interaction<size_particle_positions,size_particle_rhs>(
                                                             inout_index_particles[(intervals[idx_1].first+idx_p1)],
@@ -800,7 +811,7 @@ public:
                                                             inout_index_particles[(intervals[idx_2].first+idx_p2)],
                                                             &particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions],
                                                             &particles_current_rhs[0][(intervals[idx_2].first+idx_p2)*size_particle_rhs],
-                                                            dist_r2, cutoff_radius_compute, 0, 0, 0);
+                                                            dist_r2, cutoff_radius_compute, diff_x, diff_y, diff_z);
                                     }
                                 }
                             }
@@ -821,13 +832,15 @@ public:
                                 for(size_t idx_2 = 0 ; idx_2 < (*neighbors[idx_neighbor]).size() ; ++idx_2){
                                     for(partsize_t idx_p1 = 0 ; idx_p1 < intervals[idx_1].second ; ++idx_p1){
                                         for(partsize_t idx_p2 = 0 ; idx_p2 < (*neighbors[idx_neighbor])[idx_2].second ; ++idx_p2){
+                                            real_number diff_x, diff_y, diff_z;
                                             const real_number dist_r2 = compute_distance_r2(particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDXC_X],
                                                                                             particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDXC_Y],
                                                                                             particles_positions[(intervals[idx_1].first+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]);
+                                                                                            shift[idx_neighbor][IDXC_X], shift[idx_neighbor][IDXC_Y], shift[idx_neighbor][IDXC_Z],
+                                                                                            diff_x, diff_y, diff_z);
                                             if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
                                                 computer_thread.template compute_interaction<size_particle_positions,size_particle_rhs>(
                                                                     inout_index_particles[(intervals[idx_1].first+idx_p1)],
@@ -836,7 +849,8 @@ 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[0][((*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,
+                                                                    diff_x, diff_y, diff_z);
                                             }
                                         }
                                     }
diff --git a/cpp/particles/p2p/p2p_ghost_collisions.hpp b/cpp/particles/p2p/p2p_ghost_collisions.hpp
index b89fe8574d20620cfca8320b4f9fbebc0ef85405..13a7fb842b010f29eceb74be803738f35cf27f3b 100644
--- a/cpp/particles/p2p/p2p_ghost_collisions.hpp
+++ b/cpp/particles/p2p/p2p_ghost_collisions.hpp
@@ -216,11 +216,11 @@ public:
                              const partsize_t idx_part2,
                              const real_number pos_part2[],
                              real_number /*rhs_part2*/[],
-                             const real_number /*dist_pow2*/,
+                             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 xseparation,
+                             const real_number yseparation,
+                             const real_number zseparation){
         switch(this->current_particle_shape)
         {
             case SPHERE:
@@ -233,19 +233,14 @@ public:
                 break;
             case CYLINDER:
                 {
-                    double x, y, z, pq, xp, xq, t, s, x_dist, y_dist, z_dist ,min_distance;
+                    double pq, xp, xq, t, s, x_dist, y_dist, z_dist ,min_distance;
 
-                    /* Relative position vector of the two particles (x,y,z)^T */
-                    /* complicated usage of fmod, fabs and sgn because C++ fmod does not do what it should. */
-                    x = std::fmod(pos_part2[0],pi2)-std::fmod(pos_part1[0],pi2);
-                    y = std::fmod(pos_part2[1],pi2)-std::fmod(pos_part1[1],pi2);
-                    z = std::fmod(pos_part2[2],pi2)-std::fmod(pos_part1[2],pi2);
+                    const double x = xseparation;
+                    const double y = yseparation;
+                    const double z = zseparation;
+                    const double r = sqrt(dist_pow2);
 
-                    x = ( std::fmod( std::fabs(x) +pi ,pi2) - pi ) * sgn(x) ;
-                    y = ( std::fmod( std::fabs(y) +pi ,pi2) - pi ) * sgn(y) ;
-                    z = ( std::fmod( std::fabs(z) +pi ,pi2) - pi ) * sgn(z) ;
-
-                    if( sqrt(x*x+y*y+z*z) <= this->cylinder_length )
+                    if( r <= this->cylinder_length )
                     {
 
                         /* p and q are the orientation vectors of the first and second particles. */
@@ -384,7 +379,9 @@ public:
                     D1 = sqrt(x1x1+x0x0+t*t-2.0*x0x1-2.0*x1p0*t);
                     D2 = sqrt(x2x2+x0x0+t*t-2.0*x0x2-2.0*x2p0*t);
                     /* Check if a collision is possible. Exit if true. */
-                    if( D1<=this->disk_width and D2<=this->disk_width ){
+                    if ((D1 <= (this->disk_width/2)) and
+                        (D2 <= (this->disk_width/2)))
+                    {
                         std::pair <partsize_t, partsize_t> single_collision_pair(idx_part1, idx_part2);
                         this->collision_pairs_local.insert(single_collision_pair);
                         assert(idx_part1!=idx_part2);
@@ -397,7 +394,9 @@ public:
                     D1 = sqrt(x1x1+x0x0+t*t-2.0*x0x1-2.0*x1p0*t);
                     D2 = sqrt(x2x2+x0x0+t*t-2.0*x0x2-2.0*x2p0*t);
                     /* Check if a collision is possible. Exit if true. */
-                    if( D1<=this->disk_width and D2<=this->disk_width ){
+                    if ((D1 <= (this->disk_width/2)) and
+                        (D2 <= (this->disk_width/2)))
+                    {
                         std::pair <partsize_t, partsize_t> single_collision_pair(idx_part1, idx_part2);
                         this->collision_pairs_local.insert(single_collision_pair);
                         assert(idx_part1!=idx_part2);
diff --git a/cpp/particles/p2p/p2p_merge_collisions.hpp b/cpp/particles/p2p/p2p_merge_collisions.hpp
index 96fff9c2d7f6fa9f23d21ff346f1b8b819c9c1d5..cf93f3e8001699b3cb8de2db4c043409bae4861b 100644
--- a/cpp/particles/p2p/p2p_merge_collisions.hpp
+++ b/cpp/particles/p2p/p2p_merge_collisions.hpp
@@ -47,11 +47,16 @@ 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 /*xseparation*/,
+                             const real_number /*yseparation*/,
+                             const real_number /*zseparation*/){
         mergedParticles.emplace_back(std::max(idx_part1,idx_part2));
     }