Require Seed

doc/tutorials/01-lennard_jones/scripts/lj_tutorial.py

@@ -66,6 +66,8 @@
 # System setup
 #############################################################
 system              = espressomd.System(box_l=[1.0, 1.0, 1.0])
+system.seed         = system.cell_system.get_state()['n_nodes'] * [1234]
+np.random.seed(system.seed)
 
 if not os.path.exists('data') :
     os.mkdir('data')

doc/tutorials/01-lennard_jones/scripts/msd.py

@@ -73,6 +73,8 @@
 # System setup
 #############################################################
 system              = espressomd.System(box_l=[1.0, 1.0, 1.0])
+system.seed         = system.cell_system.get_state()['n_nodes'] * [1234]
+np.random.seed(system.seed)
 
 if not os.path.exists('data') :
     os.mkdir('data')

doc/tutorials/01-lennard_jones/scripts/two-component-visualization.py

@@ -74,6 +74,8 @@
 # System setup
 #############################################################
 system              = espressomd.System(box_l=[1.0, 1.0, 1.0])
+system.seed         = system.cell_system.get_state()['n_nodes'] * [1234]
+np.random.seed(system.seed)
 
 if not os.path.exists('data') :
     os.mkdir('data')

doc/tutorials/01-lennard_jones/scripts/two-component.py

@@ -70,6 +70,8 @@
 # System setup
 #############################################################
 system              = espressomd.System(box_l=[1.0, 1.0, 1.0])
+system.seed         = system.cell_system.get_state()['n_nodes'] * [1234]
+np.random.seed(system.seed)
 
 if not os.path.exists('data') :
     os.mkdir('data')

doc/tutorials/02-charged_system/scripts/nacl.py

@@ -52,6 +52,7 @@
 # Setup System
 box_l = (n_part / density)**(1. / 3.)
 system = espressomd.System(box_l=[box_l] * 3)
+system.seed  = system.cell_system.get_state()['n_nodes'] * [1234]
 system.periodicity = [1, 1, 1]
 system.time_step = time_step
 system.cell_system.skin = 0.3

doc/tutorials/02-charged_system/scripts/nacl_units.py

@@ -24,6 +24,8 @@
 assert_features(["ELECTROSTATICS", "MASS", "LENNARD_JONES"])
 
 system = espressomd.System(box_l=[1.0, 1.0, 1.0])
+system.seed  = system.cell_system.get_state()['n_nodes'] * [1234]
+numpy.random.seed(system.seed)
 
 print("\n--->Setup system")
 
@@ -172,3 +174,5 @@ def combination_rule_sigma(rule, sig1, sig2):
 rdf_fp.close()
 print("\n--->Written rdf.data")
 print("\n--->Done")
+
+

doc/tutorials/02-charged_system/scripts/nacl_units_confined.py

@@ -25,6 +25,8 @@
 assert_features(["ELECTROSTATICS", "CONSTRAINTS", "MASS", "LENNARD_JONES"])
 
 system = espressomd.System(box_l=[1.0, 1.0, 1.0])
+system.seed  = system.cell_system.get_state()['n_nodes'] * [1234]
+numpy.random.seed(system.seed)
 
 print("\n--->Setup system")
 

doc/tutorials/02-charged_system/scripts/nacl_units_confined_vis.py

@@ -28,6 +28,8 @@
 assert_features(["ELECTROSTATICS", "CONSTRAINTS", "MASS", "LENNARD_JONES"])
 
 system = espressomd.System(box_l=[1.0, 1.0, 1.0])
+system.seed  = system.cell_system.get_state()['n_nodes'] * [1234]
+numpy.random.seed(system.seed)
 
 print("\n--->Setup system")
 

doc/tutorials/02-charged_system/scripts/nacl_units_vis.py

@@ -27,6 +27,9 @@
 assert_features(["ELECTROSTATICS", "MASS", "LENNARD_JONES"])
 
 system = espressomd.System(box_l=[1.0, 1.0, 1.0])
