diff --git a/cpp/base.hpp b/cpp/base.hpp
index 203d707c53e267bd1555ef5fcab3989810408ee3..be05de0490194d34148f7afc888281a918169d8f 100644
--- a/cpp/base.hpp
+++ b/cpp/base.hpp
@@ -132,6 +132,33 @@ public:
     }
 };
 
+template <>
+class mpi_real_type<int>
+{
+    public:
+        static constexpr MPI_Datatype real(){
+            return MPI_INT;
+        }
+};
+
+template <>
+class mpi_real_type<long int>
+{
+    public:
+        static constexpr MPI_Datatype real(){
+            return MPI_LONG;
+        }
+};
+
+template <>
+class mpi_real_type<long long int>
+{
+    public:
+        static constexpr MPI_Datatype real(){
+            return MPI_LONG_LONG;
+        }
+};
+
 /////////////////////////////////////////////////////////////
 /////////////////////////////////////////////////////////////
 
diff --git a/cpp/hdf5_tools.cpp b/cpp/hdf5_tools.cpp
index 44da17c3789b5d14c8c21d6e626c2e7edc641fcb..57545fc75a3467865a01a8511ae19afad6604e50 100644
--- a/cpp/hdf5_tools.cpp
+++ b/cpp/hdf5_tools.cpp
@@ -412,6 +412,68 @@ int hdf5_tools::write_vector_with_single_rank(
     return EXIT_SUCCESS;
 }
 
+template <typename number>
+int hdf5_tools::gather_and_write_with_single_rank(
+        const int myrank,
+        const int rank_to_use,
+        const MPI_Comm COMM,
+        const number *data,
+        const int number_of_items,
+        const hid_t group,
+        const std::string dset_name)
+{
+    int nprocs;
+    MPI_Comm_size(COMM, &nprocs);
+    // first get number of items for each rank:
+    int *recvcounts = NULL;
+    if (myrank == rank_to_use)
+        recvcounts = new int[nprocs];
+    MPI_Gather(
+                &number_of_items,
+                1,
+                MPI_INT,
+                recvcounts,
+                1,
+                MPI_INT,
+                rank_to_use,
+                COMM);
+    // now prepare to gather data on rank_to_use
+    int total_number_of_items = 0;
+    int *displacement = NULL;
+    if (myrank == rank_to_use)
+    {
+        displacement = new int[nprocs];
+        displacement[0] = 0;
+        for(int ii = 1; ii < nprocs; ii++)
+        {
+            displacement[ii] = displacement[ii-1] + recvcounts[ii-1];
+        }
+        total_number_of_items = displacement[nprocs-1] + recvcounts[nprocs-1];
+    }
+    std::vector<number> data_to_write;
+    if (myrank == rank_to_use)
+        data_to_write.resize(total_number_of_items);
+    MPI_Gatherv(
+                data,
+                number_of_items,
+                mpi_real_type<number>::real(),
+                &data_to_write.front(),
+                recvcounts,
+                displacement,
+                mpi_real_type<number>::real(),
+                rank_to_use,
+                COMM);
+    if (myrank == rank_to_use)
+    {
+        write_vector_with_single_rank<number>(
+            data_to_write,
+            group,
+            dset_name);
+        delete[] recvcounts;
+    }
+    return EXIT_SUCCESS;
+}
+
 template
 std::vector<int> hdf5_tools::read_vector<int>(
         const hid_t,
@@ -477,3 +539,23 @@ int hdf5_tools::write_vector_with_single_rank(
         const hid_t group,
         const std::string dset_name);
 
+template
+int hdf5_tools::gather_and_write_with_single_rank<double>(
+        const int myrank,
+        const int rank_to_use,
+        const MPI_Comm COMM,
+        const double *data,
+        const int number_of_items,
+        const hid_t group,
+        const std::string dset_name);
+
+template
+int hdf5_tools::gather_and_write_with_single_rank<long long int>(
+        const int myrank,
+        const int rank_to_use,
+        const MPI_Comm COMM,
+        const long long int *data,
+        const int number_of_items,
+        const hid_t group,
+        const std::string dset_name);
+
diff --git a/cpp/hdf5_tools.hpp b/cpp/hdf5_tools.hpp
index 09f578283b7ad4ce185c027b8a825d62048eef19..b6068de6e4686e3fb478e28080f8235e4169d956 100644
--- a/cpp/hdf5_tools.hpp
+++ b/cpp/hdf5_tools.hpp
@@ -105,6 +105,16 @@ namespace hdf5_tools
             const std::vector<real_number> v,
             const hid_t group,
             const std::string dset_name);
