diff --git a/bfps/cpp/full_code/NSVEparticles.cpp b/bfps/cpp/full_code/NSVEparticles.cpp
index 953ad9f9316ed9e515be60f821ca0d49a262b63d..5f9f480de262753e132e0994fffd747dd6241e48 100644
--- a/bfps/cpp/full_code/NSVEparticles.cpp
+++ b/bfps/cpp/full_code/NSVEparticles.cpp
@@ -21,13 +21,13 @@ int NSVEparticles<rnumber>::initialize(void)
                 this->comm,
                 this->fs->iteration+1);
     this->particles_output_writer_mpi = new particles_output_hdf5<
-        long long int, double, 3, 3>(
+        long long int, double, 3>(
                 MPI_COMM_WORLD,
                 "tracers0",
                 nparticles,
                 tracers0_integration_steps);
     this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
-        long long int, double, 3, 3>(
+        long long int, double, 3>(
                 MPI_COMM_WORLD,
                 this->ps->getGlobalNbParticles(),
                 (this->simname + "_particles.h5"),
@@ -51,7 +51,7 @@ int NSVEparticles<rnumber>::write_checkpoint(void)
 {
     this->NSVE<rnumber>::write_checkpoint();
     this->particles_output_writer_mpi->open_file(this->fs->get_current_fname());
-    this->particles_output_writer_mpi->save(
+    this->particles_output_writer_mpi->template save<3>(
             this->ps->getParticlesState(),
             this->ps->getParticlesRhs(),
             this->ps->getParticlesIndexes(),
@@ -93,7 +93,7 @@ int NSVEparticles<rnumber>::do_stats()
     std::copy(this->ps->getParticlesState(),
               this->ps->getParticlesState()+3*this->ps->getLocalNbParticles(),
               pdata.get());
-    this->particles_sample_writer_mpi->save_dataset(
+    this->particles_sample_writer_mpi->template save_dataset<3>(
             "tracers0",
             "position",
             this->ps->getParticlesState(),
@@ -104,7 +104,7 @@ int NSVEparticles<rnumber>::do_stats()
 
     /// sample velocity
     this->ps->sample_compute_field(*this->tmp_vec_field, pdata.get());
-    this->particles_sample_writer_mpi->save_dataset(
+    this->particles_sample_writer_mpi->template save_dataset<3>(
             "tracers0",
             "velocity",
             this->ps->getParticlesState(),
@@ -117,7 +117,7 @@ int NSVEparticles<rnumber>::do_stats()
     this->fs->compute_Lagrangian_acceleration(this->tmp_vec_field);
     this->tmp_vec_field->ift();
     this->ps->sample_compute_field(*this->tmp_vec_field, pdata.get());
-    this->particles_sample_writer_mpi->save_dataset(
+    this->particles_sample_writer_mpi->template save_dataset<3>(
             "tracers0",
             "acceleration",
             this->ps->getParticlesState(),
diff --git a/bfps/cpp/full_code/NSVEparticles.hpp b/bfps/cpp/full_code/NSVEparticles.hpp
index 97ea5c8429924d0d4c20113bc8edd1cf3153e20d..1a2a65351c48ef7f56cb5cc7415538ed083bff67 100644
--- a/bfps/cpp/full_code/NSVEparticles.hpp
+++ b/bfps/cpp/full_code/NSVEparticles.hpp
@@ -59,8 +59,8 @@ class NSVEparticles: public NSVE<rnumber>
         /* other stuff */
         std::unique_ptr<abstract_particles_system<long long int, double>> ps;
 
-        particles_output_hdf5<long long int, double,3,3> *particles_output_writer_mpi;
-        particles_output_sampling_hdf5<long long int, double, 3, 3> *particles_sample_writer_mpi;
+        particles_output_hdf5<long long int, double,3> *particles_output_writer_mpi;
+        particles_output_sampling_hdf5<long long int, double, 3> *particles_sample_writer_mpi;
 
 
         NSVEparticles(
diff --git a/bfps/cpp/full_code/NSVEparticlesP2P.cpp b/bfps/cpp/full_code/NSVEparticlesP2P.cpp
index ff1a594ffcd450a975ab3cc0627a3d1182d11995..a6b7082e6b2ffa16e456f32d439e8b4cfecf451d 100644
--- a/bfps/cpp/full_code/NSVEparticlesP2P.cpp
+++ b/bfps/cpp/full_code/NSVEparticlesP2P.cpp
@@ -38,13 +38,13 @@ int NSVEparticlesP2P<rnumber>::initialize(void)
                 cutoff);
 
     this->particles_output_writer_mpi = new particles_output_hdf5<
-        long long int, double, 6, 6>(
+        long long int, double, 6>(
                 MPI_COMM_WORLD,
                 "tracers0",
                 nparticles,
                 tracers0_integration_steps);
     this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
-        long long int, double, 3, 3>(
+        long long int, double, 3>(
                 MPI_COMM_WORLD,
                 this->ps->getGlobalNbParticles(),
                 (this->simname + "_particles.h5"),
@@ -79,7 +79,7 @@ int NSVEparticlesP2P<rnumber>::write_checkpoint(void)
     this->NSVE<rnumber>::write_checkpoint();
     this->particles_output_writer_mpi->open_file(this->fs->get_current_fname());
     // TODO P2P write particle data too
-    this->particles_output_writer_mpi->save(
+    this->particles_output_writer_mpi->template save<6>(
             this->ps->getParticlesState(),
             this->ps->getParticlesRhs(),
             this->ps->getParticlesIndexes(),
@@ -117,7 +117,7 @@ int NSVEparticlesP2P<rnumber>::do_stats()
     std::unique_ptr<double[]> pdata1 = this->ps->extractParticlesState(3, 6);
 
     /// sample position
-    this->particles_sample_writer_mpi->save_dataset(
+    this->particles_sample_writer_mpi->template save_dataset<3>(
             "tracers0",
             "position",
             pdata0.get(),
@@ -127,7 +127,7 @@ int NSVEparticlesP2P<rnumber>::do_stats()
             this->ps->get_step_idx()-1);
 
     /// sample orientation
-    this->particles_sample_writer_mpi->save_dataset(
+    this->particles_sample_writer_mpi->template save_dataset<3>(
             "tracers0",
             "orientation",
             pdata0.get(),
@@ -138,7 +138,7 @@ int NSVEparticlesP2P<rnumber>::do_stats()
 
     /// sample velocity
     this->ps->sample_compute_field(*this->tmp_vec_field, pdata1.get());
-    this->particles_sample_writer_mpi->save_dataset(
+    this->particles_sample_writer_mpi->template save_dataset<3>(
             "tracers0",
             "velocity",
             pdata0.get(),
@@ -151,7 +151,7 @@ int NSVEparticlesP2P<rnumber>::do_stats()
     *this->tmp_vec_field = this->fs->cvorticity->get_cdata();
     this->tmp_vec_field->ift();
     this->ps->sample_compute_field(*this->tmp_vec_field, pdata1.get());
-    this->particles_sample_writer_mpi->save_dataset(
+    this->particles_sample_writer_mpi->template save_dataset<3>(
             "tracers0",
             "vorticity",
             pdata0.get(),
@@ -164,7 +164,7 @@ int NSVEparticlesP2P<rnumber>::do_stats()
     this->fs->compute_Lagrangian_acceleration(this->tmp_vec_field);
     this->tmp_vec_field->ift();
     this->ps->sample_compute_field(*this->tmp_vec_field, pdata1.get());
-    this->particles_sample_writer_mpi->save_dataset(
+    this->particles_sample_writer_mpi->template save_dataset<3>(
             "tracers0",
             "acceleration",
             pdata0.get(),
diff --git a/bfps/cpp/full_code/NSVEparticlesP2P.hpp b/bfps/cpp/full_code/NSVEparticlesP2P.hpp
index 8f435116d86f03bb2d805b0934a4a10371a754c0..b74169f5b35dabf66433cbab953e290fb6ceb943 100644
--- a/bfps/cpp/full_code/NSVEparticlesP2P.hpp
+++ b/bfps/cpp/full_code/NSVEparticlesP2P.hpp
@@ -65,8 +65,8 @@ class NSVEparticlesP2P: public NSVE<rnumber>
         /* other stuff */
         std::unique_ptr<abstract_particles_system<long long int, double>> ps;
         // TODO P2P use a reader with particle data
-        particles_output_hdf5<long long int, double,6,6> *particles_output_writer_mpi;
-        particles_output_sampling_hdf5<long long int, double,3,3> *particles_sample_writer_mpi;
+        particles_output_hdf5<long long int, double,6> *particles_output_writer_mpi;
+        particles_output_sampling_hdf5<long long int, double, 3> *particles_sample_writer_mpi;
 
 
         NSVEparticlesP2P(
diff --git a/bfps/cpp/particles/abstract_particles_output.hpp b/bfps/cpp/particles/abstract_particles_output.hpp
index 5285c90fe156ac92335141f3c75f02a057e9bbea..7510bc6e9990e8cb714580567da7e9ae9842bf87 100644
--- a/bfps/cpp/particles/abstract_particles_output.hpp
+++ b/bfps/cpp/particles/abstract_particles_output.hpp
@@ -13,7 +13,7 @@
 #include "scope_timer.hpp"
 #include "env_utils.hpp"
 
-template <class partsize_t, class real_number, int size_particle_positions, int size_particle_rhs>
+template <class partsize_t, class real_number, int size_particle_positions>
 class abstract_particles_output {
     MPI_Comm mpi_com;
     MPI_Comm mpi_com_writer;
@@ -145,6 +145,7 @@ public:
         }
     }
 
+    template <int size_particle_rhs>
     void save(
             const real_number input_particles_positions[],
             const std::unique_ptr<real_number[]> input_particles_rhs[],
@@ -264,11 +265,11 @@ public:
         }
 
         write(idx_time_step, buffer_particles_positions_send.get(), buffer_particles_rhs_send.data(),
-              nb_to_receive, particles_chunk_current_offset);
+              nb_to_receive, particles_chunk_current_offset, size_particle_rhs);
     }
 
     virtual void write(const int idx_time_step, const real_number* positions, const std::unique_ptr<real_number[]>* rhs,
-                       const partsize_t nb_particles, const partsize_t particles_idx_offset) = 0;
+                       const partsize_t nb_particles, const partsize_t particles_idx_offset, const int size_particle_rhs) = 0;
 };
 
 #endif
diff --git a/bfps/cpp/particles/particles_output_hdf5.hpp b/bfps/cpp/particles/particles_output_hdf5.hpp
index 647103ca9445dadb05c76d73710baea4e47cbec2..22d2fa851b0092f3eecfec456dc90d0e8e009169 100644
--- a/bfps/cpp/particles/particles_output_hdf5.hpp
+++ b/bfps/cpp/particles/particles_output_hdf5.hpp
@@ -10,16 +10,13 @@
 
 template <class partsize_t,
           class real_number,
-          int size_particle_positions,
-          int size_particle_rhs>
+          int size_particle_positions>
 class particles_output_hdf5 : public abstract_particles_output<partsize_t,
                                                                real_number,
-                                                               size_particle_positions,
-                                                               size_particle_rhs>{
+                                                               size_particle_positions>{
     using Parent = abstract_particles_output<partsize_t,
                                              real_number,
-                                             size_particle_positions,
-                                             size_particle_rhs>;
+                                             size_particle_positions>;
 
     std::string particle_species_name;
 
@@ -39,8 +36,7 @@ public:
                           const bool in_use_collective_io = false)
             : abstract_particles_output<partsize_t,
                                         real_number,
-                                        size_particle_positions,
-                                        size_particle_rhs>(
+                                        size_particle_positions>(
                                                 in_mpi_com,
                                                 inTotalNbParticles,
                                                 in_nb_rhs),
@@ -183,7 +179,8 @@ public:
             const real_number* particles_positions,
             const std::unique_ptr<real_number[]>* particles_rhs,
             const partsize_t nb_particles,
-            const partsize_t particles_idx_offset) final{
+            const partsize_t particles_idx_offset,
+            const int size_particle_rhs) final{
         assert(Parent::isInvolved());
 
         TIMEZONE("particles_output_hdf5::write");
diff --git a/bfps/cpp/particles/particles_output_mpiio.hpp b/bfps/cpp/particles/particles_output_mpiio.hpp
index 77dae6ca2f9441948ccf04f8a72e4a53d249894b..5810c4a06e795a32db13b6701ed4e0256db40eb3 100644
--- a/bfps/cpp/particles/particles_output_mpiio.hpp
+++ b/bfps/cpp/particles/particles_output_mpiio.hpp
@@ -11,8 +11,8 @@
 #include "particles_utils.hpp"
 
 template <class partsize_t, class real_number, int size_particle_positions, int size_particle_rhs>
-class particles_output_mpiio : public abstract_particles_output<partsize_t, real_number, size_particle_positions, size_particle_rhs>{
-    using Parent = abstract_particles_output<partsize_t, real_number, size_particle_positions, size_particle_rhs>;
+class particles_output_mpiio : public abstract_particles_output<partsize_t, real_number, size_particle_positions>{
+    using Parent = abstract_particles_output<partsize_t, real_number, size_particle_positions>;
 
     const std::string filename;
     const int nb_step_prealloc;
@@ -24,7 +24,7 @@ class particles_output_mpiio : public abstract_particles_output<partsize_t, real
 public:
     particles_output_mpiio(MPI_Comm in_mpi_com, const std::string in_filename, const partsize_t inTotalNbParticles,
                            const int in_nb_rhs, const int in_nb_step_prealloc = -1)
-            : abstract_particles_output<partsize_t, real_number, size_particle_positions, size_particle_rhs>(in_mpi_com, inTotalNbParticles, in_nb_rhs),
+            : abstract_particles_output<partsize_t, real_number, size_particle_positions>(in_mpi_com, inTotalNbParticles, in_nb_rhs),
               filename(in_filename), nb_step_prealloc(in_nb_step_prealloc), current_step_in_file(0){
         if(Parent::isInvolved()){
             {
diff --git a/bfps/cpp/particles/particles_output_sampling_hdf5.hpp b/bfps/cpp/particles/particles_output_sampling_hdf5.hpp
index 823754a5c7f880c215396553f7b0ebe6395be700..dc21322668420ac661a74c656b47ddb9d2e67e12 100644
--- a/bfps/cpp/particles/particles_output_sampling_hdf5.hpp
+++ b/bfps/cpp/particles/particles_output_sampling_hdf5.hpp
@@ -7,17 +7,14 @@
 
 template <class partsize_t,
           class real_number,
-          int size_particle_positions,
-          int size_particle_rhs>
+          int size_particle_positions>
 class particles_output_sampling_hdf5 : public abstract_particles_output<
                                        partsize_t,
                                        real_number,
-                                       size_particle_positions,
-                                       size_particle_rhs>{
+                                       size_particle_positions>{
     using Parent = abstract_particles_output<partsize_t,
                                              real_number,
-                                             size_particle_positions,
-                                             size_particle_rhs>;
+                                             size_particle_positions>;
 
     hid_t file_id, pgroup_id;
 
@@ -121,6 +118,7 @@ public:
         return EXIT_SUCCESS;
     }
 
+    template <int size_particle_rhs>
     int save_dataset(
             const std::string& in_groupname,
             const std::string& in_dataset_name,
@@ -144,7 +142,7 @@ public:
                 H5P_DEFAULT);
         AssertMpi(MPI_Bcast(&dataset_exists, 1, MPI_INT, 0, this->getCom()));
         if (dataset_exists == 0)
-            this->save(
+            this->template save<size_particle_rhs>(
                 input_particles_positions,
                 input_particles_rhs,
                 index_particles,
@@ -158,7 +156,8 @@ public:
             const real_number* /*particles_positions*/,
             const std::unique_ptr<real_number[]>* particles_rhs,
             const partsize_t nb_particles,
-            const partsize_t particles_idx_offset) final{
+            const partsize_t particles_idx_offset,
+            const int size_particle_rhs) final{
         assert(Parent::isInvolved());
 
         TIMEZONE("particles_output_hdf5::write");
diff --git a/bfps/cpp/particles/particles_sampling.hpp b/bfps/cpp/particles/particles_sampling.hpp
index a8927591c12f26a3a674d52271c4ac1a6ba51455..8baff633c155adb161c3764eaac8c9010c5fe086 100644
--- a/bfps/cpp/particles/particles_sampling.hpp
+++ b/bfps/cpp/particles/particles_sampling.hpp
@@ -21,10 +21,10 @@ void sample_from_particles_system(const field<rnumber, be, fc>& in_field, // a p
     const int size_particle_rhs = ncomp(fc);
 
     // Stop here if already exists
-    if(particles_output_sampling_hdf5<partsize_t, particles_rnumber, 3, size_particle_rhs>::DatasetExistsCol(MPI_COMM_WORLD,
-                                                                                                             filename,
-                                                                                                             parent_groupname,
-                                                                                                             datasetname)){
+    if(particles_output_sampling_hdf5<partsize_t, particles_rnumber, 3>::DatasetExistsCol(MPI_COMM_WORLD,
+                                                                                          filename,
+                                                                                          parent_groupname,
+                                                                                          datasetname)){
         return;
     }
 
@@ -36,12 +36,12 @@ void sample_from_particles_system(const field<rnumber, be, fc>& in_field, // a p
 
 
 
-    particles_output_sampling_hdf5<partsize_t, particles_rnumber, 3, size_particle_rhs> outputclass(MPI_COMM_WORLD,
-                                                                                                    ps->getGlobalNbParticles(),
-                                                                                                    filename,
-                                                                                                    parent_groupname,
-                                                                                                    datasetname);
-    outputclass.save(ps->getParticlesState(),
+    particles_output_sampling_hdf5<partsize_t, particles_rnumber, 3> outputclass(MPI_COMM_WORLD,
+                                                                                 ps->getGlobalNbParticles(),
+                                                                                 filename,
+                                                                                 parent_groupname,
+                                                                                 datasetname);
+    outputclass.template save<size_particle_rhs>(ps->getParticlesState(),
                      &sample_rhs,
                      ps->getParticlesIndexes(),
                      ps->getLocalNbParticles(),
@@ -57,10 +57,10 @@ void sample_particles_system_position(
     const std::string datasetname = fname + std::string("/") + std::to_string(ps->get_step_idx());
 
     // Stop here if already exists
-    if(particles_output_sampling_hdf5<partsize_t, particles_rnumber, 3, 3>::DatasetExistsCol(MPI_COMM_WORLD,
-                                                                                             filename,
-                                                                                             parent_groupname,
-                                                                                             datasetname)){
+    if(particles_output_sampling_hdf5<partsize_t, particles_rnumber, 3>::DatasetExistsCol(MPI_COMM_WORLD,
+                                                                                          filename,
+                                                                                          parent_groupname,
+                                                                                          datasetname)){
         return;
     }
 
@@ -68,12 +68,12 @@ void sample_particles_system_position(
     std::unique_ptr<particles_rnumber[]> sample_rhs(new particles_rnumber[3*nb_particles]);
     std::copy(ps->getParticlesState(), ps->getParticlesState() + 3*nb_particles, sample_rhs.get());
 
-    particles_output_sampling_hdf5<partsize_t, particles_rnumber, 3, 3> outputclass(MPI_COMM_WORLD,
-                                                                                    ps->getGlobalNbParticles(),
-                                                                                    filename,
-                                                                                    parent_groupname,
-                                                                                    datasetname);
-    outputclass.save(ps->getParticlesState(),
+    particles_output_sampling_hdf5<partsize_t, particles_rnumber, 3> outputclass(MPI_COMM_WORLD,
+                                                                                 ps->getGlobalNbParticles(),
+                                                                                 filename,
+                                                                                 parent_groupname,
+                                                                                 datasetname);
+    outputclass.template save<3>(ps->getParticlesState(),
                      &sample_rhs,
                      ps->getParticlesIndexes(),
                      ps->getLocalNbParticles(),