diff --git a/applications/test/BinSum/Test-BinSum.C b/applications/test/BinSum/Test-BinSum.C
index 873e7eb893..6d0abfff28 100644
--- a/applications/test/BinSum/Test-BinSum.C
+++ b/applications/test/BinSum/Test-BinSum.C
@@ -46,7 +46,7 @@ int main(int argc, char *argv[])
     scalarField samples(10000000);
     forAll(samples, i)
     {
-        samples[i] = rndGen.scalar01();
+        samples[i] = rndGen.sample01<scalar>();
     }
 
     const scalar min = 0;
diff --git a/applications/test/Distribution/Test-Distribution.C b/applications/test/Distribution/Test-Distribution.C
index f893e1e4bc..647c7a7d42 100644
--- a/applications/test/Distribution/Test-Distribution.C
+++ b/applications/test/Distribution/Test-Distribution.C
@@ -71,7 +71,7 @@ int main(int argc, char *argv[])
 
         for (label i = 0; i < randomDistributionTestSize; i++)
         {
-            dS.add(2.5*R.GaussNormal() + 8.5);
+            dS.add(2.5*R.GaussNormal<scalar>() + 8.5);
         }
 
         Info<< "Mean " << dS.mean() << nl
@@ -90,7 +90,7 @@ int main(int argc, char *argv[])
 
         for (label i = 0; i < randomDistributionTestSize; i++)
         {
-            dS2.add(1.5*R.GaussNormal() -6.0);
+            dS2.add(1.5*R.GaussNormal<scalar>() -6.0);
         }
 
         Info<< "Mean " << dS2.mean() << nl
@@ -121,7 +121,7 @@ int main(int argc, char *argv[])
 
         for (label i = 0; i < randomDistributionTestSize; i++)
         {
-            dS.add(R.scalar01() + 10*Pstream::myProcNo());
+            dS.add(R.sample01<scalar>() + 10*Pstream::myProcNo());
         }
 
         Pout<< "Mean " << dS.mean() << nl
@@ -155,7 +155,7 @@ int main(int argc, char *argv[])
 
         for (label i = 0; i < randomDistributionTestSize; i++)
         {
-            dV.add(R.vector01());
+            dV.add(R.sample01<vector>());
 
             // Adding separate GaussNormal components with component
             // weights
@@ -164,9 +164,9 @@ int main(int argc, char *argv[])
             (
                 vector
                 (
-                    R.GaussNormal()*3.0 + 1.5,
-                    R.GaussNormal()*0.25 + 4.0,
-                    R.GaussNormal()*3.0 - 1.5
+                    R.GaussNormal<scalar>()*3.0 + 1.5,
+                    R.GaussNormal<scalar>()*0.25 + 4.0,
+                    R.GaussNormal<scalar>()*3.0 - 1.5
                 ),
                 vector(1.0, 2.0, 5.0)
             );
@@ -225,7 +225,7 @@ int main(int argc, char *argv[])
 
         for (label i = 0; i < randomDistributionTestSize; i++)
         {
-            dT.add(R.tensor01());
+            dT.add(R.sample01<tensor>());
         }
 
         Info<< "Mean " << dT.mean() << nl
@@ -249,7 +249,7 @@ int main(int argc, char *argv[])
 
         for (label i = 0; i < randomDistributionTestSize; i++)
         {
-            dSyT.add(R.symmTensor01());
+            dSyT.add(R.sample01<symmTensor>());
         }
 
         Info<< "Mean " << dSyT.mean() << nl
@@ -273,7 +273,7 @@ int main(int argc, char *argv[])
 
         for (label i = 0; i < randomDistributionTestSize; i++)
         {
-            dSpT.add(R.sphericalTensor01());
+            dSpT.add(R.sample01<sphericalTensor>());
         }
 
         Info<< "Mean " << dSpT.mean() << nl
diff --git a/applications/test/Polynomial/Test-Polynomial.C b/applications/test/Polynomial/Test-Polynomial.C
index 352d15f96d..6a6fb8bb45 100644
--- a/applications/test/Polynomial/Test-Polynomial.C
+++ b/applications/test/Polynomial/Test-Polynomial.C
@@ -152,7 +152,7 @@ int main(int argc, char *argv[])
     Random rnd(123456);
     for (int i=0; i<10; i++)
     {
-        scalar x = rnd.scalar01()*100;
+        scalar x = rnd.sample01<scalar>()*100;
 
         scalar px = polyValue(x);
         scalar ipx = intPolyValue(x);
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
index 190e493322..929eeb306b 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
@@ -293,9 +293,9 @@ inline Foam::point Foam::conformalVoronoiMesh::perturbPoint
 
     scalar pert = 1e-12*defaultCellSize();
 
-    perturbedPt.x() += pert*(rndGen_.scalar01() - 0.5);
-    perturbedPt.y() += pert*(rndGen_.scalar01() - 0.5);
-    perturbedPt.z() += pert*(rndGen_.scalar01() - 0.5);
+    perturbedPt.x() += pert*(rndGen_.sample01<scalar>() - 0.5);
+    perturbedPt.y() += pert*(rndGen_.sample01<scalar>() - 0.5);
+    perturbedPt.z() += pert*(rndGen_.sample01<scalar>() - 0.5);
 
     return perturbedPt;
 }
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C
index e2d1001504..1acd01d9bf 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C
@@ -553,11 +553,11 @@ bool Foam::autoDensity::fillBox
                       + vector
                         (
                             delta.x()
-                           *(i + 0.5 + 0.1*(rndGen().scalar01() - 0.5)),
+                           *(i + 0.5 + 0.1*(rndGen().sample01<scalar>() - 0.5)),
                             delta.y()
-                           *(j + 0.5 + 0.1*(rndGen().scalar01() - 0.5)),
+                           *(j + 0.5 + 0.1*(rndGen().sample01<scalar>() - 0.5)),
                             delta.z()
-                           *(k + 0.5 + 0.1*(rndGen().scalar01() - 0.5))
+                           *(k + 0.5 + 0.1*(rndGen().sample01<scalar>() - 0.5))
                         );
                 }
             }
@@ -662,7 +662,7 @@ bool Foam::autoDensity::fillBox
                 // TODO - is there a lot of cost in the 1/density calc?  Could
                 // assess on
                 //    (1/maxDensity)/(1/localDensity) = minVolume/localVolume
