diff --git a/bfps/DNS.py b/bfps/DNS.py
index 57dc879a9db8e2afcc8e3082d8cf0a9c58185efe..faf756c704144f2e771de78eb72fd0c4bb52f12d 100644
--- a/bfps/DNS.py
+++ b/bfps/DNS.py
@@ -865,6 +865,7 @@ class DNS(_code):
                     rseed = opt.particle_rand_seed)
             if not os.path.exists(self.get_particle_file_name()):
                 with h5py.File(self.get_particle_file_name(), 'w') as particle_file:
+                    particle_file.create_group('tracers0/position')
                     particle_file.create_group('tracers0/velocity')
                     particle_file.create_group('tracers0/acceleration')
         return None
diff --git a/bfps/cpp/full_code/NSVEparticles.cpp b/bfps/cpp/full_code/NSVEparticles.cpp
index 90b948b68756788b7a21007ff3c3ae602afaf0bc..2384e3f1abbc7fa979287d4f9e89b0d0b0718cbe 100644
--- a/bfps/cpp/full_code/NSVEparticles.cpp
+++ b/bfps/cpp/full_code/NSVEparticles.cpp
@@ -77,6 +77,14 @@ int NSVEparticles<rnumber>::do_stats()
     if (!(this->iteration % this->niter_part == 0))
         return EXIT_SUCCESS;
 
+    /// sample position
+    sample_particles_system_position(
+                                 this->ps,
+                                 (this->simname + "_particles.h5"), // filename
+                                 "tracers0",                        // hdf5 parent group
+                                 "position"                         // dataset basename TODO
+                                 );
+
     /// sample velocity
     sample_from_particles_system(*this->tmp_vec_field,              // field to save
                                  this->ps,
diff --git a/bfps/cpp/particles/particles_sampling.hpp b/bfps/cpp/particles/particles_sampling.hpp
index 3adc255341f3ca879d5cae1445124091f31b4394..c6b7e295b0465c15da9014d417da0623b80c597f 100644
--- a/bfps/cpp/particles/particles_sampling.hpp
+++ b/bfps/cpp/particles/particles_sampling.hpp
@@ -48,5 +48,37 @@ void sample_from_particles_system(const field<rnumber, be, fc>& in_field, // a p
                      ps->get_step_idx());
 }
 
-#endif
+template <class partsize_t, class particles_rnumber>
+void sample_particles_system_position(
+        std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>& ps, // a pointer to an particles_system<double>
+                                  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());
+
+    // Stop here if already exists
+    if(particles_output_sampling_hdf5<partsize_t, particles_rnumber, 3, 3>::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[3*nb_particles]);
+    std::copy(ps->getParticlesPositions(), ps->getParticlesPositions() + 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->getParticlesPositions(),
+                     &sample_rhs,
+                     ps->getParticlesIndexes(),
+                     ps->getLocalNbParticles(),
+                     ps->get_step_idx());
+}
+
+#endif//PARTICLES_SAMPLING_HPP