7.2 Algoritmo de Metropolis

O primeiro algoritmo de Monte Carlo que veremos é o mais famoso e amplamente utilizado algoritmo de todos eles, o algoritmo de Metropolis, que foi introduzido por Nicolas Metropolis e seus colaboradores em um artigo de 1953 em simulações de gases de esferas rígidas (Metropolis et al. 1953). Utilizarremos esse algoritmo para ilustrar alguns conceitos gerais envolvidos em um cálculo de Monte Carlo, incluindo o equilíbrio, medida do valor esperado e o cálculo de erros. Entretanto, vamos detalhar o algoritmo e estudar como implementá-lo em um código computacional.

O algoritmo de Metropolis segue o seguinte esquema: escolhemos um conjunto de probabilidades de seleção $g(\mu \rightarrow \nu)$, uma para cada possível transição de um estado para outro, $\mu \rightarrow \nu$, então escolhemos um conjunto de probabilidades de aceitação $A(\mu \rightarrow \nu)$ tal que

$$\frac{P(\mu \rightarrow \nu )}{P(\nu \rightarrow \mu )}=\frac{g(\... ...w \nu )A(\mu \rightarrow \nu )}{g(\nu \rightarrow \mu )A(\nu \rightarrow \mu )}$$ (7.2.1)

em que

$\frac{A(\mu \rightarrow \nu )}{A(\nu \rightarrow \mu )}$ pode assumir qualquer valors entre 0 e $\infty$;

e

$g(\mu \rightarrow \nu)$ e $g(\nu \rightarrow \mu )$ pode assumir qualquer valor desejado.

A Eq. 2 satisfaz a condição de balanço detalhado dado pela equação

$$\frac{P(\mu \rightarrow \nu )}{P(\nu \rightarrow \mu )}=\frac{p_{\nu }}{p_{\mu }}=e^{-\beta (E_{\nu }-E_{\mu })}$$ (7.2.2)

O algoritmo trabalha repetindo a escolha de um novo estado $\nu$ aleatoriamente, aceitando ou rejeitando esse novo estado de acordo com a escolha de aceitação probabilística. Se o estado é aceito, o estado que antes era $\mu$ passa agora a ser o estado $\mu$, caso contrário permanece como está. Assim o processo é repetido por várias vezes.

As probabilidades de seleção $g(\mu \rightarrow \nu)$ devem ser escolhidas de modo que a condição de ergodicidade (exigência de que cada estado estará acessível a partir de todos os outros em um número finito de passos) seja cumprida (A condição de ergodicidade é a exigência que deve ser possível para os processos de Markov alcançar qualquer estado do sistema a partir de qualquer outro estado, se executado em um tempo suficiente longo. Isso é necessário para alcançar o objetivo de gerar estados com corretas probabilidades Boltzmann. Cada estado aparece $\nu$ com alguma probabilidade $p_{\nu}$, diferente de zero, na distribuição de Boltzmann, e se esse estado estiver inacessível a partir de outro estado $\mu$, não importando quanto tempo ainda continuamos o processo para, em seguida, iniciarmos do estado $\mu$: a probabilidade de encontrar $\nu$ nos estados das cadeias de Markov será zero, e não $p_{\nu}$ A condição de ergodicidade nos diz que estamos autorizados a fazer algumas transições de probabilidades a partir do processo de Markov zero, mas deve haver pelo menos um caminho diferente de zero nas transições de probabilidades entre dois estados escolhidos. Na prática, a maioria dos algoritmos de Monte Carlo define quase todas as probabilidades de transição como sendo zero, e devemos ter cuidado para que ao fazê-lo não criarmos um algoritmo que viola a ergodicidade. Para os algoritmos desenvolvidos deve-se provar explícitamente que ergodicidade está sendo satisfeita antes de uso do algoritmo na producao de resultados serios). Isto ainda nos deixa uma grande margem sobre a forma de como a passagem de um estado para outro será escolhido e dado um estado inicial $\mu$ podemos gerar qualquer número de estados candidatos $\nu$ simplesmente girando diferentes grupos de spins da rede. As energias dos sistemas em equilíbrio térmico permanecem dentro de um intervalo muito estreito de energia, as flutuações de energia são pequenas em comparação com a energia de todo o sistema. Em outras palavras, o sistema real passa a maior parte de seu tempo em um subconjunto de estados com uma estreita faixa de energias e raramente faz transições que mudam a energia do sistema de forma drástica. Isso diz que provavelmente não queremos gastar muito tempo em nossa simulação considerando as transições para os estados cuja energia é muito diferente da energia do estado atual. A forma mais simples de conseguir isto no modelo de Ising é considerar apenas os estados que diferem de um dos presentes pelo giro de um único spin. Um algoritmo que realiza este tipo de giro de spin é dito como tendo um dinâmica single-spin-flip. O algoritmo que descreveremos tem dinâmica single-spin-flip, embora isto não seja realizado no algoritmo de Metropolis. Como discutido abaixo, é a escolha particular da relação de aceitação que caracteriza o algoritmo de Metropolis. Nosso algoritmo continuaria a ser um algoritmo de Metropolis, mesmo que girasse várias vezes os vários spins de uma só vez.

