C     MODELO DE ISING EM UM CAMPO MAGNÉTICO NULO
C
C     Simulação de Monte Carlo, com o algoritmo de Metropolis, para o modelo 
C     de Ising na rede quadrada contendo L por L spins.
C
C     A simulação utiliza MCDIS passos de Monte Carlo por spin na termalização
C     e MCSTEP passos de Monte Carlo por spin no cálculo das médias. 
C     O estado de cada spin é guardado no vetor IS, de L*L componentes, cada
C     uma das quais pode assumir os valores -1 (spin para baixo) ou +1 (spin 
C     para cima).
C     Utilizamos condições de contorno periódicas, impostas com auxílio de 
C     uma matriz de primeiros vizinhos NN, de ordem 4x(L*L).
C     A simulação percorre temperaturas entre T0 e TMAX com incremento TINC,
C     calculando a energia, a magnetização, a suscetibilidade e o
C     cumulante de Binder.
C
C     Declarações de variáveis
      IMPLICIT NONE ! Força a declaração de todas as variáveis
      INTEGER L,L2,MCDIS,MCSTEP,IS,NN
      INTEGER I,I0,I1,I2,I3,I4,IH,IAUX,MC,MAG,IU,ISO
      PARAMETER (L=10,L2=L*L)
      PARAMETER (MCDIS=1000000,MCSTEP=500000) 
      DIMENSION IS(L2),NN(4,L2)
      DOUBLE PRECISION T,PV,R,DE,RMAG,AV,AV2,AV3,AV4,T0,TINC,TMAX
      DOUBLE PRECISION PB(-1:1,-4:4)
      CHARACTER CL*2
C
C     Construindo a matriz de primeiros vizinhos
      DO I=1,L2
         NN(1,I)=I+1
         IF (MOD(I,L).EQ.0) NN(1,I)=NN(1,I)-L
         NN(2,I)=I+L
         IF (NN(2,I).GT.L2) NN(2,I)=NN(2,I)-L2
         NN(3,I)=I-1
         IF (MOD(I,L).EQ.1) NN(3,I)=NN(3,I)+L
         NN(4,I)=I-L
         IF (NN(4,I).LE.0) NN(4,I)=NN(4,I)+L2
      END DO
C     Parâmetros relacionados às temperaturas
      T0=2.2 
      TMAX=2.4
      TINC=0.005
