Search code examples
randomfortrangfortranrandom-seed

Exact way of defining seed for random number generation in fortran


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.

I just want to fix a particular seed for my problem, start generating random numbers, and concentrate on the other physical aspects!

The questions are -

  • How exactly to:
    (a) call random_seed(),
    (b) defining the exact integral value of seed (like e.g. "123456789") for initialisation,
    (c) putting call random_number(),
    (d) what argument to put in random_seed(),
    (e) SIZE, GET, PUT ... what to use?

Solution

  • 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