+system.seed  = system.cell_system.get_state()['n_nodes'] * [1234]
+numpy.random.seed(system.seed)
+
 visualizer = openGLLive(system, drag_force=5 * 298,
                         background_color=[1, 1, 1], light_pos=[30, 30, 30])
 

doc/tutorials/04-lattice_boltzmann/scripts/polymer_diffusion.py

@@ -12,6 +12,7 @@
 
 # System setup
 system = espressomd.System(box_l=[32, 32, 32])
+system.seed = system.cell_system.get_state()['n_nodes'] * [1234]
 system.cell_system.skin = 0.4
 
 try:

doc/tutorials/04-lattice_boltzmann/scripts/single_particle_diffusion.py

@@ -19,6 +19,7 @@
 
 # System setup
 system = espressomd.System(box_l=[box_l, box_l, box_l])
+system.seed  = system.cell_system.get_state()['n_nodes'] * [1234]
 system.time_step = time_step
 system.cell_system.skin = 0.4
 

doc/tutorials/05-raspberry_electrophoresis/scripts/simulate_raspberry_electrophoresis.py

@@ -37,6 +37,8 @@
 # System setup
 #############################################################
 system = espressomd.System(box_l=[1.0, 1.0, 1.0])
+system.seed  = system.cell_system.get_state()['n_nodes'] * [1234]
+np.random.seed(seed=system.seed)
 system.box_l = [box_l, box_l, box_l]
 system.cell_system.skin = skin
 system.periodicity = [1, 1, 1]

doc/tutorials/06-active_matter/EXERCISES/enhanced_diffusion.py

@@ -62,7 +62,7 @@
 # Why can we get away with such a small box?
 # Could it be even smaller?
 system = espressomd.System(box_l=[10.0, 10.0, 10.0])
-
+system.seed  = system.cell_system.get_state()['n_nodes'] * [1234]
 
 system.cell_system.skin = 0.3
 system.time_step = tstep

doc/tutorials/06-active_matter/EXERCISES/flow_field.py

@@ -63,6 +63,7 @@
 dt = 0.01
 
 system = espressomd.System(box_l=[length, length, length])
+system.seed  = system.cell_system.get_state()['n_nodes'] * [1234]
 system.cell_system.skin = 0.3
 system.time_step = dt
 system.min_global_cut = 1.0

doc/tutorials/06-active_matter/EXERCISES/rectification_geometry.py

@@ -53,6 +53,7 @@
 # Setup the MD parameters
 
 system = espressomd.System(box_l=[length, diameter + 4, diameter + 4])
+system.seed  = system.cell_system.get_state()['n_nodes'] * [1234]
 system.cell_system.skin = 0.1
 system.time_step = dt
 system.min_global_cut = 0.5

doc/tutorials/06-active_matter/EXERCISES/rectification_simulation.py

@@ -88,6 +88,7 @@ def a2quat(phi, theta):
 # Setup the MD parameters
 
 system = System(box_l=[length, diameter + 4, diameter + 4])
+system.seed  = system.cell_system.get_state()['n_nodes'] * [1234]
 system.cell_system.skin = 0.1
 system.time_step = dt
 system.min_global_cut = 0.5

doc/tutorials/08-visualization/scripts/plotting.py

@@ -25,6 +25,7 @@
 # Integration parameters
 #############################################################
 system = espressomd.System(box_l=[box_l, box_l, box_l])
+system.seed  = system.cell_system.get_state()['n_nodes'] * [1234]
 system.time_step = 0.01
 system.cell_system.skin = 0.4
 system.thermostat.set_langevin(kT=1.0, gamma=1.0)

doc/tutorials/08-visualization/scripts/simulation.py

@@ -25,6 +25,7 @@
 # Integration parameters
 #############################################################
 system = espressomd.System(box_l=[box_l, box_l, box_l])
