diff --git a/bfps/cpp/particles/particles_field_computer.hpp b/bfps/cpp/particles/particles_field_computer.hpp
index 2a835c0e7a064b3ed6a666b2e7b803af8bc86845..0a155367b48dee14e02bc157e926088acff06f3b 100644
--- a/bfps/cpp/particles/particles_field_computer.hpp
+++ b/bfps/cpp/particles/particles_field_computer.hpp
@@ -13,9 +13,10 @@ template <class partsize_t,
           class interpolator_class,
           class field_class,
           int interp_neighbours,
-          class positions_updater_class >
-class particles_field_computer : public abstract_particles_distr<partsize_t, real_number, 3,3,1> {
-    using Parent = abstract_particles_distr<partsize_t, real_number, 3,3,1>;
+          class positions_updater_class,
+          class size_particle_rhs>
+class particles_field_computer : public abstract_particles_distr<partsize_t, real_number, 3,size_particle_rhs,1> {
+    using Parent = abstract_particles_distr<partsize_t, real_number, 3,size_particle_rhs,1>;
 
     const std::array<int,3> field_grid_dim;
     const std::pair<int,int> current_partition_interval;
@@ -39,7 +40,7 @@ class particles_field_computer : public abstract_particles_distr<partsize_t, rea
                                    const partsize_t nb_particles) const final{
         // Set values to zero initialy
         std::fill(particles_current_rhs,
-                  particles_current_rhs+nb_particles*3,
+                  particles_current_rhs+nb_particles*size_particle_rhs,
                   0);
     }
 
@@ -135,9 +136,9 @@ class particles_field_computer : public abstract_particles_distr<partsize_t, rea
                             const ptrdiff_t tindex = field.get_rindex_from_global(idx_x_pbc, idx_y_pbc, idx_z_pbc);
 
                             // getValue does not necessary return real_number
-                            particles_current_rhs[idxPart*3+IDX_X] += real_number(field.rval(tindex,IDX_X))*coef;
-                            particles_current_rhs[idxPart*3+IDX_Y] += real_number(field.rval(tindex,IDX_Y))*coef;
-                            particles_current_rhs[idxPart*3+IDX_Z] += real_number(field.rval(tindex,IDX_Z))*coef;
+                            for(int idx_rhs_val = 0 ; idx_rhs_val < size_particle_rhs ; ++idx_rhs_val){
+                                particles_current_rhs[idxPart*size_particle_rhs+idx_rhs_val] += real_number(field.rval(tindex,idx_rhs_val))*coef;
+                            }
                         }
                     }
                 }
@@ -151,9 +152,9 @@ class particles_field_computer : public abstract_particles_distr<partsize_t, rea
         TIMEZONE("particles_field_computer::reduce_particles");
         // Simply sum values
         for(partsize_t idxPart = 0 ; idxPart < nb_particles ; ++idxPart){
-            particles_current_rhs[idxPart*3+IDX_X] += extra_particles_current_rhs[idxPart*3+IDX_X];
-            particles_current_rhs[idxPart*3+IDX_Y] += extra_particles_current_rhs[idxPart*3+IDX_Y];
-            particles_current_rhs[idxPart*3+IDX_Z] += extra_particles_current_rhs[idxPart*3+IDX_Z];
+            for(int idx_rhs_val = 0 ; idx_rhs_val < size_particle_rhs ; ++idx_rhs_val){
+                particles_current_rhs[idxPart*size_particle_rhs+idx_rhs_val] += extra_particles_current_rhs[idxPart*size_particle_rhs+idx_rhs_val];
+            }
         }
     }
 
diff --git a/bfps/cpp/particles/particles_system.hpp b/bfps/cpp/particles/particles_system.hpp
index 6769cb89f379d5f9bb431c9217af53e901549d3b..444f0cc766bbbd3d06214117b40afe2d83590da8 100644
--- a/bfps/cpp/particles/particles_system.hpp
+++ b/bfps/cpp/particles/particles_system.hpp
@@ -11,7 +11,8 @@
 #include "particles_adams_bashforth.hpp"
 #include "scope_timer.hpp"
 
-template <class partsize_t, class real_number, class field_rnumber, class field_class, class interpolator_class, int interp_neighbours>
+template <class partsize_t, class real_number, class field_rnumber, class field_class, class interpolator_class, int interp_neighbours,
+          class size_particle_rhs>
 class particles_system : public abstract_particles_system<partsize_t, real_number> {
     MPI_Comm mpi_com;
 
@@ -20,7 +21,7 @@ class particles_system : public abstract_particles_system<partsize_t, real_numbe
 
     interpolator_class interpolator;
 
-    particles_field_computer<partsize_t, real_number, interpolator_class, field_class, interp_neighbours, particles_adams_bashforth<partsize_t, real_number, 3,3>> computer;
+    particles_field_computer<partsize_t, real_number, interpolator_class, field_class, interp_neighbours, particles_adams_bashforth<partsize_t, real_number, 3, size_particle_rhs>, size_particle_rhs> computer;
 
     std::unique_ptr<partsize_t[]> current_my_nb_particles_per_partition;
     std::unique_ptr<partsize_t[]> current_offset_particles_for_partition;
@@ -88,9 +89,9 @@ public:
         [&](const partsize_t idx1, const partsize_t idx2){
             std::swap(my_particles_positions_indexes[idx1], my_particles_positions_indexes[idx2]);
             for(int idx_rhs = 0 ; idx_rhs < int(my_particles_rhs.size()) ; ++idx_rhs){
-                for(int idx_val = 0 ; idx_val < 3 ; ++idx_val){
-                    std::swap(my_particles_rhs[idx_rhs][idx1*3 + idx_val],
-                              my_particles_rhs[idx_rhs][idx2*3 + idx_val]);
+                for(int idx_val = 0 ; idx_val < size_particle_rhs ; ++idx_val){
+                    std::swap(my_particles_rhs[idx_rhs][idx1*size_particle_rhs + idx_val],
+                              my_particles_rhs[idx_rhs][idx2*size_particle_rhs + idx_val]);
                 }
             }
         });
@@ -142,7 +143,7 @@ public:
                 my_particles_rhs[idx_rhs] = std::move(my_particles_rhs[idx_rhs-1]);
             }
             my_particles_rhs[0] = std::move(next_current);
