diff --git a/bfps/cpp/particles/p2p_computer.hpp b/bfps/cpp/particles/p2p_computer.hpp
index 46328c3a045911ed1af0345bb6a0a5772739bda0..922d65d1526f2491e3f3bbc6ec3b875d48cbc268 100644
--- a/bfps/cpp/particles/p2p_computer.hpp
+++ b/bfps/cpp/particles/p2p_computer.hpp
@@ -62,15 +62,15 @@ public:
         ///     (4 / \tau) \sum_j W_\ell ( | x^i - x^j | ) (p^i \cdot p^j)p^j
         /// \f]
         ///
-        const double dot_product = (pos_part1[3+IDX_X]*pos_part2[3+IDX_X] +
-                              pos_part1[3+IDX_Y]*pos_part2[3+IDX_Y] +
-                              pos_part1[3+IDX_Z]*pos_part2[3+IDX_Z]);
-        rhs_part1[3+IDX_X] += pos_part2[3+IDX_X] * 4 * ww * dot_product;
-        rhs_part1[3+IDX_Y] += pos_part2[3+IDX_Y] * 4 * ww * dot_product;
-        rhs_part1[3+IDX_Z] += pos_part2[3+IDX_Z] * 4 * ww * dot_product;
-        rhs_part2[3+IDX_X] += pos_part1[3+IDX_X] * 4 * ww * dot_product;
-        rhs_part2[3+IDX_Y] += pos_part1[3+IDX_Y] * 4 * ww * dot_product;
-        rhs_part2[3+IDX_Z] += pos_part1[3+IDX_Z] * 4 * ww * dot_product;
+        const double dot_product = (pos_part1[3+IDXC_X]*pos_part2[3+IDXC_X] +
+                              pos_part1[3+IDXC_Y]*pos_part2[3+IDXC_Y] +
+                              pos_part1[3+IDXC_Z]*pos_part2[3+IDXC_Z]);
+        rhs_part1[3+IDXC_X] += pos_part2[3+IDXC_X] * 4 * ww * dot_product;
+        rhs_part1[3+IDXC_Y] += pos_part2[3+IDXC_Y] * 4 * ww * dot_product;
+        rhs_part1[3+IDXC_Z] += pos_part2[3+IDXC_Z] * 4 * ww * dot_product;
+        rhs_part2[3+IDXC_X] += pos_part1[3+IDXC_X] * 4 * ww * dot_product;
+        rhs_part2[3+IDXC_Y] += pos_part1[3+IDXC_Y] * 4 * ww * dot_product;
+        rhs_part2[3+IDXC_Z] += pos_part1[3+IDXC_Z] * 4 * ww * dot_product;
     }
 
     bool isEnable() const {
diff --git a/bfps/cpp/particles/p2p_distr_mpi.hpp b/bfps/cpp/particles/p2p_distr_mpi.hpp
index 740aba18a0d353a8fb96ca699481e748bdb5d5b8..2758aecf42c62828057c03757f0a75d645d35640 100644
--- a/bfps/cpp/particles/p2p_distr_mpi.hpp
+++ b/bfps/cpp/particles/p2p_distr_mpi.hpp
@@ -103,7 +103,7 @@ protected:
 
     static int foundGridFactor(const real_number in_cutoff_radius, const std::array<real_number,3>& in_spatial_box_width){
         int idx_factor = 1;
-        while(in_cutoff_radius <= in_spatial_box_width[IDX_Z]/real_number(idx_factor+1)){
+        while(in_cutoff_radius <= in_spatial_box_width[IDXC_Z]/real_number(idx_factor+1)){
             idx_factor += 1;
         }
         return idx_factor;
@@ -126,7 +126,7 @@ public:
             spatial_box_width(in_spatial_box_width), spatial_box_offset(in_spatial_box_offset),
             cutoff_radius_compute(in_cutoff_radius),
             nb_cells_factor(foundGridFactor(in_cutoff_radius, in_spatial_box_width)),
-            cutoff_radius(in_spatial_box_width[IDX_Z]/real_number(nb_cells_factor)){
+            cutoff_radius(in_spatial_box_width[IDXC_Z]/real_number(nb_cells_factor)){
 
         AssertMpi(MPI_Comm_rank(current_com, &my_rank));
         AssertMpi(MPI_Comm_size(current_com, &nb_processes));
@@ -154,11 +154,11 @@ public:
             assert(partition_interval_size_per_proc[idx_proc_involved] != 0);
         }
 
-        assert(int(field_grid_dim[IDX_Z]) == partition_interval_offset_per_proc[nb_processes_involved]);
+        assert(int(field_grid_dim[IDXC_Z]) == partition_interval_offset_per_proc[nb_processes_involved]);
 
-        nb_cell_levels[IDX_X] = nb_cells_factor;
-        nb_cell_levels[IDX_Y] = nb_cells_factor;
-        nb_cell_levels[IDX_Z] = nb_cells_factor;
+        nb_cell_levels[IDXC_X] = nb_cells_factor;
+        nb_cell_levels[IDXC_Y] = nb_cells_factor;
+        nb_cell_levels[IDXC_Z] = nb_cells_factor;
     }
 
     virtual ~p2p_distr_mpi(){}
@@ -174,25 +174,25 @@ public:
     }
 
     long int get_cell_coord_x_from_index(const long int index) const{
-        return index % nb_cell_levels[IDX_X];
+        return index % nb_cell_levels[IDXC_X];
     }
 
     long int get_cell_coord_y_from_index(const long int index) const{
-        return (index % (nb_cell_levels[IDX_X]*nb_cell_levels[IDX_Y]))
-                / nb_cell_levels[IDX_X];
+        return (index % (nb_cell_levels[IDXC_X]*nb_cell_levels[IDXC_Y]))
+                / nb_cell_levels[IDXC_X];
     }
 
     long int get_cell_coord_z_from_index(const long int index) const{
-        return index / (nb_cell_levels[IDX_X]*nb_cell_levels[IDX_Y]);
+        return index / (nb_cell_levels[IDXC_X]*nb_cell_levels[IDXC_Y]);
     }
 
     long int first_cell_level_proc(const int dest_proc) const{
-        const real_number field_section_width_z = spatial_box_width[IDX_Z]/real_number(field_grid_dim[IDX_Z]);
+        const real_number field_section_width_z = spatial_box_width[IDXC_Z]/real_number(field_grid_dim[IDXC_Z]);
         return static_cast<long int>((field_section_width_z*real_number(partition_interval_offset_per_proc[dest_proc]))/cutoff_radius);
     }
 
     long int last_cell_level_proc(const int dest_proc) const{
-        const real_number field_section_width_z = spatial_box_width[IDX_Z]/real_number(field_grid_dim[IDX_Z]);
+        const real_number field_section_width_z = spatial_box_width[IDXC_Z]/real_number(field_grid_dim[IDXC_Z]);
         const long int limite = static_cast<long int>((field_section_width_z*real_number(partition_interval_offset_per_proc[dest_proc+1])
                                      - std::numeric_limits<real_number>::epsilon())/cutoff_radius);
         if(static_cast<real_number>(limite)*cutoff_radius
@@ -202,7 +202,7 @@ public:
         return limite;
     }
 
-    real_number apply_pbc(real_number pos, IDXS_3D dim) const{
+    real_number apply_pbc(real_number pos, COMPONENT_3D dim) const{
         while( pos < spatial_box_offset[dim] ){
             pos += spatial_box_width[dim];
         }
@@ -214,32 +214,32 @@ public:
 
     std::array<long int,3> get_cell_coordinate(const real_number pos_x, const real_number pos_y,
                                                const real_number pos_z) const {
-        const real_number diff_x = apply_pbc(pos_x,IDX_X) - spatial_box_offset[IDX_X];
-        const real_number diff_y = apply_pbc(pos_y,IDX_Y) - spatial_box_offset[IDX_Y];
-        const real_number diff_z = apply_pbc(pos_z,IDX_Z) - spatial_box_offset[IDX_Z];
+        const real_number diff_x = apply_pbc(pos_x,IDXC_X) - spatial_box_offset[IDXC_X];
+        const real_number diff_y = apply_pbc(pos_y,IDXC_Y) - spatial_box_offset[IDXC_Y];
+        const real_number diff_z = apply_pbc(pos_z,IDXC_Z) - spatial_box_offset[IDXC_Z];
         std::array<long int,3> coord;
-        coord[IDX_X] = static_cast<long int>(diff_x/cutoff_radius);
-        coord[IDX_Y] = static_cast<long int>(diff_y/cutoff_radius);
-        coord[IDX_Z] = static_cast<long int>(diff_z/cutoff_radius);
+        coord[IDXC_X] = static_cast<long int>(diff_x/cutoff_radius);
+        coord[IDXC_Y] = static_cast<long int>(diff_y/cutoff_radius);
+        coord[IDXC_Z] = static_cast<long int>(diff_z/cutoff_radius);
         return coord;
     }
 
     long int get_cell_idx(const real_number pos_x, const real_number pos_y,
                           const real_number pos_z) const {
         std::array<long int,3> coord = get_cell_coordinate(pos_x, pos_y, pos_z);
-        return ((coord[IDX_Z]*nb_cell_levels[IDX_Y])+coord[IDX_Y])*nb_cell_levels[IDX_X]+coord[IDX_X];
+        return ((coord[IDXC_Z]*nb_cell_levels[IDXC_Y])+coord[IDXC_Y])*nb_cell_levels[IDXC_X]+coord[IDXC_X];
     }
 
     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,IDX_X)-apply_pbc(x2,IDX_X)+xshift_coef*spatial_box_width[IDX_X]);
+        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);
 
-        real_number diff_y = std::abs(apply_pbc(y1,IDX_X)-apply_pbc(y2,IDX_X)+yshift_coef*spatial_box_width[IDX_Y]);
+        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);
 
-        real_number diff_z = std::abs(apply_pbc(z1,IDX_X)-apply_pbc(z2,IDX_X)+zshift_coef*spatial_box_width[IDX_Z]);
+        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);
 
         return (diff_x*diff_x) + (diff_y*diff_y) + (diff_z*diff_z);
