!
! Copyright (C) 2001 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
! USE:
! 
! call once subroutine set_rndm_seed(iseed) with an integer
! value for iseed in order to initialize the Random Number Generator
!
! subsequent calls to the real function rndm() will return 
! random numbers uniformely distributed in the [0,1) interval
!
!
function rndm ()
  !
  !   RANDOM NUMBER GENERATOR equivalent to ran1 of Num.Rec.
  !
  implicit none
  integer :: irand
  common/random_number/ irand
  real(kind(0.0d0)) :: rndm, shuffle (32)
  real(kind(0.0d0)),external :: rndx
  integer :: i
  logical :: first
  data first / .true. /
  save first, shuffle, i
  !
  if (first.or.irand.lt.0) then
     irand = - irand
     if (first) irand=1 ! starting seed, must be not be 0
     do i = 32 + 8, 1, - 1
        shuffle (min (i, 32) ) = rndx (irand)
     enddo
     i = 32 * shuffle (1) + 1
     first = .false.
  endif
  rndm = shuffle (i)
  shuffle (i) = rndx (irand)

  i = 32 * rndm + 1
  return

end function rndm

function rndx (irand)
  !
  !   RANDOM NUMBER GENERATOR equivalent to ran0 of Num.Rec.
  !
  implicit none
  integer :: im, ia, iq, ir, irand, is, it
  real(kind(1.d0)) :: rndx, obm
  logical :: first
  data first / .true. /
  save im, ia, iq, ir, obm, first
  !
  if (first) then
     ! this is 2**31-1 avoiding overflow
     im = 2 * (2**30 - 1) + 1
     obm = 1.0 / im
     ia = 7**5
     iq = im / ia
     ir = im - ia * iq
     ! starting seed, must be not be 0
!     irand = 1
     first = .false.
  endif
  is = irand / iq
  it = irand-is * iq
  irand = ia * it - is * ir
  if (irand.lt.0) irand = irand+im 
  rndx = irand * obm
  return
end function rndx


subroutine set_rndm_seed(iseed)
!
! this subroutine initialize the random number with the given seed
!
  implicit none
  integer :: irand,iseed
  common/random_number/ irand
  real(kind(1.d0)) :: dummy, rndm
  if (iseed.le.0) stop 'set_rndm_seed: seed should be a positive integer !'
! make sure rndm() has been called already once !
  dummy = rndm()
  irand = - iseed

  return
end subroutine set_rndm_seed

