diff --git a/bfps/cpp/particles/particles_output_sampling_hdf5.hpp b/bfps/cpp/particles/particles_output_sampling_hdf5.hpp
index 05a0aee6a46c996a9549cd7468d2208758a3a074..09b09a939a77d44cf461aacb988368fc5add7052 100644
--- a/bfps/cpp/particles/particles_output_sampling_hdf5.hpp
+++ b/bfps/cpp/particles/particles_output_sampling_hdf5.hpp
@@ -24,6 +24,36 @@ class particles_output_sampling_hdf5 : public abstract_particles_output<partsize
     const bool use_collective_io;
 
 public:
+    static bool DatasetExistsCol(MPI_Comm in_mpi_com,
+                                  const std::string& in_filename,
+                                  const std::string& in_groupname,
+                                 const std::string& in_dataset_name){
+        int my_rank;
+        AssertMpi(MPI_Comm_rank(in_mpi_com, &my_rank));
+
+        int dataset_exists = -1;
+
+        if(my_rank == 0){
+            // Parallel HDF5 write
+            hid_t file_id = H5Fopen(
+                    in_filename.c_str(),
+                    H5F_ACC_RDWR | H5F_ACC_DEBUG,
+                    H5P_DEFAULT);
+            assert(file_id >= 0);
+
+            dataset_exists = H5Lexists(
+                    file_id,
+                    (in_groupname + "/" + in_dataset_name).c_str(),
+                    H5P_DEFAULT);
+
+            int retTest = H5Fclose(file_id);
+            assert(retTest >= 0);
+        }
+
+        AssertMpi(MPI_Bcast( &dataset_exists, 1, MPI_INT, 0, in_mpi_com ));
+        return dataset_exists;
+    }
+
     particles_output_sampling_hdf5(MPI_Comm in_mpi_com,
                           const partsize_t inTotalNbParticles,
                                    const std::string& in_filename,
diff --git a/bfps/cpp/particles/particles_sampling.hpp b/bfps/cpp/particles/particles_sampling.hpp
index 8ec52c9a81f753d368db4d678030adee3bb7a992..0833b45f2b65d88966fc240584fcd4a45d68ef5a 100644
--- a/bfps/cpp/particles/particles_sampling.hpp
+++ b/bfps/cpp/particles/particles_sampling.hpp
@@ -17,13 +17,23 @@ void sample_from_particles_system(const field<rnumber, be, fc>& in_field, // a p
                                   const std::string& filename,
                                   const std::string& parent_groupname,
                                   const std::string& fname){
+    const std::string datasetname = fname + std::string("/") + std::to_string(ps->get_step_idx());
     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)){
+        return;
+    }
+
     const partsize_t nb_particles = ps->getLocalNbParticles();
     std::unique_ptr<particles_rnumber[]> sample_rhs(new particles_rnumber[size_particle_rhs*nb_particles]);
 
     ps->sample_compute_field(in_field, sample_rhs.get());
 
-    const std::string datasetname = fname + std::string("/") + std::to_string(ps->get_step_idx());
+
 
     particles_output_sampling_hdf5<partsize_t, particles_rnumber, 3, size_particle_rhs> outputclass(MPI_COMM_WORLD,
                                                                                                     ps->getGlobalNbParticles(),