1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     #############################################################
9c     ##                                                         ##
10c     ##  function random  --  portable random number generator  ##
11c     ##                                                         ##
12c     #############################################################
13c
14c
15c     "random" generates a random number on [0,1] via a long
16c     period generator due to L'Ecuyer with Bays-Durham shuffle
17c
18c     literature references:
19c
20c     P. L'Ecuyer, Communications of the ACM, 31, 742-774 (1988)
21c
22c     W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P.
23c     Flannery, Numerical Recipes (Fortran), 2nd Ed., Cambridge
24c     University Press, 1992, Section 7.1
25c
26c
27      function random ()
28      use inform
29      use iounit
30      use keys
31      implicit none
32      integer im1,ia1,iq1,ir1
33      integer im2,ia2,iq2,ir2
34      integer big,nshuffle
35      integer imm1,ndiv
36      real*8 factor
37      parameter (im1=2147483563)
38      parameter (ia1=40014)
39      parameter (iq1=53668)
40      parameter (ir1=12211)
41      parameter (im2=2147483399)
42      parameter (ia2=40692)
43      parameter (iq2=52774)
44      parameter (ir2=3791)
45      parameter (big=141803398)
46      parameter (nshuffle=32)
47      parameter (imm1=im1-1)
48      parameter (ndiv=1+imm1/nshuffle)
49      parameter (factor=1.0d0/im1)
50      integer i,k,iy,next
51      integer seed,seed2
52      integer year,month,day
53      integer hour,minute,second
54      integer ishuffle(nshuffle)
55      real*8 random
56      logical first
57      character*20 keyword
58      character*240 record
59      character*240 string
60      save first
61      save seed,seed2
62      save iy,ishuffle
63      data first  / .true. /
64c
65c
66c     random number seed is first set to a big number,
67c     then incremented by the seconds elapsed this decade
68c
69      if (first) then
70         first = .false.
71         seed = big
72         call calendar (year,month,day,hour,minute,second)
73         year = mod(year,10)
74         seed = seed + 32140800*year + 2678400*(month-1)
75         seed = seed + 86400*(day-1) + 3600*hour
76         seed = seed + 60*minute + second
77c
78c     search the keywords for a random number seed
79c
80         do i = 1, nkey
81            next = 1
82            record = keyline(i)
83            call gettext (record,keyword,next)
84            call upcase (keyword)
85            if (keyword(1:11) .eq. 'RANDOMSEED ') then
86               string = record(next:240)
87               read (string,*,err=10)  seed
88               seed = max(1,seed)
89            end if
90   10       continue
91         end do
92c
93c     print the value used for the random number seed
94c
95         if (verbose) then
96            write (iout,20)  seed
97   20       format (/,' Random Number Generator Initialized',
98     &                 ' with SEED :',3x,i12)
99         end if
100c
101c     warm up and then load the shuffling table
102c
103         seed2 = seed
104         do i = nshuffle+8, 1, -1
105            k = seed / iq1
106            seed = ia1 * (seed-k*iq1) - k*ir1
107            if (seed .lt. 0)  seed = seed + im1
108            if (i .le. nshuffle)  ishuffle(i) = seed
109         end do
110         iy = ishuffle(1)
111      end if
112c
113c     get a new random number value each call
114c
115      k = seed / iq1
116      seed = ia1*(seed-k*iq1) - k*ir1
117      if (seed .lt. 0)  seed = seed + im1
118      k = seed2 / iq2
119      seed2 = ia2*(seed2-k*iq2) - k*ir2
120      if (seed2 .lt. 0)  seed2 = seed2 + im2
121      i = 1 + iy/ndiv
122      iy = ishuffle(i) - seed2
123      ishuffle(i) = seed
124      if (iy .lt. 1)  iy = iy + imm1
125      random = factor * iy
126c
127c     print the value of the current random number
128c
129c     if (debug) then
130c        write (iout,30)  random
131c  30    format (' RANDOM  --  The Random Number Value is',f12.8)
132c     end if
133      return
134      end
135c
136c
137c     ############################################################
138c     ##                                                        ##
139c     ##  function normal  --  random number from normal curve  ##
140c     ##                                                        ##
141c     ############################################################
142c
143c
144c     "normal" generates a random number from a normal Gaussian
145c     distribution with a mean of zero and a variance of one
146c
147c
148      function normal ()
149      use inform
150      use iounit
151      implicit none
152      real*8 v1,v2,rsq
153      real*8 factor,store
154      real*8 normal,random
155      logical compute
156      external random
157      save compute,store
158      data compute  / .true. /
159c
160c
161c     get a pair of random values from the distribution
162c
163      if (compute) then
164   10    continue
165         v1 = 2.0d0*random() - 1.0d0
166         v2 = 2.0d0*random() - 1.0d0
167         rsq = v1**2 + v2**2
168         if (rsq .ge. 1.0d0)  goto 10
169         factor = sqrt(-2.0d0*log(rsq)/rsq)
170         store = v1 * factor
171         normal = v2 * factor
172         compute = .false.
173c
174c     use the second random value computed at the last call
175c
176      else
177         normal = store
178         compute = .true.
179      end if
180c
181c     print the value of the current random number
182c
183c     if (debug) then
184c        write (iout,20)  normal
185c  20    format (' NORMAL  --  The Random Number Value is',f12.8)
186c     end if
187      return
188      end
189c
190c
191c     ##############################################################
192c     ##                                                          ##
193c     ##  subroutine ranvec  --  unit vector in random direction  ##
194c     ##                                                          ##
195c     ##############################################################
196c
197c
198c     "ranvec" generates a unit vector in 3-dimensional
199c     space with uniformly distributed random orientation
200c
201c     literature references:
202c
203c     G. Marsaglia, Ann. Math. Stat., 43, 645 (1972)
204c
205c     R. C. Rapaport, The Art of Molecular Dynamics Simulation,
206c     2nd Edition, Cambridge University Press, 2004, Section 18.4
207c
208c
209      subroutine ranvec (vector)
210      use inform
211      use iounit
212      implicit none
213      real*8 x,y,s
214      real*8 random
215      real*8 vector(3)
216      external random
217c
218c
219c     get a pair of appropriate components in the plane
220c
221      s = 2.0d0
222      do while (s .ge. 1.0d0)
223         x = 2.0d0*random() - 1.0d0
224         y = 2.0d0*random() - 1.0d0
225         s = x**2 + y**2
226      end do
227c
228c     construct the 3-dimensional random unit vector
229c
230      vector(3) = 1.0d0 - 2.0d0*s
231      s = 2.0d0 * sqrt(1.0d0 - s)
232      vector(2) = s * y
233      vector(1) = s * x
234c
235c     print the components of the random unit vector
236c
237c     if (debug) then
238c        write (iout,10)  vector(1),vector(2),vector(3)
239c  10    format (' RANVEC  --  The Random Vector is',3f10.4)
240c     end if
241      return
242      end
243c
244c
245c     ##############################################################
246c     ##                                                          ##
247c     ##  subroutine sphere  --  uniform set of points on sphere  ##
248c     ##                                                          ##
249c     ##############################################################
250c
251c
252c     "sphere" finds a specified number of uniformly distributed
253c     points on a sphere of unit radius centered at the origin
254c
255c     literature reference:
256c
257c     E. B. Saff and A. B. J. Kuijlaars, "Distributing Many
258c     Points on a Sphere", The Mathematical Intelligencer,
259c     19, 5-11 (1997)
260c
261c
262      subroutine sphere (ndot,dot)
263      use math
264      implicit none
265      integer i,ndot
266      real*8 theta,phi
267      real*8 h,phiold
268      real*8 tot,tot1
269      real*8 dot(3,*)
270c
271c
272c     find spherical coordinates then convert to Cartesian
273c
274      tot = dble(ndot)
275      tot1 = dble(ndot-1)
276      do i = 1, ndot
277         h = -1.0d0 + 2.0d0*dble(i-1)/tot1
278         h = min(1.0d0,h)
279         theta = acos(h)
280         if (i.eq.1 .or. i.eq.ndot) then
281            phi = 0.0d0
282         else
283            phi = mod(phiold+3.6d0/sqrt(tot*(1.0d0-h*h)),2.0d0*pi)
284         end if
285         dot(1,i) = sin(theta) * cos(phi)
286         dot(2,i) = sin(theta) * sin(phi)
287         dot(3,i) = cos(theta)
288         phiold = phi
289      end do
290      return
291      end
292