program saw implicit none logical, allocatable :: space(:,:,:) integer :: N,ntimes integer :: iwalk,istep,nfail,dir integer :: i,j,k,m,x,y,z,seed double precision, allocatable :: xi(:),xsqi(:) double precision, allocatable :: yi(:),ysqi(:) double precision, allocatable :: zi(:),zsqi(:) double precision, allocatable :: xsum(:),xsqsum(:) double precision, allocatable :: ysum(:),ysqsum(:) double precision, allocatable :: zsum(:),zsqsum(:) double precision, allocatable :: weighti(:) double precision, allocatable :: weightsum(:) double precision, allocatable :: rsqsum(:) character :: buf*80 integer, external :: iargc double precision, external :: grnd seed=235689 if (iargc() < 3) then print *,'Usage: selfavoidingwalk N ntimes seed' STOP endif call getarg(1,buf) read (unit=buf,fmt=*) N call getarg(2,buf) read (unit=buf,fmt=*) ntimes call getarg(3,buf) read (unit=buf,fmt=*) seed print *,'Doing self-avoidingwalk to',N,'steps',ntimes,' times' call sgrnd(seed) ! Allocate array which is guaranteed to fit the walk print *,'Allocating array with',(2*N+1)**3,' elements' allocate(space(-N:N,-N:N,-N:N)) allocate (xsum(N)) allocate (xsqsum(N)) allocate (ysum(N)) allocate (ysqsum(N)) allocate (zsum(N)) allocate (zsqsum(N)) allocate (xi(N)) allocate (xsqi(N)) allocate (yi(N)) allocate (ysqi(N)) allocate (zi(N)) allocate (zsqi(N)) allocate (weighti(0:N)) allocate (weightsum(0:N)) allocate (rsqsum(0:N)) do i=1,N xsum(i)=0.0; ysum(i)=0.0; zsum(i)=0.0; xsqsum(i)=0.0; ysqsum(i)=0.0; zsqsum(i)=0.0; weightsum(i)=0.0 rsqsum(i)=0.0 enddo iwalk=0 nfail=0 do ! Loop over walks ! New walk do i=-N,N do j=-N,N do k=-N,N space(i,j,k)=.false. enddo enddo enddo x=0; y=0; z=0 space(x,y,z)=.true. ! First step to the north generated here x=1 space(x,y,z)=.true. weighti(0)=1.0; weighti(1)=1.0; do istep=2,N ! Loop over steps ! Count possible directions for Rosenbluth-Rosenbluth m=0 if (.not.space(x+1,y,z)) m=m+1 if (.not.space(x-1,y,z)) m=m+1 if (.not.space(x,y+1,z)) m=m+1 if (.not.space(x,y-1,z)) m=m+1 if (.not.space(x,y,z-1)) m=m+1 if (.not.space(x,y,z+1)) m=m+1 if (m==0) then ! We fail self-avoidingness, need to redo ! Same as setting weight to 0 exit endif weighti(istep)=(m/5.0)*weighti(istep-1) do ! Loop until valid direction found (one must exist since m>0) ! Select random number direction if it is valid dir=int(grnd()*6) if (dir==0) then if (space(x+1,y,z)) cycle x=x+1 else if (dir==1) then if (space(x-1,y,z)) cycle x=x-1 else if (dir==2) then if (space(x,y+1,z)) cycle y=y+1 else if (dir==3) then if (space(x,y-1,z)) cycle y=y-1 else if (dir==4) then if (space(x,y,z+1)) cycle z=z+1 else if (dir==5) then if (space(x,y,z-1)) cycle z=z-1 endif ! If we got this far direction found exit enddo space(x,y,z)=.true. ! Store results xi(istep)=x xsqi(istep)=x*x yi(istep)=y ysqi(istep)=y*y zi(istep)=z zsqi(istep)=z*z enddo ! Check whether we made it if (istep >= N) then iwalk=iwalk+1 ! We made it so include result do i=1,N xsum(i)=xsum(i)+xi(i) xsqsum(i)=xsqsum(i)+xsqi(i) ysum(i)=ysum(i)+yi(i) ysqsum(i)=ysqsum(i)+ysqi(i) zsum(i)=zsum(i)+zi(i) zsqsum(i)=zsqsum(i)+zsqi(i) weightsum(i)=weightsum(i)+weighti(i) rsqsum(i)=rsqsum(i)+weighti(i)*(xsqi(i)+ysqi(i)+zsqi(i)) enddo else nfail=nfail+1 if (mod(nfail,10000)==0) then print *,'So far',iwalk,'walks and ',nfail,'failures' endif endif if (iwalk >= ntimes) exit enddo print *,'In total',iwalk,'walks and ',nfail,'failures' open(10,file='saw3d_rosenbluth.stat',status='unknown') do i=1,N write (10,*) i,rsqsum(i)/weightsum(i) enddo close(10) ! Write out top end only for fitting to asymptotic part open(10,file='saw3d_rosenbluth_top.stat',status='unknown') do i=10,N write (10,*) i,rsqsum(i)/weightsum(i) enddo close(10) end program saw !************************************************************************ ! Mersenne twister random number generator !************************************************************************ subroutine sgrnd(seed) ! implicit integer(a-z) ! ! Period parameters parameter(N = 624) ! dimension mt(0:N-1) ! the array for the state particle common /block/mti,mt save /block/ ! ! setting initial seeds to mt[N] using ! the generator Line 25 of Table 1 in ! [KNUTH 1981, The Art of Computer Programming ! Vol. 2 (2nd Ed.), pp102] ! mt(0)= iand(seed,-1) do 1000 mti=1,N-1 mt(mti) = iand(69069 * mt(mti-1),-1) 1000 continue ! return end subroutine sgrnd !************************************************************************ double precision function grnd() ! implicit integer(a-z) ! ! Period parameters parameter(N = 624) parameter(N1 = N+1) parameter(M = 397) parameter(MATA = -1727483681) ! constant vector a parameter(UMASK = -2147483648) ! most significant w-r bits parameter(LMASK = 2147483647) ! least significant r bits ! Tempering parameters parameter(TMASKB= -1658038656) parameter(TMASKC= -272236544) ! dimension mt(0:N-1) ! the array for the state vector common /block/mti,mt save /block/ data mti/N1/ ! mti==N+1 means mt[N] is not initialized ! dimension mag01(0:1) data mag01/0, MATA/ save mag01 ! mag01(x) = x * MATA for x=0,1 ! TSHFTU(y)=ishft(y,-11) TSHFTS(y)=ishft(y,7) TSHFTT(y)=ishft(y,15) TSHFTL(y)=ishft(y,-18) ! if(mti.ge.N) then ! generate N words at one time if(mti.eq.N+1) then ! if sgrnd() has not been called, call sgrnd(4357) ! a default initial seed is used endif ! do 1000 kk=0,N-M-1 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1))) 1000 continue do 1100 kk=N-M,N-2 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1))) 1100 continue y=ior(iand(mt(N-1),UMASK),iand(mt(0),LMASK)) mt(N-1)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1))) mti = 0 endif ! y=mt(mti) mti=mti+1 y=ieor(y,TSHFTU(y)) y=ieor(y,iand(TSHFTS(y),TMASKB)) y=ieor(y,iand(TSHFTT(y),TMASKC)) y=ieor(y,TSHFTL(y)) ! if(y.lt.0) then grnd=(dble(y)+2.0d0**32)/(2.0d0**32-1.0d0) else grnd=dble(y)/(2.0d0**32-1.0d0) endif ! return end