Usando a dinâmica single-spin-flip garante que o novo estado $\nu$ vai ter uma energia $E_{\nu}$ diferenciando-se da energia corrente $E_{\mu}$ de até, no máximo, $2J$ para cada ligação entre o spin flipado e seus vizinhos. Por exemplo, em uma rede quadrada em duas dimensões cada spin tem quatro vizinhos, assim a máxima diferença de energia seria $8J$. A expressão geral $2zJ$, onde $z$ é o número de coordenação da rede, ou seja, o número de vizinhos em que cada ponto da rede possui (isto não é a mesma coisa que o ``o número de coordenação de spin''. O número de coordenação de spin é o número de spins $j$ na vizinhança de $i$ que possui o mesmo valor do spin $i$: $s_i = s_j$). Usando uma única dinâmica single-spin-flip também garante que o algoritmo obedece a ergodicidade, desde que esteja claro que seja possível pegar um estado a partir do outro na rede finita flipando os spins um por um, onde há diferença.

No algoritmo de Metropolis as probabilidades de seleção $g(\mu \rightarrow \nu)$ para cada um dos possíveis estados $\mu$ são escolhidos para serem iguais. As probabilidades de seleção de todos os outros estados são definidos como sendo zero. Suponha que há $N$ spins no sistema que desejamos simular. Com uma única dinâmica single-spin-flip teremos então $N$ diferentes spins que poderão ser flipados e, portanto, $N$ estados $\nu$ que podem alcançar um determinado estado $\mu$. Assim, há $N$ probabilidades de seleção $g(\mu \rightarrow \nu)$ que serão diferentes de zero, e cada um deles assumirá o valor

$$g(\mu \rightarrow \nu )=\frac{1}{N}$$

Com estas probabilidades de seleção, a condição de balanço detalhado, a Eq. 3, assume a forma

$$\frac{P(\mu \rightarrow \nu )}{P(\nu \rightarrow \mu )}=\frac{g(\... ...(\mu \rightarrow \nu )}{A(\nu \rightarrow \mu )}=e^{-\beta (E_{\nu }-E_{\mu })}$$

(3a)

Agora temos que escolher as taxas de aceitação $A(\mu \rightarrow \nu)$ para satisfazer a equação. Uma possibilidade seria escolher

$$A(\mu \rightarrow \nu )=A_{0}e^{\frac{-{1}}{2}\beta (E_{\nu }-E_{\mu})}$$ (7.2.3)

A constante de proporcionalidade $A_0$ não está presente na Eq. (4), podemos escolher qualquer valor para esta constante, exceto que $A(\mu \rightarrow \nu)$, sendo uma probabilidade, nunca deve se tornar maior do que $A_0$. A maior diferença de energia $E_{\nu}-E_{\mu}$ que podemos ter entre dois estados é $2zJ$, onde $z$ é o número de coordenação da rede. Isso significa que o maior valor de $e^{\frac{-{1}}{2}\beta (E_{\nu }-E_{\mu })}$ é $e^{\beta \mathit{zJ}}$. Assim, a fim de certificar-se de que $A(\mu \rightarrow \nu) \leq 1$, escolhemos

$$A_{0}\le e^{-\beta \mathit{zJ}}$$ (7.2.4)

Para tornar o algoritmo mais eficiente possível, fazemos as probabilidades de aceitação tão grande quanto possível, de modo que $A_0$ seja tão grande quanto é permitido ser, o que nos dá

$$A(\mu \rightarrow \nu )=A_{0}e^{\frac{-{1}}{2}\beta (E_{\nu }-E_{\mu }+2\mathit{zJ})}$$ (7.2.5)

Este ainda não é o algoritmo de Metropolis, mas utilizando essa probabilidade de aceitação, podemos realizar uma simulação de Monte Carlo do modelo de Ising, sendo uma amostragem correta da distribuição de Boltzmann. No entanto, a simulação será muito ineficiente, pois a relação de aceitação, a Eq. (6), será muito pequena para quase todos os movimentos.

Na Eq. 4, foi assumido uma forma funcional particular para o índice de aceitação, mas a condição do balanço detalhado, Eq. 3a, na verdade não requer que seja dessa forma. A Eq. 3a especifica apenas a relação entre pares de probabilidades de aceitação, deixando muitas possibilidades de manobra. Quando há uma restrição do tipo da Eq. 3a a maneira de maximizar as taxas de aceitação (e portanto, produzir um algoritmo mais eficiente) é sempre dar para o maior dos dois índices, o maior valor possível, ou seja, 1, e em seguida, ajustar o outro para satisfazer a restrição. Para ver como isso funciona, suponha que os dois estados $\mu$ e $\nu$, $\mu$ tenha a menor energia e $\nu$ a maior energia: $E_{\nu} < E_{mu}$. Assim a maior das duas taxas de aceitação $A(\nu \rightarrow \mu)$, então ajustamos o conjunto para ser igual a $1$. A fim de satisfazer a Eq. 3a, $A(\mu \rightarrow \nu)$ deve assumir o valor $e^{-\beta (E_{\nu }-E_{\mu })}$. Assim, o algoritmo otimizado será:

$$A(\mu \rightarrow \nu )=\genfrac{}{}{0pt}{0}{e^{-\beta (E_{\nu }-... ..._{\nu }-E_{\mu }>0}{\rightarrow \mathit{para}\;\mathit{outros}\;\mathit{casos}}$$ (7.2.6)

Em outras palavras, se selecionar um novo estado que tem uma energia inferior ou igual ao atual, sempre devemos aceitar a transição para esse estado. Se o estado tiver uma energia mais elevada, talvez pode ser aceita com a probabilidade dada acima. Este é o algoritmo de Metropolis para o modelo de Ising com uma única dinâmica single-spin-flip. Esta é a parte pioneira que abriu caminho para Metropolis e colaboradores em seu artigo gases de esferas rígidas, podendo ser aplicado a qualquer modelo de acordo com a Eq. 7.