diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc
index 95bd75d..6d814bd 100644
--- a/models/ecolab_model.cc
+++ b/models/ecolab_model.cc
@@ -32,26 +32,6 @@ using namespace classdesc;
 // TODO - move this into main library
 namespace
 {
-  /* Rounding function, randomly round up or down, in the range 0..INT_MAX */
-  inline int ROUND(double x) 
-  {
-    double dum;
-    if (x<0) x=0;
-    if (x>INT_MAX-1) x=INT_MAX-1;
-    return std::fabs(std::modf(x,&dum)) > array_urand.rand() ?
-      (int)x+1 : (int)x;
-  }
-
-  inline int ROUND(float x) {return ROUND(double(x));}
-  
-  template <class E>
-  inline array<int> ROUND(const E& x)
-  {
-    array<int> r(x.size());
-    for (size_t i=0; i<x.size(); i++)
-      r[i]=ROUND(x[i]);
-    return r;
-  }
 }
 
 
@@ -77,6 +57,27 @@ void PanmicticModel::generate(unsigned niter)
   tstep+=niter;
 } 
 
+  /* Rounding function, randomly round up or down, in the range 0..INT_MAX */
+  int EcolabPoint::ROUND(double x) 
+  {
+    double dum;
+    if (x<0) x=0;
+    if (x>std::numeric_limits<int>::max()-1) x=std::numeric_limits<int>::max()-1;
+    return std::fabs(std::modf(x,&dum))*(rand.max()-rand.min()) > (rand()-rand.min()) ?
+      (int)x+1 : (int)x;
+  }
+
+//inline int ROUND(float x) {return ROUND(double(x));}
+  
+  template <class E>
+  inline array<int> EcolabPoint::ROUND(const E& x)
+  {
+    array<int> r(x.size());
+    for (size_t i=0; i<x.size(); i++)
+      r[i]=ROUND(x[i]);
+    return r;
+  }
+
 void EcolabPoint::generate(unsigned niter, const ModelData& model)
 { 
   array<double> n(density);
diff --git a/models/ecolab_model.h b/models/ecolab_model.h
index 5741dc3..f918e2b 100644
--- a/models/ecolab_model.h
+++ b/models/ecolab_model.h
@@ -11,6 +11,7 @@
 #include <ecolab.h>
 using classdesc::Object;
 
+#include <random>
 #include <vector>
 #include <pack_stl.h>
 
@@ -65,6 +66,10 @@ struct EcolabPoint
   void condense(const array<bool>& mask, size_t mask_true);
   array<int> mutate(const array<double>&);
   unsigned nsp() const; ///< number of living species in this cell
+  /// Rounding function, randomly round up or down, in the range 0..INT_MAX
+  int ROUND(double x); 
+  template <class E> array<int> ROUND(const E& x);
+  Exclude<std::mt19937> rand; // random number generator
 };
 
 struct PanmicticModel: public ModelData, public EcolabPoint, public ecolab::Model<PanmicticModel>
@@ -72,6 +77,7 @@ struct PanmicticModel: public ModelData, public EcolabPoint, public ecolab::Mode
   ConnectionPlot connectionPlot;
   void updateConnectionPlot() {connectionPlot.update(density,interaction);}
   void makeConsistent() {ModelData::makeConsistent(density.size());}
+  void seed(double x) {rand.seed(x);}
   void generate(unsigned niter);
   void generate() {generate(1);}
   void condense();
@@ -93,6 +99,7 @@ class SpatialModel: public ModelData, public EcolabGraph<EcolabCell>, public eco
   }
   array<unsigned> nsp() const;
   void makeConsistent();
+  void seed(double x) {forAll([=](EcolabCell& cell){cell.rand.seed(x);});}
   void generate(unsigned niter);
   void generate() {generate(1);}
   void condense();