-            particles_utils::memzero(my_particles_rhs[0], 3*my_nb_particles);
+            particles_utils::memzero(my_particles_rhs[0], size_particle_rhs*my_nb_particles);
         }
     }
 
@@ -182,9 +183,9 @@ public:
             assert(std::isnan(my_particles_positions[idx_part*3+IDX_Z]) == false);
 
             for(int idx_rhs = 0 ; idx_rhs < my_particles_rhs.size() ; ++idx_rhs){
-                assert(std::isnan(my_particles_rhs[idx_rhs][idx_part*3+IDX_X]) == false);
-                assert(std::isnan(my_particles_rhs[idx_rhs][idx_part*3+IDX_Y]) == false);
-                assert(std::isnan(my_particles_rhs[idx_rhs][idx_part*3+IDX_Z]) == false);
+                for(int idx_rhs_val = 0 ; idx_rhs_val < size_particle_rhs ; ++idx_rhs_val){
+                    assert(std::isnan(my_particles_rhs[idx_rhs][idx_part*size_particle_rhs+idx_rhs_val]) == false);
+                }
             }
         }
     }
diff --git a/bfps/cpp/particles/particles_system_builder.hpp b/bfps/cpp/particles/particles_system_builder.hpp
index 6a61a5b340588fb22074ffa559f5947b7fbb9840..a3c0d09433fdd5fdd14d1b313b8d735c81a40aba 100644
--- a/bfps/cpp/particles/particles_system_builder.hpp
+++ b/bfps/cpp/particles/particles_system_builder.hpp
@@ -108,11 +108,11 @@ inline RetType evaluate(IterType1 value1, IterType2 value2, Args... args){
 ///
 //////////////////////////////////////////////////////////////////////////////
 
-template <class partsize_t, class field_rnumber, field_backend be, class particles_rnumber>
+template <class partsize_t, class field_rnumber, field_backend be, field_components fc, class particles_rnumber>
 struct particles_system_build_container {
     template <const int interpolation_size, const int spline_mode>
     static std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>> instanciate(
-             const field<field_rnumber, be, THREE>* fs_field, // (field object)
+             const field<field_rnumber, be, fc>* fs_field, // (field object)
              const kspace<be, SMOOTH>* fs_kk, // (kspace object, contains dkx, dky, dkz)
              const int nsteps, // to check coherency between parameters and hdf input file (nb rhs)
              const partsize_t nparticles, // to check coherency between parameters and hdf input file
@@ -192,8 +192,8 @@ struct particles_system_build_container {
         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];
 
         // Create the particles system
-        particles_system<partsize_t, particles_rnumber, field_rnumber, field<field_rnumber, be, THREE>, particles_generic_interp<particles_rnumber, interpolation_size,spline_mode>, interpolation_size>* part_sys
-         = new particles_system<partsize_t, particles_rnumber, field_rnumber, field<field_rnumber, be, THREE>, particles_generic_interp<particles_rnumber, interpolation_size,spline_mode>, interpolation_size>(field_grid_dim,
+        particles_system<partsize_t, particles_rnumber, field_rnumber, field<field_rnumber, be, fc>, particles_generic_interp<particles_rnumber, interpolation_size,spline_mode>, interpolation_size, ncomp(fc)>>* part_sys
+         = new particles_system<partsize_t, particles_rnumber, field_rnumber, field<field_rnumber, be, fc>, particles_generic_interp<particles_rnumber, interpolation_size,spline_mode, interpolation_size, ncomp(fc)>>(field_grid_dim,
                                                                                                    spatial_box_width,
                                                                                                    spatial_box_offset,
                                                                                                    spatial_partition_width,
@@ -231,9 +231,9 @@ struct particles_system_build_container {
 };
 
 
-template <class partsize_t, class field_rnumber, field_backend be, class particles_rnumber = double>
+template <class partsize_t, class field_rnumber, field_backend be, field_components fc, class particles_rnumber = double>
 inline std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>> particles_system_builder(
-        const field<field_rnumber, be, THREE>* fs_field, // (field object)
+        const field<field_rnumber, be, fc>* fs_field, // (field object)
         const kspace<be, SMOOTH>* fs_kk, // (kspace object, contains dkx, dky, dkz)
         const int nsteps, // to check coherency between parameters and hdf input file (nb rhs)
         const partsize_t nparticles, // to check coherency between parameters and hdf input file
@@ -246,7 +246,7 @@ inline std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>
     return Template_double_for_if::evaluate<std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>,
                        int, 1, 11, 1, // interpolation_size
                        int, 0, 3, 1, // spline_mode
-                       particles_system_build_container<partsize_t, field_rnumber,be,particles_rnumber>>(
+                       particles_system_build_container<partsize_t, field_rnumber,be,fc,particles_rnumber>>(
                            interpolation_size, // template iterator 1
                            spline_mode, // template iterator 2
                            fs_field,fs_kk, nsteps, nparticles, fname_input, inDatanameState, inDatanameRhs, mpi_comm, in_current_iteration);