diff --git a/TurTLE/DNS.py b/TurTLE/DNS.py
index c0ac8328553e9c8a5e6785fef24a2150de288cb4..686d74a0f577dbf74db61baf8bf031408d80a19d 100644
--- a/TurTLE/DNS.py
+++ b/TurTLE/DNS.py
@@ -140,6 +140,7 @@ class DNS(_code):
         self.NSVEp_extra_parameters['niter_part_fine_period'] = int(10)
         self.NSVEp_extra_parameters['niter_part_fine_duration'] = int(0)
         self.NSVEp_extra_parameters['nparticles'] = int(10)
+        self.NSVEp_extra_parameters['cpp_random_particles'] = int(0)
         self.NSVEp_extra_parameters['tracers0_integration_steps'] = int(4)
         self.NSVEp_extra_parameters['tracers0_neighbours'] = int(1)
         self.NSVEp_extra_parameters['tracers0_smoothness'] = int(1)
@@ -914,7 +915,7 @@ class DNS(_code):
                             if self.dns_type == 'NSVE_Stokes_particles':
                                 dset[cc*batch_size:(cc+1)*batch_size, 3:] = self.parameters['initial_particle_vel']*get_random_versors(batch_size)
                             else:
-                                dset[cc*batch_size:(cc+1)*batch_size, 3:] = get_random_versors(batch_size)                             
+                                dset[cc*batch_size:(cc+1)*batch_size, 3:] = get_random_versors(batch_size)
                         nn -= batch_size
                     else:
                         dset[cc*batch_size:cc*batch_size+nn, :3] = get_random_phases(nn)
@@ -1057,10 +1058,11 @@ class DNS(_code):
     def generate_particle_data(
             self,
             opt = None):
-        if self.parameters['nparticles'] > 0:
-            self.generate_tracer_state(
-                    species = 0,
-                    rseed = opt.particle_rand_seed)
+        if (self.parameters['nparticles'] > 0):
+            if (self.parameters['cpp_random_particles'] == 0):
+                self.generate_tracer_state(
+                        species = 0,
+                        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')
@@ -1136,7 +1138,7 @@ class DNS(_code):
                     data = self.generate_vector_field(
                        write_to_file = False,
                        spectra_slope = 2.0,
-                       amplitude = 0.05)                
+                       amplitude = 0.05)
                 f['vorticity/complex/{0}'.format(0)] = data
             f.close()
         if self.dns_type == 'kraichnan_field':
diff --git a/cpp/full_code/NSVEparticles.cpp b/cpp/full_code/NSVEparticles.cpp
index 81e946771bccca1046e60ece2c568bf131ce0db7..ecae213e3f43f90fbba91cbffe0e3e97bdb02f98 100644
--- a/cpp/full_code/NSVEparticles.cpp
+++ b/cpp/full_code/NSVEparticles.cpp
@@ -28,6 +28,8 @@
 #include <cmath>
 #include "NSVEparticles.hpp"
 #include "scope_timer.hpp"
+#include "particles/interpolation/particle_set.hpp"
+#include "particles/particles_input_random.hpp"
 
 template <typename rnumber>
 int NSVEparticles<rnumber>::initialize(void)
@@ -41,6 +43,63 @@ int NSVEparticles<rnumber>::initialize(void)
             this->fs->cvelocity->rlayout->comm,
             this->fs->cvelocity->fftw_plan_rigor);
 
+    // initialize particle output object
+    // this may be needed before the particle system,
+    // if cpp_random_particles is true
+    this->particles_output_writer_mpi = new particles_output_hdf5<
+        long long int, double, 3>(
+                MPI_COMM_WORLD,
+                "tracers0",
+                nparticles,
+                tracers0_integration_steps);
+
+    // if this is the first iteration, and we are supposed to generate our own random particles,
+    // we do that here, and then we write them to file for future reference.
+    if ((this->cpp_random_particles > 0) &&
+         this->fs->iteration == 0)
+    {
+        // temporary particle set
+        // interpolation setup is irrelevant for output
+        particle_set<3, 1, 1> pset(
+                this->fs->cvelocity->rlayout,
+                this->fs->kk->dkx,
+                this->fs->kk->dky,
+                this->fs->kk->dkz);
+        particles_input_random<long long int, double, 3> pinput(
+                this->comm,
+                this->nparticles,
+                this->cpp_random_particles, // seed
+                pset.getSpatialLowLimitZ(),
+                pset.getSpatialUpLimitZ());
+        // initialize particle set
+        pset.init(pinput);
+        // allocate dummy rhs data
+        // we must fill these arrays with 0
+        // otherwise we get floating point exceptions because data is being written/read
+        std::vector<std::unique_ptr<double[]>> rhs_data;
+        rhs_data.resize(tracers0_integration_steps);
+        for (int counter = 0; counter < tracers0_integration_steps; counter++)
+        {
+            rhs_data[counter].reset(new double[pset.getLocalNumberOfParticles()*3]);
+            std::fill_n(rhs_data[counter].get(), pset.getLocalNumberOfParticles()*3, 0);
+        }
+        // create dummy file_layout object
+        std::vector<hsize_t> file_layout;
+        file_layout.resize(1);
+        file_layout[0] = this->nparticles;
+        // write data
+        this->particles_output_writer_mpi->open_file(this->fs->get_current_fname());
+        this->particles_output_writer_mpi->update_particle_species_name("tracers0");
+        this->particles_output_writer_mpi->setParticleFileLayout(file_layout);
+        this->particles_output_writer_mpi->template save<3>(
+                pset.getParticleState(),
+                rhs_data.data(),
+                pset.getParticleIndex(),
+                pset.getLocalNumberOfParticles(),
+                0);
+        this->particles_output_writer_mpi->close_file();
+    }
+
     DEBUG_MSG_WAIT(MPI_COMM_WORLD, "about to call particles_system_builder\n");
     this->ps = particles_system_builder(
                 this->fs->cvelocity,              // (field object)
@@ -55,13 +114,8 @@ int NSVEparticles<rnumber>::initialize(void)
                 this->comm,
                 this->fs->iteration+1);
     DEBUG_MSG_WAIT(MPI_COMM_WORLD, "after call to particles_system_builder\n");
