!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Fortran90 program, Monte Carlo simulation of the Ising model!programmed by Hans-Marc Erkinger, erkinger@itp.tu-graz.ac.at!Final version, 21/08/2000!Compiled with Fortran 90 Lahey-Compiler on RedHat Linux!f95-lah --wide --mod MOD_lahey --pca --tpp --warn --trace -O!NEEDS (in same directory):!condor.dat!input_pl.dat!--------------------------------------------------------------!Purpose of program:!Simulates the Ising model with a local update algorithm in the!high-temperature representation.!Algorithm shows very large correlations, comparable to a !Metropolis-algorithm with z >= 2 for all dimensions observed!(2d, 3d, 4d)!Further information can be found in my diploma thesis:!"A new cluster algorithm for the Ising model", !TU Graz, August 2000!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Module for creating random numbers!---------------------------------------------------MODULE random_mod  !---------------------------------------------------  IMPLICIT NONE  INTEGER iseedCONTAINS  !---------------------------------------------------  SUBROUTINE ranset(iseed)    !---------------------------------------------------    INTEGER iseed,i(1)    i(1)=iseed    CALL RANDOM_SEED(put=i(1:1))  END SUBROUTINE ranset  !---------------------------------------------------  SUBROUTINE ranget(iseed)    !---------------------------------------------------    INTEGER iseed,i(1)    CALL RANDOM_SEED(get=i(1:1))    iseed=i(1)  END SUBROUTINE ranget  !---------------------------------------------------  FUNCTION random()    !---------------------------------------------------    REAL*8 random,x(1)    CALL RANDOM_NUMBER(x(1:1)) !internal rng in F90    random = x(1)  END FUNCTION randomEND MODULE random_mod!---------------------------------------------------!Utility module for creating next neighbours (nn),!list of plaquettes(p) and corresponding bonds (b)!belonging to a plaquette!---------------------------------------------------MODULE utility  !---------------------------------------------------  IMPLICIT NONECONTAINS  !---------------------------------------------------  SUBROUTINE create_nn_2d(nx,ny,v,b,nn,p)    INTEGER, INTENT(in) :: nx, ny, v    INTEGER, DIMENSION(1:v,1:4), INTENT(inout) :: nn    INTEGER,DIMENSION(1:v,1:2), INTENT(in) :: b    INTEGER,DIMENSION(1:v,1:4), INTENT(inout) :: p    INTEGER, DIMENSION(:,:), ALLOCATABLE :: sites    INTEGER :: i,j    ALLOCATE(sites(1:ny+2,1:nx+2))    sites=0    DO i=1,ny       DO j=1,nx          sites(i+1,j+1)=(i-1)*nx+j       ENDDO    ENDDO    sites(1,:)=sites(ny+1,:)    sites(ny+2,:)=sites(2,:)    sites(:,1)=sites(:,nx+1)    sites(:,nx+2)=sites(:,2)    DO i=1,ny       DO j=1,nx          nn((i-1)*nx+j,1)=sites(i+1,j)       !-x dir          nn((i-1)*nx+j,2)=sites(i+1,j+2)     !+x dir          nn((i-1)*nx+j,3)=sites(i,j+1)       !-y dir           nn((i-1)*nx+j,4)=sites(i+2,j+1)     !+y dir        ENDDO    ENDDO    !create a list of plaquettes and corresponding bonds    DO i = 1,v       p(i,1)=b(i,1)       p(i,2)=b(i,2)       p(i,3)=b(nn(i,2),2)       p(i,4)=b(nn(i,4),1)    ENDDO    RETURN;  END SUBROUTINE create_nn_2d  !---------------------------------------------------  SUBROUTINE create_nn_3d(nx,ny,nz,v,b,nn,p)    INTEGER, INTENT(in) :: nx, ny, nz, v    INTEGER, DIMENSION(1:v,1:6), INTENT(inout) :: nn    INTEGER,DIMENSION(1:v,1:3), INTENT(in) :: b    INTEGER,DIMENSION(1:3*v,1:4), INTENT(inout) :: p    INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: sites    INTEGER :: i,j,k,csite    ALLOCATE(sites(1:nz+2,1:ny+2,1:nx+2))    sites=0    DO i=1,nz       DO j=1,ny          DO k=1,nx             sites(i+1,j+1,k+1) = (i-1)*ny*nx + (j-1)*nx + k          ENDDO       ENDDO    ENDDO    sites(1,:,:)=sites(nz+1,:,:)    sites(nz+2,:,:)=sites(2,:,:)    sites(:,1,:)=sites(:,ny+1,:)    sites(:,ny+2,:)=sites(:,2,:)    sites(:,:,1)=sites(:,:,nx+1)    sites(:,:,nx+2)=sites(:,:,2)    DO i=1,nz       DO j=1,ny          DO k=1,nx             csite=(i-1)*nx*ny+(j-1)*nx+k             nn(csite,1)=sites(i+1,j+1,k)      !-x dir             nn(csite,2)=sites(i+1,j+1,k+2)    !+x dir              nn(csite,3)=sites(i+1,j,k+1)      !-y dir             nn(csite,4)=sites(i+1,j+2,k+1)    !+y dir             nn(csite,5)=sites(i,j+1,k+1)      !-z dir             nn(csite,6)=sites(i+2,j+1,k+1)    !+z dir          ENDDO       ENDDO    ENDDO    DO i = 1,v       !the 1,2-case       p(3*(i-1)+1,1)=b(i,1)       p(3*(i-1)+1,2)=b(i,2)       p(3*(i-1)+1,3)=b(nn(i,2),2)       p(3*(i-1)+1,4)=b(nn(i,4),1)       !the 1,3-case       p(3*(i-1)+2,1)=b(i,1)       p(3*(i-1)+2,2)=b(i,3)       p(3*(i-1)+2,3)=b(nn(i,2),3)       p(3*(i-1)+2,4)=b(nn(i,6),1)       !the 2,3-case       p(3*(i-1)+3,1)=b(i,2)       p(3*(i-1)+3,2)=b(i,3)       p(3*(i-1)+3,3)=b(nn(i,4),3)       p(3*(i-1)+3,4)=b(nn(i,6),2)    ENDDO    RETURN;  END SUBROUTINE create_nn_3d  !---------------------------------------------------  SUBROUTINE create_nn_4d(nx,ny,nz,nw,v,b,nn,p)    INTEGER, INTENT(in) :: nx, ny, nz, nw, v    INTEGER, DIMENSION(1:v,1:8), INTENT(inout) :: nn    INTEGER,DIMENSION(1:v,1:4), INTENT(in) :: b    INTEGER,DIMENSION(1:6*v,1:4), INTENT(inout) :: p    INTEGER, DIMENSION(:,:,:,:), ALLOCATABLE :: sites    INTEGER :: i,j,k,l,csite    ALLOCATE(sites(1:nw+2,1:nz+2,1:ny+2,1:nx+2))    sites=0    DO i=1,nw       DO j=1,nz          DO k=1,ny             DO l=1,nx                sites(i+1,j+1,k+1,l+1) = (i-1)*nz*ny*nx + (j-1)*ny*nx + (k-1)*nx + l             ENDDO          ENDDO       ENDDO    ENDDO    sites(1,:,:,:)=sites(nw+1,:,:,:)    sites(nw+2,:,:,:)=sites(2,:,:,:)    sites(:,1,:,:)=sites(:,nz+1,:,:)    sites(:,nz+2,:,:)=sites(:,2,:,:)    sites(:,:,1,:)=sites(:,:,ny+1,:)    sites(:,:,ny+2,:)=sites(:,:,2,:)    sites(:,:,:,1)=sites(:,:,:,nx+1)    sites(:,:,:,nx+2)=sites(:,:,:,2)    DO i=1,nw       DO j=1,nz          DO k=1,ny             DO l=1,nx                csite=(i-1)*nz*ny*nx + (j-1)*ny*nx + (k-1)*nx + l                nn(csite,1)=sites(i+1,j+1,k+1,l)      !-x dir                nn(csite,2)=sites(i+1,j+1,k+1,l+2)    !+x dir                 nn(csite,3)=sites(i+1,j+1,k,l+1)      !-y dir                nn(csite,4)=sites(i+1,j+1,k+2,l+1)    !+y dir                nn(csite,5)=sites(i+1,j,k+1,l+1)      !-z dir                nn(csite,6)=sites(i+1,j+2,k+1,l+1)    !+z dir                nn(csite,7)=sites(i,j+1,k+1,l+1)      !-w dir                nn(csite,8)=sites(i+2,j+1,k+1,l+1)    !+w dir             ENDDO          ENDDO       ENDDO    ENDDO    DO i = 1,v       !the 1,2-case       p(6*(i-1)+1,1)=b(i,1)       p(6*(i-1)+1,2)=b(i,2)       p(6*(i-1)+1,3)=b(nn(i,2),2)       p(6*(i-1)+1,4)=b(nn(i,4),1)       !the 1,3-case       p(6*(i-1)+2,1)=b(i,1)       p(6*(i-1)+2,2)=b(i,3)       p(6*(i-1)+2,3)=b(nn(i,2),3)       p(6*(i-1)+2,4)=b(nn(i,6),1)       !the 1,4-case       p(6*(i-1)+3,1)=b(i,1)       p(6*(i-1)+3,2)=b(i,4)       p(6*(i-1)+3,3)=b(nn(i,2),4)       p(6*(i-1)+3,4)=b(nn(i,8),1)       !the 2,3-case       p(6*(i-1)+4,1)=b(i,2)       p(6*(i-1)+4,2)=b(i,3)       p(6*(i-1)+4,3)=b(nn(i,4),3)       p(6*(i-1)+4,4)=b(nn(i,6),2)       !the 2,4-case       p(6*(i-1)+5,1)=b(i,2)       p(6*(i-1)+5,2)=b(i,4)       p(6*(i-1)+5,3)=b(nn(i,4),4)       p(6*(i-1)+5,4)=b(nn(i,8),2)       !the 3,4-case       p(6*(i-1)+6,1)=b(i,3)       p(6*(i-1)+6,2)=b(i,4)       p(6*(i-1)+6,3)=b(nn(i,6),4)       p(6*(i-1)+6,4)=b(nn(i,8),3)    ENDDO    RETURN  END SUBROUTINE create_nn_4d  !---------------------------------------------------END MODULE utility!---------------------------------------------------PROGRAM plaquettes  USE random_mod  USE utility  IMPLICIT NONE  INTEGER :: nx,ny,nz,nw,flag, counter  INTEGER, DIMENSION(:,:), ALLOCATABLE :: nn, b, p  INTEGER, DIMENSION(:), ALLOCATABLE :: bval  INTEGER :: i, j, i1, j1, c, c2, n, k  REAL*8  :: rand, beta, thbeta, chbeta, shbeta, E, Esum  REAL*8  :: propup  INTEGER :: d, v, dataf=3   INTEGER :: u1=1, u2=2, condor, mc, mc2  REAL :: iterreal, thermreal  INTEGER :: iter, therm, nbnew, nbold, randpl, bcount  CHARACTER*40 :: data_format, data_file  INTEGER :: s, n2, c1  REAL*8 ::  time2, time1  !needs file 'condor.dat' in same directory  !condor.dat consists of 1 line only  !line = either 0 or 1   OPEN(unit=u1,file='condor.dat')  READ(u1,*) condor;  CLOSE(unit=u1)  IF (condor.EQ.0) THEN     OPEN(unit=u1,file='input_pl.dat')     READ(u1,*) nx;      READ(u1,*) d;     READ(u1,*) iterreal ! iterations for MC;      READ(u1,*) thermreal ! number of thermalization sweeps     READ(u1,*) beta;      CLOSE(unit=u1)  ENDIF  IF (condor.EQ.1) THEN     READ(5,*) nx;      READ(5,*) d;     READ(5,*) iterreal ! iterations for MC;      READ(5,*) thermreal ! number of thermalization sweeps     READ(5,*) beta;   ENDIF  iter = INT(iterreal)   !necessary for reading 1E3-style input  therm = INT(thermreal)  WRITE(*,*) nx   WRITE(*,*) d  WRITE(*,*) iter    WRITE(*,*) therm  WRITE(*,*) beta  WRITE(*,*) '-----------'  IF (MOD(nx,2).EQ.1) THEN     WRITE(*,*) 'nx has to be even. Exiting.'     STOP  ENDIF  IF (d.EQ.2) THEN     ny = nx;     v = nx*ny;  ELSEIF (d.EQ.3) THEN     ny = nx; nz = nx;     v = nx*ny*nz;  ELSEIF (d.EQ.4) THEN     ny = nx; nz = nx; nw = nx;     v = nx*ny*nz*nw  ENDIF  ALLOCATE(nn(1:v,1:2*d))  ALLOCATE(b(1:v,1:d),p(1:d*(d-1)/2*v,1:4))  ALLOCATE(bval(1:d*v))  !create the filename for datafile  IF (nx .LT. 1000) THEN     data_format = '(a11,i3,a1,i1,a2,f6.4,a4)'  ENDIF  IF (nx .LT. 100) THEN     data_format = '(a11,i2,a1,i1,a2,f6.4,a4)'  ENDIF  IF (nx .LT. 10) THEN     data_format = '(a11,i1,a1,i1,a2,f6.4,a4)'  ENDIF  WRITE(data_file,data_format) 'dataseries_', nx, '_', d, 'd_', beta, '.dat'  OPEN(unit=dataf,file=data_file)  !label the bonds linearly  c = 0;  DO i = 1,v     DO j = 1,d        c = c+1;        b(i,j) = c;     ENDDO  ENDDO  IF (d.EQ.2) THEN     CALL create_nn_2d(nx,ny,v,b,nn,p)  ELSEIF (d.EQ.3) THEN     CALL create_nn_3d(nx,ny,nz,v,b,nn,p)  ELSEIF (d.EQ.4) THEN     CALL create_nn_4d(nx,ny,nz,nw,v,b,nn,p)  ENDIF  !b=-1 means bond set, b=+1 means no bond set  bval = +1     !initially no bonds are set, cold start  IF (ALL(bval.EQ.+1)) THEN     nbnew = 0     nbold = 0  ELSE                nbnew = v*d  !only valid, if all bonds = -1     nbold = v*d  END IF  thbeta = TANH(beta);  shbeta = SINH(beta);  chbeta = COSH(beta);  k = 0  E = 0; Esum = 0;  mc = 0; mc2 = 0;  n = 0;  c1 = 1;  CALL cpu_time(time1) !start timer  DO s = 1,iter+therm !start main MC loop     DO n2 = 1,v        n = n + 1;        !select a random plaquette        randpl = INT(d*(d-1)/2*v*random()) + 1          bcount = 0        !count the number of bonds set at random plaqu.        DO j = 1,4           IF (bval(p(randpl,j)).EQ.-1) THEN              bcount = bcount+1           ENDIF        ENDDO        nbnew = nbold + 4 - 2*bcount        !calculate update propability, Metropolis        propup = MIN(1,EXP((nbnew-nbold)*LOG(thbeta)))        rand = random()            IF (rand.LE.propup) THEN   !Move accepted           mc = mc+1;           DO j = 1,4              bval(p(randpl,j))=bval(p(randpl,j))*(-1)           ENDDO           nbold = nbnew;        ELSE                 !Move not accepted           nbnew = nbold           mc2 = mc2+1;        ENDIF        IF (s.GT.therm) THEN           E = -d*thbeta - 1.0/v * nbnew * 1.0/(shbeta*chbeta)  !v-normed energy           Esum=Esum+E           k=k+1        ENDIF     ENDDO ! end main MC loop     IF (s.GT.therm) THEN             WRITE(dataf,*) E !only write every v-th value to file     ENDIF     ! print how many iterations done, how long sill to go     IF (s .EQ. (iter+therm)/10*c1) THEN        IF (c1.EQ.1) CALL cpu_time(time2)        WRITE(*,'(i5,a6,i5,a8,f10.2,a12)') s, ' done,', iter+therm-s, ' to do, ', (time2-time1)*(10-c1), ' sec. to go.'        c1=c1+1     ENDIF  ENDDO  WRITE(*,*) Esum/k  WRITE(*,*) 'Number of accepted Moves: ', mc  WRITE(*,*) '# NOT accepted: ', mc2  WRITE(*,*) 'Sum: ', mc+mc2  CLOSE(unit=dataf)END PROGRAM plaquettes