+system.seed  = system.cell_system.get_state()['n_nodes'] * [1234]
 system.time_step = 0.01
 system.cell_system.skin = 0.4
 system.thermostat.set_langevin(kT=1.0, gamma=1.0)

doc/tutorials/08-visualization/scripts/visualization.py

@@ -25,6 +25,7 @@
 # Integration parameters
 #############################################################
 system = espressomd.System(box_l=[box_l, box_l, box_l])
+system.seed  = system.cell_system.get_state()['n_nodes'] * [1234]
 system.time_step = 0.01
 system.cell_system.skin = 0.4
 system.thermostat.set_langevin(kT=1.0, gamma=1.0)

src/core/cuda_common_cuda.cu

@@ -25,13 +25,13 @@
 #include "cuda_utils.hpp"
 #include "errorhandling.hpp"
 #include "interaction_data.hpp"
-#include "random.hpp"
+
+#include <random>
 
 #if defined(OMPI_MPI_H) || defined(_MPI_H)
 #error CU-file includes mpi.h! This should not happen!
 #endif
 
-static int max_ran = 1000000;
 static CUDA_global_part_vars global_part_vars_host = {0, 0, 0};
 static __device__ __constant__ CUDA_global_part_vars global_part_vars_device;
 
@@ -197,7 +197,7 @@ void gpu_change_number_of_part_to_comm() {
   if (global_part_vars_host.number_of_particles != n_part &&
       global_part_vars_host.communication_enabled == 1 && this_node == 0) {
 
-    global_part_vars_host.seed = (unsigned int)i_random(max_ran);
+    global_part_vars_host.seed = (unsigned int)std::random_device{}();
     global_part_vars_host.number_of_particles = n_part;
 
     cuda_safe_mem(cudaMemcpyToSymbol(global_part_vars_device,

src/core/dpd.cpp

@@ -135,7 +135,12 @@ Vector3d dpd_pair_force(const Particle *p1, const Particle *p2, IA_parameters *i
       vel12_dot_d12 += (p1->m.v[j] - p2->m.v[j]) * d[j];
     auto const friction = ia_params->dpd_pref1 * omega2 * vel12_dot_d12 * time_step;
     // random force prefactor
-    auto const noise = ia_params->dpd_pref2 * omega * (d_random() - 0.5);
+    double noise;
+    if (ia_params->dpd_pref2 > 0.0) {
+      noise = ia_params->dpd_pref2 * omega * (d_random() - 0.5);
+    } else {
+      noise = 0.0;
+    }
     for (int j = 0; j < 3; j++) {
        f[j] += (noise - friction) * d[j];
     }
@@ -149,7 +154,11 @@ Vector3d dpd_pair_force(const Particle *p1, const Particle *p2, IA_parameters *i
             noise_vec[3];
     for (int i = 0; i < 3; i++) {
       // noise vector
-      noise_vec[i] = d_random() - 0.5;
+      if (ia_params->dpd_pref2 > 0.0) {
+        noise_vec[i] = d_random() - 0.5;
+      } else {
+        noise_vec[i] = 0.0;
+      }
       // Projection Matrix
       for (int j = 0; j < 3; j++) {
         P_times_dist_sqr[i][j] -= d[i] * d[j];

src/core/fene.hpp

@@ -61,7 +61,7 @@ inline int calc_fene_pair_force(Particle *p1, Particle *p2,
   double fac = -iaparams->p.fene.k * dr / ((1.0 - dr*dr*iaparams->p.fene.drmax2i));
   if (fabs(dr) > ROUND_ERROR_PREC) {
      if(len > ROUND_ERROR_PREC) {  /* Regular case */
-	fac /= len ; 
+       fac /= len;
      } else { /* dx[] == 0: the force is undefined. Let's use a random direction */
         for(int i = 0;i < 3;i++) dx[i] = d_random()-0.5;
         fac /= sqrt(sqrlen(dx));

src/core/ghmc.cpp

@@ -266,12 +266,16 @@ void simple_momentum_update() {
     sigmat = sqrt(temperature / (p).p.mass);
 #endif
     for (int j = 0; j < 3; j++) {
-      p.m.v[j] = sigmat * gaussian_random() * time_step;
+      if (sigmat > 0.0) {
+        p.m.v[j] = sigmat * gaussian_random() * time_step;
+      }
 #ifdef ROTATION
 #ifdef ROTATIONAL_INERTIA
       sigmar = sqrt(temperature / p.p.rinertia[j]);
 #endif
-      p.m.omega[j] = sigmar * gaussian_random();
+      if (sigmar > 0.0) {
+        p.m.omega[j] = sigmar * gaussian_random();
+      }
 #endif
     }
   }
@@ -289,14 +293,22 @@ void partial_momentum_update() {
     sigmat = sqrt(temperature / (p).p.mass);
 #endif
     for (int j = 0; j < 3; j++) {
-      p.m.v[j] =
+      if (sigmat > 0.0) {
+        p.m.v[j] =
           cosp * (p.m.v[j]) + sinp * (sigmat * gaussian_random() * time_step);
+      } else {
+        p.m.v[j] = cosp * p.m.v[j];
+      }
 #ifdef ROTATION
 #ifdef ROTATIONAL_INERTIA
       sigmar = sqrt(temperature / p.p.rinertia[j]);
 #endif
-      p.m.omega[j] =
+      if (sigmar > 0.0) {
+        p.m.omega[j] =
           cosp * (p.m.omega[j]) + sinp * (sigmar * gaussian_random());
+      } else {
+        p.m.omega[j] = cosp * p.m.omega[j];
+      }
 #endif
     }
   }
@@ -379,9 +391,6 @@ void ghmc_mc() {
     else
       boltzmann = exp(-beta * boltzmann);
 
-    // fprintf(stderr,"old hamiltonian : %f, new hamiltonian: % f, boltzmann
-    // factor: %f\n",ghmcdata.hmlt_old,ghmcdata.hmlt_new,boltzmann);
-
     if (d_random() < boltzmann) {
       ghmcdata.acc++;
       ghmc_mc_res = GHMC_MOVE_ACCEPT;

src/core/lb.cpp

@@ -3013,18 +3013,22 @@ void calc_particle_lattice_ia() {
 
     /* draw random numbers for local particles */
     for (auto &p : local_cells.particles()) {
+      if (lb_coupl_pref2 > 0.0) {
 #ifdef GAUSSRANDOM
-      p.lc.f_random[0] = lb_coupl_pref2 * gaussian_random();
-      p.lc.f_random[1] = lb_coupl_pref2 * gaussian_random();
-      p.lc.f_random[2] = lb_coupl_pref2 * gaussian_random();
+        p.lc.f_random[0] = lb_coupl_pref2 * gaussian_random();
+        p.lc.f_random[1] = lb_coupl_pref2 * gaussian_random();
+        p.lc.f_random[2] = lb_coupl_pref2 * gaussian_random();
 #elif defined(GAUSSRANDOMCUT)
-      p.lc.f_random[0] = lb_coupl_pref2 * gaussian_random_cut();
-      p.lc.f_random[1] = lb_coupl_pref2 * gaussian_random_cut();
-      p.lc.f_random[2] = lb_coupl_pref2 * gaussian_random_cut();
+        p.lc.f_random[0] = lb_coupl_pref2 * gaussian_random_cut();
+        p.lc.f_random[1] = lb_coupl_pref2 * gaussian_random_cut();
+        p.lc.f_random[2] = lb_coupl_pref2 * gaussian_random_cut();
 #elif defined(FLATNOISE)
-      p.lc.f_random[0] = lb_coupl_pref * (d_random() - 0.5);
-      p.lc.f_random[1] = lb_coupl_pref * (d_random() - 0.5);
-      p.lc.f_random[2] = lb_coupl_pref * (d_random() - 0.5);
+        p.lc.f_random[0] = lb_coupl_pref * (d_random() - 0.5);
+        p.lc.f_random[1] = lb_coupl_pref * (d_random() - 0.5);
+        p.lc.f_random[2] = lb_coupl_pref * (d_random() - 0.5);
+      } else {
+        p.lc.f_random = {0.0, 0.0, 0.0};
+      }
 #else // GAUSSRANDOM
 #error No noise type defined for the CPU LB
 #endif // GAUSSRANDOM

src/core/lbgpu.cpp

@@ -47,6 +47,7 @@
 #include <cstdlib>
 #include <ctime>
 #include <cstring>
+#include <random>
 
 #ifndef D3Q19
 #error The implementation only works for D3Q19 so far!
@@ -218,11 +219,7 @@ void lb_realloc_particles_gpu(){
 
   lbpar_gpu.number_of_particles = n_part;
   LB_TRACE (printf("#particles realloc\t %u \n", lbpar_gpu.number_of_particles));
-  //fprintf(stderr, "%u \t \n", lbpar_gpu.number_of_particles);
-  /**-----------------------------------------------------*/
-  /** allocating of the needed memory for several structs */
-  /**-----------------------------------------------------*/
-  lbpar_gpu.your_seed = (unsigned int)i_random(max_ran);
+  lbpar_gpu.your_seed = (unsigned int)std::random_device{}();
 
   LB_TRACE (fprintf(stderr,"test your_seed %u \n", lbpar_gpu.your_seed));
 
@@ -232,12 +229,7 @@ void lb_realloc_particles_gpu(){
 /** (Re-)initializes the fluid according to the given value of rho. */
 void lb_reinit_fluid_gpu() {
 
-  //lbpar_gpu.your_seed = (unsigned int)i_random(max_ran);
   lb_reinit_parameters_gpu();
-//#ifdef SHANCHEN
-//  lb_calc_particle_lattice_ia_gpu();
-//  copy_forces_from_GPU();
-//#endif 
   if(lbpar_gpu.number_of_nodes != 0){
     lb_reinit_GPU(&lbpar_gpu);
     lbpar_gpu.reinit = 1;

src/core/random.cpp

@@ -38,6 +38,9 @@ std::mt19937 generator;
 std::normal_distribution<double> normal_distribution(0, 1);
 std::uniform_real_distribution<double> uniform_real_distribution(0, 1);
 
+bool user_has_seeded = false;
+bool unseeded_error_thrown = false;
+
 /** Local functions */
 
 /**
@@ -70,7 +73,7 @@ int get_state_size_of_generator() {
 /** Communication */
 
 void mpi_random_seed_slave(int pnode, int cnt) {
-  int this_seed;
+  int this_seed;  user_has_seeded=true;
 
   MPI_Scatter(nullptr, 1, MPI_INT, &this_seed, 1, MPI_INT, 0, comm_cart);
 
@@ -80,6 +83,7 @@ void mpi_random_seed_slave(int pnode, int cnt) {
 
 void mpi_random_seed(int cnt, vector<int> &seeds) {
   int this_seed;
+  user_has_seeded = true;
   mpi_call(mpi_random_seed_slave, -1, cnt);
 
   MPI_Scatter(&seeds[0], 1, MPI_INT, &this_seed, 1, MPI_INT, 0, comm_cart);
@@ -90,13 +94,15 @@ void mpi_random_seed(int cnt, vector<int> &seeds) {
 }
 
 void mpi_random_set_stat_slave(int, int) {
+  user_has_seeded = true;
   string msg;
   mpiCallbacks().comm().recv(0, SOME_TAG, msg);
 
   set_state(msg);
 }
 
 void mpi_random_set_stat(const vector<string> &stat) {
+  user_has_seeded = true;
   mpi_call(mpi_random_set_stat_slave, 0, 0);
 
   for (int i = 1; i < n_nodes; i++) {