comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 14:15:19 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

adapt to OO-style RNG type and drop unnecessary use of pointers

Browse source 
Fortran does pass-by-reference, but `POINTER` was needed for the
procedures acting on the `rng_stream_type`. With those refactored as
type-bound procedures, almost all uses of pointers can be avoided.

One of the two remaining cases (in `helium_methods.F`) could likely be
avoided as well by using an allocatable and `move_alloc`, but this would
need another nested type like `rng_stream_p_type` but with an
allocatable member instead of a pointer (or a complete redesign of the
RNG setup in that part).
This commit is contained in:
Tiziano Müller 2019-07-18 15:47:06 +02:00 • committed by Tiziano Müller
parent 9579ec3659
commit b8d84a3028
56 changed files with 553 additions and 972 deletions
Download patch file Download diff file Expand all files Collapse all files

7
src/common/distribution_1d_types.F
View file

@ -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)

20
src/csvr_system_types.F
View file

@ -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

11
src/csvr_system_utils.F
View file

@ -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

21
src/distribution_methods.F
View file

@ -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

23
src/environment.F
View file

@ -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")

20
src/ewalds_multipole.F
View file

@ -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))

19
src/fm/cp_fm_types.F
View file

@ -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)

13
src/gle_system_types.F
View file

@ -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

11
src/global_types.F
View file

@ -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

24
src/hfx_load_balance_methods.F
View file

@ -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)

8
src/input_restart_force_eval.F
View file

@ -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

135
src/library_tests.F
View file

@ -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

3
src/metadynamics.F
View file

@ -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)

20
src/metadynamics_types.F
View file

@ -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)

3
src/mode_selective.F
View file

@ -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)))

3
src/motion/dimer_types.F
View file

@ -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

5
src/motion/helium_common.F
View file

@ -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

100
src/motion/helium_methods.F
View file

@ -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

47
src/motion/helium_sampling.F
View file

@ -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

173
src/motion/helium_worm.F
View file

@ -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

24
src/motion/input_cp2k_restarts.F
View file

@ -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

28
src/motion/integrator.F
View file

@ -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

8
src/motion/mc/mc_control.F
View file

@ -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

33
src/motion/mc/mc_coordinates.F
View file

@ -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

33
src/motion/mc/mc_ensembles.F
View file

@ -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

25
src/motion/mc/mc_ge_moves.F
View file

@ -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

70
src/motion/mc/mc_moves.F
View file

@ -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

17
src/motion/mc/mc_run.F
View file

@ -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
! **************************************************************************************************

40
src/motion/mc/tamc_run.F
View file

@ -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)

40
src/motion/md_vel_utils.F
View file

@ -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

8
src/motion/neb_md_utils.F
View file

@ -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)

13
src/motion/pint_gle.F
View file

@ -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

46
src/motion/pint_methods.F
View file

@ -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

25
src/motion/pint_piglet.F
View file

@ -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)

24
src/motion/pint_pile.F
View file

@ -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

13
src/motion/pint_public.F
View file

@ -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) + &

42
src/motion/pint_qtb.F
View file

@ -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

8
src/motion/pint_types.F
View file

@ -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

10
src/motion/thermostat/al_system_dynamics.F
View file

@ -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

7
src/motion/thermostat/csvr_system_init.F
View file

@ -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__, &

7
src/motion/thermostat/extended_system_init.F
View file

@ -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

20
src/motion/thermostat/gle_system_dynamics.F
View file

@ -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

28
src/motion/wiener_process.F
View file

@ -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)

11
src/optimize_input.F
View file

@ -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

14
src/pao_ml_neuralnet.F
View file

@ -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))

8
src/start/input_cp2k.F
View file

@ -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. "// &

6
src/statistical_methods.F
View file

@ -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

14
src/swarm/glbopt_mincrawl.F
View file

@ -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

4
src/swarm/glbopt_worker.F
View file

@ -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

19
src/tmc/tmc_calculations.F
View file

@ -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

83
src/tmc/tmc_moves.F
View file

@ -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

20
src/tmc/tmc_setup.F
View file

@ -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)

74
src/tmc/tmc_tree_build.F
View file

@ -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

3
src/tmc/tmc_types.F
View file

@ -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

33
src/tmc/tmc_worker.F
View file

@ -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, &

1
tests/LIBTEST/test_01.inp
View file

@ -12,5 +12,4 @@
FFT 10
CLEBSCH_GORDON 10
MPI 4
RANDOM_NUMBER_GENERATOR 1000
&END