diff --git a/src/common/distribution_1d_types.F b/src/common/distribution_1d_types.F
index 2be0c1ad78..54851be90d 100644
--- a/src/common/distribution_1d_types.F
+++ b/src/common/distribution_1d_types.F
@@ -25,8 +25,7 @@ MODULE distribution_1d_types
    USE cp_para_env,                     ONLY: cp_para_env_release,&
                                               cp_para_env_retain
    USE cp_para_types,                   ONLY: cp_para_env_type
-   USE parallel_rng_types,              ONLY: delete_rng_stream,&
-                                              rng_stream_p_type
+   USE parallel_rng_types,              ONLY: rng_stream_p_type
 #include "../base/base_uses.f90"
 
    IMPLICIT NONE
@@ -205,8 +204,8 @@ CONTAINS
                      DO iparticle_local = 1, nparticle_local
                         IF (ASSOCIATED(local_particle_set(iparticle_kind)% &
                                        rng(iparticle_local)%stream)) THEN
-                           CALL delete_rng_stream(local_particle_set(iparticle_kind)% &
-                                                  rng(iparticle_local)%stream)
+                           DEALLOCATE (local_particle_set(iparticle_kind)% &
+                                       rng(iparticle_local)%stream)
                         END IF
                      END DO
                      DEALLOCATE (local_particle_set(iparticle_kind)%rng)
diff --git a/src/csvr_system_types.F b/src/csvr_system_types.F
index a1998eac18..d7726689be 100644
--- a/src/csvr_system_types.F
+++ b/src/csvr_system_types.F
@@ -17,8 +17,6 @@ MODULE csvr_system_types
                                               section_vals_val_get
    USE kinds,                           ONLY: dp
    USE parallel_rng_types,              ONLY: GAUSSIAN,&
-                                              create_rng_stream,&
-                                              delete_rng_stream,&
                                               next_rng_seed,&
                                               rng_stream_type
    USE simpar_types,                    ONLY: simpar_type
@@ -39,7 +37,7 @@ MODULE csvr_system_types
       REAL(KIND=dp)                           :: nkt
       REAL(KIND=dp)                           :: thermostat_energy
       REAL(KIND=dp)                           :: region_kin_energy
-      TYPE(rng_stream_type), POINTER          :: gaussian_rng_stream
+      TYPE(rng_stream_type)                   :: gaussian_rng_stream
    END TYPE csvr_thermo_type
 
 ! **************************************************************************************************
@@ -106,7 +104,6 @@ CONTAINS
       ALLOCATE (csvr%nvt(csvr%loc_num_csvr))
       DO i = 1, csvr%loc_num_csvr
          csvr%nvt(i)%thermostat_energy = 0.0_dp
-         NULLIFY (csvr%nvt(i)%gaussian_rng_stream)
       END DO
       ! Initialize the gaussian stream random number
       ALLOCATE (seed(3, 2, csvr%glob_num_csvr))
@@ -123,9 +120,8 @@ CONTAINS
          my_seed = seed(:, :, my_index)
          WRITE (UNIT=name, FMT="(A,I8)") "Wiener process for Thermostat #", my_index
          CALL compress(name)
-         CALL create_rng_stream(rng_stream=csvr%nvt(ithermo)%gaussian_rng_stream, &
-                                name=name, distribution_type=GAUSSIAN, extended_precision=.TRUE., &
-                                seed=my_seed)
+         csvr%nvt(ithermo)%gaussian_rng_stream = rng_stream_type( &
+                                                 name=name, distribution_type=GAUSSIAN, extended_precision=.TRUE., seed=my_seed)
       END DO
       DEALLOCATE (seed)
 
@@ -160,16 +156,8 @@ CONTAINS
       CHARACTER(LEN=*), PARAMETER :: routineN = 'csvr_thermo_dealloc', &
          routineP = moduleN//':'//routineN
 
-      INTEGER                                            :: i
-
-      IF (ASSOCIATED(nvt)) THEN
-         DO i = 1, SIZE(nvt)
-            IF (ASSOCIATED(nvt(i)%gaussian_rng_stream)) THEN
-               CALL delete_rng_stream(nvt(i)%gaussian_rng_stream)
-            ENDIF
-         END DO
+      IF (ASSOCIATED(nvt)) &
          DEALLOCATE (nvt)
-      ENDIF
    END SUBROUTINE csvr_thermo_dealloc
 
 END MODULE csvr_system_types
diff --git a/src/csvr_system_utils.F b/src/csvr_system_utils.F
index 61563bc31c..20c083e52b 100644
--- a/src/csvr_system_utils.F
+++ b/src/csvr_system_utils.F
@@ -11,8 +11,7 @@
 MODULE csvr_system_utils
 
    USE kinds,                           ONLY: dp
-   USE parallel_rng_types,              ONLY: next_random_number,&
-                                              rng_stream_type
+   USE parallel_rng_types,              ONLY: rng_stream_type
 #include "./base/base_uses.f90"
 
    IMPLICIT NONE
@@ -80,7 +79,7 @@ CONTAINS
       REAL(KIND=dp), INTENT(IN)                          :: kk, sigma
       INTEGER, INTENT(IN)                                :: ndeg
       REAL(KIND=dp), INTENT(IN)                          :: taut
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
       REAL(KIND=dp)                                      :: my_res
 
       CHARACTER(len=*), PARAMETER :: routineN = 'rescaling_factor', &
@@ -95,7 +94,7 @@ CONTAINS
          ELSE
             factor = 0.0_dp
          END IF
-         rr = next_random_number(rng_stream)
+         rr = rng_stream%next()
          reverse = 1.0_dp
          ! reverse of momentum is implemented to have the correct limit to Langevin dynamics for ndeg=1
          ! condition: rr < -SQRT(ndeg*kk*factor/(sigma*(1.0_dp-factor)))
@@ -123,7 +122,7 @@ CONTAINS
 ! **************************************************************************************************
    FUNCTION sumnoises(nn, rng_stream) RESULT(sum_gauss)
       INTEGER, INTENT(IN)                                :: nn
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
       REAL(KIND=dp)                                      :: sum_gauss
 
       CHARACTER(len=*), PARAMETER :: routineN = 'sumnoises', routineP = moduleN//':'//routineN
@@ -132,7 +131,7 @@ CONTAINS
 
       sum_gauss = 0.0_dp
       DO i = 1, nn
-         sum_gauss = sum_gauss + next_random_number(rng_stream)**2
+         sum_gauss = sum_gauss + rng_stream%next()**2
       END DO
 
    END FUNCTION sumnoises
diff --git a/src/distribution_methods.F b/src/distribution_methods.F
index 5d989756fe..69311117d0 100644
--- a/src/distribution_methods.F
+++ b/src/distribution_methods.F
@@ -60,9 +60,6 @@ MODULE distribution_methods
                                               molecule_kind_type
    USE molecule_types,                  ONLY: molecule_type
    USE parallel_rng_types,              ONLY: UNIFORM,&
-                                              create_rng_stream,&
-                                              delete_rng_stream,&
-                                              next_random_number,&
                                               rng_stream_type
    USE particle_types,                  ONLY: particle_type
    USE qs_kind_types,                   ONLY: get_qs_kind,&
@@ -1310,25 +1307,23 @@ CONTAINS
          dvec(3), old_var, rn, scaled_cent(3, ncent), var_cl(ncent)
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: dmat
       REAL(KIND=dp), DIMENSION(3, 2)                     :: initial_seed
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type)                              :: rng_stream
 
       CALL timeset(routineN, handle)
 
       initial_seed = REAL(seed, dp); nat = SIZE(coord, 2)
-      NULLIFY (rng_stream)
       ALLOCATE (dmat(ncent, nat))
 
-      CALL create_rng_stream(rng_stream=rng_stream, &
-                             name="kmeans uniform distribution [0,1]", &
-                             distribution_type=UNIFORM, seed=initial_seed)
+      rng_stream = rng_stream_type(name="kmeans uniform distribution [0,1]", &
+                                   distribution_type=UNIFORM, seed=initial_seed)
 
 ! try to find a clever initial guess with centers being somewhat distributed
-      rn = next_random_number(rng_stream)
+      rn = rng_stream%next()
       ind = CEILING(rn*nat)
       cent_coord(:, 1) = coord(:, ind)
       DO i = 2, ncent
          DO
-            rn = next_random_number(rng_stream)
+            rn = rng_stream%next()
             ind = CEILING(rn*nat)
             cent_coord(:, i) = coord(:, ind)
             devi = HUGE(1.0_dp)
@@ -1337,7 +1332,7 @@ CONTAINS
                dist = SQRT(DOT_PRODUCT(dvec, dvec))
                IF (dist .LT. devi) devi = dist
             END DO
-            rn = next_random_number(rng_stream)
+            rn = rng_stream%next()
             IF (rn .LT. devi**2/169.0) EXIT
          END DO
       END DO
@@ -1381,7 +1376,7 @@ CONTAINS
 
             DO i = 1, ncent
                IF (nat_cl(i) == 0) THEN
-                  rn = next_random_number(rng_stream)
+                  rn = rng_stream%next()
                   scaled_cent(:, i) = scaled_coord(:, CEILING(rn*nat))
                ELSE
                   average(:, i, 1) = average(:, i, 1)/REAL(nat_cl(i), dp)
@@ -1395,8 +1390,6 @@ CONTAINS
          END IF
       END DO
 
-      CALL delete_rng_stream(rng_stream)
-
       CALL timestop(handle)
 
    END SUBROUTINE kmeans
diff --git a/src/environment.F b/src/environment.F
index 81a0767779..6d91484e60 100644
--- a/src/environment.F
+++ b/src/environment.F
@@ -82,10 +82,8 @@ MODULE environment
                                               init_spherical_harmonics
    USE parallel_rng_types,              ONLY: GAUSSIAN,&
                                               check_rng,&
-                                              create_rng_stream,&
-                                              init_rng,&
-                                              write_rng_matrices,&
-                                              write_rng_stream
+                                              rng_stream_type,&
+                                              write_rng_matrices
    USE physcon,                         ONLY: write_physcon
    USE reference_manager,               ONLY: collect_citations_from_ranks,&
                                               print_all_references,&
@@ -411,7 +409,6 @@ CONTAINS
 
       ! Initialize the parallel random number generator
 
-      CALL init_rng()
       iw = cp_print_key_unit_nr(logger, root_section, "GLOBAL%PRINT/RNG_MATRICES", &
                                 extension=".Log")
       IF (iw > 0) THEN
@@ -432,11 +429,11 @@ CONTAINS
          CPABORT("Supply exactly 1 or 6 arguments for SEED in &GLOBAL only!")
       END IF
 
-      CALL create_rng_stream(rng_stream=globenv%gaussian_rng_stream, &
-                             name="Global Gaussian random numbers", &
-                             distribution_type=GAUSSIAN, &
-                             seed=initial_seed, &
-                             extended_precision=.TRUE.)
+      globenv%gaussian_rng_stream = rng_stream_type( &
+                                    name="Global Gaussian random numbers", &
+                                    distribution_type=GAUSSIAN, &
+                                    seed=initial_seed, &
+                                    extended_precision=.TRUE.)
 
       iw = cp_print_key_unit_nr(logger, root_section, "GLOBAL%PRINT/RNG_CHECK", &
                                 extension=".Log")
@@ -449,9 +446,9 @@ CONTAINS
 
       iw = cp_print_key_unit_nr(logger, root_section, "GLOBAL%PRINT/GLOBAL_GAUSSIAN_RNG", &
                                 extension=".Log")
-      IF (iw > 0) THEN
-         CALL write_rng_stream(globenv%gaussian_rng_stream, iw, write_all=.TRUE.)
-      END IF
+      IF (iw > 0) &
+         CALL globenv%gaussian_rng_stream%write(iw, write_all=.TRUE.)
+
       CALL cp_print_key_finished_output(iw, logger, root_section, &
                                         "GLOBAL%PRINT/GLOBAL_GAUSSIAN_RNG")
 
diff --git a/src/ewalds_multipole.F b/src/ewalds_multipole.F
index beb5c01618..cc3fb1dec8 100644
--- a/src/ewalds_multipole.F
+++ b/src/ewalds_multipole.F
@@ -44,9 +44,6 @@ MODULE ewalds_multipole
    USE mathlib,                         ONLY: matvec_3x3
    USE message_passing,                 ONLY: mp_sum
    USE parallel_rng_types,              ONLY: UNIFORM,&
-                                              create_rng_stream,&
-                                              delete_rng_stream,&
-                                              random_numbers,&
                                               rng_stream_type
    USE particle_types,                  ONLY: particle_type
    USE pw_grid_types,                   ONLY: pw_grid_type