-    this->particles_output_writer_mpi = new particles_output_hdf5<
-        long long int, double, 3>(
-                MPI_COMM_WORLD,
-                "tracers0",
-                nparticles,
-                tracers0_integration_steps);
-    this->particles_output_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
+
+    // initialize sample object
     this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
         long long int, double, 3>(
                 MPI_COMM_WORLD,
@@ -69,6 +123,9 @@ int NSVEparticles<rnumber>::initialize(void)
                 (this->simname + "_particles.h5"),
                 "tracers0",
                 "position/0");
+
+    // set particle file layout for output objects
+    this->particles_output_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
     this->particles_sample_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
     return EXIT_SUCCESS;
 }
@@ -198,6 +255,7 @@ int NSVEparticles<rnumber>::read_parameters(void)
     this->niter_part_fine_period = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part_fine_period");
     this->niter_part_fine_duration = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part_fine_duration");
     this->nparticles = hdf5_tools::read_value<int>(parameter_file, "parameters/nparticles");
+    this->cpp_random_particles = hdf5_tools::read_value<int>(parameter_file, "parameters/cpp_random_particles");
     this->tracers0_integration_steps = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_integration_steps");
     this->tracers0_neighbours = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_neighbours");
     this->tracers0_smoothness = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_smoothness");
diff --git a/cpp/full_code/NSVEparticles.hpp b/cpp/full_code/NSVEparticles.hpp
index 8b70ead9b084aa4f7693c243b2f04e5b06c2d572..6dc3ca8e75011f6b942be1e4ad0b529a44b585ae 100644
--- a/cpp/full_code/NSVEparticles.hpp
+++ b/cpp/full_code/NSVEparticles.hpp
@@ -57,6 +57,7 @@ class NSVEparticles: public NSVE<rnumber>
         int tracers0_integration_steps;
         int tracers0_neighbours;
         int tracers0_smoothness;
+        int cpp_random_particles;
 
         /* other stuff */
         std::unique_ptr<abstract_particles_system<long long int, double>> ps;
diff --git a/cpp/particles/particles_system_builder.hpp b/cpp/particles/particles_system_builder.hpp
index 038c5774c1fc8fe93d2d2b4bb9476dc3ce6b386c..a44b2e63a329df664903fc6f1fa96f281cf85b36 100644
--- a/cpp/particles/particles_system_builder.hpp
+++ b/cpp/particles/particles_system_builder.hpp
@@ -40,7 +40,7 @@
 #include "kspace.hpp"
 
 const int MAX_INTERPOLATION_NEIGHBOURS=5;
-const int MAX_INTERPOLATION_SMOOTHNESS=2;
+const int MAX_INTERPOLATION_SMOOTHNESS=3;
 
 //////////////////////////////////////////////////////////////////////////////
 ///
diff --git a/tests/DNS/test_scaling.py b/tests/DNS/test_scaling.py
index 6a287492dd70b7e53b8d1a9dd04be5796fd01e69..ce16072e0356ca3e7e4a1baeccf1d48ac76b2ef3 100644
--- a/tests/DNS/test_scaling.py
+++ b/tests/DNS/test_scaling.py
@@ -57,6 +57,7 @@ def get_DNS_parameters(
         DNS_parameters += [
                 '--precision', 'double']
     # check that source sim exists
+    print('looking for ', os.path.join(src_dirname, src_simname + '.h5'))
     assert(os.path.exists(os.path.join(src_dirname, src_simname + '.h5')))
     # check that source checkpoint exists
     dns_src = TurTLE.DNS(simname = src_simname, work_dir = src_dirname)
@@ -87,7 +88,7 @@ def get_DNS_parameters(
         DNS_parameters += [
                 '--tracers0_neighbours', '{0}'.format(nneighbours),
                 '--tracers0_smoothness', '{0}'.format(smoothness),
-                '--particle-rand-seed', '2']
+                '--cpp_random_particles', '2']
     if no_submit:
         DNS_parameters += ['--no-submit']
     DNS_parameters += ['--environment', environment,