X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Frandom.f90;fp=source%2Funres%2Frandom.f90;h=0000000000000000000000000000000000000000;hb=2611e37f82b576a1366c2d78fce87c1a55852037;hp=fa14312fbd6d9b13a7fa7f79b1fd50a550e033d3;hpb=a1e9d5a6b704c90ebd10f86a655780212a880dce;p=unres4.git diff --git a/source/unres/random.f90 b/source/unres/random.f90 deleted file mode 100644 index fa14312..0000000 --- a/source/unres/random.f90 +++ /dev/null @@ -1,577 +0,0 @@ - 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