Algoritmo de Wang-Landau

Nomes: Rafael Abel da Silveira e William Machado Pantaleão

Introdução

Simulações computacionais, como o método de Monte Carlo, são vastamente utilizadas para estudar transições de fase e fenômenos críticos. O método padrão para simulações de Monte Carlo é o algoritmo de Metropolis, entretanto, algoritmos novos e mais eficientes são usados em simulações modernas, como o algoritmo de Wang-Landau. Ao contrário dos métodos convencionais de Monte Carlo, que geram diretamente uma distribuição canônica $g(E)e^{-E/k_{B}T}$ a uma dada temperatura $T$, a abordagem de Wang-Landau estima a densidade de estados $g(E)$ diretamente por meio de um passeio aleatório, que produz um histograma plano no espaço de energia $H(E)$.

Mesmo para modelos que podem ser resolvidos analiticamente, a densidade de estados $g(E)$ não pode ser determinada para sistemas maiores . Com o algoritmo de Wang-Landau, é possível obter a $g(E)$ a partir de um passeio aleatório. A estimativa para $g(E)$ é melhorada a cada etapa do passeio aleatório, usando um fator de modificação $f$ cuidadosamente controlado, para produzir um resultado que converge para o valor real rapidamente.

Amostragem de Wang-Landau

No início da simulação, $g(E)$ é desconhecido e fazemos uma estimativa inicial para ele. A abordagem mais simples é definir $g(E)=1$ para todas as energias possíveis $E$. A configuração de spin inicial para toda a rede pode ser escolhida arbitrariamente. Então, uma caminhada aleatória no espaço de energia é iniciada pela formação de estados de teste, cada um dos quais é produzido escolhendo aleatoriamente um spin e alterando seu estado.

Cada vez que uma energia $E$ é visitada, o histograma $H(E)$ é incrementado em 1. A estimativa de $g(E)$ é então modificada por um fator multiplicativo $f$, e o valor atualizado realiza um passeio aleatório adicional no espaço de $E$.

Se $E_{1}$ e $E_{2}$ são as energias antes e depois de um valor de spin ser alterado, a probabilidade de transição da energia $E_{1}$ para $E_{2}$ é dada por

$p(E_{1}\rightarrow E_{2})=min\left({\frac {g(E_{1})}{g(E_{2})}},1\right).$

A razão das probabilidades de transição de $E_{1}$ para $E_{2}$ e de $E_{2}$ a $E_{1}$ podem ser calculados como

${\frac {p(E_{1}\rightarrow E_{2})}{p(E_{2}\rightarrow E_{1})}}={\frac {g(E_{1})}{g(E_{2})}}.$

Logo, o algoritmo de passeio aleatório satisfaz o equilíbrio detalhado:

${\frac {1}{g(E_{1})}}p(E_{1}\rightarrow E_{2})={\frac {1}{g(E_{2})}}p(E_{2}\rightarrow E_{1}),$

onde $1/g(E1)$ é a probabilidade na energia $E_{1}$ e $p(E_{1}\rightarrow E_{2})$ é a probabilidade de transição de $E_{1}$ para $E_{2}$.

Se o estado de energia $E_{2}$ é aceito, a densidade de estados $g(E_{2})$ é multiplicada pelo fator de modificação $f>1$ de maneira que $g(E_{2})\to f\times g(E_{2})$ e a entrada no histograma para $H(E_{2})$ é atualizada de forma $H(E_{2})\to H(E_{2})+1$. Se o estado de energia não é aceito, a densidade de estados $g(E_{1})$ é multiplicada pelo fator de modificação, $g(E_{1})\to f\times g(E_{1})$ e $H(E_{1})$ é atualizada de forma $H(E_{1})\to H(E_{1})+1$.

Flatness

O procedimento de passeio aleatório é seguido até o histograma $H(E)$ estar reto (do inglês, "flat"), e para determinar isso, a cada $n$ iterações verificamos se todos valores possíveis de $E$ estão a uma distância, no máximo, $x\%$ de $\langle H(E)\rangle$. A variável $x$ é denominada "flatness". Quando o histograma está reto, todos estados de energia foram visitados aproximadamente igualmente.

O número de passos, </math> n </math> que devemos realizar antes de checar deve ser maior que $L*L$ onde $L$ indica o tamanho da rede, para que o algoritmo tenha a oportunidade de visitar cada estado pelo menos uma vez.

Para sistemas simples, podemos utilizar um valor tão alto quanto 95%, entretanto, para este trabalho foi escolhido o valor de 80%.

Fator de modificação

Em geral, como $g(E)$ se torna muito grande, trabalhamos com o logaritmo natural dessas quantidades, $\ln g(E)$. Portanto, cada atualização da densidade de estados é dada por $\ln g(E)\to \ln g(E)+\ln f$. O valor comumente utilizado para o fator de modificação é $f=f_{0}=e$.

Quando o histograma é considerado reto, pelas condições descritas acima, reduzimos o valor de $f$ de forma que o novo valor será $f_{1}=f_{0}/2$, resetamos o histograma $H(E)$ e recomeçamos o passeio aleatório.

