diff --git a/cpp/full_code/ornstein_uhlenbeck_process.cpp b/cpp/full_code/ornstein_uhlenbeck_process.cpp
index d12c09bc14fcb50eacd1b1ac55ff5d31b36314f8..f8266fdb391d253026abaa2fe97bc263c49da0d2 100644
--- a/cpp/full_code/ornstein_uhlenbeck_process.cpp
+++ b/cpp/full_code/ornstein_uhlenbeck_process.cpp
@@ -48,7 +48,17 @@ ornstein_uhlenbeck_process<rnumber,be>::ornstein_uhlenbeck_process(
 
     assert(this->kk->kM2 >= this->ou_kmax_squ); //I can avoid one if statement later in the step method
 
+    gen.resize(omp_get_max_threads());
 
+    for(int thread=0;thread<omp_get_max_threads();thread++)
+    {
+            int current_seed = 
+                    this->ou_field->clayout->myrank*omp_get_max_threads() +
+                    thread;
+            gen[thread].seed(current_seed);
+    }
+
+    
     this->initialize_B();
 
 }
@@ -78,10 +88,11 @@ void ornstein_uhlenbeck_process<rnumber,be>::step(
         if (k2 <= this->ou_kmax_squ && k2 >= this->ou_kmin_squ)
         {
             double g = this->gamma(sqrt(k2));
+            int tr = omp_get_thread_num();
             double random[6] =
             {
-                this->d(this->gen),this->d(this->gen),this->d(this->gen),
-                this->d(this->gen),this->d(this->gen),this->d(this->gen)
+                this->d(this->gen[tr]),this->d(this->gen[tr]),this->d(this->gen[tr]),
+                this->d(this->gen[tr]),this->d(this->gen[tr]),this->d(this->gen[tr])
             };
             for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
                 this->ou_field->cval(cindex,cc,i) =
diff --git a/cpp/full_code/ornstein_uhlenbeck_process.hpp b/cpp/full_code/ornstein_uhlenbeck_process.hpp
index d20fddd986df652ea8d47a9707ca59bd0f5037f7..7710b46702ead3502ab93a7a6056be20d62ed542 100644
--- a/cpp/full_code/ornstein_uhlenbeck_process.hpp
+++ b/cpp/full_code/ornstein_uhlenbeck_process.hpp
@@ -31,8 +31,8 @@ class ornstein_uhlenbeck_process{
         std::vector<double> nom; //number of modes
         std::vector<double> aom; //absolute k value of modes
 
-        std::random_device rd{};
-        std::mt19937 gen{rd()}; //32-bit, maybe fix for double prec TODO
+        //std::random_device rd{};
+        std::vector<std::mt19937_64> gen; //32-bit, maybe fix for double prec TODO
         std::normal_distribution<double> d{0,1};
 
         ornstein_uhlenbeck_process(