@@ -276,9 +276,9 @@ public:
             for(int idxPartition = 0 ; idxPartition < current_partition_size ; ++idxPartition){
                 #pragma omp parallel for schedule(static)
                 for(partsize_t idxPart = current_offset_particles_for_partition[idxPartition] ; idxPart < current_offset_particles_for_partition[idxPartition+1] ; ++idxPart ){
-                    particles_coord[idxPart] = get_cell_idx(particles_positions[(idxPart)*size_particle_positions + IDX_X],
-                                                                              particles_positions[(idxPart)*size_particle_positions + IDX_Y],
-                                                                              particles_positions[(idxPart)*size_particle_positions + IDX_Z]);
+                    particles_coord[idxPart] = get_cell_idx(particles_positions[(idxPart)*size_particle_positions + IDXC_X],
+                                                                              particles_positions[(idxPart)*size_particle_positions + IDXC_Y],
+                                                                              particles_positions[(idxPart)*size_particle_positions + IDXC_Z]);
                     assert(my_down_z_cell_level <= get_cell_coord_z_from_index(particles_coord[idxPart]));
                     assert(get_cell_coord_z_from_index(particles_coord[idxPart]) <= my_top_z_cell_level);
                 }
@@ -378,11 +378,11 @@ public:
             int dest_proc = (my_rank+1)%nb_processes_involved;
             while(dest_proc != my_rank
                   && (my_top_z_cell_level == first_cell_level_proc(dest_proc)
-                      || (my_top_z_cell_level+1)%nb_cell_levels[IDX_Z] == first_cell_level_proc(dest_proc))){
+                      || (my_top_z_cell_level+1)%nb_cell_levels[IDXC_Z] == first_cell_level_proc(dest_proc))){
                 // Find if we have to send 1 or 2 cell levels
                 int nb_levels_to_send = 1;
                 if(my_nb_cell_levels > 1 // I have more than one level
-                        && (my_top_z_cell_level-1+2)%nb_cell_levels[IDX_Z] <= last_cell_level_proc(dest_proc)){
+                        && (my_top_z_cell_level-1+2)%nb_cell_levels[IDXC_Z] <= last_cell_level_proc(dest_proc)){
                     nb_levels_to_send += 1;
                 }
 
@@ -400,11 +400,11 @@ public:
             int src_proc = (my_rank-1+nb_processes_involved)%nb_processes_involved;
             while(src_proc != my_rank
                   && (last_cell_level_proc(src_proc) == my_down_z_cell_level
-                      || (last_cell_level_proc(src_proc)+1)%nb_cell_levels[IDX_Z] == my_down_z_cell_level)){
+                      || (last_cell_level_proc(src_proc)+1)%nb_cell_levels[IDXC_Z] == my_down_z_cell_level)){
                 // Find if we have to send 1 or 2 cell levels
                 int nb_levels_to_recv = 1;
                 if(my_nb_cell_levels > 1 // I have more than one level
-                        && first_cell_level_proc(src_proc) <= (my_down_z_cell_level-1+2)%nb_cell_levels[IDX_Z]){
+                        && first_cell_level_proc(src_proc) <= (my_down_z_cell_level-1+2)%nb_cell_levels[IDXC_Z]){
                     nb_levels_to_recv += 1;
                 }
 
@@ -564,14 +564,14 @@ public:
                         // Compute
                         partsize_t idxPart = 0;
                         while(idxPart != NbParticlesToReceive){
-                            const long int current_cell_idx = get_cell_idx(descriptor.toCompute[idxPart*size_particle_positions + IDX_X],
-                                                                           descriptor.toCompute[idxPart*size_particle_positions + IDX_Y],
-                                                                           descriptor.toCompute[idxPart*size_particle_positions + IDX_Z]);
+                            const long int current_cell_idx = get_cell_idx(descriptor.toCompute[idxPart*size_particle_positions + IDXC_X],
+                                                                           descriptor.toCompute[idxPart*size_particle_positions + IDXC_Y],
+                                                                           descriptor.toCompute[idxPart*size_particle_positions + IDXC_Z]);
                             partsize_t nb_parts_in_cell = 1;
                             while(idxPart+nb_parts_in_cell != NbParticlesToReceive
-                                  && current_cell_idx == get_cell_idx(descriptor.toCompute[(idxPart+nb_parts_in_cell)*size_particle_positions + IDX_X],
-                                                                     descriptor.toCompute[(idxPart+nb_parts_in_cell)*size_particle_positions + IDX_Y],
-                                                                     descriptor.toCompute[(idxPart+nb_parts_in_cell)*size_particle_positions + IDX_Z])){
+                                  && current_cell_idx == get_cell_idx(descriptor.toCompute[(idxPart+nb_parts_in_cell)*size_particle_positions + IDXC_X],
+                                                                     descriptor.toCompute[(idxPart+nb_parts_in_cell)*size_particle_positions + IDXC_Y],
+                                                                     descriptor.toCompute[(idxPart+nb_parts_in_cell)*size_particle_positions + IDXC_Z])){
                                 nb_parts_in_cell += 1;
                             }
 
@@ -589,20 +589,20 @@ 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 + IDX_X],
-                                                                                                descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions + IDX_Y],
-                                                                                                descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions + IDX_Z],
-                                                                                                particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_X],
-                                                                                                particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Y],
-                                                                                                particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Z],
-                                                                                                shift[idx_neighbor][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_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){
                                                     in_computer.template compute_interaction<size_particle_positions, size_particle_rhs>(
                                                                         &descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions],
                                                                         &descriptor.results[(idxPart+idx_p1)*size_particle_rhs],
                                                                         &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][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_Z]);
