کد: انتخاب همه
program mont
integer,parameter::lengh=31,particle=961,montstep=75,time=10000
integer::vicinity1,vicinity2,vicinity3,i1x,i2x,j1y,j2y,h(lengh,lengh),inew,jnew,inew2,jnew2,K,l,t,x,y,i,j
real,parameter::E0=1,E1=0.3,KT=0.025 !....... tempreture 300 k..... energy write with ev
real(8)::random1,random2,random3,random4,r,varians(time),hmean,summ,sigma,E2,E3,dE
real(4)::ran2,seed,corelation1,corelation2,alfa
!***************************** initil condition ***************************
seed=.3d0
h(1:lengh,1:lengh)=0.0
!***************************** time loop **********************************
do t=1,time !(1 open)
!***************************** loop of the distribute particle ************
do k=1,particle !(2 open)
random1=ran2(seed)
random2=ran2(seed)
random1=random1*(30)+1
random2=random2*(30)+1
i=int(random1)
j=int(random2)
h(i,j)=(h(i,j)+1)
!***************************** calculte of first energy ********************
call vicinity(i,j,h,vicinity2)
E2=-(E0+vicinity2*(E1))
inew=i
jnew=j
!****************************** start monte carlo step *********************
do l=1,montstep !(3 open)
30 call random_number (random3)
if (0<random3.and.random3<0.25) then !(4 open)
i=i+1
if (i>lengh) then
i=1
end if
else if (0.25<=random3.and.random3<0.5) then
i=i-1
if (i<1) then
i=lengh
end if
else if (0.5<=random3.and.random3<0.75) then
j=j+1
if (j>lengh) then
j=1
end if
else if (0.75<=random3.and.random3<1) then
j=j-1
if (j<1) then
j=lengh
end if
!******************************* calculte of second energy *****************
call vicinity(i,j,h,vicinity3)
E3=-(E0+(vicinity3)*(E1))
dE=E3-E2
!***************** checking for accept or reject particle n new home **********
call random_number(random4)
if (dE<0.and.random4<exp(-dE/KT)) then
inew2=i
jnew2=j
else
Goto 30
end if
end if !(4 close)
end do !(3 close)
h(inew,jnew)=h(inew,jnew)-1
h(inew2,jnew2)=(h(inew2,jnew2)+1)
end do !(2 close)
!************************ calculate of varians ***********************************
summ=0
sigma=0
do i=1,lengh
do j=1,lengh
summ= summ+h(i,j)
end do
end do
hmean=summ/(lengh*lengh)
do i=1,lengh
do j=1,lengh
sigma=sigma+((h(i,j)-hmean)**2)
end do
end do
varians=sqrt((1.d0/(lengh*lengh))*sigma)
write(46,*) log(real(t)), log(varians(t))
write(47,*) varians(t)
end do !(1 close)
!************************************calculate corolation*****************
corelation1=0
do x=1,lengh
do y=1,lengh
do i=1,lengh
do j=1,lengh
corelation1=corelation1+abs(h(i,j)-h(x,y))
end do
end do
end do
end do
corelation2=(corelation1/(lengh**4))/varians(time)
alfa=log(real(corelation2))/log(real(lengh))
write (*,*) "corelation2=", corelation2
end program mont
!===================================================================================
!=====================subroutine vicinity for calculate number of naiberhood=========
subroutine vicinity(i,j,h,vicinity1)
integer,parameter::lengh=31
integer::n,i1x,i2x,j1y,j2y,j,i,h(lengh,lengh)
vicinity1=1
i1x=i+1
if (i1x>lengh) then
i1x=1
end if
if (h(i,j)<= h(i1x,j)) then
vicinity1=vicinity1+1
i2x=i-1
if (i2x<1) then
i2x=lengh
end if
if (h(i,j)<= h(i2x,j)) then
vicinity1=vicinity1+1
j1y=j+1
if (j1y>lengh) then
j1y=1
end if
if (h(i,j)<= h(i,j1y)) then
vicinity1=vicinity1+1
j2y=j-1
if (j2y<1) then
j2y=lengh
end if
if (h(i,j)<= h(i,j2y)) then
vicinity1=vicinity1+1
end if
end if
end if
end if
end subroutine vicinity
!=====================================================================================
!================================= functional of ran2 for generate random number======
FUNCTION ran2(idum)
INTEGER IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
REAL idum,ran2,AM,EPS,RNMX
PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1, &
IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791, &
NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
INTEGER idum2,j,k,iv(NTAB),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/NTAB*0/, iy/0/
if (idum.le.0) then
idum=max(-idum,1.0)
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