Exemplo de código para a amostragem de Wang-Landau em Python3
Ir para navegação
Ir para pesquisar
O código foi feito para execução em notebook, por isso as bibliotecas foram importadas apenas uma vez e os código está dividido em seções.
Amostragem Wang-Landau
import numpy as np
import matplotlib.pyplot as plt
import random as R
from numba import jit # Usamos o numba para fazer cálculos mais rápidos
# Variáveis inicias
J = -1
T = 2.3
beta = 1/T
L = 16
N = L**2
Qmax = 2 # 2 = [0,1]
grid = np.arange(0,L**2).reshape((L,L))
# Gerador de uma grade LxL com spins aleatórios
def gentable(grid): # Gera uma tabela com os vizinhos de todas as posições
n = len(grid)
table = np.empty((0,4), int)
for i in range(0,n):
for j in range(0,n):
table = np.vstack([table,neighbor(grid,n,i,j)])
return table
def neighbor(grid,n,row,column): # Devolve os vizinhos do item i,j dado uma grid
return [grid[row][column+1] if column<n-1 else grid[row][0],
grid[row-1][column] if row>0 else grid[n-1][column],
grid[row][column-1] if column>0 else grid[row][n-1],
grid[row+1][column] if row<n-1 else grid[0][column]]
@jit
def Round(N,J,table,states,He,lnGe,ene,newEne,logf):
for i in range(10000): # Fazemos os 10000 MC step
for j in range(N): # MC step
selectedSpin = R.randint(0,N-1) # Selecionamos um spin para ser flipado
delE = -2*J*states[selectedSpin]*sum(states[table[selectedSpin]]) # Calulamos a variação de energia
newEne = ene + delE # Soma dos estados ao redor do selecionado, se são todos -1, e o spin selecionado é 1 e o resultado é -4 * 2 * J
histogramIndex = int(ene/4 + N/2) # Descobrimos qual o index de H(E) que o sistema está
newHistogramIndex = int(newEne/4 + N/2) # Descobrimos qual o index de H(E) que o sistema quer ir
if(R.random() < np.exp(lnGe[histogramIndex]-lnGe[newHistogramIndex])): # Condição para aceitar o spin
# Se aceitar o flip
states[selectedSpin] *= -1 # flipa-se o spin
ene = newEne # atualiza-se a energia
lnGe[newHistogramIndex] += logf # Aumenta o g(E) da nova energia pelo parâmetro f
He[newHistogramIndex] += 1 # Aumenta o H(E) da nova energia
else:
# Se rejeitar o flip
lnGe[histogramIndex] += logf # Aumenta o g(E) da mesma energia
He[histogramIndex] += 1 # Aumenta o H(E) da mesma energia
return states,ene,He,lnGe
def sim(N,J,T,table,states):
lnGe = np.zeros(N+1) # Inicializamos g(E)
ene = 0
newEne = 0
f = np.e # Inicializamos f como e
logf = np.log(f)
for i in range(len(states)): # Calcuclamos a energia pela primeira vez
ene += states[i]*states[table[i][0]]+states[i]*states[table[i][1]]
ene *= J
while(f - 1 > 1e-6): # Enquanto o parâmetro F nao chega ao limite de precisão desejado
He = np.zeros(N+1) # Se inicia ou reinicia os histogramas as serem "planificados"
He[int(ene/4 + N/2)] += 1 # Adiciona no H(E) a energia inicial
while(True): # condição de parada no (*), feito desta maneira pois queremos que o código rode pelo menos uma vez antes da condição de parada.
states,ene,He,lnGe = Round(N,J,table,states,He,lnGe,ene,newEne,logf) #Um 'round' de 10000 MC steps
HeMask = He[He != 0] #Necessário por causa das energias proibidas em -2N + 8 e 2N - 8
if(not np.any(HeMask/np.mean(HeMask) < 0.8)): # (*) Vemos se todos as medidas de H(E) estão por volta de 80% da média.
#se estiver por volta de 80%, saimos do while, se não vamos pra outro round.
break
f = np.sqrt(f) # Após a planificação do histograma, tiramos a raiz do parâmetro f.
logf = np.log(f)
FinalGe = np.exp(lnGe - lnGe[0] + np.log(2)) # Acabada a simulação, normalizamos a g(E).
FinalGe[1] = 1 # O algoritmo não da conta que no modelo de Ising estas energias não tem densidade
FinalGe[N-1] = 1
return FinalGe
table = gentable(grid)
states = np.random.choice([1,-1], size=N)
Ge = sim(N,J,T,table,states) # Inicia a simulação e recebe g(E)
Plot Distribuição Canônica
#xlin e modificações no GePlot são feitas para a função ter 0 e a curva ser fechada no plot
xlin = np.arange(-2-4/N,2+5/N,4/N)
GePlot = np.insert(Ge, 0, 1, axis=0)
GePlot = np.insert(GePlot, len(Ge), 1, axis=0)
plt.plot(xlin,np.log(GePlot), color='black')
plt.title('Densidade de Estados')
plt.ylabel(r'$\ln(g(e))$')
plt.xlabel(r'$E/N$')
plt.axhline(color = 'gray', linewidth = 0.7)
plt.axvline(color = 'gray', linewidth = 0.7)
plt.show()
Plot Densidade de Estados
T = 2.3
T2 = 2.4
T3 = 2.2
beta = 1/T
beta2 = 1/T2
beta3 = 1/T3
lnGe = np.log(Ge)
xlin = np.arange(-2,2+1/N,4/N)
E = np.arange(-2*N,2*N+1,4)
E[1] = 0 #Não há densidade com essa energia
E[N-1] = 0
Pe = Ge*np.exp(-E*beta)
Pe2 = Ge*np.exp(-E*beta2)
Pe3 = Ge*np.exp(-E*beta3)
plt.plot(xlin,Pe/np.max(Pe),label = r"$K_bT_c \approx 2.3$")
plt.plot(xlin,Pe2/np.max(Pe2),label = r"$K_bT = 2.4$")
plt.plot(xlin,Pe3/np.max(Pe3),label = r"$K_bT = 2.2$")
plt.legend()
plt.title('Distribuição Canônica')
plt.xlim((-2,0))
plt.ylabel(r'$g(E)e^{E/K_bT_c}$')
plt.xlabel(r'$E/N$')
plt.axhline(color = 'gray', linewidth = 0.7)
plt.axvline(color = 'gray', linewidth = 0.7)
plt.show()
Plot Energia Interna
Tlin = np.arange(0.05,8,0.01, np.float128) # KbT de 0.05 até 8.
betaLin = 1/Tlin
beta = 1/0.8
E = np.arange(-2*N,2*N+1,4, np.float128)
E[1] = 0 #Não há densidade com essa energia
E[N-1] = 0
def U(beta,E,Ge):
return np.sum(E*Ge*np.exp(-E*beta))/np.sum(Ge*np.exp(-E*beta))
Ulin = np.zeros(len(Tlin))
for i in range(len(betaLin)):
Ulin[i] = U(betaLin[i],E,Ge)
plt.plot(Tlin,Ulin/N, color='black')
plt.title('U - Energia Interna')
plt.xlim((0,8))
plt.ylabel(r'$U(T)/N$')
plt.xlabel(r'$k_BT$')
plt.axhline(color = 'gray', linewidth = 0.7)
plt.axvline(color = 'gray', linewidth = 0.7)
plt.show()
Plot Calor Específico
Tlin = np.arange(0.05,8,0.01, np.float128) # KbT de 0.05 até 8.
betaLin = 1/Tlin
beta = 1/0.8
E = np.arange(-2*N,2*N+1,4, np.float128)
E[1] = 0 #Não há densidade com essa energia
E[N-1] = 0
def E2(beta,E,Ge):
return np.sum(np.power(E,2)*Ge*np.exp(-E*beta))/np.sum(Ge*np.exp(-E*beta)) #Calcula o segundo momento da distribuição de energia
def E1(beta,E,Ge):
return np.sum(E*Ge*np.exp(-E*beta))/np.sum(Ge*np.exp(-E*beta)) #Calcula o primeiro momento momento da distribuição de energia
E2lin = np.zeros(len(Tlin))
E1lin = np.zeros(len(Tlin))
for i in range(len(betaLin)):
E2lin[i] = E2(betaLin[i],E,Ge)
E1lin[i] = E1(betaLin[i],E,Ge)
Clin = (E2lin - E1lin**2)/Tlin**2
plt.plot(Tlin,Clin/N, color='black')
plt.title('C - Calor Específico')
plt.xlim((0,8))
plt.ylabel(r'$C(T)/N$')
plt.xlabel(r'$k_BT$')
plt.axhline(color = 'gray', linewidth = 0.7)
plt.axvline(color = 'gray', linewidth = 0.7)
plt.show()
Plot Energia Livre de Helmholtz
Tlin = np.arange(0.05,8,0.01, np.float128) # KbT de 0.05 até 8.
betaLin = 1/Tlin
beta = 1/0.8
E = np.arange(-2*N,2*N+1,4, np.float128)
E[1] = 0 #Não há densidade com essa energia
E[N-1] = 0
def Z(beta,E,Ge):
return np.sum(Ge*np.exp(-E*beta)) # Calcula a função de partição para a temperatura T
def F(beta,E,Ge):
return -np.log(Z(beta,E,Ge))/beta#Calcula o primeiro momento momento da distribuição de energia
Flin = np.zeros(len(Tlin))
for i in range(len(betaLin)):
Flin[i] = F(betaLin[i],E,Ge)
plt.plot(Tlin,Flin/N, color='black')
plt.title('F - Energia Livre de Helmholtz')
plt.xlim((0,8))
plt.ylim((-6,-1.8))
plt.ylabel(r'$F(T)/N$')
plt.xlabel(r'$k_BT$')
plt.axhline(color = 'gray', linewidth = 0.7)
plt.axvline(color = 'gray', linewidth = 0.7)
plt.show()
Plot Entropia
Tlin = np.arange(0.05,8,0.01, np.float128) # KbT de 0.05 até 8.
betaLin = 1/Tlin
beta = 1/0.8
E = np.arange(-2*N,2*N+1,4, np.float128)
E[1] = 0 #Não há densidade com essa energia
E[N-1] = 0
Slin = (Ulin - Flin)/Tlin
plt.plot(Tlin,Slin/N, color='black')
plt.title('S - Entropia', color='black')
plt.xlim((0,8))
plt.ylabel(r'$F(T)/N$')
plt.xlabel(r'$k_BT$')
plt.axhline(color = 'gray', linewidth = 0.7)
plt.axvline(color = 'gray', linewidth = 0.7)
plt.show()