1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
|
!
! Copyright (C) 2001-2012 Quantum ESPRESSO 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 .
!
!----------------------------------------------------------------------------
MODULE random_numbers
!----------------------------------------------------------------------------
!
USE kinds, ONLY : DP
!
IMPLICIT NONE
!
INTERFACE gauss_dist
!
MODULE PROCEDURE gauss_dist_scal, gauss_dist_vect
!
END INTERFACE
!
CONTAINS
!
!------------------------------------------------------------------------
FUNCTION randy ( irand )
!------------------------------------------------------------------------
!
! x=randy(n): reseed with initial seed idum=n ( 0 <= n <= ic, see below)
! if randy is not explicitly initialized, it will be
! initialized with seed idum=0 the first time it is called
! x=randy() : generate uniform real(DP) numbers x in [0,1]
!
REAL(DP) :: randy
INTEGER, optional :: irand
!
INTEGER , PARAMETER :: m = 714025, &
ia = 1366, &
ic = 150889, &
ntab = 97
REAL(DP), PARAMETER :: rm = 1.0_DP / m
INTEGER :: j
INTEGER, SAVE :: ir(ntab), iy, idum=0
LOGICAL, SAVE :: first=.true.
!
IF ( present(irand) ) THEN
idum = MIN( ABS(irand), ic)
first=.true.
END IF
IF ( first ) THEN
!
first = .false.
idum = MOD( ic - idum, m )
!
DO j=1,ntab
idum=mod(ia*idum+ic,m)
ir(j)=idum
END DO
idum=mod(ia*idum+ic,m)
iy=idum
END IF
j=1+(ntab*iy)/m
IF( j > ntab .OR. j < 1 ) call errore('randy','j out of range',ABS(j)+1)
iy=ir(j)
randy=iy*rm
idum=mod(ia*idum+ic,m)
ir(j)=idum
!
RETURN
!
END FUNCTION randy
!
!------------------------------------------------------------------------
SUBROUTINE set_random_seed ( )
!------------------------------------------------------------------------
!
! poor-man random seed for randy
!
INTEGER, DIMENSION (8) :: itime
INTEGER :: iseed, irand
!
CALL date_and_time ( values = itime )
! itime contains: year, month, day, time difference in minutes, hours,
! minutes, seconds and milliseconds.
iseed = ( itime(8) + itime(6) ) * ( itime(7) + itime(4) )
irand = randy ( iseed )
!
END SUBROUTINE set_random_seed
!
!-----------------------------------------------------------------------
FUNCTION gauss_dist_scal( mu, sigma )
!-----------------------------------------------------------------------
!
! ... this function generates a number taken from a normal
! ... distribution of mean value \mu and variance \sigma
!
IMPLICIT NONE
!
REAL(DP), INTENT(IN) :: mu
REAL(DP), INTENT(IN) :: sigma
REAL(DP) :: gauss_dist_scal
!
REAL(DP) :: x1, x2, w
!
!
gaussian_loop: DO
!
x1 = 2.0_DP * randy() - 1.0_DP
x2 = 2.0_DP * randy() - 1.0_DP
!
w = x1 * x1 + x2 * x2
!
IF ( w < 1.0_DP ) EXIT gaussian_loop
!
END DO gaussian_loop
!
w = SQRT( ( - 2.0_DP * LOG( w ) ) / w )
!
gauss_dist_scal = x1 * w * sigma + mu
!
RETURN
!
END FUNCTION gauss_dist_scal
!
!-----------------------------------------------------------------------
FUNCTION gauss_dist_vect( mu, sigma, dim )
!-----------------------------------------------------------------------
!
! ... this function generates an array of numbers taken from a normal
! ... distribution of mean value \mu and variance \sigma
!
IMPLICIT NONE
!
REAL(DP), INTENT(IN) :: mu
REAL(DP), INTENT(IN) :: sigma
INTEGER, INTENT(IN) :: dim
REAL(DP) :: gauss_dist_vect( dim )
!
REAL(DP) :: x1, x2, w
INTEGER :: i
!
!
DO i = 1, dim, 2
!
gaussian_loop: DO
!
x1 = 2.0_DP * randy() - 1.0_DP
x2 = 2.0_DP * randy() - 1.0_DP
!
w = x1 * x1 + x2 * x2
!
IF ( w < 1.0_DP ) EXIT gaussian_loop
!
END DO gaussian_loop
!
w = SQRT( ( - 2.0_DP * LOG( w ) ) / w )
!
gauss_dist_vect(i) = x1 * w * sigma
!
IF ( i >= dim ) EXIT
!
gauss_dist_vect(i+1) = x2 * w * sigma
!
END DO
!
gauss_dist_vect(:) = gauss_dist_vect(:) + mu
!
RETURN
!
END FUNCTION gauss_dist_vect
!
END MODULE random_numbers
|
