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()