Metropolis - Potts 2D: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
Sem resumo de edição
Sem resumo de edição
 
(2 revisões intermediárias pelo mesmo usuário não estão sendo mostradas)
Linha 22: Linha 22:
       MCS = 2000
       MCS = 2000
       Temp = 1.
       Temp = 1.
       Q_state = 100
       Q_state = 7
        
        
!      DO v=1,10
!      DO v=1,10

Edição atual tal como às 19h52min de 18 de outubro de 2022

  
   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