diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index 8eea54a1d6d71db562badf85280f34f20099906c..a5a0abd2350241b2604f22725c30aceb3b44b827 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -536,10 +536,11 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
                 self.particle_loop += 'ps{0}->synchronize();\n'.format(self.particle_species)
             elif particle_class == 'rFFTW_particles':
                 self.particle_loop += 'ps{0}->step();\n'.format(self.particle_species)
-        self.particle_loop += (('if (ps{0}->iteration % niter_part == 0)\n' +
-                                '{{\n' +
-                                'ps{0}->write(stat_file, false);\n').format(self.particle_species) +
-                               output_vel_acc + '}\n')
+        self.particle_stat_src += (
+                ('if (ps{0}->iteration % niter_part == 0)\n' +
+                 '{{\n' +
+                 'ps{0}->write(stat_file, false);\n').format(self.particle_species) +
+                output_vel_acc + '}\n')
         self.particle_species += 1
         return None
     def get_data_file(self):
diff --git a/bfps/cpp/rFFTW_particles.cpp b/bfps/cpp/rFFTW_particles.cpp
index 8a8f2ffdf7173a3e9ec2b75d84cfa4b2112fddc9..151d057a7958260c50302f08c28b708eb58ab822 100644
--- a/bfps/cpp/rFFTW_particles.cpp
+++ b/bfps/cpp/rFFTW_particles.cpp
@@ -112,7 +112,7 @@ template <int particle_type, class rnumber, int interp_neighbours>
 void rFFTW_particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth(int nsteps)
 {
     ptrdiff_t ii;
-    this->get_rhs(0, this->state, this->rhs[0]);
+    this->get_rhs(1, this->state, this->rhs[0]);
     switch(nsteps)
     {
         case 1:
diff --git a/bfps/fluid_base.py b/bfps/fluid_base.py
index 62dfcfd264e2b6f0987ac0cfe547873e840fb686..ba44f306bce344d73ef39c51618e51b48bf98a6c 100644
--- a/bfps/fluid_base.py
+++ b/bfps/fluid_base.py
@@ -90,13 +90,14 @@ class fluid_particle_base(bfps.code):
         self.fluid_loop = ''
         self.fluid_end  = ''
         self.fluid_output = ''
+        self.stat_src = ''
         self.particle_includes = ''
         self.particle_variables = ''
         self.particle_definitions = ''
         self.particle_start = ''
         self.particle_loop = ''
         self.particle_end  = ''
-        self.stat_src = ''
+        self.particle_stat_src = ''
         self.file_datasets_grow   = ''
         return None
     def finalize_code(self):
@@ -144,6 +145,7 @@ class fluid_particle_base(bfps.code):
                              'return file_problems;\n'
                              '}\n')
         self.definitions += 'void do_stats()\n{\n' + self.stat_src + '}\n'
+        self.definitions += 'void do_particle_stats()\n{\n' + self.particle_stat_src + '}\n'
         # take care of wisdom
         if self.use_fftw_wisdom:
             if self.dtype == np.float32:
@@ -191,6 +193,7 @@ class fluid_particle_base(bfps.code):
                                return EXIT_SUCCESS;
                            }
                            do_stats();
+                           do_particle_stats();
                            //endcpp
                            """
         output_time_difference = ('time1 = clock();\n' +
@@ -202,16 +205,21 @@ class fluid_particle_base(bfps.code):
                                       '<< iteration << " took " ' +
                                       '<< time_difference/nprocs << " seconds" << std::endl;\n' +
                                   'time0 = time1;\n')
-        self.main       += 'for (int max_iter = iteration+niter_todo; iteration < max_iter; iteration++)\n{\n'
+        self.main       += 'for (int max_iter = iteration+niter_todo; iteration < max_iter; iteration++)\n'
+        self.main       += '{\n'
         self.main       += output_time_difference
-        self.main       += self.fluid_loop
         if self.particle_species > 0:
             self.main   += self.particle_loop
-        self.main       += 'if (iteration % niter_stat == 0) do_stats();\n}\n'
+        self.main       += self.fluid_loop
+        self.main       += 'if (iteration % niter_stat == 0) do_stats();\n'
+        if self.particle_species > 0:
+            self.main       += 'if (iteration % niter_part == 0) do_particle_stats();\n'
+        self.main       += '}\n'
         self.main       += output_time_difference
+        self.main       += 'do_stats();\n'
+        self.main       += 'do_particle_stats();\n'
         if self.particle_species > 0:
             self.main   += self.particle_end
-        self.main       += 'do_stats();\n'
         self.main       += self.fluid_end
         return None
     def read_rfield(