! ! ...1.........2.........3.........4.........5.........6.........6.........7.........8 program main implicit none integer(8),parameter :: n1x=50 integer(8),parameter :: n2x=50 integer(8),parameter :: nclasses=10 real(8) ,parameter :: jint=1.d0 ! spin-spin interaction real(8) ,parameter :: hext=0.d0 ! external field real(8) ,parameter :: kbt=2.3d0 integer(8),parameter :: nstep=10**7 ! number of mc steps real(8) ,parameter :: tslice=1.d+5 integer(8) :: nmem(nclasses,n2x) integer(8) :: nmem1 integer(8) :: state(n1x,n2x) integer(8) :: class(n1x,n2x) integer(8) :: i,i1,i2,istep,icount integer(8) :: i1p,i1m,i2p,i2m,nn,imem,isvar integer(8) :: iclass,thisclass integer(8) :: nslice integer(8) :: isteplast real(8) :: x(nclasses) real(8) :: p1(nclasses) real(8) :: b real(8) :: deltae(nclasses) real(8) :: svar real(8) :: acceptanceratio real(8) :: etot,etotav real(8) :: mag,magav real(8) :: ran real(8) :: time real(8) :: deltat real(8) :: profile(nclasses) real(8) :: weight ! ************************************************************************************ weight=1.d0/real(n1x*n2x) ! ! ==================================================================================== ! == set up process table == ! == each process class refers to a spin flip in a given environment == ! == each process class is fully defined by the central state and the == ! == average orientation of the nearest neigbors == ! ==================================================================================== ! == class= 1: central state 0; all neighbors 0 == ! == class= 2: central state 0; 1 neighbor with state 1 == ! == class= 3: central state 0; 2 neighbors with state 1 == ! == class= 4: central state 0; 3 neighbors with state 1 == ! == class= 5: central state 0; 4 neighbors with state 1 == ! == class= 6: central state 1; all neighbors 0 == ! == class= 7: central state 1; 1 neighbor with state 1 == ! == class= 8: central state 1; 2 neighbors with state 1 == ! == class= 9: central state 1; 3 neighbors with state 1 == ! == class=10: central state 1; 4 neighbors with state 1 == ! ==================================================================================== do i=1,nclasses isvar=int(real(i-1)/5.d0) ! state of central atom nn=i-1-5*isvar ! sum of spin up states on nearest neighbors deltae(i)=-2.d0*real(2*isvar-1)*(-hext-jint*real(2*nn-4)) p1(i)=min(1,exp(-deltae(i)/kbt)) ! probability for a spin flip enddo p1(:)=p1(:)/sum(p1) ! relative jump rates for the different process classes ! ! ==================================================================================== ! == prepare initial state == ! ==================================================================================== do i1=1,n1x do i2=1,n2x call random_number(ran) state(i1,i2)=nint(ran) enddo enddo ! ! ==================================================================================== ! == calculate initial properties == ! ==================================================================================== call totalenergy(n1x,n2x,hext,jint,state,etot,mag) print*,'magnetization',mag ! ! ==================================================================================== ! == set up process table == ! == class assigns each site a process class == ! == nmem counts the number of sites in each class in a vertical line == ! ==================================================================================== nmem(:,:)=0 do i1=1,n1x do i2=1,n2x i1p=1+modulo(i1,n1x) i1m=1+modulo(i1-2,n1x) i2p=1+modulo(i2,n2x) i2m=1+modulo(i2-2,n2x) nn=state(i1p,i2)+state(i1m,i2)+state(i1,i2p)+state(i1,i2m) iclass=1+nn+5*state(i1,i2) nmem(iclass,i2)=nmem(iclass,i2)+1 class(i1,i2)=iclass enddo enddo ! ! ==================================================================================== ! == monte carlo loop == ! ==================================================================================== open(10,file='nfold.dat') open(11,file='nfold.out') time=0.d0 magav=0.d0 etotav=0.d0 nslice=0 isteplast=0 do istep=1,nstep if(modulo(istep,nstep/100_8).eq.0)print*,'percent finished ',100*istep/nstep ! ! ================================================================================== ! == select a class == ! ================================================================================== svar=0.d0 do i=1,nclasses svar=svar+p1(i)*real(sum(nmem(i,:))) x(i)=svar enddo x=x/svar ! == the first class corresponds to the interval [0,x(1)], ! == the second class corresponds to the interval [x(1),x(2)], ! == the last class corresponds to the interval [x(n-1),x(n)], where x(n)=1. call random_number(ran) do i=1,nclasses if(x(i).gt.ran) then thisclass=i ! process class "thisclass" is selected exit end if enddo ! ! ================================================================================== ! == select a random site from this class == ! ================================================================================== nmem1=sum(nmem(thisclass,:)) call random_number(ran) isvar=nint(0.5d0+ran*real(nmem1-1)) ! process isvar in this class is selected isvar=min(isvar,nmem1) ! atom number in this class isvar=max(isvar,1) ! atom number in this class ! icount=0 do i2=1,n2x if(icount+nmem(thisclass,i2).lt.isvar) then icount=icount+nmem(thisclass,i2) else do i1=1,n1x if(class(i1,i2).eq.thisclass) then icount=icount+1 if(icount.eq.isvar) then goto 1000 end if end if enddo end if enddo 1000 continue ! ! ================================================================================== ! == increment the time == ! ================================================================================== call random_number(ran) deltat=-log(ran)/p1(i) time=time+deltat ! ! ================================================================================== ! == update average values == ! ================================================================================== magav=magav+mag*deltat etotav=etotav+etot*deltat ! ! ================================================================================== ! == flip spin == ! ================================================================================== state(i1,i2)=1-state(i1,i2) mag=mag+2*(2*state(i1,i2)-1)*weight etot=etot+deltae(thisclass) ! ! ================================================================================== ! == update process list == ! ================================================================================== i1p=1+modulo(i1,n1x) i1m=1+modulo(i1-2,n1x) i2p=1+modulo(i2,n2x) i2m=1+modulo(i2-2,n2x) nmem(class(i1,i2),i2)=nmem(class(i1,i2),i2)-1 nmem(class(i1m,i2),i2)=nmem(class(i1m,i2),i2)-1 nmem(class(i1p,i2),i2)=nmem(class(i1p,i2),i2)-1 nmem(class(i1,i2m),i2m)=nmem(class(i1,i2m),i2m)-1 nmem(class(i1,i2p),i2p)=nmem(class(i1,i2p),i2p)-1 if(state(i1,i2).eq.1) then class(i1,i2)=class(i1,i2)+5 class(i1p,i2)=class(i1p,i2)+1 class(i1m,i2)=class(i1m,i2)+1 class(i1,i2p)=class(i1,i2p)+1 class(i1,i2m)=class(i1,i2m)+1 else class(i1,i2)=class(i1,i2)-5 class(i1p,i2)=class(i1p,i2)-1 class(i1m,i2)=class(i1m,i2)-1 class(i1,i2p)=class(i1,i2p)-1 class(i1,i2m)=class(i1,i2m)-1 end if nmem(class(i1,i2),i2)=nmem(class(i1,i2),i2)+1 nmem(class(i1m,i2),i2)=nmem(class(i1m,i2),i2)+1 nmem(class(i1p,i2),i2)=nmem(class(i1p,i2),i2)+1 nmem(class(i1,i2m),i2m)=nmem(class(i1,i2m),i2m)+1 nmem(class(i1,i2p),i2p)=nmem(class(i1,i2p),i2p)+1 ! ! ================================================================================== ! == write result == ! ================================================================================== if(time.gt.tslice) then svar=1.d0/tslice etotav=etotav*svar magav=magav*svar acceptanceratio=time/real(istep-isteplast) write(10,*)nslice*tslice,etot*weight,etotav*weight,mag,magav,acceptanceratio write(11,*)nslice*tslice,etot*weight,mag call plotstate(11,n1x,n2x,state) isteplast=istep nslice=nslice+1 time=time-tslice magav=0.d0 etotav=0.d0 end if enddo close(10) ! ! ==================================================================================== ! == analyze result == ! ==================================================================================== close(11) print*,'magnetization',2.d0*sum(state)/real(n1x*n2x)-1.d0 stop end ! ! ...1.........2.........3.........4.........5.........6.........6.........7.........8 subroutine totalenergy(n1x,n2x,hext,jint,state,etot,mag) implicit none integer(8),intent(in) :: n1x integer(8),intent(in) :: n2x real(8) ,intent(in) :: hext real(8) ,intent(in) :: jint integer(8),intent(in) :: state(n1x,n2x) real(8) ,intent(out):: etot real(8) ,intent(out):: mag integer(8) :: i1,i2 integer(8) :: i1p,i1m,i2p,i2m integer(8) :: nn real(8) :: hint ! ************************************************************************************ etot=0.d0 mag=0.d0 do i1=1,n1x do i2=1,n2x i1p=1+modulo(i1,n1x) i1m=1+modulo(i1-2,n1x) i2p=1+modulo(i2,n2x) i2m=1+modulo(i2-2,n2x) nn=state(i1p,i2)+state(i1m,i2)+state(i1,i2p)+state(i1,i2m) hint=jint*real(2*nn-4) ! sum_k jint*sigma_k etot=etot-(hext+0.5d0*hint)*real(2*state(i1,i2)-1) mag=mag+real(2*state(i1,i2)-1) enddo enddo mag=mag/real(n1x*n2x) return end ! ! ...1.........2.........3.........4.........5.........6.........6.........7.........8 subroutine plotstate(nfil,n1x,n2x,state) ! ************************************************************************************ implicit none integer ,intent(in) :: nfil integer(8),intent(in) :: n1x integer(8),intent(in) :: n2x integer(8),intent(in) :: state(n1x,n2x) integer(8) :: i1,i2 character(n2x) :: string ! ************************************************************************************ do i2=1,n2x string(i2:i2)='-' enddo write(nfil,*)'|'//string//'|' do i1=1,n1x string='' do i2=1,n2x if(state(i1,i2).eq.1)string(i2:i2)='o' enddo write(nfil,*)'|'//string//'|' enddo do i2=1,n2x string(i2:i2)='-' enddo write(nfil,*)'|'//string//'|' return end