C
C     Criamos um arquivo de dados para salvar os resultados, com a informação
C     sobre o tamanho do sistema
!      CALL CHARINT (L,2,CL)                    
!     OPEN(9,file='ising_testeL'//CL//'.dat')
      OPEN(9,file='ising_testeL10_f.dat') 
C     
      T=T0                 ! Começamos da temperatura mínima
      DO WHILE (T.LT.TMAX) ! Varremos os valores de temperatura
         PRINT*,T
C        Zerando as variáveis que armazenam as grandezas de interesse
         AV=0D0
         AV2=0D0
         AV3=0D0
         AV4=0D0
C        Calculando a matriz de pesos de Bolzmann para a temperatura T.
C        Isso evita o cálculo da exponencial a cada tentativa de inversão de 
C        um spin, acelerando bastante a simulação.
         DO I0=-1,1,2
            DO IAUX=-4,4,2
C              Variação da energia ao invertermos um spin do estado I0 para
C              -I0 quando a soma dos estados de seus 4 vizinhos é IAUX
               DE=DBLE(2*I0*IAUX)
C              Fator de Bolzmann correspondente. Note que variações negativas
C              de energia dão origem a fatores maiores do que 1.
               PB(I0,IAUX)=EXP(-DE/T)
            END DO
         END DO
C        Na temperatura mais baixa, impomos que os spins apontam inicialmente
C        todos para cima. Em outras temperaturas, mantemos a configuração do
C        final da simulação na temperatura anterior.
         IF (T.EQ.T0) THEN
            DO I=1,L2    
               IS(I)=1
            END DO       
         END IF
C        Loop da termalização
         DO 12 MC=1,MCDIS 
C           Damos a cada spin, sequencialmente, a chance de ser invertido. 
            DO IH=1,L2
               ISO=IS(IH)       ! Estado do spin candidato à inversão
               I1=IS(NN(1,IH))  ! Estado do vizinho à direita
               I2=IS(NN(2,IH))  ! Estado do vizinho abaixo
               I3=IS(NN(3,IH))  ! Estado do vizinho à esquerda
               I4=IS(NN(4,IH))  ! Estado do vizinho acima
               IAUX=I1+I2+I3+I4 
               PV=PB(ISO,IAUX)  ! Fator de Boltzmann PV associado à inversão
               R=1D0
               IF (PV.LT.R) R=RAND()  ! Geramos um número aleatório se PV<1
C              A inversão do spin ocorre de fato se o fator de Boltzmann for
C              maior que o número aleatório gerado. Isso sempre ocorrerá se
C              o fator for maior que 1 (variação negativa de energia).
               IF (PV.GT.R) IS(IH)=-ISO
            END DO
12       CONTINUE ! Fim do loop de termalização
C        Loop de medição
         DO 13 MC=1,MCSTEP
            DO IH=1,L2
               ISO=IS(IH)      
               I1=IS(NN(1,IH))
               I2=IS(NN(2,IH))
               I3=IS(NN(3,IH))
               I4=IS(NN(4,IH))
               IAUX=I1+I2+I3+I4
               PV=PB(ISO,IAUX)
               R=1D0
               IF (PV.LT.R) R=RAND()  ! Geramos um número aleatório se PV<1
               IF (PV.GT.R) IS(IH)=-ISO
            END DO
C           Ao final de cada varredura da rede, tomamos as medidas
            MAG=0
            IU=0
            DO I=1,L2 ! Somamos as magnetizações e as energias de interação
               MAG=MAG+IS(I)
               IU=IU-IS(I)*(IS(NN(1,I))+IS(NN(2,I)))
            END DO
            RMAG=ABS(DBLE(MAG))/L2  ! Módulo da magnetização por spin
            AV=AV+RMAG              ! Acumulamos a magnetização
            AV2=AV2+DBLE(IU)/L2     ! Acumulamos a energia por spin
            AV3=AV3+RMAG**2         ! Acumulamos a magnetização ao quadrado
            AV4=AV4+RMAG**4         ! Acumulamos a magnetização à quarta
13       CONTINUE ! Fim do loop de medição
C        Cálculo das médias
         AV=AV/MCSTEP
         AV2=AV2/MCSTEP
         AV3=AV3/MCSTEP
         AV4=AV4/MCSTEP
         AV4=1D0-AV4/(3*AV3*AV3) ! Cálculo do cumulante de Binder
         AV3=L2*(AV3-AV*AV)/T    ! Cálculo da suscetibilidade magnética
         WRITE(9,99) T,AV,AV3,AV2,AV4  ! Salvando os resultados no arquivo
         ! 1a coluna guarda a temperatura
         ! 2a coluna guarda a magnetização
         ! 3a coluna guarda a suscetibilidade magnética
         ! 4a coluna guarda a energia média
         ! 5a coluna guarda o cumulante de Binder
C        Incrementamos a temperatura
         T=T+TINC
      END DO ! Fim do loop das temperaturas
      CLOSE(9)
99    FORMAT (10F16.6)
      END
C
C
C
      SUBROUTINE CHARINT (i,n,ci) 
C     Produz uma variável de caractere 'ci' a partir de um número inteiro 'i',
C     registrando 'n' algarismos
      INTEGER i,n
      CHARACTER*10 ci
C
      DO j=n,1,-1
         ci((n+1-j):(n+1-j))=CHAR(48+MOD(i,10**j)/10**(j-1))
      END DO
      RETURN
      END