+                                                                        dist_r2, cutoff_radius_compute, shift[idx_neighbor][IDXC_X], shift[idx_neighbor][IDXC_Y], shift[idx_neighbor][IDXC_Z]);
                                                 }
                                             }
                                         }
@@ -671,12 +671,12 @@ 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){
-                            const real_number dist_r2 = compute_distance_r2(particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_X],
-                                                                            particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Y],
-                                                                            particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Z],
-                                                                            particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions + IDX_X],
-                                                                            particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions + IDX_Y],
-                                                                            particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions + IDX_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);
                             if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
                                 in_computer.template compute_interaction<size_particle_positions,size_particle_rhs>(
@@ -693,12 +693,12 @@ 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){
-                                const real_number dist_r2 = compute_distance_r2(particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_X],
-                                                                                particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Y],
-                                                                                particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Z],
-                                                                                particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions + IDX_X],
-                                                                                particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions + IDX_Y],
-                                                                                particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions + IDX_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);
                                 if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
                                     in_computer.template compute_interaction<size_particle_positions,size_particle_rhs>(
@@ -727,20 +727,20 @@ 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){
-                                        const real_number dist_r2 = compute_distance_r2(particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_X],
-                                                                                        particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Y],
-                                                                                        particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Z],
-                                                                                        particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_X],
-                                                                                        particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Y],
-                                                                                        particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Z],
-                                                                                        shift[idx_neighbor][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_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]);
                                         if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
                                             in_computer.template compute_interaction<size_particle_positions,size_particle_rhs>(
                                                                 &particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions],
                                                                 &particles_current_rhs[(intervals[idx_1].first+idx_p1)*size_particle_rhs],
                                                                 &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][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_Z]);
+                                                                dist_r2, cutoff_radius_compute, shift[idx_neighbor][IDXC_X], shift[idx_neighbor][IDXC_Y], shift[idx_neighbor][IDXC_Z]);
                                         }
                                     }
                                 }