-                if (localDensity/maxDensity > rndGen().scalar01())
+                if (localDensity/maxDensity > rndGen().sample01<scalar>())
                 {
                     scalar localVolume = 1/localDensity;
 
@@ -675,7 +675,7 @@ bool Foam::autoDensity::fillBox
                         scalar addProbability =
                            (totalVolume - volumeAdded)/localVolume;
 
-                        scalar r = rndGen().scalar01();
+                        scalar r = rndGen().sample01<scalar>();
 
                         if (debug)
                         {
@@ -729,7 +729,7 @@ bool Foam::autoDensity::fillBox
         {
             trialPoints++;
 
-            point p = min + cmptMultiply(span, rndGen().vector01());
+            point p = min + cmptMultiply(span, rndGen().sample01<vector>());
 
             scalar localSize = cellShapeControls().cellSize(p);
 
@@ -785,7 +785,7 @@ bool Foam::autoDensity::fillBox
 
                 // Accept possible placements proportional to the relative local
                 // density
-                if (localDensity/maxDensity > rndGen().scalar01())
+                if (localDensity/maxDensity > rndGen().sample01<scalar>())
                 {
                     scalar localVolume = 1/localDensity;
 
@@ -798,7 +798,7 @@ bool Foam::autoDensity::fillBox
                         scalar addProbability =
                             (totalVolume - volumeAdded)/localVolume;
 
-                        scalar r = rndGen().scalar01();
+                        scalar r = rndGen().sample01<scalar>();
 
                         if (debug)
                         {
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/bodyCentredCubic/bodyCentredCubic.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/bodyCentredCubic/bodyCentredCubic.C
index 17a4285d86..8991e7ae9d 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/bodyCentredCubic/bodyCentredCubic.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/bodyCentredCubic/bodyCentredCubic.C
@@ -129,9 +129,9 @@ List<Vb::Point> bodyCentredCubic::initialPoints() const
 
                 if (randomiseInitialGrid_)
                 {
-                    pA.x() += pert*(rndGen().scalar01() - 0.5);
-                    pA.y() += pert*(rndGen().scalar01() - 0.5);
-                    pA.z() += pert*(rndGen().scalar01() - 0.5);
+                    pA.x() += pert*(rndGen().sample01<scalar>() - 0.5);
+                    pA.y() += pert*(rndGen().sample01<scalar>() - 0.5);
+                    pA.z() += pert*(rndGen().sample01<scalar>() - 0.5);
                 }
 
                 if (Pstream::parRun())
@@ -150,9 +150,9 @@ List<Vb::Point> bodyCentredCubic::initialPoints() const
 
                 if (randomiseInitialGrid_)
                 {
-                    pB.x() += pert*(rndGen().scalar01() - 0.5);
-                    pB.y() += pert*(rndGen().scalar01() - 0.5);
-                    pB.z() += pert*(rndGen().scalar01() - 0.5);
+                    pB.x() += pert*(rndGen().sample01<scalar>() - 0.5);
+                    pB.y() += pert*(rndGen().sample01<scalar>() - 0.5);
+                    pB.z() += pert*(rndGen().sample01<scalar>() - 0.5);
                 }
 
                 if (Pstream::parRun())
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/faceCentredCubic/faceCentredCubic.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/faceCentredCubic/faceCentredCubic.C
index 9fd36aee3a..36052a3bbb 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/faceCentredCubic/faceCentredCubic.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/faceCentredCubic/faceCentredCubic.C
@@ -127,9 +127,9 @@ List<Vb::Point> faceCentredCubic::initialPoints() const
 
                 if (randomiseInitialGrid_)
                 {
-                    p.x() += pert*(rndGen().scalar01() - 0.5);
-                    p.y() += pert*(rndGen().scalar01() - 0.5);
-                    p.z() += pert*(rndGen().scalar01() - 0.5);
+                    p.x() += pert*(rndGen().sample01<scalar>() - 0.5);
+                    p.y() += pert*(rndGen().sample01<scalar>() - 0.5);
+                    p.z() += pert*(rndGen().sample01<scalar>() - 0.5);
                 }
 
                 if (Pstream::parRun())
@@ -155,9 +155,9 @@ List<Vb::Point> faceCentredCubic::initialPoints() const
 
                 if (randomiseInitialGrid_)
                 {
-                    p.x() += pert*(rndGen().scalar01() - 0.5);
-                    p.y() += pert*(rndGen().scalar01() - 0.5);
-                    p.z() += pert*(rndGen().scalar01() - 0.5);
+                    p.x() += pert*(rndGen().sample01<scalar>() - 0.5);
+                    p.y() += pert*(rndGen().sample01<scalar>() - 0.5);
+                    p.z() += pert*(rndGen().sample01<scalar>() - 0.5);
                 }
 
                 if (Pstream::parRun())
@@ -183,9 +183,9 @@ List<Vb::Point> faceCentredCubic::initialPoints() const
 
                 if (randomiseInitialGrid_)
                 {
-                    p.x() += pert*(rndGen().scalar01() - 0.5);
-                    p.y() += pert*(rndGen().scalar01() - 0.5);
-                    p.z() += pert*(rndGen().scalar01() - 0.5);
+                    p.x() += pert*(rndGen().sample01<scalar>() - 0.5);
+                    p.y() += pert*(rndGen().sample01<scalar>() - 0.5);
+                    p.z() += pert*(rndGen().sample01<scalar>() - 0.5);
                 }
 
                 if (Pstream::parRun())
@@ -211,9 +211,9 @@ List<Vb::Point> faceCentredCubic::initialPoints() const
 
                 if (randomiseInitialGrid_)
                 {
-                    p.x() += pert*(rndGen().scalar01() - 0.5);
-                    p.y() += pert*(rndGen().scalar01() - 0.5);
-                    p.z() += pert*(rndGen().scalar01() - 0.5);
+                    p.x() += pert*(rndGen().sample01<scalar>() - 0.5);
+                    p.y() += pert*(rndGen().sample01<scalar>() - 0.5);
+                    p.z() += pert*(rndGen().sample01<scalar>() - 0.5);
                 }
 
                 if (Pstream::parRun())
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
index 38742a6203..a09153cf20 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
@@ -198,9 +198,15 @@ List<Vb::Point> pointFile::initialPoints() const
 
             if (randomiseInitialGrid_)
             {
-                p.x() += randomPerturbationCoeff_*(rndGen().scalar01() - 0.5);
-                p.y() += randomPerturbationCoeff_*(rndGen().scalar01() - 0.5);
-                p.z() += randomPerturbationCoeff_*(rndGen().scalar01() - 0.5);
+                p.x() +=
+                    randomPerturbationCoeff_
+                   *(rndGen().sample01<scalar>() - 0.5);
+                p.y() +=
+                    randomPerturbationCoeff_
+                   *(rndGen().sample01<scalar>() - 0.5);
+                p.z() +=
+                    randomPerturbationCoeff_
+                   *(rndGen().sample01<scalar>() - 0.5);
             }
 
             initialPoints.append(Vb::Point(p.x(), p.y(), p.z()));
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C
index 103a09c23a..7727153e5b 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C
@@ -86,9 +86,9 @@ void rayShooting::splitLine
         {
             Foam::point newPt
             (
-                midPoint.x() + pert*(rndGen().scalar01() - 0.5),
-                midPoint.y() + pert*(rndGen().scalar01() - 0.5),
-                midPoint.z() + pert*(rndGen().scalar01() - 0.5)
+                midPoint.x() + pert*(rndGen().sample01<scalar>() - 0.5),
+                midPoint.y() + pert*(rndGen().sample01<scalar>() - 0.5),
+                midPoint.z() + pert*(rndGen().sample01<scalar>() - 0.5)
             );
 
             if
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C
index 319db49293..e7d8d8dcf1 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C
@@ -128,9 +128,9 @@ List<Vb::Point> uniformGrid::initialPoints() const
 
                 if (randomiseInitialGrid_)
                 {
-                    p.x() += pert*(rndGen().scalar01() - 0.5);
-                    p.y() += pert*(rndGen().scalar01() - 0.5);
-                    p.z() += pert*(rndGen().scalar01() - 0.5);
+                    p.x() += pert*(rndGen().sample01<scalar>() - 0.5);
+                    p.y() += pert*(rndGen().sample01<scalar>() - 0.5);
+                    p.z() += pert*(rndGen().sample01<scalar>() - 0.5);
                 }
 
                 if
diff --git a/applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/CV2D.C b/applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/CV2D.C
index 92d08f29d6..1a540eaf8d 100644
--- a/applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/CV2D.C
+++ b/applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/CV2D.C
@@ -264,8 +264,8 @@ void Foam::CV2D::insertGrid()
 
             if (meshControls().randomiseInitialGrid())
             {
-                p.x() += pert*(rndGen.scalar01() - 0.5);
-                p.y() += pert*(rndGen.scalar01() - 0.5);
+                p.x() += pert*(rndGen.sample01<scalar>() - 0.5);
+                p.y() += pert*(rndGen.sample01<scalar>() - 0.5);
             }
 
             if (qSurf_.wellInside(p, 0.5*meshControls().minCellSize2()))
diff --git a/applications/utilities/postProcessing/miscellaneous/pdfPlot/createFields.H b/applications/utilities/postProcessing/miscellaneous/pdfPlot/createFields.H
index 85ec143ca9..68229a7a93 100644
--- a/applications/utilities/postProcessing/miscellaneous/pdfPlot/createFields.H
+++ b/applications/utilities/postProcessing/miscellaneous/pdfPlot/createFields.H
@@ -20,7 +20,7 @@
     const fileName pdfPath = runTime.path()/"pdf";
     mkDir(pdfPath);
 
-    cachedRandom rndGen(label(0), -1);
+    Random rndGen;
 
     autoPtr<distributionModels::distributionModel> p
     (
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C b/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C
index abc3d321d2..02d0b0c8e6 100644
--- a/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C
+++ b/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C
@@ -783,7 +783,7 @@ labelPair edgeIntersectionsAndShuffleCGAL
                     const edge& e = edges[edgeI];
                     forAll(e, eI)
                     {
-                        vector d = rndGen.vector01()-p05;
+                        vector d = rndGen.sample01<vector>() - p05;
                         surf1Points[mp[e[eI]]] += surf1PointTol[e[eI]]*d;
                     }
                 }
diff --git a/applications/utilities/surface/surfaceInertia/surfaceInertia.C b/applications/utilities/surface/surfaceInertia/surfaceInertia.C
index b0808d239d..fb45390688 100644
--- a/applications/utilities/surface/surfaceInertia/surfaceInertia.C
+++ b/applications/utilities/surface/surfaceInertia/surfaceInertia.C
@@ -125,9 +125,9 @@ int main(int argc, char *argv[])
             << "No eigenValues found, shape may have symmetry, "
             << "perturbing inertia tensor diagonal" << endl;
 
-        J.xx() *= 1.0 + SMALL*rand.scalar01();
-        J.yy() *= 1.0 + SMALL*rand.scalar01();
-        J.zz() *= 1.0 + SMALL*rand.scalar01();
+        J.xx() *= 1.0 + SMALL*rand.sample01<scalar>();
+        J.yy() *= 1.0 + SMALL*rand.sample01<scalar>();
+        J.zz() *= 1.0 + SMALL*rand.sample01<scalar>();
 
         eVal = eigenValues(J);
 
diff --git a/src/OSspecific/POSIX/POSIX.C b/src/OSspecific/POSIX/POSIX.C
index 5594acdee8..d444b83a0b 100644
--- a/src/OSspecific/POSIX/POSIX.C
+++ b/src/OSspecific/POSIX/POSIX.C
@@ -59,7 +59,6 @@ Description
 #include <link.h>
 
 #include <netinet/in.h>
-
 #ifdef USE_RANDOM
     #include <climits>
     #if INT_MAX    != 2147483647
@@ -1430,6 +1429,15 @@ Foam::fileNameList Foam::dlLoaded()
     return libs;
 }
 
+Foam::label Foam::osRandomBufferSize()
+{
+    #ifdef USE_RANDOM
+    return sizeof(random_data);
+    #else
+    return sizeof(drand48_data);
+    #endif
+}
+
 
 void Foam::osRandomSeed(const label seed)
 {
@@ -1461,4 +1469,42 @@ Foam::scalar Foam::osRandomDouble()
 }
 
 
+void Foam::osRandomSeed(const label seed, List<char>& buffer)
+{
+    #ifdef USE_RANDOM
+    srandom_r((unsigned int)seed, reinterpret_cast<random_data*>(buffer.begin()));
+    #else
+    srand48_r(seed, reinterpret_cast<drand48_data*>(buffer.begin()));
+    #endif
+}
+
+
+Foam::label Foam::osRandomInteger(List<char>& buffer)
+{
+    #ifdef USE_RANDOM
+    int32_t result;
+    random_r(reinterpret_cast<random_data*>(buffer.begin()), &result);
+    return result;
+    #else
+    long result;
+    lrand48_r(reinterpret_cast<drand48_data*>(buffer.begin()), &result);
+    return result;
+    #endif
+}
+
+
+Foam::scalar Foam::osRandomDouble(List<char>& buffer)
+{
+    #ifdef USE_RANDOM
+    int32_t result;
+    random_r(reinterpret_cast<random_data*>(buffer.begin()), &result);
+    return (scalar)result/INT_MAX;
+    #else
+    double result;
+    drand48_r(reinterpret_cast<drand48_data*>(buffer.begin()), &result);
+    return result;
+    #endif
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index 14af54effd..5636f7c937 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -111,7 +111,6 @@ $(sha1)/SHA1.C
 $(sha1)/SHA1Digest.C
 
 primitives/random/Random/Random.C
-primitives/random/cachedRandom/cachedRandom.C
 
 ranges = primitives/ranges
 $(ranges)/labelRange/labelRange.C
diff --git a/src/OpenFOAM/include/OSspecific.H b/src/OpenFOAM/include/OSspecific.H
index 3d5e3bbc7b..6dee7af27b 100644
--- a/src/OpenFOAM/include/OSspecific.H
+++ b/src/OpenFOAM/include/OSspecific.H
@@ -245,6 +245,18 @@ label osRandomInteger();
 //- Return random double precision (uniform distribution between 0 and 1)
 scalar osRandomDouble();
 
+//- Return the size of the buffer for re-entrant random number generator
+label osRandomBufferSize();
+
+//- Seed random number generator.
+void osRandomSeed(const label seed, List<char>& buffer);
+
+//- Return random integer (uniform distribution between 0 and 2^31)
+label osRandomInteger(List<char>& buffer);
+
+//- Return random double precision (uniform distribution between 0 and 1)
+scalar osRandomDouble(List<char>& buffer);
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H
index badea19542..2369e77ca8 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H
@@ -42,7 +42,6 @@ SourceFiles
 #include "point.H"
 #include "primitiveFieldsFwd.H"
 #include "pointHit.H"
-#include "cachedRandom.H"
 #include "Random.H"
 #include "FixedList.H"
 #include "UList.H"
@@ -239,10 +238,6 @@ public:
             //  uniform distribution
             inline Point randomPoint(Random& rndGen) const;
 
-            //- Return a random point in the tetrahedron from a
-            //  uniform distribution
-            inline Point randomPoint(cachedRandom& rndGen) const;
-
             //- Calculate the barycentric coordinates of the given
             //  point, in the same order as a, b, c, d.  Returns the
             //  determinant of the solution.
diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
index 2c14d16121..d95f42dc1e 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
@@ -249,42 +249,6 @@ inline Point Foam::tetrahedron<Point, PointRef>::randomPoint
     // Adapted from
     // http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html
 
-    scalar s = rndGen.scalar01();
-    scalar t = rndGen.scalar01();
-    scalar u = rndGen.scalar01();
-
-    if (s + t > 1.0)
-    {
-        s = 1.0 - s;
-        t = 1.0 - t;
-    }
-
-    if (t + u > 1.0)
-    {
-        scalar tmp = u;
-        u = 1.0 - s - t;
-        t = 1.0 - tmp;
-    }
-    else if (s + t + u > 1.0)
-    {
-        scalar tmp = u;
-        u = s + t + u - 1.0;
-        s = 1.0 - t - tmp;
-    }
-
-    return (1 - s - t - u)*a_ + s*b_ + t*c_ + u*d_;
-}
-
-
-template<class Point, class PointRef>
-inline Point Foam::tetrahedron<Point, PointRef>::randomPoint
-(
-    cachedRandom& rndGen
-) const
-{
-    // Adapted from
-    // http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html
-
     scalar s = rndGen.sample01<scalar>();
     scalar t = rndGen.sample01<scalar>();
     scalar u = rndGen.sample01<scalar>();
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
index ea15519517..b97483ca27 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
@@ -40,7 +40,6 @@ SourceFiles
 #include "tensor.H"
 #include "pointHit.H"
 #include "Random.H"
-#include "cachedRandom.H"
 #include "FixedList.H"
 #include "UList.H"
 #include "linePointRef.H"
@@ -238,10 +237,6 @@ public:
             //  distribution
             inline Point randomPoint(Random& rndGen) const;
 
-            //- Return a random point on the triangle from a uniform
-            //  distribution
-            inline Point randomPoint(cachedRandom& rndGen) const;
-
             //- Calculate the barycentric coordinates of the given
             //  point, in the same order as a, b, c.  Returns the
             //  determinant of the solution.
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
index eee54fd92c..605e10589b 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
@@ -249,25 +249,6 @@ inline Point Foam::triangle<Point, PointRef>::randomPoint(Random& rndGen) const
     // from "Graphics Gems", Academic Press, 1990
     // http://tog.acm.org/GraphicsGems/gems/TriPoints.c
 
-    scalar s = rndGen.scalar01();
-
-    scalar t = sqrt(rndGen.scalar01());
-
-    return (1 - t)*a_ + (1 - s)*t*b_ + s*t*c_;
-}
-
-
-template<class Point, class PointRef>
-inline Point Foam::triangle<Point, PointRef>::randomPoint
-(
-    cachedRandom& rndGen
-) const
-{
-    // Generating Random Points in Triangles
-    // by Greg Turk
-    // from "Graphics Gems", Academic Press, 1990
-    // http://tog.acm.org/GraphicsGems/gems/TriPoints.c
-
     scalar s = rndGen.sample01<scalar>();
 
     scalar t = sqrt(rndGen.sample01<scalar>());
diff --git a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H
index edc42983df..7c27f9fc26 100644
--- a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H
+++ b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H
@@ -335,7 +335,7 @@ public:
             label distanceCmp(const point& pt, const treeBoundBox& other) const;
 
             //- Return slightly wider bounding box
-            //  Extends all dimensions with s*span*Random::scalar01()
+            //  Extends all dimensions with s*span*Random::sample01<scalar>()
             //  and guarantees in any direction s*mag(span) minimum width
             inline treeBoundBox extend(Random& rndGen, const scalar s) const;
 
diff --git a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H
index 4b073fd975..750a6737ff 100644
--- a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H
+++ b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H
@@ -336,8 +336,8 @@ inline Foam::treeBoundBox Foam::treeBoundBox::extend
         newSpan[dir] = Foam::max(newSpan[dir], minSpan);
     }
 
-    bb.min() -= cmptMultiply(s * rndGen.vector01(), newSpan);
-    bb.max() += cmptMultiply(s * rndGen.vector01(), newSpan);
+    bb.min() -= cmptMultiply(s*rndGen.sample01<vector>(), newSpan);
+    bb.max() += cmptMultiply(s*rndGen.sample01<vector>(), newSpan);
 
     return bb;
 }
diff --git a/src/OpenFOAM/primitives/random/Random/Random.C b/src/OpenFOAM/primitives/random/Random/Random.C
index 3fa14012c1..c17f5c2668 100644
--- a/src/OpenFOAM/primitives/random/Random/Random.C
+++ b/src/OpenFOAM/primitives/random/Random/Random.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,173 +25,237 @@ License
 
 #include "Random.H"
 #include "OSspecific.H"
+#include "PstreamReduceOps.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * //
 
-#if INT_MAX    != 2147483647
-#    error "INT_MAX    != 2147483647"
-#    error "The random number generator may not work!"
-#endif
+Foam::scalar Foam::Random::scalar01()
+{
+    return osRandomDouble(buffer_);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::Random::Random(const label seed)
+:
+    buffer_(osRandomBufferSize()),
+    seed_(seed),
+    hasGaussSample_(false),
+    gaussSample_(0)
+{
+    // Initialise the random number generator
+    osRandomSeed(seed_, buffer_);
+}
+
+
+Foam::Random::Random(const Random& r, const bool reset)
+:
+    buffer_(r.buffer_),
+    seed_(r.seed_),
+    hasGaussSample_(r.hasGaussSample_),
+    gaussSample_(r.gaussSample_)
+{
+    if (reset)
+    {
+        hasGaussSample_ = false;
+        gaussSample_ = 0;
+
+        // Re-initialise the samples
+        osRandomSeed(seed_, buffer_);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::Random::~Random()
+{}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::Random::Random(const label seed)
+void Foam::Random::reset(const label seed)
 {
-    if (seed > 1)
-    {
-        Seed = seed;
-    }
-    else
-    {
-        Seed = 1;
-    }
-
-    osRandomSeed(Seed);
+    seed_ = seed;
+    osRandomSeed(seed_, buffer_);
 }
 
 
 int Foam::Random::bit()
 {
-    if (osRandomInteger() > INT_MAX/2)
+    return osRandomInteger(buffer_) & 1;
+}
+
+
+template<>
+Foam::scalar Foam::Random::sample01()
+{
+    return scalar01();
+}
+
+
+template<>
+Foam::label Foam::Random::sample01()
+{
+    return round(scalar01());
+}
+
+
+template<>
+Foam::scalar Foam::Random::GaussNormal()
+{
+    if (hasGaussSample_)
     {
-        return 1;
+        hasGaussSample_ = false;
+        return gaussSample_;
     }
     else
     {
-        return 0;
-    }
-}
-
-
-Foam::scalar Foam::Random::scalar01()
-{
-    return osRandomDouble();
-}
-
-
-Foam::vector Foam::Random::vector01()
-{
-    vector rndVec;
-    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
-    {
-        rndVec.component(cmpt) = scalar01();
-    }
-
-    return rndVec;
-}
-
-
-Foam::sphericalTensor Foam::Random::sphericalTensor01()
-{
-    sphericalTensor rndTen;
-    rndTen.ii() = scalar01();
-
-    return rndTen;
-}
-
-
-Foam::symmTensor Foam::Random::symmTensor01()
-{
-    symmTensor rndTen;
-    for (direction cmpt=0; cmpt<symmTensor::nComponents; cmpt++)
-    {
-        rndTen.component(cmpt) = scalar01();
-    }
-
-    return rndTen;
-}
-
-
-Foam::tensor Foam::Random::tensor01()
-{
-    tensor rndTen;
-    for (direction cmpt=0; cmpt<tensor::nComponents; cmpt++)
-    {
-        rndTen.component(cmpt) = scalar01();
-    }
-
-    return rndTen;
-}
-
-
-Foam::label Foam::Random::integer(const label lower, const label upper)
-{
-    return lower + (osRandomInteger() % (upper+1-lower));
-}
-
-
-Foam::vector Foam::Random::position(const vector& start, const vector& end)
-{
-    vector rndVec(start);
-
-    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
-    {
-        rndVec.component(cmpt) +=
-            scalar01()*(end.component(cmpt) - start.component(cmpt));
-    }
-
-    return rndVec;
-}
-
-
-void Foam::Random::randomise(scalar& s)
-{
-     s = scalar01();
-}
-
-
-void Foam::Random::randomise(vector& v)
-{
-    v = vector01();
-}
-
-
-void Foam::Random::randomise(sphericalTensor& st)
-{
-    st = sphericalTensor01();
-}
-
-
-void Foam::Random::randomise(symmTensor& st)
-{
-    st = symmTensor01();
-}
-
-
-void Foam::Random::randomise(tensor& t)
-{
-    t = tensor01();
-}
-
-
-Foam::scalar Foam::Random::GaussNormal()
-{
-    static int iset = 0;
-    static scalar gset;
-    scalar fac, rsq, v1, v2;
-
-    if (iset == 0)
-    {
+        scalar rsq, v1, v2;
         do
         {
-            v1 = 2.0*scalar01() - 1.0;
-            v2 = 2.0*scalar01() - 1.0;
+            v1 = 2*scalar01() - 1;
+            v2 = 2*scalar01() - 1;
             rsq = sqr(v1) + sqr(v2);
-        } while (rsq >= 1.0 || rsq == 0.0);
+        } while (rsq >= 1 || rsq == 0);
 
-        fac = sqrt(-2.0*log(rsq)/rsq);
-        gset = v1*fac;
-        iset = 1;
+        scalar fac = sqrt(-2*log(rsq)/rsq);
+
+        gaussSample_ = v1*fac;
+        hasGaussSample_ = true;
 
         return v2*fac;
     }
-    else
-    {
-        iset = 0;
+}
 
-        return gset;
+
+template<>
+Foam::label Foam::Random::GaussNormal()
+{
+    return round(GaussNormal<scalar>());
+}
+
+
+template<>
+Foam::scalar Foam::Random::position
+(
+    const scalar& start,
+    const scalar& end
+)
+{
+    return start + scalar01()*(end - start);
+}
+
+
+template<>
+Foam::label Foam::Random::position(const label& start, const label& end)
+{
+    return start + round(scalar01()*(end - start));
+}
+
+
+template<>
+Foam::scalar Foam::Random::globalSample01()
+{
+    scalar value = -GREAT;
+
+    if (Pstream::master())
+    {
+        value = scalar01();
     }
+
+    Pstream::scatter(value);
+
+    return value;
+}
+
+
+template<>
+Foam::label Foam::Random::globalSample01()
+{
+    scalar value = -GREAT;
+
+    if (Pstream::master())
+    {
+        value = scalar01();
+    }
+
+    Pstream::scatter(value);
+
+    return round(value);
+}
+
+
+template<>
+Foam::scalar Foam::Random::globalGaussNormal()
+{
+    scalar value = -GREAT;
+
+    if (Pstream::master())
+    {
+        value = GaussNormal<scalar>();
+    }
+
+    Pstream::scatter(value);
+
+    return value;
+}
+
+
+template<>
+Foam::label Foam::Random::globalGaussNormal()
+{
+    scalar value = -GREAT;
+
+    if (Pstream::master())
+    {
+        value = GaussNormal<scalar>();
+    }
+
+    Pstream::scatter(value);
+
+    return round(value);
+}
+
+
+template<>
+Foam::scalar Foam::Random::globalPosition
+(
+    const scalar& start,
+    const scalar& end
+)
+{
+    scalar value = -GREAT;
+
+    if (Pstream::master())
+    {
+        value = scalar01()*(end - start);
+    }
+
+    Pstream::scatter(value);
+
+    return start + value;
+}
+
+
+template<>
+Foam::label Foam::Random::globalPosition
+(
+    const label& start,
+    const label& end
+)
+{
+    label value = labelMin;
+
+    if (Pstream::master())
+    {
+        value = round(scalar01()*(end - start));
+    }
+
+    Pstream::scatter(value);
+
+    return start + value;
 }
 
 
diff --git a/src/OpenFOAM/primitives/random/Random/Random.H b/src/OpenFOAM/primitives/random/Random/Random.H
index 1dd109a235..fb348fdafa 100644
--- a/src/OpenFOAM/primitives/random/Random/Random.H
+++ b/src/OpenFOAM/primitives/random/Random/Random.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,24 +25,29 @@ Class
     Foam::Random
 
 Description
-    Simple random number generator.
+    Random number generator.
+
 
 SourceFiles
+    RandomI.H
     Random.C
+    RandomTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef Random_H
 #define Random_H
 
-#include "vector.H"
-#include "tensor.H"
+#include "List.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+// Forward declaration of classes
+class Random;
+
 /*---------------------------------------------------------------------------*\
                            Class Random Declaration
 \*---------------------------------------------------------------------------*/
@@ -51,60 +56,165 @@ class Random
 {
     // Private data
 
-        label Seed;
+        //- Buffer used by re-entrant random number generator
+        List<char> buffer_;
+
+        //- Initial random number seed
+        label seed_;
+
+        //- Indicator, which tells if there is a stored gaussian sample
+        bool hasGaussSample_;
+
+        //- Stored sample value
+        scalar gaussSample_;
+
+
+    // Private Member Functions
+
+        //- Returns the current sample
+        scalar scalar01();
 
 
 public:
 
-
     // Constructors
 
-        //- Construct given seed
-        Random(const label);
+        //- Construct given seed and sample count
+        Random(const label seed = 123456);
+
+        //- Copy constructor with optional reset of sampleI
+        Random(const Random& r, const bool reset = false);
+
+
+    // Destructor
+    ~Random();
 
 
     // Member functions
 
-        int bit();
+        // Access
 
-        //- Scalar [0..1] (so including 0,1)
-        scalar scalar01();
+            //- Return const access to the initial random number seed
+            inline label seed() const;
 
-        //- Vector with every component scalar01
-        vector vector01();
+            //- Reset the random number generator seed
+            void reset(const label seed);
 
-        //- sphericalTensor with every component scalar01
-        sphericalTensor sphericalTensor01();
 
-        //- symmTensor with every component scalar01
-        symmTensor symmTensor01();
+        // Evaluation
 
-        //- Tensor with every component scalar01
-        tensor tensor01();
+            // Random numbers
 
-        //- Label [lower..upper]
-        label integer(const label lower, const label upper);
+                //- Return a random bit
+                int bit();
 
-        vector position(const vector&, const vector&);
+                //- Return a sample whose components lie in the range 0-1
+                template<class Type>
+                Type sample01();
 
-        void randomise(scalar&);
-        void randomise(vector&);
-        void randomise(sphericalTensor&);
-        void randomise(symmTensor&);
-        void randomise(tensor&);
+                //- Return a sample whose components are normally distributed
+                //  with zero mean and unity variance N(0, 1)
+                template<class Type>
+                Type GaussNormal();
 
-        //- Return a normal Gaussian randon number
-        //  with zero mean and unity variance N(0, 1)
-        scalar GaussNormal();
+                //- Return a sample between start and end
+                template<class Type>
+                Type position(const Type& start, const Type& end);
+
+                //- Randomise value in the range 0-1
+                template<class Type>
+                void randomise01(Type& value);
+
+                //- Shuffle the values in the list
+                template<class Type>
+                void shuffle(UList<Type>& values);
+
+
+            // Global random numbers - consistent across all processors
+
+                //- Return a sample whose components lie in the range 0-1
+                template<class Type>
+                Type globalSample01();
+
+                //- Return a sample whose components are normally distributed
+                //  with zero mean and unity variance N(0, 1)
+                template<class Type>
+                Type globalGaussNormal();
+
+                //- Return a sample between start and end
+                template<class Type>
+                Type globalPosition(const Type& start, const Type& end);
+
+                //- Randomise value in the range 0-1
+                template<class Type>
+                void globalRandomise01(Type& value);
 };
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Template specialisations
+
+template<>
+scalar Random::sample01<scalar>();
+
+template<>
+label Random::sample01<label>();
+
+template<>
+scalar Random::GaussNormal<scalar>();
+
+template<>
+label Random::GaussNormal<label>();
+
+template<>
+scalar Random::position<scalar>
+(
+    const scalar& start,
+    const scalar& end
+);
+
+template<>
+label Random::position<label>(const label& start, const label& end);
+
+template<>
+scalar Random::globalSample01<scalar>();
+
+template<>
+label Random::globalSample01<label>();
+
+template<>
+scalar Random::globalGaussNormal<scalar>();
+
+template<>
+label Random::globalGaussNormal<label>();
+
+template<>
+scalar Random::globalPosition<scalar>
+(
+    const scalar& start,
+    const scalar& end
+);
+
+template<>
+label Random::globalPosition<label>(const label& start, const label& end);
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "RandomI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "RandomTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandomI.H b/src/OpenFOAM/primitives/random/Random/RandomI.H
similarity index 80%
rename from src/OpenFOAM/primitives/random/cachedRandom/cachedRandomI.H
rename to src/OpenFOAM/primitives/random/Random/RandomI.H
index 8d23cbcede..ef5734dd63 100644
--- a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandomI.H
+++ b/src/OpenFOAM/primitives/random/Random/RandomI.H
@@ -23,32 +23,14 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "cachedRandom.H"
+#include "Random.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-inline Foam::label Foam::cachedRandom::seed() const
+inline Foam::label Foam::Random::seed() const
 {
     return seed_;
 }
 
 
-inline const Foam::scalarList& Foam::cachedRandom::samples() const
-{
-    return samples_;
-}
-
-
-inline Foam::label Foam::cachedRandom::sampleI() const
-{
-    return sampleI_;
-}
-
-
-inline Foam::label& Foam::cachedRandom::sampleI()
-{
-    return sampleI_;
-}
-
-
 // ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandomTemplates.C b/src/OpenFOAM/primitives/random/Random/RandomTemplates.C
similarity index 86%
rename from src/OpenFOAM/primitives/random/cachedRandom/cachedRandomTemplates.C
rename to src/OpenFOAM/primitives/random/Random/RandomTemplates.C
index 66612e765b..461c46e242 100644
--- a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandomTemplates.C
+++ b/src/OpenFOAM/primitives/random/Random/RandomTemplates.C
@@ -23,13 +23,13 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "cachedRandom.H"
+#include "Random.H"
 #include "Pstream.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-Type Foam::cachedRandom::sample01()
+Type Foam::Random::sample01()
 {
     Type value;
     for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
@@ -42,7 +42,7 @@ Type Foam::cachedRandom::sample01()
 
 
 template<class Type>
-Type Foam::cachedRandom::GaussNormal()
+Type Foam::Random::GaussNormal()
 {
     Type value;
     for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
@@ -55,7 +55,7 @@ Type Foam::cachedRandom::GaussNormal()
 
 
 template<class Type>
-Type Foam::cachedRandom::position(const Type& start, const Type& end)
+Type Foam::Random::position(const Type& start, const Type& end)
 {
     Type value(start);
     for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
@@ -69,14 +69,14 @@ Type Foam::cachedRandom::position(const Type& start, const Type& end)
 
 
 template<class Type>
-void Foam::cachedRandom::randomise01(Type& value)
+void Foam::Random::randomise01(Type& value)
 {
     value = sample01<Type>();
 }
 
 
 template<class Type>
-void Foam::cachedRandom::shuffle(UList<Type>& values)
+void Foam::Random::shuffle(UList<Type>& values)
 {
     const label nSample = values.size();
     label posI = nSample - 1;
@@ -93,7 +93,7 @@ void Foam::cachedRandom::shuffle(UList<Type>& values)
 
 
 template<class Type>
-Type Foam::cachedRandom::globalSample01()
+Type Foam::Random::globalSample01()
 {
     Type value = -GREAT*pTraits<Type>::one;
 
@@ -109,7 +109,7 @@ Type Foam::cachedRandom::globalSample01()
 
 
 template<class Type>
-Type Foam::cachedRandom::globalGaussNormal()
+Type Foam::Random::globalGaussNormal()
 {
     Type value = -GREAT*pTraits<Type>::one;
 
@@ -125,7 +125,7 @@ Type Foam::cachedRandom::globalGaussNormal()
 
 
 template<class Type>
-Type Foam::cachedRandom::globalPosition(const Type& start, const Type& end)
+Type Foam::Random::globalPosition(const Type& start, const Type& end)
 {
     Type value = -GREAT*pTraits<Type>::one;
 
@@ -141,7 +141,7 @@ Type Foam::cachedRandom::globalPosition(const Type& start, const Type& end)
 
 
 template<class Type>
-void Foam::cachedRandom::globalRandomise01(Type& value)
+void Foam::Random::globalRandomise01(Type& value)
 {
     value = -GREAT*pTraits<Type>::one;
 
diff --git a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.C b/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.C
deleted file mode 100644
index e2fec45161..0000000000
--- a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.C
+++ /dev/null
@@ -1,299 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "cachedRandom.H"
-#include "OSspecific.H"
-#include "PstreamReduceOps.H"
-
-// * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * //
-
-Foam::scalar Foam::cachedRandom::scalar01()
-{
-    if (sampleI_ < 0)
-    {
-        return osRandomDouble();
-    }
-
-    if (sampleI_ == samples_.size() - 1)
-    {
-        scalar s = samples_[sampleI_];
-        sampleI_ = 0;
-        return s;
-    }
-    else
-    {
-        scalar s = samples_[sampleI_];
-        sampleI_++;
-        return s;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::cachedRandom::cachedRandom(const label seed, const label count)
-:
-    seed_(1),
-    samples_(0),
-    sampleI_(-1),
-    hasGaussSample_(false),
-    gaussSample_(0)
-{
-    if (seed > 1)
-    {
-        seed_ = seed;
-    }
-
-    // Initialise samples
-    osRandomSeed(seed_);
-
-    // Samples will be cached if count > 0
-    if (count > 0)
-    {
-        samples_.setSize(count);
-        forAll(samples_, i)
-        {
-            samples_[i] = osRandomDouble();
-        }
-
-        sampleI_ = 0;
-    }
-}
-
-
-Foam::cachedRandom::cachedRandom(const cachedRandom& cr, const bool reset)
-:
-    seed_(cr.seed_),
-    samples_(cr.samples_),
-    sampleI_(cr.sampleI_),
-    hasGaussSample_(cr.hasGaussSample_),
-    gaussSample_(cr.gaussSample_)
-{
-    //if (sampleI_ == -1)
-    //{
-    //    WarningInFunction
-    //        << "Copy constructor called, but samples not being cached. "
-    //        << "This may lead to non-repeatable behaviour" << endl;
-    //
-    //}
-
-    if (reset)
-    {
-        hasGaussSample_ = false;
-        gaussSample_ = 0;
-
-        if (samples_.size())
-        {
-            sampleI_ = 0;
-        }
-        else
-        {
-            // Re-initialise the samples
-            osRandomSeed(seed_);
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::cachedRandom::~cachedRandom()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<>
-Foam::scalar Foam::cachedRandom::sample01()
-{
-    return scalar01();
-}
-
-
-template<>
-Foam::label Foam::cachedRandom::sample01()
-{
-    return round(scalar01());
-}
-
-
-template<>
-Foam::scalar Foam::cachedRandom::GaussNormal()
-{
-    if (hasGaussSample_)
-    {
-        hasGaussSample_ = false;
-        return gaussSample_;
-    }
-    else
-    {
-        scalar rsq, v1, v2;
-        do
-        {
-            v1 = 2*scalar01() - 1;
-            v2 = 2*scalar01() - 1;
-            rsq = sqr(v1) + sqr(v2);
-        } while (rsq >= 1 || rsq == 0);
-
-        scalar fac = sqrt(-2*log(rsq)/rsq);
-
-        gaussSample_ = v1*fac;
-        hasGaussSample_ = true;
-
-        return v2*fac;
-    }
-}
-
-
-template<>
-Foam::label Foam::cachedRandom::GaussNormal()
-{
-    return round(GaussNormal<scalar>());
-}
-
-
-template<>
-Foam::scalar Foam::cachedRandom::position
-(
-    const scalar& start,
-    const scalar& end
-)
-{
-    return start + scalar01()*(end - start);
-}
-
-
-template<>
-Foam::label Foam::cachedRandom::position(const label& start, const label& end)
-{
-    return start + round(scalar01()*(end - start));
-}
-
-
-template<>
-Foam::scalar Foam::cachedRandom::globalSample01()
-{
-    scalar value = -GREAT;
-
-    if (Pstream::master())
-    {
-        value = scalar01();
-    }
-
-    Pstream::scatter(value);
-
-    return value;
-}
-
-
-template<>
-Foam::label Foam::cachedRandom::globalSample01()
-{
-    scalar value = -GREAT;
-
-    if (Pstream::master())
-    {
-        value = scalar01();
-    }
-
-    Pstream::scatter(value);
-
-    return round(value);
-}
-
-
-template<>
-Foam::scalar Foam::cachedRandom::globalGaussNormal()
-{
-    scalar value = -GREAT;
-
-    if (Pstream::master())
-    {
-        value = GaussNormal<scalar>();
-    }
-
-    Pstream::scatter(value);
-
-    return value;
-}
-
-
-template<>
-Foam::label Foam::cachedRandom::globalGaussNormal()
-{
-    scalar value = -GREAT;
-
-    if (Pstream::master())
-    {
-        value = GaussNormal<scalar>();
-    }
-
-    Pstream::scatter(value);
-
-    return round(value);
-}
-
-
-template<>
-Foam::scalar Foam::cachedRandom::globalPosition
-(
-    const scalar& start,
-    const scalar& end
-)
-{
-    scalar value = -GREAT;
-
-    if (Pstream::master())
-    {
-        value = scalar01()*(end - start);
-    }
-
-    Pstream::scatter(value);
-
-    return start + value;
-}
-
-
-template<>
-Foam::label Foam::cachedRandom::globalPosition
-(
-    const label& start,
-    const label& end
-)
-{
-    label value = labelMin;
-
-    if (Pstream::master())
-    {
-        value = round(scalar01()*(end - start));
-    }
-
-    Pstream::scatter(value);
-
-    return start + value;
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.H b/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.H
deleted file mode 100644
index a42034284a..0000000000
--- a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.H
+++ /dev/null
@@ -1,238 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::cachedRandom
-
-Description
-    Random number generator.
-
-    Pre-computes and caches samples on construction, so that when sample01()
-    is called, the function simply returns the next (pre-computed) sample. On
-    reaching the last sample, the sample sequence is repeated.
-
-    Constructed using a seed and sample count. If the supplied count is
-    negative, no caching is performed, and a new sample is generated on each
-    call to sample01().
-
-    Note: the copy constructor cannot be used if count = -1.
-
-SourceFiles
-    cachedRandomI.H
-    cachedRandom.C
-    cachedRandomTemplates.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef cachedRandom_H
-#define cachedRandom_H
-
-#include "scalarList.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-class cachedRandom;
-
-/*---------------------------------------------------------------------------*\
-                       Class cachedRandom Declaration
-\*---------------------------------------------------------------------------*/
-
-class cachedRandom
-{
-    // Private data
-
-        //- Initial random number seed
-        label seed_;
-
-        //- List of scalar samples
-        scalarList samples_;
-
-        //- Current sample marker
-        label sampleI_;
-
-        //- Indicator, which tells if there is a stored gaussian sample
-        bool hasGaussSample_;
-
-        //- Stored sample value
-        scalar gaussSample_;
-
-
-    // Private Member Functions
-
-        //- Returns the current sample
-        scalar scalar01();
-
-
-public:
-
-    // Constructors
-
-        //- Construct given seed and sample count
-        cachedRandom(const label seed, const label count);
-
-        //- Copy constructor with optional reset of sampleI
-        cachedRandom(const cachedRandom& cr, const bool reset = false);
-
-
-    // Destructor
-    ~cachedRandom();
-
-
-    // Member functions
-
-        // Access
-
-            //- Return const access to the initial random number seed
-            inline label seed() const;
-
-            //- Return const access to the list of samples
-            inline const scalarList& samples() const;
-
-            //- Return the current sample marker
-            inline label sampleI() const;
-
-
-        // Manipulation
-
-            //- Return non-const access to the sample marker
-            inline label& sampleI();
-
-
-        // Evaluation
-
-            // Random numbers
-
-                //- Return a sample whose components lie in the range 0-1
-                template<class Type>
-                Type sample01();
-
-                //- Return a sample whose components are normally distributed
-                //  with zero mean and unity variance N(0, 1)
-                template<class Type>
-                Type GaussNormal();
-
-                //- Return a sample between start and end
-                template<class Type>
-                Type position(const Type& start, const Type& end);
-
-                //- Randomise value in the range 0-1
-                template<class Type>
-                void randomise01(Type& value);
-
-                //- Shuffle the values in the list
-                template<class Type>
-                void shuffle(UList<Type>& values);
-
-
-            // Global random numbers - consistent across all processors
-
-                //- Return a sample whose components lie in the range 0-1
-                template<class Type>
-                Type globalSample01();
-
-                //- Return a sample whose components are normally distributed
-                //  with zero mean and unity variance N(0, 1)
-                template<class Type>
-                Type globalGaussNormal();
-
-                //- Return a sample between start and end
-                template<class Type>
-                Type globalPosition(const Type& start, const Type& end);
-
-                //- Randomise value in the range 0-1
-                template<class Type>
-                void globalRandomise01(Type& value);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-// Template specialisations
-
-template<>
-scalar cachedRandom::sample01<scalar>();
-
-template<>
-label cachedRandom::sample01<label>();
-
-template<>
-scalar cachedRandom::GaussNormal<scalar>();
-
-template<>
-label cachedRandom::GaussNormal<label>();
-
-template<>
-scalar cachedRandom::position<scalar>
-(
-    const scalar& start,
-    const scalar& end
-);
-
-template<>
-label cachedRandom::position<label>(const label& start, const label& end);
-
-template<>
-scalar cachedRandom::globalSample01<scalar>();
-
-template<>
-label cachedRandom::globalSample01<label>();
-
-template<>
-scalar cachedRandom::globalGaussNormal<scalar>();
-
-template<>
-label cachedRandom::globalGaussNormal<label>();
-
-template<>
-scalar cachedRandom::globalPosition<scalar>
-(
-    const scalar& start,
-    const scalar& end
-);
-
-template<>
-label cachedRandom::globalPosition<label>(const label& start, const label& end);
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "cachedRandomI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "cachedRandomTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddy.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddy.C
index b4a73f3b91..1883fac278 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddy.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddy.C
@@ -129,7 +129,7 @@ Foam::eddy::eddy
     const scalar x,
     const scalar sigmaX,
     const symmTensor& R,
-    cachedRandom& rndGen
+    Random& rndGen
 )
 :
     patchFaceI_(patchFaceI),
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddy.H b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddy.H
index 5bd0feb0f1..be527caa2f 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddy.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddy.H
@@ -41,7 +41,7 @@ SourceFiles
 #include "vector.H"
 #include "point.H"
 #include "tensor.H"
-#include "cachedRandom.H"
+#include "Random.H"
 #include "coordinateSystem.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -112,7 +112,7 @@ class eddy
         ) const;
 
         //- Return a number with zero mean and unit variance
-        inline scalar epsi(cachedRandom& rndGen) const;
+        inline scalar epsi(Random& rndGen) const;
 
 
 public:
@@ -133,7 +133,7 @@ public:
             const scalar x,             // distance from reference position
             const scalar sigmaX,        // length scale
             const symmTensor& R,        // Stress tensor
-            cachedRandom& rndGen
+            Random& rndGen
         );
 
         //- Construct copy
@@ -175,7 +175,7 @@ public:
             inline label dir1() const;
 
             //- Return random vector of -1 and 1's
-            inline vector epsilon(cachedRandom& rndGen) const;
+            inline vector epsilon(Random& rndGen) const;
 
 
         // Helper functions
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddyI.H b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddyI.H
index 82c579fcc5..91531bf13d 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddyI.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddyI.H
@@ -29,7 +29,7 @@ using namespace Foam::constant;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-Foam::scalar Foam::eddy::epsi(cachedRandom& rndGen) const
+Foam::scalar Foam::eddy::epsi(Random& rndGen) const
 {
     // Random number with zero mean and unit variance
     return rndGen.sample01<scalar>() > 0.5 ? 1 : -1;
@@ -84,7 +84,7 @@ inline Foam::scalar Foam::eddy::c1() const
 }
 
 
-Foam::vector Foam::eddy::epsilon(cachedRandom& rndGen) const
+Foam::vector Foam::eddy::epsilon(Random& rndGen) const
 {
     return vector(epsi(rndGen), epsi(rndGen), epsi(rndGen));
 }
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
index c68231649a..5449a6f5e7 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
@@ -746,7 +746,7 @@ turbulentDFSEMInletFvPatchVectorField
     nCellPerEddy_(5),
     patchNormal_(vector::zero),
     v0_(0),
-    rndGen_(0, -1),
+    rndGen_(Pstream::myProcNo()),
     sigmax_(size(), 0),
     maxSigmaX_(0),
     nEddy_(0),
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H
index b1d0afe62a..d14d82266f 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H
@@ -81,7 +81,7 @@ SourceFiles
 #define turbulentDFSEMInletFvPatchVectorField_H
 
 #include "fixedValueFvPatchFields.H"
-#include "cachedRandom.H"
+#include "Random.H"
 #include "eddy.H"
 #include "pointIndexHit.H"
 #include "instantList.H"
@@ -184,7 +184,7 @@ class turbulentDFSEMInletFvPatchVectorField
         scalar v0_;
 
         //- Random number generator
-        cachedRandom rndGen_;
+        Random rndGen_;
 
         //- Length scale per patch face
         scalarField sigmax_;
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentInlet/turbulentInletFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentInlet/turbulentInletFvPatchField.C
index ec71541bd0..b457ecd85a 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentInlet/turbulentInletFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentInlet/turbulentInletFvPatchField.C
@@ -166,7 +166,7 @@ void Foam::turbulentInletFvPatchField<Type>::updateCoeffs()
 
         forAll(patchField, facei)
         {
-            ranGen_.randomise(randomField[facei]);
+            ranGen_.randomise01<Type>(randomField[facei]);
         }
 
         // Correction-factor to compensate for the loss of RMS fluctuation
diff --git a/src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.H b/src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.H
index effb0904fc..580babf826 100644
--- a/src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.H
+++ b/src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.H
@@ -71,7 +71,6 @@ SourceFiles
 #include "surfaceFieldsFwd.H"
 #include "vectorList.H"
 #include "globalIndex.H"
-#include "cachedRandom.H"
 #include "eulerianParticle.H"
 #include "IOdictionary.H"
 #include "injectedParticleCloud.H"
diff --git a/src/functionObjects/field/particleDistribution/particleDistribution.C b/src/functionObjects/field/particleDistribution/particleDistribution.C
index 71919b0933..b8cf35122b 100644
--- a/src/functionObjects/field/particleDistribution/particleDistribution.C
+++ b/src/functionObjects/field/particleDistribution/particleDistribution.C
@@ -61,7 +61,7 @@ Foam::functionObjects::particleDistribution::particleDistribution
     cloudName_("unknown-cloudName"),
     nameVsBinWidth_(),
     tagFieldName_("none"),
-    rndGen_(1234, -1),
+    rndGen_(),
     writerPtr_(nullptr)
 {
     read(dict);
diff --git a/src/functionObjects/field/particleDistribution/particleDistribution.H b/src/functionObjects/field/particleDistribution/particleDistribution.H
index fec579f193..4e7765e1b1 100644
--- a/src/functionObjects/field/particleDistribution/particleDistribution.H
+++ b/src/functionObjects/field/particleDistribution/particleDistribution.H
@@ -73,7 +73,7 @@ SourceFiles
 #include "fvMeshFunctionObject.H"
 #include "writeFile.H"
 #include "scalarField.H"
-#include "cachedRandom.H"
+#include "Random.H"
 #include "Tuple2.H"
 #include "writer.H"
 
@@ -107,7 +107,7 @@ protected:
         word tagFieldName_;
 
         //- Random number generator - used by distribution models
-        cachedRandom rndGen_;
+        Random rndGen_;
 
         //- Writer
         autoPtr<writer<scalar>> writerPtr_;
diff --git a/src/functionObjects/field/randomise/randomiseTemplates.C b/src/functionObjects/field/randomise/randomiseTemplates.C
index a466ef3551..f0c691eb39 100644
--- a/src/functionObjects/field/randomise/randomiseTemplates.C
+++ b/src/functionObjects/field/randomise/randomiseTemplates.C
@@ -47,7 +47,7 @@ bool Foam::functionObjects::randomise::calcRandomised()
         forAll(field, celli)
         {
             Type rndPert;
-            rand.randomise(rndPert);
+            rand.randomise01(rndPert);
             rndPert = 2.0*rndPert - pTraits<Type>::one;
             rndPert /= mag(rndPert);
             rfield[celli] += magPerturbation_*rndPert;
diff --git a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
index 6da1c1aa8e..fadaf872bc 100644
--- a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
+++ b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
@@ -152,7 +152,7 @@ void Foam::DSMCCloud<ParcelType>::initialise
                 if
                 (
                     (particlesRequired - nParticlesToInsert)
-                  > rndGen_.scalar01()
+                  > rndGen_.sample01<scalar>()
                 )
                 {
                     nParticlesToInsert++;
@@ -282,7 +282,7 @@ void Foam::DSMCCloud<ParcelType>::collisions()
                 // subCell candidate selection procedure
 
                 // Select the first collision candidate
-                label candidateP = rndGen_.integer(0, nC - 1);
+                label candidateP = rndGen_.position<label>(0, nC - 1);
 
                 // Declare the second collision candidate
                 label candidateQ = -1;
@@ -298,7 +298,8 @@ void Foam::DSMCCloud<ParcelType>::collisions()
 
                     do
                     {
-                        candidateQ = subCellPs[rndGen_.integer(0, nSC - 1)];
+                        label i = rndGen_.position<label>(0, nSC - 1);
+                        candidateQ = subCellPs[i];
                     } while (candidateP == candidateQ);
                 }
                 else
@@ -309,7 +310,7 @@ void Foam::DSMCCloud<ParcelType>::collisions()
 
                     do
                     {
-                        candidateQ = rndGen_.integer(0, nC - 1);
+                        candidateQ = rndGen_.position<label>(0, nC - 1);
                     } while (candidateP == candidateQ);
                 }
 
@@ -317,15 +318,15 @@ void Foam::DSMCCloud<ParcelType>::collisions()
                 // uniform candidate selection procedure
 
                 // // Select the first collision candidate
-                // label candidateP = rndGen_.integer(0, nC-1);
+                // label candidateP = rndGen_.position<label>(0, nC-1);
 
                 // // Select a possible second collision candidate
-                // label candidateQ = rndGen_.integer(0, nC-1);
+                // label candidateQ = rndGen_.position<label>(0, nC-1);
 
                 // // If the same candidate is chosen, choose again
                 // while (candidateP == candidateQ)
                 // {
-                //     candidateQ = rndGen_.integer(0, nC-1);
+                //     candidateQ = rndGen_.position<label>(0, nC-1);
                 // }
 
                 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -348,7 +349,7 @@ void Foam::DSMCCloud<ParcelType>::collisions()
                     sigmaTcRMax_[celli] = sigmaTcR;
                 }
 
-                if ((sigmaTcR/sigmaTcRMax) > rndGen_.scalar01())
+                if ((sigmaTcR/sigmaTcRMax) > rndGen_.sample01<scalar>())
                 {
                     binaryCollision().collide
                     (
@@ -648,7 +649,7 @@ Foam::DSMCCloud<ParcelType>::DSMCCloud
         mesh_
     ),
     constProps_(),
-    rndGen_(label(149382906) + 7183*Pstream::myProcNo()),
+    rndGen_(Pstream::myProcNo()),
     boundaryT_
     (
         volScalarField
@@ -711,7 +712,7 @@ Foam::DSMCCloud<ParcelType>::DSMCCloud
     // and 1.
     forAll(collisionSelectionRemainder_, i)
     {
-        collisionSelectionRemainder_[i] = rndGen_.scalar01();
+        collisionSelectionRemainder_[i] = rndGen_.sample01<scalar>();
     }
 
     if (readFields)
@@ -900,7 +901,7 @@ Foam::DSMCCloud<ParcelType>::DSMCCloud
         )
     ),
     constProps_(),
-    rndGen_(label(971501) + 1526*Pstream::myProcNo()),
+    rndGen_(Pstream::myProcNo()),
     boundaryT_
     (
         volScalarField
@@ -1039,12 +1040,7 @@ Foam::vector Foam::DSMCCloud<ParcelType>::equipartitionLinearVelocity
 {
     return
         sqrt(physicoChemical::k.value()*temperature/mass)
-       *vector
-        (
-            rndGen_.GaussNormal(),
-            rndGen_.GaussNormal(),
-            rndGen_.GaussNormal()
-        );
+       *rndGen_.GaussNormal<vector>();
 }
 
 
@@ -1064,7 +1060,9 @@ Foam::scalar Foam::DSMCCloud<ParcelType>::equipartitionInternalEnergy
     else if (iDof < 2.0 + SMALL && iDof > 2.0 - SMALL)
     {
         // Special case for iDof = 2, i.e. diatomics;
-        Ei = -log(rndGen_.scalar01())*physicoChemical::k.value()*temperature;
+        Ei =
+           -log(rndGen_.sample01<scalar>())
+           *physicoChemical::k.value()*temperature;
     }
     else
     {
@@ -1074,9 +1072,9 @@ Foam::scalar Foam::DSMCCloud<ParcelType>::equipartitionInternalEnergy
 
         do
         {
-            energyRatio = 10*rndGen_.scalar01();
+            energyRatio = 10*rndGen_.sample01<scalar>();
             P = pow((energyRatio/a), a)*exp(a - energyRatio);
-        } while (P < rndGen_.scalar01());
+        } while (P < rndGen_.sample01<scalar>());
 
         Ei = energyRatio*physicoChemical::k.value()*temperature;
     }
diff --git a/src/lagrangian/DSMC/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C b/src/lagrangian/DSMC/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C
index 9139e09be3..d81ce400f6 100644
--- a/src/lagrangian/DSMC/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C
+++ b/src/lagrangian/DSMC/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C
@@ -39,14 +39,14 @@ Foam::scalar Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::energyRatio
 {
     CloudType& cloud(this->owner());
 
-    Random& rndGen(cloud.rndGen());
+    Random& rndGen = cloud.rndGen();
 
     scalar ChiAMinusOne = ChiA - 1;
     scalar ChiBMinusOne = ChiB - 1;
 
     if (ChiAMinusOne < SMALL && ChiBMinusOne < SMALL)
     {
-        return rndGen.scalar01();
+        return rndGen.sample01<scalar>();
     }
 
     scalar energyRatio;
@@ -56,7 +56,7 @@ Foam::scalar Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::energyRatio
     {
         P = 0;
 
-        energyRatio = rndGen.scalar01();
+        energyRatio = rndGen.sample01<scalar>();
 
         if (ChiAMinusOne < SMALL)
         {
@@ -81,7 +81,7 @@ Foam::scalar Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::energyRatio
                     ChiBMinusOne
                 );
         }
-    } while (P < rndGen.scalar01());
+    } while (P < rndGen.sample01<scalar>());
 
     return energyRatio;
 }
@@ -186,7 +186,7 @@ void Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::collide
     scalar& EiP = pP.Ei();
     scalar& EiQ = pQ.Ei();
 
-    Random& rndGen(cloud.rndGen());
+    Random& rndGen = cloud.rndGen();
 
     scalar inverseCollisionNumber = 1/relaxationCollisionNumber_;
 
@@ -217,13 +217,14 @@ void Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::collide
 
     if (iDofP > 0)
     {
-        if (inverseCollisionNumber > rndGen.scalar01())
+        if (inverseCollisionNumber > rndGen.sample01<scalar>())
         {
             availableEnergy += preCollisionEiP;
 
             if (iDofP == 2)
             {
-                scalar energyRatio = 1.0 - pow(rndGen.scalar01(), (1.0/ChiB));
+                scalar energyRatio =
+                    1.0 - pow(rndGen.sample01<scalar>(), (1.0/ChiB));
                 EiP = energyRatio*availableEnergy;
             }
             else
@@ -238,13 +239,13 @@ void Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::collide
 
     if (iDofQ > 0)
     {
-        if (inverseCollisionNumber > rndGen.scalar01())
+        if (inverseCollisionNumber > rndGen.sample01<scalar>())
         {
             availableEnergy += preCollisionEiQ;
 
             if (iDofQ == 2)
             {
-                scalar energyRatio = 1.0 - pow(rndGen.scalar01(), (1.0/ChiB));
+                scalar energyRatio = 1.0 - pow(rndGen.sample01<scalar>(), (1.0/ChiB));
                 EiQ = energyRatio*availableEnergy;
             }
             else
@@ -261,9 +262,9 @@ void Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::collide
     scalar cR = sqrt(2.0*availableEnergy/mR);
 
     // Variable Hard Sphere collision part
-    scalar cosTheta = 2.0*rndGen.scalar01() - 1.0;
+    scalar cosTheta = 2.0*rndGen.sample01<scalar>() - 1.0;
     scalar sinTheta = sqrt(1.0 - cosTheta*cosTheta);
-    scalar phi = twoPi*rndGen.scalar01();
+    scalar phi = twoPi*rndGen.sample01<scalar>();
 
     vector postCollisionRelU =
         cR
diff --git a/src/lagrangian/DSMC/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C b/src/lagrangian/DSMC/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C
index 19d1ed61d8..3f20bd3bf9 100644
--- a/src/lagrangian/DSMC/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C
+++ b/src/lagrangian/DSMC/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C
@@ -121,7 +121,7 @@ void Foam::VariableHardSphere<CloudType>::collide
     vector& UP = pP.U();
     vector& UQ = pQ.U();
 
-    Random& rndGen(cloud.rndGen());
+    Random& rndGen = cloud.rndGen();
 
     scalar mP = cloud.constProps(typeIdP).mass();
 
@@ -131,11 +131,11 @@ void Foam::VariableHardSphere<CloudType>::collide
 
     scalar cR = mag(UP - UQ);
 
-    scalar cosTheta = 2.0*rndGen.scalar01() - 1.0;
+    scalar cosTheta = 2.0*rndGen.sample01<scalar>() - 1.0;
 
     scalar sinTheta = sqrt(1.0 - cosTheta*cosTheta);
 
-    scalar phi = twoPi*rndGen.scalar01();
+    scalar phi = twoPi*rndGen.sample01<scalar>();
 
     vector postCollisionRelU =
         cR
diff --git a/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
index 314bd25977..4ed53b5cf6 100644
--- a/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
+++ b/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
@@ -147,7 +147,7 @@ void Foam::FreeStream<CloudType>::inflow()
 
     const scalar deltaT = mesh.time().deltaTValue();
 
-    Random& rndGen(cloud.rndGen());
+    Random& rndGen = cloud.rndGen();
 
     scalar sqrtPi = sqrt(pi);
 
@@ -291,7 +291,7 @@ void Foam::FreeStream<CloudType>::inflow()
 
                 // Add another particle with a probability proportional to the
                 // remainder of taking the integer part of faceAccumulator
-                if ((faceAccumulator - nI) > rndGen.scalar01())
+                if ((faceAccumulator - nI) > rndGen.sample01<scalar>())
                 {
                     nI++;
                 }
@@ -307,7 +307,7 @@ void Foam::FreeStream<CloudType>::inflow()
                     // Choose a triangle to insert on, based on their relative
                     // area
 
-                    scalar triSelection = rndGen.scalar01();
+                    scalar triSelection = rndGen.sample01<scalar>();
 
                     // Selected triangle
                     label selectedTriI = -1;
@@ -371,7 +371,7 @@ void Foam::FreeStream<CloudType>::inflow()
                     do
                     {
                         uNormalThermal =
-                            randomScaling*(2.0*rndGen.scalar01() - 1);
+                            randomScaling*(2.0*rndGen.sample01<scalar>() - 1);
 
                         uNormal = uNormalThermal + sCosTheta;
 
@@ -385,13 +385,13 @@ void Foam::FreeStream<CloudType>::inflow()
                                *exp(uNormProbCoeffB - sqr(uNormalThermal));
                         }
 
-                    } while (P < rndGen.scalar01());
+                    } while (P < rndGen.sample01<scalar>());
 
                     vector U =
                         sqrt(physicoChemical::k.value()*faceTemperature/mass)
                        *(
-                            rndGen.GaussNormal()*t1
-                          + rndGen.GaussNormal()*t2
+                            rndGen.GaussNormal<scalar>()*t1
+                          + rndGen.GaussNormal<scalar>()*t2
                         )
                       + (t1 & faceVelocity)*t1
                       + (t2 & faceVelocity)*t2
diff --git a/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C b/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
index 3bc25401ed..9f4ff92bde 100644
--- a/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
+++ b/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
@@ -78,19 +78,18 @@ void Foam::MaxwellianThermal<CloudType>::correct
 
     CloudType& cloud(this->owner());
 
-    Random& rndGen(cloud.rndGen());
+    Random& rndGen = cloud.rndGen();
 
     while (mag(Ut) < SMALL)
     {
         // If the incident velocity is parallel to the face normal, no
         // tangential direction can be chosen.  Add a perturbation to the
         // incoming velocity and recalculate.
-
         U = vector
         (
-            U.x()*(0.8 + 0.2*rndGen.scalar01()),
-            U.y()*(0.8 + 0.2*rndGen.scalar01()),
-            U.z()*(0.8 + 0.2*rndGen.scalar01())
+            U.x()*(0.8 + 0.2*rndGen.sample01<scalar>()),
+            U.y()*(0.8 + 0.2*rndGen.sample01<scalar>()),
+            U.z()*(0.8 + 0.2*rndGen.sample01<scalar>())
         );
 
         U_dot_nw = U & nw;
@@ -113,9 +112,9 @@ void Foam::MaxwellianThermal<CloudType>::correct
     U =
         sqrt(physicoChemical::k.value()*T/mass)
        *(
-            rndGen.GaussNormal()*tw1
-          + rndGen.GaussNormal()*tw2
-          - sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw
+            rndGen.GaussNormal<scalar>()*tw1
+          + rndGen.GaussNormal<scalar>()*tw2
+          - sqrt(-2.0*log(max(1 - rndGen.sample01<scalar>(), VSMALL)))*nw
         );
 
     U += cloud.boundaryU().boundaryField()[wppIndex][wppLocalFace];
diff --git a/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C b/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C
index 06e5a7079c..056251c533 100644
--- a/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C
+++ b/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C
@@ -73,9 +73,9 @@ void Foam::MixedDiffuseSpecular<CloudType>::correct
 
     CloudType& cloud(this->owner());
 
-    Random& rndGen(cloud.rndGen());
+    Random& rndGen = cloud.rndGen();
 
-    if (diffuseFraction_ > rndGen.scalar01())
+    if (diffuseFraction_ > rndGen.sample01<scalar>())
     {
         // Diffuse reflection
 
@@ -87,12 +87,11 @@ void Foam::MixedDiffuseSpecular<CloudType>::correct
             // If the incident velocity is parallel to the face normal, no
             // tangential direction can be chosen.  Add a perturbation to the
             // incoming velocity and recalculate.
-
             U = vector
             (
-                U.x()*(0.8 + 0.2*rndGen.scalar01()),
-                U.y()*(0.8 + 0.2*rndGen.scalar01()),
-                U.z()*(0.8 + 0.2*rndGen.scalar01())
+                U.x()*(0.8 + 0.2*rndGen.sample01<scalar>()),
+                U.y()*(0.8 + 0.2*rndGen.sample01<scalar>()),
+                U.z()*(0.8 + 0.2*rndGen.sample01<scalar>())
             );
 
             U_dot_nw = U & nw;
@@ -115,9 +114,9 @@ void Foam::MixedDiffuseSpecular<CloudType>::correct
         U =
             sqrt(physicoChemical::k.value()*T/mass)
            *(
-                rndGen.GaussNormal()*tw1
-              + rndGen.GaussNormal()*tw2
-              - sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw
+                rndGen.GaussNormal<scalar>()*tw1
+              + rndGen.GaussNormal<scalar>()*tw2
+              - sqrt(-2.0*log(max(1 - rndGen.sample01<scalar>(), VSMALL)))*nw
             );
 
         U += cloud.boundaryU().boundaryField()[wppIndex][wppLocalFace];
diff --git a/src/lagrangian/distributionModels/RosinRammler/RosinRammler.C b/src/lagrangian/distributionModels/RosinRammler/RosinRammler.C
index 2e5a045619..0ed45f62a2 100644
--- a/src/lagrangian/distributionModels/RosinRammler/RosinRammler.C
+++ b/src/lagrangian/distributionModels/RosinRammler/RosinRammler.C
@@ -42,7 +42,7 @@ namespace distributionModels
 Foam::distributionModels::RosinRammler::RosinRammler
 (
     const dictionary& dict,
-    cachedRandom& rndGen
+    Random& rndGen
 )
 :
     distributionModel(typeName, dict, rndGen),
diff --git a/src/lagrangian/distributionModels/RosinRammler/RosinRammler.H b/src/lagrangian/distributionModels/RosinRammler/RosinRammler.H
index 448165a703..f7f7183d38 100644
--- a/src/lagrangian/distributionModels/RosinRammler/RosinRammler.H
+++ b/src/lagrangian/distributionModels/RosinRammler/RosinRammler.H
@@ -82,7 +82,7 @@ public:
     // Constructors
 
         //- Construct from components
-        RosinRammler(const dictionary& dict, cachedRandom& rndGen);
+        RosinRammler(const dictionary& dict, Random& rndGen);
 
         //- Construct copy
         RosinRammler(const RosinRammler& p);
diff --git a/src/lagrangian/distributionModels/binned/binned.C b/src/lagrangian/distributionModels/binned/binned.C
index b9ed2813f7..57530fa649 100644
--- a/src/lagrangian/distributionModels/binned/binned.C
+++ b/src/lagrangian/distributionModels/binned/binned.C
@@ -78,7 +78,7 @@ void Foam::distributionModels::binned::initialise()
 Foam::distributionModels::binned::binned
 (
     const dictionary& dict,
-    cachedRandom& rndGen
+    Random& rndGen
 )
 :
     distributionModel(typeName, dict, rndGen),
@@ -102,7 +102,7 @@ Foam::distributionModels::binned::binned
 (
     const UList<scalar>& sampleData,
     const scalar binWidth,
-    cachedRandom& rndGen
+    Random& rndGen
 )
 :
     distributionModel(typeName, dictionary::null, rndGen),
diff --git a/src/lagrangian/distributionModels/binned/binned.H b/src/lagrangian/distributionModels/binned/binned.H
index a579a4a40b..4e719a1ef0 100644
--- a/src/lagrangian/distributionModels/binned/binned.H
+++ b/src/lagrangian/distributionModels/binned/binned.H
@@ -96,14 +96,14 @@ public:
     // Constructors
 
         //- Construct from dictionary
-        binned(const dictionary& dict, cachedRandom& rndGen);
+        binned(const dictionary& dict, Random& rndGen);
 
         //- Construct from components
         binned
         (
             const UList<scalar>& sampleData,
             const scalar binWidth,
-            cachedRandom& rndGen
+            Random& rndGen
         );
 
         //- Construct copy
diff --git a/src/lagrangian/distributionModels/distributionModel/distributionModel.C b/src/lagrangian/distributionModels/distributionModel/distributionModel.C
index 20a5cf82a4..6343ea3757 100644
--- a/src/lagrangian/distributionModels/distributionModel/distributionModel.C
+++ b/src/lagrangian/distributionModels/distributionModel/distributionModel.C
@@ -66,7 +66,7 @@ Foam::distributionModels::distributionModel::distributionModel
 (
     const word& name,
     const dictionary& dict,
-    cachedRandom& rndGen
+    Random& rndGen
 )
 :
     distributionModelDict_(dict),
diff --git a/src/lagrangian/distributionModels/distributionModel/distributionModel.H b/src/lagrangian/distributionModels/distributionModel/distributionModel.H
index 5e712e10ec..0be1489929 100644
--- a/src/lagrangian/distributionModels/distributionModel/distributionModel.H
+++ b/src/lagrangian/distributionModels/distributionModel/distributionModel.H
@@ -55,7 +55,7 @@ SourceFiles
 
 #include "IOdictionary.H"
 #include "autoPtr.H"
-#include "cachedRandom.H"
+#include "Random.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -79,7 +79,7 @@ protected:
         const dictionary distributionModelDict_;
 
         //- Reference to the random number generator
-        cachedRandom& rndGen_;
+        Random& rndGen_;
 
 
     // Protected Member Functions
@@ -102,7 +102,7 @@ public:
         dictionary,
         (
             const dictionary& dict,
-            cachedRandom& rndGen
+            Random& rndGen
         ),
         (dict, rndGen)
     );
@@ -115,7 +115,7 @@ public:
         (
             const word& name,
             const dictionary& dict,
-            cachedRandom& rndGen
+            Random& rndGen
         );
 
         //- Construct copy
@@ -129,7 +129,7 @@ public:
     static autoPtr<distributionModel> New
     (
         const dictionary& dict,
-        cachedRandom& rndGen
+        Random& rndGen
     );
 
 
diff --git a/src/lagrangian/distributionModels/distributionModel/distributionModelNew.C b/src/lagrangian/distributionModels/distributionModel/distributionModelNew.C
index 3e3f4ba249..e458085a43 100644
--- a/src/lagrangian/distributionModels/distributionModel/distributionModelNew.C
+++ b/src/lagrangian/distributionModels/distributionModel/distributionModelNew.C
@@ -31,7 +31,7 @@ Foam::autoPtr<Foam::distributionModels::distributionModel>
 Foam::distributionModels::distributionModel::New
 (
     const dictionary& dict,
-    cachedRandom& rndGen
+    Random& rndGen
 )
 {
     const word modelType(dict.lookup("type"));
diff --git a/src/lagrangian/distributionModels/exponential/exponential.C b/src/lagrangian/distributionModels/exponential/exponential.C
index 013e5c39bb..3253151351 100644
--- a/src/lagrangian/distributionModels/exponential/exponential.C
+++ b/src/lagrangian/distributionModels/exponential/exponential.C
@@ -42,7 +42,7 @@ namespace distributionModels
 Foam::distributionModels::exponential::exponential
 (
     const dictionary& dict,
-    cachedRandom& rndGen
+    Random& rndGen
 )
 :
     distributionModel(typeName, dict, rndGen),
diff --git a/src/lagrangian/distributionModels/exponential/exponential.H b/src/lagrangian/distributionModels/exponential/exponential.H
index 76f62acbb8..a68bc4f804 100644
--- a/src/lagrangian/distributionModels/exponential/exponential.H
+++ b/src/lagrangian/distributionModels/exponential/exponential.H
@@ -75,7 +75,7 @@ public:
     // Constructors
 
         //- Construct from components
-        exponential(const dictionary& dict, cachedRandom& rndGen);
+        exponential(const dictionary& dict, Random& rndGen);
 
         //- Construct copy
         exponential(const exponential& p);
diff --git a/src/lagrangian/distributionModels/fixedValue/fixedValue.C b/src/lagrangian/distributionModels/fixedValue/fixedValue.C
index f8f0c234a9..0275a404ad 100644
--- a/src/lagrangian/distributionModels/fixedValue/fixedValue.C
+++ b/src/lagrangian/distributionModels/fixedValue/fixedValue.C
@@ -42,7 +42,7 @@ namespace distributionModels
 Foam::distributionModels::fixedValue::fixedValue
 (
     const dictionary& dict,
-    cachedRandom& rndGen
+    Random& rndGen
 )
 :
     distributionModel(typeName, dict, rndGen),
diff --git a/src/lagrangian/distributionModels/fixedValue/fixedValue.H b/src/lagrangian/distributionModels/fixedValue/fixedValue.H
index aa08973031..875d172e42 100644
--- a/src/lagrangian/distributionModels/fixedValue/fixedValue.H
+++ b/src/lagrangian/distributionModels/fixedValue/fixedValue.H
@@ -66,7 +66,7 @@ public:
     // Constructors
 
         //- Construct from components
-        fixedValue(const dictionary& dict, cachedRandom& rndGen);
+        fixedValue(const dictionary& dict, Random& rndGen);
 
         //- Construct copy
         fixedValue(const fixedValue& p);
diff --git a/src/lagrangian/distributionModels/general/general.C b/src/lagrangian/distributionModels/general/general.C
index 87339642bf..9f96ff7b80 100644
--- a/src/lagrangian/distributionModels/general/general.C
+++ b/src/lagrangian/distributionModels/general/general.C
@@ -78,7 +78,7 @@ void Foam::distributionModels::general::initialise()
 Foam::distributionModels::general::general
 (
     const dictionary& dict,
-    cachedRandom& rndGen
+    Random& rndGen
 )
 :
     distributionModel(typeName, dict, rndGen),
@@ -96,7 +96,7 @@ Foam::distributionModels::general::general
 (
     const UList<scalar>& sampleData,
     const scalar binWidth,
-    cachedRandom& rndGen
+    Random& rndGen
 )
 :
     distributionModel(typeName, dictionary::null, rndGen),
diff --git a/src/lagrangian/distributionModels/general/general.H b/src/lagrangian/distributionModels/general/general.H
index db8a18194a..689cfa6cc7 100644
--- a/src/lagrangian/distributionModels/general/general.H
+++ b/src/lagrangian/distributionModels/general/general.H
@@ -96,14 +96,14 @@ public:
     // Constructors
 
         //- Construct from components
-        general(const dictionary& dict, cachedRandom& rndGen);
+        general(const dictionary& dict, Random& rndGen);
 
         //- Construct from components
         general
         (
             const UList<scalar>& sampleData,
             const scalar binWidth,
-            cachedRandom& rndGen
+            Random& rndGen
         );
 
         //- Construct copy
diff --git a/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.C b/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.C
index 015343b711..23339eca0f 100644
--- a/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.C
+++ b/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.C
@@ -42,7 +42,7 @@ namespace distributionModels
 Foam::distributionModels::massRosinRammler::massRosinRammler
 (
     const dictionary& dict,
-    cachedRandom& rndGen
+    Random& rndGen
 )
 :
     distributionModel(typeName, dict, rndGen),
diff --git a/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.H b/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.H
index f1745326b6..1f180f5bd6 100644
--- a/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.H
+++ b/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.H
@@ -93,7 +93,7 @@ public:
     // Constructors
 
         //- Construct from components
-        massRosinRammler(const dictionary& dict, cachedRandom& rndGen);
+        massRosinRammler(const dictionary& dict, Random& rndGen);
 
         //- Construct copy
         massRosinRammler(const massRosinRammler& p);
diff --git a/src/lagrangian/distributionModels/multiNormal/multiNormal.C b/src/lagrangian/distributionModels/multiNormal/multiNormal.C
index 1cbb5a5856..8779281fc0 100644
--- a/src/lagrangian/distributionModels/multiNormal/multiNormal.C
+++ b/src/lagrangian/distributionModels/multiNormal/multiNormal.C
@@ -42,7 +42,7 @@ namespace distributionModels
 Foam::distributionModels::multiNormal::multiNormal
 (
     const dictionary& dict,
-    cachedRandom& rndGen
+    Random& rndGen
 )
 :
     distributionModel(typeName, dict, rndGen),
diff --git a/src/lagrangian/distributionModels/multiNormal/multiNormal.H b/src/lagrangian/distributionModels/multiNormal/multiNormal.H
index 03bb0e270d..a58a96bba1 100644
--- a/src/lagrangian/distributionModels/multiNormal/multiNormal.H
+++ b/src/lagrangian/distributionModels/multiNormal/multiNormal.H
@@ -85,7 +85,7 @@ public:
     // Constructors
 
         //- Construct from components
-        multiNormal(const dictionary& dict, cachedRandom& rndGen);
+        multiNormal(const dictionary& dict, Random& rndGen);
 
         //- Construct copy
         multiNormal(const multiNormal& p);
diff --git a/src/lagrangian/distributionModels/normal/normal.C b/src/lagrangian/distributionModels/normal/normal.C
index 1fc16376a9..a6ebd9d594 100644
--- a/src/lagrangian/distributionModels/normal/normal.C
+++ b/src/lagrangian/distributionModels/normal/normal.C
@@ -43,7 +43,7 @@ namespace distributionModels
 Foam::distributionModels::normal::normal
 (
     const dictionary& dict,
-    cachedRandom& rndGen
+    Random& rndGen
 )
 :
     distributionModel(typeName, dict, rndGen),
diff --git a/src/lagrangian/distributionModels/normal/normal.H b/src/lagrangian/distributionModels/normal/normal.H
index d78ec5a86a..77d11e426a 100644
--- a/src/lagrangian/distributionModels/normal/normal.H
+++ b/src/lagrangian/distributionModels/normal/normal.H
@@ -85,7 +85,7 @@ public:
     // Constructors
 
         //- Construct from components
-        normal(const dictionary& dict, cachedRandom& rndGen);
+        normal(const dictionary& dict, Random& rndGen);
 
         //- Construct copy
         normal(const normal& p);
diff --git a/src/lagrangian/distributionModels/uniform/uniform.C b/src/lagrangian/distributionModels/uniform/uniform.C
index 72a0d2fc3d..ccf60544b5 100644
--- a/src/lagrangian/distributionModels/uniform/uniform.C
+++ b/src/lagrangian/distributionModels/uniform/uniform.C
@@ -42,7 +42,7 @@ namespace distributionModels
 Foam::distributionModels::uniform::uniform
 (
     const dictionary& dict,
-    cachedRandom& rndGen
+    Random& rndGen
 )
 :
     distributionModel(typeName, dict, rndGen),
diff --git a/src/lagrangian/distributionModels/uniform/uniform.H b/src/lagrangian/distributionModels/uniform/uniform.H
index 3efaee8219..a47bb8495f 100644
--- a/src/lagrangian/distributionModels/uniform/uniform.H
+++ b/src/lagrangian/distributionModels/uniform/uniform.H
@@ -70,7 +70,7 @@ public:
     // Constructors
 
         //- Construct from components
-        uniform(const dictionary& dict, cachedRandom& rndGen);
+        uniform(const dictionary& dict, Random& rndGen);
 
         //- Construct copy
         uniform(const uniform& p);
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index d82b1bba67..0dfb75094d 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -316,13 +316,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
     (
         particleProperties_.subOrEmptyDict("subModels", solution_.active())
     ),
-    rndGen_
-    (
-        label(0),
-        solution_.steadyState() ?
-        particleProperties_.lookupOrDefault<label>("randomSampleSize", 100000)
-      : -1
-    ),
+    rndGen_(Pstream::myProcNo()),
     cellOccupancyPtr_(),
     cellLengthScale_(cbrt(mesh_.V())),
     rho_(rho),
@@ -515,7 +509,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
     solution_(mesh),
     constProps_(),
     subModelProperties_(dictionary::null),
-    rndGen_(0, 0),
+    rndGen_(),
     cellOccupancyPtr_(nullptr),
     cellLengthScale_(c.cellLengthScale_),
     rho_(c.rho_),
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
index d0c8b57603..15abbe9ac9 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
@@ -59,7 +59,7 @@ SourceFiles
 #include "kinematicCloud.H"
 #include "IOdictionary.H"
 #include "autoPtr.H"
-#include "cachedRandom.H"
+#include "Random.H"
 #include "fvMesh.H"
 #include "volFields.H"
 #include "fvMatrices.H"
@@ -163,7 +163,7 @@ protected:
         const dictionary subModelProperties_;
 
         //- Random number generator - used by some injection routines
-        mutable cachedRandom rndGen_;
+        mutable Random rndGen_;
 
         //- Cell occupancy information for each parcel, (demand driven)
         autoPtr<List<DynamicList<parcelType*>>> cellOccupancyPtr_;
@@ -366,7 +366,7 @@ public:
             // Cloud data
 
                 //- Return reference to the random object
-                inline cachedRandom& rndGen() const;
+                inline Random& rndGen() const;
 
                 //- Return the cell occupancy information for each
                 //  parcel, non-const access, the caller is
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
index ce92455169..f5ac343730 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
@@ -356,7 +356,7 @@ inline Foam::scalar Foam::KinematicCloud<CloudType>::Dmax() const
 
 
 template<class CloudType>
-inline Foam::cachedRandom& Foam::KinematicCloud<CloudType>::rndGen() const
+inline Foam::Random& Foam::KinematicCloud<CloudType>::rndGen() const
 {
     return rndGen_;
 }
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
index c559e3f5df..5cc2c3a767 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
@@ -171,16 +171,16 @@ void Foam::ParticleCollector<CloudType>::initConcentricCircles()
     }
     else
     {
-        // set 4 quadrants for single sector cases
+        // Set 4 quadrants for single sector cases
         nS = 4;
 
         vector tangent = Zero;
         scalar magTangent = 0.0;
 
-        Random rnd(1234);
+        Random& rnd = this->owner().rndGen();
         while (magTangent < SMALL)
         {
-            vector v = rnd.vector01();
+            vector v = rnd.sample01<vector>();
 
             tangent = v - (v & normal_[0])*normal_[0];
             magTangent = mag(tangent);
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/CellZoneInjection/CellZoneInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/CellZoneInjection/CellZoneInjection.C
index ffa8163dff..06ed2d542c 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/CellZoneInjection/CellZoneInjection.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/CellZoneInjection/CellZoneInjection.C
@@ -40,7 +40,7 @@ void Foam::CellZoneInjection<CloudType>::setPositions
     const fvMesh& mesh = this->owner().mesh();
     const scalarField& V = mesh.V();
     const label nCells = cellZoneCells.size();
-    cachedRandom& rnd = this->owner().rndGen();
+    Random& rnd = this->owner().rndGen();
 
     DynamicList<vector> positions(nCells);          // initial size only
     DynamicList<label> injectorCells(nCells);       // initial size only
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjection/ConeInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjection/ConeInjection.C
index 9406a3d9dc..11300c4273 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjection/ConeInjection.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjection/ConeInjection.C
@@ -110,7 +110,7 @@ Foam::ConeInjection<CloudType>::ConeInjection
         vector tangent = Zero;
         scalar magTangent = 0.0;
 
-        cachedRandom& rnd = this->owner().rndGen();
+        Random& rnd = this->owner().rndGen();
         while (magTangent < SMALL)
         {
             vector v = rnd.sample01<vector>();
@@ -260,7 +260,7 @@ void Foam::ConeInjection<CloudType>::setProperties
     typename CloudType::parcelType& parcel
 )
 {
-    cachedRandom& rnd = this->owner().rndGen();
+    Random& rnd = this->owner().rndGen();
 
     // Set particle velocity
     const label i = parcelI % positionAxis_.size();
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.C
index 208bc3511b..b8bcd7464e 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.C
@@ -175,7 +175,7 @@ Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection
 
     setFlowType();
 
-    cachedRandom& rndGen = this->owner().rndGen();
+    Random& rndGen = this->owner().rndGen();
 
     // Normalise direction vector
     direction_ /= mag(direction_);
@@ -321,7 +321,7 @@ void Foam::ConeNozzleInjection<CloudType>::setPositionAndCell
     label& tetPti
 )
 {
-    cachedRandom& rndGen = this->owner().rndGen();
+    Random& rndGen = this->owner().rndGen();
 
     scalar beta = mathematical::twoPi*rndGen.sample01<scalar>();
     normal_ = tanVec1_*cos(beta) + tanVec2_*sin(beta);
@@ -373,7 +373,7 @@ void Foam::ConeNozzleInjection<CloudType>::setProperties
     typename CloudType::parcelType& parcel
 )
 {
-    cachedRandom& rndGen = this->owner().rndGen();
+    Random& rndGen = this->owner().rndGen();
 
     // Set particle velocity
     const scalar deg2Rad = mathematical::pi/180.0;
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InflationInjection/InflationInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InflationInjection/InflationInjection.C
index d3700001f8..41cf67a56d 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InflationInjection/InflationInjection.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InflationInjection/InflationInjection.C
@@ -205,7 +205,7 @@ Foam::label Foam::InflationInjection<CloudType>::parcelsToInject
 
     newParticles_.clear();
 
-    cachedRandom& rnd = this->owner().rndGen();
+    Random& rnd = this->owner().rndGen();
 
     // Diameter factor, when splitting particles into 4, this is the
     // factor that modifies the diameter.
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectedParticleDistributionInjection/InjectedParticleDistributionInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectedParticleDistributionInjection/InjectedParticleDistributionInjection.C
index 3bef483fc3..2ae5e8031b 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectedParticleDistributionInjection/InjectedParticleDistributionInjection.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectedParticleDistributionInjection/InjectedParticleDistributionInjection.C
@@ -148,7 +148,7 @@ void Foam::InjectedParticleDistributionInjection<CloudType>::initialise()
     scalar minTime = GREAT;
 
     // Populate injector properties, filtering out invalid entries
-    cachedRandom& rnd = this->owner().rndGen();
+    Random& rnd = this->owner().rndGen();
     label injectori = 0;
     forAll(injDiameter, i)
     {
@@ -448,7 +448,7 @@ void Foam::InjectedParticleDistributionInjection<CloudType>::setPositionAndCell
     label& tetPtI
 )
 {
-    cachedRandom& rnd = this->owner().rndGen();
+    Random& rnd = this->owner().rndGen();
     currentInjectori_ = rnd.globalPosition<label>(0, position_.size() - 1);
     currentSamplei_ = rnd.globalPosition<label>(0, resampleSize_ - 1);
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/KinematicLookupTableInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/KinematicLookupTableInjection.C
index d081be437f..739599c153 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/KinematicLookupTableInjection.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/KinematicLookupTableInjection.C
@@ -182,7 +182,7 @@ void Foam::KinematicLookupTableInjection<CloudType>::setPositionAndCell
     label injectorI = 0;
     if (randomise_)
     {
-        cachedRandom& rnd = this->owner().rndGen();
+        Random& rnd = this->owner().rndGen();
         injectorI = rnd.position<label>(0, injectorCells_.size() - 1);
     }
     else
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.C
index c1fd0f678e..1244c85dcb 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.C
@@ -162,7 +162,7 @@ Foam::label Foam::PatchFlowRateInjection<CloudType>::parcelsToInject
 
         scalar nParcels = parcelConcentration_*c*flowRate()*dt;
 
-        cachedRandom& rnd = this->owner().rndGen();
+        Random& rnd = this->owner().rndGen();
 
         label nParcelsToInject = floor(nParcels);
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/PatchInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/PatchInjection.C
index 13227e34ef..c0a898c156 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/PatchInjection.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/PatchInjection.C
@@ -121,7 +121,7 @@ Foam::label Foam::PatchInjection<CloudType>::parcelsToInject
     if ((time0 >= 0.0) && (time0 < duration_))
     {
         scalar nParcels = (time1 - time0)*parcelsPerSecond_;
-        cachedRandom& rnd = this->owner().rndGen();
+        Random& rnd = this->owner().rndGen();
         scalar rndPos = rnd.globalPosition(scalar(0), scalar(1));
         label nParcelsToInject = floor(nParcels);
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.C
index 522725fdeb..713c994902 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.C
@@ -26,7 +26,7 @@ License
 #include "patchInjectionBase.H"
 #include "polyMesh.H"
 #include "SubField.H"
-#include "cachedRandom.H"
+#include "Random.H"
 #include "triPointRef.H"
 #include "volFields.H"
 #include "polyMeshTetDecomposition.H"
@@ -151,7 +151,7 @@ void Foam::patchInjectionBase::updateMesh(const polyMesh& mesh)
 void Foam::patchInjectionBase::setPositionAndCell
 (
     const fvMesh& mesh,
-    cachedRandom& rnd,
+    Random& rnd,
     vector& position,
     label& cellOwner,
     label& tetFacei,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.H b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.H
index 2124b3acfc..8b1d837dee 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.H
@@ -54,7 +54,7 @@ namespace Foam
 // Forward class declarations
 class polyMesh;
 class fvMesh;
-class cachedRandom;
+class Random;
 
 /*---------------------------------------------------------------------------*\
                      Class patchInjectionBase Declaration
@@ -118,7 +118,7 @@ public:
         virtual void setPositionAndCell
         (
             const fvMesh& mesh,
-            cachedRandom& rnd,
+            Random& rnd,
             vector& position,
             label& cellOwner,
             label& tetFacei,
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C
index 6aa5354455..068568a1b5 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C
+++ b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C
@@ -71,7 +71,7 @@ Foam::scalar Foam::IsotropyModels::Stochastic<CloudType>::sampleGauss()
     }
     else
     {
-        cachedRandom& rndGen = this->owner().rndGen();
+        Random& rndGen = this->owner().rndGen();
 
         scalar f, m, x, y;
 
@@ -98,7 +98,7 @@ void Foam::IsotropyModels::Stochastic<CloudType>::calculate()
 {
     const fvMesh& mesh = this->owner().mesh();
     const scalar deltaT(this->owner().db().time().deltaTValue());
-    cachedRandom& rndGen = this->owner().rndGen();
+    Random& rndGen = this->owner().rndGen();
 
     const scalar oneBySqrtThree = sqrt(1.0/3.0);
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.C b/src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.C
index 40290f9a85..7b0982e6a6 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.C
@@ -181,7 +181,7 @@ void Foam::ReactingLookupTableInjection<CloudType>::setPositionAndCell
     label injectorI = 0;
     if (randomise_)
     {
-        cachedRandom& rnd = this->owner().rndGen();
+        Random& rnd = this->owner().rndGen();
         injectorI = rnd.position<label>(0, injectorCells_.size() - 1);
     }
     else
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.C b/src/lagrangian/intermediate/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.C
index 7361402d21..79423e612e 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.C
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.C
@@ -187,7 +187,7 @@ void Foam::ReactingMultiphaseLookupTableInjection<CloudType>::setPositionAndCell
     label injectorI = 0;
     if (randomise_)
     {
-        cachedRandom& rnd = this->owner().rndGen();
+        Random& rnd = this->owner().rndGen();
         injectorI = rnd.position<label>(0, injectorCells_.size() - 1);
     }
     else
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.C b/src/lagrangian/intermediate/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.C
index b93166d4d4..b746e3fd3f 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.C
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.C
@@ -182,7 +182,7 @@ void Foam::ThermoLookupTableInjection<CloudType>::setPositionAndCell
     label injectorI = 0;
     if (randomise_)
     {
-        cachedRandom& rnd = this->owner().rndGen();
+        Random& rnd = this->owner().rndGen();
         injectorI = rnd.position<label>(0, injectorCells_.size() - 1);
     }
     else
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H
index 3faf3102e1..0de4186d09 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H
@@ -101,7 +101,7 @@ protected:
         typedef typename CloudType::parcelType parcelType;
 
         //- Reference to the cloud random number generator
-        cachedRandom& rndGen_;
+        Random& rndGen_;
 
         //- Reference to the cloud thermo package
         const SLGThermo& thermo_;
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
index 761dfc8d94..7db8ded122 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
@@ -1041,11 +1041,11 @@ void Foam::moleculeCloud::createMolecule
     {
         pi = equipartitionAngularMomentum(temperature, cP);
 
-        scalar phi(rndGen_.scalar01()*twoPi);
+        scalar phi(rndGen_.sample01<scalar>()*twoPi);
 
-        scalar theta(rndGen_.scalar01()*twoPi);
+        scalar theta(rndGen_.sample01<scalar>()*twoPi);
 
-        scalar psi(rndGen_.scalar01()*twoPi);
+        scalar psi(rndGen_.sample01<scalar>()*twoPi);
 
         Q = tensor
         (
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
index 63f89a69bd..ebc738cab2 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
@@ -301,12 +301,9 @@ inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
     scalar mass
 )
 {
-    return sqrt(physicoChemical::k.value()*temperature/mass)*vector
-    (
-        rndGen_.GaussNormal(),
-        rndGen_.GaussNormal(),
-        rndGen_.GaussNormal()
-    );
+    return
+        sqrt(physicoChemical::k.value()*temperature/mass)
+       *rndGen_.GaussNormal<vector>();
 }
 
 
@@ -323,17 +320,17 @@ inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum
         return sqrtKbT*vector
         (
             0.0,
-            sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
-            sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
+            sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal<scalar>(),
+            sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal<scalar>()
         );
     }
     else
     {
         return sqrtKbT*vector
         (
-            sqrt(cP.momentOfInertia().xx())*rndGen_.GaussNormal(),
-            sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
-            sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
+            sqrt(cP.momentOfInertia().xx())*rndGen_.GaussNormal<scalar>(),
+            sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal<scalar>(),
+            sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal<scalar>()
         );
     }
 }
diff --git a/src/lagrangian/spray/submodels/AtomizationModel/AtomizationModel/AtomizationModel.H b/src/lagrangian/spray/submodels/AtomizationModel/AtomizationModel/AtomizationModel.H
index 216990ea91..1630aafdb2 100644
--- a/src/lagrangian/spray/submodels/AtomizationModel/AtomizationModel/AtomizationModel.H
+++ b/src/lagrangian/spray/submodels/AtomizationModel/AtomizationModel/AtomizationModel.H
@@ -136,7 +136,7 @@ public:
             const vector& injectionPos,
             const scalar pAmbient,
             const scalar chi,
-            cachedRandom& rndGen
+            Random& rndGen
         ) const = 0;
 };
 
diff --git a/src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.C b/src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.C
index 63b5d4ff5c..63670b0845 100644
--- a/src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.C
+++ b/src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.C
@@ -92,7 +92,7 @@ void Foam::BlobsSheetAtomization<CloudType>::update
     const vector& injectionPos,
     const scalar pAmbient,
     const scalar chi,
-    cachedRandom& rndGen
+    Random& rndGen
 ) const
 {
     scalar lBU =
diff --git a/src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.H b/src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.H
index 56175bf095..e95e930c5a 100644
--- a/src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.H
+++ b/src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.H
@@ -122,7 +122,7 @@ public:
             const vector& injectionPos,
             const scalar pAmbient,
             const scalar chi,
-            cachedRandom& rndGen
+            Random& rndGen
         ) const;
 };
 
diff --git a/src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.C b/src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.C
index c21fdfd583..2cc52829da 100644
--- a/src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.C
+++ b/src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.C
@@ -118,7 +118,7 @@ void Foam::LISAAtomization<CloudType>::update
     const vector& injectionPos,
     const scalar pAmbient,
     const scalar chi,
-    cachedRandom& rndGen
+    Random& rndGen
 ) const
 {
     if (volFlowRate < SMALL)
diff --git a/src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.H b/src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.H
index 1c8ed5c897..2a0f1a5ef0 100644
--- a/src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.H
+++ b/src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.H
@@ -142,7 +142,7 @@ public:
             const vector& injectionPos,
             const scalar pAmbient,
             const scalar chi,
-            cachedRandom& rndGen
+            Random& rndGen
         ) const;
 };
 
diff --git a/src/lagrangian/spray/submodels/AtomizationModel/NoAtomization/NoAtomization.C b/src/lagrangian/spray/submodels/AtomizationModel/NoAtomization/NoAtomization.C
index 038e90c5b4..ed5fed5ec6 100644
--- a/src/lagrangian/spray/submodels/AtomizationModel/NoAtomization/NoAtomization.C
+++ b/src/lagrangian/spray/submodels/AtomizationModel/NoAtomization/NoAtomization.C
@@ -95,7 +95,7 @@ void Foam::NoAtomization<CloudType>::update
     const vector& injectionPos,
     const scalar pAmbient,
     const scalar chi,
-    cachedRandom& rndGen
+    Random& rndGen
 ) const
 {}
 
diff --git a/src/lagrangian/spray/submodels/AtomizationModel/NoAtomization/NoAtomization.H b/src/lagrangian/spray/submodels/AtomizationModel/NoAtomization/NoAtomization.H
index 975f511810..24133879f3 100644
--- a/src/lagrangian/spray/submodels/AtomizationModel/NoAtomization/NoAtomization.H
+++ b/src/lagrangian/spray/submodels/AtomizationModel/NoAtomization/NoAtomization.H
@@ -105,7 +105,7 @@ public:
             const vector& injectionPos,
             const scalar pAmbient,
             const scalar chi,
-            cachedRandom& rndGen
+            Random& rndGen
         ) const;
 };
 
diff --git a/src/lagrangian/spray/submodels/BreakupModel/SHF/SHF.C b/src/lagrangian/spray/submodels/BreakupModel/SHF/SHF.C
index 41e3ac307a..642df18a86 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/SHF/SHF.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/SHF/SHF.C
@@ -138,7 +138,7 @@ bool Foam::SHF<CloudType>::update
     scalar& massChild
 )
 {
-    cachedRandom& rndGen = this->owner().rndGen();
+    Random& rndGen = this->owner().rndGen();
 
     bool addChild = false;
 
diff --git a/src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C b/src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C
index 4f3607b6e9..5817c8c441 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C
@@ -110,7 +110,7 @@ bool Foam::TAB<CloudType>::update
     scalar& massChild
 )
 {
-    cachedRandom& rndGen = this->owner().rndGen();
+    Random& rndGen = this->owner().rndGen();
 
     scalar r = 0.5*d;
     scalar r2 = r*r;
diff --git a/src/lagrangian/turbulence/submodels/Kinematic/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C b/src/lagrangian/turbulence/submodels/Kinematic/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C
index ef93f6a3cf..cdba049f9e 100644
--- a/src/lagrangian/turbulence/submodels/Kinematic/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C
+++ b/src/lagrangian/turbulence/submodels/Kinematic/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C
@@ -100,7 +100,7 @@ Foam::vector Foam::GradientDispersionRAS<CloudType>::update
     scalar& tTurb
 )
 {
-    cachedRandom& rnd = this->owner().rndGen();
+    Random& rnd = this->owner().rndGen();
 
     const scalar cps = 0.16432;
 
diff --git a/src/lagrangian/turbulence/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C b/src/lagrangian/turbulence/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C
index e78fcf9ed9..56ee5d5054 100644
--- a/src/lagrangian/turbulence/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C
+++ b/src/lagrangian/turbulence/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C
@@ -71,7 +71,7 @@ Foam::vector Foam::StochasticDispersionRAS<CloudType>::update
     scalar& tTurb
 )
 {
-    cachedRandom& rnd = this->owner().rndGen();
+    Random& rnd = this->owner().rndGen();
 
     const scalar cps = 0.16432;
 
diff --git a/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C b/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C
index a2228fa5e3..86a24ff4c1 100644
--- a/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C
+++ b/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C
@@ -206,7 +206,7 @@ Foam::forceSuSp Foam::BrownianMotionForce<CloudType>::calcCoupled
 
     // To generate a spherical distribution:
 
-    cachedRandom& rnd = this->owner().rndGen();
+    Random& rnd = this->owner().rndGen();
 
     const scalar theta = rnd.sample01<scalar>()*twoPi;
     const scalar u = 2*rnd.sample01<scalar>() - 1;
diff --git a/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H b/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H
index 0c7a5e89f7..2990a0f253 100644
--- a/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H
+++ b/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H
@@ -46,7 +46,7 @@ SourceFiles
 #define BrownianMotionForce_H
 
 #include "ParticleForce.H"
-#include "cachedRandom.H"
+#include "Random.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,7 +65,7 @@ class BrownianMotionForce
     // Private data
 
         //- Reference to the cloud random number generator
-        cachedRandom& rndGen_;
+        Random& rndGen_;
 
         //- Molecular free path length [m]
         const scalar lambda_;
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
index 8d15aebb2c..8799d2f6ae 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
@@ -259,7 +259,7 @@ bool Foam::edgeIntersections::inlinePerturb
             if (perturbStart)
             {
                 // Perturb with something (hopefully) larger than tolerance.
-                scalar t = 4.0*(rndGen.scalar01() - 0.5);
+                scalar t = 4.0*(rndGen.sample01<scalar>() - 0.5);
                 points1[v0] += t*surf1PointTol[e[0]]*n;
 
                 const labelList& pEdges = surf1.pointEdges()[e[0]];
@@ -272,7 +272,7 @@ bool Foam::edgeIntersections::inlinePerturb
             if (perturbEnd)
             {
                 // Perturb with something larger than tolerance.
-                scalar t = 4.0*(rndGen.scalar01() - 0.5);
+                scalar t = 4.0*(rndGen.sample01<scalar>() - 0.5);
                 points1[v1] += t*surf1PointTol[e[1]]*n;
 
                 const labelList& pEdges = surf1.pointEdges()[e[1]];
@@ -319,7 +319,7 @@ bool Foam::edgeIntersections::rotatePerturb
             //label pointi = e[0];
 
             // Generate random vector slightly larger than tolerance.
-            vector rndVec = rndGen.vector01() - vector(0.5, 0.5, 0.5);
+            vector rndVec = rndGen.sample01<vector>() - vector(0.5, 0.5, 0.5);
 
             // Make sure rndVec only perp to edge
             vector n(points1[meshPoints[e[1]]] - points1[meshPoints[e[0]]]);
@@ -404,7 +404,8 @@ bool Foam::edgeIntersections::offsetPerturb
         if (nearType == triPointRef::POINT || nearType == triPointRef::EDGE)
         {
             // Shift edge towards tri centre
-            vector offset = 0.01*rndGen.scalar01()*(ctr - pHit.hitPoint());
+            vector offset =
+                0.01*rndGen.sample01<scalar>()*(ctr - pHit.hitPoint());
 
             // shift e[0]
             points1[meshPoints[e[0]]] += offset;
diff --git a/src/randomProcesses/processes/UOprocess/UOprocess.C b/src/randomProcesses/processes/UOprocess/UOprocess.C
index 631ae7da52..b15e8ff84d 100644
--- a/src/randomProcesses/processes/UOprocess/UOprocess.C
+++ b/src/randomProcesses/processes/UOprocess/UOprocess.C
@@ -40,9 +40,9 @@ complexVector UOprocess::WeinerProcess()
 {
     return RootDeltaT*complexVector
     (
-        complex(GaussGen.GaussNormal(), GaussGen.GaussNormal()),
-        complex(GaussGen.GaussNormal(), GaussGen.GaussNormal()),
-        complex(GaussGen.GaussNormal(), GaussGen.GaussNormal())
+        complex(GaussGen.GaussNormal<scalar>(), GaussGen.GaussNormal<scalar>()),
+        complex(GaussGen.GaussNormal<scalar>(), GaussGen.GaussNormal<scalar>()),
+        complex(GaussGen.GaussNormal<scalar>(), GaussGen.GaussNormal<scalar>())
     );
 }
 
@@ -57,7 +57,7 @@ UOprocess::UOprocess
     const dictionary& UOdict
 )
 :
-    GaussGen(label(0)),
+    GaussGen(),
     Mesh(kmesh),
     DeltaT(deltaT),
     RootDeltaT(sqrt(DeltaT)),
diff --git a/src/randomProcesses/turbulence/turbGen.C b/src/randomProcesses/turbulence/turbGen.C
index 66162c9a13..b942ee49b6 100644
--- a/src/randomProcesses/turbulence/turbGen.C
+++ b/src/randomProcesses/turbulence/turbGen.C
@@ -51,8 +51,8 @@ Foam::vectorField Foam::turbGen::U()
 
     forAll(K, i)
     {
-        s[i] = RanGen.vector01();
-        rndPhases[i] = RanGen.scalar01();
+        s[i] = RanGen.sample01<vector>();
+        rndPhases[i] = RanGen.sample01<scalar>();
     }
 
     s = K ^ s;
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.C
index 9e6c790540..ef46144e4b 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.C
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.C
@@ -99,7 +99,7 @@ contactAngleForce::contactAngleForce
 :
     force(typeName, owner, dict),
     Ccf_(readScalar(coeffDict_.lookup("Ccf"))),
-    rndGen_(label(0), -1),
+    rndGen_(),
     distribution_
     (
         distributionModels::distributionModel::New
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.H b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.H
index 577c861bbf..44e5efc176 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.H
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.H
@@ -40,7 +40,7 @@ SourceFiles
 
 #include "force.H"
 #include "distributionModel.H"
-#include "cachedRandom.H"
+#include "Random.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -67,7 +67,7 @@ private:
         scalar Ccf_;
 
         //- Random number generator
-        cachedRandom rndGen_;
+        Random rndGen_;
 
         //- Parcel size PDF model
         const autoPtr<distributionModels::distributionModel> distribution_;
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/drippingInjection/drippingInjection.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/drippingInjection/drippingInjection.C
index a429af6bec..3a1e40b494 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/drippingInjection/drippingInjection.C
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/drippingInjection/drippingInjection.C
@@ -57,7 +57,7 @@ drippingInjection::drippingInjection
     injectionModel(type(), owner, dict),
     deltaStable_(readScalar(coeffDict_.lookup("deltaStable"))),
     particlesPerParcel_(readScalar(coeffDict_.lookup("particlesPerParcel"))),
-    rndGen_(label(0), -1),
+    rndGen_(),
     parcelDistribution_
     (
         distributionModels::distributionModel::New
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/drippingInjection/drippingInjection.H b/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/drippingInjection/drippingInjection.H
index e2e2d607e7..277006fcaf 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/drippingInjection/drippingInjection.H
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/drippingInjection/drippingInjection.H
@@ -42,7 +42,7 @@ SourceFiles
 
 #include "injectionModel.H"
 #include "distributionModel.H"
-#include "cachedRandom.H"
+#include "Random.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -84,7 +84,7 @@ protected:
         scalar particlesPerParcel_;
 
         //- Random number generator
-        cachedRandom rndGen_;
+        Random rndGen_;
 
         //- Parcel size PDF model
         const autoPtr<distributionModels::distributionModel>
diff --git a/src/renumber/renumberMethods/randomRenumber/randomRenumber.C b/src/renumber/renumberMethods/randomRenumber/randomRenumber.C
index e5268f7ed7..ecc8d32dce 100644
--- a/src/renumber/renumberMethods/randomRenumber/randomRenumber.C
+++ b/src/renumber/renumberMethods/randomRenumber/randomRenumber.C
@@ -65,7 +65,7 @@ Foam::labelList Foam::randomRenumber::renumber
     {
         forAll(newToOld, i)
         {
-            label j = rndGen.integer(0, newToOld.size()-1);
+            label j = rndGen.position<label>(0, newToOld.size()-1);
             Swap(newToOld[i], newToOld[j]);
         }
     }
diff --git a/src/sampling/sampledSet/patchSeed/patchSeedSet.C b/src/sampling/sampledSet/patchSeed/patchSeedSet.C
index db106e4336..9060f8aeb7 100644
--- a/src/sampling/sampledSet/patchSeed/patchSeedSet.C
+++ b/src/sampling/sampledSet/patchSeed/patchSeedSet.C
@@ -227,7 +227,7 @@ void Foam::patchSeedSet::calcSamples
         {
             forAll(subset, i)
             {
-                label j = rndGen.integer(0, subset.size()-1);
+                label j = rndGen.position<label>(0, subset.size()-1);
                 Swap(subset[i], subset[j]);
             }
         }