@@ -1322,14 +1319,14 @@ $: ewalds_multipole_sr_macro(mode="SCREENED_COULOMB_ERF", store_energy=True, sto
                                                 e_field2
       REAL(KIND=dp), POINTER, &
          DIMENSION(:, :, :)                  :: quadrupoles
-      TYPE(rng_stream_type), POINTER           :: random_stream
+      TYPE(rng_stream_type)                    :: random_stream
       TYPE(multi_charge_type), DIMENSION(:), &
          POINTER                             :: multipoles
 
-      NULLIFY (random_stream, multipoles, charges, dipoles, g_forces, g_pv, &
+      NULLIFY (multipoles, charges, dipoles, g_forces, g_pv, &
                r_forces, r_pv, e_field1, e_field2)
-      CALL create_rng_stream(random_stream, name="DEBUG_EWALD_MULTIPOLE", &
-                             distribution_type=UNIFORM)
+      random_stream = rng_stream_type(name="DEBUG_EWALD_MULTIPOLE", &
+                                      distribution_type=UNIFORM)
       ! check:  charge - charge
       task = .FALSE.
       nparticles = SIZE(particle_set)
@@ -1553,7 +1550,6 @@ $: ewalds_multipole_sr_macro(mode="SCREENED_COULOMB_ERF", store_energy=True, sto
       DEALLOCATE (e_field2)
       DEALLOCATE (g_pv)
       DEALLOCATE (r_pv)
-      CALL delete_rng_stream(random_stream)
 
    CONTAINS
 ! **************************************************************************************************
@@ -1704,7 +1700,7 @@ $: ewalds_multipole_sr_macro(mode="SCREENED_COULOMB_ERF", store_energy=True, sto
       INTEGER, INTENT(IN)                                :: idim, istart, iend
       CHARACTER(LEN=*), INTENT(IN)                       :: label
       REAL(KIND=dp), INTENT(IN)                          :: echarge
-      TYPE(rng_stream_type), POINTER                     :: random_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: random_stream
       REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: charges
       REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: dipoles
       REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
@@ -1748,7 +1744,7 @@ $: ewalds_multipole_sr_macro(mode="SCREENED_COULOMB_ERF", store_energy=True, sto
                CPASSERT(ASSOCIATED(dipoles))
                ALLOCATE (multipoles(i)%charge_typ(isize)%charge(2))
                ALLOCATE (multipoles(i)%charge_typ(isize)%pos(3, 2))
-               CALL random_numbers(rvec, random_stream)
+               CALL random_stream%fill(rvec)
                rvec = rvec/(2.0_dp*SQRT(DOT_PRODUCT(rvec, rvec)))*dx
                multipoles(i)%charge_typ(isize)%charge(1) = echarge
                multipoles(i)%charge_typ(isize)%pos(1:3, 1) = rvec
@@ -1762,8 +1758,8 @@ $: ewalds_multipole_sr_macro(mode="SCREENED_COULOMB_ERF", store_energy=True, sto
                CPASSERT(ASSOCIATED(quadrupoles))
                ALLOCATE (multipoles(i)%charge_typ(isize)%charge(4))
                ALLOCATE (multipoles(i)%charge_typ(isize)%pos(3, 4))
-               CALL random_numbers(rvec1, random_stream)
-               CALL random_numbers(rvec2, random_stream)
+               CALL random_stream%fill(rvec1)
+               CALL random_stream%fill(rvec2)
                rvec1 = rvec1/SQRT(DOT_PRODUCT(rvec1, rvec1))
                rvec2 = rvec2-DOT_PRODUCT(rvec2, rvec1)*rvec1
                rvec2 = rvec2/SQRT(DOT_PRODUCT(rvec2, rvec2))
diff --git a/src/fm/cp_fm_types.F b/src/fm/cp_fm_types.F
index 5ac3626ae5..93c30ac3d0 100644
--- a/src/fm/cp_fm_types.F
+++ b/src/fm/cp_fm_types.F
@@ -31,11 +31,6 @@ MODULE cp_fm_types
         cp2k_is_parallel, mp_allgather, mp_any_source, mp_bcast, mp_irecv, mp_isend, mp_max, &
         mp_proc_null, mp_recv, mp_request_null, mp_send, mp_sum, mp_waitall
    USE parallel_rng_types,              ONLY: UNIFORM,&
-                                              create_rng_stream,&
-                                              delete_rng_stream,&
-                                              get_rng_stream,&
-                                              random_numbers,&
-                                              reset_to_next_rng_substream,&
                                               rng_stream_type
 #include "../base/base_uses.f90"
 
@@ -299,7 +294,7 @@ CONTAINS
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: buff
       REAL(KIND=dp), DIMENSION(3, 2), SAVE               :: seed
       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: local_data
-      TYPE(rng_stream_type), POINTER                     :: rng
+      TYPE(rng_stream_type)                              :: rng
 
       CALL timeset(routineN, handle)
 
@@ -311,9 +306,8 @@ CONTAINS
       ! guarantee same seed over all tasks
       CALL mp_bcast(seed, 0, matrix%matrix_struct%para_env%group)
 
-      NULLIFY (rng)
-      CALL create_rng_stream(rng, "cp_fm_init_random_stream", distribution_type=UNIFORM, &
-                             extended_precision=.TRUE., seed=seed)
+      rng = rng_stream_type("cp_fm_init_random_stream", distribution_type=UNIFORM, &
+                            extended_precision=.TRUE., seed=seed)
 
       CPASSERT(.NOT. matrix%use_sp)
 
@@ -339,11 +333,11 @@ CONTAINS
       DO icol_local = 1, ncol_local
          CPASSERT(col_indices(icol_local) > icol_global)
          DO
-            CALL reset_to_next_rng_substream(rng)
+            CALL rng%reset_to_next_substream()
             icol_global = icol_global + 1
             IF (icol_global == col_indices(icol_local)) EXIT
          ENDDO
-         CALL random_numbers(buff, rng)
+         CALL rng%fill(buff)
          DO irow_local = 1, nrow_local
             local_data(irow_local, icol_local) = buff(row_indices(irow_local))
          ENDDO
@@ -352,8 +346,7 @@ CONTAINS
       DEALLOCATE (buff)
 
       ! store seed before deletion (unclear if this is the proper seed)
-      CALL get_rng_stream(rng, ig=seed)
-      CALL delete_rng_stream(rng)
+      CALL rng%get(ig=seed)
 
       CALL timestop(handle)
 
diff --git a/src/gle_system_types.F b/src/gle_system_types.F
index f124346672..c453366be2 100644
--- a/src/gle_system_types.F
+++ b/src/gle_system_types.F
@@ -19,8 +19,6 @@ MODULE gle_system_types
                                               section_vals_val_get
    USE kinds,                           ONLY: dp
    USE parallel_rng_types,              ONLY: GAUSSIAN,&
-                                              create_rng_stream,&
-                                              delete_rng_stream,&
                                               next_rng_seed,&
                                               rng_stream_type
    USE string_utilities,                ONLY: compress
@@ -38,7 +36,7 @@ MODULE gle_system_types
       INTEGER                                 :: degrees_of_freedom
       REAL(KIND=dp)                           :: nkt, kin_energy, thermostat_energy
       REAL(KIND=dp), DIMENSION(:), POINTER    :: s
-      TYPE(rng_stream_type), POINTER          :: gaussian_rng_stream
+      TYPE(rng_stream_type)                   :: gaussian_rng_stream
    END TYPE gle_thermo_type
 
 ! **************************************************************************************************
@@ -202,14 +200,12 @@ CONTAINS
       ! Update initial seed
       initial_seed = next_rng_seed(seed(:, :, gle%glob_num_gle))
       DO ithermo = 1, gle%loc_num_gle
-         NULLIFY (gle%nvt(ithermo)%gaussian_rng_stream)
          my_index = gle%map_info%index(ithermo)
          my_seed = seed(:, :, my_index)
          WRITE (UNIT=name, FMT="(A,I8)") "Wiener process for Thermostat #", my_index
          CALL compress(name)
-         CALL create_rng_stream(rng_stream=gle%nvt(ithermo)%gaussian_rng_stream, &
-                                name=name, distribution_type=GAUSSIAN, extended_precision=.TRUE., &
-                                seed=my_seed)
+         gle%nvt(ithermo)%gaussian_rng_stream = rng_stream_type( &
+                                                name=name, distribution_type=GAUSSIAN, extended_precision=.TRUE., seed=my_seed)
       END DO
 
       DEALLOCATE (seed)
@@ -243,9 +239,6 @@ CONTAINS
          IF (ASSOCIATED(gle%nvt)) THEN
             DO i = 1, SIZE(gle%nvt)
                DEALLOCATE (gle%nvt(i)%s)
-               IF (ASSOCIATED(gle%nvt(i)%gaussian_rng_stream)) THEN
-                  CALL delete_rng_stream(gle%nvt(i)%gaussian_rng_stream)
-               END IF
             END DO
             DEALLOCATE (gle%nvt)
          ENDIF
diff --git a/src/global_types.F b/src/global_types.F
index 67bd4c0f58..25037c6b86 100644
--- a/src/global_types.F
+++ b/src/global_types.F
@@ -23,8 +23,7 @@ MODULE global_types
                                               default_string_length,&
                                               dp
    USE machine,                         ONLY: m_walltime
-   USE parallel_rng_types,              ONLY: delete_rng_stream,&
-                                              rng_stream_type
+   USE parallel_rng_types,              ONLY: rng_stream_type
 #include "./base/base_uses.f90"
 
    IMPLICIT NONE
@@ -65,7 +64,7 @@ MODULE global_types
 ! **************************************************************************************************
    TYPE global_environment_type
       INTEGER :: id_nr, ref_count
-      TYPE(rng_stream_type), POINTER          :: gaussian_rng_stream
+      TYPE(rng_stream_type), ALLOCATABLE      :: gaussian_rng_stream
       CHARACTER(LEN=default_string_length)    :: diag_library
       CHARACTER(LEN=default_string_length)    :: default_fft_library
       CHARACTER(LEN=default_path_length)      :: fftw_wisdom_file_name
@@ -111,7 +110,6 @@ CONTAINS
       globenv%idum = 0 !! random number seed
       globenv%blacs_grid_layout = BLACS_GRID_SQUARE
       globenv%cp2k_start_time = m_walltime()
-      NULLIFY (globenv%gaussian_rng_stream)
    END SUBROUTINE globenv_create
 
 ! **************************************************************************************************
@@ -144,9 +142,8 @@ CONTAINS
          CPASSERT(globenv%ref_count > 0)
          globenv%ref_count = globenv%ref_count - 1
          IF (globenv%ref_count == 0) THEN
-            IF (ASSOCIATED(globenv%gaussian_rng_stream)) THEN
-               CALL delete_rng_stream(globenv%gaussian_rng_stream)
-            END IF
+            IF (ALLOCATED(globenv%gaussian_rng_stream)) &
+               DEALLOCATE (globenv%gaussian_rng_stream)
             DEALLOCATE (globenv)
          END IF
       END IF
diff --git a/src/hfx_load_balance_methods.F b/src/hfx_load_balance_methods.F
index 9c16226518..4350ffe054 100644
--- a/src/hfx_load_balance_methods.F
+++ b/src/hfx_load_balance_methods.F
@@ -32,9 +32,6 @@ MODULE hfx_load_balance_methods
                                               mp_sync,&
                                               mp_waitall
    USE parallel_rng_types,              ONLY: UNIFORM,&
-                                              create_rng_stream,&
-                                              delete_rng_stream,&
-                                              next_random_number,&
                                               rng_stream_type
    USE particle_types,                  ONLY: particle_type
    USE util,                            ONLY: sort
@@ -1849,7 +1846,7 @@ CONTAINS
       INTEGER                                            :: i, itmp, j, nstep
       INTEGER(int_8), DIMENSION(:), POINTER              :: my_cost_cpu, tmp_cost, tmp_cpu_cost
       INTEGER, DIMENSION(:), POINTER                     :: tmp_cpu_index, tmp_index
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), ALLOCATABLE                 :: rng_stream
 
       nstep = MAX(1, INT(number_of_processes)/2)
 
@@ -1870,12 +1867,9 @@ CONTAINS
       ! it also avoids degenerate cases where thousands of zero sized tasks
       ! are assigned to the same (least loaded) cpu
       !
-      IF (do_randomize) THEN
-         NULLIFY (rng_stream)
-         CALL create_rng_stream(rng_stream=rng_stream, &
-                                name="uniform_rng", &
-                                distribution_type=UNIFORM)
-      END IF
+      IF (do_randomize) &
+         rng_stream = rng_stream_type(name="uniform_rng", &
+                                      distribution_type=UNIFORM)
 
       DO i = total_number_of_bins, 1, -nstep
          tmp_cpu_cost = my_cost_cpu
@@ -1890,10 +1884,6 @@ CONTAINS
          ENDDO
       ENDDO
 
-      IF (do_randomize) THEN
-         CALL delete_rng_stream(rng_stream)
-      END IF
-
       DEALLOCATE (tmp_cost, tmp_index, tmp_cpu_cost)
       DEALLOCATE (tmp_cpu_index, my_cost_cpu)
    END SUBROUTINE optimize_distribution
@@ -2538,15 +2528,15 @@ CONTAINS
    SUBROUTINE reshuffle(size, array, rng_stream)
       INTEGER                                            :: size
       INTEGER, DIMENSION(size)                           :: array
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type)                              :: rng_stream
 
       INTEGER                                            :: i, idx1, idx2, tmp
       REAL(dp)                                           :: x
 
       DO i = 1, size*10
-         x = next_random_number(rng_stream)
+         x = rng_stream%next()
          idx1 = INT(x*(size + 1 - 1)) + 1
-         x = next_random_number(rng_stream)
+         x = rng_stream%next()
          idx2 = INT(x*(size + 1 - 1)) + 1
 
          tmp = array(idx1)
diff --git a/src/input_restart_force_eval.F b/src/input_restart_force_eval.F
index 923471f214..6c47a3fc8a 100644
--- a/src/input_restart_force_eval.F
+++ b/src/input_restart_force_eval.F
@@ -49,8 +49,7 @@ MODULE input_restart_force_eval
    USE molecule_types,                  ONLY: get_molecule,&
                                               molecule_type
    USE multipole_types,                 ONLY: multipole_type
-   USE parallel_rng_types,              ONLY: dump_rng_stream,&
-                                              rng_record_length
+   USE parallel_rng_types,              ONLY: rng_record_length
    USE particle_list_types,             ONLY: particle_list_type
    USE qmmm_ff_fist,                    ONLY: qmmm_ff_precond_only_qm
    USE qs_environment_types,            ONLY: get_qs_env
@@ -265,9 +264,8 @@ CONTAINS
                nparticle_local = local_particles%n_el(iparticle_kind)
                DO iparticle_local = 1, nparticle_local
                   IF (iparticle == local_particles%list(iparticle_kind)%array(iparticle_local)) THEN
-                     CALL dump_rng_stream(rng_stream=local_particles%local_particle_set(iparticle_kind)% &
-                                          rng(iparticle_local)%stream, &
-                                          rng_record=rng_record)
+                     CALL local_particles%local_particle_set(iparticle_kind)% &
+                        rng(iparticle_local)%stream%dump(rng_record=rng_record)
                      CALL string_to_ascii(rng_record, ascii(:, iparticle))
                   END IF
                END DO
diff --git a/src/library_tests.F b/src/library_tests.F
index 83fb32a160..c7b46cd9e5 100644
--- a/src/library_tests.F
+++ b/src/library_tests.F
@@ -78,14 +78,8 @@ MODULE library_tests
                                               mpi_perf_test
    USE minimax_exp,                     ONLY: validate_exp_minimax
    USE mp2_weights,                     ONLY: test_least_square_ft
-   USE parallel_rng_types,              ONLY: GAUSSIAN,&
-                                              UNIFORM,&
-                                              check_rng,&
-                                              create_rng_stream,&
-                                              delete_rng_stream,&
-                                              next_random_number,&
-                                              rng_stream_type,&
-                                              write_rng_stream
+   USE parallel_rng_types,              ONLY: UNIFORM,&
+                                              rng_stream_type
    USE pw_grid_types,                   ONLY: FULLSPACE,&
                                               HALFSPACE,&
                                               pw_grid_type
@@ -174,8 +168,6 @@ CONTAINS
       !
       IF (runtest(8) /= 0) CALL mpi_perf_test(para_env%group, runtest(8), iw)
       !
-      IF (runtest(9) /= 0) CALL rng_test(para_env, iw)
-      !
       IF (runtest(10) /= 0) CALL validate_exp_minimax(runtest(10), iw)
       !
       IF (runtest(11) /= 0) CALL test_least_square_ft(runtest(11), iw)
@@ -275,7 +267,6 @@ CONTAINS
       CALL section_vals_val_get(test_section, 'ERI', i_val=runtest(4))
       CALL section_vals_val_get(test_section, 'CLEBSCH_GORDON', i_val=runtest(6))
       CALL section_vals_val_get(test_section, 'MPI', i_val=runtest(8))
-      CALL section_vals_val_get(test_section, 'RNG', i_val=runtest(9))
       CALL section_vals_val_get(test_section, 'MINIMAX', i_val=runtest(10))
       CALL section_vals_val_get(test_section, 'LEAST_SQ_FT', i_val=runtest(11))
 
@@ -1041,113 +1032,6 @@ CONTAINS
 
    END SUBROUTINE pw_fft_test
 
-! **************************************************************************************************
-!> \brief Test the parallel (pseudo)random number generator (RNG).
-!> \param para_env ...
-!> \param output_unit ...
-!> \par History
-!>      JGH  6-Feb-2001 : Test and performance code
-!> \author JGH  1-JAN-2001
-! **************************************************************************************************
-   SUBROUTINE rng_test(para_env, output_unit)
-      TYPE(cp_para_env_type), POINTER                    :: para_env
-      INTEGER                                            :: output_unit
-
-      CHARACTER(LEN=*), PARAMETER :: routineN = 'rng_test', routineP = moduleN//':'//routineN
-
-      INTEGER                                            :: i, n
-      LOGICAL                                            :: ionode
-      REAL(KIND=dp)                                      :: t, tend, tmax, tmin, tstart, tsum, tsum2
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
-
-      ionode = para_env%ionode
-      n = runtest(9)
-      NULLIFY (rng_stream)
-
-      ! Check correctness
-
-      CALL check_rng(output_unit, ionode)
-
-      ! Check performance
-
-      IF (ionode) THEN
-         WRITE (UNIT=output_unit, FMT="(/,/,T2,A,I10,A)") &
-            "Check distributions using", n, " random numbers:"
-      END IF
-
-      ! Test uniform distribution [0,1]
-
-      CALL create_rng_stream(rng_stream=rng_stream, &
-                             name="Test uniform distribution [0,1]", &
-                             distribution_type=UNIFORM, &
-                             extended_precision=.TRUE.)
-
-      IF (ionode) THEN
-         CALL write_rng_stream(rng_stream, output_unit, write_all=.TRUE.)
-      END IF
-
-      tmax = -HUGE(0.0_dp)
-      tmin = +HUGE(0.0_dp)
-      tsum = 0.0_dp
-      tsum2 = 0.0_dp
-
-      tstart = m_walltime()
-      DO i = 1, n
-         t = next_random_number(rng_stream)
-         tsum = tsum + t
-         tsum2 = tsum2 + t*t
-         IF (t > tmax) tmax = t
-         IF (t < tmin) tmin = t
-      END DO
-      tend = m_walltime()
-
-      IF (ionode) THEN
-         WRITE (UNIT=output_unit, FMT="(/,(T4,A,F12.6))") &
-            "Minimum: ", tmin, &
-            "Maximum: ", tmax, &
-            "Average: ", tsum/REAL(n, KIND=dp), &
-            "Variance:", tsum2/REAL(n, KIND=dp), &
-            "Time [s]:", tend - tstart
-      END IF
-
-      CALL delete_rng_stream(rng_stream)
-
-      ! Test normal Gaussian distribution
-
-      CALL create_rng_stream(rng_stream=rng_stream, &
-                             name="Test normal Gaussian distribution", &
-                             distribution_type=GAUSSIAN, &
-                             extended_precision=.TRUE.)
-
-      tmax = -HUGE(0.0_dp)
-      tmin = +HUGE(0.0_dp)
-      tsum = 0.0_dp
-      tsum2 = 0.0_dp
-
-      tstart = m_walltime()
-      DO i = 1, n
-         t = next_random_number(rng_stream)
-         tsum = tsum + t
-         tsum2 = tsum2 + t*t
-         IF (t > tmax) tmax = t
-         IF (t < tmin) tmin = t
-      END DO
-      tend = m_walltime()
-
-      IF (ionode) THEN
-         CALL write_rng_stream(rng_stream, output_unit)
-         WRITE (UNIT=output_unit, FMT="(/,(T4,A,F12.6))") &
-            "Minimum: ", tmin, &
-            "Maximum: ", tmax, &
-            "Average: ", tsum/REAL(n, KIND=dp), &
-            "Variance:", tsum2/REAL(n, KIND=dp), &
-            "Time [s]:", tend - tstart
-      END IF
-
-      CALL delete_rng_stream(rng_stream)
-
-   END SUBROUTINE rng_test
-
 ! **************************************************************************************************
 !> \brief Tests the eigensolver library routines
 !> \param para_env ...
@@ -1175,7 +1059,7 @@ CONTAINS
       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
       TYPE(cp_fm_struct_type), POINTER                   :: fmstruct
       TYPE(cp_fm_type), POINTER                          :: eigenvectors, matrix, work
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), ALLOCATABLE                 :: rng_stream
 
       group = para_env%group
       source = para_env%source
@@ -1263,11 +1147,10 @@ CONTAINS
          IF (para_env%ionode) THEN
             SELECT CASE (init_method)
             CASE (do_mat_random)
-               NULLIFY (rng_stream)
-               CALL create_rng_stream(rng_stream=rng_stream, &
-                                      name="rng_stream", &
-                                      distribution_type=UNIFORM, &
-                                      extended_precision=.TRUE.)
+               rng_stream = rng_stream_type( &
+                            name="rng_stream", &
+                            distribution_type=UNIFORM, &
+                            extended_precision=.TRUE.)
             CASE (do_mat_read)
                CALL open_file(file_name="MATRIX", &
                               file_action="READ", &
@@ -1282,7 +1165,7 @@ CONTAINS
                SELECT CASE (init_method)
                CASE (do_mat_random)
                   DO j = i, n
-                     buffer(1, j) = next_random_number(rng_stream) - 0.5_dp
+                     buffer(1, j) = rng_stream%next() - 0.5_dp
                   END DO
                   !MK activate/modify for a diagonal dominant symmetric matrix:
                   !MK buffer(1,i) = 10.0_dp*buffer(1,i)
@@ -1328,8 +1211,6 @@ CONTAINS
 
          IF (para_env%ionode) THEN
             SELECT CASE (init_method)
-            CASE (do_mat_random)
-               CALL delete_rng_stream(rng_stream=rng_stream)
             CASE (do_mat_read)
                CALL close_file(unit_number=unit_number)
             END SELECT
diff --git a/src/metadynamics.F b/src/metadynamics.F
index 2fc5e8b53a..0a589d0131 100644
--- a/src/metadynamics.F
+++ b/src/metadynamics.F
@@ -50,7 +50,6 @@ MODULE metadynamics
                                  meta_walls, &
                                  restart_hills, &
                                  synchronize_multiple_walkers
-   USE parallel_rng_types, ONLY: next_random_number
    USE particle_list_types, ONLY: particle_list_type
 #if defined (__PLUMED2)
    USE physcon, ONLY: angstrom, &
@@ -401,7 +400,7 @@ CONTAINS
          meta_env%ekin_s = 0.0_dp
          DO i_c = 1, meta_env%n_colvar
             cv => meta_env%metavar(i_c)
-            cv%vvp = next_random_number(force_env%globenv%gaussian_rng_stream)
+            cv%vvp = force_env%globenv%gaussian_rng_stream%next()
             meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
          END DO
          ekin_w = 0.5_dp*meta_env%temp_wanted*REAL(meta_env%n_colvar, KIND=dp)
diff --git a/src/metadynamics_types.F b/src/metadynamics_types.F
index 5a5bfd40cc..d46fb4987b 100644
--- a/src/metadynamics_types.F
+++ b/src/metadynamics_types.F
@@ -16,8 +16,7 @@ MODULE metadynamics_types
                                               section_vals_val_get
    USE kinds,                           ONLY: default_path_length,&
                                               dp
-   USE parallel_rng_types,              ONLY: delete_rng_stream,&
-                                              rng_stream_p_type
+   USE parallel_rng_types,              ONLY: rng_stream_type
 #include "./base/base_uses.f90"
 
    IMPLICIT NONE
@@ -118,8 +117,8 @@ MODULE metadynamics_types
       TYPE(multiple_walkers_type), POINTER      :: multiple_walkers
       TYPE(cp_para_env_type), POINTER           :: para_env
       TYPE(section_vals_type), POINTER          :: metadyn_section
-      TYPE(rng_stream_p_type), DIMENSION(:), &
-         POINTER :: rng
+      TYPE(rng_stream_type), DIMENSION(:), &
+         ALLOCATABLE :: rng
       INTEGER                                   :: TAMCSteps
       REAL(kind=dp)                             :: zdt
    END TYPE meta_env_type
@@ -160,7 +159,6 @@ CONTAINS
       NULLIFY (meta_env%multiple_walkers, &
                meta_env%metadyn_section, &
                meta_env%time, &
-               meta_env%rng, &
                meta_env%hills_env)
 
       meta_env%use_plumed = .FALSE.
@@ -227,9 +225,6 @@ CONTAINS
       CALL section_vals_val_get(metadyn_section, "LANGEVIN", l_val=do_langevin)
       IF (do_langevin) THEN
          ALLOCATE (meta_env%rng(meta_env%n_colvar))
-         DO i = 1, meta_env%n_colvar
-            NULLIFY (meta_env%rng(meta_env%n_colvar)%stream)
-         END DO
       ENDIF
    END SUBROUTINE metadyn_create
 
@@ -310,14 +305,9 @@ CONTAINS
             END IF
 
             ! Langevin on COLVARS
-            IF (meta_env%langevin) THEN
-               DO i = 1, SIZE(meta_env%rng)
-                  IF (ASSOCIATED(meta_env%rng(i)%stream)) THEN
-                     CALL delete_rng_stream(meta_env%rng(i)%stream)
-                  END IF
-               END DO
+            IF (meta_env%langevin) &
                DEALLOCATE (meta_env%rng)
-            ENDIF
+
             NULLIFY (meta_env%time)
             NULLIFY (meta_env%metadyn_section)
             DEALLOCATE (meta_env)
diff --git a/src/mode_selective.F b/src/mode_selective.F
index 922b58b6ca..f851e5c3bc 100644
--- a/src/mode_selective.F
+++ b/src/mode_selective.F
@@ -42,7 +42,6 @@ MODULE mode_selective
    USE mathlib,                         ONLY: diamat_all
    USE message_passing,                 ONLY: mp_bcast
    USE molden_utils,                    ONLY: write_vibrations_molden
-   USE parallel_rng_types,              ONLY: next_random_number
    USE particle_types,                  ONLY: particle_type
    USE physcon,                         ONLY: massunit,&
                                               vibfac
@@ -370,7 +369,7 @@ CONTAINS
             DO j = 1, natoms
                DO k = 1, 3
                   jj = (map_atoms(j) - 1)*3 + k
-                  ms_vib%b_vec(jj, i) = ABS(next_random_number(globenv%gaussian_rng_stream))
+                  ms_vib%b_vec(jj, i) = ABS(globenv%gaussian_rng_stream%next())
                END DO
             END DO
             norm = SQRT(DOT_PRODUCT(ms_vib%b_vec(:, i), ms_vib%b_vec(:, i)))
diff --git a/src/motion/dimer_types.F b/src/motion/dimer_types.F
index 954d7d47e9..0a4e4b8b78 100644
--- a/src/motion/dimer_types.F
+++ b/src/motion/dimer_types.F
@@ -32,7 +32,6 @@ MODULE dimer_types
    USE molecule_kind_types,             ONLY: fixd_constraint_type,&
                                               get_molecule_kind,&
                                               molecule_kind_type
-   USE parallel_rng_types,              ONLY: random_numbers
 #include "../base/base_uses.f90"
 
    IMPLICIT NONE
@@ -162,7 +161,7 @@ CONTAINS
          END DO
          CPASSERT(isize == SIZE(dimer_env%nvec))
       ELSE
-         CALL random_numbers(dimer_env%nvec, globenv%gaussian_rng_stream)
+         CALL globenv%gaussian_rng_stream%fill(dimer_env%nvec)
       END IF
       ! Check for translation in the dimer vector and remove them
       IF (natom > 1) THEN
diff --git a/src/motion/helium_common.F b/src/motion/helium_common.F
index 63b184d542..6fac62e953 100644
--- a/src/motion/helium_common.F
+++ b/src/motion/helium_common.F
@@ -27,7 +27,6 @@ MODULE helium_common
                                               dp
    USE mathconstants,                   ONLY: pi
    USE memory_utilities,                ONLY: reallocate
-   USE parallel_rng_types,              ONLY: next_random_number
    USE physcon,                         ONLY: angstrom,&
                                               bohr
    USE pint_public,                     ONLY: pint_com_pos
@@ -1339,7 +1338,7 @@ CONTAINS
          ! number of random numbers to generate:
          c = 1000000000
          DO j = 1, c
-            v = next_random_number(helium%rng_stream_uniform)
+            v = helium%rng_stream_uniform%next()
             ! walk down the search tree:
             k = nb - 1
             DO
@@ -1418,7 +1417,7 @@ CONTAINS
       ! (should not be taken, but just in case it does we have something valid)
 
       helium%pweight = 0.0_dp
-      t = next_random_number(helium%rng_stream_uniform)
+      t = helium%rng_stream_uniform%next()
       helium%ptable(1) = 1 + INT(t*nb)
       helium%ptable(2) = -1
 
diff --git a/src/motion/helium_methods.F b/src/motion/helium_methods.F
index e308dbfd2f..da524dd5c7 100644
--- a/src/motion/helium_methods.F
+++ b/src/motion/helium_methods.F
@@ -62,12 +62,8 @@ MODULE helium_methods
                                               mp_comm_split_direct
    USE parallel_rng_types,              ONLY: GAUSSIAN,&
                                               UNIFORM,&
-                                              create_rng_stream,&
-                                              delete_rng_stream,&
-                                              next_random_number,&
                                               rng_stream_p_type,&
-                                              rng_stream_type,&
-                                              set_rng_stream
+                                              rng_stream_type
    USE particle_list_types,             ONLY: particle_list_type
    USE physcon,                         ONLY: a_mass,&
                                               angstrom,&
@@ -967,8 +963,8 @@ CONTAINS
                         helium_env(k)%helium%u0, &
                         helium_env(k)%helium%e0)
 
-               CALL delete_rng_stream(helium_env(k)%helium%rng_stream_uniform)
-               CALL delete_rng_stream(helium_env(k)%helium%rng_stream_gaussian)
+               DEALLOCATE (helium_env(k)%helium%rng_stream_uniform)
+               DEALLOCATE (helium_env(k)%helium%rng_stream_gaussian)
 
                ! deallocate solute-related arrays
                IF (helium_env(k)%helium%solute_present) THEN
@@ -1222,7 +1218,7 @@ CONTAINS
                   !minHeHedsttmp = 0.90_dp**(iter/100)*minHeHedst
                   minHeHedsttmp = 0.90_dp**MIN(0, iter - 2)*minHeHedst
                   DO ic = 1, 3
-                     r1 = next_random_number(helium_env(k)%helium%rng_stream_uniform)
+                     r1 = helium_env(k)%helium%rng_stream_uniform%next()
                      r1 = 2.0_dp*r1 - 1.0_dp
                      r1 = r1*helium_env(k)%helium%cell_size
                      centroids(ic, ia) = r1
@@ -1292,7 +1288,7 @@ CONTAINS
                   ! if sampling fails to often, reduce he he criterion
                   minHeHedsttmp = 0.90_dp**MIN(0, iter - 2)*minHeHedst
                   DO ic = 1, 3
-                     rvek(ic) = next_random_number(helium_env(k)%helium%rng_stream_uniform)
+                     rvek(ic) = helium_env(k)%helium%rng_stream_uniform%next()
                      rvek(ic) = 2.0_dp*rvek(ic) - 1.0_dp
                      rvek(ic) = rvek(ic)*helium_env(k)%helium%droplet_radius
                   END DO
@@ -1419,7 +1415,7 @@ CONTAINS
             DO imode = 2, p
                omega = 2.0_dp*p*kbT*SIN((imode - 1)*pip)
                variance = kbT*p/(helium_env%he_mass_au*omega**2)
-               rand = next_random_number(helium_env%rng_stream_gaussian)
+               rand = helium_env%rng_stream_gaussian%next()
                nmhecoords(imode) = rand*SQRT(variance)
             END DO
             helium_env%pos(idim, iatom, 1:p) = MATMUL(u2x, nmhecoords)
@@ -1885,7 +1881,6 @@ CONTAINS
       REAL(KIND=dp), DIMENSION(3, 2)                     :: initial_seed
       TYPE(cp_logger_type), POINTER                      :: logger
       TYPE(rng_stream_p_type), DIMENSION(:), POINTER     :: gaussian_array, uniform_array
-      TYPE(rng_stream_type), POINTER                     :: next_rngs, prev_rngs
 
       NULLIFY (logger)
       logger => cp_get_default_logger()
@@ -1902,73 +1897,52 @@ CONTAINS
       ALLOCATE (uniform_array(helium_env(1)%helium%num_env), &
                 gaussian_array(helium_env(1)%helium%num_env))
       DO i = 1, helium_env(1)%helium%num_env
-         NULLIFY (uniform_array(i)%stream, gaussian_array(i)%stream)
+         ALLOCATE (uniform_array(i)%stream, gaussian_array(i)%stream)
       END DO
-      NULLIFY (prev_rngs, next_rngs)
 
       ! Create num_env RNG streams on processor all processors
       ! and distribute them so, that each processor gets unique
       ! RN sequences for his helium environments
       ! COMMENT: rng_stream can not be used with mp_bcast
 
-      CALL create_rng_stream(prev_rngs, &
-                             name="helium_rns_uniform", &
-                             distribution_type=UNIFORM, &
-                             extended_precision=.TRUE., &
-                             seed=initial_seed)
-      uniform_array(1)%stream => prev_rngs
-
-      CALL create_rng_stream(next_rngs, &
-                             name="helium_rns_gaussian", &
-                             last_rng_stream=prev_rngs, &
-                             distribution_type=GAUSSIAN, &
-                             extended_precision=.TRUE.)
-      gaussian_array(1)%stream => next_rngs
-      NULLIFY (prev_rngs)
-      prev_rngs => next_rngs
-      NULLIFY (next_rngs)
+      uniform_array(1)%stream = rng_stream_type(name="helium_rns_uniform", &
+                                                distribution_type=UNIFORM, &
+                                                extended_precision=.TRUE., &
+                                                seed=initial_seed)
 
+      gaussian_array(1)%stream = rng_stream_type(name="helium_rns_gaussian", &
+                                                 distribution_type=GAUSSIAN, &
+                                                 extended_precision=.TRUE., &
+                                                 last_rng_stream=uniform_array(1)%stream)
       DO i = 2, helium_env(1)%helium%num_env
-         CALL create_rng_stream(next_rngs, &
-                                name="helium_rns_uniform", &
-                                last_rng_stream=prev_rngs, &
-                                distribution_type=UNIFORM, &
-                                extended_precision=.TRUE.)
-         uniform_array(i)%stream => next_rngs
-         prev_rngs => next_rngs
-         NULLIFY (next_rngs)
-
-         CALL create_rng_stream(next_rngs, &
-                                name="helium_rns_gaussian", &
-                                last_rng_stream=prev_rngs, &
-                                distribution_type=GAUSSIAN, &
-                                extended_precision=.TRUE.)
-         gaussian_array(i)%stream => next_rngs
-         prev_rngs => next_rngs
-         NULLIFY (next_rngs)
+         uniform_array(i)%stream = rng_stream_type(name="helium_rns_uniform", &
+                                                   distribution_type=UNIFORM, &
+                                                   extended_precision=.TRUE., &
+                                                   last_rng_stream=gaussian_array(i - 1)%stream)
 
+         gaussian_array(i)%stream = rng_stream_type(name="helium_rns_uniform", &
+                                                    distribution_type=GAUSSIAN, &
+                                                    extended_precision=.TRUE., &
+                                                    last_rng_stream=uniform_array(i)%stream)
       END DO
 
-      NULLIFY (prev_rngs)
-
       offset = 0
       DO i = 1, logger%para_env%mepos
          offset = offset + helium_env(1)%env_all(i)
       END DO
 
-      IF (ASSOCIATED(helium_env)) THEN
-         DO i = 1, SIZE(helium_env)
-            NULLIFY (helium_env(i)%helium%rng_stream_uniform, &
-                     helium_env(i)%helium%rng_stream_gaussian)
-            helium_env(i)%helium%rng_stream_uniform => uniform_array(offset + i)%stream
-            helium_env(i)%helium%rng_stream_gaussian => gaussian_array(offset + i)%stream
-         END DO
-      END IF
+      DO i = 1, SIZE(helium_env)
+         NULLIFY (helium_env(i)%helium%rng_stream_uniform, &
+                  helium_env(i)%helium%rng_stream_gaussian)
+         helium_env(i)%helium%rng_stream_uniform => uniform_array(offset + i)%stream
+         helium_env(i)%helium%rng_stream_gaussian => gaussian_array(offset + i)%stream
+      END DO
 
       DO i = 1, helium_env(1)%helium%num_env
          IF (i .LE. offset .OR. i .GT. offset + SIZE(helium_env)) THEN
-            CALL delete_rng_stream(uniform_array(i)%stream)
-            CALL delete_rng_stream(gaussian_array(i)%stream)
+            ! only deallocate pointers here which were not passed on to helium_env(*)%helium
+            DEALLOCATE (uniform_array(i)%stream)
+            DEALLOCATE (gaussian_array(i)%stream)
          END IF
          NULLIFY (uniform_array(i)%stream)
          NULLIFY (gaussian_array(i)%stream)
@@ -1976,8 +1950,6 @@ CONTAINS
 
       DEALLOCATE (uniform_array)
       DEALLOCATE (gaussian_array)
-
-      RETURN
    END SUBROUTINE helium_rng_init
 
 ! ***************************************************************************
@@ -2060,8 +2032,8 @@ CONTAINS
             ELSE
                lbf = .FALSE.
             END IF
-            CALL set_rng_stream(helium_env(k)%helium%rng_stream_uniform, bg=bg, cg=cg, ig=ig, &
-                                buffer=bu, buffer_filled=lbf)
+            CALL helium_env(k)%helium%rng_stream_uniform%set(bg=bg, cg=cg, ig=ig, &
+                                                             buffer=bu, buffer_filled=lbf)
             bg(:, :) = UNPACK(message(off + 21:off + 26), MASK=m, FIELD=f)
             cg(:, :) = UNPACK(message(off + 27:off + 32), MASK=m, FIELD=f)
             ig(:, :) = UNPACK(message(off + 33:off + 38), MASK=m, FIELD=f)
@@ -2072,8 +2044,8 @@ CONTAINS
             ELSE
                lbf = .FALSE.
             END IF
-            CALL set_rng_stream(helium_env(k)%helium%rng_stream_gaussian, bg=bg, cg=cg, ig=ig, &
-                                buffer=bu, buffer_filled=lbf)
+            CALL helium_env(k)%helium%rng_stream_gaussian%set(bg=bg, cg=cg, ig=ig, &
+                                                              buffer=bu, buffer_filled=lbf)
          END DO
       END IF
 
diff --git a/src/motion/helium_sampling.F b/src/motion/helium_sampling.F
index e7e2ccaa51..f161acefd2 100644
--- a/src/motion/helium_sampling.F
+++ b/src/motion/helium_sampling.F
@@ -43,7 +43,6 @@ MODULE helium_sampling
    USE machine,                         ONLY: m_walltime
    USE message_passing,                 ONLY: mp_bcast,&
                                               mp_sum
-   USE parallel_rng_types,              ONLY: next_random_number
    USE physcon,                         ONLY: angstrom
    USE pint_public,                     ONLY: pint_com_pos
    USE pint_types,                      ONLY: pint_env_type
@@ -228,7 +227,7 @@ CONTAINS
                ! 'rotation state' in helium%relrot wich is within (0, helium%beads-1)
                ! (this is needed to sample different fragments of the permutation
                ! paths in try_permutations)
-               rnd = next_random_number(helium_env(k)%helium%rng_stream_uniform)
+               rnd = helium_env(k)%helium%rng_stream_uniform%next()
                nslices = INT(rnd*helium_env(k)%helium%beads)
                CALL helium_rotate(helium_env(k)%helium, nslices)
 
@@ -373,7 +372,7 @@ CONTAINS
             END DO
          ELSE IF (helium_env(1)%helium%get_helium_forces == helium_forces_last) THEN
             IF (logger%para_env%ionode) THEN
-               sel_mp_source = INT(next_random_number(helium_env(1)%helium%rng_stream_uniform) &
+               sel_mp_source = INT(helium_env(1)%helium%rng_stream_uniform%next() &
                                    *REAL(helium_env(1)%helium%num_env, dp))
             END IF
             CALL mp_bcast(sel_mp_source, logger%para_env%source, helium_env(1)%comm)
@@ -599,7 +598,7 @@ CONTAINS
          SELECT CASE (helium%m_dist_type)
 
          CASE (helium_mdist_singlev)
-            x = next_random_number(helium%rng_stream_uniform)
+            x = helium%rng_stream_uniform%next()
             IF (x .LT. r) THEN
                cyclen = 1
             ELSE
@@ -607,24 +606,24 @@ CONTAINS
             END IF
 
          CASE (helium_mdist_uniform)
-            x = next_random_number(helium%rng_stream_uniform)
+            x = helium%rng_stream_uniform%next()
             IF (x .LT. r) THEN
                cyclen = helium%m_value
             ELSE
                DO
-                  x = next_random_number(helium%rng_stream_uniform)
+                  x = helium%rng_stream_uniform%next()
                   cyclen = INT(helium%maxcycle*x) + 1
                   IF (cyclen .NE. helium%m_value) EXIT
                END DO
             END IF
 
          CASE (helium_mdist_linear)
-            x = next_random_number(helium%rng_stream_uniform)
+            x = helium%rng_stream_uniform%next()
             IF (x .LT. r) THEN
                cyclen = helium%m_value
             ELSE
                DO
-                  x = next_random_number(helium%rng_stream_uniform)
+                  x = helium%rng_stream_uniform%next()
                   y = SQRT(2.0_dp*x)
                   cyclen = INT(helium%maxcycle*y/SQRT(2.0_dp)) + 1
                   IF (cyclen .NE. helium%m_value) EXIT
@@ -632,12 +631,12 @@ CONTAINS
             END IF
 
          CASE (helium_mdist_quadratic)
-            x = next_random_number(helium%rng_stream_uniform)
+            x = helium%rng_stream_uniform%next()
             IF (x .LT. r) THEN
                cyclen = helium%m_value
             ELSE
                DO
-                  x = next_random_number(helium%rng_stream_uniform)
+                  x = helium%rng_stream_uniform%next()
                   y = (3.0_dp*x)**(1.0_dp/3.0_dp)
                   z = 3.0_dp**(1.0_dp/3.0_dp)
                   cyclen = INT(helium%maxcycle*y/z) + 1
@@ -646,13 +645,13 @@ CONTAINS
             END IF
 
          CASE (helium_mdist_exponential)
-            x = next_random_number(helium%rng_stream_uniform)
+            x = helium%rng_stream_uniform%next()
             IF (x .LT. r) THEN
                cyclen = helium%m_value
             ELSE
                DO
                   DO
-                     x = next_random_number(helium%rng_stream_uniform)
+                     x = helium%rng_stream_uniform%next()
                      IF (x .GE. 0.01_dp) EXIT
                   END DO
                   z = -LOG(0.01_dp)
@@ -663,12 +662,12 @@ CONTAINS
             END IF
 
          CASE (helium_mdist_gaussian)
-            x = next_random_number(helium%rng_stream_uniform)
+            x = helium%rng_stream_uniform%next()
             IF (x .LT. r) THEN
                cyclen = 1
             ELSE
                DO
-                  x = next_random_number(helium%rng_stream_gaussian)
+                  x = helium%rng_stream_gaussian%next()
                   cyclen = INT(x*0.75_dp + helium%m_value - 0.5_dp) + 1
                   IF (cyclen .NE. 1) EXIT
                END DO
@@ -688,7 +687,7 @@ CONTAINS
 
          ! check, if permutation of this length can be constructed
          IF (cyclen == 1) THEN
-            rnd = next_random_number(helium%rng_stream_uniform)
+            rnd = helium%rng_stream_uniform%next()
             helium%ptable(1) = 1 + INT(rnd*helium%atoms)
             helium%ptable(2) = -1
             helium%pweight = 0.0_dp
@@ -846,7 +845,7 @@ CONTAINS
             DO k = 1, cyclen
                CALL helium_boxmean_3d(helium, work(:, p(k), pk1), work(:, p(k), pk2), tmp1)
                DO c = 1, 3
-                  x = next_random_number(rng_stream=helium%rng_stream_gaussian, variance=1.0_dp)
+                  x = helium%rng_stream_gaussian%next(variance=1.0_dp)
                   x = sigma*x
                   tmp1(c) = tmp1(c) + x
                   tmp2(c) = x
@@ -872,7 +871,7 @@ CONTAINS
          DO k = 1, cyclen
             CALL helium_boxmean_3d(helium, work(:, p(k), pk1), work(:, perm(p(1 + MOD(k, cyclen))), 1), tmp1)
             DO c = 1, 3
-               x = next_random_number(rng_stream=helium%rng_stream_gaussian, variance=1.0_dp)
+               x = helium%rng_stream_gaussian%next(variance=1.0_dp)
                x = sigma*x
                tmp1(c) = tmp1(c) + x
                tmp2(c) = x
@@ -954,7 +953,7 @@ CONTAINS
             ds = ds - x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
          END DO
          ! ok now accept or reject:
-         rtmp = next_random_number(helium%rng_stream_uniform)
+         rtmp = helium%rng_stream_uniform%next()
 !      IF ((dtk+ds-pds < 0.0_dp).AND.(EXP(dtk+ds-pds)<rtmp)) THEN
          IF (dtk + ds - pds < 0.0_dp) THEN
             IF (EXP(dtk + ds - pds) < rtmp) THEN
@@ -993,7 +992,7 @@ CONTAINS
          DO k = 1, cyclen
             CALL helium_boxmean_3d(helium, work(:, p(k), pk1), work(:, p(k), pk2), tmp1)
             DO c = 1, 3
-               x = next_random_number(rng_stream=helium%rng_stream_gaussian, variance=1.0_dp)
+               x = helium%rng_stream_gaussian%next(variance=1.0_dp)
                x = sigma*x
                tmp1(c) = tmp1(c) + x
                tmp2(c) = x
@@ -1017,7 +1016,7 @@ CONTAINS
       DO k = 1, cyclen
          CALL helium_boxmean_3d(helium, work(:, p(k), pk1), work(:, perm(p(1 + MOD(k, cyclen))), 1), tmp1)
          DO c = 1, 3
-            x = next_random_number(rng_stream=helium%rng_stream_gaussian, variance=1.0_dp)
+            x = helium%rng_stream_gaussian%next(variance=1.0_dp)
             x = sigma*x
             tmp1(c) = tmp1(c) + x
             tmp2(c) = x
@@ -1127,7 +1126,7 @@ CONTAINS
       END DO
       DEALLOCATE (work3)
       ! ok now accept or reject:
-      rtmp = next_random_number(helium%rng_stream_uniform)
+      rtmp = helium%rng_stream_uniform%next()
 !   IF ((dtk+ds-pds<0.0_dp).AND.(EXP(dtk+ds-pds)<rtmp)) THEN
       IF (dtk + ds - pds < 0.0_dp) THEN
          IF (EXP(dtk + ds - pds) < rtmp) THEN
@@ -1256,10 +1255,10 @@ CONTAINS
       nmatrix => helium%nmatrix
 
       p(len + 1) = -1
-      rnd = next_random_number(helium%rng_stream_uniform)
+      rnd = helium%rng_stream_uniform%next()
       p(1) = INT(n*rnd) + 1
       DO k = 1, len - 1
-         t = next_random_number(helium%rng_stream_uniform)
+         t = helium%rng_stream_uniform%next()
          ! find the corresponding path to connect to
          ! using the precalculated optimal decision tree:
          i = n - 1
@@ -1287,7 +1286,7 @@ CONTAINS
       s1 = s1 + ipmatrix(p(len), perm(p(1)))
       s2 = s2 + ipmatrix(p(len), perm(p(len)))
       ! final accept/reject:
-      rnd = next_random_number(helium%rng_stream_uniform)
+      rnd = helium%rng_stream_uniform%next()
       t = s1*rnd
       IF (t > s2) RETURN
       ! ok, we have accepted the permutation
diff --git a/src/motion/helium_worm.F b/src/motion/helium_worm.F
index 2e7516360f..b0c870ffa0 100644
--- a/src/motion/helium_worm.F
+++ b/src/motion/helium_worm.F
@@ -25,7 +25,6 @@ MODULE helium_worm
    USE kinds,                           ONLY: default_string_length,&
                                               dp
    USE mathconstants,                   ONLY: pi
-   USE parallel_rng_types,              ONLY: next_random_number
    USE pint_types,                      ONLY: pint_env_type
 #include "../base/base_uses.f90"
 
@@ -96,11 +95,11 @@ CONTAINS
       IF (helium%worm_allow_open) THEN
          DO ! Exit criterion at the end of the loop
             DO iMC = 1, nMC
-               imove = next_random_number(helium%rng_stream_uniform, 1, helium%worm_all_limit)
+               imove = helium%rng_stream_uniform%next(1, helium%worm_all_limit)
                IF (helium%worm_is_closed) THEN
                   IF ((imove >= helium%worm_centroid_min) .AND. (imove <= helium%worm_centroid_max)) THEN
                      ! centroid move
-                     iatom = next_random_number(helium%rng_stream_uniform, 1, helium%atoms)
+                     iatom = helium%rng_stream_uniform%next(1, helium%atoms)
                      CALL worm_centroid_move(pint_env, helium, &
                                              iatom, helium%worm_centroid_drmax, ac)
                      ncentratt = ncentratt + 1
@@ -109,18 +108,18 @@ CONTAINS
                      ! staging is adjusted to conserve these weights
                   ELSE IF ((imove >= helium%worm_centroid_max + 1) .AND. (imove <= helium%worm_open_close_min - 1)) THEN
                      ! staging move
-                     iatom = next_random_number(helium%rng_stream_uniform, 1, helium%atoms)
-                     ibead = next_random_number(helium%rng_stream_uniform, 1, helium%beads)
-                     staging_l = next_random_number(helium%rng_stream_uniform, 2, helium%worm_staging_l)
+                     iatom = helium%rng_stream_uniform%next(1, helium%atoms)
+                     ibead = helium%rng_stream_uniform%next(1, helium%beads)
+                     staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
                      CALL worm_staging_move(pint_env, helium, &
                                             iatom, ibead, staging_l, ac)
                      nstagatt = nstagatt + 1
                      nstagacc = nstagacc + ac
                   ELSE IF ((imove >= helium%worm_open_close_min) .AND. (imove <= helium%worm_open_close_max)) THEN
                      ! attempt opening of worm
-                     iatom = next_random_number(helium%rng_stream_uniform, 1, helium%atoms)
-                     ibead = next_random_number(helium%rng_stream_uniform, 1, helium%beads)
-                     staging_l = next_random_number(helium%rng_stream_uniform, 2, helium%worm_staging_l)
+                     iatom = helium%rng_stream_uniform%next(1, helium%atoms)
+                     ibead = helium%rng_stream_uniform%next(1, helium%beads)
+                     staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
                      CALL worm_open_move(pint_env, helium, &
                                          iatom, ibead, staging_l, ac)
                      nopenatt = nopenatt + 1
@@ -132,16 +131,16 @@ CONTAINS
                ELSE ! worm is open
                   IF ((imove >= helium%worm_centroid_min) .AND. (imove <= helium%worm_centroid_max)) THEN
                      ! centroid move
-                     iatom = next_random_number(helium%rng_stream_uniform, 1, helium%atoms)
+                     iatom = helium%rng_stream_uniform%next(1, helium%atoms)
                      CALL worm_centroid_move(pint_env, helium, &
                                              iatom, helium%worm_centroid_drmax, ac)
                      ncentratt = ncentratt + 1
                      ncentracc = ncentracc + ac
                   ELSE IF ((imove >= helium%worm_staging_min) .AND. (imove <= helium%worm_staging_max)) THEN
                      ! staging move
-                     iatom = next_random_number(helium%rng_stream_uniform, 1, helium%atoms)
-                     ibead = next_random_number(helium%rng_stream_uniform, 1, helium%beads)
-                     staging_l = next_random_number(helium%rng_stream_uniform, 2, helium%worm_staging_l)
+                     iatom = helium%rng_stream_uniform%next(1, helium%atoms)
+                     ibead = helium%rng_stream_uniform%next(1, helium%beads)
+                     staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
                      CALL worm_staging_move(pint_env, helium, &
                                             iatom, ibead, staging_l, ac)
                      nstagatt = nstagatt + 1
@@ -149,7 +148,7 @@ CONTAINS
                   ELSE IF ((imove >= helium%worm_fcrawl_min) .AND. (imove <= helium%worm_fcrawl_max)) THEN
                      ! crawl forward
                      DO icrawl = 1, helium%worm_repeat_crawl
-                        staging_l = next_random_number(helium%rng_stream_uniform, 2, helium%worm_staging_l)
+                        staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
                         CALL worm_crawl_move_forward(pint_env, helium, &
                                                      staging_l, ac)
                         ncrawlfwdatt = ncrawlfwdatt + 1
@@ -158,7 +157,7 @@ CONTAINS
                   ELSE IF ((imove >= helium%worm_bcrawl_min) .AND. (imove <= helium%worm_bcrawl_max)) THEN
                      ! crawl backward
                      DO icrawl = 1, helium%worm_repeat_crawl
-                        staging_l = next_random_number(helium%rng_stream_uniform, 2, helium%worm_staging_l)
+                        staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
                         CALL worm_crawl_move_backward(pint_env, helium, &
                                                       staging_l, ac)
                         ncrawlbwdatt = ncrawlbwdatt + 1
@@ -166,20 +165,20 @@ CONTAINS
                      END DO
                   ELSE IF ((imove >= helium%worm_head_min) .AND. (imove <= helium%worm_head_max)) THEN
                      ! move head
-                     staging_l = next_random_number(helium%rng_stream_uniform, 2, helium%worm_staging_l)
+                     staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
                      CALL worm_head_move(pint_env, helium, &
                                          staging_l, ac)
                      nmoveheadatt = nmoveheadatt + 1
                      nmoveheadacc = nmoveheadacc + ac
                   ELSE IF ((imove >= helium%worm_tail_min) .AND. (imove <= helium%worm_tail_max)) THEN
                      ! move tail
-                     staging_l = next_random_number(helium%rng_stream_uniform, 2, helium%worm_staging_l)
+                     staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
                      CALL worm_tail_move(pint_env, helium, &
                                          staging_l, ac)
                      nmovetailatt = nmovetailatt + 1
                      nmovetailacc = nmovetailacc + ac
                   ELSE IF ((imove >= helium%worm_swap_min) .AND. (imove <= helium%worm_swap_max)) THEN
-                     staging_l = next_random_number(helium%rng_stream_uniform, 2, helium%worm_staging_l)
+                     staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
                      CALL worm_swap_move(pint_env, helium, &
                                          helium%atoms, staging_l, ac)
                      npswapacc = npswapacc + ac
@@ -187,7 +186,7 @@ CONTAINS
                      nswapatt = nswapatt + 1
                   ELSE IF ((imove >= helium%worm_open_close_min) .AND. (imove <= helium%worm_open_close_max)) THEN
                      ! attempt closing of worm
-                     staging_l = next_random_number(helium%rng_stream_uniform, 2, helium%worm_staging_l)
+                     staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
                      CALL worm_close_move(pint_env, helium, &
                                           staging_l, ac)
                      ncloseatt = ncloseatt + 1
@@ -220,19 +219,19 @@ CONTAINS
          END DO !attempts loop
       ELSE ! only closed configurations allowed
          DO iMC = 1, nMC
-            imove = next_random_number(helium%rng_stream_uniform, 1, helium%worm_all_limit)
+            imove = helium%rng_stream_uniform%next(1, helium%worm_all_limit)
 
             IF ((imove >= helium%worm_centroid_min) .AND. (imove <= helium%worm_centroid_max)) THEN
                ! centroid move
-               iatom = next_random_number(helium%rng_stream_uniform, 1, helium%atoms)
+               iatom = helium%rng_stream_uniform%next(1, helium%atoms)
                CALL worm_centroid_move(pint_env, helium, &
                                        iatom, helium%worm_centroid_drmax, ac)
                ncentratt = ncentratt + 1
                ncentracc = ncentracc + ac
             ELSE IF ((imove >= helium%worm_staging_min) .AND. (imove <= helium%worm_staging_max)) THEN
                ! staging move
-               iatom = next_random_number(helium%rng_stream_uniform, 1, helium%atoms)
-               ibead = next_random_number(helium%rng_stream_uniform, 1, helium%beads)
+               iatom = helium%rng_stream_uniform%next(1, helium%atoms)
+               ibead = helium%rng_stream_uniform%next(1, helium%beads)
                CALL worm_staging_move(pint_env, helium, &
                                       iatom, ibead, helium%worm_staging_l, ac)
                nstagatt = nstagatt + 1
@@ -355,7 +354,7 @@ CONTAINS
       REAL(KIND=dp), DIMENSION(3)                        :: dr, dro, new_com, old_com
 
       DO ic = 1, 3
-         rtmp = next_random_number(helium%rng_stream_uniform)
+         rtmp = helium%rng_stream_uniform%next()
          dr(ic) = (2.0_dp*rtmp - 1.0_dp)*drmax
       END DO
 
@@ -417,7 +416,7 @@ CONTAINS
          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
             should_reject = .TRUE.
          ELSE
-            rtmp = next_random_number(helium%rng_stream_uniform)
+            rtmp = helium%rng_stream_uniform%next()
             IF (EXP(sdiff) < rtmp) THEN
                should_reject = .TRUE.
             END IF
@@ -846,8 +845,6 @@ CONTAINS
 
       partaction = partaction*helium%tau
 
-      RETURN
-
    END FUNCTION worm_centroid_move_inter_action
 
 ! **************************************************************************************************
@@ -886,8 +883,7 @@ CONTAINS
       invstagemass = rk*weight*imass
       ! proposing new positions
       DO idim = 1, 3
-         new_path(idim, 1) = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                                variance=helium%tau*invstagemass)
+         new_path(idim, 1) = helium%rng_stream_gaussian%next(variance=helium%tau*invstagemass)
       END DO
       new_path(:, 1) = new_path(:, 1) + weight*(re(:) + rk*rs(:))
 
@@ -899,8 +895,7 @@ CONTAINS
          invstagemass = rk*weight*imass
          ! proposing new positions
          DO idim = 1, 3
-            new_path(idim, istage) = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                                        variance=helium%tau*invstagemass)
+            new_path(idim, istage) = helium%rng_stream_gaussian%next(variance=helium%tau*invstagemass)
          END DO
          new_path(:, istage) = new_path(:, istage) + weight*(rk*new_path(:, istage - 1) + re(:))
       END DO
@@ -1025,7 +1020,7 @@ CONTAINS
          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
             should_reject = .TRUE.
          ELSE
-            rtmp = next_random_number(helium%rng_stream_uniform)
+            rtmp = helium%rng_stream_uniform%next()
             IF (EXP(sdiff) < rtmp) THEN
                should_reject = .TRUE.
             END IF
@@ -1134,8 +1129,6 @@ CONTAINS
          END DO
       END IF
 
-      RETURN
-
    END SUBROUTINE worm_staging_move
 
 ! **************************************************************************************************
@@ -1196,15 +1189,13 @@ CONTAINS
          ! gro head from startbead
          DO kbead = startbead + 1, endbead - 1
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
             END DO
          END DO
          ! last grow head bead
          DO idim = 1, 3
-            xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                    variance=helium%hb2m*helium%tau)
+            xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
             helium%worm_xtra_bead_work(idim) = helium%work(idim, startatom, endbead - 1) + xr
          END DO
       ELSE IF (endbead /= 1) THEN
@@ -1212,28 +1203,24 @@ CONTAINS
          ! grow from startbead
          DO kbead = startbead + 1, helium%beads
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
             END DO
          END DO
          ! bead one of endatom relative to last on startatom
          DO idim = 1, 3
-            xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                    variance=helium%hb2m*helium%tau)
+            xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
             helium%work(idim, endatom, 1) = helium%work(idim, startatom, helium%beads) + xr
          END DO
          ! everything on endatom
          DO kbead = 2, endbead - 1
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, endatom, kbead) = helium%work(idim, endatom, kbead - 1) + xr
             END DO
          END DO
          DO idim = 1, 3
