diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index 6f345644a3a5933ecb93d918b09601cf864675b8..0b26266ef9ed99c9aef3490b85c7ceaf2e8d2096 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -240,7 +240,10 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
             smoothness = 1,
             name = 'particle_field'):
         self.particle_variables += 'interpolator<{0}, {1}> *vel_{2}, *acc_{2};\n'.format(self.C_dtype, neighbours, name)
+        self.parameters[name + '_type'] = interp_type
+        self.parameters[name + '_neighbours'] = neighbours
         if interp_type == 'spline':
+            self.parameters[name + '_smoothness'] = smoothness
             beta_name = 'beta_n{0}_m{1}'.format(neighbours, smoothness)
         elif interp_type == 'Lagrange':
             beta_name = 'beta_Lagrange_n{0}'.format(neighbours)
@@ -263,9 +266,6 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
         return None
     def add_particles(
             self,
-            neighbours = 1,
-            smoothness = 1,
-            interp_type = 'spline',
             integration_method = 'AdamsBashforth',
             integration_steps = 2,
             kcut = 'fs->kM',
@@ -275,12 +275,11 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
             integration_steps = 4
         elif integration_method == 'Heun':
             integration_steps = 2
-        self.parameters['integration_method{0}'.format(self.particle_species)] = integration_method
-        self.parameters['interp_type{0}'.format(self.particle_species)] = interp_type
-        self.parameters['neighbours{0}'.format(self.particle_species)] = neighbours
-        self.parameters['smoothness{0}'.format(self.particle_species)] = smoothness
-        self.parameters['kcut{0}'.format(self.particle_species)] = kcut
-        self.parameters['integration_steps{0}'.format(self.particle_species)] = integration_steps
+        neighbours = self.parameters[fields_name + '_neighbours']
+        self.parameters['tracers{0}_field'.format(self.particle_species)] = fields_name
+        self.parameters['tracers{0}_integration_method'.format(self.particle_species)] = integration_method
+        self.parameters['tracers{0}_kcut'.format(self.particle_species)] = kcut
+        self.parameters['tracers{0}_integration_steps'.format(self.particle_species)] = integration_steps
         self.file_datasets_grow += """
                         //begincpp
                         temp_string = (std::string("/particles/") +
@@ -341,10 +340,6 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
                           }}
                           //endcpp
                           """.format(self.particle_species, update_field, compute_acc)
-        if interp_type == 'spline':
-            beta_name = 'beta_n{0}_m{1}'.format(neighbours, smoothness)
-        elif interp_type == 'Lagrange':
-            beta_name = 'beta_Lagrange_n{0}'.format(neighbours)
         self.particle_start += 'sprintf(fname, "tracers{0}");\n'.format(self.particle_species)
         self.particle_end += ('ps{0}->write(stat_file);\n' +
                               'delete ps{0};\n').format(self.particle_species)
@@ -361,7 +356,7 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
         self.particle_start += ('ps{0} = new particles<VELOCITY_TRACER, {1}, {2},{3}>(\n' +
                                     'fname, fs, vel_{4},\n' +
                                     'nparticles,\n' +
-                                    'niter_part, integration_steps{0});\n').format(
+                                    'niter_part, tracers{0}_integration_steps);\n').format(
                                             self.particle_species, self.C_dtype, multistep, neighbours, fields_name)
         self.particle_start += ('ps{0}->dt = dt;\n' +
                                 'ps{0}->iteration = iteration;\n' +
diff --git a/bfps/cpp/particles.cpp b/bfps/cpp/particles.cpp
index 7109f409c60d73798a717a9b74bfb7ce532fbea1..8a4848a78d93a6bd5d3c1d8c54a5e17f14f0f79e 100644
--- a/bfps/cpp/particles.cpp
+++ b/bfps/cpp/particles.cpp
@@ -289,7 +289,7 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::synchroniz
             MPI_DOUBLE,
             MPI_SUM,
             this->comm);
-    if (this->integration_steps >= 1)
+    if (this->integration_steps >= 1 && multistep)
     {
         for (int i=0; i<this->integration_steps; i++)
         {
@@ -586,7 +586,7 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::read(hid_t
         H5Sclose(mspace);
         H5Sclose(rspace);
         H5Dclose(dset);
-        if (this->iteration > 0)
+        if (this->iteration > 0 && multistep)
         {
             temp_string = (std::string("/particles/") +
                            std::string(this->name) +
@@ -617,15 +617,13 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::read(hid_t
             MPI_DOUBLE,
             0,
             this->comm);
-    for (int i = 0; i<this->integration_steps; i++)
-    {
+    if (multistep) for (int i = 0; i<this->integration_steps; i++)
         MPI_Bcast(
                 this->rhs[i],
                 this->array_size,
                 MPI_DOUBLE,
                 0,
                 this->comm);
-    }
     // initial assignment of particles
     for (int p=0; p<this->nparticles; p++)
     {
@@ -659,7 +657,7 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::write(hid_
         H5Sclose(mspace);
         H5Sclose(wspace);
         H5Dclose(dset);
-        if (write_rhs)
+        if (write_rhs && multistep)
         {
             temp_string = (std::string("/particles/") +
                            std::string(this->name) +
diff --git a/bfps/fluid_base.py b/bfps/fluid_base.py
index 8d06010328c1fba5368154ef0a01d8913d7e7c93..999e7e372e3385f53a491384972e07e4f78ca5fb 100644
--- a/bfps/fluid_base.py
+++ b/bfps/fluid_base.py
@@ -129,9 +129,10 @@ class fluid_particle_base(bfps.code):
                                  'grow_single_dataset(dset, niter_todo/niter_part);\n' +
                                  'H5Dclose(dset);\n')
         self.definitions += ('full_name = (std::string(name) + std::string("/rhs"));\n' +
+                             'if (H5Lexists(g_id, full_name.c_str(), H5P_DEFAULT))\n{\n' +
                              'dset = H5Dopen(g_id, full_name.c_str(), H5P_DEFAULT);\n' +
                              'grow_single_dataset(dset, 1);\n' +
-                             'H5Dclose(dset);\n' +
+                             'H5Dclose(dset);\n}\n' +
                              'return 0;\n}\n')
         self.definitions += ('int grow_file_datasets()\n{\n' +
                              'int file_problems = 0;\n' +
@@ -429,24 +430,25 @@ class fluid_particle_base(bfps.code):
                                      dtype = np.int64,
                                      compression = 'gzip')
             for s in range(self.particle_species):
-                time_chunk = 2**20 // (8*3*
-                                       self.parameters['nparticles']*
-                                       self.parameters['integration_steps{0}'.format(s)])
-                time_chunk = max(time_chunk, 1)
-                ofile.create_dataset('particles/tracers{0}/rhs'.format(s),
-                                     (1,
-                                      self.parameters['integration_steps{0}'.format(s)],
-                                      self.parameters['nparticles'],
-                                      3),
-                                     maxshape = (None,
-                                                 self.parameters['integration_steps{0}'.format(s)],
-                                                 self.parameters['nparticles'],
-                                                 3),
-                                     chunks =  (time_chunk,
-                                                self.parameters['integration_steps{0}'.format(s)],
-                                                self.parameters['nparticles'],
-                                                3),
-                                     dtype = np.float64)
+                if self.parameters['tracers{0}_integration_method'.format(s)] == 'AdamsBashforth':
+                    time_chunk = 2**20 // (8*3*
+                                           self.parameters['nparticles']*
+                                           self.parameters['tracers{0}_integration_steps'.format(s)])
+                    time_chunk = max(time_chunk, 1)
+                    ofile.create_dataset('particles/tracers{0}/rhs'.format(s),
+                                         (1,
+                                          self.parameters['tracers{0}_integration_steps'.format(s)],
+                                          self.parameters['nparticles'],
+                                          3),
+                                         maxshape = (None,
+                                                     self.parameters['tracers{0}_integration_steps'.format(s)],
+                                                     self.parameters['nparticles'],
+                                                     3),
+                                         chunks =  (time_chunk,
+                                                    self.parameters['tracers{0}_integration_steps'.format(s)],
+                                                    self.parameters['nparticles'],
+                                                    3),
+                                         dtype = np.float64)
                 time_chunk = 2**20 // (8*3*self.parameters['nparticles'])
                 time_chunk = max(time_chunk, 1)
                 ofile.create_dataset(
diff --git a/tests/test_base.py b/tests/test_base.py
index 6ec5aa63a2109585940c941b0667d975acb8b180..aafd41e6ce8c7ac129e34758e958fd9366c8b098 100755
--- a/tests/test_base.py
+++ b/tests/test_base.py
@@ -97,7 +97,7 @@ def launch(
         c.add_particle_fields(kcut = 'fs->kM/2', name = 'filtered', neighbours = opt.neighbours)
         c.add_particles(
                 kcut = 'fs->kM/2',
-                integration_steps = 1, neighbours = opt.neighbours, smoothness = opt.smoothness,
+                integration_steps = 1,
                 fields_name = 'filtered')
         #for integr_steps in range(1, 7):
         #    c.add_particles(
@@ -105,16 +105,11 @@ def launch(
         #            neighbours = opt.neighbours,
         #            smoothness = opt.smoothness,
         #            fields_name = 'regular')
-        for info in [(2, 1, 'spline', 'Heun'),
-                     (4, 1, 'spline', 'cRK4'),
-                     (2, 1, 'Lagrange', 'Heun'),
-                     (4, 1, 'Lagrange', 'cRK4')]:
+        for info in [(2, 'Heun'),
+                     (4, 'cRK4')]:
             c.add_particles(
                     integration_steps = info[0],
-                    neighbours = opt.neighbours,
-                    smoothness = info[1],
-                    interp_type = info[2],
-                    integration_method = info[3],
+                    integration_method = info[1],
                     fields_name = 'regular')
     c.fill_up_fluid_code()
     c.finalize_code()