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 = 7 
       
!      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