-            xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                    variance=helium%hb2m*helium%tau)
+            xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
             helium%worm_xtra_bead_work(idim) = helium%work(idim, endatom, endbead - 1) + xr
          END DO
       ELSE ! imagtimewrap and headbead = 1
@@ -1241,15 +1228,13 @@ CONTAINS
          ! grow from startbead
          DO kbead = startbead + 1, helium%beads
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
             END DO
          END DO
          ! bead one of endatom relative to last on startatom
          DO idim = 1, 3
-            xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                    variance=helium%hb2m*helium%tau)
+            xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
             helium%worm_xtra_bead_work(idim) = helium%work(idim, startatom, helium%beads) + xr
          END DO
       END IF
@@ -1280,7 +1265,7 @@ CONTAINS
          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
             should_reject = .TRUE.
          ELSE
-            rtmp = next_random_number(helium%rng_stream_uniform)
+            rtmp = helium%rng_stream_uniform%next()
             IF (EXP(sdiff) < rtmp) THEN
                should_reject = .TRUE.
             END IF
@@ -1520,7 +1505,7 @@ CONTAINS
          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
             should_reject = .TRUE.
          ELSE
-            rtmp = next_random_number(helium%rng_stream_uniform)
+            rtmp = helium%rng_stream_uniform%next()
             IF (EXP(sdiff) < rtmp) THEN
                should_reject = .TRUE.
             END IF
@@ -1701,15 +1686,13 @@ CONTAINS
          ! gro head from startbead
          DO kbead = startbead + 1, endbead - 1
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
             END DO
          END DO
          ! last grow head bead
          DO idim = 1, 3
-            xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                    variance=helium%hb2m*helium%tau)
+            xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
             helium%worm_xtra_bead_work(idim) = helium%work(idim, startatom, endbead - 1) + xr
          END DO
       ELSE IF (endbead /= 1) THEN
@@ -1717,28 +1700,24 @@ CONTAINS
          ! grow from startbead
          DO kbead = startbead + 1, helium%beads
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
             END DO
          END DO
          ! bead one of endatom relative to last on startatom
          DO idim = 1, 3
-            xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                    variance=helium%hb2m*helium%tau)
+            xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
             helium%work(idim, endatom, 1) = helium%work(idim, startatom, helium%beads) + xr
          END DO
          ! everything on endatom
          DO kbead = 2, endbead - 1
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, endatom, kbead) = helium%work(idim, endatom, kbead - 1) + xr
             END DO
          END DO
          DO idim = 1, 3
-            xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                    variance=helium%hb2m*helium%tau)
+            xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
             helium%worm_xtra_bead_work(idim) = helium%work(idim, endatom, endbead - 1) + xr
          END DO
       ELSE ! imagtimewrap and headbead = 1
@@ -1746,15 +1725,13 @@ CONTAINS
          ! grow from startbead
          DO kbead = startbead + 1, helium%beads
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
             END DO
          END DO
          ! bead one of endatom relative to last on startatom
          DO idim = 1, 3
-            xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                    variance=helium%hb2m*helium%tau)
+            xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
             helium%worm_xtra_bead_work(idim) = helium%work(idim, startatom, helium%beads) + xr
          END DO
       END IF
@@ -1778,7 +1755,7 @@ CONTAINS
          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
             should_reject = .TRUE.
          ELSE
-            rtmp = next_random_number(helium%rng_stream_uniform)
+            rtmp = helium%rng_stream_uniform%next()
             IF (EXP(sdiff) < rtmp) THEN
                should_reject = .TRUE.
             END IF
@@ -1961,8 +1938,7 @@ CONTAINS
          ! gro tail from endbead to startbead (confusing eh?)
          DO kbead = endbead - 1, startbead, -1
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead + 1) + xr
             END DO
          END DO
@@ -1971,24 +1947,21 @@ CONTAINS
          ! grow from endbead
          DO kbead = endbead - 1, 1, -1
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, endatom, kbead) = helium%work(idim, endatom, kbead + 1) + xr
             END DO
          END DO
 
          ! over imaginary time boundary
          DO idim = 1, 3
-            xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                    variance=helium%hb2m*helium%tau)
+            xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
             helium%work(idim, startatom, helium%beads) = helium%work(idim, endatom, 1) + xr
          END DO
 
          ! rest on startatom
          DO kbead = helium%beads - 1, startbead, -1
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead + 1) + xr
             END DO
          END DO
@@ -2013,7 +1986,7 @@ CONTAINS
          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
             should_reject = .TRUE.
          ELSE
-            rtmp = next_random_number(helium%rng_stream_uniform)
+            rtmp = helium%rng_stream_uniform%next()
             IF (EXP(sdiff) < rtmp) THEN
                should_reject = .TRUE.
             END IF
@@ -2222,15 +2195,13 @@ CONTAINS
          ! gro head from startbead
          DO kbead = helium%worm_bead_idx + 1, helium%worm_bead_idx_work - 1
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, helium%worm_atom_idx, kbead) = helium%work(idim, helium%worm_atom_idx, kbead - 1) + xr
             END DO
          END DO
          ! last grow head bead
          DO idim = 1, 3
-            xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                    variance=helium%hb2m*helium%tau)
+            xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
             helium%worm_xtra_bead_work(idim) = helium%work(idim, helium%worm_atom_idx, helium%worm_bead_idx_work - 1) + xr
          END DO
       ELSE IF (helium%worm_bead_idx_work /= 1) THEN
@@ -2238,28 +2209,24 @@ CONTAINS
          ! grow from startbead
          DO kbead = helium%worm_bead_idx + 1, helium%beads
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, helium%worm_atom_idx, kbead) = helium%work(idim, helium%worm_atom_idx, kbead - 1) + xr
             END DO
          END DO
          ! bead one of endatom relative to last on helium%worm_atom_idx
          DO idim = 1, 3
-            xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                    variance=helium%hb2m*helium%tau)
+            xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
             helium%work(idim, helium%worm_atom_idx_work, 1) = helium%work(idim, helium%worm_atom_idx, helium%beads) + xr
          END DO
          ! everything on endatom
          DO kbead = 2, helium%worm_bead_idx_work - 1
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, helium%worm_atom_idx_work, kbead) = helium%work(idim, helium%worm_atom_idx_work, kbead - 1) + xr
             END DO
          END DO
          DO idim = 1, 3
-            xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                    variance=helium%hb2m*helium%tau)
+            xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
             helium%worm_xtra_bead_work(idim) = helium%work(idim, helium%worm_atom_idx_work, helium%worm_bead_idx_work - 1) + xr
          END DO
       ELSE ! imagtimewrap and headbead = 1
@@ -2267,15 +2234,13 @@ CONTAINS
          ! grow from startbead
          DO kbead = helium%worm_bead_idx + 1, helium%beads
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, helium%worm_atom_idx, kbead) = helium%work(idim, helium%worm_atom_idx, kbead - 1) + xr
             END DO
          END DO
          ! bead one of endatom relative to last on helium%worm_atom_idx
          DO idim = 1, 3
-            xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                    variance=helium%hb2m*helium%tau)
+            xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
             helium%worm_xtra_bead_work(idim) = helium%work(idim, helium%worm_atom_idx, helium%beads) + xr
          END DO
       END IF
@@ -2301,7 +2266,7 @@ CONTAINS
          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
             should_reject = .TRUE.
          ELSE
-            rtmp = next_random_number(helium%rng_stream_uniform)
+            rtmp = helium%rng_stream_uniform%next()
             IF (EXP(sdiff) < rtmp) THEN
                should_reject = .TRUE.
             END IF
@@ -2519,8 +2484,7 @@ CONTAINS
          ! gro tail from endbead to startbead (confusing eh?)
          DO kbead = helium%worm_bead_idx - 1, helium%worm_bead_idx_work, -1
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, helium%worm_atom_idx, kbead) = helium%work(idim, helium%worm_atom_idx, kbead + 1) + xr
             END DO
 
@@ -2530,24 +2494,21 @@ CONTAINS
          ! grow from endbead
          DO kbead = helium%worm_bead_idx - 1, 1, -1
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, helium%worm_atom_idx, kbead) = helium%work(idim, helium%worm_atom_idx, kbead + 1) + xr
             END DO
          END DO
 
          ! over imaginary time boundary
          DO idim = 1, 3
-            xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                    variance=helium%hb2m*helium%tau)
+            xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
             helium%work(idim, helium%worm_atom_idx_work, helium%beads) = helium%work(idim, helium%worm_atom_idx, 1) + xr
          END DO
 
          ! rest on startatom
          DO kbead = helium%beads - 1, helium%worm_bead_idx_work, -1
             DO idim = 1, 3
-               xr = next_random_number(rng_stream=helium%rng_stream_gaussian, &
-                                       variance=helium%hb2m*helium%tau)
+               xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
                helium%work(idim, helium%worm_atom_idx_work, kbead) = helium%work(idim, helium%worm_atom_idx_work, kbead + 1) + xr
             END DO
          END DO
@@ -2573,7 +2534,7 @@ CONTAINS
          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
             should_reject = .TRUE.
          ELSE
-            rtmp = next_random_number(helium%rng_stream_uniform)
+            rtmp = helium%rng_stream_uniform%next()
             IF (EXP(sdiff) < rtmp) THEN
                should_reject = .TRUE.
             END IF
@@ -2795,7 +2756,7 @@ CONTAINS
       forwarddensmatsum = SUM(forwarddensmat)
 
       ! Select an atom with its corresponding probability
-      rtmp = next_random_number(helium%rng_stream_uniform)*forwarddensmatsum
+      rtmp = helium%rng_stream_uniform%next()*forwarddensmatsum
       fendatom = 1
       DO WHILE (rtmp >= forwarddensmat(fendatom))
          rtmp = rtmp - forwarddensmat(fendatom)
@@ -2941,7 +2902,7 @@ CONTAINS
          IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
             should_reject = .TRUE.
          ELSE
