plaquettes.f90

用fortran编写的程序.主要用于研究统计物理里的Ising模型.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Fortran90 program, Monte Carlo simulation of the Ising model
!programmed by Hans-Marc Erkinger, [e-mail]
!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 iseed
CONTAINS
  !---------------------------------------------------
  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 random
END 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 NONE
CONTAINS
  !---------------------------------------------------
  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