2 !-----------------------------------------------------------------------------
4 use prng ! prng.f90 or prng_32.f90
8 !-----------------------------------------------------------------------------
10 !-----------------------------------------------------------------------------
12 !-----------------------------------------------------------------------------
13 real(kind=8) function ran_number(x1,x2)
14 ! Calculate a random real number from the range (x1,x2).
16 use MPI_data ! include 'COMMON.SETUP'
20 ! implicit real*8 (a-h,o-z)
21 ! include 'DIMENSIONS'
22 real(kind=8) :: x1,x2,fctor
24 integer :: ix(1) !el ix ---> ix(1)
25 data fctor /2147483647.0D0/
27 ran_number=x1+(x2-x1)*prng_next(me)
30 ran_number=x1+(x2-x1)*ix(1)/fctor
33 end function ran_number
34 !-----------------------------------------------------------------------------
35 integer function iran_num(n1,n2)
36 ! Calculate a random integer number from the range (n1,n2).
38 use MPI_data ! include 'COMMON.SETUP'
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 integer :: n1,n2,ix(1) !el ix ---> ix(1)
44 real(kind=4) :: fctor=2147483647.0
46 ix(1)=n1+(n2-n1+1)*prng_next(me)
47 if (ix(1).lt.n1) ix(1)=n1
48 if (ix(1).gt.n2) ix(1)=n2
52 ix(1)=n1+(n2-n1+1)*(ix(1)/fctor)
53 if (ix(1).gt.n2) ix(1)=n2
58 !-----------------------------------------------------------------------------
59 real(kind=8) function binorm(x1,x2,sigma1,sigma2,ak)
60 ! implicit real*8 (a-h,o-z)
62 real(kind=8) :: x1,x2,sigma1,sigma2,ak,alowb,aupb,seg,alen
63 ! print '(a)','Enter BINORM.'
64 alowb=dmin1(x1-3.0D0*sigma1,x2-3.0D0*sigma2)
65 aupb=dmax1(x1+3.0D0*sigma1,x2+3.0D0*sigma2)
66 seg=sigma1/(sigma1+ak*sigma2)
67 alen=ran_number(0.0D0,1.0D0)
69 binorm=anorm_distr(x1,sigma1,alowb,aupb)
71 binorm=anorm_distr(x2,sigma2,alowb,aupb)
73 ! print '(a)','Exiting BINORM.'
76 !-----------------------------------------------------------------------------
77 real(kind=8) function anorm_distr(x,sigma,alowb,aupb)
78 ! implicit real*8 (a-h,o-z)
79 ! to make a normally distributed deviate with zero mean and unit variance
82 real(kind=8) :: fac,gset,rsq,v1,v2,ran1
83 real(kind=8) :: x,sigma,alowb,aupb,gaussdev
86 !elwrite(iout,*) "anorm distr start",x,sigma,alowb,aupb
88 1 v1=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
89 v2=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
91 if(rsq.ge.1.d0.or.rsq.eq.0.0d0) goto 1
92 fac=sqrt(-2.0d0*log(rsq)/rsq)
100 anorm_distr=x+gaussdev*sigma
101 !elwrite(iout,*) "anorm distr end",x,sigma,alowb,aupb,anorm_distr
103 end function anorm_distr
104 !-----------------------------------------------------------------------------
105 subroutine mult_norm(lda,n,a,x,fail)
107 ! Generate the vector X whose elements obey the multiple-normal distribution
108 ! from exp(-0.5*X'AX). LDA is the leading dimension of the moment matrix A,
109 ! n is the dimension of the problem. FAIL is set at .TRUE., if the smallest
110 ! eigenvalue of the matrix A is close to 0.
112 ! implicit double precision (a-h,o-z)
114 real(kind=8),dimension(lda,n) :: a
115 real(kind=8),dimension(n) :: x
116 real(kind=8),dimension(100) :: eig,work
117 real(kind=8),dimension(3,3) :: vec
118 real(kind=8) :: eig_limit=1.0D-8
120 real(kind=8) :: sigma,alim,xi
122 ! print '(a)','Enter MULT_NORM.'
124 ! Find the smallest eigenvalue of the matrix A.
127 ! print '(8f10.5)',(a(i,j),j=1,n)
130 call f02faf('V','U',2,a,lda,eig,work,100,ifail)
132 call djacob(2,lda,10000,1.0d-10,a,vec,eig)
134 ! print '(8f10.5)',(eig(i),i=1,n)
137 ! print '(8f10.5)',(a(i,j),j=1,n)
139 if (eig(1).lt.eig_limit) then
140 print *,'From Mult_Norm: Eigenvalues of A are too small.'
145 ! Generate points following the normal distributions along the principal
146 ! axes of the moment matrix. Store in WORK.
149 sigma=1.0D0/dsqrt(eig(i))
151 work(i)=anorm_distr(0.0D0,sigma,-alim,alim)
154 ! Transform the vector of normal variables back to the original basis.
164 end subroutine mult_norm
165 !-----------------------------------------------------------------------------
166 subroutine mult_norm1(lda,n,a,z,box,x,fail)
168 ! Generate the vector X whose elements obey the multi-gaussian multi-dimensional
169 ! distribution from sum_{i=1}^m W(i)exp[-0.5*X'(i)A(i)X(i)]. LDA is the
170 ! leading dimension of the moment matrix A, n is the dimension of the
171 ! distribution, nlob is the number of lobes. FAIL is set at .TRUE., if the
172 ! smallest eigenvalue of the matrix A is close to 0.
174 ! implicit real*8 (a-h,o-z)
175 ! include 'DIMENSIONS'
177 use MPI_data !include 'COMMON.SETUP'
182 integer :: lda,n,i,j,istep
183 real(kind=8),dimension(lda,n) :: a
184 real(kind=8),dimension(n) :: z,x
185 real(kind=8),dimension(n,n) :: box
186 real(kind=8) :: etmp,ww,xi,dec
187 ! include 'COMMON.IOUNITS'
190 ! Generate points following the normal distributions along the principal
191 ! axes of the moment matrix. Store in WORK.
193 !d print *,'CG Processor',me,' entered MultNorm1.'
194 !d print '(2(1pe14.4),3x,1pe14.4)',((a(i,j),j=1,2),z(i),i=1,2)
196 !d print *,i,box(1,i),box(2,i)
200 if (istep.gt.10000) then
201 ! write (iout,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
203 ! write (*,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
205 ! write (iout,*) 'box',box
206 ! write (iout,*) 'a',a
207 ! write (iout,*) 'z',z
212 x(i)=ran_number(box(1,i),box(2,i))
217 ww=ww+0.5D0*a(i,i)*xi*xi
219 ww=ww+a(i,j)*xi*pinorm(x(j)-z(j))
222 dec=ran_number(0.0D0,1.0D0)
223 ! print *,(x(i),i=1,n),ww,dexp(-ww),dec
224 !rc if (dec.gt.dexp(-ww)) goto 10
230 if (dec.gt.etmp) goto 10
231 !d print *,'CG Processor',me,' exitting MultNorm1.'
233 end subroutine mult_norm1
234 !-----------------------------------------------------------------------------
236 !-----------------------------------------------------------------------------
239 integer,parameter :: IA=16807,IM=2147483647,IQ=127773,IR=2836
240 integer,parameter :: MASK=123459876
242 real(kind=4),parameter :: AM=1./IM
246 idum=IA*(idum-k*IQ)-IR*k
247 if (idum.lt.0) idum=idum+IM
252 ! (C) Copr. 1986-92 Numerical Recipes Software *11915
253 !-----------------------------------------------------------------------------
256 integer,parameter :: IA=16807,IM=2147483647,IQ=127773,IR=2836
257 integer,parameter :: NTAB=32,NDIV=1+(IM-1)/NTAB
259 real(kind=4),parameter :: AM=1./IM,EPS=1.2e-7,RNMX=1.-EPS
261 integer,dimension(NTAB) :: iv
263 DATA iv /NTAB*0/, iy /0/
264 if (idum.le.0.or.iy.eq.0) then
268 idum=IA*(idum-k*IQ)-IR*k
269 if (idum.lt.0) idum=idum+IM
270 if (j.le.NTAB) iv(j)=idum
275 idum=IA*(idum-k*IQ)-IR*k
276 if (idum.lt.0) idum=idum+IM
283 ! (C) Copr. 1986-92 Numerical Recipes Software *11915
284 !-----------------------------------------------------------------------------
287 integer,parameter :: IM1=2147483563,IM2=2147483399,IMM1=IM1-1,&
288 IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,&
289 NTAB=32,NDIV=1+IMM1/NTAB
291 real(kind=4),parameter :: AM=1./IM1,EPS=1.2e-7,RNMX=1.-EPS
292 integer :: idum2,j,k,iy
293 integer,dimension(NTAB) :: iv
295 DATA idum2/123456789/, iv/NTAB*0/, iy/0/
301 idum=IA1*(idum-k*IQ1)-k*IR1
302 if (idum.lt.0) idum=idum+IM1
303 if (j.le.NTAB) iv(j)=idum
308 idum=IA1*(idum-k*IQ1)-k*IR1
309 if (idum.lt.0) idum=idum+IM1
311 idum2=IA2*(idum2-k*IQ2)-k*IR2
312 if (idum2.lt.0) idum2=idum2+IM2
316 if(iy.lt.1)iy=iy+IMM1
320 ! (C) Copr. 1986-92 Numerical Recipes Software *11915
321 !-----------------------------------------------------------------------------
324 integer,parameter :: MBIG=1000000000,MSEED=161803398,MZ=0
327 real(kind=4),parameter :: FAC=1./MBIG
328 ! PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG)
329 integer :: i,iff,ii,inext,inextp,k
331 integer,dimension(55) :: ma
333 SAVE iff,inext,inextp,ma
335 if(idum.lt.0.or.iff.eq.0)then
345 if(mk.lt.MZ)mk=mk+MBIG
350 ma(i)=ma(i)-ma(1+mod(i+30,55))
351 if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG
359 if(inext.eq.56)inext=1
361 if(inextp.eq.56)inextp=1
362 mj=ma(inext)-ma(inextp)
363 if(mj.lt.MZ)mj=mj+MBIG
368 ! (C) Copr. 1986-92 Numerical Recipes Software *11915
369 !-----------------------------------------------------------------------------
371 !-----------------------------------------------------------------------------
372 ! $Date: 1994/10/04 16:19:52 $
376 ! See help for RANDOMV on the PSFSHARE disk to understand these
377 ! subroutines. This is the VS Fortran version of this code.
380 subroutine VRND(VEC,N)
383 !el integer,dimension(250) :: A
384 integer :: LOOP,N !el,I,I147
385 integer,dimension(N) :: VEC
386 !el COMMON /VRANDD/ A, I, I147
389 IF(.NOT.(I.GE.251))GOTO 23002
393 IF(.NOT.(I147.GE.251))GOTO 23004
396 A(I)=IEOR(A(I147),A(I))
401 !-----------------------------------------------------------------------------
402 real(kind=8) function RNDV(IDUM)
404 real(kind=8) :: RM1,RM2,R(99)
405 INTEGER :: IA1,IC1,M1, IA2,IC2,M2, IA3,IC3,M3, IDUM
406 INTEGER :: IX1,IX2,IX3,J
408 DATA IA1,IC1,M1/1279,351762,1664557/
409 DATA IA2,IC2,M2/2011,221592,1048583/
410 DATA IA3,IC3,M3/15551,6150,29101/
411 IF(.NOT.(IDUM.LT.0))GOTO 23006
413 IX1 = MOD(IA1*IX1+IC1,M1)
415 IX1 = MOD(IA1*IX1+IC1,M1)
420 IX1 = MOD(IA1*IX1+IC1,M1)
421 IX2 = MOD(IA2*IX2+IC2,M2)
422 R(J) = (DBLE(IX1)+DBLE(IX2)*RM2)*RM1
425 IX1 = MOD(IA1*IX1+IC1,M1)
426 IX2 = MOD(IA2*IX2+IC2,M2)
427 IX3 = MOD(IA3*IX3+IC3,M3)
430 R(J) = (DBLE(IX1)+DBLE(IX2)*RM2)*RM1
434 !-----------------------------------------------------------------------------
435 subroutine VRNDST(SEED)
438 !el INTEGER :: A(250)
439 INTEGER :: LOOP,IDUM !el,I,I147
441 !el real(kind=8) :: RNDV
442 !el COMMON /VRANDD/ A, I, I147
447 A(LOOP)=INT(RNDV(IDUM)*2147483647)
450 end subroutine VRNDST
451 !-----------------------------------------------------------------------------
452 subroutine VRNDIN(IODEV)
454 INTEGER :: IODEV !el, A(250),I, I147
455 !el COMMON/VRANDD/ A, I, I147
456 READ(IODEV) A, I, I147
458 end subroutine VRNDIN
459 !-----------------------------------------------------------------------------
460 subroutine VRNDOU(IODEV)
461 ! This corresponds to VRNDOUT in the APFTN64 version
463 INTEGER :: IODEV !el, A(250),I, I147
464 !el COMMON/VRANDD/ A, I, I147
465 WRITE(IODEV) A, I, I147
467 end subroutine VRNDOU
468 !-----------------------------------------------------------------------------
469 real(kind=8) function RNUNF(N)
470 INTEGER :: IRAN1(2000),N
471 real(kind=8) :: FCTOR
472 DATA FCTOR /2147483647.0D0/
473 ! We get only one random number, here! DR 9/1/92
475 RNUNF= DBLE( IRAN1(1) ) / FCTOR
476 !******************************
477 ! write(6,*) 'rnunf in rnunf = ',rnunf
480 !-----------------------------------------------------------------------------
482 !-----------------------------------------------------------------------------
483 subroutine random_init(seed)
485 ! Initialize random number generator
487 ! implicit real*8 (a-h,o-z)
488 ! include 'DIMENSIONS'
492 logical :: OKRandom !, prng_restart
494 integer :: iseed_array(4)
495 integer(kind=8) :: iseed
500 integer :: ierr,error_msg
501 ! include 'COMMON.IOUNITS'
502 ! include 'COMMON.TIME1'
503 ! include 'COMMON.THREAD'
504 ! include 'COMMON.SBRIDGE'
505 ! include 'COMMON.CONTROL'
506 ! include 'COMMON.MCM'
507 ! include 'COMMON.MAP'
508 ! include 'COMMON.HEADER'
509 ! include 'COMMON.CSA'
510 ! include 'COMMON.CHAIN'
511 ! include 'COMMON.MUCA'
512 ! include 'COMMON.MD'
513 ! include 'COMMON.FFIELD'
514 ! include 'COMMON.SETUP'
515 iseed=-dint(dabs(seed))
517 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
518 'Random seed undefined. The program will stop.'
519 write (*,'(/80(1h*)/20x,a/80(1h*))') &
520 'Random seed undefined. The program will stop.'
522 call mpi_finalize(mpi_comm_world,ierr)
524 stop 'Bad random seed.'
527 if (fg_rank.eq.0) then
531 write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
532 OKRandom = prng_restart(me,iseed)
536 iseed_array(i) = dint(seed/tmp)
537 seed=seed-iseed_array(i)*tmp
540 write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',&
541 (iseed_array(i),i=1,4)
542 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',&
543 (iseed_array(i),i=1,4)
544 OKRandom = prng_restart(me,iseed_array)
548 r1=ran_number(0.0D0,1.0D0)
550 write (iout,*) 'ran_num',r1
551 if (r1.lt.0.0d0) OKRandom=.false.
553 if (.not.OKRandom) then
554 write (iout,*) 'PRNG IS NOT WORKING!!!'
555 print *,'PRNG IS NOT WORKING!!!'
558 call mpi_abort(mpi_comm_world,error_msg,ierr)
561 write (iout,*) 'too many processors for parallel prng'
562 write (*,*) 'too many processors for parallel prng'
570 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
573 end subroutine random_init
574 !-----------------------------------------------------------------------------
575 !-----------------------------------------------------------------------------