-            rtmp = next_random_number(helium%rng_stream_uniform)
+            rtmp = helium%rng_stream_uniform%next()
             IF (EXP(sdiff) < rtmp) THEN
                should_reject = .TRUE.
             END IF
diff --git a/src/motion/input_cp2k_restarts.F b/src/motion/input_cp2k_restarts.F
index 4de232dcc3..9ce782002c 100644
--- a/src/motion/input_cp2k_restarts.F
+++ b/src/motion/input_cp2k_restarts.F
@@ -64,9 +64,7 @@ MODULE input_cp2k_restarts
    USE molecule_kind_list_types,        ONLY: molecule_kind_list_type
    USE molecule_list_types,             ONLY: molecule_list_type
    USE neb_types,                       ONLY: neb_var_type
-   USE parallel_rng_types,              ONLY: dump_rng_stream,&
-                                              get_rng_stream,&
-                                              rng_record_length
+   USE parallel_rng_types,              ONLY: rng_record_length
    USE particle_list_types,             ONLY: particle_list_type
    USE particle_types,                  ONLY: get_particle_pos_or_vel,&
                                               particle_type
@@ -1008,8 +1006,7 @@ CONTAINS
       ELSE IF (pint_env%pimd_thermostat == thermostat_pile) THEN
          tmpsec => section_vals_get_subs_vals(pint_section, &
                                               "PILE%RNG_INIT")
-         CALL dump_rng_stream(rng_stream=pint_env%pile_therm%gaussian_rng_stream, &
-                              rng_record=rng_record)
+         CALL pint_env%pile_therm%gaussian_rng_stream%dump(rng_record)
          CALL string_to_ascii(rng_record, ascii(:, 1))
          CALL section_rng_val_set(rng_section=tmpsec, nsize=1, &
                                   ascii=ascii)
@@ -1023,8 +1020,7 @@ CONTAINS
       ELSE IF (pint_env%pimd_thermostat == thermostat_piglet) THEN
          tmpsec => section_vals_get_subs_vals(pint_section, &
                                               "PIGLET%RNG_INIT")
-         CALL dump_rng_stream(rng_stream=pint_env%piglet_therm%gaussian_rng_stream, &
-                              rng_record=rng_record)
+         CALL pint_env%piglet_therm%gaussian_rng_stream%dump(rng_record)
          CALL string_to_ascii(rng_record, ascii(:, 1))
          CALL section_rng_val_set(rng_section=tmpsec, nsize=1, &
                                   ascii=ascii)
@@ -1261,8 +1257,8 @@ CONTAINS
          real_msg_gather(:) = 0.0_dp
 
          DO i = 1, SIZE(helium_env)
-            CALL get_rng_stream(helium_env(i)%helium%rng_stream_uniform, bg=bg, cg=cg, ig=ig, &
-                                buffer=bu, buffer_filled=lbf)
+            CALL helium_env(i)%helium%rng_stream_uniform%get(bg=bg, cg=cg, ig=ig, &
+                                                             buffer=bu, buffer_filled=lbf)
             off = 0
             real_msg(off + 1:off + 6) = PACK(bg, .TRUE.)
             real_msg(off + 7:off + 12) = PACK(cg, .TRUE.)
@@ -1274,8 +1270,8 @@ CONTAINS
             END IF
             real_msg(off + 19) = bf
             real_msg(off + 20) = bu
-            CALL get_rng_stream(helium_env(i)%helium%rng_stream_gaussian, bg=bg, cg=cg, ig=ig, &
-                                buffer=bu, buffer_filled=lbf)
+            CALL helium_env(i)%helium%rng_stream_gaussian%get(bg=bg, cg=cg, ig=ig, &
+                                                              buffer=bu, buffer_filled=lbf)
             off = 20
             real_msg(off + 1:off + 6) = PACK(bg, .TRUE.)
             real_msg(off + 7:off + 12) = PACK(cg, .TRUE.)
@@ -1553,8 +1549,7 @@ CONTAINS
       dwork = 0
       DO i = 1, csvr%loc_num_csvr
          my_index = csvr%map_info%index(i)
-         CALL dump_rng_stream(rng_stream=csvr%nvt(i)%gaussian_rng_stream, &
-                              rng_record=rng_record)
+         CALL csvr%nvt(i)%gaussian_rng_stream%dump(rng_record)
          CALL string_to_ascii(rng_record, dwork(:, my_index))
       END DO
 
@@ -1681,8 +1676,7 @@ CONTAINS
       dwork = 0
       DO i = 1, loc_num
          j = gle%map_info%index(i)
-         CALL dump_rng_stream(rng_stream=gle%nvt(i)%gaussian_rng_stream, &
-                              rng_record=rng_record)
+         CALL gle%nvt(i)%gaussian_rng_stream%dump(rng_record)
          CALL string_to_ascii(rng_record, dwork(:, j))
       END DO
 
diff --git a/src/motion/integrator.F b/src/motion/integrator.F
index c08a20ecc1..2b4f6487ed 100644
--- a/src/motion/integrator.F
+++ b/src/motion/integrator.F
@@ -73,8 +73,6 @@ MODULE integrator
    USE molecule_list_types,             ONLY: molecule_list_type
    USE molecule_types,                  ONLY: global_constraint_type,&
                                               molecule_type
-   USE parallel_rng_types,              ONLY: next_random_number,&
-                                              rng_stream_type
    USE particle_list_types,             ONLY: particle_list_type
    USE particle_types,                  ONLY: particle_type,&
                                               update_particle_set
@@ -154,7 +152,6 @@ CONTAINS
       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
       TYPE(particle_list_type), POINTER                  :: particles
       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
       TYPE(simpar_type), POINTER                         :: simpar
       TYPE(thermal_region_type), POINTER                 :: thermal_region
       TYPE(thermal_regions_type), POINTER                :: thermal_regions
@@ -162,7 +159,7 @@ CONTAINS
 
       NULLIFY (cell, para_env, gci, force_env)
       NULLIFY (atomic_kinds, local_particles, subsys, local_molecules, molecule_kinds, molecules)
-      NULLIFY (molecule_kind_set, molecule_set, particles, particle_set, rng_stream, simpar, virial)
+      NULLIFY (molecule_kind_set, molecule_set, particles, particle_set, simpar, virial)
       NULLIFY (thermal_region, thermal_regions, itimes)
 
       CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, &
@@ -253,11 +250,12 @@ CONTAINS
             iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
             IF (do_langevin(iparticle)) THEN
                sigma = var_w(iparticle)*mass
-               rng_stream => local_particles%local_particle_set(iparticle_kind)% &
-                             rng(iparticle_local)%stream
-               w(1, iparticle) = next_random_number(rng_stream, variance=sigma)
-               w(2, iparticle) = next_random_number(rng_stream, variance=sigma)
-               w(3, iparticle) = next_random_number(rng_stream, variance=sigma)
+               ASSOCIATE (rng_stream=>local_particles%local_particle_set(iparticle_kind)% &
+                          rng(iparticle_local))
+                  w(1, iparticle) = rng_stream%stream%next(variance=sigma)
+                  w(2, iparticle) = rng_stream%stream%next(variance=sigma)
+                  w(3, iparticle) = rng_stream%stream%next(variance=sigma)
+               END ASSOCIATE
             END IF
          END DO
       END DO
@@ -744,7 +742,6 @@ CONTAINS
                                                             shell_particles
       TYPE(particle_type), DIMENSION(:), POINTER         :: core_particle_set, particle_set, &
                                                             shell_particle_set
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
       TYPE(simpar_type), POINTER                         :: simpar
       TYPE(thermostat_type), POINTER                     :: thermostat_coeff, thermostat_fast, &
                                                             thermostat_shell, thermostat_slow
@@ -831,8 +828,7 @@ CONTAINS
       IF (ASSOCIATED(force_env%meta_env)) THEN
          IF (force_env%meta_env%langevin) THEN
             DO ivar = 1, force_env%meta_env%n_colvar
-               rng_stream => force_env%meta_env%rng(ivar)%stream
-               rand(ivar) = next_random_number(rng_stream)
+               rand(ivar) = force_env%meta_env%rng(ivar)%next()
             ENDDO
             CALL metadyn_velocities_colvar(force_env, rand)
          ENDIF
@@ -956,7 +952,6 @@ CONTAINS
                                                             shell_particles
       TYPE(particle_type), DIMENSION(:), POINTER         :: core_particle_set, particle_set, &
                                                             shell_particle_set
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
       TYPE(simpar_type), POINTER                         :: simpar
       TYPE(thermostat_type), POINTER                     :: thermostat_coeff, thermostat_part, &
                                                             thermostat_shell
@@ -1040,8 +1035,7 @@ CONTAINS
       IF (ASSOCIATED(force_env%meta_env)) THEN
          IF (force_env%meta_env%langevin) THEN
             DO ivar = 1, force_env%meta_env%n_colvar
-               rng_stream => force_env%meta_env%rng(ivar)%stream
-               rand(ivar) = next_random_number(rng_stream)
+               rand(ivar) = force_env%meta_env%rng(ivar)%next()
             ENDDO
             CALL metadyn_velocities_colvar(force_env, rand)
          ENDIF
@@ -1208,7 +1202,6 @@ CONTAINS
                                                             shell_particles
       TYPE(particle_type), DIMENSION(:), POINTER         :: core_particle_set, particle_set, &
                                                             shell_particle_set
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
       TYPE(simpar_type), POINTER                         :: simpar
       TYPE(thermostat_type), POINTER                     :: thermostat_baro, thermostat_part, &
                                                             thermostat_shell
@@ -1324,8 +1317,7 @@ CONTAINS
       IF (ASSOCIATED(force_env%meta_env)) THEN
          IF (force_env%meta_env%langevin) THEN
             DO ivar = 1, force_env%meta_env%n_colvar
-               rng_stream => force_env%meta_env%rng(ivar)%stream
-               rand(ivar) = next_random_number(rng_stream)
+               rand(ivar) = force_env%meta_env%rng(ivar)%next()
             ENDDO
             CALL metadyn_velocities_colvar(force_env, rand)
          ENDIF
diff --git a/src/motion/mc/mc_control.F b/src/motion/mc/mc_control.F
index aa1522ab7b..4b1e9a3504 100644
--- a/src/motion/mc/mc_control.F
+++ b/src/motion/mc/mc_control.F
@@ -45,8 +45,7 @@ MODULE mc_control
    USE molecule_kind_types,             ONLY: atom_type,&
                                               get_molecule_kind,&
                                               molecule_kind_type
-   USE parallel_rng_types,              ONLY: next_random_number,&
-                                              rng_stream_type
+   USE parallel_rng_types,              ONLY: rng_stream_type
    USE particle_list_types,             ONLY: particle_list_type
    USE physcon,                         ONLY: angstrom
 #include "../../base/base_uses.f90"
@@ -172,7 +171,7 @@ CONTAINS
       TYPE(force_env_type), POINTER                      :: force_env
       INTEGER, INTENT(IN)                                :: iw
       INTEGER, INTENT(INOUT)                             :: mc_nunits_tot
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(len=*), PARAMETER :: routineN = 'read_mc_restart', &
          routineP = moduleN//':'//routineN
@@ -318,7 +317,7 @@ CONTAINS
       ! advance the random number sequence based on the restart step
       IF (ionode) THEN
          DO i = 1, nstart + 1
-            rand = next_random_number(rng_stream)
+            rand = rng_stream%next()
          ENDDO
       ENDIF
 
@@ -420,4 +419,3 @@ CONTAINS
    END SUBROUTINE mc_create_bias_force_env
 
 END MODULE mc_control
-
diff --git a/src/motion/mc/mc_coordinates.F b/src/motion/mc/mc_coordinates.F
index b2da2102b1..732c80efdd 100644
--- a/src/motion/mc/mc_coordinates.F
+++ b/src/motion/mc/mc_coordinates.F
@@ -24,8 +24,7 @@ MODULE mc_coordinates
                                               mc_simpar_type
    USE message_passing,                 ONLY: mp_bcast
    USE molecule_types,                  ONLY: molecule_type
-   USE parallel_rng_types,              ONLY: next_random_number,&
-                                              rng_stream_type
+   USE parallel_rng_types,              ONLY: rng_stream_type
    USE particle_list_types,             ONLY: particle_list_type
    USE physcon,                         ONLY: angstrom
 #include "../../base/base_uses.f90"
@@ -345,7 +344,7 @@ CONTAINS
       LOGICAL, INTENT(IN)                                :: ionode, lremove
       INTEGER, DIMENSION(:), INTENT(IN)                  :: mol_type, nchains
       INTEGER, INTENT(IN)                                :: source, group
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
       INTEGER, INTENT(IN), OPTIONAL                      :: avbmc_atom
       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: rmin, rmax
       CHARACTER(LEN=*), INTENT(IN), OPTIONAL             :: move_type
@@ -463,7 +462,7 @@ CONTAINS
             ELSE
 ! find a new insertion point somewhere in the box
                DO i = 1, 3
-                  rand = next_random_number(rng_stream)
+                  rand = rng_stream%next()
                   r_insert(i) = rand*abc(i)
                ENDDO
 
@@ -548,7 +547,7 @@ CONTAINS
       total_running_weight = 0.0E0_dp
       choosen = 0
       IF (ionode) THEN
-         rand = next_random_number(rng_stream)
+         rand = rng_stream%next()
 !         CALL random_number(rand)
       ENDIF
       CALL mp_bcast(rand, source, group)
@@ -631,7 +630,7 @@ CONTAINS
       REAL(KIND=dp), DIMENSION(1:natoms), INTENT(IN)     :: mass
       REAL(KIND=dp), DIMENSION(1:3, 1:natoms), &
          INTENT(INOUT)                                   :: r
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(len=*), PARAMETER :: routineN = 'rotate_molecule', &
          routineP = moduleN//':'//routineN
@@ -649,7 +648,7 @@ CONTAINS
       CALL get_center_of_mass(r(:, :), natoms, center_of_mass(:), mass(:))
 
 ! call a random number to figure out how far we're moving
-      rand = next_random_number(rng_stream)
+      rand = rng_stream%next()
       dgamma = pi*(rand - 0.5E0_dp)*2.0E0_dp
 
 ! *** set up the rotation matrix ***
@@ -719,7 +718,7 @@ CONTAINS
 
       TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info
       INTEGER, INTENT(OUT)                               :: start_atom, box_number, molecule_type
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
       INTEGER, INTENT(IN), OPTIONAL                      :: box, molecule_type_old
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'find_mc_test_molecule', &
@@ -747,7 +746,7 @@ CONTAINS
 
       IF (PRESENT(box) .AND. PRESENT(molecule_type_old)) THEN
 ! only need to find the atom number the molecule starts on
-         rand = next_random_number(rng_stream)
+         rand = rng_stream%next()
          molecule_number = CEILING(rand*REAL(nchains(molecule_type_old, box), KIND=dp))
 
          start_mol = 1
@@ -767,7 +766,7 @@ CONTAINS
 
       ELSEIF (PRESENT(box)) THEN
 ! any molecule in box...need to find molecule type and start atom
-         rand = next_random_number(rng_stream)
+         rand = rng_stream%next()
          molecule_number = CEILING(rand*REAL(SUM(nchains(:, box)), KIND=dp))
 
          start_mol = 1
@@ -785,7 +784,7 @@ CONTAINS
 
       ELSEIF (PRESENT(molecule_type_old)) THEN
 ! any molecule of type molecule_type_old...need to find box number and start atom
-         rand = next_random_number(rng_stream)
+         rand = rng_stream%next()
          molecule_number = CEILING(rand*REAL(SUM(nchains(molecule_type_old, :)), KIND=dp))
 
 ! find which box it's in
@@ -818,7 +817,7 @@ CONTAINS
          DO ibox = 1, SIZE(nchains(1, :))
             nchains_tot = nchains_tot + SUM(nchains(:, ibox))
          ENDDO
-         rand = next_random_number(rng_stream)
+         rand = rng_stream%next()
          molecule_number = CEILING(rand*REAL(nchains_tot, KIND=dp))
 
          molecule_type = mol_type(molecule_number)
@@ -939,7 +938,7 @@ CONTAINS
       CHARACTER(LEN=*), INTENT(IN)                       :: move_type
       REAL(KIND=dp), DIMENSION(1:3), INTENT(OUT)         :: r_insert
       REAL(KIND=dp), DIMENSION(1:3), INTENT(IN)          :: abc
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       INTEGER                                            :: i
       REAL(dp)                                           :: dist, eta_1, eta_2, eta_sq, rand
@@ -950,8 +949,8 @@ CONTAINS
       IF (move_type == 'in') THEN
 ! generate a random unit vector, from Allen and Tildesley
          DO
-            eta_1 = next_random_number(rng_stream)
-            eta_2 = next_random_number(rng_stream)
+            eta_1 = rng_stream%next()
+            eta_2 = rng_stream%next()
             eta_sq = eta_1**2 + eta_2**2
             IF (eta_sq .LT. 1.0_dp) THEN
                r_insert(1) = 2.0_dp*eta_1*SQRT(1.0_dp - eta_sq)
@@ -962,7 +961,7 @@ CONTAINS
          ENDDO
 
 ! now scale that vector to be within the "in" region
-         rand = next_random_number(rng_stream)
+         rand = rng_stream%next()
          r_insert(1:3) = r_insert(1:3)*(rand*(rmax**3 - rmin**3) + rmin**3)** &
                          (1.0_dp/3.0_dp)
 
@@ -972,7 +971,7 @@ CONTAINS
 ! find a new insertion point somewhere in the box
          DO
             DO i = 1, 3
-               rand = next_random_number(rng_stream)
+               rand = rng_stream%next()
                r_insert(i) = rand*abc(i)
             ENDDO
 
diff --git a/src/motion/mc/mc_ensembles.F b/src/motion/mc/mc_ensembles.F
index 42cc153303..a6b0e5384a 100644
--- a/src/motion/mc/mc_ensembles.F
+++ b/src/motion/mc/mc_ensembles.F
@@ -73,8 +73,7 @@ MODULE mc_ensembles
                                               mc_simulation_parameters_p_type,&
                                               set_mc_par
    USE message_passing,                 ONLY: mp_bcast
-   USE parallel_rng_types,              ONLY: next_random_number,&
-                                              rng_stream_type
+   USE parallel_rng_types,              ONLY: rng_stream_type
    USE particle_list_types,             ONLY: particle_list_p_type,&
                                               particle_list_type
    USE particle_methods,                ONLY: write_particle_coordinates
@@ -117,7 +116,7 @@ CONTAINS
       TYPE(global_environment_type), POINTER             :: globenv
       TYPE(section_type), POINTER                        :: input_declaration
       INTEGER, INTENT(IN)                                :: nboxes
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(len=*), PARAMETER :: routineN = 'mc_run_ensemble', &
          routineP = moduleN//':'//routineN
@@ -444,7 +443,7 @@ CONTAINS
             WRITE (iw, *) "------- On Monte Carlo Step ", nnstep
          ENDIF
 
-         IF (ionode) rand = next_random_number(rng_stream)
+         IF (ionode) rand = rng_stream%next()
          ! broadcast the random number, to make sure we're on the same move
          CALL mp_bcast(rand, source, group)
 
@@ -469,7 +468,7 @@ CONTAINS
                                       r_old, rng_stream)
             CASE ("GEMC_NPT")
                ! we need to select a box based on the probability given in the input file
-               IF (ionode) rand = next_random_number(rng_stream)
+               IF (ionode) rand = rng_stream%next()
                CALL mp_bcast(rand, source, group)
 
                DO ibox = 1, nboxes
@@ -574,7 +573,7 @@ CONTAINS
             ENDIF
 
 ! pick a box at random
-            IF (ionode) rand = next_random_number(rng_stream)
+            IF (ionode) rand = rng_stream%next()
             CALL mp_bcast(rand, source, group)
 
             DO ibox = 1, nboxes
@@ -600,7 +599,7 @@ CONTAINS
             ENDIF
 
             ! first, pick a box to do it for
-            IF (ionode) rand = next_random_number(rng_stream)
+            IF (ionode) rand = rng_stream%next()
             CALL mp_bcast(rand, source, group)
 
             IF (nboxes .EQ. 2) THEN
@@ -614,7 +613,7 @@ CONTAINS
             ENDIF
 
             ! now pick a molecule type to do it for
-            IF (ionode) rand = next_random_number(rng_stream)
+            IF (ionode) rand = rng_stream%next()
             CALL mp_bcast(rand, source, group)
             molecule_type_swap = 0
             DO imol_type = 1, nmol_types
@@ -654,7 +653,7 @@ CONTAINS
 
                   ! choose if we're swapping into the bonded region of mol_target, or
                   ! into the nonbonded region
-                  rand = next_random_number(rng_stream)
+                  rand = rng_stream%next()
 
                ENDIF
                CALL mp_bcast(start_atom_swap, source, group)
@@ -690,12 +689,12 @@ CONTAINS
 
             DO imove = 1, nmoves
 
-               IF (ionode) rand = next_random_number(rng_stream)
+               IF (ionode) rand = rng_stream%next()
                CALL mp_bcast(rand, source, group)
                IF (rand .LT. pmtraion) THEN
                   ! change molecular conformation
                   ! first, pick a box to do it for
-                  IF (ionode) rand = next_random_number(rng_stream)
+                  IF (ionode) rand = rng_stream%next()
                   CALL mp_bcast(rand, source, group)
                   IF (nboxes .EQ. 2) THEN
                      IF (rand .LT. 0.75E0_dp) THEN
@@ -708,7 +707,7 @@ CONTAINS
                   ENDIF
 
                   ! figure out which molecule type we're looking for
-                  IF (ionode) rand = next_random_number(rng_stream)
+                  IF (ionode) rand = rng_stream%next()
                   CALL mp_bcast(rand, source, group)
                   molecule_type = 0
                   DO imol_type = 1, nmol_types
@@ -736,7 +735,7 @@ CONTAINS
                                                    box=box_number, molecule_type_old=molecule_type)
 
                         ! choose if we're changing a bond length or an angle
-                        rand = next_random_number(rng_stream)
+                        rand = rng_stream%next()
                      ENDIF
                      CALL mp_bcast(rand, source, group)
                      CALL mp_bcast(start_atom, source, group)
@@ -766,7 +765,7 @@ CONTAINS
                ELSEIF (rand .LT. pmtrans) THEN
                   ! translate a whole molecule in the system
                   ! pick a molecule type
-                  IF (ionode) rand = next_random_number(rng_stream)
+                  IF (ionode) rand = rng_stream%next()
                   CALL mp_bcast(rand, source, group)
                   molecule_type = 0
                   DO imol_type = 1, nmol_types
@@ -798,7 +797,7 @@ CONTAINS
                ELSEIF (rand .LT. pmcltrans) THEN
                   ! translate a whole cluster in the system
                   ! first, pick a box to do it for
-                  IF (ionode) rand = next_random_number(rng_stream)
+                  IF (ionode) rand = rng_stream%next()
                   CALL mp_bcast(rand, source, group)
 
                   DO ibox = 1, nboxes
@@ -819,7 +818,7 @@ CONTAINS
                ELSE
                   !     rotate a whole molecule in the system
                   ! pick a molecule type
-                  IF (ionode) rand = next_random_number(rng_stream)
+                  IF (ionode) rand = rng_stream%next()
                   CALL mp_bcast(rand, source, group)
                   molecule_type = 0
                   DO imol_type = 1, nmol_types
@@ -1174,7 +1173,7 @@ CONTAINS
    SUBROUTINE mc_compute_virial(mc_env, rng_stream)
 
       TYPE(mc_environment_p_type), DIMENSION(:), POINTER :: mc_env
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(len=*), PARAMETER :: routineN = 'mc_compute_virial', &
          routineP = moduleN//':'//routineN
diff --git a/src/motion/mc/mc_ge_moves.F b/src/motion/mc/mc_ge_moves.F
index bd5cbdcd31..3b57927865 100644
--- a/src/motion/mc/mc_ge_moves.F
+++ b/src/motion/mc/mc_ge_moves.F
@@ -45,8 +45,7 @@ MODULE mc_ge_moves
         mc_molecule_info_destroy, mc_molecule_info_type, mc_moves_p_type, &
         mc_simulation_parameters_p_type, set_mc_par
    USE message_passing,                 ONLY: mp_bcast
-   USE parallel_rng_types,              ONLY: next_random_number,&
-                                              rng_stream_type
+   USE parallel_rng_types,              ONLY: rng_stream_type
    USE particle_list_types,             ONLY: particle_list_p_type,&
                                               particle_list_type
    USE particle_methods,                ONLY: write_particle_coordinates
@@ -117,7 +116,7 @@ CONTAINS
       INTEGER, DIMENSION(:), INTENT(IN)                  :: box_flag
       TYPE(cp_subsys_p_type), DIMENSION(:), POINTER      :: subsys
       TYPE(particle_list_p_type), DIMENSION(:), POINTER  :: particles
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
       REAL(KIND=dp), INTENT(IN)                          :: unit_conv
 
       CHARACTER(len=*), PARAMETER :: routineN = 'mc_Quickstep_move', &
@@ -284,7 +283,7 @@ CONTAINS
          w = 1.0E0_dp
          rand = 0.0E0_dp
       ELSE
-         IF (ionode) rand = next_random_number(rng_stream)
+         IF (ionode) rand = rng_stream%next()
          CALL mp_bcast(rand, source, group)
       ENDIF
 
@@ -435,7 +434,7 @@ CONTAINS
       TYPE(section_type), POINTER                        :: input_declaration
       TYPE(cp_para_env_type), POINTER                    :: para_env
       REAL(KIND=dp), DIMENSION(1:2), INTENT(INOUT)       :: bias_energy_old, last_bias_energy
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(len=*), PARAMETER :: routineN = 'mc_ge_swap_move', &
          routineP = moduleN//':'//routineN
@@ -511,7 +510,7 @@ CONTAINS
       ENDDO
 
 !     choose a direction to swap
-      IF (ionode) rand = next_random_number(rng_stream)
+      IF (ionode) rand = rng_stream%next()
       CALL mp_bcast(rand, source, group)
 
       IF (rand .LE. 0.50E0_dp) THEN
@@ -531,7 +530,7 @@ CONTAINS
 !      insert_box=1
 
 ! now choose a molecule type at random
-      IF (ionode) rand = next_random_number(rng_stream)
+      IF (ionode) rand = rng_stream%next()
       CALL mp_bcast(rand, source, group)
       DO itype = 1, nmol_types
          IF (rand .LT. pmswap_mol(itype)) THEN
@@ -552,7 +551,7 @@ CONTAINS
             moves(molecule_type, insert_box)%moves%empty + 1
       ELSE
 
-         IF (ionode) rand = next_random_number(rng_stream)
+         IF (ionode) rand = rng_stream%next()
          CALL mp_bcast(rand, source, group)
          imolecule = CEILING(rand*nchains(molecule_type, remove_box))
 ! figure out the atom number this molecule starts on
@@ -627,7 +626,7 @@ CONTAINS
       IF (ionode) THEN
 ! choose an insertion point
          DO idim = 1, 3
-            rand = next_random_number(rng_stream)
+            rand = rng_stream%next()
             pos_insert(idim) = rand*abc_insert(idim)
          ENDDO
       ENDIF
@@ -1008,7 +1007,7 @@ CONTAINS
       IF (w .GE. 1.0E0_dp) THEN
          rand = 0.0E0_dp
       ELSE
-         IF (ionode) rand = next_random_number(rng_stream)
+         IF (ionode) rand = rng_stream%next()
          CALL mp_bcast(rand, source, group)
       ENDIF
 
@@ -1138,7 +1137,7 @@ CONTAINS
       INTEGER, INTENT(IN)                                :: nnstep
       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: old_energy, energy_check
       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: r_old
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(len=*), PARAMETER :: routineN = 'mc_ge_volume_move', &
          routineP = moduleN//':'//routineN
@@ -1237,7 +1236,7 @@ CONTAINS
       ENDDO
 
 ! call a random number to figure out how far we're moving