A simulação é parada para um valor de $f$ predeterminado. No caso, usamos $f_{final}=\exp(10^{-8})=1.00000001$.

Aplicação ao Modelo de Ising 2D

Modelo de Ising

O modelo de Ising é uma rede 2D, de tamanho $L\times L$ que consiste de uma variável discreta em cada sítio que pode ser usada para representar o momento de dipolo magnético de um átomo Cada sítio pode ter o valor de spin $+1$ ou $-1$.

Para este trabalho, o hamiltoniano de interação pode ser calculado por ${\mathcal {H}}=-\sum _{\langle i,j\rangle }\sigma _{i}\sigma _{j}$ onde $\langle i,j$ indica pares distintos de vizinhos-mais-próximos.

Com a densidade de estados, podemos calcular as seguintes quantidades termodinâmicas:

Energia interna:

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

Calor específico:

$C(T)={\frac {\partial U(T)}{\partial T}}={\frac {\langle E^{2}\rangle -\langle E\rangle ^{2}}{k_{B}T^{2}}}$

Energia livre de Helmoltz:

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

Entropia:

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

Finalmente, podemos também calcular a distribuição canônica usando:

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

Algoritmo

Resumindo, o passo a passo do algoritmo pode ser escrito como:

1. Defino $g(E)=1$ para todos $E$ e o fator de modificação inicial $f_{0}=e$;

2. Aleatoriamente, escolho um spin e troco o seu valor. Aceito a transição com probabilidade $p(E_{1}\to E_{2})=min(g(E_{1})/g(E_{2}),1)$;

3. Modifico a densidade de estados $g(E)\to g(E)\times f$ e atualizo o histograma $H(E)$;

4. Continuo até o histograma estar reto, então diminuo o valor de $\ln f$ pela metade e reseto o histograma $H(E)$;

5. Repito os passos 2-4 até $f<1.00000001$.

6. Obtendo a $g(E)$, posso calcular as quantidades termodinâmicas descritas anteriormente.

Resultados

Densidade de estados

A estimativa da densidade de estados para $L=16$ usando a amostragem de Wang-Landau é mostrada na Fig. 1, junto com os resultados exatos de Beale .

Os fatores de modificação inicial e final para os passeios aleatórios foram Falhou ao verificar gramática (erro de sintaxe): {\displaystyle ln(f_0) = 1 e ln(f_{final}) = 10^{−8}} .

O histograma $H(E)$ foi considerado plano quando todas as entradas não eram inferiores a 80% da média Falhou ao verificar gramática (erro de sintaxe): {\displaystyle \left<H(E)\right>} .

A densidade absoluta de estados na Fig. 1 é obtida pela condição de que o número de estados fundamentais seja 2 para o modelo de Ising 2D.

Com a escala logarítmica usada na Fig. 1, os dados simulados e a solução exata se sobrepõem perfeitamente. Na inserção da Fig. 1, vemos que o erro relativo é, de fato, muito pequeno.

Figura 1: Logaritmo da densidade de estados $\ln g(E)$ do modelo de Ising 2D com $L=16$.

Distribuição canônica

Podemos calcular a distribuição canônica usando a equação Falhou ao verificar gramática (erro de sintaxe): {\displaystyle P(E, T) = g(E)e^{−E/k_BT}} a qualquer temperatura, sem a necessidade de realizar várias simulações.

Na Fig. 2, mostramos a distribuição canônica resultante na temperatura crítica $T_{c}$, que exibe um único pico.

Figura 2: Distribuição canônica $P(E,T_{c})=g(E)e^{-E/k_{B}T_{c}}$ na temperatura de transição do modelo de Ising 2D com $L=16$ e $k_{B}T_{c}=2.3$.

As distribuições em temperaturas acima e abaixo de $T_{c}$ também apresentam pico único, conforme ilustrado no detalhe da Fig. 3.

Figura 3: Distribuição canônica $P(E,T)=g(E)e^{-E/k_{B}T}$ nas temperaturas $T=2.2$ e $T=2.4$ do modelo de Ising 2D com $L=16$.

Quantidades termodinâmicas

Os resultados simulados, calculados diretamente com as equações descritas no capítulo anterior, e as soluções exatas se sobrepõem quase perfeitamente em uma ampla região de temperatura de $k_{B}T=0-8$.

Testes mais rigorosos de precisão são fornecidos a partir dos erros relativos para as respectivas quantidades termodinâmicas. Os erros relativos são muito pequenos para toda a região de temperatura de $k_{B}T=0-8$.

Como o sistema tem uma transição de fase de segunda ordem, a primeira derivada da energia livre é uma função contínua da temperatura. Não há saltos na energia interna ou na entropia, mesmo no limite, pois o tamanho do sistema vai para o infinito.

Figura 4: Quantidades termodinâmicas do modelo de Ising 2D com $L=16$ calculado a partir da densidade de estados $g(E)$. Na figura, estão mostrados: (a) energia interna, (b) calor específico, (c) energia livre de Helmholtz e (d) entropia. Os gráficos inset são os erros relativos.

Referências

<references>

Código

asdfasdf