diff --git a/bfps/cpp/particles/p2p_tree.hpp b/bfps/cpp/particles/p2p_tree.hpp
index 3d92c4e5666a1b1d94e71a6ed920a5bc063cfcba..a4441543e4d4711b8f1f311431e6e982fb34913f 100644
--- a/bfps/cpp/particles/p2p_tree.hpp
+++ b/bfps/cpp/particles/p2p_tree.hpp
@@ -11,21 +11,21 @@ class p2p_tree{
     std::array<long int,3> nb_cell_levels;
 
     long int get_cell_coord_x_from_index(const long int index) const{
-        return index % nb_cell_levels[IDX_X];
+        return index % nb_cell_levels[IDXC_X];
     }
 
     long int get_cell_coord_y_from_index(const long int index) const{
-        return (index % (nb_cell_levels[IDX_X]*nb_cell_levels[IDX_Y]))
-                / nb_cell_levels[IDX_X];
+        return (index % (nb_cell_levels[IDXC_X]*nb_cell_levels[IDXC_Y]))
+                / nb_cell_levels[IDXC_X];
     }
 
     long int get_cell_coord_z_from_index(const long int index) const{
-        return index / (nb_cell_levels[IDX_X]*nb_cell_levels[IDX_Y]);
+        return index / (nb_cell_levels[IDXC_X]*nb_cell_levels[IDXC_Y]);
     }
 
     long int get_cell_idx(const long int idx_x, const long int idx_y,
                           const long int idx_z) const {
-        return (((idx_z*nb_cell_levels[IDX_Y])+idx_y)*nb_cell_levels[IDX_X])+idx_x;
+        return (((idx_z*nb_cell_levels[IDXC_Y])+idx_y)*nb_cell_levels[IDXC_X])+idx_x;
     }
 
 public:
@@ -61,11 +61,11 @@ public:
             long int neigh_x_pbc = neigh_x+idx_x;
             ShiftType shift_x = 0;
             if(neigh_x_pbc < 0){
-                neigh_x_pbc += nb_cell_levels[IDX_X];
+                neigh_x_pbc += nb_cell_levels[IDXC_X];
                 shift_x = 1;
             }
-            else if(nb_cell_levels[IDX_X] <= neigh_x_pbc){
-                neigh_x_pbc -= nb_cell_levels[IDX_X];
+            else if(nb_cell_levels[IDXC_X] <= neigh_x_pbc){
+                neigh_x_pbc -= nb_cell_levels[IDXC_X];
                 shift_x = -1;
             }
 
@@ -73,11 +73,11 @@ public:
                 long int neigh_y_pbc = neigh_y+idx_y;
                 ShiftType shift_y = 0;
                 if(neigh_y_pbc < 0){
-                    neigh_y_pbc += nb_cell_levels[IDX_Y];
+                    neigh_y_pbc += nb_cell_levels[IDXC_Y];
                     shift_y = 1;
                 }
-                else if(nb_cell_levels[IDX_Y] <= neigh_y_pbc){
-                    neigh_y_pbc -= nb_cell_levels[IDX_Y];
+                else if(nb_cell_levels[IDXC_Y] <= neigh_y_pbc){
+                    neigh_y_pbc -= nb_cell_levels[IDXC_Y];
                     shift_y = -1;
                 }
 
@@ -85,11 +85,11 @@ public:
                     long int neigh_z_pbc = neigh_z+idx_z;
                     ShiftType shift_z = 0;
                     if(neigh_z_pbc < 0){
-                        neigh_z_pbc += nb_cell_levels[IDX_Z];
+                        neigh_z_pbc += nb_cell_levels[IDXC_Z];
                         shift_z = 1;
                     }
-                    else if(nb_cell_levels[IDX_Z] <= neigh_z_pbc){
-                        neigh_z_pbc -= nb_cell_levels[IDX_Z];
+                    else if(nb_cell_levels[IDXC_Z] <= neigh_z_pbc){
+                        neigh_z_pbc -= nb_cell_levels[IDXC_Z];
                         shift_z = -1;
                     }
 
@@ -102,9 +102,9 @@ public:
                             output[nbNeighbors] = &(iter->second);
                             output_indexes[nbNeighbors] = idx_neigh;
 
-                            shift[nbNeighbors][IDX_X] = shift_x;
-                            shift[nbNeighbors][IDX_Y] = shift_y;
-                            shift[nbNeighbors][IDX_Z] = shift_z;
+                            shift[nbNeighbors][IDXC_X] = shift_x;
+                            shift[nbNeighbors][IDXC_Y] = shift_y;
+                            shift[nbNeighbors][IDXC_Z] = shift_z;
 
                             nbNeighbors += 1;
                         }
diff --git a/bfps/cpp/particles/particles_distr_mpi.hpp b/bfps/cpp/particles/particles_distr_mpi.hpp
index 251119be950ad8219ee611d858ca8ac6c9986dbe..b03749fdf9b8fb05e17648f3165e26b2e6bc4267 100644
--- a/bfps/cpp/particles/particles_distr_mpi.hpp
+++ b/bfps/cpp/particles/particles_distr_mpi.hpp
@@ -127,7 +127,7 @@ public:
             assert(partition_interval_size_per_proc[idx_proc_involved] != 0);
         }
 
-        assert(int(field_grid_dim[IDX_Z]) == partition_interval_offset_per_proc[nb_processes_involved]);
+        assert(int(field_grid_dim[IDXC_Z]) == partition_interval_offset_per_proc[nb_processes_involved]);
     }
 
     virtual ~particles_distr_mpi(){}
@@ -522,11 +522,11 @@ public:
             partsize_t partOffset = 0;
             for(int idxPartition = 0 ; idxPartition < current_partition_size ; ++idxPartition){
                 for(partsize_t idx = 0 ; idx < current_my_nb_particles_per_partition[idxPartition] ; ++idx){
-                    const int partition_level = in_computer.pbc_field_layer((*inout_positions_particles)[(idx+partOffset)*size_particle_positions+IDX_Z], IDX_Z);
+                    const int partition_level = in_computer.pbc_field_layer((*inout_positions_particles)[(idx+partOffset)*size_particle_positions+IDXC_Z], IDXC_Z);
                     variable_used_only_in_assert(partition_level);
                     assert(partition_level == current_partition_interval.first + idxPartition
-                           || partition_level == (current_partition_interval.first + idxPartition-1+int(field_grid_dim[IDX_Z]))%int(field_grid_dim[IDX_Z])
-                           || partition_level == (current_partition_interval.first + idxPartition+1)%int(field_grid_dim[IDX_Z]));
+                           || partition_level == (current_partition_interval.first + idxPartition-1+int(field_grid_dim[IDXC_Z]))%int(field_grid_dim[IDXC_Z])
+                           || partition_level == (current_partition_interval.first + idxPartition+1)%int(field_grid_dim[IDXC_Z]));
                 }
                 partOffset += current_my_nb_particles_per_partition[idxPartition];
             }
@@ -543,11 +543,11 @@ public:
         // Find particles outside my interval
         const partsize_t nbOutLower = particles_utils::partition_extra<partsize_t, size_particle_positions>(&(*inout_positions_particles)[0], current_my_nb_particles_per_partition[0],
                     [&](const real_number val[]){
-            const int partition_level = in_computer.pbc_field_layer(val[IDX_Z], IDX_Z);
+            const int partition_level = in_computer.pbc_field_layer(val[IDXC_Z], IDXC_Z);
             assert(partition_level == current_partition_interval.first
-                   || partition_level == (current_partition_interval.first-1+int(field_grid_dim[IDX_Z]))%int(field_grid_dim[IDX_Z])
-                   || partition_level == (current_partition_interval.first+1)%int(field_grid_dim[IDX_Z]));
-            const bool isLower = partition_level == (current_partition_interval.first-1+int(field_grid_dim[IDX_Z]))%int(field_grid_dim[IDX_Z]);
+                   || partition_level == (current_partition_interval.first-1+int(field_grid_dim[IDXC_Z]))%int(field_grid_dim[IDXC_Z])
+                   || partition_level == (current_partition_interval.first+1)%int(field_grid_dim[IDXC_Z]));
+            const bool isLower = partition_level == (current_partition_interval.first-1+int(field_grid_dim[IDXC_Z]))%int(field_grid_dim[IDXC_Z]);
             return isLower;
         },
                     [&](const partsize_t idx1, const partsize_t idx2){
@@ -569,11 +569,11 @@ public:
                     &(*inout_positions_particles)[(current_offset_particles_for_partition[current_partition_size-1]+offesetOutLow)*size_particle_positions],
                     myTotalNbParticles - (current_offset_particles_for_partition[current_partition_size-1]+offesetOutLow),
                     [&](const real_number val[]){
-            const int partition_level = in_computer.pbc_field_layer(val[IDX_Z], IDX_Z);
+            const int partition_level = in_computer.pbc_field_layer(val[IDXC_Z], IDXC_Z);
             assert(partition_level == (current_partition_interval.second-1)
-                   || partition_level == ((current_partition_interval.second-1)-1+int(field_grid_dim[IDX_Z]))%int(field_grid_dim[IDX_Z])
-                   || partition_level == ((current_partition_interval.second-1)+1)%int(field_grid_dim[IDX_Z]));
-            const bool isUpper = (partition_level == ((current_partition_interval.second-1)+1)%int(field_grid_dim[IDX_Z]));
+                   || partition_level == ((current_partition_interval.second-1)-1+int(field_grid_dim[IDXC_Z]))%int(field_grid_dim[IDXC_Z])
+                   || partition_level == ((current_partition_interval.second-1)+1)%int(field_grid_dim[IDXC_Z]));
+            const bool isUpper = (partition_level == ((current_partition_interval.second-1)+1)%int(field_grid_dim[IDXC_Z]));
             return !isUpper;
         },
                     [&](const partsize_t idx1, const partsize_t idx2){
@@ -822,7 +822,7 @@ public:
                                              myTotalNbParticles,current_partition_size,
                                              current_my_nb_particles_per_partition, current_offset_particles_for_partition.get(),
                                              [&](const real_number& z_pos){
-                const int partition_level = in_computer.pbc_field_layer(z_pos, IDX_Z);
+                const int partition_level = in_computer.pbc_field_layer(z_pos, IDXC_Z);
                 assert(current_partition_interval.first <= partition_level && partition_level < current_partition_interval.second);
                 return partition_level - current_partition_interval.first;
             },
@@ -845,7 +845,7 @@ public:
                     assert(current_my_nb_particles_per_partition[idxPartition] ==
                            current_offset_particles_for_partition[idxPartition+1] - current_offset_particles_for_partition[idxPartition]);
                     for(partsize_t idx = current_offset_particles_for_partition[idxPartition] ; idx < current_offset_particles_for_partition[idxPartition+1] ; ++idx){
-                        assert(in_computer.pbc_field_layer((*inout_positions_particles)[idx*size_particle_positions+IDX_Z], IDX_Z)-current_partition_interval.first == idxPartition);
+                        assert(in_computer.pbc_field_layer((*inout_positions_particles)[idx*size_particle_positions+IDXC_Z], IDXC_Z)-current_partition_interval.first == idxPartition);
                     }
                 }
             }
