Amostragem de Wang-Landau: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
Sem resumo de edição
 
(52 revisões intermediárias pelo mesmo usuário não estão sendo mostradas)
Linha 1: Linha 1:
O algorítmo de amostragem de Wang-Landau é um método de amostragem para simulações de Monte Carlo introduzido por F.Wang e D.P Landau em 2001 que apresenta diversas vantagens sobre outros métodos para sistemas de spins. Dentre eles podemos citar o algoritmo de Metropolis, o algoritmo de ''clustering'' de Wolff, ou em um modelamento de ''ensamble'' multi-canônico. Nestes dois últimos métodos são utilizados métodos de repesagem de histogramas que são limitados em sistemas grandes devido as baixa qualidade estatística nas asas do histograma.  
O [https://pt.wikipedia.org/wiki/Algoritmo algoritmo] de amostragem de Wang-Landau é um método de amostragem para simulações de [https://pt.wikipedia.org/wiki/Monte_Carlo Monte Carlo] introduzido por F. Wang and D. P. Landau em 2001
Dentro deste contexto, o algoritmo Wang-Landau promete resolver problemas encontrados em amostragens convencionais como o ''critical slowing down'' para temperaturas próximas da temperatura critica <math>T_c</math> utilizando-se de caminhadas aleatórias controladas para mapear a densidade de estados <math>g(E)</math> de um sistema, sem fazer uso de qualquer repesagem de histogramas.
<ref>F. Wang and D. P. Landau, ‘‘Efficient, multiple-range random walk algorithm to calculate the density of states,’’ Phys. Rev. Lett. 86, 2050–2053
!2001"; ‘‘Determining the density of states for classical statistical models:
A random walk algorithm to produce a flat histogram,’’ Phys. Rev. E 64,
056101 !2001".
</ref>
(com uma versão didática publicada em 2004 por D.P Landau, Shan-Ho Tsai e M. Exlerem),
<ref>Landau, D. P., Shan-Ho Tsai, and M. Exler. "A new approach to Monte Carlo simulations in statistical physics: Wang-Landau sampling." American Journal of Physics 72.10 (2004): 1294-1302.</ref>
que apresenta diversas vantagens sobre outros métodos para sistemas de [https://pt.wikipedia.org/wiki/Spin spins]. Dentre eles podemos citar o algoritmo de Metropolis
<ref>N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E.
Teller, ‘‘Equation of state calculations by fast computing machines,’’ J.
Chem. Phys. 21, 1087–1092 !1953".</ref>
, o algoritmo de ''clustering'' de Wolff
<ref>U. Wolff, ‘‘Collective Monte Carlo updating for spin systems,’’ Phys. Rev.
Lett. 62, 361–364 !1989".</ref>
, ou em um modelamento de ''ensamble'' multi-canônico
<ref>B. A. Berg and T. Neuhaus, ‘‘Multicanonical ensemble: A new approach to
simulate first-order phase transitions,’’ Phys. Rev. Lett. 68, 9–12 !1992"</ref>
. Nestes dois últimos métodos são utilizados métodos de repesagem de [https://pt.wikipedia.org/wiki/Histograma histogramas]<ref>
A. M. Ferrenberg and R. H. Swendsen, ‘‘New Monte Carlo technique for
studying phase transitions,’’ Phys. Rev. Lett. 61, 2635–2638 !1988"; ‘‘Optimized Monte Carlo data analysis,’’ ibid. 63, 1195–1198 !1989".</ref>
que são limitados em sistemas grandes devido as baixa qualidade estatística nas asas do histograma.  
Dentro deste contexto, o algoritmo Wang-Landau promete resolver problemas encontrados em amostragens convencionais como o ''critical slowing down'' para temperaturas próximas da temperatura crítica <math>T_c</math> utilizando-se de caminhadas aleatórias controladas para mapear a densidade de estados <math>g(E)</math> de um sistema, sem depender do método de repesagem de histogramas.


A maioria dos algoritmos de amostragem convencionais geram uma distribuição canônica não normalizada
----


<math> P(E,T) = g(E)e^{-E/K_BT}, </math>
A maioria dos algoritmos de amostragem convencionais geram uma [https://pt.wikipedia.org/wiki/Conjunto_can%C3%B3nico distribuição canônica] não normalizada


para uma determinada temperatura $T$, Geralmente estas distribuições são estreitas e se faz necessário múltiplas simulações para obter algum parâmetro termodinâmico para uma distribuição significantemente grande de temperaturas.
<div id="equação de distribuição de probabilidade canônica"> <math> P(E,T) = g(E)e^{-E/K_BT}, </math> </div>
Como <math>g(E)</math>não depende de temperatura do sistema, se pudermos encontrá-lo para todo <math>E</math>, podemos encontrar a função de partição


para uma determinada temperatura <math>T</math>, Geralmente estas distribuições são estreitas e se faz necessário múltiplas simulações para obter algum parâmetro [https://pt.wikipedia.org/wiki/Termodin%C3%A2mica termodinâmico] para uma distribuição significantemente grande de temperaturas, na Fig. 2 vemos exemplos destas distribuições.
Como <math>g(E)</math> não depende de temperatura do sistema, se pudermos encontrá-lo para todo <math>E</math>, podemos encontrar a [https://pt.wikipedia.org/wiki/Fun%C3%A7%C3%A3o_de_parti%C3%A7%C3%A3o_(mec%C3%A2nica_estat%C3%ADstica) função de partição]


<math> Z(T) = \sum_E g(E)e^{-E/K_BT}, </math>  
<div id="pf"> <math> Z(T) = \sum_E g(E)e^{-E/K_BT}, </math> </div>


e o sistema esta essencialmente resolvido, uma vez que grande parte dos parâmetros termodinâmicos podem ser derivados de $Z$. Além disso, a amostragem de Wang-Landau é provada ser útil em diversas aplicações como o antiferromagneto de Potts, sistemas de spins aleatórios, sistemas quânticos, etc... .
e o sistema esta essencialmente resolvido, uma vez que grande parte dos parâmetros termodinâmicos podem ser derivados de <math>Z</math>. Além disso, a amostragem de Wang-Landau é provada ser útil em diversas aplicações como o [https://pt.wikipedia.org/wiki/Antiferromagnetismo antiferromagneto] de Potts
<ref>C. Yamaguchi and Y. Okabe, ‘‘Three-dimensional antiferromagnetic
q-state Potts models: application of the Wang-Landau algorithm,’’ J. Phys.
A 34, 8781– 8794 !2001"</ref>
, sistemas de spins aleatórios<ref>Y. Okabe, Y. Tomita, and C. Yamaguchi, ‘‘Application of new Monte Carlo
algorithms to random spin systems,’’ Comput. Phys. Commun. 146, 63– 68
!2002".</ref>
, sistemas [https://pt.wikipedia.org/wiki/Mec%C3%A2nica_qu%C3%A2ntica quânticos]<ref>
2M. Troyer, S. Wessel, and F. Alet, ‘‘Flat histogram methods for quantum
systems: Algorithms to overcome tunneling problems and calculate the
free energy,’’ Phys. Rev. Lett. 90, 120201 !2003". </ref>, etc... .


==Descrição do algoritmo de Wang Landau==
==Descrição do algoritmo de Wang Landau==
Descreveremos o funcionamento do algoritmo de Wang-Landau num sistema de spins clássicos de 2 estados com valores discretos de energia e sem campo magnético. Portanto quando nos referirmos a <math>g(E)</math> como densidade de estados, interpretamos como o número de estados com energia E.
Descreveremos o funcionamento do algoritmo de Wang-Landau num sistema de spins clássicos de 2 estados numa grade bidimensional com valores discretos de energia e sem [https://pt.wikipedia.org/wiki/Campo_magn%C3%A9tico campo magnético]([https://pt.wikipedia.org/wiki/Modelo_de_Ising Modelo de Ising]2D). Portanto quando nos referirmos a <math>g(E)</math> como densidade de estados, interpretamos como o número de estados com energia E.
A amostragem de Wang-Landau faz caminhadas aleatórias no espaço de energia mudando os estados de spins aleatoriamente selecionados, porém esta mudança só é aceita com probabilidade proporcional a reciproca da densidade de estados. Durante a caminhada também se acumula o número de vezes que uma energia é visitada durante a caminhada <math>H(E)</math>, isto é, ao visitarmos a energia <math>E_1</math> faz-se a atualização da variável <math> H(E_1) \rightarrow H(E_1) + 1 </math>. Por outro lado, a atualização da densidade de estados <math>g(E_1)</math> se da por um fator multiplicativo <math>f</math> (<math> g(E_1) \rightarrow g(E_1) \times f </math>) controlado ao longo da simulação para que seja muito próximo de 1 ao final das caminhadas.  
A amostragem de Wang-Landau faz caminhadas aleatórias no espaço de energia mudando os estados de spins aleatoriamente selecionados, porém esta mudança só é aceita com probabilidade proporcional a reciproca da densidade de estados. Durante a caminhada também se acumula o número de vezes que uma energia é visitada durante a caminhada <math>H(E)</math>, isto é, ao visitarmos a energia <math>E_1</math> faz-se a atualização da variável <math> H(E_1) \rightarrow H(E_1) + 1 </math>. Por outro lado, a atualização da densidade de estados <math>g(E_1)</math> se da por um fator multiplicativo <math>f</math> (<math> g(E_1) \rightarrow g(E_1) \times f </math>) controlado ao longo da simulação para que seja muito próximo de 1 ao final das caminhadas.  


Linha 27: Linha 59:
#*Se não aceitarmos a mudança de estado do spin, fazemos as atualizações de <math>g(E)</math> e <math>H(E)</math> como <math>g(E_1) \rightarrow g(E_1) \times f</math> e <math>H(E_1) \rightarrow H(E_1) +1</math> respectivamente, de maneira a recontar o estado <math>E_1</math>.  
#*Se não aceitarmos a mudança de estado do spin, fazemos as atualizações de <math>g(E)</math> e <math>H(E)</math> como <math>g(E_1) \rightarrow g(E_1) \times f</math> e <math>H(E_1) \rightarrow H(E_1) +1</math> respectivamente, de maneira a recontar o estado <math>E_1</math>.  
#*Destaca-se que em ambos os casos usamos <math>\ln(g(E)) = \rightarrow \ln(g(E)) + \ln(f)</math>, pois ao longo da simulação acabamos usando números muito grandes.  
#*Destaca-se que em ambos os casos usamos <math>\ln(g(E)) = \rightarrow \ln(g(E)) + \ln(f)</math>, pois ao longo da simulação acabamos usando números muito grandes.  
#Faz-se esta caminhada aleatória nos diferentes estados do sistema até que o histograma $H(E)$ esteja aproximadamente plano. O critério para decidir se um histograma está plano é arbitrário. Para um hamiltoniano Ising 2D este critério pode ser definido tão alto quanto 95\%(i.e. todos os valores de $H(E)$ devem ser pelo menos 95\% de $\left<H(E)\right>$), porém valores mais altos que isso podem resultar no programa nunca identificando o histograma como plano.
#Faz-se esta caminhada aleatória nos diferentes estados do sistema até que o histograma <math>H(E)</math> esteja aproximadamente plano.  
#Checa-se se $H(E)$ está plano a cada 10000 passos MC. Quando está plano, podemos dizer que todos os estados foram visitados uma quantidade de vezes aproximadamente igual e a densidade de estados $g(E)$ converge ao valor real com precisão da ordem de $f$.
#*O critério para decidir se um histograma está plano é arbitrário. Para um [https://pt.wikipedia.org/wiki/Hamiltoniano_(mec%C3%A2nica_qu%C3%A2ntica) hamiltoniano] de Ising 2D este critério pode ser definido tão alto quanto 95% (i.e. todos os valores de <math>H(E)</math> devem ser pelo menos 95% de <math> \langle H(E) \rangle </math>), porém valores mais altos que isso podem resultar no programa nunca identificando o histograma como plano.
#Reduz-se o fator $f$ da seguinte maneira $f_1 \rightarrow \sqrt{f_0}$, reinicia-se o histograma $H(E) = 0$ e recomeça-se a caminhada aleatória com este novo fator $f$.(Todos os parâmetros não mencionados neste passo permanecem intocados.)
#Checa-se se <math>H(E)</math> está plano a cada 10000 passos MC. Quando está plano, podemos dizer que todos os estados foram visitados uma quantidade de vezes aproximadamente igual e a densidade de estados <math>g(E)</math> converge ao valor real com precisão da ordem de <math>f</math>. (O critério de 10000 passos é arbitrário e pode ser modificado de acordo com o propósito da simulação).
#Continuamos executando os passos VI-VIII reduzindo o fator $f$ segundo a seguinte expressão $f_{i+1} = \sqrt{f_{i}}$
#*Um passo MC corresponde a <math>L^2 = N</math> passos.
#Encerra-se a simulação quando $f_{final}$ estiver da ordem do erro desejado. Claro que $f_{final}$ pode ser escolhido arbitrariamente pequeno, mas sempre com um certo limite razoável ($10^{-6}-10^{-8}$), ou a simulação pode tomar tempos não razoáveis para completar.
#Reduz-se o fator <math>f</math> da seguinte maneira <math>f_1 \rightarrow \sqrt{f_0}</math>, reinicia-se o histograma <math>H(E) = 0</math> e recomeça-se a caminhada aleatória com este novo fator <math>f</math>. (Todos os parâmetros não mencionados neste passo permanecem intocados).
#Continuamos executando os passos 5-7, reduzindo o fator <math>f</math> segundo a seguinte expressão <math>f_{i+1} = \sqrt{f_{i}}</math>
#Encerra-se a simulação quando <math>f_{final}</math> estiver da ordem do erro desejado.  
#*Claro que <math>f_{final}</math> pode ser escolhido arbitrariamente pequeno, mas sempre com um certo limite razoável <math>(10^{-6}-10^{-8})</math>, ou a simulação pode tomar tempos não razoáveis para completar.
 
==Observações sobre o algoritmo==
Nesta seção discutimos algumas observações importantes a se levar em conta na implementação do método de amostragem de Wang-Landau
===Fator de modificação f===
Quando tratamos da atualização do fator <math>f</math>, a expressão <math>f_{i+1} = \sqrt{f_{i}}</math> é apenas uma recomendação, uma vez que outros valores de <math>n>1</math> podem ser escolhidos para uma atualização do tipo <math>f_{i+1} = f_{i}^{(1/n)}</math>. Não obstante, <math>n = 2</math> é adequado para o sistema estudado (Ising 2D), resultando em boa acurácia em relativamente pouco tempo de simulação.
 
===Implementação paralela===
A simulação pode ser melhorada ainda fazendo múltiplas caminhadas aleatórias paralelamente no espaço de energias. Restringindo o alcance das caminhadas proporcionalmente com o número de caminhantes em paralelo (e.g. No caso de 2 caminhantes simultâneos, dividimos o espaço de energias em 2 e restringimos um caminhante para a metade inferior das energias, e o outro para a parte superior) e depois juntando as densidades de estados resultantes.
<ref>Vogel, Thomas, et al. "Generic, hierarchical framework for massively parallel Wang-Landau sampling." Physical review letters 110.21 (2013): 210603.</ref>
 
===Balanço detalhado===
A condição de balanço detalhado inicialmente não é satisfeita uma vez que <math>g(E)</math> é constantemente modificada durante a caminhada aleatória. Porém após várias iterações, a condição é satisfeita a medida que <math>f</math> se aproxima de 1. Observa-se que se <math>P(E_1 \rightarrow E_2)</math> é a probabilidade de transição da energia <math>E_1</math> para a energia <math>E_2</math>, utilizando a equação do passo 4 do algoritmo  temos que:
 
<math>
  \frac{P(E_1 \rightarrow E_2)}{P(E_2 \rightarrow E_1)} = \frac{g(E_1)}{g(E_2)}.
</math>
 
Podemos reescrever a equação de uma forma mais familiar
 
<math>
  \frac{1}{g(E_1)}P(E_1 \rightarrow E_2) = \frac{1}{g(E_2)}P(E_2 \rightarrow E_1)
</math>
 
que é a condição de balanceamento detalhado, uma vez que interpretamos que <math>1/g(E_1)</math> como a probabilidade do sistema possuir a energia <math>E_1</math> e analogamente para <math>E_2</math>. Concluímos então que a condição de balanço detalhado é satisfeita com precisão proporcional a <math>\ln(f)</math>.
===Escalabilidade===
Quando analisamos um modelo de Ising 2D de tamanho <math>L\times L</math>, temos que o número de configurações <math>(2^N, N = L^2)</math> aumenta exponencialmente com <math>N</math>, enquanto o número de possíveis energias é por volta de <math>2N</math> e aumenta linearmente com N. Implicando que temos uma escalabilidade muito boa para as caminhadas no espaço de energia quando o objetivo é estimar <math>g(E)</math> uma vez que um aumento no tamanho da grade não implica em um aumento exponencial, mas sim linear, no tempo de execução.
===Normalização===
É necessário ressaltar que após a simulação completa, o algoritmo de Wang-Landau nos fornece apenas a densidade de estados relativa. Para extrairmos a real densidade de estados <math>g_n(E)</math> é necessário que utilizemos uma das duas condições: Que o número total de estados possíveis é <math>\sum_Eg_n(E) = Q^N</math> ou que o numero de estados fundamentais é <math>Q</math> (onde <math>Q = 2</math> para o modelo de Ising 2D pois os spins possuem apenas dois estados).
 
Pela primeira condição, podemos obter a densidade de estados normalizada através da equação
 
<math>\ln[g_n(E)] = \ln[g(E)] - \ln[\sum_Eg(E)] + N\ln(2)</math>,
 
enquanto pela segunda condição temos que 
 
<math>\ln[g_n(E)] = \ln[g(E)] - \ln[g(E = -2N)] + \ln(2)</math>.
 
A segunda normalização é preferível pois garante precisão para estados de menor energia, o que é necessário para o calculo de parâmetros termodinâmicos a baixas temperaturas.
 
===Grandezas Termodinâmicos===
Uma vez que temos a densidade de estados, podemos calcular diversas grandezas termodinâmicos, como a [https://pt.wikipedia.org/wiki/Energia_interna energia interna] <math>U(T)</math>, [https://pt.wikipedia.org/wiki/Calor_espec%C3%ADfico calor especifico] <math>C(T)</math>, [https://pt.wikipedia.org/wiki/Energia_livre_de_Helmholtz energia livre de Helmholtz] <math>F(T)</math> e [https://pt.wikipedia.org/wiki/Entropia entropia] <math>S(T)</math> através das seguintes equações:
<div id="equações para parâmetros termodinâmicos">
      <math>U(T) = \frac{\sum_EEg(E)e^{-E/k_BT}}{\sum_Eg(E)e^{-E/k_BT}}</math>
 
      <math>C(T) = \frac{\langle E^2\rangle - \langle E\rangle^2}{k_BT^2}</math>
 
      <math>F(T) = -k_BT\ln\left(\sum_Eg(E)e^{-E/k_BT}\right)</math>
 
      <math>S(T) = \frac{U(T) - F(T)}{T}</math>
</div>
Estas grandezas dependem apenas da temperatura uma vez que já foi encontrado a densidade de estados <math>g(E)</math>, o que contorna os problemas de ''critical slowing down'' na temperatura crítica pois a simulação não precisa ser refeita para cada uma das temperaturas, por consequência também dispensa a necessidade do método de repesagem de histogramas, contornando o problema de estatística fraca nas asas dos histogramas.
 
==Exemplo: Modelo de Ising 2D por amostragem de Wang-Landau==
Podemos aplicar o algoritmo de Wang-Landau para um sistema [https://pt.wikipedia.org/wiki/Ferromagnetismo#:~:text=Ferromagnetismo%20%C3%A9%20o%20mecanismo%20b%C3%A1sico,diferentes%20de%20magnetismo%20s%C3%A3o%20distinguidos. ferromagnético] 2D de Ising com interação de primeiros vizinhos e uma grade quadrada <math>L\times L</math> com [https://pt.wikipedia.org/wiki/Condi%C3%A7%C3%B5es_de_fronteira_peri%C3%B3dicas condições de contorno periódicas] no qual o hamiltoniano é dado por
 
<math>\mathcal{H} = -\sum_{\langle i,j\rangle}\sigma_{i}\sigma_j</math>,
 
onde <math>\sigma_{i} = +1</math> para spin para cima, <math>\sigma_{i} = -1</math> para baixo e <math>\langle i,j\rangle</math> indica a soma entre os primeiros vizinhos. Podemos calcular os erros entre os valores obtidos através da simulação e os valores exatos conhecidos para esse sistema através da seguinte expressão
 
<math>\varepsilon(X) = \frac{|X_{sim}-X_{exato}|}{X_{exato}}</math>,
 
onde <math>X_{sim}</math> é o valor obtido na simulação para o parâmetro e <math>X_{exato}</math> é o valor exato para o mesmo parâmetro. Na Fig. 1 temos a densidade de estados <math>g(E)</math> para um sistema <math>16 \times 16</math> com <math>f_{final} \sim 10^{-6}</math> e com critério para um histograma <math>H(E)</math> plano de 80%. Com esses parâmetros, a ''Google Compute Engine'' padrão para [https://pt.wikipedia.org/wiki/Python Python3] do Google Colab realiza a simulação em <math>\sim 40s</math>. Para as Figs. 1,2,3,4,5,6 os valores exatos estão sobrepostos aos valores obtidos em simulação como uma linha pontilhada.
 
[[Arquivo: erickStateDensity2.png|400px|thumb|center| Figura 1: Logaritmo da densidade de estados <math>\ln(g(E))</math> para o modelo de Ising 2D com <math>L = 16</math>, há uma brusca queda na densidade de estados perto das pontas uma vez que as energias <math>-2L^2 + 4</math> e <math>2L^2 - 4</math> são inalcançáveis pelo sistema.]]
 
Fazendo uso da [[#equação de distribuição de probabilidade canônica]], podemos encontrar a distribuição de probabilidade das energias para uma determinada temperatura <math>T</math>. Devido a natureza do algoritmo, estas distribuições que por meios convencionais levariam tempo demasiadamente grande para ser calculado ou necessitariam de repesagem de histograma são instantaneamente calculadas com o algoritmo de Wang-Landau independentemente da temperatura. Na Fig.2 temos um exemplo das distribuições, incluindo uma próxima da temperatura critica <math>T_c</math>.
 
[[Arquivo: erickProbs2.png|400px|thumb|center| Figura 2: A distribuição canônica <math>P(E,T) = g(E)e^{-E/K_BT}</math> em 3 temperaturas, incluindo a temperatura crítica <math>T_c</math> para o modelo de Ising 2D com <math>L = 16</math>.]]
 
Fazendo uso das [[#equações para parâmetros termodinâmicos]], podemos calculá-los para este sistema em função da temperatura quase que instantaneamente. Estes são ilustrados nas Figs. 3,4,5,6.
<div><ul>
<li style="display: inline-block;">[[Arquivo: erickU2.png|400px|thumb|center| Figura 3: A Energia interna do sistema em função da temperatura para o modelo de Ising 2D com <math>L = 16</math>. Perceba que a inflexão se da aproximadamente na temperatura crítica <math>k_BT_c \approx 2.3</math>.]]</li>
<li style="display: inline-block;">[[Arquivo: erickC2.png|400px|thumb|center| Figura 4: O Calor Específico do sistema em função da temperatura para o modelo de Ising 2D com <math>L = 16</math>. Perceba que o pico se da aproximadamente na temperatura crítica <math>k_BT_c \approx 2.3</math>.]]</li>
<li style="display: inline-block;">[[Arquivo: erickF2.png|400px|thumb|center| Figura 5: A Energia Livre do sistema em função da temperatura para o modelo de Ising 2D com <math>L = 16</math>. Perceba que o começo da queda se da aproximadamente na temperatura crítica <math>k_BT_c \approx 2.3</math>.]]</li>
<li style="display: inline-block;">[[Arquivo: erickS2.png|400px|thumb|center| Figura 6: A Entropia do sistema em função da temperatura para o modelo de Ising 2D com <math>L = 16</math>. Perceba que a inflexão se da aproximadamente na temperatura crítica <math>k_BT_c \approx 2.3</math>.]]</li>
</ul></div>
 
==Código do Exemplo==
[[Exemplo de código para a amostragem de Wang-Landau em Python3]]
 
==Referências==
<references />

Edição atual tal como às 01h25min de 19 de outubro de 2022

O algoritmo de amostragem de Wang-Landau é um método de amostragem para simulações de Monte Carlo introduzido por F. Wang and D. P. Landau em 2001 [1] (com uma versão didática publicada em 2004 por D.P Landau, Shan-Ho Tsai e M. Exlerem), [2] que apresenta diversas vantagens sobre outros métodos para sistemas de spins. Dentre eles podemos citar o algoritmo de Metropolis [3] , o algoritmo de clustering de Wolff [4] , ou em um modelamento de ensamble multi-canônico [5] . Nestes dois últimos métodos são utilizados métodos de repesagem de histogramas[6] que são limitados em sistemas grandes devido as baixa qualidade estatística nas asas do histograma. Dentro deste contexto, o algoritmo Wang-Landau promete resolver problemas encontrados em amostragens convencionais como o critical slowing down para temperaturas próximas da temperatura crítica utilizando-se de caminhadas aleatórias controladas para mapear a densidade de estados de um sistema, sem depender do método de repesagem de histogramas.


A maioria dos algoritmos de amostragem convencionais geram uma distribuição canônica não normalizada

para uma determinada temperatura , Geralmente estas distribuições são estreitas e se faz necessário múltiplas simulações para obter algum parâmetro termodinâmico para uma distribuição significantemente grande de temperaturas, na Fig. 2 vemos exemplos destas distribuições. Como não depende de temperatura do sistema, se pudermos encontrá-lo para todo , podemos encontrar a função de partição

e o sistema esta essencialmente resolvido, uma vez que grande parte dos parâmetros termodinâmicos podem ser derivados de . Além disso, a amostragem de Wang-Landau é provada ser útil em diversas aplicações como o antiferromagneto de Potts [7] , sistemas de spins aleatórios[8] , sistemas quânticos[9], etc... .

Descrição do algoritmo de Wang Landau

Descreveremos o funcionamento do algoritmo de Wang-Landau num sistema de spins clássicos de 2 estados numa grade bidimensional com valores discretos de energia e sem campo magnético(Modelo de Ising2D). Portanto quando nos referirmos a como densidade de estados, interpretamos como o número de estados com energia E. A amostragem de Wang-Landau faz caminhadas aleatórias no espaço de energia mudando os estados de spins aleatoriamente selecionados, porém esta mudança só é aceita com probabilidade proporcional a reciproca da densidade de estados. Durante a caminhada também se acumula o número de vezes que uma energia é visitada durante a caminhada , isto é, ao visitarmos a energia faz-se a atualização da variável . Por outro lado, a atualização da densidade de estados se da por um fator multiplicativo () controlado ao longo da simulação para que seja muito próximo de 1 ao final das caminhadas.

Podemos descrever os passos do algoritmo da seguinte maneira:

  1. Inicializamos as densidades de energias com para todo , da mesma forma para todo .
  2. Inicializamos e um sistema de spins de valor 1 e -1 aleatoriamente distribuídos.
    • O valor de é arbitrário e deve ser escolhido não muito pequeno, pois irá fazer com que a simulação demore muito tempo para explorar diversas energias, por outro lado se escolhido muito grande, levará a erros estatísticos significativos.
  3. Começamos a caminhada inicial escolhendo aleatoriamente um dos spins e mudando o seu estado.
  4. Se denotamos como a energia antes da mudança de estado do spin selecionado e como a energia após, aceitamos este novo estado com a seguinte probabilidade:
    • Se aceitarmos a mudança de estado do spin, fazemos as atualizações de e como e respectivamente.
    • Se não aceitarmos a mudança de estado do spin, fazemos as atualizações de e como e respectivamente, de maneira a recontar o estado .
    • Destaca-se que em ambos os casos usamos , pois ao longo da simulação acabamos usando números muito grandes.
  5. Faz-se esta caminhada aleatória nos diferentes estados do sistema até que o histograma esteja aproximadamente plano.
    • O critério para decidir se um histograma está plano é arbitrário. Para um hamiltoniano de Ising 2D este critério pode ser definido tão alto quanto 95% (i.e. todos os valores de devem ser pelo menos 95% de ), porém valores mais altos que isso podem resultar no programa nunca identificando o histograma como plano.
  6. Checa-se se está plano a cada 10000 passos MC. Quando está plano, podemos dizer que todos os estados foram visitados uma quantidade de vezes aproximadamente igual e a densidade de estados converge ao valor real com precisão da ordem de . (O critério de 10000 passos é arbitrário e pode ser modificado de acordo com o propósito da simulação).
    • Um passo MC corresponde a passos.
  7. Reduz-se o fator da seguinte maneira , reinicia-se o histograma e recomeça-se a caminhada aleatória com este novo fator . (Todos os parâmetros não mencionados neste passo permanecem intocados).
  8. Continuamos executando os passos 5-7, reduzindo o fator segundo a seguinte expressão
  9. Encerra-se a simulação quando estiver da ordem do erro desejado.
    • Claro que pode ser escolhido arbitrariamente pequeno, mas sempre com um certo limite razoável , ou a simulação pode tomar tempos não razoáveis para completar.

Observações sobre o algoritmo

Nesta seção discutimos algumas observações importantes a se levar em conta na implementação do método de amostragem de Wang-Landau

Fator de modificação f

Quando tratamos da atualização do fator , a expressão é apenas uma recomendação, uma vez que outros valores de podem ser escolhidos para uma atualização do tipo . Não obstante, é adequado para o sistema estudado (Ising 2D), resultando em boa acurácia em relativamente pouco tempo de simulação.

Implementação paralela

A simulação pode ser melhorada ainda fazendo múltiplas caminhadas aleatórias paralelamente no espaço de energias. Restringindo o alcance das caminhadas proporcionalmente com o número de caminhantes em paralelo (e.g. No caso de 2 caminhantes simultâneos, dividimos o espaço de energias em 2 e restringimos um caminhante para a metade inferior das energias, e o outro para a parte superior) e depois juntando as densidades de estados resultantes. [10]

Balanço detalhado

A condição de balanço detalhado inicialmente não é satisfeita uma vez que é constantemente modificada durante a caminhada aleatória. Porém após várias iterações, a condição é satisfeita a medida que se aproxima de 1. Observa-se que se é a probabilidade de transição da energia para a energia , utilizando a equação do passo 4 do algoritmo temos que:

Podemos reescrever a equação de uma forma mais familiar

que é a condição de balanceamento detalhado, uma vez que interpretamos que como a probabilidade do sistema possuir a energia e analogamente para . Concluímos então que a condição de balanço detalhado é satisfeita com precisão proporcional a .

Escalabilidade

Quando analisamos um modelo de Ising 2D de tamanho , temos que o número de configurações aumenta exponencialmente com , enquanto o número de possíveis energias é por volta de e aumenta linearmente com N. Implicando que temos uma escalabilidade muito boa para as caminhadas no espaço de energia quando o objetivo é estimar uma vez que um aumento no tamanho da grade não implica em um aumento exponencial, mas sim linear, no tempo de execução.

Normalização

É necessário ressaltar que após a simulação completa, o algoritmo de Wang-Landau nos fornece apenas a densidade de estados relativa. Para extrairmos a real densidade de estados é necessário que utilizemos uma das duas condições: Que o número total de estados possíveis é ou que o numero de estados fundamentais é (onde para o modelo de Ising 2D pois os spins possuem apenas dois estados).

Pela primeira condição, podemos obter a densidade de estados normalizada através da equação

,

enquanto pela segunda condição temos que

.

A segunda normalização é preferível pois garante precisão para estados de menor energia, o que é necessário para o calculo de parâmetros termodinâmicos a baixas temperaturas.

Grandezas Termodinâmicos

Uma vez que temos a densidade de estados, podemos calcular diversas grandezas termodinâmicos, como a energia interna , calor especifico , energia livre de Helmholtz e entropia através das seguintes equações:

      
      
      
      

Estas grandezas dependem apenas da temperatura uma vez que já foi encontrado a densidade de estados , o que contorna os problemas de critical slowing down na temperatura crítica pois a simulação não precisa ser refeita para cada uma das temperaturas, por consequência também dispensa a necessidade do método de repesagem de histogramas, contornando o problema de estatística fraca nas asas dos histogramas.

Exemplo: Modelo de Ising 2D por amostragem de Wang-Landau

Podemos aplicar o algoritmo de Wang-Landau para um sistema ferromagnético 2D de Ising com interação de primeiros vizinhos e uma grade quadrada com condições de contorno periódicas no qual o hamiltoniano é dado por

,

onde para spin para cima, para baixo e indica a soma entre os primeiros vizinhos. Podemos calcular os erros entre os valores obtidos através da simulação e os valores exatos conhecidos para esse sistema através da seguinte expressão

,

onde é o valor obtido na simulação para o parâmetro e é o valor exato para o mesmo parâmetro. Na Fig. 1 temos a densidade de estados para um sistema com e com critério para um histograma plano de 80%. Com esses parâmetros, a Google Compute Engine padrão para Python3 do Google Colab realiza a simulação em . Para as Figs. 1,2,3,4,5,6 os valores exatos estão sobrepostos aos valores obtidos em simulação como uma linha pontilhada.

Figura 1: Logaritmo da densidade de estados para o modelo de Ising 2D com , há uma brusca queda na densidade de estados perto das pontas uma vez que as energias e são inalcançáveis pelo sistema.

Fazendo uso da #equação de distribuição de probabilidade canônica, podemos encontrar a distribuição de probabilidade das energias para uma determinada temperatura . Devido a natureza do algoritmo, estas distribuições que por meios convencionais levariam tempo demasiadamente grande para ser calculado ou necessitariam de repesagem de histograma são instantaneamente calculadas com o algoritmo de Wang-Landau independentemente da temperatura. Na Fig.2 temos um exemplo das distribuições, incluindo uma próxima da temperatura critica .

Figura 2: A distribuição canônica em 3 temperaturas, incluindo a temperatura crítica para o modelo de Ising 2D com .

Fazendo uso das #equações para parâmetros termodinâmicos, podemos calculá-los para este sistema em função da temperatura quase que instantaneamente. Estes são ilustrados nas Figs. 3,4,5,6.

  • Figura 3: A Energia interna do sistema em função da temperatura para o modelo de Ising 2D com . Perceba que a inflexão se da aproximadamente na temperatura crítica .
  • Figura 4: O Calor Específico do sistema em função da temperatura para o modelo de Ising 2D com . Perceba que o pico se da aproximadamente na temperatura crítica .
  • Figura 5: A Energia Livre do sistema em função da temperatura para o modelo de Ising 2D com . Perceba que o começo da queda se da aproximadamente na temperatura crítica .
  • Figura 6: A Entropia do sistema em função da temperatura para o modelo de Ising 2D com . Perceba que a inflexão se da aproximadamente na temperatura crítica .

Código do Exemplo

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

Referências

  1. F. Wang and D. P. Landau, ‘‘Efficient, multiple-range random walk algorithm to calculate the density of states,’’ Phys. Rev. Lett. 86, 2050–2053 !2001"; ‘‘Determining the density of states for classical statistical models: A random walk algorithm to produce a flat histogram,’’ Phys. Rev. E 64, 056101 !2001".
  2. Landau, D. P., Shan-Ho Tsai, and M. Exler. "A new approach to Monte Carlo simulations in statistical physics: Wang-Landau sampling." American Journal of Physics 72.10 (2004): 1294-1302.
  3. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, ‘‘Equation of state calculations by fast computing machines,’’ J. Chem. Phys. 21, 1087–1092 !1953".
  4. U. Wolff, ‘‘Collective Monte Carlo updating for spin systems,’’ Phys. Rev. Lett. 62, 361–364 !1989".
  5. B. A. Berg and T. Neuhaus, ‘‘Multicanonical ensemble: A new approach to simulate first-order phase transitions,’’ Phys. Rev. Lett. 68, 9–12 !1992"
  6. A. M. Ferrenberg and R. H. Swendsen, ‘‘New Monte Carlo technique for studying phase transitions,’’ Phys. Rev. Lett. 61, 2635–2638 !1988"; ‘‘Optimized Monte Carlo data analysis,’’ ibid. 63, 1195–1198 !1989".
  7. C. Yamaguchi and Y. Okabe, ‘‘Three-dimensional antiferromagnetic q-state Potts models: application of the Wang-Landau algorithm,’’ J. Phys. A 34, 8781– 8794 !2001"
  8. Y. Okabe, Y. Tomita, and C. Yamaguchi, ‘‘Application of new Monte Carlo algorithms to random spin systems,’’ Comput. Phys. Commun. 146, 63– 68 !2002".
  9. 2M. Troyer, S. Wessel, and F. Alet, ‘‘Flat histogram methods for quantum systems: Algorithms to overcome tunneling problems and calculate the free energy,’’ Phys. Rev. Lett. 90, 120201 !2003".
  10. Vogel, Thomas, et al. "Generic, hierarchical framework for massively parallel Wang-Landau sampling." Physical review letters 110.21 (2013): 210603.