Exemplo de código para a amostragem de Wang-Landau em Python3

De Física Computacional
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()