diff --git a/bfps/cpp/particles/particles_field_computer.hpp b/bfps/cpp/particles/particles_field_computer.hpp
index 330763c85cc15ddddd08c2b2f99ff57d905ed859..50d4df7818bba1df96e10f0704f4ad995c65a133 100644
--- a/bfps/cpp/particles/particles_field_computer.hpp
+++ b/bfps/cpp/particles/particles_field_computer.hpp
@@ -44,9 +44,9 @@ public:
         : field_grid_dim({{int(in_field_grid_dim[0]),int(in_field_grid_dim[1]),int(in_field_grid_dim[2])}}), current_partition_interval(in_current_partitions),
           interpolator(in_interpolator),
           spatial_box_width(in_spatial_box_width), spatial_box_offset(in_spatial_box_offset), box_step_width(in_box_step_width){
-        deriv[IDX_X] = 0;
-        deriv[IDX_Y] = 0;
-        deriv[IDX_Z] = 0;
+        deriv[IDXC_X] = 0;
+        deriv[IDXC_Y] = 0;
+        deriv[IDXC_Z] = 0;
     }
 
     ////////////////////////////////////////////////////////////////////////
@@ -82,25 +82,25 @@ public:
         TIMEZONE("particles_field_computer::apply_computation");
 
         for(partsize_t idxPart = 0 ; idxPart < nb_particles ; ++idxPart){
-            const real_number reltv_x = get_norm_pos_in_cell(particles_positions[idxPart*size_particle_positions+IDX_X], IDX_X);
-            const real_number reltv_y = get_norm_pos_in_cell(particles_positions[idxPart*size_particle_positions+IDX_Y], IDX_Y);
-            const real_number reltv_z = get_norm_pos_in_cell(particles_positions[idxPart*size_particle_positions+IDX_Z], IDX_Z);
+            const real_number reltv_x = get_norm_pos_in_cell(particles_positions[idxPart*size_particle_positions+IDXC_X], IDXC_X);
+            const real_number reltv_y = get_norm_pos_in_cell(particles_positions[idxPart*size_particle_positions+IDXC_Y], IDXC_Y);
+            const real_number reltv_z = get_norm_pos_in_cell(particles_positions[idxPart*size_particle_positions+IDXC_Z], IDXC_Z);
 
             typename interpolator_class::real_number
                 bx[interp_neighbours*2+2],
                 by[interp_neighbours*2+2],
                 bz[interp_neighbours*2+2];
-            interpolator.compute_beta(deriv[IDX_X], reltv_x, bx);
-            interpolator.compute_beta(deriv[IDX_Y], reltv_y, by);
-            interpolator.compute_beta(deriv[IDX_Z], reltv_z, bz);
+            interpolator.compute_beta(deriv[IDXC_X], reltv_x, bx);
+            interpolator.compute_beta(deriv[IDXC_Y], reltv_y, by);
+            interpolator.compute_beta(deriv[IDXC_Z], reltv_z, bz);
 
-            const int partGridIdx_x = pbc_field_layer(particles_positions[idxPart*size_particle_positions+IDX_X], IDX_X);
-            const int partGridIdx_y = pbc_field_layer(particles_positions[idxPart*size_particle_positions+IDX_Y], IDX_Y);
-            const int partGridIdx_z = pbc_field_layer(particles_positions[idxPart*size_particle_positions+IDX_Z], IDX_Z);
+            const int partGridIdx_x = pbc_field_layer(particles_positions[idxPart*size_particle_positions+IDXC_X], IDXC_X);
+            const int partGridIdx_y = pbc_field_layer(particles_positions[idxPart*size_particle_positions+IDXC_Y], IDXC_Y);
+            const int partGridIdx_z = pbc_field_layer(particles_positions[idxPart*size_particle_positions+IDXC_Z], IDXC_Z);
 
-            assert(0 <= partGridIdx_x && partGridIdx_x < int(field_grid_dim[IDX_X]));
-            assert(0 <= partGridIdx_y && partGridIdx_y < int(field_grid_dim[IDX_Y]));
-            assert(0 <= partGridIdx_z && partGridIdx_z < int(field_grid_dim[IDX_Z]));
+            assert(0 <= partGridIdx_x && partGridIdx_x < int(field_grid_dim[IDXC_X]));
+            assert(0 <= partGridIdx_y && partGridIdx_y < int(field_grid_dim[IDXC_Y]));
+            assert(0 <= partGridIdx_z && partGridIdx_z < int(field_grid_dim[IDXC_Z]));
 
             const int interp_limit_mx = partGridIdx_x-interp_neighbours;
             const int interp_limit_x = partGridIdx_x+interp_neighbours+1;
@@ -113,8 +113,8 @@ public:
             int nb_z_intervals;
 
             if((partGridIdx_z-interp_neighbours) < 0){
-                assert(partGridIdx_z+interp_neighbours+1 < int(field_grid_dim[IDX_Z]));
-                interp_limit_mz[0] = std::max(current_partition_interval.first, partGridIdx_z-interp_neighbours+int(field_grid_dim[IDX_Z]));
+                assert(partGridIdx_z+interp_neighbours+1 < int(field_grid_dim[IDXC_Z]));
+                interp_limit_mz[0] = std::max(current_partition_interval.first, partGridIdx_z-interp_neighbours+int(field_grid_dim[IDXC_Z]));
                 interp_limit_z[0] = current_partition_interval.second-1;
 
                 interp_limit_mz[1] = std::max(0, current_partition_interval.first);
@@ -122,12 +122,12 @@ public:
 
                 nb_z_intervals = 2;
             }
-            else if(int(field_grid_dim[IDX_Z]) <= (partGridIdx_z+interp_neighbours+1)){
+            else if(int(field_grid_dim[IDXC_Z]) <= (partGridIdx_z+interp_neighbours+1)){
                 interp_limit_mz[0] = std::max(current_partition_interval.first, partGridIdx_z-interp_neighbours);
-                interp_limit_z[0] = std::min(int(field_grid_dim[IDX_Z])-1,current_partition_interval.second-1);
+                interp_limit_z[0] = std::min(int(field_grid_dim[IDXC_Z])-1,current_partition_interval.second-1);
 
                 interp_limit_mz[1] = std::max(0, current_partition_interval.first);
-                interp_limit_z[1] = std::min(partGridIdx_z+interp_neighbours+1-int(field_grid_dim[IDX_Z]), current_partition_interval.second-1);
+                interp_limit_z[1] = std::min(partGridIdx_z+interp_neighbours+1-int(field_grid_dim[IDXC_Z]), current_partition_interval.second-1);
 
                 nb_z_intervals = 2;
             }
@@ -139,19 +139,19 @@ public:
 
             for(int idx_inter = 0 ; idx_inter < nb_z_intervals ; ++idx_inter){
                 for(int idx_z = interp_limit_mz[idx_inter] ; idx_z <= interp_limit_z[idx_inter] ; ++idx_z ){
-                    const int idx_z_pbc = (idx_z + field_grid_dim[IDX_Z])%field_grid_dim[IDX_Z];
+                    const int idx_z_pbc = (idx_z + field_grid_dim[IDXC_Z])%field_grid_dim[IDXC_Z];
                     assert(current_partition_interval.first <= idx_z_pbc && idx_z_pbc < current_partition_interval.second);
-                    assert(((idx_z+field_grid_dim[IDX_Z]-interp_limit_mz_bz)%field_grid_dim[IDX_Z]) < interp_neighbours*2+2);
+                    assert(((idx_z+field_grid_dim[IDXC_Z]-interp_limit_mz_bz)%field_grid_dim[IDXC_Z]) < interp_neighbours*2+2);
 
                     for(int idx_x = interp_limit_mx ; idx_x <= interp_limit_x ; ++idx_x ){
-                        const int idx_x_pbc = (idx_x + field_grid_dim[IDX_X])%field_grid_dim[IDX_X];
+                        const int idx_x_pbc = (idx_x + field_grid_dim[IDXC_X])%field_grid_dim[IDXC_X];
                         assert(idx_x-interp_limit_mx < interp_neighbours*2+2);
 
                         for(int idx_y = interp_limit_my ; idx_y <= interp_limit_y ; ++idx_y ){
-                            const int idx_y_pbc = (idx_y + field_grid_dim[IDX_Y])%field_grid_dim[IDX_Y];
+                            const int idx_y_pbc = (idx_y + field_grid_dim[IDXC_Y])%field_grid_dim[IDXC_Y];
                             assert(idx_y-interp_limit_my < interp_neighbours*2+2);
 
-                            const real_number coef = (bz[((idx_z+field_grid_dim[IDX_Z]-interp_limit_mz_bz)%field_grid_dim[IDX_Z])]
+                            const real_number coef = (bz[((idx_z+field_grid_dim[IDXC_Z]-interp_limit_mz_bz)%field_grid_dim[IDXC_Z])]
                                                     * by[idx_y-interp_limit_my]
                                                     * bx[idx_x-interp_limit_mx]);
 
diff --git a/bfps/cpp/particles/particles_input_hdf5.hpp b/bfps/cpp/particles/particles_input_hdf5.hpp
index 40fef3c48a734dfc2b449382da9831b6f5feeebb..5231872d2f22dce1e8422cac05a0f7b6b72f9d10 100644
--- a/bfps/cpp/particles/particles_input_hdf5.hpp
+++ b/bfps/cpp/particles/particles_input_hdf5.hpp
@@ -240,7 +240,7 @@ public:
                                                 &split_particles_positions[previousOffset*size_particle_positions],
                                                  partsize_t(load_splitter.getMySize())-previousOffset,
                                                  [&](const real_number val[]){
-                    const real_number shiftPos = val[IDX_Z]-spatial_box_offset;
+                    const real_number shiftPos = val[IDXC_Z]-spatial_box_offset;
                     const real_number nbRepeat = floor(shiftPos/spatial_box_width);
                     const real_number posInBox = shiftPos - (spatial_box_width*nbRepeat);
                     return posInBox < limitPartitionShifted;
diff --git a/bfps/cpp/particles/particles_system.hpp b/bfps/cpp/particles/particles_system.hpp
index f8688f61ab71d07058618cf13b16b62a6ea6c019..460383c382076bbd7fb685bf78188954249ac839 100644
--- a/bfps/cpp/particles/particles_system.hpp
+++ b/bfps/cpp/particles/particles_system.hpp
@@ -70,7 +70,7 @@ public:
                      particles_inner_computer_class in_computer_particules_inner,
                      const int in_current_iteration = 1)
         : mpi_com(in_mpi_com),
-          current_partition_interval({in_local_field_offset[IDX_Z], in_local_field_offset[IDX_Z] + in_local_field_dims[IDX_Z]}),
+          current_partition_interval({in_local_field_offset[IDXC_Z], in_local_field_offset[IDXC_Z] + in_local_field_dims[IDXC_Z]}),
           partition_interval_size(current_partition_interval.second - current_partition_interval.first),
           interpolator(),
           particles_distr(in_mpi_com, current_partition_interval,field_grid_dim),
@@ -100,7 +100,7 @@ public:
         my_nb_particles = particles_input.getLocalNbParticles();
 
         for(partsize_t idx_part = 0 ; idx_part < my_nb_particles ; ++idx_part){ // TODO remove me
-            const int partition_level = computer.pbc_field_layer(my_particles_positions[idx_part*size_particle_positions+IDX_Z], IDX_Z);
+            const int partition_level = computer.pbc_field_layer(my_particles_positions[idx_part*size_particle_positions+IDXC_Z], IDXC_Z);
             variable_used_only_in_assert(partition_level);
             assert(partition_level >= current_partition_interval.first);
             assert(partition_level < current_partition_interval.second);
@@ -109,7 +109,7 @@ public:
         particles_utils::partition_extra_z<partsize_t, size_particle_positions>(&my_particles_positions[0], my_nb_particles, partition_interval_size,
                                               current_my_nb_particles_per_partition.get(), current_offset_particles_for_partition.get(),
         [&](const real_number& z_pos){
-            const int partition_level = computer.pbc_field_layer(z_pos, IDX_Z);
+            const int partition_level = computer.pbc_field_layer(z_pos, IDXC_Z);
             assert(current_partition_interval.first <= partition_level && partition_level < current_partition_interval.second);
             return partition_level - current_partition_interval.first;
         },
@@ -128,7 +128,7 @@ public:
                 assert(current_my_nb_particles_per_partition[idxPartition] ==
                        current_offset_particles_for_partition[idxPartition+1] - current_offset_particles_for_partition[idxPartition]);
                 for(partsize_t idx = current_offset_particles_for_partition[idxPartition] ; idx < current_offset_particles_for_partition[idxPartition+1] ; ++idx){
-                    assert(computer.pbc_field_layer(my_particles_positions[idx*size_particle_positions+IDX_Z], IDX_Z)-current_partition_interval.first == idxPartition);
+                    assert(computer.pbc_field_layer(my_particles_positions[idx*size_particle_positions+IDXC_Z], IDXC_Z)-current_partition_interval.first == idxPartition);
                 }
             }
         }
@@ -331,9 +331,9 @@ public:
 
     void checkNan() const { // TODO remove
         for(partsize_t idx_part = 0 ; idx_part < my_nb_particles ; ++idx_part){ // TODO remove me
-            assert(std::isnan(my_particles_positions[idx_part*size_particle_positions+IDX_X]) == false);
-            assert(std::isnan(my_particles_positions[idx_part*size_particle_positions+IDX_Y]) == false);
-            assert(std::isnan(my_particles_positions[idx_part*size_particle_positions+IDX_Z]) == false);
+            assert(std::isnan(my_particles_positions[idx_part*size_particle_positions+IDXC_X]) == false);
+            assert(std::isnan(my_particles_positions[idx_part*size_particle_positions+IDXC_Y]) == false);
+            assert(std::isnan(my_particles_positions[idx_part*size_particle_positions+IDXC_Z]) == false);
 
             for(int idx_rhs = 0 ; idx_rhs < my_particles_rhs.size() ; ++idx_rhs){
                 for(int idx_rhs_val = 0 ; idx_rhs_val < size_particle_rhs ; ++idx_rhs_val){
diff --git a/bfps/cpp/particles/particles_system_builder.hpp b/bfps/cpp/particles/particles_system_builder.hpp
index 3e72ceaac18c68da3bef767738ba54d4b4a66d96..a9a140ac9921905513c1adcf558e071093b5afbb 100644
--- a/bfps/cpp/particles/particles_system_builder.hpp
+++ b/bfps/cpp/particles/particles_system_builder.hpp
@@ -130,21 +130,21 @@ struct particles_system_build_container {
 
         // The size of the field grid (global size) all_size seems
         std::array<size_t,3> field_grid_dim;
-        field_grid_dim[IDX_X] = fs_field->rlayout->sizes[FIELD_IDX_X];// nx
-        field_grid_dim[IDX_Y] = fs_field->rlayout->sizes[FIELD_IDX_Y];// nx
-        field_grid_dim[IDX_Z] = fs_field->rlayout->sizes[FIELD_IDX_Z];// nz
+        field_grid_dim[IDXC_X] = fs_field->rlayout->sizes[IDXV_X];// nx
+        field_grid_dim[IDXC_Y] = fs_field->rlayout->sizes[IDXV_Y];// nx
+        field_grid_dim[IDXC_Z] = fs_field->rlayout->sizes[IDXV_Z];// nz
 
         // The size of the local field grid (the field nodes that belong to current process)
         std::array<size_t,3> local_field_dims;
-        local_field_dims[IDX_X] = fs_field->rlayout->subsizes[FIELD_IDX_X];
-        local_field_dims[IDX_Y] = fs_field->rlayout->subsizes[FIELD_IDX_Y];
-        local_field_dims[IDX_Z] = fs_field->rlayout->subsizes[FIELD_IDX_Z];
+        local_field_dims[IDXC_X] = fs_field->rlayout->subsizes[IDXV_X];
+        local_field_dims[IDXC_Y] = fs_field->rlayout->subsizes[IDXV_Y];
+        local_field_dims[IDXC_Z] = fs_field->rlayout->subsizes[IDXV_Z];
 
         // The offset of the local field grid
         std::array<size_t,3> local_field_offset;
-        local_field_offset[IDX_X] = fs_field->rlayout->starts[FIELD_IDX_X];
-        local_field_offset[IDX_Y] = fs_field->rlayout->starts[FIELD_IDX_Y];
-        local_field_offset[IDX_Z] = fs_field->rlayout->starts[FIELD_IDX_Z];
+        local_field_offset[IDXC_X] = fs_field->rlayout->starts[IDXV_X];
+        local_field_offset[IDXC_Y] = fs_field->rlayout->starts[IDXV_Y];
+        local_field_offset[IDXC_Z] = fs_field->rlayout->starts[IDXV_Z];
 
 
         // Retreive split from fftw to know processes that have no work
@@ -152,51 +152,51 @@ struct particles_system_build_container {
         AssertMpi(MPI_Comm_rank(mpi_comm, &my_rank));
         AssertMpi(MPI_Comm_size(mpi_comm, &nb_processes));
 
-        const int split_step = (int(field_grid_dim[IDX_Z])+nb_processes-1)/nb_processes;
-        const int nb_processes_involved = (int(field_grid_dim[IDX_Z])+split_step-1)/split_step;
+        const int split_step = (int(field_grid_dim[IDXC_Z])+nb_processes-1)/nb_processes;
+        const int nb_processes_involved = (int(field_grid_dim[IDXC_Z])+split_step-1)/split_step;
 
-        assert((my_rank < nb_processes_involved && local_field_dims[IDX_Z] != 0)
-               || (nb_processes_involved <= my_rank && local_field_dims[IDX_Z] == 0));
-        assert(nb_processes_involved <= int(field_grid_dim[IDX_Z]));
+        assert((my_rank < nb_processes_involved && local_field_dims[IDXC_Z] != 0)
+               || (nb_processes_involved <= my_rank && local_field_dims[IDXC_Z] == 0));
+        assert(nb_processes_involved <= int(field_grid_dim[IDXC_Z]));
 
         // Make the idle processes starting from the limit (and not 0 as set by fftw)
         if(nb_processes_involved <= my_rank){
-            local_field_offset[IDX_Z] = field_grid_dim[IDX_Z];
+            local_field_offset[IDXC_Z] = field_grid_dim[IDXC_Z];
         }
 
         // Ensure that 1D partitioning is used
         {
-            assert(local_field_offset[IDX_X] == 0);
-            assert(local_field_offset[IDX_Y] == 0);
-            assert(local_field_dims[IDX_X] == field_grid_dim[IDX_X]);
-            assert(local_field_dims[IDX_Y] == field_grid_dim[IDX_Y]);
-
-            assert(my_rank >= nb_processes_involved || ((my_rank == 0 && local_field_offset[IDX_Z] == 0)
-                   || (my_rank != 0 && local_field_offset[IDX_Z] != 0)));
-            assert(my_rank >= nb_processes_involved || ((my_rank == nb_processes_involved-1 && local_field_offset[IDX_Z]+local_field_dims[IDX_Z] == field_grid_dim[IDX_Z])
-                   || (my_rank != nb_processes_involved-1 && local_field_offset[IDX_Z]+local_field_dims[IDX_Z] != field_grid_dim[IDX_Z])));
+            assert(local_field_offset[IDXC_X] == 0);
+            assert(local_field_offset[IDXC_Y] == 0);
+            assert(local_field_dims[IDXC_X] == field_grid_dim[IDXC_X]);
+            assert(local_field_dims[IDXC_Y] == field_grid_dim[IDXC_Y]);
+
+            assert(my_rank >= nb_processes_involved || ((my_rank == 0 && local_field_offset[IDXC_Z] == 0)
+                   || (my_rank != 0 && local_field_offset[IDXC_Z] != 0)));
+            assert(my_rank >= nb_processes_involved || ((my_rank == nb_processes_involved-1 && local_field_offset[IDXC_Z]+local_field_dims[IDXC_Z] == field_grid_dim[IDXC_Z])
+                   || (my_rank != nb_processes_involved-1 && local_field_offset[IDXC_Z]+local_field_dims[IDXC_Z] != field_grid_dim[IDXC_Z])));
         }
 
         // The spatial box size (all particles should be included inside)
         std::array<particles_rnumber,3> spatial_box_width;
-        spatial_box_width[IDX_X] = 4 * acos(0) / (fs_kk->dkx);
-        spatial_box_width[IDX_Y] = 4 * acos(0) / (fs_kk->dky);
-        spatial_box_width[IDX_Z] = 4 * acos(0) / (fs_kk->dkz);
+        spatial_box_width[IDXC_X] = 4 * acos(0) / (fs_kk->dkx);
+        spatial_box_width[IDXC_Y] = 4 * acos(0) / (fs_kk->dky);
+        spatial_box_width[IDXC_Z] = 4 * acos(0) / (fs_kk->dkz);
 
         // Box is in the corner
         std::array<particles_rnumber,3> spatial_box_offset;
-        spatial_box_offset[IDX_X] = 0;
-        spatial_box_offset[IDX_Y] = 0;
-        spatial_box_offset[IDX_Z] = 0;
+        spatial_box_offset[IDXC_X] = 0;
+        spatial_box_offset[IDXC_Y] = 0;
+        spatial_box_offset[IDXC_Z] = 0;
 
         // The distance between two field nodes in z
         std::array<particles_rnumber,3> spatial_partition_width;
-        spatial_partition_width[IDX_X] = spatial_box_width[IDX_X]/particles_rnumber(field_grid_dim[IDX_X]);
-        spatial_partition_width[IDX_Y] = spatial_box_width[IDX_Y]/particles_rnumber(field_grid_dim[IDX_Y]);
-        spatial_partition_width[IDX_Z] = spatial_box_width[IDX_Z]/particles_rnumber(field_grid_dim[IDX_Z]);
+        spatial_partition_width[IDXC_X] = spatial_box_width[IDXC_X]/particles_rnumber(field_grid_dim[IDXC_X]);
+        spatial_partition_width[IDXC_Y] = spatial_box_width[IDXC_Y]/particles_rnumber(field_grid_dim[IDXC_Y]);
+        spatial_partition_width[IDXC_Z] = spatial_box_width[IDXC_Z]/particles_rnumber(field_grid_dim[IDXC_Z]);
         // The spatial interval of the current process
-        const particles_rnumber my_spatial_low_limit_z = particles_rnumber(local_field_offset[IDX_Z])*spatial_partition_width[IDX_Z];
-        const particles_rnumber my_spatial_up_limit_z = particles_rnumber(local_field_offset[IDX_Z]+local_field_dims[IDX_Z])*spatial_partition_width[IDX_Z];
+        const particles_rnumber my_spatial_low_limit_z = particles_rnumber(local_field_offset[IDXC_Z])*spatial_partition_width[IDXC_Z];
+        const particles_rnumber my_spatial_up_limit_z = particles_rnumber(local_field_offset[IDXC_Z]+local_field_dims[IDXC_Z])*spatial_partition_width[IDXC_Z];
 
         // Create the particles system
         using particles_system_type = particles_system<partsize_t, particles_rnumber, field_rnumber,
diff --git a/bfps/cpp/particles/particles_utils.hpp b/bfps/cpp/particles/particles_utils.hpp
index 146dc4399477b72c30329edff587d35d7b44d69d..76371e64f7bebcbc5d54004f417461c6d211a558 100644
--- a/bfps/cpp/particles/particles_utils.hpp
+++ b/bfps/cpp/particles/particles_utils.hpp
@@ -19,16 +19,28 @@
 #define AssertMpi(X) if(MPI_SUCCESS != (X)) { printf("MPI Error at line %d\n",__LINE__); fflush(stdout) ; throw std::runtime_error("Stop from from mpi erro"); }
 #endif
 
-enum IDXS_3D {
-    IDX_X = 0,
-    IDX_Y = 1,
-    IDX_Z = 2
+enum IDX_COMPONENT_3D {
+    IDXC_X = 0,
+    IDXC_Y = 1,
+    IDXC_Z = 2
 };
 
-enum FIELD_IDXS_3D {
-    FIELD_IDX_X = 2,
-    FIELD_IDX_Y = 1,
-    FIELD_IDX_Z = 0
+enum IDX_COMPONENT_DEL_3D {
+    IDXC_DX_X = 0,
+    IDXC_DX_Y = 1,
+    IDXC_DX_Z = 2,
+    IDXC_DY_X = 3,
+    IDXC_DY_Y = 4,
+    IDXC_DY_Z = 5,
+    IDXC_DZ_X = 6,
+    IDXC_DZ_Y = 7,
+    IDXC_DZ_Z = 8,
+};
+
+enum IDX_VARIABLE_3D {
+    IDXV_X = 2,
+    IDXV_Y = 1,
+    IDXV_Z = 0
 };
 
 namespace particles_utils {
@@ -123,7 +135,7 @@ inline void partition_extra_z(real_number* array, const partsize_t size, const i
     if(nb_partitions == 2){
         const partsize_t size_current = partition_extra<partsize_t, nb_values>(array, size,
                 [&](const real_number inval[]){
-            return partitions_levels(inval[IDX_Z]) == 0;
+            return partitions_levels(inval[IDXC_Z]) == 0;
         }, pdcswap);
         partitions_size[0] = size_current;
         partitions_size[1] = size-size_current;
@@ -152,7 +164,7 @@ inline void partition_extra_z(real_number* array, const partsize_t size, const i
             const partsize_t size_current = partition_extra<partsize_t, nb_values>(&array[partitions_offset[current_part.first]*nb_values],
                                                      size_unpart,
                     [&](const real_number inval[]){
-                return partitions_levels(inval[IDX_Z]) <= idx_middle;
+                return partitions_levels(inval[IDXC_Z]) <= idx_middle;
             }, pdcswap, partitions_offset[current_part.first]);
 
             partitions_offset[idx_middle+1] = size_current + partitions_offset[current_part.first];