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