-      IF (ionode) rand = next_random_number(rng_stream)
+      IF (ionode) rand = rng_stream%next()
       CALL mp_bcast(rand, source, group)
 
       vol_dis = rmvolume*(rand - 0.5E0_dp)*2.0E0_dp
@@ -1363,7 +1362,7 @@ CONTAINS
          w = 1.0E0_dp
          rand = 0.0E0_dp
       ELSE
-         IF (ionode) rand = next_random_number(rng_stream)
+         IF (ionode) rand = rng_stream%next()
          CALL mp_bcast(rand, source, group)
       ENDIF
 
diff --git a/src/motion/mc/mc_moves.F b/src/motion/mc/mc_moves.F
index 07f680a108..230b179a46 100644
--- a/src/motion/mc/mc_moves.F
+++ b/src/motion/mc/mc_moves.F
@@ -51,8 +51,7 @@ MODULE mc_moves
                                               get_molecule_kind,&
                                               molecule_kind_type,&
                                               torsion_type
-   USE parallel_rng_types,              ONLY: next_random_number,&
-                                              rng_stream_type
+   USE parallel_rng_types,              ONLY: rng_stream_type
    USE particle_list_types,             ONLY: particle_list_type
    USE physcon,                         ONLY: angstrom
 #include "../../base/base_uses.f90"
@@ -144,7 +143,7 @@ CONTAINS
       REAL(KIND=dp), INTENT(INOUT)                       :: bias_energy
       CHARACTER(LEN=*), INTENT(IN)                       :: move_type
       LOGICAL, INTENT(OUT)                               :: lreject
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(len=*), PARAMETER :: routineN = 'mc_conformation_change', &
          routineP = moduleN//':'//routineN
@@ -350,7 +349,7 @@ CONTAINS
             rand = 0.0E0_dp
          ELSE
             IF (ionode) THEN
-               rand = next_random_number(rng_stream)
+               rand = rng_stream%next()
             ENDIF
             CALL mp_bcast(rand, source, group)
          ENDIF
@@ -440,7 +439,7 @@ CONTAINS
       REAL(KIND=dp), INTENT(INOUT)                       :: bias_energy
       INTEGER, INTENT(IN)                                :: molecule_type
       LOGICAL, INTENT(OUT)                               :: lreject
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(len=*), PARAMETER :: routineN = 'mc_molecule_translation', &
          routineP = moduleN//':'//routineN
@@ -536,13 +535,13 @@ CONTAINS
 ! move one molecule in the system
 
 ! call a random number to figure out which direction we're moving
-      IF (ionode) rand = next_random_number(rng_stream)
+      IF (ionode) rand = rng_stream%next()
       CALL mp_bcast(rand, source, group)
       ! 1,2,3 with equal prob
       move_direction = INT(3*rand) + 1
 
 ! call a random number to figure out how far we're moving
-      IF (ionode) rand = next_random_number(rng_stream)
+      IF (ionode) rand = rng_stream%next()
       CALL mp_bcast(rand, source, group)
       dis_mol = rmtrans(molecule_type)*(rand - 0.5E0_dp)*2.0E0_dp
 
@@ -592,7 +591,7 @@ CONTAINS
             w = 1.0E0_dp
             rand = 0.0E0_dp
          ELSE
-            IF (ionode) rand = next_random_number(rng_stream)
+            IF (ionode) rand = rng_stream%next()
             CALL mp_bcast(rand, source, group)
          ENDIF
 
@@ -663,7 +662,7 @@ CONTAINS
       INTEGER, INTENT(IN)                                :: box_number, start_atom, molecule_type
       REAL(KIND=dp), INTENT(INOUT)                       :: bias_energy
       LOGICAL, INTENT(OUT)                               :: lreject
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(len=*), PARAMETER :: routineN = 'mc_molecule_rotation', &
          routineP = moduleN//':'//routineN
@@ -771,7 +770,7 @@ CONTAINS
 ! rotate one molecule in the system
 
 ! call a random number to figure out which direction we're moving
-      IF (ionode) rand = next_random_number(rng_stream)
+      IF (ionode) rand = rng_stream%next()
 !      CALL RANDOM_NUMBER(rand)
       CALL mp_bcast(rand, source, group)
       ! 1,2,3 with equal prob
@@ -798,7 +797,7 @@ CONTAINS
       nzcm = nzcm/masstot
 
 ! call a random number to figure out how far we're moving
-      IF (ionode) rand = next_random_number(rng_stream)
+      IF (ionode) rand = rng_stream%next()
       CALL mp_bcast(rand, source, group)
       dgamma = rmrot(molecule_type)*(rand - 0.5E0_dp)*2.0E0_dp
 
@@ -895,7 +894,7 @@ CONTAINS
             w = 1.0E0_dp
             rand = 0.0E0_dp
          ELSE
-            IF (ionode) rand = next_random_number(rng_stream)
+            IF (ionode) rand = rng_stream%next()
             CALL mp_bcast(rand, source, group)
          ENDIF
 
@@ -968,7 +967,7 @@ CONTAINS
       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: r_old
       INTEGER, INTENT(IN)                                :: iw
       INTEGER, DIMENSION(1:3, 1:2), INTENT(INOUT)        :: discrete_array
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'mc_volume_move', routineP = moduleN//':'//routineN
 
@@ -1047,7 +1046,7 @@ CONTAINS
 ! now do the move
 
 ! call a random number to figure out how far we're moving
-      IF (ionode) rand = next_random_number(rng_stream)
+      IF (ionode) rand = rng_stream%next()
       CALL mp_bcast(rand, source, group)
 
 ! find the test cell lenghts for the discrete volume move
@@ -1063,7 +1062,7 @@ CONTAINS
 ! if we're increasing the volume, we need to find a side we can increase
          IF (lincrease) THEN
             DO
-               IF (ionode) rand = next_random_number(rng_stream)
+               IF (ionode) rand = rng_stream%next()
                CALL mp_bcast(rand, source, group)
                iside_change = CEILING(3.0_dp*rand)
                IF (discrete_array(iside_change, 1) .EQ. 1) THEN
@@ -1074,7 +1073,7 @@ CONTAINS
             ENDDO
          ELSE
             DO
-               IF (ionode) rand = next_random_number(rng_stream)
+               IF (ionode) rand = rng_stream%next()
                CALL mp_bcast(rand, source, group)
                iside_change = CEILING(3.0_dp*rand)
                IF (discrete_array(iside_change, 2) .EQ. 1) THEN
@@ -1240,7 +1239,7 @@ CONTAINS
          w = 1.0E0_dp
          rand = 0.0E0_dp
       ELSE
-         IF (ionode) rand = next_random_number(rng_stream)
+         IF (ionode) rand = rng_stream%next()
          CALL mp_bcast(rand, source, group)
       ENDIF
 
@@ -1311,7 +1310,7 @@ CONTAINS
       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
       REAL(KIND=dp), INTENT(OUT)                         :: dis_length
       TYPE(particle_list_type), POINTER                  :: particles
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(len=*), PARAMETER :: routineN = 'change_bond_length', &
          routineP = moduleN//':'//routineN
@@ -1347,7 +1346,7 @@ CONTAINS
 
 ! pick which bond in the molecule at random
       IF (ionode) THEN
-         rand = next_random_number(rng_stream)
+         rand = rng_stream%next()
       ENDIF
       CALL mp_bcast(rand, source, group)
       CALL get_molecule_kind(molecule_kind, natom=natom, nbond=nbond, &
@@ -1406,7 +1405,7 @@ CONTAINS
       ENDDO
 
 ! choose a displacement
-      IF (ionode) rand = next_random_number(rng_stream)
+      IF (ionode) rand = rng_stream%next()
       CALL mp_bcast(rand, source, group)
 
       dis_length = rmbond(molecule_type)*2.0E0_dp*(rand - 0.5E0_dp)
@@ -1478,7 +1477,7 @@ CONTAINS
       INTEGER, INTENT(IN)                                :: molecule_type
       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
       TYPE(particle_list_type), POINTER                  :: particles
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(len=*), PARAMETER :: routineN = 'change_bond_angle', &
          routineP = moduleN//':'//routineN
@@ -1517,7 +1516,7 @@ CONTAINS
 
 ! pick which bond in the molecule at random
       IF (ionode) THEN
-         rand = next_random_number(rng_stream)
+         rand = rng_stream%next()
       ENDIF
       CALL mp_bcast(rand, source, group)
       CALL get_molecule_kind(molecule_kind, natom=natom, nbend=nbend, &
@@ -1573,7 +1572,7 @@ CONTAINS
       ENDDO
 
 ! choose a displacement
-      IF (ionode) rand = next_random_number(rng_stream)
+      IF (ionode) rand = rng_stream%next()
       CALL mp_bcast(rand, source, group)
 
       dis_angle = rmangle(molecule_type)*2.0E0_dp*(rand - 0.5E0_dp)
@@ -1738,7 +1737,7 @@ CONTAINS
       INTEGER, INTENT(IN)                                :: molecule_type
       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
       TYPE(particle_list_type), POINTER                  :: particles
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(len=*), PARAMETER :: routineN = 'change_dihedral', &
          routineP = moduleN//':'//routineN
@@ -1778,7 +1777,7 @@ CONTAINS
 
 ! pick which bond in the molecule at random
       IF (ionode) THEN
-         rand = next_random_number(rng_stream)
+         rand = rng_stream%next()
 !      CALL RANDOM_NUMBER(rand)
       ENDIF
       CALL mp_bcast(rand, source, group)
@@ -1837,7 +1836,7 @@ CONTAINS
       ENDDO
 
 ! choose a displacement
-      IF (ionode) rand = next_random_number(rng_stream)
+      IF (ionode) rand = rng_stream%next()
       CALL mp_bcast(rand, source, group)
 
       dis_angle = rmdihedral(molecule_type)*2.0E0_dp*(rand - 0.5E0_dp)
@@ -1959,7 +1958,7 @@ CONTAINS
                                                             molecule_type, box_number
       REAL(KIND=dp), INTENT(INOUT)                       :: bias_energy_old, last_bias_energy
       CHARACTER(LEN=*), INTENT(IN)                       :: move_type
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(len=*), PARAMETER :: routineN = 'mc_avbmc_move', routineP = moduleN//':'//routineN
 
@@ -2247,7 +2246,7 @@ CONTAINS
       IF (w .GE. 1.0E0_dp) THEN
          rand = 0.0E0_dp
       ELSE
-         IF (ionode) rand = next_random_number(rng_stream)
+         IF (ionode) rand = rng_stream%next()
          CALL mp_bcast(rand, source, group)
       ENDIF
 
@@ -2355,7 +2354,7 @@ CONTAINS
       INTEGER, INTENT(IN)                                :: box_number
       REAL(KIND=dp), INTENT(INOUT)                       :: energy_check
       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: r_old
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'mc_hmc_move', routineP = moduleN//':'//routineN
 
@@ -2425,7 +2424,7 @@ CONTAINS
          w = 1.0E0_dp
          rand = 0.0E0_dp
       ELSE
-         IF (ionode) rand = next_random_number(rng_stream)
+         IF (ionode) rand = rng_stream%next()
          CALL mp_bcast(rand, source, group)
       ENDIF
 
@@ -2490,7 +2489,7 @@ CONTAINS
       INTEGER, INTENT(IN)                                :: box_number
       REAL(KIND=dp), INTENT(INOUT)                       :: bias_energy
       LOGICAL, INTENT(OUT)                               :: lreject
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(len=*), PARAMETER :: routineN = 'mc_cluster_translation', &
          routineP = moduleN//':'//routineN
@@ -2586,17 +2585,17 @@ CONTAINS
       ENDIF
 
 ! call a random number to figure out which direction we're moving
-      IF (ionode) rand = next_random_number(rng_stream)
+      IF (ionode) rand = rng_stream%next()
       CALL mp_bcast(rand, source, group)
       move_direction = INT(3*rand) + 1
 
 ! call a random number to figure out how far we're moving
-      IF (ionode) rand = next_random_number(rng_stream)
+      IF (ionode) rand = rng_stream%next()
       CALL mp_bcast(rand, source, group)
       dis_mol = rmcltrans*(rand - 0.5E0_dp)*2.0E0_dp
 
 ! choosing cluster
-      IF (ionode) rand = next_random_number(rng_stream)
+      IF (ionode) rand = rng_stream%next()
       CALL mp_bcast(rand, source, group)
       jpart = INT(1 + rand*total_clus)
 
@@ -2685,7 +2684,7 @@ CONTAINS
             w = 1.0E0_dp
             rand = 0.0E0_dp
          ELSE
-            IF (ionode) rand = next_random_number(rng_stream)
+            IF (ionode) rand = rng_stream%next()
             CALL mp_bcast(rand, source, group)
          ENDIF
          IF (rand .LT. w) THEN
@@ -2725,4 +2724,3 @@ CONTAINS
    END SUBROUTINE mc_cluster_translation
 
 END MODULE mc_moves
-
diff --git a/src/motion/mc/mc_run.F b/src/motion/mc/mc_run.F
index c7ad9d6f3b..d936bdeba7 100644
--- a/src/motion/mc/mc_run.F
+++ b/src/motion/mc/mc_run.F
@@ -66,8 +66,6 @@ MODULE mc_run
         mc_sim_par_create, mc_sim_par_destroy, mc_simulation_parameters_p_type, read_mc_section, &
         set_mc_par
    USE parallel_rng_types,              ONLY: UNIFORM,&
-                                              create_rng_stream,&
-                                              delete_rng_stream,&
                                               rng_stream_type
    USE physcon,                         ONLY: angstrom
 #include "../../base/base_uses.f90"
@@ -126,11 +124,11 @@ CONTAINS
       TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info
       TYPE(mc_simulation_parameters_p_type), &
          DIMENSION(:), POINTER                           :: mc_par
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type)                              :: rng_stream
       TYPE(section_vals_type), POINTER                   :: force_env_section, mc_section, &
                                                             root_section
 
-      NULLIFY (mc_env, mc_par, force_env_2, rng_stream, force_env_section, &
+      NULLIFY (mc_env, mc_par, force_env_2, force_env_section, &
                root_section, mc_molecule_info)
 
       CALL force_env_retain(force_env_1)
@@ -144,10 +142,10 @@ CONTAINS
 ! set some values...will use get_globenv if that ever comes around
 
 ! initialize the random numbers
-      CALL create_rng_stream(rng_stream=rng_stream, &
-                             last_rng_stream=force_env_1%globenv%gaussian_rng_stream, &
-                             name="first", &
-                             distribution_type=UNIFORM)
+      rng_stream = rng_stream_type( &
+                   last_rng_stream=force_env_1%globenv%gaussian_rng_stream, &
+                   name="first", &
+                   distribution_type=UNIFORM)
 
 ! need to figure out how many boxes we have, based on the value
 ! of mc_par % ensemble
@@ -338,9 +336,6 @@ CONTAINS
       DEALLOCATE (mc_env)
       DEALLOCATE (force_env)
 
-! delete the random numbers
-      CALL delete_rng_stream(rng_stream)
-
    END SUBROUTINE do_mon_car
 
 ! **************************************************************************************************
diff --git a/src/motion/mc/tamc_run.F b/src/motion/mc/tamc_run.F
index 9c0250f30b..4b956419f3 100644
--- a/src/motion/mc/tamc_run.F
+++ b/src/motion/mc/tamc_run.F
@@ -99,9 +99,6 @@ MODULE tamc_run
    USE molecule_types,                  ONLY: global_constraint_type,&
                                               molecule_type
    USE parallel_rng_types,              ONLY: UNIFORM,&
-                                              create_rng_stream,&
-                                              delete_rng_stream,&
-                                              next_random_number,&
                                               rng_stream_type
    USE particle_list_types,             ONLY: particle_list_type
    USE particle_types,                  ONLY: particle_type
@@ -190,7 +187,7 @@ CONTAINS
       TYPE(meta_env_type), POINTER                       :: meta_env_saved
       TYPE(particle_list_type), POINTER                  :: particles
       TYPE(reftraj_type), POINTER                        :: reftraj
-      TYPE(rng_stream_type), POINTER                     :: rng_stream, rng_stream_mc
+      TYPE(rng_stream_type)                              :: rng_stream_mc
       TYPE(section_vals_type), POINTER :: constraint_section, force_env_section, &
          free_energy_section, fs_section, global_section, mc_section, md_section, motion_section, &
          reftraj_section, subsys_section, work_section
@@ -308,7 +305,7 @@ CONTAINS
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! MC setup up
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-      NULLIFY (mc_env, mc_par, rng_stream_mc, MCaverages)
+      NULLIFY (mc_env, mc_par, MCaverages)
 
       CALL section_vals_get(force_env_section, n_repetition=isos)
       CPASSERT(isos == 1)
@@ -316,9 +313,8 @@ CONTAINS
 
 ! initialize the random numbers
 !       IF (para_env%ionode) THEN
-      CALL create_rng_stream(rng_stream=rng_stream_mc, &
-                             name="Random numbers for monte carlo acc/rej", &
-                             distribution_type=UNIFORM)
+      rng_stream_mc = rng_stream_type(name="Random numbers for monte carlo acc/rej", &
+                                      distribution_type=UNIFORM)
 !       ENDIF
 !!!!! this shoudl go in a routine hmc_read
 
@@ -361,9 +357,9 @@ CONTAINS
                          particles=particles, virial=virial)
 
       DO i = 1, rand2skip
-         auxRandom = next_random_number(rng_stream_mc)
+         auxRandom = rng_stream_mc%next()
          DO j = 1, 3*SIZE(particles%els)
-            auxRandom = next_random_number(globenv%gaussian_rng_stream)
+            auxRandom = globenv%gaussian_rng_stream%next()
          ENDDO
       ENDDO
 
@@ -435,12 +431,10 @@ CONTAINS
          IF (ASSOCIATED(force_env%meta_env)) THEN
             IF (force_env%meta_env%langevin) THEN
                CALL create_wiener_process_cv(force_env%meta_env)
-               NULLIFY (rng_stream)
                DO j = 1, (rand2skip - 1)/nmccycles
                   DO i = 1, force_env%meta_env%n_colvar
-                     rng_stream => force_env%meta_env%rng(i)%stream
-                     auxRandom = next_random_number(rng_stream)
-                     auxRandom = next_random_number(rng_stream)
+                     auxRandom = force_env%meta_env%rng(i)%next()
+                     auxRandom = force_env%meta_env%rng(i)%next()
                   ENDDO
                ENDDO
             ENDIF
@@ -652,7 +646,6 @@ CONTAINS
       ! Clean restartable sections..
       IF (my_rm_restart_info) CALL remove_restart_info(force_env%root_section)
 !     IF (para_env%ionode) THEN
-      CALL delete_rng_stream(rng_stream_mc)
 !     ENDIF
       CALL MC_ENV_RELEASE(mc_env)
       DEALLOCATE (mc_par)
@@ -884,7 +877,7 @@ CONTAINS
       TYPE(mc_environment_type), POINTER                 :: mc_env
       TYPE(mc_moves_type), POINTER                       :: moves, gmoves
       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: r
-      TYPE(rng_stream_type), POINTER                     :: rng_stream_mc
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream_mc
       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: xieta, An, fz
       TYPE(mc_averages_type), INTENT(INOUT), POINTER     :: averages
       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: zbuff
@@ -910,7 +903,6 @@ CONTAINS
       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
       TYPE(particle_list_type), POINTER                  :: particles
       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
       TYPE(simpar_type), POINTER                         :: simpar
       TYPE(virial_type), POINTER                         :: virial
 
@@ -918,7 +910,6 @@ CONTAINS
       logger => cp_get_default_logger()
       output_unit = cp_logger_get_default_io_unit(logger)
 
-      NULLIFY (rng_stream)
 ! quantitites to be nullified for the get_md_env
       NULLIFY (simpar, force_env, para_env)
 ! quantities to be nullified for the force_env_get environment
@@ -955,9 +946,8 @@ CONTAINS
       !    *** Velocity Verlet for Langevin *** v(t)--> v(t+1/2)
       !!!!!! noise xi is in the first half, eta in the second half
       DO ivar = 1, force_env%meta_env%n_colvar
-         rng_stream => force_env%meta_env%rng(ivar)%stream
-         xieta(ivar) = next_random_number(rng_stream)
-         xieta(ivar + force_env%meta_env%n_colvar) = next_random_number(rng_stream)
+         xieta(ivar) = force_env%meta_env%rng(ivar)%next()
+         xieta(ivar + force_env%meta_env%n_colvar) = force_env%meta_env%rng(ivar)%next()
          gamma = force_env%meta_env%metavar(ivar)%gamma
          mass = force_env%meta_env%metavar(ivar)%mass
          sigma = SQRT((force_env%meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*gamma/mass)
@@ -1010,7 +1000,7 @@ CONTAINS
       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: r
       TYPE(mc_simpar_type), POINTER                      :: mc_par
       TYPE(mc_moves_type), POINTER                       :: moves, gmoves
-      TYPE(rng_stream_type), POINTER                     :: rng_stream_mc
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream_mc
       INTEGER, INTENT(IN)                                :: output_unit
       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: fz, zbuff
       INTEGER, INTENT(IN), OPTIONAL                      :: nskip
@@ -1146,7 +1136,7 @@ CONTAINS
       REAL(KIND=dp), INTENT(INOUT)                       :: old_epx, old_epz, energy_check
       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: r
       INTEGER, INTENT(IN)                                :: output_unit
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: zbuff
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'mc_hmc_move', routineP = moduleN//':'//routineN
@@ -1240,7 +1230,7 @@ CONTAINS
          w = EXP(value)
       ENDIF
 
-      rand = next_random_number(rng_stream)
+      rand = rng_stream%next()
       IF (rand < w) THEN
 ! accept the move
          moves%hmc%successes = moves%hmc%successes + 1
@@ -1444,7 +1434,7 @@ CONTAINS
       meta_env%ekin_s = 0.0_dp
       DO i_c = 1, meta_env%n_colvar
          cv => meta_env%metavar(i_c)
-         cv%vvp = next_random_number(force_env%globenv%gaussian_rng_stream)
+         cv%vvp = force_env%globenv%gaussian_rng_stream%next()
          meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
       END DO
       ekin_w = 0.5_dp*meta_env%temp_wanted*REAL(meta_env%n_colvar, KIND=dp)
diff --git a/src/motion/md_vel_utils.F b/src/motion/md_vel_utils.F
index e253c8672c..9b9d27f214 100644
--- a/src/motion/md_vel_utils.F
+++ b/src/motion/md_vel_utils.F
@@ -66,10 +66,6 @@ MODULE md_vel_utils
                                               get_molecule_kind_set,&
                                               molecule_kind_type
    USE parallel_rng_types,              ONLY: UNIFORM,&
-                                              create_rng_stream,&
-                                              delete_rng_stream,&
-                                              next_random_number,&
-                                              random_numbers,&
                                               rng_stream_type
    USE particle_list_types,             ONLY: particle_list_type
    USE particle_types,                  ONLY: particle_type
@@ -920,9 +916,8 @@ CONTAINS
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: dr, eigenvectors
       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
       TYPE(cp_para_env_type), POINTER                    :: para_env
-      TYPE(rng_stream_type), POINTER                     :: random_stream
+      TYPE(rng_stream_type), ALLOCATABLE                 :: random_stream
 
-      NULLIFY (random_stream)
       CALL cite_reference(West2006)
       natoms = SIZE(particles)
       temperature = simpar%temp_ext
@@ -949,14 +944,15 @@ CONTAINS
       CALL section_vals_val_get(md_section, "INITIAL_VIBRATION%PHASE", r_val=my_phase)
       my_phase = MIN(1.0_dp, my_phase)
       ! generate random numbers
-      CALL create_rng_stream(random_stream, name="MD_INIT_VIB", distribution_type=UNIFORM)
-      CALL random_numbers(random, random_stream)
+      random_stream = rng_stream_type(name="MD_INIT_VIB", distribution_type=UNIFORM)
+      CALL random_stream%fill(random)
       IF (my_phase .LT. 0.0_dp) THEN
-         CALL random_numbers(phase, random_stream)
+         CALL random_stream%fill(phase)
       ELSE
          phase = my_phase
       END IF
-      CALL delete_rng_stream(random_stream)
+      DEALLOCATE (random_stream)
+
       ! the first three modes are acoustic with zero frequencies,
       ! exclude these from considerations
       my_dof = dof - 3
@@ -1234,24 +1230,24 @@ CONTAINS
          IF (mass .NE. 0.0) THEN
             SELECT CASE (is_fixed(i))
             CASE (use_perd_x)
-               part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream)/SQRT(mass)
-               part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream)/SQRT(mass)
+               part(i)%v(2) = globenv%gaussian_rng_stream%next()/SQRT(mass)
+               part(i)%v(3) = globenv%gaussian_rng_stream%next()/SQRT(mass)
             CASE (use_perd_y)
-               part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream)/SQRT(mass)
-               part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream)/SQRT(mass)
+               part(i)%v(1) = globenv%gaussian_rng_stream%next()/SQRT(mass)
+               part(i)%v(3) = globenv%gaussian_rng_stream%next()/SQRT(mass)
             CASE (use_perd_z)
-               part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream)/SQRT(mass)
-               part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream)/SQRT(mass)
+               part(i)%v(1) = globenv%gaussian_rng_stream%next()/SQRT(mass)
+               part(i)%v(2) = globenv%gaussian_rng_stream%next()/SQRT(mass)
             CASE (use_perd_xy)
-               part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream)/SQRT(mass)
+               part(i)%v(3) = globenv%gaussian_rng_stream%next()/SQRT(mass)
             CASE (use_perd_xz)
-               part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream)/SQRT(mass)
+               part(i)%v(2) = globenv%gaussian_rng_stream%next()/SQRT(mass)
             CASE (use_perd_yz)
-               part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream)/SQRT(mass)
+               part(i)%v(1) = globenv%gaussian_rng_stream%next()/SQRT(mass)
             CASE (use_perd_none)
-               part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream)/SQRT(mass)
-               part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream)/SQRT(mass)
-               part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream)/SQRT(mass)
+               part(i)%v(1) = globenv%gaussian_rng_stream%next()/SQRT(mass)
+               part(i)%v(2) = globenv%gaussian_rng_stream%next()/SQRT(mass)
+               part(i)%v(3) = globenv%gaussian_rng_stream%next()/SQRT(mass)
             END SELECT
          END IF
       END DO
diff --git a/src/motion/neb_md_utils.F b/src/motion/neb_md_utils.F
index 3f53c1e723..0da62ab41d 100644
--- a/src/motion/neb_md_utils.F
+++ b/src/motion/neb_md_utils.F
@@ -27,7 +27,6 @@ MODULE neb_md_utils
    USE kinds,                           ONLY: dp
    USE neb_types,                       ONLY: neb_type,&
                                               neb_var_type
-   USE parallel_rng_types,              ONLY: next_random_number
    USE particle_types,                  ONLY: get_particle_pos_or_vel,&
                                               particle_type,&
                                               update_particle_pos_or_vel
@@ -73,7 +72,7 @@ CONTAINS
       CHARACTER(LEN=*), PARAMETER :: routineN = 'neb_initialize_velocity', &
          routineP = moduleN//':'//routineN
 
-      INTEGER                                            :: iatom, ivar, k, natom, nparticle, nvar
+      INTEGER                                            :: iatom, ivar, natom, nparticle, nvar
       REAL(KIND=dp)                                      :: akin, mass, mass_tot, sc, temp, &
                                                             temp_ext, tmp_r1
       REAL(KIND=dp), DIMENSION(3)                        :: v, vcom
@@ -87,10 +86,7 @@ CONTAINS
          natom = SIZE(particle_set)
          mass_tot = 0.0_dp
          vcom(1:3) = 0.0_dp
-         vels(:, i_rep) = 0.0_dp
-         DO k = 1, nparticle
-            vels(k, i_rep) = next_random_number(globenv%gaussian_rng_stream)
-         END DO
+         CALL globenv%gaussian_rng_stream%fill(vels(:, i_rep))
          ! Check always if BAND is working in Cartesian or in internal coordinates
          ! If working in cartesian coordinates let's get rid of the COM
          ! Compute also the total mass (both in Cartesian and internal)
