In [ ]:

    

##### Dependências e config iniciais #####

#%matplotlib inline #Comente/remova essa linha e reinicie o IPython/Jupyter caso queira plots em janelas separadas

import numpy as np
import matplotlib.pyplot as plt
import time


    

In [ ]:

    

##### Bloco de condições iniciais e de simulação #####

#Constante de simulação
J = 1 #Valor númerico do spin
h = 0 #Constante associada a cada sítio de spin

#Parâmetros da simulação
N = 1000   #Iterações de Metropolis a serem feitas
Nmtr = 10    #Tamanho da matriz de spins inicial

T = np.linspace(0, 5, 16)    #Temperaturas a serem iteradas. (~10s por iteração em meu PC)
#Multiplique pela constante de Boltzmann para obter a grandeza em Kelvins

#Matriz para simulação
sigma_0 = np.ones((Nmtr, Nmtr))    #Matriz com os spins iniciais


    

In [ ]:

    

##### Bloco de funções #####

###
#Obter a energia total de uma matriz de spins - extrapolação do método usado na sala
#não utilizado devido a ineficiência
def energiaIsing2(spin_grid):
    shp = spin_grid.shape
    E = 0
    if len(shp) == 1:
        line_len = spin_grid.shape[0]
        for i in range(0, line_len-1):
            E -= J * spin_grid[i] * spin_grid[i+1] + h * spin_grid[i]
        E += -(J * spin_grid[line_len-1] * spin_grid[0] + h * spin_grid[line_len-1])
    else:
        for r in spin_grid:
            E += energiaIsing(r)
            E += energiaIsing(np.transpose(r))
    return E    

###
#Obter a energia total de uma matriz de spins - método matricial mais eficiente
def energiaIsing(spin_grid):
    #Obter uma matriz com a soma de todos os vizinhos circulares em cada sítio da matriz 
    neighbors = (np.roll(spin_grid, 1, axis=0) + np.roll(spin_grid, -1, axis=0) +
    np.roll(spin_grid, 1, axis=1) + np.roll(spin_grid, -1, axis=1))
    
    #Retornar a soma de todos os elementos da matriz anterior junto com o H associado
    #Note que a divisão por dois da matriz de vizinhos é devido a cada ligação ser associada duas vezes
    return -(J * np.sum(np.multiply(spin_grid, neighbors)) / 2 + np.sum(h * spin_grid))

###
#Fazer iteração do algoritimo de Metropolis para uma rede 2D de spins
def iterMetropolis(sigma_0, t):
    E = energiaIsing(sigma_0)    #Energia inicial
    m = np.zeros(N)    #Lista contendo a magnetização média em cada iteração
    e = np.zeros(N)    #Lista contento a energia média em cada iteração
    cfg = np.zeros(N)
    t1 = time.time()
    
    #Inicializar a matriz para uso nas iterações
    sigma = np.copy(sigma_0)
    
    for n in range(0, N):
        for i in range(0, Nmtr):
            for j in range(0, Nmtr):
                #Inicializar uma matriz temporária para a flipagem do spin
                sigma_temp = np.copy(sigma)        
                sigma_temp[i, j] = -sigma[i,j]
                
                E1 = energiaIsing(sigma_temp)
                Eflip = E1 - E
                
                if Eflip <= 0:
                    E = E1
                    sigma = np.copy(sigma_temp)
                else:
                    Pflip = np.exp(-Eflip / t)
                    r = np.random.rand(1)
                    
                    if r < Pflip:
                        E = E1
                        sigma = np.copy(sigma_temp)
                        
        #Guardar a magnetização e energia média da iteração                
        m[n] = np.sum(sigma) / sigma.size
        e[n] = E / sigma.size
        cfg[n] = n
    
    t2 = time.time()
    print("T/k_b: %.2f, tempo de execução: %.2f s" % (t, t2-t1))
    output = {'Em': np.average(e), 'M': np.average(m), 'e':e, 'm':m, 'cfg':cfg}
    return output


    

In [ ]:

    

##### Bloco de simulações #####
Em = []
M = []

for t in T:    
    ite = iterMetropolis(sigma_0, t)
    Em.append(ite['Em'])
    M.append(ite['M'])

Em = np.array(Em)
M = np.array(M)


    

In [ ]:

    

##### Plotagem da magnetização média #####

plt.plot(T, M, '-*', label='Magnetização média')
plt.plot(T, np.zeros(len(T)), '-')
plt.ylim((-1.1*J, 1.1*J))
plt.xlim((0, np.max(T)))
plt.legend(loc=0)
plt.xlabel("Temperatura ($K/k_b$)")
plt.ylabel("<m>")
plt.title("Magnetização média por temperatura \n h=%s, %s x %s spins, %s iter. de Metropolis" % (h, Nmtr, Nmtr, N))
plt.show()


    

In [ ]:

    

##### Plotagem da energia média #####

plt.plot(T, Em, '-*', label='Energia total do sistema')
plt.plot(T, np.zeros(len(T)), '-')
plt.legend(loc=0)
plt.xlim((0, np.max(T)))
plt.xlabel("Temperatura [$K/k_b$]")
plt.ylabel("<E>")
plt.title("Energia média por temperatura \n h=%s, %s x %s spins, %s iter. de Metropolis" % (h, Nmtr, Nmtr, N))
plt.show()