Amostragem de Wang-Landau

O algoritmo de amostragem de Wang-Landau é um método de amostragem para simulações de Monte Carlo introduzido por D.P Landau, Shan-Ho Tsai e M. Exlerem em 2004 com suas primeiras concepções traçadas por F. Wang and 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. 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 fazer uso de qualquer 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. 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})\rightarrow 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})\rightarrow 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}\rightarrow E_{2})={\text{min}}\left({\frac {g(E_{1})}{g(E_{2})}},1\right)$
   • Se aceitarmos a mudança de estado do spin, fazemos as atualizações de $g(E)$ e $H(E)$ como $g(E_{2})\rightarrow g(E_{2})\times f$ e $H(E_{2})\rightarrow 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})\rightarrow g(E_{1})\times f$ e $H(E_{1})\rightarrow H(E_{1})+1$ respectivamente, de maneira a recontar o estado $E_{1}$.
   • Destaca-se que em ambos os casos usamos $\ln(g(E))=\rightarrow \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 $\langle H(E)\rangle$), 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$.
   • Um passo MC corresponde a $L^{2}=N$ passos.
7. 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).
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 grande parte dos sistemas estudados, 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}\rightarrow 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}\rightarrow E_{2})}{P(E_{2}\rightarrow 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}\rightarrow E_{2})={\frac {1}{g(E_{2})}}P(E_{2}\rightarrow 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 há baixas temperaturas.

Parâmetros Termodinâmicos

Uma vez que temos a densidade de estados, podemos calcular diversos parâmetros 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 {\langle E^{2}\rangle -\langle E\rangle ^{2}}{k_{B}T^{2}}}$

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

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

Estes parâmetros termodinâmicos 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 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

${\mathcal {H}}=-\sum _{\langle i,j\rangle }\sigma _{i}\sigma _{j}$

onde $\sigma _{i}=+1$ para spin para cima, $\sigma _{i}=-1$ para baixo e $\langle i,j\rangle$ indica a soma entre os primeiros vizinhos. 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$.

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}+8$ e $2L^{2}-8$ 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$.

Exemplo de código

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