Amostragem de Wang-Landau

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 (com uma versão didática publicada em 2004 por D.P Landau, Shan-Ho Tsai e M. Exlerem), 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. 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 $T_c$ utilizando-se de caminhadas aleatórias controladas para mapear a densidade de estados $g(E)$ 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

$P(E,T)=g(E)e^{-E/K_{B}T},$

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, na Fig. 2 vemos exemplos destas distribuições. Como $g(E)$ não depende de temperatura do sistema, se pudermos encontrá-lo para todo $E$, podemos encontrar a função de partição

$Z(T)=\sum _{E}g(E)e^{-E/K_{B}T},$

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... .

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 $g(E)$ 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 $H(E)$, isto é, ao visitarmos a energia $E_1$ faz-se a atualização da variável $H(E_1)\to H(E_1)+1$. Por outro lado, a atualização da densidade de estados $g(E_1)$ se da por um fator multiplicativo $f$ ($g(E_1)\to g(E_1)\times f$) 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 $g(E)=1$ para todo $E$, da mesma forma $H(E)=0$ para todo $E$.
2. Inicializamos $f=f_0=e\approx 2.71828$ e um sistema $L\times L$ de spins de valor 1 e -1 aleatoriamente distribuídos.
   • O valor de $f$ é 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 $E_1$ como a energia antes da mudança de estado do spin selecionado e $E_2$ como a energia após, aceitamos este novo estado com a seguinte probabilidade: $P(E_1\to E_2)=\text{min}(\frac{g(E_1)}{g(E_2)},1)$
   • Se aceitarmos a mudança de estado do spin, fazemos as atualizações de $g(E)$ e $H(E)$ como $g(E_2)\to g(E_2)\times f$ e $H(E_2)\to H(E_2)+1$ respectivamente.
   • Se não aceitarmos a mudança de estado do spin, fazemos as atualizações de $g(E)$ e $H(E)$ como $g(E_1)\to g(E_1)\times f$ e $H(E_1)\to H(E_1)+1$ respectivamente, de maneira a recontar o estado $E_1$.
   • Destaca-se que em ambos os casos usamos $\ln (g(E))=\to \ln (g(E))+\ln (f)$, 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 $H(E)$ 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 $H(E)$ devem ser pelo menos 95% de $⟨H(E)⟩$), porém valores mais altos que isso podem resultar no programa nunca identificando o histograma como plano.
6. 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 de 10000 passos é arbitrário e pode ser modificado de acordo com o propósito da simulação).
   • Um passo MC corresponde a $L^2=N$ passos.
7. Reduz-se o fator $f$ da seguinte maneira $f_1\to \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).
8. Continuamos executando os passos 5-7, reduzindo o fator $f$ segundo a seguinte expressão $f_{i+1}=\sqrt{f_i}$
9. 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.

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 $f$, a expressão $f_{i+1}=\sqrt{f_i}$ é apenas uma recomendação, uma vez que outros valores de $n>1$ podem ser escolhidos para uma atualização do tipo $f_{i+1}=f_i^{(1/n)}$. Não obstante, $n=2$ é 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.

Balanço detalhado

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

$\frac{P(E_1\to E_2)}{P(E_2\to E_1)}=\frac{g(E_1)}{g(E_2)}.$

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

$\frac{1}{g(E_1)}P(E_1\to E_2)=\frac{1}{g(E_2)}P(E_2\to E_1)$

que é a condição de balanceamento detalhado, uma vez que interpretamos que $1/g(E_1)$ como a probabilidade do sistema possuir a energia $E_1$ e analogamente para $E_2$. Concluímos então que a condição de balanço detalhado é satisfeita com precisão proporcional a $\ln (f)$.

Escalabilidade

Quando analisamos um modelo de Ising 2D de tamanho $L\times L$, temos que o número de configurações $(2^N,N=L^2)$ aumenta exponencialmente com $N$, enquanto o número de possíveis energias é por volta de $2N$ e aumenta linearmente com N. Implicando que temos uma escalabilidade muito boa para as caminhadas no espaço de energia quando o objetivo é estimar $g(E)$ 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 $g_n(E)$ é necessário que utilizemos uma das duas condições: Que o número total de estados possíveis é $\sum _E{g}_n(E)=Q^N$ ou que o numero de estados fundamentais é $Q$ (onde $Q=2$ 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

$\ln [g_n(E)]=\ln [g(E)]-\ln [\sum _{E}g(E)]+N\ln (2)$,

enquanto pela segunda condição temos que

$\ln [g_n(E)]=\ln [g(E)]-\ln [g(E=-2N)]+\ln (2)$.

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 $U(T)$, calor especifico $C(T)$, energia livre de Helmholtz $F(T)$ e entropia $S(T)$ através das seguintes equações:

      $U(T)=\frac{\sum _{E}Eg(E)e^{-E/k_{B}T}}{\sum _{E}g(E)e^{-E/k_{B}T}}$

      $C(T)=\frac{⟨E^2⟩-⟨E{⟩}^2}{k_B{T}^2}$

      $F(T)=-k_{B}T\ln (\sum _{E}g(E)e^{-E/k_{B}T})$

      $S(T)=\frac{U(T)-F(T)}{T}$

Estas grandezas dependem apenas da temperatura uma vez que já foi encontrado a densidade de estados $g(E)$, 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 $L\times L$ com condições de contorno periódicas no qual o hamiltoniano é dado por

${\mathscr{H}}=-\sum _{⟨i,j⟩}{\sigma }_i{\sigma }_j$,

onde ${\sigma }_i=+1$ para spin para cima, ${\sigma }_i=-1$ para baixo e $⟨i,j⟩$ 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

$\epsilon (X)=\frac{|X_{sim}-X_{exato}|}{X_{exato}}$,

onde $X_{sim}$ é o valor obtido na simulação para o parâmetro e $X_{exato}$ é o valor exato para o mesmo parâmetro. Na Fig. 1 temos a densidade de estados $g(E)$ para um sistema $16\times 16$ com $f_{final}\sim 10^{-6}$ e com critério para um histograma $H(E)$ plano de 80%. Com esses parâmetros, a Google Compute Engine padrão para Python3 do Google Colab realiza a simulação em $\sim 40s$. 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 $\ln (g(E))$ para o modelo de Ising 2D com $L=16$, há uma brusca queda na densidade de estados perto das pontas uma vez que as energias $-2L^2+4$ e $2L^2-4$ 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 $T$. 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 $T_c$.

Figura 2: A distribuição canônica $P(E,T)=g(E)e^{-E/K_{B}T}$ em 3 temperaturas, incluindo a temperatura crítica $T_c$ para o modelo de Ising 2D com $L=16$.

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 $L=16$. Perceba que a inflexão se da aproximadamente na temperatura crítica $k_B{T}_c\approx 2.3$.

• 

  Figura 4: O Calor Específico do sistema em função da temperatura para o modelo de Ising 2D com $L=16$. Perceba que o pico se da aproximadamente na temperatura crítica $k_B{T}_c\approx 2.3$.

• 

  Figura 5: A Energia Livre do sistema em função da temperatura para o modelo de Ising 2D com $L=16$. Perceba que o começo da queda se da aproximadamente na temperatura crítica $k_B{T}_c\approx 2.3$.

• 

  Figura 6: A Entropia do sistema em função da temperatura para o modelo de Ising 2D com $L=16$. Perceba que a inflexão se da aproximadamente na temperatura crítica $k_B{T}_c\approx 2.3$.

Código do Exemplo

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