diff --git a/src/motion/pint_gle.F b/src/motion/pint_gle.F
index f6ae04d292..c252231233 100644
--- a/src/motion/pint_gle.F
+++ b/src/motion/pint_gle.F
@@ -15,8 +15,6 @@ MODULE pint_gle
    USE gle_system_dynamics,             ONLY: gle_cholesky_stab
    USE gle_system_types,                ONLY: gle_type
    USE kinds,                           ONLY: dp
-   USE parallel_rng_types,              ONLY: next_random_number,&
-                                              rng_stream_type
    USE pint_types,                      ONLY: pint_env_type
 #include "../base/base_uses.f90"
 
@@ -55,7 +53,6 @@ CONTAINS
 
       INTEGER                                            :: i, ib, idim, imap, j
       REAL(dp) :: mf, rr(pint_env%gle%ndim), cc(pint_env%gle%ndim, pint_env%gle%ndim)
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
 
       CALL gle_cholesky_stab(pint_env%gle%c_mat, cc, pint_env%gle%ndim)
       DO i = 1, pint_env%gle%loc_num_gle
@@ -63,9 +60,8 @@ CONTAINS
          ib = 1 + (imap - 1)/pint_env%ndim
          idim = 1 + MOD(imap - 1, pint_env%ndim)
          mf = 1.0_dp/SQRT(pint_env%mass_fict(ib, idim))
-         rng_stream => pint_env%gle%nvt(i)%gaussian_rng_stream
          DO j = 1, pint_env%gle%ndim
-            rr(j) = next_random_number(rng_stream)*mf
+            rr(j) = pint_env%gle%nvt(i)%gaussian_rng_stream%next()*mf
          END DO
          pint_env%gle%nvt(i)%s = MATMUL(cc, rr)
       END DO
@@ -86,12 +82,10 @@ CONTAINS
       REAL(dp)                                           :: alpha, beta, mf, rr
       REAL(dp), DIMENSION(:, :), POINTER                 :: a_mat, e_tmp, h_tmp, s_tmp
       TYPE(gle_type), POINTER                            :: gle
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
 
       CALL timeset(routineN, handle)
 
       gle => pint_env%gle
-      NULLIFY (rng_stream)
       ndim = gle%ndim
 
       ALLOCATE (s_tmp(ndim, gle%loc_num_gle))
@@ -108,13 +102,12 @@ CONTAINS
          gle%nvt(ideg)%thermostat_energy = gle%nvt(ideg)%thermostat_energy &
                                            + 0.5_dp*pint_env%mass_fict(ib, idim)*gle%nvt(ideg)%s(1)**2
          s_tmp(1, imap) = gle%nvt(ideg)%s(1)
-         rng_stream => gle%nvt(ideg)%gaussian_rng_stream
-         rr = next_random_number(rng_stream)
+         rr = gle%nvt(ideg)%gaussian_rng_stream%next()
          mf = 1.0_dp/SQRT(pint_env%mass_fict(ib, idim))
          e_tmp(1, imap) = rr*mf
          DO iadd = 2, ndim
             s_tmp(iadd, imap) = gle%nvt(ideg)%s(iadd)
-            rr = next_random_number(rng_stream)
+            rr = gle%nvt(ideg)%gaussian_rng_stream%next()
             e_tmp(iadd, imap) = rr*mf
          END DO
       END DO
diff --git a/src/motion/pint_methods.F b/src/motion/pint_methods.F
index a4b4aed7e4..ad6d0d24e5 100644
--- a/src/motion/pint_methods.F
+++ b/src/motion/pint_methods.F
@@ -89,10 +89,6 @@ MODULE pint_methods
    USE molecule_types,                  ONLY: global_constraint_type,&
                                               molecule_type
    USE parallel_rng_types,              ONLY: GAUSSIAN,&
-                                              create_rng_stream,&
-                                              delete_rng_stream,&
-                                              next_random_number,&
-                                              next_rng_seed,&
                                               rng_stream_type
    USE particle_list_types,             ONLY: particle_list_type
    USE particle_types,                  ONLY: particle_type
@@ -185,7 +181,6 @@ CONTAINS
                                                             ierr, ig, itmp, nrep, prep, stat
       LOGICAL                                            :: explicit, ltmp
       REAL(kind=dp)                                      :: dt, mass, omega
-      REAL(kind=dp), DIMENSION(3, 2)                     :: seed
       TYPE(cp_subsys_type), POINTER                      :: subsys
       TYPE(f_env_type), POINTER                          :: f_env
       TYPE(global_constraint_type), POINTER              :: gci
@@ -338,14 +333,10 @@ CONTAINS
 !MK ... but we have to initialise v_tol
       pint_env%v_tol = 0.0_dp ! to be fixed
 
-      NULLIFY (pint_env%randomG)
-
-      seed(:, :) = next_rng_seed()
-      CALL create_rng_stream(pint_env%randomG, &
-                             name="pint_randomG", &
-                             distribution_type=GAUSSIAN, &
-                             extended_precision=.TRUE., &
-                             seed=seed)
+      pint_env%randomG = rng_stream_type( &
+                         name="pint_randomG", &
+                         distribution_type=GAUSSIAN, &
+                         extended_precision=.TRUE.)
 
       ALLOCATE (pint_env%e_pot_bead(pint_env%p))
       pint_env%e_pot_bead = 0._dp
@@ -731,7 +722,6 @@ CONTAINS
             IF (ASSOCIATED(pint_env%normalmode_env)) THEN
                CALL normalmode_release(pint_env%normalmode_env)
             END IF
-            CALL delete_rng_stream(pint_env%randomG)
 
             DEALLOCATE (pint_env%mass)
             DEALLOCATE (pint_env%e_pot_bead)
@@ -1052,7 +1042,7 @@ CONTAINS
       REAL(kind=dp), DIMENSION(3)                        :: x0
       REAL(kind=dp), DIMENSION(3, 2)                     :: seed
       REAL(kind=dp), DIMENSION(:), POINTER               :: bx, r_vals
-      TYPE(rng_stream_type), POINTER                     :: rng_gaussian
+      TYPE(rng_stream_type), ALLOCATABLE                 :: rng_gaussian
       TYPE(section_vals_type), POINTER                   :: input_section
 
       CPASSERT(ASSOCIATED(pint_env))
@@ -1075,16 +1065,15 @@ CONTAINS
 
          NULLIFY (bx)
          ALLOCATE (bx(3*pint_env%p))
-         NULLIFY (rng_gaussian)
          CALL section_vals_val_get(pint_env%input, &
                                    "MOTION%PINT%INIT%LEVY_SEED", i_val=input_seed)
          seed(:, :) = REAL(input_seed, KIND=dp)
 !      seed(:,:) = next_rng_seed()
-         CALL create_rng_stream(rng_gaussian, &
-                                name="tmp_rng_gaussian", &
-                                distribution_type=GAUSSIAN, &
-                                extended_precision=.TRUE., &
-                                seed=seed)
+         rng_gaussian = rng_stream_type( &
+                        name="tmp_rng_gaussian", &
+                        distribution_type=GAUSSIAN, &
+                        extended_precision=.TRUE., &
+                        seed=seed)
 
          CALL section_vals_val_get(pint_env%input, &
                                    "MOTION%PINT%INIT%LEVY_CORRELATED", &
@@ -1126,7 +1115,6 @@ CONTAINS
 
          END IF
 
-         CALL delete_rng_stream(rng_gaussian)
          DEALLOCATE (bx)
          done_levy = .TRUE.
       END IF
@@ -1164,9 +1152,8 @@ CONTAINS
          DO idim = 1, pint_env%ndim
             DO ib = 1, pint_env%p
                pint_env%x(ib, idim) = pint_env%x(ib, idim) + &
-                                      next_random_number(rng_stream=pint_env%randomG, &
-                                                         variance=pint_env%beta/ &
-                                                         SQRT(12.0_dp*pint_env%mass(idim)))
+                                      pint_env%randomG%next(variance=pint_env%beta/ &
+                                                            SQRT(12.0_dp*pint_env%mass(idim)))
             END DO
          END DO
          done_rand = .TRUE.
@@ -1386,8 +1373,7 @@ CONTAINS
       DO idim = 1, SIZE(pint_env%uv, 2)
          DO ib = first_mode, SIZE(pint_env%uv, 1)
             pint_env%uv(ib, idim) = &
-               next_random_number(rng_stream=pint_env%randomG, &
-                                  variance=target_t/pint_env%mass_fict(ib, idim))
+               pint_env%randomG%next(variance=target_t/pint_env%mass_fict(ib, idim))
          END DO
       END DO
 
@@ -1399,8 +1385,7 @@ CONTAINS
       IF (ltmp) THEN
          CALL pint_u2x(pint_env, ux=pint_env%uv, x=pint_env%v)
          DO idim = 1, pint_env%ndim
-            rtmp = next_random_number(rng_stream=pint_env%randomG, &
-                                      variance=pint_env%mass(idim)*pint_env%kT) &
+            rtmp = pint_env%randomG%next(variance=pint_env%mass(idim)*pint_env%kT) &
                    /pint_env%mass(idim)
             DO ib = 1, pint_env%p
                pint_env%v(ib, idim) = pint_env%v(ib, idim) + rtmp
@@ -1569,8 +1554,7 @@ CONTAINS
             DO ib = 1, SIZE(pint_env%tv, 2)
                DO inos = 1, SIZE(pint_env%tv, 1)
                   pint_env%tv(inos, ib, idim) = &
-                     next_random_number(rng_stream=pint_env%randomG, &
-                                        variance=mykt/pint_env%Q(ib))
+                     pint_env%randomG%next(variance=mykt/pint_env%Q(ib))
                END DO
             END DO
          END DO
diff --git a/src/motion/pint_piglet.F b/src/motion/pint_piglet.F
index c9a390f8dd..037707fe95 100644
--- a/src/motion/pint_piglet.F
+++ b/src/motion/pint_piglet.F
@@ -29,11 +29,9 @@ MODULE pint_piglet
                                               dp
    USE message_passing,                 ONLY: mp_bcast
    USE parallel_rng_types,              ONLY: GAUSSIAN,&
-                                              create_rng_stream,&
-                                              delete_rng_stream,&
-                                              next_random_number,&
-                                              read_rng_stream,&
-                                              rng_record_length
+                                              rng_record_length,&
+                                              rng_stream_type,&
+                                              rng_stream_type_from_record
    USE pint_io,                         ONLY: pint_write_line
    USE pint_types,                      ONLY: piglet_therm_type,&
                                               pint_env_type
@@ -277,21 +275,19 @@ CONTAINS
 
       !prepare Random number generator
       NULLIFY (rng_section)
-      NULLIFY (piglet_therm%gaussian_rng_stream)
       rng_section => section_vals_get_subs_vals(section, &
                                                 subsection_name="RNG_INIT")
       CALL section_vals_get(rng_section, explicit=explicit)
       IF (explicit) THEN
          CALL section_vals_val_get(rng_section, "_DEFAULT_KEYWORD_", &
                                    i_rep_val=1, c_val=rng_record)
-         CALL read_rng_stream(rng_stream=piglet_therm%gaussian_rng_stream, &
-                              rng_record=rng_record)
+         piglet_therm%gaussian_rng_stream = rng_stream_type_from_record(rng_record)
       ELSE
          initial_seed(:, :) = REAL(pint_env%thermostat_rng_seed, dp)
-         CALL create_rng_stream(rng_stream=piglet_therm%gaussian_rng_stream, &
-                                name="piglet_rng_gaussian", distribution_type=GAUSSIAN, &
-                                extended_precision=.TRUE., &
-                                seed=initial_seed)
+         piglet_therm%gaussian_rng_stream = rng_stream_type( &
+                                            name="piglet_rng_gaussian", distribution_type=GAUSSIAN, &
+                                            extended_precision=.TRUE., &
+                                            seed=initial_seed)
       END IF
 
       !Compute the T and S matrices on every mpi process
@@ -385,7 +381,7 @@ CONTAINS
             ! Fill a vector with random numbers
             DO idim = 1, piglet_therm%ndim
                DO j = 1, piglet_therm%nsp1
-                  piglet_therm%temp2(j, idim) = next_random_number(piglet_therm%gaussian_rng_stream)
+                  piglet_therm%temp2(j, idim) = piglet_therm%gaussian_rng_stream%next()
                   !piglet_therm%temp2(j,idim) = 1.0_dp
                END DO
             END DO
@@ -469,7 +465,7 @@ CONTAINS
          !fill temp2 with gaussian random noise
          DO j = 1, nsp1
             DO idim = 1, ndim
-               piglet_therm%temp2(j, idim) = next_random_number(piglet_therm%gaussian_rng_stream)
+               piglet_therm%temp2(j, idim) = piglet_therm%gaussian_rng_stream%next()
                !piglet_therm%temp2(j,idim) = 1.0_dp
             END DO
          END DO
@@ -550,7 +546,6 @@ CONTAINS
             DEALLOCATE (piglet_therm%c_mat)
             DEALLOCATE (piglet_therm%gle_t)
             DEALLOCATE (piglet_therm%gle_s)
-            CALL delete_rng_stream(piglet_therm%gaussian_rng_stream)
             DEALLOCATE (piglet_therm%smalls)
             DEALLOCATE (piglet_therm%temp1)
             DEALLOCATE (piglet_therm%temp2)
diff --git a/src/motion/pint_pile.F b/src/motion/pint_pile.F
index 4ba4b48809..5edf169151 100644
--- a/src/motion/pint_pile.F
+++ b/src/motion/pint_pile.F
@@ -17,11 +17,9 @@ MODULE pint_pile
                                               section_vals_val_get
    USE kinds,                           ONLY: dp
    USE parallel_rng_types,              ONLY: GAUSSIAN,&
-                                              create_rng_stream,&
-                                              delete_rng_stream,&
-                                              next_random_number,&
-                                              read_rng_stream,&
-                                              rng_record_length
+                                              rng_record_length,&
+                                              rng_stream_type,&
+                                              rng_stream_type_from_record
    USE pint_types,                      ONLY: normalmode_env_type,&
                                               pile_therm_type,&
                                               pint_env_type
@@ -105,21 +103,20 @@ CONTAINS
 
       !prepare Random number generator
       NULLIFY (rng_section)
-      NULLIFY (pile_therm%gaussian_rng_stream)
       rng_section => section_vals_get_subs_vals(section, &
                                                 subsection_name="RNG_INIT")
       CALL section_vals_get(rng_section, explicit=explicit)
       IF (explicit) THEN
          CALL section_vals_val_get(rng_section, "_DEFAULT_KEYWORD_", &
                                    i_rep_val=1, c_val=rng_record)
-         CALL read_rng_stream(rng_stream=pile_therm%gaussian_rng_stream, &
-                              rng_record=rng_record)
+
+         pile_therm%gaussian_rng_stream = rng_stream_type_from_record(rng_record)
       ELSE
          initial_seed(:, :) = REAL(pint_env%thermostat_rng_seed, dp)
-         CALL create_rng_stream(rng_stream=pile_therm%gaussian_rng_stream, &
-                                name="pile_rng_gaussian", distribution_type=GAUSSIAN, &
-                                extended_precision=.TRUE., &
-                                seed=initial_seed)
+         pile_therm%gaussian_rng_stream = rng_stream_type( &
+                                          name="pile_rng_gaussian", distribution_type=GAUSSIAN, &
+                                          extended_precision=.TRUE., &
+                                          seed=initial_seed)
       END IF
 
    END SUBROUTINE pint_pile_init
@@ -151,7 +148,7 @@ CONTAINS
          DO ibead = first_mode, p
             vnew(ibead, idim) = pile_therm%c1(ibead)*vold(ibead, idim) + &
                                 pile_therm%massfact(ibead, idim)*pile_therm%c2(ibead)* &
-                                next_random_number(pile_therm%gaussian_rng_stream)
+                                pile_therm%gaussian_rng_stream%next()
             delta_ekin = delta_ekin + masses(ibead, idim)*( &
                          vnew(ibead, idim)*vnew(ibead, idim) - &
                          vold(ibead, idim)*vold(ibead, idim))
@@ -185,7 +182,6 @@ CONTAINS
             DEALLOCATE (pile_therm%c2)
             DEALLOCATE (pile_therm%g_fric)
             DEALLOCATE (pile_therm%massfact)
-            CALL delete_rng_stream(pile_therm%gaussian_rng_stream)
             DEALLOCATE (pile_therm)
          END IF
       END IF
diff --git a/src/motion/pint_public.F b/src/motion/pint_public.F
index d23794e0d3..da224844fa 100644
--- a/src/motion/pint_public.F
+++ b/src/motion/pint_public.F
@@ -13,8 +13,7 @@
 MODULE pint_public
 
    USE kinds,                           ONLY: dp
-   USE parallel_rng_types,              ONLY: next_random_number,&
-                                              rng_stream_type
+   USE parallel_rng_types,              ONLY: rng_stream_type
    USE pint_types,                      ONLY: pint_env_type
 #include "../base/base_uses.f90"
 
@@ -122,7 +121,7 @@ CONTAINS
 !
       INTEGER, INTENT(IN)                                :: n
       REAL(kind=dp), INTENT(IN)                          :: t
-      TYPE(rng_stream_type), POINTER                     :: rng_gaussian
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_gaussian
       REAL(kind=dp), DIMENSION(:), POINTER               :: x
       INTEGER, INTENT(OUT)                               :: nout
 
@@ -177,8 +176,7 @@ CONTAINS
             ! generate new point and save it under j
             DO ic = 1, 3
                xc = (x(3*i1 + ic) + x(3*i2 + ic))/2.0
-               xc = xc + next_random_number(rng_stream=rng_gaussian, &
-                                            variance=vrnc)
+               xc = xc + rng_gaussian%next(variance=vrnc)
                x(3*j + ic) = xc
             END DO
             nout = nout + 1
@@ -220,7 +218,7 @@ CONTAINS
       INTEGER, INTENT(IN)                                :: n
       REAL(kind=dp), INTENT(IN)                          :: v
       REAL(kind=dp), DIMENSION(:), POINTER               :: x
-      TYPE(rng_stream_type), POINTER                     :: rng_gaussian
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_gaussian
 
       CHARACTER(len=*), PARAMETER :: routineN = 'pint_levy_walk', routineP = moduleN//':'//routineN
 
@@ -233,8 +231,7 @@ CONTAINS
       x(3) = x0(3)
       DO ib = 1, n - 1
          DO ic = 1, 3
-            r = next_random_number(rng_stream=rng_gaussian, &
-                                   variance=1.0_dp)
+            r = rng_gaussian%next(variance=1.0_dp)
             tau_i = (REAL(ib, dp) - 1.0_dp)/REAL(n, dp)
             tau_i1 = (REAL(ib + 1, dp) - 1.0_dp)/REAL(n, dp)
             x(ib*3 + ic) = (x((ib - 1)*3 + ic)*(1.0_dp - tau_i1) + &
diff --git a/src/motion/pint_qtb.F b/src/motion/pint_qtb.F
index fdd21c4f39..de0854109f 100644
--- a/src/motion/pint_qtb.F
+++ b/src/motion/pint_qtb.F
@@ -34,12 +34,9 @@ MODULE pint_qtb
                                               twopi
    USE message_passing,                 ONLY: mp_bcast
    USE parallel_rng_types,              ONLY: GAUSSIAN,&
-                                              create_rng_stream,&
-                                              delete_rng_stream,&
-                                              dump_rng_stream,&
-                                              next_random_number,&
-                                              read_rng_stream,&
-                                              rng_record_length
+                                              rng_record_length,&
+                                              rng_stream_type,&
+                                              rng_stream_type_from_record
    USE pint_io,                         ONLY: pint_write_line
    USE pint_types,                      ONLY: normalmode_env_type,&
                                               pint_env_type,&
@@ -128,21 +125,19 @@ CONTAINS
 
       !prepare Random number generator
       NULLIFY (rng_section)
-      NULLIFY (qtb_therm%gaussian_rng_stream)
       rng_section => section_vals_get_subs_vals(section, &
                                                 subsection_name="RNG_INIT")
       CALL section_vals_get(rng_section, explicit=restart)
       IF (restart) THEN
          CALL section_vals_val_get(rng_section, "_DEFAULT_KEYWORD_", &
                                    i_rep_val=1, c_val=rng_record)
-         CALL read_rng_stream(rng_stream=qtb_therm%gaussian_rng_stream, &
-                              rng_record=rng_record)
+         qtb_therm%gaussian_rng_stream = rng_stream_type_from_record(rng_record)
       ELSE
          initial_seed(:, :) = REAL(pint_env%thermostat_rng_seed, dp)
-         CALL create_rng_stream(rng_stream=qtb_therm%gaussian_rng_stream, &
-                                name="qtb_rng_gaussian", distribution_type=GAUSSIAN, &
-                                extended_precision=.TRUE., &
-                                seed=initial_seed)
+         qtb_therm%gaussian_rng_stream = rng_stream_type( &
+                                         name="qtb_rng_gaussian", distribution_type=GAUSSIAN, &
+                                         extended_precision=.TRUE., &
+                                         seed=initial_seed)
       END IF
 
       !Initialization of the QTB random forces
@@ -183,15 +178,14 @@ CONTAINS
                DO i = 1, qtb_therm%nf - 1
                   qtb_therm%rng_status(i) = qtb_therm%rng_status(i + 1)
                END DO
-               CALL dump_rng_stream(rng_stream=qtb_therm%gaussian_rng_stream, &
-                                    rng_record=qtb_therm%rng_status(qtb_therm%nf))
+               CALL qtb_therm%gaussian_rng_stream%dump(qtb_therm%rng_status(qtb_therm%nf))
             END IF
             DO idim = 1, ndim
                !update random numbers
                DO i = 1, qtb_therm%nf - 1
                   qtb_therm%r(i, ibead, idim) = qtb_therm%r(i + 1, ibead, idim)
                END DO
-               qtb_therm%r(qtb_therm%nf, ibead, idim) = next_random_number(qtb_therm%gaussian_rng_stream)
+               qtb_therm%r(qtb_therm%nf, ibead, idim) = qtb_therm%gaussian_rng_stream%next()
                !compute new random force through the convolution product
                qtb_therm%rf(ibead, idim) = 0.0_dp
                DO i = 1, qtb_therm%nf
@@ -248,7 +242,6 @@ CONTAINS
             DEALLOCATE (qtb_therm%cpt)
             DEALLOCATE (qtb_therm%step)
             DEALLOCATE (qtb_therm%rng_status)
-            CALL delete_rng_stream(qtb_therm%gaussian_rng_stream)
             DEALLOCATE (qtb_therm)
          END IF
       END IF
@@ -435,15 +428,14 @@ CONTAINS
       ELSE
          !update the rng status
          DO i = 1, qtb_therm%nf
-            CALL dump_rng_stream(rng_stream=qtb_therm%gaussian_rng_stream, &
-                                 rng_record=qtb_therm%rng_status(i))
+            CALL qtb_therm%gaussian_rng_stream%dump(qtb_therm%rng_status(i))
          END DO
          !if no restart then initialize random numbers from scratch
          qtb_therm%cpt = 0
          DO idim = 1, ndim
             DO ibead = 1, p
                DO i = 1, nf
-                  qtb_therm%r(i, ibead, idim) = next_random_number(qtb_therm%gaussian_rng_stream)
+                  qtb_therm%r(i, ibead, idim) = qtb_therm%gaussian_rng_stream%next()
                END DO
             END DO
          END DO
@@ -486,14 +478,13 @@ CONTAINS
          qtb_therm%cpt = 0
          !update the rng status
          DO i = 1, qtb_therm%nf
-            CALL dump_rng_stream(rng_stream=qtb_therm%gaussian_rng_stream, &
-                                 rng_record=qtb_therm%rng_status(i))
+            CALL qtb_therm%gaussian_rng_stream%dump(qtb_therm%rng_status(i))
          END DO
          !first random numbers initialized from scratch
          DO idim = 1, pint_env%ndim
             DO ibead = 1, pint_env%p
                DO i = 1, qtb_therm%nf
-                  qtb_therm%r(i, ibead, idim) = next_random_number(qtb_therm%gaussian_rng_stream)
+                  qtb_therm%r(i, ibead, idim) = qtb_therm%gaussian_rng_stream%next()
                END DO
             END DO
          END DO
@@ -518,15 +509,14 @@ CONTAINS
                   DO i = 1, qtb_therm%nf - 1
                      qtb_therm%rng_status(i) = qtb_therm%rng_status(i + 1)
                   END DO
-                  CALL dump_rng_stream(rng_stream=qtb_therm%gaussian_rng_stream, &
-                                       rng_record=qtb_therm%rng_status(qtb_therm%nf))
+                  CALL qtb_therm%gaussian_rng_stream%dump(qtb_therm%rng_status(qtb_therm%nf))
                END IF
                DO idim = 1, pint_env%ndim
                   !update random numbers
                   DO i = 1, qtb_therm%nf - 1
                      qtb_therm%r(i, ibead, idim) = qtb_therm%r(i + 1, ibead, idim)
                   END DO
-                  qtb_therm%r(qtb_therm%nf, ibead, idim) = next_random_number(qtb_therm%gaussian_rng_stream)
+                  qtb_therm%r(qtb_therm%nf, ibead, idim) = qtb_therm%gaussian_rng_stream%next()
                END DO
                qtb_therm%cpt(ibead) = 0
             END IF
diff --git a/src/motion/pint_types.F b/src/motion/pint_types.F
index 63a54a88d1..40c4008e40 100644
--- a/src/motion/pint_types.F
+++ b/src/motion/pint_types.F
@@ -121,7 +121,7 @@ MODULE pint_types
       TYPE(section_vals_type), POINTER :: input
       TYPE(staging_env_type), POINTER :: staging_env
       TYPE(normalmode_env_type), POINTER :: normalmode_env
-      TYPE(rng_stream_type), POINTER :: randomG
+      TYPE(rng_stream_type) :: randomG
       TYPE(gle_type), POINTER        :: gle
       REAL(KIND=dp), DIMENSION(e_num_ids) :: energy
       REAL(KIND=dp), DIMENSION(:), POINTER :: mass, e_pot_bead
@@ -206,7 +206,7 @@ MODULE pint_types
       REAL(KIND=dp), DIMENSION(:), POINTER       :: c2
       REAL(KIND=dp), DIMENSION(:), POINTER       :: g_fric
       REAL(KIND=dp), DIMENSION(:, :), POINTER    :: massfact
-      TYPE(rng_stream_type), POINTER             :: gaussian_rng_stream
+      TYPE(rng_stream_type)                      :: gaussian_rng_stream
    END TYPE pile_therm_type
 
    ! ***************************************************************************
@@ -240,7 +240,7 @@ MODULE pint_types
       REAL(KIND=dp), DIMENSION(:, :), POINTER :: temp1
       REAL(KIND=dp), DIMENSION(:, :), POINTER :: temp2
       REAL(KIND=dp), DIMENSION(:, :), POINTER :: sqrtmass
-      TYPE(rng_stream_type), POINTER       :: gaussian_rng_stream
+      TYPE(rng_stream_type)                :: gaussian_rng_stream
    END TYPE piglet_therm_type
 
    ! ***************************************************************************
@@ -283,7 +283,7 @@ MODULE pint_types
       INTEGER                                    :: fp
       INTEGER                                    :: nf
       REAL(KIND=dp)                              :: thermostat_energy
-      TYPE(rng_stream_type), POINTER             :: gaussian_rng_stream
+      TYPE(rng_stream_type)                      :: gaussian_rng_stream
       CHARACTER(LEN=rng_record_length), DIMENSION(:), POINTER :: rng_status
    END TYPE qtb_therm_type
 
diff --git a/src/motion/thermostat/al_system_dynamics.F b/src/motion/thermostat/al_system_dynamics.F
index 2c7f2519cc..344ae06055 100644
--- a/src/motion/thermostat/al_system_dynamics.F
+++ b/src/motion/thermostat/al_system_dynamics.F
@@ -19,8 +19,6 @@ MODULE al_system_dynamics
    USE molecule_kind_types,             ONLY: molecule_kind_type
    USE molecule_types,                  ONLY: get_molecule,&
                                               molecule_type
-   USE parallel_rng_types,              ONLY: next_random_number,&
-                                              rng_stream_type
    USE particle_types,                  ONLY: particle_type
    USE thermostat_utils,                ONLY: ke_region_particles,&
                                               vel_rescale_particles
@@ -190,7 +188,6 @@ CONTAINS
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: w(:, :)
       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
       TYPE(molecule_type), POINTER                       :: molecule
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
 
       present_vel = PRESENT(vel)
 
@@ -222,10 +219,9 @@ CONTAINS
          CPASSERT(check)
          DO iparticle_local = 1, nparticle_local
             ipart = local_particles%list(iparticle_kind)%array(iparticle_local)
-            rng_stream => local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream
-            w(1, ipart) = next_random_number(rng_stream, variance=1.0_dp)
-            w(2, ipart) = next_random_number(rng_stream, variance=1.0_dp)
-            w(3, ipart) = next_random_number(rng_stream, variance=1.0_dp)
+            w(1, ipart) = local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream%next(variance=1.0_dp)
+            w(2, ipart) = local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream%next(variance=1.0_dp)
+            w(3, ipart) = local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream%next(variance=1.0_dp)
          END DO
       END DO
 
diff --git a/src/motion/thermostat/csvr_system_init.F b/src/motion/thermostat/csvr_system_init.F
index 4a19e4eedb..14ae4bf251 100644
--- a/src/motion/thermostat/csvr_system_init.F
+++ b/src/motion/thermostat/csvr_system_init.F
@@ -21,8 +21,8 @@ MODULE csvr_system_init
    USE molecule_kind_types,             ONLY: molecule_kind_type
    USE molecule_types,                  ONLY: global_constraint_type,&
                                               molecule_type
-   USE parallel_rng_types,              ONLY: read_rng_stream,&
-                                              rng_record_length
+   USE parallel_rng_types,              ONLY: rng_record_length,&
+                                              rng_stream_type_from_record
    USE simpar_types,                    ONLY: simpar_type
    USE thermostat_types,                ONLY: thermostat_info_type
 #include "../../base/base_uses.f90"
@@ -182,8 +182,7 @@ CONTAINS
                my_index = csvr%map_info%index(i)
                CALL section_vals_val_get(section_vals=work_section, keyword_name="_DEFAULT_KEYWORD_", &
                                          i_rep_val=my_index, c_val=rng_record)
-               CALL read_rng_stream(rng_stream=csvr%nvt(i)%gaussian_rng_stream, &
-                                    rng_record=rng_record)
+               csvr%nvt(i)%gaussian_rng_stream = rng_stream_type_from_record(rng_record)
             END DO
          ELSE
             CALL cp_abort(__LOCATION__, &
diff --git a/src/motion/thermostat/extended_system_init.F b/src/motion/thermostat/extended_system_init.F
index 7e31a4bd2d..42954ffbe2 100644
--- a/src/motion/thermostat/extended_system_init.F
+++ b/src/motion/thermostat/extended_system_init.F
@@ -42,7 +42,6 @@ MODULE extended_system_init
    USE molecule_kind_types,             ONLY: molecule_kind_type
    USE molecule_types,                  ONLY: global_constraint_type,&
                                               molecule_type
-   USE parallel_rng_types,              ONLY: next_random_number
    USE simpar_types,                    ONLY: simpar_type
    USE thermostat_types,                ONLY: thermostat_info_type
    USE thermostat_utils,                ONLY: get_nhc_energies
@@ -750,7 +749,7 @@ CONTAINS
          ! Map deterministically determined random number to nhc % v
          DO i = 1, nhc%loc_num_nhc
             DO j = 1, nhc%nhc_len
-               nhc%nvt(j, i)%v = next_random_number(globenv%gaussian_rng_stream)
+               nhc%nvt(j, i)%v = globenv%gaussian_rng_stream%next()
             END DO
          END DO
 
@@ -787,7 +786,7 @@ CONTAINS
 
       CASE DEFAULT
          DO i = 1, tot_rn
-            array_of_rn(i) = next_random_number(globenv%gaussian_rng_stream)
+            array_of_rn(i) = globenv%gaussian_rng_stream%next()
          END DO
          ! Map deterministically determined random number to nhc % v
          DO i = 1, nhc%loc_num_nhc
@@ -884,7 +883,7 @@ CONTAINS
       ! initializing velocities
       DO i = 1, SIZE(npt, 1)
          DO j = i, SIZE(npt, 2)
-            v = next_random_number(globenv%gaussian_rng_stream)
+            v = globenv%gaussian_rng_stream%next()
             ! Symmetrizing the initial barostat velocities to ensure
             ! no rotation of the cell under NPT_F
             npt(j, i)%v = v
diff --git a/src/motion/thermostat/gle_system_dynamics.F b/src/motion/thermostat/gle_system_dynamics.F
index 7fa00ae115..ba7c7d2e70 100644
--- a/src/motion/thermostat/gle_system_dynamics.F
+++ b/src/motion/thermostat/gle_system_dynamics.F
@@ -28,10 +28,8 @@ MODULE gle_system_dynamics
    USE molecule_kind_types,             ONLY: molecule_kind_type
    USE molecule_types,                  ONLY: global_constraint_type,&
                                               molecule_type
-   USE parallel_rng_types,              ONLY: next_random_number,&
-                                              read_rng_stream,&
-                                              rng_record_length,&
-                                              rng_stream_type
+   USE parallel_rng_types,              ONLY: rng_record_length,&
+                                              rng_stream_type_from_record
    USE particle_types,                  ONLY: particle_type
    USE simpar_types,                    ONLY: simpar_type
    USE thermostat_mapping,              ONLY: thermostat_mapping_region
@@ -93,12 +91,10 @@ CONTAINS
       REAL(dp)                                           :: alpha, beta, rr
       REAL(dp), DIMENSION(:, :), POINTER                 :: a_mat, e_tmp, h_tmp, s_tmp
       TYPE(map_info_type), POINTER                       :: map_info
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
 
       CALL timeset(routineN, handle)
       my_shell_adiabatic = .FALSE.
       IF (PRESENT(shell_adiabatic)) my_shell_adiabatic = shell_adiabatic
-      NULLIFY (rng_stream)
       present_vel = PRESENT(vel)
       ndim = gle%ndim
       ALLOCATE (s_tmp(ndim, gle%loc_num_gle))
@@ -122,12 +118,11 @@ CONTAINS
          IF (gle%nvt(ideg)%nkt == 0.0_dp) CYCLE
          gle%nvt(ideg)%s(1) = map_info%s_kin(imap)
          s_tmp(1, imap) = map_info%s_kin(imap)
-         rng_stream => gle%nvt(ideg)%gaussian_rng_stream
-         rr = next_random_number(rng_stream)
+         rr = gle%nvt(ideg)%gaussian_rng_stream%next()
          e_tmp(1, imap) = rr
          DO iadd = 2, ndim
             s_tmp(iadd, imap) = gle%nvt(ideg)%s(iadd)
-            rr = next_random_number(rng_stream)
+            rr = gle%nvt(ideg)%gaussian_rng_stream%next()
             e_tmp(iadd, imap) = rr
          END DO
       END DO
@@ -502,8 +497,7 @@ CONTAINS
                   ind = map_info%index(i)
                   CALL section_vals_val_get(section_vals=work_section, keyword_name="_DEFAULT_KEYWORD_", &
                                             i_rep_val=ind, c_val=rng_record)
-                  CALL read_rng_stream(rng_stream=gle%nvt(i)%gaussian_rng_stream, &
-                                       rng_record=rng_record)
+                  gle%nvt(i)%gaussian_rng_stream = rng_stream_type_from_record(rng_record)
                END DO
             ELSE
                CALL cp_abort(__LOCATION__, &
@@ -528,14 +522,12 @@ CONTAINS
 
       INTEGER                                            :: i, j
       REAL(dp)                                           :: rr(gle%ndim), cc(gle%ndim, gle%ndim)
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
 
       CALL gle_cholesky_stab(gle%c_mat, cc, gle%ndim)
       DO i = 1, gle%loc_num_gle
-         rng_stream => gle%nvt(i)%gaussian_rng_stream
          DO j = 1, gle%ndim
             ! here s should be properly initialized, when it is not read from restart
-            rr(j) = next_random_number(rng_stream)
+            rr(j) = gle%nvt(i)%gaussian_rng_stream%next()
          END DO
          gle%nvt(i)%s = MATMUL(cc, rr)
       END DO
diff --git a/src/motion/wiener_process.F b/src/motion/wiener_process.F
index b169626d95..11d8a7af02 100644
--- a/src/motion/wiener_process.F
+++ b/src/motion/wiener_process.F
@@ -29,10 +29,10 @@ MODULE wiener_process
                                               need_per_atom_wiener_process
    USE metadynamics_types,              ONLY: meta_env_type
    USE parallel_rng_types,              ONLY: GAUSSIAN,&
-                                              create_rng_stream,&
                                               next_rng_seed,&
-                                              read_rng_stream,&
-                                              rng_record_length
+                                              rng_record_length,&
+                                              rng_stream_type,&
+                                              rng_stream_type_from_record
    USE particle_list_types,             ONLY: particle_list_type
    USE simpar_types,                    ONLY: simpar_type
    USE string_utilities,                ONLY: compress
@@ -107,7 +107,7 @@ CONTAINS
             nparticle_local = local_particles%n_el(iparticle_kind)
             ALLOCATE (local_particles%local_particle_set(iparticle_kind)%rng(nparticle_local))
             DO iparticle_local = 1, nparticle_local
-               NULLIFY (local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream)
+               ALLOCATE (local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream)
             END DO
          END DO
 
@@ -132,9 +132,9 @@ CONTAINS
                iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
                WRITE (UNIT=name, FMT="(A,I8)") "Wiener process for particle", iparticle
                CALL compress(name)
-               CALL create_rng_stream(rng_stream=local_particles%local_particle_set(iparticle_kind)% &
-                                      rng(iparticle_local)%stream, name=name, distribution_type=GAUSSIAN, &
-                                      extended_precision=.TRUE., seed=seed(:, :, iparticle))
+               local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)% &
+                  stream = rng_stream_type(name=name, distribution_type=GAUSSIAN, &
+                                           extended_precision=.TRUE., seed=seed(:, :, iparticle))
             END DO
          END DO
 
@@ -190,10 +190,8 @@ CONTAINS
                                                keyword_name="_DEFAULT_KEYWORD_", &
                                                i_rep_val=iparticle, &
                                                c_val=rng_record)
-                     CALL read_rng_stream(rng_stream=distribution_1d% &
-                                          local_particle_set(iparticle_kind)% &
-                                          rng(iparticle_local)%stream, &
-                                          rng_record=rng_record)
+                     distribution_1d%local_particle_set(iparticle_kind)%rng(iparticle_local)% &
+                        stream = rng_stream_type_from_record(rng_record)
                   END IF
                END DO
             END DO
@@ -227,10 +225,6 @@ CONTAINS
 
       initial_seed = next_rng_seed()
 
-      DO i_c = 1, meta_env%n_colvar
-         NULLIFY (meta_env%rng(i_c)%stream)
-      END DO
-
       ! Each process generates all seeds. The seed generation should be
       ! quite fast and in this way a broadcast is avoided.
 
@@ -248,8 +242,8 @@ CONTAINS
       DO i_c = 1, meta_env%n_colvar
          WRITE (UNIT=name, FMT="(A,I8)") "Wiener process for COLVAR", i_c
          CALL compress(name)
-         CALL create_rng_stream(rng_stream=meta_env%rng(i_c)%stream, name=name, distribution_type=GAUSSIAN, &
-                                extended_precision=.TRUE., seed=seed(:, :, i_c))
+         meta_env%rng(i_c) = rng_stream_type(name=name, distribution_type=GAUSSIAN, &
+                                             extended_precision=.TRUE., seed=seed(:, :, i_c))
       END DO
       DEALLOCATE (seed)
 
diff --git a/src/optimize_input.F b/src/optimize_input.F
index a8186f3689..243a36b2ef 100644
--- a/src/optimize_input.F
+++ b/src/optimize_input.F
@@ -45,9 +45,6 @@ MODULE optimize_input
                                               mp_environ,&
                                               mp_sum
    USE parallel_rng_types,              ONLY: UNIFORM,&
-                                              create_rng_stream,&
-                                              delete_rng_stream,&
-                                              next_random_number,&
                                               rng_stream_type
    USE physcon,                         ONLY: bohr
    USE powell,                          ONLY: opt_state_type,&
@@ -116,7 +113,7 @@ CONTAINS
       INTEGER                                            :: handle, i_var
       REAL(KIND=dp)                                      :: random_number, seed(3, 2)
       TYPE(oi_env_type)                                  :: oi_env
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), ALLOCATABLE                 :: rng_stream
 
       CALL timeset(routineN, handle)
 
@@ -126,18 +123,16 @@ CONTAINS
 
       ! if we have been asked to randomize the variables, we do this.
       IF (oi_env%randomize_variables .NE. 0.0_dp) THEN
-         NULLIFY (rng_stream)
          seed = REAL(oi_env%seed, KIND=dp)
-         CALL create_rng_stream(rng_stream, "run_optimize_input", distribution_type=UNIFORM, seed=seed)
+         rng_stream = rng_stream_type("run_optimize_input", distribution_type=UNIFORM, seed=seed)
          DO i_var = 1, SIZE(oi_env%variables, 1)
             IF (.NOT. oi_env%variables(i_var)%fixed) THEN
                ! change with a random percentage the variable
-               random_number = next_random_number(rng_stream)
+               random_number = rng_stream%next()
                oi_env%variables(i_var)%value = oi_env%variables(i_var)%value* &
                                                (1.0_dp + (2*random_number - 1.0_dp)*oi_env%randomize_variables/100.0_dp)
             ENDIF
          ENDDO
-         CALL delete_rng_stream(rng_stream)
       ENDIF
 
       ! proceed to actual methods
diff --git a/src/pao_ml_neuralnet.F b/src/pao_ml_neuralnet.F
index 8d21f64d94..596dae32f5 100644
--- a/src/pao_ml_neuralnet.F
+++ b/src/pao_ml_neuralnet.F
@@ -11,10 +11,7 @@ MODULE pao_ml_neuralnet
    USE kinds,                           ONLY: dp
    USE pao_types,                       ONLY: pao_env_type,&
                                               training_matrix_type
-   USE parallel_rng_types,              ONLY: create_rng_stream,&
-                                              delete_rng_stream,&
-                                              next_random_number,&
-                                              rng_stream_type
+   USE parallel_rng_types,              ONLY: rng_stream_type
 #include "./base/base_uses.f90"
 
    IMPLICIT NONE
@@ -126,11 +123,9 @@ CONTAINS
       REAL(dp)                                           :: bak, eps, error, error1, error2, num_grad
       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: prediction
       REAL(dp), ALLOCATABLE, DIMENSION(:, :, :)          :: gradient
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type)                              :: rng_stream
       TYPE(training_matrix_type), POINTER                :: training_matrix
 