+
+    template <typename number>
+    int gather_and_write_with_single_rank(
+            const int myrank,
+            const int rank_to_use,
+            const MPI_Comm COMM,
+            const number *data,
+            const int number_of_items,
+            const hid_t group,
+            const std::string dset_name);
 }
 
 #endif//HDF5_TOOLS_HPP
diff --git a/cpp/particles/interpolation/abstract_particle_set.hpp b/cpp/particles/interpolation/abstract_particle_set.hpp
index a5a8442e5b333414d933eb4e53c54d3691d3da67..621787b3d1fd9f8473e369cd320fcf3fd73fba4b 100644
--- a/cpp/particles/interpolation/abstract_particle_set.hpp
+++ b/cpp/particles/interpolation/abstract_particle_set.hpp
@@ -222,6 +222,38 @@ class abstract_particle_set
             return EXIT_SUCCESS;
         }
 
+        template <typename field_rnumber,
+                  field_backend be,
+                  field_components fc>
+        int writeSample(
+                field<field_rnumber, be, fc> *field_to_sample,
+                const std::string file_name,
+                const std::string species_name,
+                const std::string field_name,
+                const int iteration)
+        {
+            TIMEZONE("abstract_particle_set::writeSample_direct_to_file");
+            // allocate temporary array
+            std::unique_ptr<particle_rnumber[]> pdata(new particle_rnumber[ncomp(fc)*this->getLocalNumberOfParticles()]);
+            // clean up temporary array
+            set_particle_data_to_zero<partsize_t, particle_rnumber, ncomp(fc)>(
+                    pdata.get(), this->getLocalNumberOfParticles());
+            this->sample(*field_to_sample, pdata.get());
+            hid_t file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
+            hdf5_tools::gather_and_write_with_single_rank(
+                    field_to_sample->myrank,
+                    0,
+                    field_to_sample->comm,
+                    pdata.get(),
+                    this->getLocalNumberOfParticles(),
+                    file_id,
+                    species_name + std::string("/") + field_name + std::string("/") + std::to_string(iteration));
+            H5Fclose(file_id);
+            // deallocate temporary array
+            delete[] pdata.release();
+            return EXIT_SUCCESS;
+        }
+
         int writeStateChunk(
                 const int i0,
                 const int contiguous_state_chunk,
diff --git a/cpp/particles/interpolation/particle_set.hpp b/cpp/particles/interpolation/particle_set.hpp
index cccfd5764158421af2e3f82785c819eb18d0bbfc..206e9da519432a1226c441517f97d26c5f03c123 100644
--- a/cpp/particles/interpolation/particle_set.hpp
+++ b/cpp/particles/interpolation/particle_set.hpp
@@ -374,6 +374,7 @@ class particle_set: public abstract_particle_set
                     if (src_index == indices_to_copy[iii])
                     {
                         tmp_local_index[this->local_number_of_particles] = src_index;
+                        //tmp_local_index[this->local_number_of_particles] = iii;
                         std::copy(src.getParticlesState()+state_size*ii,
                                   src.getParticlesState()+state_size*(ii+1),
                                   tmp_local_state + state_size*this->local_number_of_particles);
diff --git a/cpp/particles/particles_output_sampling_hdf5.hpp b/cpp/particles/particles_output_sampling_hdf5.hpp
index 04174b15f46730de1f42d6cee15e0d36b255b58e..11a7dcc13ad2b292a9f8215b669bdcf4c97addb9 100644
--- a/cpp/particles/particles_output_sampling_hdf5.hpp
+++ b/cpp/particles/particles_output_sampling_hdf5.hpp
@@ -190,7 +190,7 @@ public:
             const int size_particle_rhs) final{
         assert(Parent::isInvolved());
 
-        TIMEZONE("particles_output_hdf5::write");
+        TIMEZONE("particles_output_sampling_hdf5::write");
 
         assert(particles_idx_offset < Parent::getTotalNbParticles() ||
                (particles_idx_offset == Parent::getTotalNbParticles() &&