I am trying to simulate a physical problem (Brownian dynamics) for which I need to generate random numbers.
Previously, I was just using the following to choose random numbers from a Gaussian distribution:
subroutine box_muller(l)
implicit none
real :: l, u1, u2
real, parameter :: pi = acos(-1.0)
call random_seed()
call random_number(u1)
call random_number(u2)
l = sqrt(-2*log(u1))*cos(2*pi*u2)
end subroutine
which is the Box-Muller algorithm.
I was calling this subroutine at each instant (time-step).
But I think I have not yet understood how to properly use the call random_seed()
stuff.
For example, I was trying the following after some internet searches:
program seeding
implicit none
integer:: n
INTEGER :: seed
REAL :: num
CALL RANDOM_SEED (size=seed)
seed = 123456789
DO n = 1, 5
CALL RANDOM_NUMBER(num)
PRINT *, num
END DO
end program
It gives the following output:
4.80166674E-02
0.797537327
0.294925630
0.107502699
0.987043023
Though I don't know if I have did it correctly.
I also tried the following in my main problem:
subroutine box_muller(x)
implicit none
integer :: n
integer,allocatable :: seed(:)
real,intent(out) :: x
real,parameter :: pi=3.14159265
real :: u1,u2
call random_seed(size=n)
allocate(seed(n))
seed = 12 ! putting arbitrary seed to all elements
call random_number(u1)
call random_number(u2)
x = sqrt(-2*log(u1))*cos(2*pi*u2)
end subroutine
But BELIEVE ME, I UNDERSTOOD NOTHING OF WHAT I AM DOING.
The questions are -
call random_seed()
,call random_number()
,Do you want/need the sequence of random numbers to be repeatable from run to run? If no, you don't have to call random_seed()
at all. If yes, your last piece of code is almost correct, but it misses a call:
integer :: n
integer,allocatable :: seed(:)
call random_seed(size=n) ! n is processor-dependent
allocate(seed(n))
seed = 12 ! putting arbitrary seed to all elements
call random_seed(put=seed) ! effectively initializing the RNG with the seed