diff --git a/src/YM/YMfields.jl b/src/YM/YMfields.jl
index c6e5433..d853c49 100644
--- a/src/YM/YMfields.jl
+++ b/src/YM/YMfields.jl
@@ -9,11 +9,11 @@
 ### created: Thu Jul 15 15:16:47 2021
 ###                               
 
-function randomize!(f, lp::SpaceParm, ymws::YMworkspace) 
+function randomize!(f, lp::SpaceParm, ymws::YMworkspace; curng=CUDA.default_rng()) 
         
     if ymws.ALG == SU2alg
         @timeit "Randomize SU(2) algebra field" begin
-            m = CUDA.randn(ymws.PRC, lp.bsz,lp.ndim,3,lp.rsz)
+            m = Random.randn(curng, ymws.PRC, lp.bsz,lp.ndim,3,lp.rsz)
             CUDA.@sync begin
                 CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_assign_SU2!(f,m,lp)
             end
@@ -23,7 +23,7 @@ function randomize!(f, lp::SpaceParm, ymws::YMworkspace)
 
     if ymws.ALG == SU3alg
         @timeit "Randomize SU(3) algebra field" begin
-            m = CUDA.randn(ymws.PRC, lp.bsz,lp.ndim,8,lp.rsz)
+            m = Random.randn(curng, ymws.PRC, lp.bsz,lp.ndim,8,lp.rsz)
             CUDA.@sync begin
                 CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_assign_SU3!(f,m,lp)
             end
diff --git a/src/YM/YMhmc.jl b/src/YM/YMhmc.jl
index c973875..81a93c0 100644
--- a/src/YM/YMhmc.jl
+++ b/src/YM/YMhmc.jl
@@ -58,13 +58,13 @@ function hamiltonian(mom, U, lp, gp, ymws)
     return K+V
 end
 
-function HMC!(U, int::IntrScheme, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}; noacc=false) where T
+function HMC!(U, int::IntrScheme, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}; noacc=false, rng=Random.default_rng(), curng=CUDA.default_rng()) where T
 
     @timeit "HMC trayectory" begin
         
         ymws.U1 .= U
         
-        randomize!(ymws.mom, lp, ymws)
+        randomize!(ymws.mom, lp, ymws; curng)
         hini = hamiltonian(ymws.mom, U, lp, gp, ymws)
         
         MD!(ymws.mom, U, int, lp, gp, ymws)
@@ -78,7 +78,7 @@ function HMC!(U, int::IntrScheme, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspac
         end
         
         if (pacc < 1.0)
-            r = rand()
+            r = rand(rng)
             if (pacc < r) 
                 U .= ymws.U1
                 acc = false
@@ -90,7 +90,7 @@ function HMC!(U, int::IntrScheme, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspac
     end
     return dh, acc
 end
-HMC!(U, eps, ns, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}; noacc=false) where T = HMC!(U, omf4(T, eps, ns), lp, gp, ymws; noacc=noacc)
+HMC!(U, eps, ns, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}; noacc=false, rng=Random.default_rng(), curng=CUDA.default_rng()) where T = HMC!(U, omf4(T, eps, ns), lp, gp, ymws; noacc=noacc, rng, curng)
 
 function MD!(mom, U, int::IntrScheme{NI, T}, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}) where {NI, T <: AbstractFloat}