Metropolis - Potts 2D

  MODULE variaveis_globais
      REAL(8) :: pi = 4*atan(1.)
      INTEGER :: L,L2,MCS,o,v,Q_state
      REAL :: Temp,B
      INTEGER,DIMENSION(33) :: seed
      REAL,DIMENSION(10) :: temperatura
  END MODULE variaveis_globais

  PROGRAM potts
      USE variaveis_globais
      IMPLICIT NONE
      DO o=1,33
          seed(o) =o 
      END DO
      CALL random_seed(put=seed)

      B = 1 
      L = 64 
      L2 = L*L
      MCS = 2000
      Temp = 1.
      Q_state = 100 
      
   !   DO v=1,10
   !       temperatura(v) = 0.2*v 
   !   END DO

      CALL metropolis
      
  END PROGRAM 
 !-------------------------------------------------------------!
 !                   Algoritmo de Metropolis                   !
 !-------------------------------------------------------------!
  SUBROUTINE metropolis
      USE variaveis_globais
      IMPLICIT NONE
      INTEGER :: i,j,t,s,r,r2,r_new,f,n
      INTEGER :: k_right,k_left,k_up,k_down
      REAL(8) :: de,rn,prob,E0,E1,E2
      INTEGER,DIMENSION(L2) :: right,left,up,down
      REAL(8),DIMENSION(L2) :: energia
      REAL,DIMENSION(L2) :: spin
      REAL(8),DIMENSION(MCS) :: e_media,mag
      REAL(8),DIMENSION(:),ALLOCATABLE :: theta
      REAL(8),DIMENSION(10) :: energia_media
       
      DO i=1,L
          DO j=1,L
              right((i-1)*L+j) = (i-1)*L+j+1
              left((i-1)*L+j) = (i-1)*L+j-1
              up((i-1)*L+j) = (i-1)*L+j-L
              down((i-1)*L+j) = (i-1)*L+j+L
              IF (i == 1) THEN
                  up((i-1)*L+j) = L2-L+j
              END IF
              IF (i == L) THEN
                  down((i-1)*L+j) = j
              END IF
              IF (j == 1) THEN
                  left((i-1)*L+j) = i*L
              END IF
              IF (j == L) THEN
                  right((i-1)*L+j) = i*L-(L-1)
              END IF
          END DO
      END DO

      spin = spins_aleat(Q_state)
      ALLOCATE(theta(Q_state))
   !   WRITE(*,*) spin
  
      E0 = 0
      DO s=1,L2
          k_right =  kronecker(spin(i),spin(right(i))) 
          k_left =  kronecker(spin(i),spin(left(i))) 
          k_up =  kronecker(spin(i),spin(up(i))) 
          k_down =  kronecker(spin(i),spin(down(i))) 
          E0 = E0 - B*(k_right + k_left + k_up + k_down)
      END DO
      WRITE(*,*) E0
      
      DO n=1,Q_state
          theta(n) = 2.*pi*n/Q_state
      END DO

      OPEN(2,file='energia.dat',status='unknown')
     ! DO f=1,10 
        DO t=1,MCS
          DO s=1,L2
              r = u(1,L2)
2             r_new = u(1,L2)
              IF (spin(r_new) == spin(r)) GO TO 2
              IF (spin(r_new) /= spin(r)) GO TO 3
              
              spin(r_new) = theta(u(1,Q_state))
3             k_right = kronecker(spin(r),spin(right(r))) 
              k_left = kronecker(spin(r),spin(left(r))) 
              k_up = kronecker(spin(r),spin(up(r))) 
              k_down = kronecker(spin(r),spin(down(r))) 
              E1 = B*(k_right + k_left + k_up + k_down)
              
              k_right = kronecker(spin(r_new),spin(right(r))) 
              k_left = kronecker(spin(r_new),spin(left(r))) 
              k_up = kronecker(spin(r_new),spin(up(r))) 
              k_down = kronecker(spin(r_new),spin(down(r))) 
              E2 = B*(k_right + k_left + k_up + k_down)
              
              de = E2 - E1 
              IF (de .le. 0) THEN
                  spin(r) = spin(r_new)
              ELSE
                  prob = exp(-de/Temp)
                  CALL random_number(rn)
                  IF (rn < prob) THEN
                      spin(r) = spin(r_new)
                  END IF
              END IF
              energia(s) = E0 + de
          END DO   
          e_media(t) = 1./L2*sum(energia)
          mag(t) = 1./L2*sum(spin)
        WRITE(2,*) t,e_media(t),mag(t)
        END DO
     !   energia_media(f) = 1./MCS*sum(e_media)
      !  WRITE(2,*) temperatura(f),energia_media(f)
     ! END DO
      CLOSE(2)

      CONTAINS 
!-------------------------------------------------------------!
!                          Funções                            !
!-------------------------------------------------------------!

!Inicializa os spins de forma aleatória
      FUNCTION spins_aleat(Q)
          USE variaveis_globais
          INTEGER :: n,k,Q,m
          REAL,DIMENSION(Q) :: theta
          REAL,DIMENSION(L2) :: spins_aleat
          DO n=1,Q
              theta(n) = 2.*pi*n/Q
          END DO
          IF (Q .le. 2) THEN
              DO k=1,L2
                  spins_aleat(k) = cos(theta(u(1,Q)))
              END DO 
          ELSE IF (Q > 2) THEN
              DO m=1,L2
                  spins_aleat(m) = theta(u(1,Q))
              END DO
          END IF
          RETURN
      END FUNCTION spins_aleat

!Inicializa os spins d forma ordenada, todos iguais
      FUNCTION spins(Q)
          USE variaveis_globais
          INTEGER :: n,k,Q,m
          REAL,DIMENSION(Q) :: theta
          REAL,DIMENSION(L2) :: spins
          DO n=1,Q
              theta(n) = 2.*pi*n/Q
          END DO
          IF (Q .le. 2) THEN
              DO k=1,L2
                  spins(k) = cos(theta(1))
              END DO 
          ELSE IF (Q > 2) THEN
              DO m=1,L2
                  spins(m) = theta(1)
              END DO
          END IF
          RETURN
      END FUNCTION spins

!Delta de Kronecker       
      FUNCTION kronecker(i,j)
          REAL :: i,j
          INTEGER :: kronecker
          IF (i == j) THEN
              kronecker = 1
          ELSE
              kronecker = 0
          END IF
          RETURN
      END FUNCTION kronecker

!Gera números inteiros aleatórios
      FUNCTION u(a,b)
          INTEGER :: u,a,b
          REAL :: r
          CALL random_number(r)
          u = int(b*r) + a
          RETURN
      END FUNCTION u
!-------------------------------------------------------------!
      END SUBROUTINE metropolis