7.3 Implementando o Algoritmo de Metropolis

Veremos agora como escrever um código computacional para realizar simulações do modelo de Ising utilizando o algoritmo de Metropolis. Para facilitar continuaremos utilizando o caso com o campo magnético $B=0$, embora a generalização para o caso $B \neq 0$ não seja difícil. Vale lembrar que quase todos os estudos anteriores do modelo de Ising, incluindo a solução exata de Onsager em duas dimensões, foram para sistemas com $B=0$.

Inicialmente, precisamos de uma rede de spins, então definimos um conjunto de $N$ variáveis, uma matriz, que pode assumir valores $\pm 1$. Assim construiremos uma matriz que possuem somente números inteiros contendo valores $+1$ e $-1$. Normalmente, é aplicada condições periódicas de contorno à matriz, isto é, impõe-se que os spins de uma extremidade da rede são vizinhos dos spins da outra ponta da rede, semelhante o teorema de Bloch. Isso garante que todas os spins têm o mesmo número de vizinhos e uma geometria local, e que não há nenhum spin com propriedades diferentes um das dos outros; todos os spins são equivalentes e o sistema é completamente invariante sobre translação. Na prática, isso melhora consideravelmente a qualidade dos resultados da simulação.

Uma variação sobre a ideia de condição periódica de contorno é usar a ``condição de contorno helicoidal'', que é ligeiramente diferente da condição periódica de contorno tradicional, possui todos os mesmos benefícios, é consideravelmente mais simples de implementar e pode tornar o código significativamente mais rápido.

A seguir deveremos decidir em qual temperatura, ou alternativamente, qual o valor de $\beta$ que desejaremos realizar a simulação, e será necessário atribuir um valor inicial para cada spin – que será o estado inicial do sistema. Em muitos casos, o estado inicial que é escolhido não é particularmente importante, embora, às vezes, uma escolha criteriosa pode reduzir o tempo necessário para chegar ao equilíbrio. Os dois estados iniciais mais utilizados são o estados de temperatura zero e o estado de temperatura infinita. Em $T = 0$ o modelo de Ising estará em seu estado fundamental. Quando a energia de interação $J$ é maior que zero e o campo externo $B$ é zero (como serão nos casos das nossas simulações), há na verdade dois estados fundamentais. Estes estados são os casos em que os spins são todos para cima ou todos para baixo. É fácil ver que estes estados devem ser os fundamentais, uma vez que nesses estados cada par de spins contribui para a menor energia possível ($-J$) no primeiro termo do Hamiltoniano da Eq. 1. Em qualquer outro estado, haverá pares de spins que contribuem com energia $+J$ para o hamiltoniano, de modo que seu valor total será maior. Se $B \neq 0$ haverá apenas um estado fundamental, o campo magnético garante que apenas um dos dois estados fundamentais seja favorecido em relação um ao outro. O outro estado inicialcomumente utilizado é o estado $T = \infty$. Quando $T = \infty$ a energia térmica disponível $kT$ para flipar o spin é infinitamente maior do que a energia de interação spin-spin $J$, então os spins são orientados para cima ou para baixo de maneira randômica de forma a ficarem não correlacionados.

As duas opções do estado inicial são populares porque correspondem a uma temperatura conhecida e bem definida e facilmente de serem geradas. Há no entanto, um outro estado inicial, que às vezes pode ser muito útil. Muitas vezes não basta realizar uma simulação em uma única temperatura, para isto é realizado um conjunto de simulações para diferentes valores de $T$, no intuito de investigar o comportamento do modelo em função da variação de temperatura. Neste caso leva-se vantagem escolher como o estado inicial do sistema o estado final, para uma simulação a uma temperatura próxima. Por exemplo, suponha que estamos interessados em investigar uma gama de temperaturas entre $T = 1,0$ e $T = 2,0$ em passos de $0,1$. (Aqui referimos-nos à temperatura em unidades de energia de modo que $k = 1$. Assim, quando dizemos que $T = 2,0$ queremos dizer $\beta^{-1} = 2,0$). Então, poderíamos iniciar a simulação em $T = 1,0$ usando o estado da temperatura zero com todos os spins alinhados como o estado inicial. Ao final da simulação, o sistema estará em equilíbrio a $T = 1,0$, assim poderemos utilizar o estado final da simulação como a entrada do estado inicial da simulação em $T = 1,1$, e assim por diante.

Iniciando a simulação, o primeiro passo é gerar um novo estado. O novo estado deve ser diferente do atual apenas pelo flip de um spin, e cada estado deve ser exatamente tão provável como qualquer outro a ser gerado. Esta é uma tarefa fácil de realizar. Pegaremos um único spin $k$ aleatoriamente na rede para ser flipado. Em seguida calcularemos a diferença de energia $E_{\nu}-E_{\mu}$ entre o estado novo e o antigo, a fim de aplicar a Eq. 7. A maneira simples de realizar o cálculo da energia é calcular $E_{\mu}$ diretamente substituindo os valores de spins ( $s^{\mu}_{i}$) do estado $\mu$ no Hamiltoniano (Eq. 1), então flipar o spin $k$ e calcular $E_{\nu}$, obter a diferença. No entanto, não é a maneira mais eficiente de realizar esta tarefa. Mesmo com $B=0$, é necessário realizar a soma do primeiro termo da Eq. 1, que tem tantos termos quantos os termos de conexões da rede, que é $1/2Nz$. Porém, a maioria desses termos não se alteram quando é realizado apenas um flip de spin. Os únicos termos que se alteram são os que envolvem o flip de spin. Os outros termos permanecem com o mesmo valor e se cancela quando tomamos a diferença $E_{\nu}-E_{\mu}$. A mudança na energia entre dois estados será:

$$E_{\nu }-E_{\mu }=-J\sum _{\langle \mathit{ij}\rangle }s_{i}^{\nu }s_{j}^{\nu }-J\sum _{\langle \mathit{ij}\rangle }s_{i}^{\mu }s_{j}^{\mu }$$

$$E_{\nu }-E_{\mu }=-J\sum _{i\mathit{n.n.}\mathit{para}k}s_{i}^{\mu }(s_{k}^{\nu }-s_{k}^{\mu })$$ (7.3.1)

Na segunda linha a soma é apenas para os spins $i$ que são os primeiros vizinhos do spin flipado $k$ e usamos o fato de que todas esses spins não se invertem, de modo que ( $s^{\nu}_{i} - s^{\mu}_{i}$). Se $s^{\mu}_{k}=+1$, então após o spin $k$ ter sido flipado deveremos ter $s^{\nu}_{k}=-1$, de modo que $s^{\nu}_{k} - s^{\mu}_{k}=-2$. Por outro lado, se $s^{\mu}_{k}=-1$ então $s^{\nu}_{k} - s^{\mu}_{k}=+2$. Assim, podemos escrever

$$s_{k}^{\nu }-s_{k}^{\mu }=-2s_{k}^{\mu }$$ (7.3.2)

e então

$$E_{\nu }-E_{\mu }=2J\sum _{i\mathit{n.n.}\mathit{para}k}s_{i}^{\mu }s_{k}^{\mu }$$

$$E_{\nu }-E_{\mu }=2\mathit{Js}_{k}^{\mu }\sum_{i\mathit{n.n.}\mathit{para}k}s_{i}^{\mu }$$ (7.3.3)

Esta expressão envolve apenas a soma sobre termos $z$, ao invés de $1/2Nz$, e não será necessário realizar multiplicações para os termos na soma, por isso é muito mais eficiente do que avaliar diretamente a variação na energia. Isto envolve apenas os valores dos spins no estado $\mu$, avaliando previamente antes de flipar realmente o spin $k$.

O algoritmo consiste em calcular $E_{\nu}-E_{\mu}$ da Eq. 10 e em seguida pela regra da Eq. 7: se $E_{\nu }-E_{\mu} \leq 0$ aceita-se a transição e flipando $s_k \rightarrow -s_k$. Se $E_{\nu }-E_{\mu} > 0$ poderíamos, ou não, flipar o spin. Com o algoritmo de Metropolis o flip do spin será dado pela probabilidade $A(\nu \rightarrow \mu) = e^{-\beta(E_{\nu}-E_{\mu})}$. Podemos fazer isso da seguinte forma: avaliamos a relação de aceitação $A(\nu \rightarrow \mu)$, usando o valor de $E_{\nu}-E_{\mu}$ da Eq. 10 , então escolhe-se um número aleatório $r$ entre zero e um. A rigor, o número pode ser igual a zero, mas deve ser inferior a um ( $0 \leq r < 1$). Se esse número for menor do que a nossa relação de aceitação, $r < A(\nu \rightarrow \mu)$, então o spin é flipado, caso contrário deixa-se como está.

Essa é sequência completa do algoritmo. Agora basta repetir os mesmos cálculos, escolhendo um flip e calculando a variação da energia, então decidindo se ocorrerá ou não o flip de acordo com a Eq. 7. Na verdade, há um outro truque que pode deixar o algoritmo um pouco mais rápido. Uma das partes mais lentas do algoritmo é o cálculo da exponencial, pois calcula-se a energia do novo estado, faz-se a comparação de energia e só então flipa ou não o spin. O cálculo de exponenciais em um computador é feito geralmente utilizando uma aproximação polinomial, que envolve a execução multiplicações de ponto flutuante e somatórios, o que pode levar um intervalo de tempo considerável. É possível contornar este esforço e tornar o código computacional mais rápido ao verificar que a quantidade da Eq. 10 que esta sendo calculada para uma pequena quantidade de valores. Cada termo da soma só pode assumir valores $+1$ e $-1$. Então, todo o somatório, que tem z termos, só pode ter valores $-z$, $-z+2$,$-z+4$ e assim por diante até $+z$, um total de $z+1$ valores possíveis. Só precisamos então calcular o exponencial quando a soma for negativa (Eq. 7), de fato temos apenas $1/2z$ valores de $E_{\nu}-E_{\mu}$ no qual é necessário o cálculo das exponenciais. Assim, utilizando o bom senso, faz todo o sentido calcular os valores dessas $1/2z$ exponenciais antes de iniciar o cálculo, e armazená-los na memória do computador, em que pode ser feito uma pesquisa durante o andamento da simulação. O cálculo é realizado apenas uma vez no início não sendo necessário calcular novamente a exponencial durante a simulação. O único cálculo necessário de ponto flutuante será na geração do número aleatório $r$, todos os outros cálculos envolvem apenas números inteiros, que é mais rápido que o tratamento com números reais.