--- /dev/null
+ module random
+!-----------------------------------------------------------------------------
+ use io_units
+ use prng ! prng.f90 or prng_32.f90
+ use math
+ implicit none
+! public :: rndv
+!
+!-----------------------------------------------------------------------------
+ contains
+!-----------------------------------------------------------------------------
+! gen_rand_conf.F
+!-----------------------------------------------------------------------------
+ real(kind=8) function ran_number(x1,x2)
+! Calculate a random real number from the range (x1,x2).
+#ifdef MPI
+ use MPI_data ! include 'COMMON.SETUP'
+ include "mpif.h"
+#endif
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ real(kind=8) :: x1,x2,fctor
+!el local variables
+ integer :: ix(1) !el ix ---> ix(1)
+ data fctor /2147483647.0D0/
+#ifdef MPI
+ ran_number=x1+(x2-x1)*prng_next(me)
+#else
+ call vrnd(ix(1),1)
+ ran_number=x1+(x2-x1)*ix(1)/fctor
+#endif
+ return
+ end function ran_number
+!-----------------------------------------------------------------------------
+ integer function iran_num(n1,n2)
+! Calculate a random integer number from the range (n1,n2).
+#ifdef MPI
+ use MPI_data ! include 'COMMON.SETUP'
+ include "mpif.h"
+#endif
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ integer :: n1,n2,ix(1) !el ix ---> ix(1)
+ real(kind=4) :: fctor=2147483647.0
+#ifdef MPI
+ ix(1)=n1+(n2-n1+1)*prng_next(me)
+ if (ix(1).lt.n1) ix(1)=n1
+ if (ix(1).gt.n2) ix(1)=n2
+ iran_num=ix(1)
+#else
+ call vrnd(ix(1),1)
+ ix(1)=n1+(n2-n1+1)*(ix(1)/fctor)
+ if (ix(1).gt.n2) ix(1)=n2
+ iran_num=ix(1)
+#endif
+ return
+ end function iran_num
+!-----------------------------------------------------------------------------
+ real(kind=8) function binorm(x1,x2,sigma1,sigma2,ak)
+! implicit real*8 (a-h,o-z)
+!el local variables
+ real(kind=8) :: x1,x2,sigma1,sigma2,ak,alowb,aupb,seg,alen
+! print '(a)','Enter BINORM.'
+ alowb=dmin1(x1-3.0D0*sigma1,x2-3.0D0*sigma2)
+ aupb=dmax1(x1+3.0D0*sigma1,x2+3.0D0*sigma2)
+ seg=sigma1/(sigma1+ak*sigma2)
+ alen=ran_number(0.0D0,1.0D0)
+ if (alen.lt.seg) then
+ binorm=anorm_distr(x1,sigma1,alowb,aupb)
+ else
+ binorm=anorm_distr(x2,sigma2,alowb,aupb)
+ endif
+! print '(a)','Exiting BINORM.'
+ return
+ end function binorm
+!-----------------------------------------------------------------------------
+ real(kind=8) function anorm_distr(x,sigma,alowb,aupb)
+! implicit real*8 (a-h,o-z)
+! to make a normally distributed deviate with zero mean and unit variance
+!
+ integer :: iset
+ real(kind=8) :: fac,gset,rsq,v1,v2,ran1
+ real(kind=8) :: x,sigma,alowb,aupb,gaussdev
+ save iset,gset
+ data iset/0/
+!elwrite(iout,*) "anorm distr start",x,sigma,alowb,aupb
+ if(iset.eq.0) then
+ 1 v1=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
+ v2=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
+ rsq=v1**2+v2**2
+ if(rsq.ge.1.d0.or.rsq.eq.0.0d0) goto 1
+ fac=sqrt(-2.0d0*log(rsq)/rsq)
+ gset=v1*fac
+ gaussdev=v2*fac
+ iset=1
+ else
+ gaussdev=gset
+ iset=0
+ endif
+ anorm_distr=x+gaussdev*sigma
+!elwrite(iout,*) "anorm distr end",x,sigma,alowb,aupb,anorm_distr
+ return
+ end function anorm_distr
+!-----------------------------------------------------------------------------
+ subroutine mult_norm(lda,n,a,x,fail)
+!
+! Generate the vector X whose elements obey the multiple-normal distribution
+! from exp(-0.5*X'AX). LDA is the leading dimension of the moment matrix A,
+! n is the dimension of the problem. FAIL is set at .TRUE., if the smallest
+! eigenvalue of the matrix A is close to 0.
+!
+! implicit double precision (a-h,o-z)
+ integer :: lda,n,i,j
+ real(kind=8),dimension(lda,n) :: a
+ real(kind=8),dimension(n) :: x
+ real(kind=8),dimension(100) :: eig,work
+ real(kind=8),dimension(3,3) :: vec
+ real(kind=8) :: eig_limit=1.0D-8
+ logical :: fail
+ real(kind=8) :: sigma,alim,xi
+ fail=.false.
+! print '(a)','Enter MULT_NORM.'
+!
+! Find the smallest eigenvalue of the matrix A.
+!
+! do i=1,n
+! print '(8f10.5)',(a(i,j),j=1,n)
+! enddo
+#ifdef NAG
+ call f02faf('V','U',2,a,lda,eig,work,100,ifail)
+#else
+ call djacob(2,lda,10000,1.0d-10,a,vec,eig)
+#endif
+! print '(8f10.5)',(eig(i),i=1,n)
+! print '(a)'
+! do i=1,n
+! print '(8f10.5)',(a(i,j),j=1,n)
+! enddo
+ if (eig(1).lt.eig_limit) then
+ print *,'From Mult_Norm: Eigenvalues of A are too small.'
+ fail=.true.
+ return
+ endif
+!
+! Generate points following the normal distributions along the principal
+! axes of the moment matrix. Store in WORK.
+!
+ do i=1,n
+ sigma=1.0D0/dsqrt(eig(i))
+ alim=-3.0D0*sigma
+ work(i)=anorm_distr(0.0D0,sigma,-alim,alim)
+ enddo
+!
+! Transform the vector of normal variables back to the original basis.
+!
+ do i=1,n
+ xi=0.0D0
+ do j=1,n
+ xi=xi+a(i,j)*work(j)
+ enddo
+ x(i)=xi
+ enddo
+ return
+ end subroutine mult_norm
+!-----------------------------------------------------------------------------
+ subroutine mult_norm1(lda,n,a,z,box,x,fail)
+!
+! Generate the vector X whose elements obey the multi-gaussian multi-dimensional
+! distribution from sum_{i=1}^m W(i)exp[-0.5*X'(i)A(i)X(i)]. LDA is the
+! leading dimension of the moment matrix A, n is the dimension of the
+! distribution, nlob is the number of lobes. FAIL is set at .TRUE., if the
+! smallest eigenvalue of the matrix A is close to 0.
+!
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+#ifdef MP
+ use MPI_data !include 'COMMON.SETUP'
+#endif
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ integer :: lda,n,i,j,istep
+ real(kind=8),dimension(lda,n) :: a
+ real(kind=8),dimension(n) :: z,x
+ real(kind=8),dimension(n,n) :: box
+ real(kind=8) :: etmp,ww,xi,dec
+! include 'COMMON.IOUNITS'
+ logical :: fail
+!
+! Generate points following the normal distributions along the principal
+! axes of the moment matrix. Store in WORK.
+!
+!d print *,'CG Processor',me,' entered MultNorm1.'
+!d print '(2(1pe14.4),3x,1pe14.4)',((a(i,j),j=1,2),z(i),i=1,2)
+!d do i=1,n
+!d print *,i,box(1,i),box(2,i)
+!d enddo
+ istep = 0
+ 10 istep = istep + 1
+ if (istep.gt.10000) then
+! write (iout,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
+! & ' in MultNorm1.'
+! write (*,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
+! & ' in MultNorm1.'
+! write (iout,*) 'box',box
+! write (iout,*) 'a',a
+! write (iout,*) 'z',z
+ fail=.true.
+ return
+ endif
+ do i=1,n
+ x(i)=ran_number(box(1,i),box(2,i))
+ enddo
+ ww=0.0D0
+ do i=1,n
+ xi=pinorm(x(i)-z(i))
+ ww=ww+0.5D0*a(i,i)*xi*xi
+ do j=i+1,n
+ ww=ww+a(i,j)*xi*pinorm(x(j)-z(j))
+ enddo
+ enddo
+ dec=ran_number(0.0D0,1.0D0)
+! print *,(x(i),i=1,n),ww,dexp(-ww),dec
+!rc if (dec.gt.dexp(-ww)) goto 10
+ if(-ww.lt.100) then
+ etmp=dexp(-ww)
+ else
+ return
+ endif
+ if (dec.gt.etmp) goto 10
+!d print *,'CG Processor',me,' exitting MultNorm1.'
+ return
+ end subroutine mult_norm1
+!-----------------------------------------------------------------------------
+! ran.f
+!-----------------------------------------------------------------------------
+ function ran0(idum)
+ integer :: idum
+ integer,parameter :: IA=16807,IM=2147483647,IQ=127773,IR=2836
+ integer,parameter :: MASK=123459876
+ real(kind=4) :: ran0
+ real(kind=4),parameter :: AM=1./IM
+ integer :: k
+ idum=ieor(idum,MASK)
+ k=idum/IQ
+ idum=IA*(idum-k*IQ)-IR*k
+ if (idum.lt.0) idum=idum+IM
+ ran0=AM*idum
+ idum=ieor(idum,MASK)
+ return
+ end function ran0
+! (C) Copr. 1986-92 Numerical Recipes Software *11915
+!-----------------------------------------------------------------------------
+ function ran1(idum)
+ integer :: idum
+ integer,parameter :: IA=16807,IM=2147483647,IQ=127773,IR=2836
+ integer,parameter :: NTAB=32,NDIV=1+(IM-1)/NTAB
+ real(kind=4) :: ran1
+ real(kind=4),parameter :: AM=1./IM,EPS=1.2e-7,RNMX=1.-EPS
+ integer :: j,k,iy
+ integer,dimension(NTAB) :: iv
+ SAVE iv,iy
+ DATA iv /NTAB*0/, iy /0/
+ if (idum.le.0.or.iy.eq.0) then
+ idum=max(-idum,1)
+ do 11 j=NTAB+8,1,-1
+ k=idum/IQ
+ idum=IA*(idum-k*IQ)-IR*k
+ if (idum.lt.0) idum=idum+IM
+ if (j.le.NTAB) iv(j)=idum
+11 continue
+ iy=iv(1)
+ endif
+ k=idum/IQ
+ idum=IA*(idum-k*IQ)-IR*k
+ if (idum.lt.0) idum=idum+IM
+ j=1+iy/NDIV
+ iy=iv(j)
+ iv(j)=idum
+ ran1=min(AM*iy,RNMX)
+ return
+ end function ran1
+! (C) Copr. 1986-92 Numerical Recipes Software *11915
+!-----------------------------------------------------------------------------
+ function ran2(idum)
+ integer :: idum
+ integer,parameter :: IM1=2147483563,IM2=2147483399,IMM1=IM1-1,&
+ IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,&
+ NTAB=32,NDIV=1+IMM1/NTAB
+ real(kind=4) :: ran2
+ real(kind=4),parameter :: AM=1./IM1,EPS=1.2e-7,RNMX=1.-EPS
+ integer :: idum2,j,k,iy
+ integer,dimension(NTAB) :: iv
+ SAVE iv,iy,idum2
+ DATA idum2/123456789/, iv/NTAB*0/, iy/0/
+ if (idum.le.0) then
+ idum=max(-idum,1)
+ idum2=idum
+ do 11 j=NTAB+8,1,-1
+ k=idum/IQ1
+ idum=IA1*(idum-k*IQ1)-k*IR1
+ if (idum.lt.0) idum=idum+IM1
+ if (j.le.NTAB) iv(j)=idum
+11 continue
+ iy=iv(1)
+ endif
+ k=idum/IQ1
+ idum=IA1*(idum-k*IQ1)-k*IR1
+ if (idum.lt.0) idum=idum+IM1
+ k=idum2/IQ2
+ idum2=IA2*(idum2-k*IQ2)-k*IR2
+ if (idum2.lt.0) idum2=idum2+IM2
+ j=1+iy/NDIV
+ iy=iv(j)-idum2
+ iv(j)=idum
+ if(iy.lt.1)iy=iy+IMM1
+ ran2=min(AM*iy,RNMX)
+ return
+ end function ran2
+! (C) Copr. 1986-92 Numerical Recipes Software *11915
+!-----------------------------------------------------------------------------
+ function ran3(idum)
+ integer :: idum
+ integer,parameter :: MBIG=1000000000,MSEED=161803398,MZ=0
+! REAL MBIG,MSEED,MZ
+ real(kind=4) :: ran3
+ real(kind=4),parameter :: FAC=1./MBIG
+! PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG)
+ integer :: i,iff,ii,inext,inextp,k
+ integer :: mj,mk
+ integer,dimension(55) :: ma
+! REAL mj,mk,ma(55)
+ SAVE iff,inext,inextp,ma
+ DATA iff /0/
+ if(idum.lt.0.or.iff.eq.0)then
+ iff=1
+ mj=MSEED-iabs(idum)
+ mj=mod(mj,MBIG)
+ ma(55)=mj
+ mk=1
+ do 11 i=1,54
+ ii=mod(21*i,55)
+ ma(ii)=mk
+ mk=mj-mk
+ if(mk.lt.MZ)mk=mk+MBIG
+ mj=ma(ii)
+11 continue
+ do 13 k=1,4
+ do 12 i=1,55
+ ma(i)=ma(i)-ma(1+mod(i+30,55))
+ if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG
+12 continue
+13 continue
+ inext=0
+ inextp=31
+ idum=1
+ endif
+ inext=inext+1
+ if(inext.eq.56)inext=1
+ inextp=inextp+1
+ if(inextp.eq.56)inextp=1
+ mj=ma(inext)-ma(inextp)
+ if(mj.lt.MZ)mj=mj+MBIG
+ ma(inext)=mj
+ ran3=mj*FAC
+ return
+ end function ran3
+! (C) Copr. 1986-92 Numerical Recipes Software *11915
+!-----------------------------------------------------------------------------
+! randgens.f
+!-----------------------------------------------------------------------------
+! $Date: 1994/10/04 16:19:52 $
+! $Revision: 2.1 $
+!
+!
+! See help for RANDOMV on the PSFSHARE disk to understand these
+! subroutines. This is the VS Fortran version of this code.
+!
+!
+ subroutine VRND(VEC,N)
+
+ use comm_vrandd
+!el integer,dimension(250) :: A
+ integer :: LOOP,N !el,I,I147
+ integer,dimension(N) :: VEC
+!el COMMON /VRANDD/ A, I, I147
+ DO 23000 LOOP=1,N
+ I=I+1
+ IF(.NOT.(I.GE.251))GOTO 23002
+ I=1
+23002 CONTINUE
+ I147=I147+1
+ IF(.NOT.(I147.GE.251))GOTO 23004
+ I147=1
+23004 CONTINUE
+ A(I)=IEOR(A(I147),A(I))
+ VEC(LOOP)=A(I)
+23000 CONTINUE
+ return
+ end subroutine VRND
+!-----------------------------------------------------------------------------
+ real(kind=8) function RNDV(IDUM)
+
+ real(kind=8) :: RM1,RM2,R(99)
+ INTEGER :: IA1,IC1,M1, IA2,IC2,M2, IA3,IC3,M3, IDUM
+ INTEGER :: IX1,IX2,IX3,J
+ SAVE
+ DATA IA1,IC1,M1/1279,351762,1664557/
+ DATA IA2,IC2,M2/2011,221592,1048583/
+ DATA IA3,IC3,M3/15551,6150,29101/
+ IF(.NOT.(IDUM.LT.0))GOTO 23006
+ IX1 = MOD(-IDUM,M1)
+ IX1 = MOD(IA1*IX1+IC1,M1)
+ IX2 = MOD(IX1,M2)
+ IX1 = MOD(IA1*IX1+IC1,M1)
+ IX3 = MOD(IX1,M3)
+ RM1 = 1./DBLE(M1)
+ RM2 = 1./DBLE(M2)
+ DO 23008 J = 1,99
+ IX1 = MOD(IA1*IX1+IC1,M1)
+ IX2 = MOD(IA2*IX2+IC2,M2)
+ R(J) = (DBLE(IX1)+DBLE(IX2)*RM2)*RM1
+23008 CONTINUE
+23006 CONTINUE
+ IX1 = MOD(IA1*IX1+IC1,M1)
+ IX2 = MOD(IA2*IX2+IC2,M2)
+ IX3 = MOD(IA3*IX3+IC3,M3)
+ J = 1+(99*IX3)/M3
+ RNDV = R(J)
+ R(J) = (DBLE(IX1)+DBLE(IX2)*RM2)*RM1
+ IDUM = IX1
+ return
+ end function RNDV
+!-----------------------------------------------------------------------------
+ subroutine VRNDST(SEED)
+
+ use comm_vrandd
+!el INTEGER :: A(250)
+ INTEGER :: LOOP,IDUM !el,I,I147
+ INTEGER :: SEED
+!el real(kind=8) :: RNDV
+!el COMMON /VRANDD/ A, I, I147
+ I=0
+ I147=103
+ IDUM=SEED
+ DO 23010 LOOP=1,250
+ A(LOOP)=INT(RNDV(IDUM)*2147483647)
+23010 CONTINUE
+ return
+ end subroutine VRNDST
+!-----------------------------------------------------------------------------
+ subroutine VRNDIN(IODEV)
+ use comm_vrandd
+ INTEGER :: IODEV !el, A(250),I, I147
+!el COMMON/VRANDD/ A, I, I147
+ READ(IODEV) A, I, I147
+ return
+ end subroutine VRNDIN
+!-----------------------------------------------------------------------------
+ subroutine VRNDOU(IODEV)
+! This corresponds to VRNDOUT in the APFTN64 version
+ use comm_vrandd
+ INTEGER :: IODEV !el, A(250),I, I147
+!el COMMON/VRANDD/ A, I, I147
+ WRITE(IODEV) A, I, I147
+ return
+ end subroutine VRNDOU
+!-----------------------------------------------------------------------------
+ real(kind=8) function RNUNF(N)
+ INTEGER :: IRAN1(2000),N
+ real(kind=8) :: FCTOR
+ DATA FCTOR /2147483647.0D0/
+! We get only one random number, here! DR 9/1/92
+ CALL VRND(IRAN1,1)
+ RNUNF= DBLE( IRAN1(1) ) / FCTOR
+!******************************
+! write(6,*) 'rnunf in rnunf = ',rnunf
+ return
+ end function RNUNF
+!-----------------------------------------------------------------------------
+! readrtns_CSA.F
+!-----------------------------------------------------------------------------
+ subroutine random_init(seed)
+!
+! Initialize random number generator
+!
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+#ifdef MPI
+ use MPI_data
+ include 'mpif.h'
+ logical :: OKRandom !, prng_restart
+ real(kind=8) :: r1
+ integer :: iseed_array(4)
+ integer(kind=8) :: iseed
+#else
+ integer :: iseed
+#endif
+ real(kind=8) :: seed
+ integer :: ierr,error_msg
+! include 'COMMON.IOUNITS'
+! include 'COMMON.TIME1'
+! include 'COMMON.THREAD'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.CONTROL'
+! include 'COMMON.MCM'
+! include 'COMMON.MAP'
+! include 'COMMON.HEADER'
+! include 'COMMON.CSA'
+! include 'COMMON.CHAIN'
+! include 'COMMON.MUCA'
+! include 'COMMON.MD'
+! include 'COMMON.FFIELD'
+! include 'COMMON.SETUP'
+ iseed=-dint(dabs(seed))
+ if (iseed.eq.0) then
+ write (iout,'(/80(1h*)/20x,a/80(1h*))') &
+ 'Random seed undefined. The program will stop.'
+ write (*,'(/80(1h*)/20x,a/80(1h*))') &
+ 'Random seed undefined. The program will stop.'
+#ifdef MPI
+ call mpi_finalize(mpi_comm_world,ierr)
+#endif
+ stop 'Bad random seed.'
+ endif
+#ifdef MPI
+ if (fg_rank.eq.0) then
+ seed=seed*(me+1)+1
+#ifdef AMD64
+ if(me.eq.king) &
+ write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
+ OKRandom = prng_restart(me,iseed)
+#else
+ do i=1,4
+ tmp=65536.0d0**(4-i)
+ iseed_array(i) = dint(seed/tmp)
+ seed=seed-iseed_array(i)*tmp
+ enddo
+ if(me.eq.king) &
+ write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',&
+ (iseed_array(i),i=1,4)
+ write (*,*) 'MPI: node= ',me, ' iseed(4)= ',&
+ (iseed_array(i),i=1,4)
+ OKRandom = prng_restart(me,iseed_array)
+#endif
+ if (OKRandom) then
+! r1 = prng_next(me)
+ r1=ran_number(0.0D0,1.0D0)
+ if(me.eq.king) &
+ write (iout,*) 'ran_num',r1
+ if (r1.lt.0.0d0) OKRandom=.false.
+ endif
+ if (.not.OKRandom) then
+ write (iout,*) 'PRNG IS NOT WORKING!!!'
+ print *,'PRNG IS NOT WORKING!!!'
+ if (me.eq.0) then
+ call flush(iout)
+ call mpi_abort(mpi_comm_world,error_msg,ierr)
+ stop
+ else
+ write (iout,*) 'too many processors for parallel prng'
+ write (*,*) 'too many processors for parallel prng'
+ call flush(iout)
+ stop
+ endif
+ endif
+ endif
+#else
+ call vrndst(iseed)
+ write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
+#endif
+ return
+ end subroutine random_init
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+ end module random
projects / unres4.git / blobdiff