-      NULLIFY (rng_stream)
-
       ! TODO this could be parallelized over ranks
       DO ikind = 1, SIZE(pao%ml_training_matrices)
          training_matrix => pao%ml_training_matrices(ikind)
@@ -150,15 +145,14 @@ CONTAINS
          ALLOCATE (training_matrix%NN(nlayers, width, width))
 
          ! initialize network with random numbers from -1.0 ... +1.0
-         CALL create_rng_stream(rng_stream, name="pao_nn")
+         rng_stream = rng_stream_type(name="pao_nn")
          DO ilayer = 1, nlayers
          DO i = 1, width
          DO j = 1, width
-            training_matrix%NN(ilayer, i, j) = -1.0_dp + 2.0_dp*next_random_number(rng_stream)
+            training_matrix%NN(ilayer, i, j) = -1.0_dp + 2.0_dp*rng_stream%next()
          ENDDO
          ENDDO
          ENDDO
-         CALL delete_rng_stream(rng_stream)
 
          ! train the network using backpropagation
          ALLOCATE (gradient(nlayers, width, width))
diff --git a/src/start/input_cp2k.F b/src/start/input_cp2k.F
index 988c85f5ef..3f2e181f0b 100644
--- a/src/start/input_cp2k.F
+++ b/src/start/input_cp2k.F
@@ -163,7 +163,7 @@ CONTAINS
 
       CALL section_create(section, __LOCATION__, name="TEST", &
                           description="Tests to perform on the supported libraries.", &
-                          n_keywords=7, n_subsections=0, repeats=.FALSE.)
+                          n_keywords=6, n_subsections=0, repeats=.FALSE.)
 
       NULLIFY (keyword, print_key)
       CALL keyword_create(keyword, __LOCATION__, name="MEMORY", &
@@ -227,12 +227,6 @@ CONTAINS
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
-      CALL keyword_create(keyword, __LOCATION__, name="RANDOM_NUMBER_GENERATOR", variants=(/"rng"/), &
-                          description=" Tests the parallel random number generator (RNG)", &
-                          usage="rng 1000000", default_i_val=0)
-      CALL section_add_keyword(section, keyword)
-      CALL keyword_release(keyword)
-
       CALL keyword_create(keyword, __LOCATION__, name="MINIMAX", &
                           description="Tests validity of minimax coefficients for approximating 1/x "// &
                           "as a sum of exponentials. "// &
diff --git a/src/statistical_methods.F b/src/statistical_methods.F
index 8b69a86216..5854e52e31 100644
--- a/src/statistical_methods.F
+++ b/src/statistical_methods.F
@@ -15,7 +15,6 @@ MODULE statistical_methods
    USE cp_log_handling,                 ONLY: cp_logger_get_default_io_unit
    USE global_types,                    ONLY: global_environment_type
    USE kinds,                           ONLY: dp
-   USE parallel_rng_types,              ONLY: next_random_number
    USE util,                            ONLY: sort
 #include "./base/base_uses.f90"
 
@@ -456,12 +455,11 @@ CONTAINS
          NULLIFY (xdata)
          ALLOCATE (xdata(n))
          DO i = 1, 10
-            xdata(i) = 5.0_dp - REAL(i, KIND=dp)/2.0_dp + 0.1* &
-                       next_random_number(globenv%gaussian_rng_stream)
+            xdata(i) = 5.0_dp - REAL(i, KIND=dp)/2.0_dp + 0.1*globenv%gaussian_rng_stream%next()
             WRITE (3, *) xdata(i)
          END DO
          DO i = 11, n
-            xdata(i) = 0.1*next_random_number(globenv%gaussian_rng_stream)
+            xdata(i) = 0.1*globenv%gaussian_rng_stream%next()
          END DO
 
          ! Test for trend
diff --git a/src/swarm/glbopt_mincrawl.F b/src/swarm/glbopt_mincrawl.F
index 569f88ba2d..c968d6ec2e 100644
--- a/src/swarm/glbopt_mincrawl.F
+++ b/src/swarm/glbopt_mincrawl.F
@@ -26,10 +26,7 @@ MODULE glbopt_mincrawl
                                               section_vals_val_get
    USE kinds,                           ONLY: default_string_length,&
                                               dp
-   USE parallel_rng_types,              ONLY: create_rng_stream,&
-                                              delete_rng_stream,&
-                                              next_random_number,&
-                                              rng_stream_type
+   USE parallel_rng_types,              ONLY: rng_stream_type
    USE particle_methods,                ONLY: write_particle_coordinates
    USE particle_types,                  ONLY: particle_type
    USE physcon,                         ONLY: kelvin
@@ -88,7 +85,7 @@ MODULE glbopt_mincrawl
       INTEGER                                             :: iw = 0
       INTEGER                                             :: minima_traj_unit = 0
       TYPE(section_vals_type), POINTER                    :: mincrawl_section => Null()
-      TYPE(rng_stream_type), POINTER                       :: rng_stream => Null()
+      TYPE(rng_stream_type)                               :: rng_stream
       TYPE(particle_type), DIMENSION(:), POINTER          :: particle_set => Null()
    END TYPE mincrawl_type
 
@@ -154,7 +151,7 @@ CONTAINS
          this%tempdist_init(i) = 1.0/(1.0 + EXP((this%tempstep_init - i)/this%tempdist_init_width))
       ENDDO
 
-      CALL create_rng_stream(this%rng_stream, name="mincrawl")
+      this%rng_stream = rng_stream_type(name="mincrawl")
    END SUBROUTINE mincrawl_init
 
 ! **************************************************************************************************
@@ -274,10 +271,10 @@ CONTAINS
       REAL(KIND=dp)                                      :: a, r
 
       DO
-         r = next_random_number(this%rng_stream)
+         r = this%rng_stream%next()
          step = INT(r*SIZE(minima%tempdist)) + 1
          a = 1.0 - 2.0*ABS(minima%tempdist(step) - 0.5)
-         r = next_random_number(this%rng_stream)
+         r = this%rng_stream%next()
          IF (r < a) EXIT
       END DO
 
@@ -516,7 +513,6 @@ CONTAINS
                                         this%mincrawl_section, "MINIMA_TRAJECTORY")
 
       CALL history_finalize(this%history)
-      CALL delete_rng_stream(this%rng_stream)
    END SUBROUTINE mincrawl_finalize
 
 END MODULE glbopt_mincrawl
diff --git a/src/swarm/glbopt_worker.F b/src/swarm/glbopt_worker.F
index f188a0733b..658f44a4d8 100644
--- a/src/swarm/glbopt_worker.F
+++ b/src/swarm/glbopt_worker.F
@@ -32,7 +32,6 @@ MODULE glbopt_worker
    USE md_run,                          ONLY: qs_mol_dyn
    USE mdctrl_types,                    ONLY: glbopt_mdctrl_data_type,&
                                               mdctrl_type
-   USE parallel_rng_types,              ONLY: reset_to_next_rng_substream
    USE physcon,                         ONLY: angstrom,&
                                               kelvin
    USE swarm_message,                   ONLY: swarm_message_add,&
@@ -116,7 +115,7 @@ CONTAINS
 
       ! We want different random-number-streams for each worker
       DO i = 1, worker_id
-         CALL reset_to_next_rng_substream(worker%globenv%gaussian_rng_stream)
+         CALL worker%globenv%gaussian_rng_stream%reset_to_next_substream()
       END DO
 
       CALL cp_subsys_get(worker%subsys, natom=worker%n_atoms)
@@ -385,4 +384,3 @@ CONTAINS
    END SUBROUTINE glbopt_worker_finalize
 
 END MODULE glbopt_worker
-
diff --git a/src/tmc/tmc_calculations.F b/src/tmc/tmc_calculations.F
index d7c234e130..a4b19c8a70 100644
--- a/src/tmc/tmc_calculations.F
+++ b/src/tmc/tmc_calculations.F
@@ -22,10 +22,7 @@ MODULE tmc_calculations
                                               set_cell
    USE kinds,                           ONLY: dp
    USE mathconstants,                   ONLY: pi
-   USE parallel_rng_types,              ONLY: get_rng_stream,&
-                                              next_random_number,&
-                                              rng_stream_type,&
-                                              set_rng_stream
+   USE parallel_rng_types,              ONLY: rng_stream_type
    USE physcon,                         ONLY: boltzmann,&
                                               joule
    USE tmc_move_types,                  ONLY: mv_type_MD
@@ -123,7 +120,7 @@ CONTAINS
       END SELECT
 
       ! ---     wait a bit
-      rnd = next_random_number(rng_stream=tmc_env%rng_stream)
+      rnd = tmc_env%rng_stream%next()
       !rnd = 0.5
 !TODO    IF(worker_random_wait.AND.exact_approx_pot)THEN
 !      CALL SYSTEM_CLOCK(time0, time_rate, time_max)
@@ -361,7 +358,7 @@ CONTAINS
       REAL(KIND=dp), DIMENSION(:), POINTER               :: vel
       TYPE(tmc_atom_type), DIMENSION(:), POINTER         :: atoms
       REAL(KIND=dp)                                      :: temerature
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
       REAL(KIND=dp), DIMENSION(3, 2, 3)                  :: rnd_seed
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'init_vel', routineP = moduleN//':'//routineN
@@ -374,19 +371,17 @@ CONTAINS
       CPASSERT(ASSOCIATED(vel))
       CPASSERT(ASSOCIATED(atoms))
 
-      CALL set_rng_stream(rng_stream=rng_stream, bg=rnd_seed(:, :, 1), &
-                          cg=rnd_seed(:, :, 2), ig=rnd_seed(:, :, 3))
+      CALL rng_stream%set(bg=rnd_seed(:, :, 1), cg=rnd_seed(:, :, 2), ig=rnd_seed(:, :, 3))
       DO i = 1, SIZE(vel)
-         rnd1 = next_random_number(rng_stream)
-         rnd2 = next_random_number(rng_stream)
+         rnd1 = rng_stream%next()
+         rnd2 = rng_stream%next()
 
          mass_tmp = atoms(INT(i/REAL(3, KIND=dp)) + 1)%mass
 
          vel(i) = SQRT(-2.0_dp*LOG(rnd1))*COS(2.0_dp*PI*rnd2)* &
                   SQRT(kB*temerature/mass_tmp)
       END DO
-      CALL get_rng_stream(rng_stream=rng_stream, bg=rnd_seed(:, :, 1), &
-                          cg=rnd_seed(:, :, 2), ig=rnd_seed(:, :, 3))
+      CALL rng_stream%get(bg=rnd_seed(:, :, 1), cg=rnd_seed(:, :, 2), ig=rnd_seed(:, :, 3))
 
    END SUBROUTINE init_vel
 
diff --git a/src/tmc/tmc_moves.F b/src/tmc/tmc_moves.F
index 9e7fc3d0c0..633ed313f8 100644
--- a/src/tmc/tmc_moves.F
+++ b/src/tmc/tmc_moves.F
@@ -19,10 +19,7 @@ MODULE tmc_moves
    USE mathconstants,                   ONLY: pi
    USE mathlib,                         ONLY: dihedral_angle,&
                                               rotate_vector
-   USE parallel_rng_types,              ONLY: get_rng_stream,&
-                                              next_random_number,&
-                                              rng_stream_type,&
-                                              set_rng_stream
+   USE parallel_rng_types,              ONLY: rng_stream_type
    USE physcon,                         ONLY: boltzmann,&
                                               joule
    USE tmc_calculations,                ONLY: center_of_mass,&
@@ -70,7 +67,7 @@ CONTAINS
                          new_subbox, move_rejected)
       TYPE(tmc_param_type), POINTER                      :: tmc_params
       TYPE(tmc_move_type), POINTER                       :: move_types
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
       TYPE(tree_type), POINTER                           :: elem
       INTEGER                                            :: mv_conf
       LOGICAL                                            :: new_subbox, move_rejected
@@ -87,12 +84,11 @@ CONTAINS
 
       CPASSERT(ASSOCIATED(tmc_params))
       CPASSERT(ASSOCIATED(move_types))
-      CPASSERT(ASSOCIATED(rng_stream))
       CPASSERT(ASSOCIATED(elem))
 
       move_rejected = .FALSE.
 
