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/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/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 d689fa0b1e..a73d7a0e5e 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/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.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/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 31054caa50..0dfb75094d 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -316,7 +316,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
     (
         particleProperties_.subOrEmptyDict("subModels", solution_.active())
     ),
-    rndGen_(),
+    rndGen_(Pstream::myProcNo()),
     cellOccupancyPtr_(),
     cellLengthScale_(cbrt(mesh_.V())),
     rho_(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 559bcfc6e6..c5cd2adf56 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
@@ -260,7 +260,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]];
@@ -273,7 +273,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]];
@@ -320,7 +320,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]]]);
@@ -405,7 +405,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.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.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]);
             }
         }