-      CALL set_rng_stream(rng_stream=rng_stream, bg=elem%rng_seed(:, :, 1), &
+      CALL rng_stream%set(bg=elem%rng_seed(:, :, 1), &
                           cg=elem%rng_seed(:, :, 2), ig=elem%rng_seed(:, :, 3))
 
       IF (new_subbox) THEN
@@ -137,7 +133,7 @@ CONTAINS
             IF (tmc_params%nr_elem_mv .EQ. 0) THEN
                ind = (i - 1)*(tmc_params%dim_per_elem) + 1
             ELSE
-               rnd = next_random_number(rng_stream) ! next random number
+               rnd = rng_stream%next()
                ind = tmc_params%dim_per_elem* &
                      INT(rnd*(SIZE(elem%pos)/tmc_params%dim_per_elem)) + 1
             END IF
@@ -145,7 +141,7 @@ CONTAINS
             IF (elem%elem_stat(ind) .EQ. status_ok) THEN
                ! displace atom
                DO d = 0, tmc_params%dim_per_elem - 1
-                  rnd = next_random_number(rng_stream) ! next random number
+                  rnd = rng_stream%next()
                   elem%pos(ind + d) = elem%pos(ind + d) + (rnd - 0.5)*2.0* &
                                       move_types%mv_size(mv_type_atom_trans, mv_conf)
                END DO
@@ -197,7 +193,7 @@ CONTAINS
                IF (tmc_params%nr_elem_mv .EQ. 0) THEN
                   m = counter
                ELSE
-                  rnd = next_random_number(rng_stream) ! next random number
+                  rnd = rng_stream%next()
                   m = INT(rnd*nr_molec) + 1
                END IF
                CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, &
@@ -210,7 +206,7 @@ CONTAINS
                IF (mol_in_sb(m) .EQ. status_ok) THEN
                   ! calculate displacement
                   DO d = 1, tmc_params%dim_per_elem
-                     rnd = next_random_number(rng_stream) ! next random number
+                     rnd = rng_stream%next()
                      direction(d) = (rnd - 0.5)*2.0_dp*move_types%mv_size( &
                                     mv_type_mol_trans, mv_conf)
                   END DO
@@ -271,7 +267,7 @@ CONTAINS
                IF (tmc_params%nr_elem_mv .EQ. 0) THEN
                   m = counter
                ELSE
-                  rnd = next_random_number(rng_stream) ! next random number
+                  rnd = rng_stream%next()
                   m = INT(rnd*nr_molec) + 1
                END IF
                CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, &
@@ -352,7 +348,7 @@ CONTAINS
                        cp_to_string(elem%move_type))
       END SELECT
 
-      CALL get_rng_stream(rng_stream=rng_stream, bg=elem%rng_seed(:, :, 1), &
+      CALL rng_stream%get(bg=elem%rng_seed(:, :, 1), &
                           cg=elem%rng_seed(:, :, 2), ig=elem%rng_seed(:, :, 3))
 
    END SUBROUTINE change_pos
@@ -472,7 +468,7 @@ CONTAINS
    SUBROUTINE elements_in_new_subbox(tmc_params, rng_stream, elem, &
                                      nr_of_sub_box_elements)
       TYPE(tmc_param_type), POINTER                      :: tmc_params
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
       TYPE(tree_type), POINTER                           :: elem
       INTEGER, INTENT(OUT)                               :: nr_of_sub_box_elements
 
@@ -501,19 +497,17 @@ CONTAINS
          ALLOCATE (atom_tmp(tmc_params%dim_per_elem))
          nr_of_sub_box_elements = 0
          ! -- define the center of the sub box
-         CALL set_rng_stream(rng_stream=rng_stream, &
-                             bg=elem%rng_seed(:, :, 1), cg=elem%rng_seed(:, :, 2), &
+         CALL rng_stream%set(bg=elem%rng_seed(:, :, 1), cg=elem%rng_seed(:, :, 2), &
                              ig=elem%rng_seed(:, :, 3))
 
          CALL get_cell(cell=tmc_params%cell, abc=box_size)
          DO i = 1, SIZE(tmc_params%sub_box_size)
-            rnd = next_random_number(rng_stream) ! next random number
+            rnd = rng_stream%next()
             center_of_sub_box(i) = rnd*box_size(i)
          END DO
          elem%subbox_center(:) = center_of_sub_box(:)
 
-         CALL get_rng_stream(rng_stream=rng_stream, &
-                             bg=elem%rng_seed(:, :, 1), cg=elem%rng_seed(:, :, 2), &
+         CALL rng_stream%get(bg=elem%rng_seed(:, :, 1), cg=elem%rng_seed(:, :, 2), &
                              ig=elem%rng_seed(:, :, 3))
 
          ! check all elements if they are in subbox
@@ -552,7 +546,7 @@ CONTAINS
       INTEGER                                            :: ind_start, ind_end
       REAL(KIND=dp)                                      :: max_angle
       TYPE(tmc_move_type), POINTER                       :: move_types
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
       INTEGER                                            :: dim_per_elem
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'do_mol_rot', routineP = moduleN//':'//routineN
@@ -571,11 +565,11 @@ CONTAINS
       CPASSERT(ASSOCIATED(move_types))
 
       ! calculate rotation matrix (using quanternions)
-      rnd = next_random_number(rng_stream) ! next random number
+      rnd = rng_stream%next()
       a1 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf)
-      rnd = next_random_number(rng_stream) ! next random number
+      rnd = rng_stream%next()
       a2 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf)
-      rnd = next_random_number(rng_stream) ! next random number
+      rnd = rng_stream%next()
       a3 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf)
       q0 = COS(a2/2)*COS((a1 + a3)/2.0_dp)
       q1 = SIN(a2/2)*COS((a1 - a3)/2.0_dp)
@@ -615,7 +609,7 @@ CONTAINS
       TYPE(tmc_atom_type)                                :: atom_kind
       REAL(KIND=dp), INTENT(IN)                          :: phi, temp
       LOGICAL                                            :: rnd_sign_change
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'vel_change', routineP = moduleN//':'//routineN
 
@@ -624,12 +618,10 @@ CONTAINS
 
       kB = boltzmann/joule
 
-      CPASSERT(ASSOCIATED(rng_stream))
-
       !phi = move_types%mv_size(mv_type_MD,1) ! TODO parallel tempering move sizes for vel_change
       ! hence first producing a gaussian random number
-      rnd1 = next_random_number(rng_stream)
-      rnd2 = next_random_number(rng_stream)
+      rnd1 = rng_stream%next()
+      rnd2 = rng_stream%next()
 
       rnd_g = SQRT(-2.0_dp*LOG(rnd1))*COS(2.0_dp*PI*rnd2)
       !we can also produce a second one in the same step:
@@ -643,7 +635,7 @@ CONTAINS
       ! can be switched of using MD_vel_invert
       ! without still the balance condition should be fullfilled
 
-      rnd3 = next_random_number(rng_stream)
+      rnd3 = rng_stream%next()
       IF (rnd3 .GE. 0.5 .AND. rnd_sign_change) THEN
          d = -1
       ELSE
@@ -669,7 +661,7 @@ CONTAINS
                                                  tmc_params)
       TYPE(tree_type), POINTER                           :: elem
       LOGICAL                                            :: short_loop
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
       TYPE(tmc_param_type), POINTER                      :: tmc_params
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'search_and_do_proton_displace_loop', &
@@ -684,7 +676,6 @@ CONTAINS
       NULLIFY (mol_arr)
 
       CPASSERT(ASSOCIATED(elem))
-      CPASSERT(ASSOCIATED(rng_stream))
       CPASSERT(ASSOCIATED(tmc_params))
 
       ! start the timing
@@ -698,7 +689,7 @@ CONTAINS
       mol_arr(:) = -1
       donor_acceptor = not_selected
       ! select randomly if neighboring molecule is donor / acceptor
-      IF (next_random_number(rng_stream) .LT. 0.5_dp) THEN
+      IF (rng_stream%next() .LT. 0.5_dp) THEN
          donor_acceptor = proton_acceptor
       ELSE
          donor_acceptor = proton_donor
@@ -706,7 +697,7 @@ CONTAINS
 
       ! first step build loop
       !  select randomly one atom
-      rnd = next_random_number(rng_stream)
+      rnd = rng_stream%next()
       ! the randomly selected first atom
       mol = INT(rnd*nr_mol) + 1
       counter = counter + 1
@@ -773,7 +764,7 @@ CONTAINS
       TYPE(tree_type), POINTER                           :: elem
       INTEGER                                            :: mol, donor_acceptor
       TYPE(tmc_param_type), POINTER                      :: tmc_params
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'find_nearest_proton_acceptor_donator', &
          routineP = moduleN//':'//routineN
@@ -885,7 +876,7 @@ CONTAINS
       END IF
 
       ! select randomly the next neighboring molecule
-      rnd = next_random_number(rng_stream)
+      rnd = rng_stream%next()
       ! the randomly selected atom: return value!
       mol_tmp = neighbor_mol(INT(rnd*SIZE(neighbor_mol(:))) + 1)
       mol = mol_tmp
@@ -1125,7 +1116,7 @@ CONTAINS
       TYPE(tree_type), POINTER                           :: conf
       INTEGER                                            :: T_ind
       TYPE(tmc_move_type), POINTER                       :: move_types
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
       TYPE(tmc_param_type), POINTER                      :: tmc_params
       LOGICAL                                            :: mv_cen_of_mass
 
@@ -1141,7 +1132,6 @@ CONTAINS
 
       CPASSERT(ASSOCIATED(conf))
       CPASSERT(ASSOCIATED(move_types))
-      CPASSERT(ASSOCIATED(rng_stream))
       CPASSERT(ASSOCIATED(tmc_params))
       CPASSERT(T_ind .GT. 0 .AND. T_ind .LE. tmc_params%nr_temp)
       CPASSERT(tmc_params%dim_per_elem .EQ. 3)
@@ -1163,15 +1153,15 @@ CONTAINS
          IF (tmc_params%v_isotropic) THEN
             CALL get_scaled_cell(cell=tmc_params%cell, box_scale=conf%box_scale, &
                                  abc=box_length_new, vol=vol)
-            rnd = next_random_number(rng_stream)
+            rnd = rng_stream%next()
             vol = vol + (rnd - 0.5_dp)*2.0_dp*move_types%mv_size(mv_type_volume_move, T_ind)
             box_length_new(:) = vol**(1/REAL(3, KIND=dp))
          ELSE
             CALL get_scaled_cell(cell=tmc_params%cell, box_scale=conf%box_scale, &
                                  abc=box_length_new, vol=vol)
-            rnd = next_random_number(rng_stream)
+            rnd = rng_stream%next()
             vol = vol + (rnd - 0.5_dp)*2.0_dp*move_types%mv_size(mv_type_volume_move, T_ind)
-            rnd = next_random_number(rng_stream)
+            rnd = rng_stream%next()
             dir = 1 + INT(rnd*3)
             box_length_new(dir) = 1.0_dp
             box_length_new(dir) = vol/PRODUCT(box_length_new(:))
@@ -1181,15 +1171,15 @@ CONTAINS
          ! increase / decrease box lenght in this direction
          ! l_n = l_o +- rnd * mv_size
          IF (tmc_params%v_isotropic) THEN
-            rnd = next_random_number(rng_stream)
+            rnd = rng_stream%next()
             box_length_new(:) = box_length_new(:) + &
                                 (rnd - 0.5_dp)*2.0_dp* &
                                 move_types%mv_size(mv_type_volume_move, T_ind)
          ELSE
             ! select a random direction
-            rnd = next_random_number(rng_stream)
+            rnd = rng_stream%next()
             dir = 1 + INT(rnd*3)
-            rnd = next_random_number(rng_stream)
+            rnd = rng_stream%next()
             box_length_new(dir) = box_length_new(dir) + &
                                   (rnd - 0.5_dp)*2.0_dp* &
                                   move_types%mv_size(mv_type_volume_move, T_ind)
@@ -1254,7 +1244,7 @@ CONTAINS
    SUBROUTINE swap_atoms(conf, move_types, rng_stream, tmc_params)
       TYPE(tree_type), POINTER                           :: conf
       TYPE(tmc_move_type), POINTER                       :: move_types
-      TYPE(rng_stream_type), POINTER                     :: rng_stream
+      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
       TYPE(tmc_param_type), POINTER                      :: tmc_params
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'swap_atoms', routineP = moduleN//':'//routineN
@@ -1265,7 +1255,6 @@ CONTAINS
 
       CPASSERT(ASSOCIATED(conf))
       CPASSERT(ASSOCIATED(move_types))
-      CPASSERT(ASSOCIATED(rng_stream))
       CPASSERT(ASSOCIATED(tmc_params))
       CPASSERT(ASSOCIATED(tmc_params%atoms))
 
@@ -1273,10 +1262,10 @@ CONTAINS
       atom_search_loop: DO
          ! select one atom randomly
          a_1 = INT(SIZE(conf%pos)/REAL(tmc_params%dim_per_elem, KIND=dp)* &
-                   next_random_number(rng_stream)) + 1
+                   rng_stream%next()) + 1
          ! select the second atom randomly
          a_2 = INT(SIZE(conf%pos)/REAL(tmc_params%dim_per_elem, KIND=dp)* &
-                   next_random_number(rng_stream)) + 1
+                   rng_stream%next()) + 1
          ! check if they have different kinds
          IF (tmc_params%atoms(a_1)%name .NE. tmc_params%atoms(a_2)%name) THEN
             ! if present, check if atoms have different type related to the specified table
diff --git a/src/tmc/tmc_setup.F b/src/tmc/tmc_setup.F
index 35e7cb4036..572f2f55ee 100644
--- a/src/tmc/tmc_setup.F
+++ b/src/tmc/tmc_setup.F
@@ -43,8 +43,7 @@ MODULE tmc_setup
                                               mp_comm_free,&
                                               mp_comm_split_direct
    USE parallel_rng_types,              ONLY: UNIFORM,&
-                                              create_rng_stream,&
-                                              delete_rng_stream
+                                              rng_stream_type
    USE physcon,                         ONLY: au2a => angstrom,&
                                               au2bar => bar
    USE tmc_analysis,                    ONLY: analysis_init,&
@@ -176,15 +175,15 @@ CONTAINS
                          tmc_env%m_env%rnd_init*10.0_dp, &
                          tmc_env%m_env%rnd_init*2.0_dp/), &
                        (/3, 2/))
-            CALL create_rng_stream(rng_stream=tmc_env%rng_stream, &
-                                   name="TMC_deterministic_rng_stream", &
-                                   seed=init_rng_seed(:, :), &
-                                   distribution_type=UNIFORM)
+            tmc_env%rng_stream = rng_stream_type( &
+                                 name="TMC_deterministic_rng_stream", &
+                                 seed=init_rng_seed(:, :), &
+                                 distribution_type=UNIFORM)
             DEALLOCATE (init_rng_seed)
          ELSE
-            CALL create_rng_stream(rng_stream=tmc_env%rng_stream, &
-                                   name="TMC_rng_stream", &
-                                   distribution_type=UNIFORM)
+            tmc_env%rng_stream = rng_stream_type( &
+                                 name="TMC_rng_stream", &
+                                 distribution_type=UNIFORM)
          END IF
 
          ! start running master and worker routines
@@ -310,8 +309,7 @@ CONTAINS
          END IF ! unused worker groups have nothing to do
 
          ! delete the random numbers
-         CPASSERT(ASSOCIATED(tmc_env%rng_stream))
-         CALL delete_rng_stream(tmc_env%rng_stream)
+         DEALLOCATE (tmc_env%rng_stream)
 
          ! deallocate the move types
          CALL finalize_mv_types(tmc_env%params)
diff --git a/src/tmc/tmc_tree_build.F b/src/tmc/tmc_tree_build.F
index 3450b945e6..68bbf7be05 100644
--- a/src/tmc/tmc_tree_build.F
+++ b/src/tmc/tmc_tree_build.F
@@ -33,10 +33,6 @@
 MODULE tmc_tree_build
    USE cp_log_handling,                 ONLY: cp_to_string
    USE kinds,                           ONLY: dp
-   USE parallel_rng_types,              ONLY: get_rng_stream,&
-                                              next_random_number,&
-                                              reset_to_next_rng_substream,&
-                                              set_rng_stream
    USE tmc_calculations,                ONLY: calc_e_kin,&
                                               init_vel
    USE tmc_dot_tree,                    ONLY: create_dot,&
@@ -290,18 +286,19 @@ CONTAINS
                                          nr_temp=tmc_env%params%nr_temp)
 
       ! use initial/default values
-      CALL get_rng_stream(rng_stream=tmc_env%rng_stream, &
-                          bg=tmc_env%m_env%gt_act%rng_seed(:, :, 1), &
-                          cg=tmc_env%m_env%gt_act%rng_seed(:, :, 2), &
-                          ig=tmc_env%m_env%gt_act%rng_seed(:, :, 3))
+      CALL tmc_env%rng_stream%get( &
+         bg=tmc_env%m_env%gt_act%rng_seed(:, :, 1), &
+         cg=tmc_env%m_env%gt_act%rng_seed(:, :, 2), &
+         ig=tmc_env%m_env%gt_act%rng_seed(:, :, 3))
 
       global_tree => tmc_env%m_env%gt_act
       tmc_env%m_env%gt_head => tmc_env%m_env%gt_act
 
       ! set global random seed
-      CALL set_rng_stream(rng_stream=tmc_env%rng_stream, bg=global_tree%rng_seed(:, :, 1), &
-                          cg=global_tree%rng_seed(:, :, 2), ig=global_tree%rng_seed(:, :, 3))
-      global_tree%rnd_nr = next_random_number(tmc_env%rng_stream)
+      CALL tmc_env%rng_stream%set(bg=global_tree%rng_seed(:, :, 1), &
+                                  cg=global_tree%rng_seed(:, :, 2), &
+                                  ig=global_tree%rng_seed(:, :, 3))
+      global_tree%rnd_nr = tmc_env%rng_stream%next()
 
       !-- SUBTREES: set initial values
       DO i = 1, SIZE(global_tree%conf)
@@ -328,10 +325,10 @@ CONTAINS
          END IF
 
          !-- different random seeds for every subtree
-         CALL reset_to_next_rng_substream(tmc_env%rng_stream)
-         CALL get_rng_stream(rng_stream=tmc_env%rng_stream, bg=global_tree%conf(i)%elem%rng_seed(:, :, 1), &
-                             cg=global_tree%conf(i)%elem%rng_seed(:, :, 2), &
-                             ig=global_tree%conf(i)%elem%rng_seed(:, :, 3))
+         CALL tmc_env%rng_stream%reset_to_next_substream()
+         CALL tmc_env%rng_stream%get(bg=global_tree%conf(i)%elem%rng_seed(:, :, 1), &
+                                     cg=global_tree%conf(i)%elem%rng_seed(:, :, 2), &
+                                     ig=global_tree%conf(i)%elem%rng_seed(:, :, 3))
 
          !-- gaussian distributed velocities
          !-- calculating the kinetic energy of the initial configuration velocity
@@ -645,16 +642,16 @@ CONTAINS
             new_elem%conf(:) = tmp_elem%conf(:)
 
             !-- set rnd nr generator and set next conf to change
-            CALL set_rng_stream(rng_stream=tmc_env%rng_stream, &
-                                bg=new_elem%parent%rng_seed(:, :, 1), &
-                                cg=new_elem%parent%rng_seed(:, :, 2), &
-                                ig=new_elem%parent%rng_seed(:, :, 3))
-            CALL reset_to_next_rng_substream(tmc_env%rng_stream)
+            CALL tmc_env%rng_stream%set( &
+               bg=new_elem%parent%rng_seed(:, :, 1), &
+               cg=new_elem%parent%rng_seed(:, :, 2), &
+               ig=new_elem%parent%rng_seed(:, :, 3))
+            CALL tmc_env%rng_stream%reset_to_next_substream()
             ! the random number for acceptance check
-            new_elem%rnd_nr = next_random_number(tmc_env%rng_stream)
+            new_elem%rnd_nr = tmc_env%rng_stream%next()
 
             ! the next configuration index to move
-            !rnd = next_random_number(tmc_env%rng_stream)
+            !rnd = tmc_env%rng_stream%next()
             !new_elem%mv_conf = 1+INT(size(new_elem%conf)*rnd)
             ! one temperature after each other
             new_elem%mv_conf = new_elem%parent%mv_next_conf
@@ -665,12 +662,11 @@ CONTAINS
             IF (n_acc) new_elem%Temp = tmp_elem%Temp*(1 - tmc_env%m_env%temp_decrease)
 
             !-- rnd for swap
-            rnd = next_random_number(tmc_env%rng_stream)
-            rnd2 = next_random_number(tmc_env%rng_stream)
-            CALL get_rng_stream(rng_stream=tmc_env%rng_stream, &
-                                bg=new_elem%rng_seed(:, :, 1), &
-                                cg=new_elem%rng_seed(:, :, 2), &
-                                ig=new_elem%rng_seed(:, :, 3))
+            rnd = tmc_env%rng_stream%next()
+            rnd2 = tmc_env%rng_stream%next()
+            CALL tmc_env%rng_stream%get(bg=new_elem%rng_seed(:, :, 1), &
+                                        cg=new_elem%rng_seed(:, :, 2), &
+                                        ig=new_elem%rng_seed(:, :, 3))
 
             ! swap moves are not part of the subtree structure,
             !  because exisiting elements from DIFFERENT subtrees are swaped
@@ -893,14 +889,14 @@ CONTAINS
       END IF
 
       ! set new substream of random number generator
-      CALL set_rng_stream(rng_stream=tmc_env%rng_stream, &
-                          bg=new_elem%rng_seed(:, :, 1), &
-                          cg=new_elem%rng_seed(:, :, 2), &
-                          ig=new_elem%rng_seed(:, :, 3))
-      CALL reset_to_next_rng_substream(tmc_env%rng_stream)
+      CALL tmc_env%rng_stream%set( &
+         bg=new_elem%rng_seed(:, :, 1), &
+         cg=new_elem%rng_seed(:, :, 2), &
+         ig=new_elem%rng_seed(:, :, 3))
+      CALL tmc_env%rng_stream%reset_to_next_substream()
 
       ! set the temperature for the NMC moves
-      rnd = next_random_number(tmc_env%rng_stream)
+      rnd = tmc_env%rng_stream%next()
       IF (tmc_env%params%NMC_inp_file .NE. "") THEN
          new_elem%temp_created = INT(tmc_env%params%nr_temp*rnd) + 1
       ELSE
@@ -908,15 +904,15 @@ CONTAINS
       END IF
 
       ! rnd nr for selecting move
-      rnd = next_random_number(tmc_env%rng_stream)
+      rnd = tmc_env%rng_stream%next()
       !-- set move type
       new_elem%move_type = select_random_move_type( &
                            move_types=tmc_env%params%move_types, &
                            rnd=rnd)
-      CALL get_rng_stream(rng_stream=tmc_env%rng_stream, &
-                          bg=new_elem%rng_seed(:, :, 1), &
-                          cg=new_elem%rng_seed(:, :, 2), &
-                          ig=new_elem%rng_seed(:, :, 3))
+      CALL tmc_env%rng_stream%get( &
+         bg=new_elem%rng_seed(:, :, 1), &
+         cg=new_elem%rng_seed(:, :, 2), &
+         ig=new_elem%rng_seed(:, :, 3))
 
       ! move is only done by the master,
       !  when standard MC moves with single potential are done
diff --git a/src/tmc/tmc_types.F b/src/tmc/tmc_types.F
index b7eb237c9f..16aed7c2be 100644
--- a/src/tmc/tmc_types.F
+++ b/src/tmc/tmc_types.F
@@ -60,7 +60,7 @@ MODULE tmc_types
    TYPE tmc_env_type
       TYPE(tmc_comp_set_type), POINTER              :: tmc_comp_set
       TYPE(tmc_param_type), POINTER                 :: params
-      TYPE(rng_stream_type), POINTER                :: rng_stream
+      TYPE(rng_stream_type), ALLOCATABLE            :: rng_stream
       TYPE(master_env_type), POINTER                :: m_env
       TYPE(worker_env_type), POINTER                :: w_env
    END TYPE tmc_env_type
@@ -199,7 +199,6 @@ CONTAINS
       NULLIFY (tmc_env%tmc_comp_set%para_env_m_ana)
       NULLIFY (tmc_env%tmc_comp_set%para_env_m_only)
 
-      NULLIFY (tmc_env%rng_stream)
       NULLIFY (tmc_env%m_env, tmc_env%w_env)
 
       ! initialize the parameter section
diff --git a/src/tmc/tmc_worker.F b/src/tmc/tmc_worker.F
index f38337d861..e3dd03b804 100644
--- a/src/tmc/tmc_worker.F
+++ b/src/tmc/tmc_worker.F
@@ -47,9 +47,6 @@ MODULE tmc_worker
    USE kinds,                           ONLY: default_string_length,&
                                               dp
    USE molecule_list_types,             ONLY: molecule_list_type
-   USE parallel_rng_types,              ONLY: get_rng_stream,&
-                                              next_random_number,&
-                                              set_rng_stream
    USE particle_list_types,             ONLY: particle_list_type
    USE tmc_analysis,                    ONLY: analysis_init,&
                                               analysis_restart_print,&
@@ -613,7 +610,7 @@ CONTAINS
       CPASSERT(ASSOCIATED(tmc_env))
       CPASSERT(ASSOCIATED(tmc_env%params))
       CPASSERT(ASSOCIATED(tmc_env%tmc_comp_set))
-      CPASSERT(ASSOCIATED(tmc_env%rng_stream))
+      CPASSERT(ALLOCATED(tmc_env%rng_stream))
       CPASSERT(ASSOCIATED(conf))
       CPASSERT(conf%temp_created .GT. 0)
       CPASSERT(conf%temp_created .LE. tmc_env%params%nr_temp)
@@ -662,15 +659,15 @@ CONTAINS
          END SELECT
 
          ! set move type
-         CALL set_rng_stream(rng_stream=tmc_env%rng_stream, &
-                             bg=conf%rng_seed(:, :, 1), cg=conf%rng_seed(:, :, 2), &
-                             ig=conf%rng_seed(:, :, 3))
+         CALL tmc_env%rng_stream%set( &
+            bg=conf%rng_seed(:, :, 1), cg=conf%rng_seed(:, :, 2), &
+            ig=conf%rng_seed(:, :, 3))
          conf%move_type = select_random_move_type( &
                           move_types=tmc_env%params%nmc_move_types, &
-                          rnd=next_random_number(tmc_env%rng_stream))
-         CALL get_rng_stream(rng_stream=tmc_env%rng_stream, &
-                             bg=conf%rng_seed(:, :, 1), cg=conf%rng_seed(:, :, 2), &
-                             ig=conf%rng_seed(:, :, 3))
+                          rnd=tmc_env%rng_stream%next())
+         CALL tmc_env%rng_stream%get( &
+            bg=conf%rng_seed(:, :, 1), cg=conf%rng_seed(:, :, 2), &
+            ig=conf%rng_seed(:, :, 3))
 
          ! do move
          CALL change_pos(tmc_params=tmc_env%params, &
@@ -704,13 +701,13 @@ CONTAINS
          END IF
 
          !check NMC step
-         CALL set_rng_stream(rng_stream=tmc_env%rng_stream, &
-                             bg=conf%rng_seed(:, :, 1), cg=conf%rng_seed(:, :, 2), &
-                             ig=conf%rng_seed(:, :, 3))
-         rnd_nr = next_random_number(tmc_env%rng_stream)
-         CALL get_rng_stream(rng_stream=tmc_env%rng_stream, &
-                             bg=conf%rng_seed(:, :, 1), cg=conf%rng_seed(:, :, 2), &
-                             ig=conf%rng_seed(:, :, 3))
+         CALL tmc_env%rng_stream%set( &
+            bg=conf%rng_seed(:, :, 1), cg=conf%rng_seed(:, :, 2), &
+            ig=conf%rng_seed(:, :, 3))
+         rnd_nr = tmc_env%rng_stream%next()
+         CALL tmc_env%rng_stream%get( &
+            bg=conf%rng_seed(:, :, 1), cg=conf%rng_seed(:, :, 2), &
+            ig=conf%rng_seed(:, :, 3))
 
          IF (.NOT. change_rejected) THEN
             CALL acceptance_check(tree_element=conf, parent_element=last_acc_conf, &
diff --git a/tests/LIBTEST/test_01.inp b/tests/LIBTEST/test_01.inp
index 7a86404617..be4afba40b 100644
--- a/tests/LIBTEST/test_01.inp
+++ b/tests/LIBTEST/test_01.inp
@@ -12,5 +12,4 @@
   FFT 10
   CLEBSCH_GORDON 10
   MPI 4
-  RANDOM_NUMBER_GENERATOR 1000
 &END