Gás de Rede 2D

Este trabalho tem o intuito de apresentar uma descrição do modelo do Gás de Rede bidimensional e traçar paralelos com o modelo de Ising, discutindo a relação entre seus Hamiltonianos.

Gás de Rede

O Modelo do Gás de Rede 2D consiste em um sistema de $N$ partículas da forma $\sigma ={\sigma _{1},\sigma _{1},\dotsc ,\sigma _{N}}$ onde cada sítio da rede pode assumir o valor $1$, ocupado por uma partícula, ou $0$, não ocupado por uma partícula. A energia total do sistema é dada pelo Hamiltoniano do Gás de Rede, descrito pela equação

${\mathcal {H}}=-\epsilon \sum _{\langle i,j\rangle }\sigma _{i}\sigma _{j}\qquad (1)$,

onde o somatório é dado entre os quatro vizinhos mais próximos e $\epsilon$ é a constante de interação entre as partículas, para $\epsilon \geq 0$ a interação é atrativa. Por se tratar de uma rede quadrada com $L^{2}=L\times L$ sítios, apenas uma parcela da rede é ocupada por partículas, ou seja, possuímos uma densidade constante $\rho$ de partículas. Podemos expressar a condição da densidade constante da forma

$\sum _{i}\sigma _{i}=\rho L^{2}$

Fazendo uma mudança de variáveis da forma $s_{i}=2\sigma _{i}-1$ saímos da situação de ocupação e não ocupação de sítios e obtemos variáveis do Modelo de Ising , spins Up e Down. A variável $s_{i}$ assume valor $+1$ (up) quando o sítio esta ocupado por uma partícula e $-1$ quando não está. Aplicando a mudança de variáveis no Hamiltoniano do Gás de Rede obtemos

${\mathcal {H}}=-{\frac {1}{4}}\epsilon \sum _{\langle i,j\rangle }(s_{i}+1)(s_{j}+1)$

$=-{\frac {1}{4}}\epsilon \sum _{\langle i,j\rangle }[s_{i}s_{j}+s_{i}+s_{j}+1]$

$=-{\frac {1}{4}}\epsilon \sum _{\langle i,j\rangle }s_{i}s_{j}-{\frac {1}{4}}\epsilon \sum _{\langle i,j\rangle }s_{i}-{\frac {1}{4}}\epsilon \sum _{\langle i,j\rangle }s_{j}-{\frac {1}{4}}\epsilon \sum _{\langle i,j\rangle }1$.

Abrindo os termos dos três últimos somatórios se encontra que

$\sum _{\langle i,j\rangle }s_{i}=\sum _{\langle i,j\rangle }s_{j}=z\sum _{i}s_{i}$, e

$\sum _{\langle i,j\rangle }1=zL^{2}$,

o que simplifica a expressão para

${\mathcal {H}}=-{\frac {1}{4}}\epsilon \sum _{\langle i,j\rangle }s_{i}s_{j}-{\frac {1}{2}}\epsilon z\sum _{i}s_{i}-{\frac {1}{4}}\epsilon zL^{2}\qquad (2)$,

onde $z$ é o número de coordenação, isto é, o número de vizinhos de cada sítio, que para a rede quadrada é 4. Note a mudança no índice do somatório de $s_{i}$. Vemos que este Hamiltoniano possui a forma do Hamiltoniano do modelo de Ising com aplicação de um campo magnético externo, a menos de uma constante. Entretanto, com a densidade constante, fazendo a mesma mudança de variáveis, obtemos

$\sum _{i}s_{i}=(2\rho -1)L^{2}\qquad (3)$.

Aplicando essa relação no Hamiltoniano obtemos

${\mathcal {H}}=-{\frac {1}{4}}\epsilon \sum _{\langle i,j\rangle }s_{i}s_{j}-z\epsilon L^{2}\rho$.

Como $z$, $\rho$ e $L^{2}$ são constantes, o segundo termo é constante. Definindo $J=\epsilon /4$ o Hamiltoniano se torna

${\mathcal {H}}=-J\sum _{\langle i,j\rangle }s_{i}s_{j}+{\text{constante}}$.

A constante não influencia nos valores da medida, pois ela se cancela na equação de mudança de estado . Vemos que, ao assumir a densidade constante, o modelo do Gás de Rede se torna o modelo de Ising sem a presença de um campo magnético. A condição da densidade constante é a magnetização total do modelo.

É interessante notar que podemos tratar o modelo do gás de rede de duas formas: (a) assumir a densidade constante e utilizar o primeiro Hamiltoniano apresentado, ou (b) não assumir a densidade constante, aplicar a mudança de variáveis discutida e utilizar o Hamiltoniano que possuí a forma de Ising com campo magnético.

Método de Monte Carlo

Princípios gerais

Para desenvolver a dinâmica do gás de rede, isto é, fazer com que ocorram alterações na posição das partículas, recorre-se ao método de Monte Carlo. Trata-se, em linhas gerais, de um método estatístico baseado em amostragem aleatória. É comum, ao introduzir-se o assunto, utilizar problemas de cálculo de áreas, como, por exemplo, a clássica história das crianças que desenharam, na areia, um círculo contido em um quadrado, de modo que o diâmetro do círculo seja igual ao lado do quadrado. Jogam, aleatoriamente, pedras no quadrado. Após uma grande quantidade de arremessos, compara-se o valor de pedras que estão dentro do círculo com o número total dentro do quadrado, fazendo a razão entre esses dois valores, encontra-se uma aproximação para a razão entre a área do círculo e a do quadrado, que equivale a $\pi /4$. Desse modo, pode ser encontrada uma aproximação do valor de $\pi$ através do método de Monte Carlo, que simula uma aproximação experimental. É possível propor um novo modelo, no qual o quadrado é grande e as pedras são atiradas aleatoriamente de dentro do quadrado, de modo que a posição de arremesso é atualizada para a posição da última pedra jogada. A principal alteração nesse método é a possibilidade de que uma pedra seja atirada para fora do quadrado. Nesse caso, a posição para o próximo arremesso não muda e as pedras que foram parar fora das bordas do quadrado são empilhadas naquela posição. Essas mudanças no método caracterizam a aplicação de Monte Carlo em uma cadeia de Markov, definida de modo que a configuração do próximo estado depende somente da configuração atual. Até aqui, vimos o método de Monte Carlo aplicado em situações nas quais a probabilidade de um evento ocorrer na configuração subsequente do sistema é zero ou um. Para novas distribuições de probabilidade, utiliza-se o algoritmo de Metrópolis.

Algoritmo de Metrópolis

Como mencionado, o método de Monte Carlo funciona através de atualizações em um sistema, fazendo com que ele passe de um estado para outro. Consideremos $\mu$ como sendo o estado atual do sistema. O próximo estado será chamado $\nu$. Chamamos a probabilidade do sistema passar do estado $\mu$ para o estado $\nu$ de $p(\mu \to \nu )$. A ideia do algoritmo de Metrópolis é aplicar essa probabilidade no método de Monte Carlo. Para isso, um método muito utilizado é calcular a probabilidade $p(\mu \to \nu )$ e compará-la a um número aleatório $r$, de modo que $0<r<1$. Se $p(\mu \to \nu )>r$, então o estado $\nu$ é adotado pelo sistema. Para aplicar essas ideias ao problema do gás de rede, o método mais simples é utilizando o algoritmo de Kawasaki, que faz com que dois sítios, selecionados aleatoriamente, tenham seus valores trocados. Isso implica que a magnetização do sistema será constante. Para calcular a probabilidade de aceitação dessa mudança, calcula-se a variação de energia entre os estados,$\Delta E=E_{\nu }-E_{\mu }$, e aplica-se o valor obtido na equação de probabilidade de aceitação de Metrópolis:

$p(\mu \to \nu )=\left\{{\begin{array}{rc}e^{\frac {-\Delta E}{k_{B}T}},\quad {\text{se}}\quad \Delta E>0\\1,\quad {\text{caso contrario}}\end{array}}\right.$

Para este trabalho, inicialmente foram feitas algumas alterações no algoritmo de Kawasaki. Limita-se a troca de valores somente entre um sítio e seus quatro vizinhos principais (isto é, acima, abaixo, à direita e à esquerda). Isso porque se tinha a intenção de recriar o comportamento difusivo de um gás no modelo do gás de rede em duas dimensões, porém, essa escolha para a dinâmica das simulações limita em muito a facilidade do sistema equilibrar, como será mostrado na parte de resultados. O algoritmo inicial consiste na seguinte ideia: No estado $\mu$, seleciona-se aleatoriamente um sítio $i$. A partir dele, é sorteado um de seus vizinhos, digamos, $j$, com o qual haverá a troca de valor. Antes que ocorra a troca, são calculados os valores de dois Hamiltonianos: o que relaciona o sítio $i$ com seus vizinhos e o que relaciona $j$ com seus vizinhos. Somados, eles representam o valor de $E_{\mu }$. Após a troca, esses valores são recalculados, formando o valor de $E_{\nu }$. Com esses valores, segue-se o cálculo da probabilidade da forma usual.

Após verificar que a dinâmica inicialmente escolhida não era apropriada, optou-se por usar o algoritmo de Kawasaki tradicional, fazendo uso de uma listagem de quais eram os estados dos sítios e escolhendo os pares a serem trocados sempre como sítios de estados diferentes. Isso faz com que a aceitação do algoritmo de Monte Carlo seja muito mais elevada que ao do nosso algoritmo inicial, pois no caso em que grandes clusters são formados as trocas eram sempre rejeitadas de acordo com a dinâmica inicial.

Considerações práticas

Já consideramos os aspectos iniciais do método de Monte Carlo, mostrando suas aplicações. Vimos, do mesmo modo, a utilização do algoritmo de Metrópolis, dentro do qual aplicamos o algoritmo de Kawasaki (modificado) para desenvolver o sistema do gás de rede. Resta elucidar a função do parâmetro temporal, isto é, como as iterações desses algoritmos influenciam o sistema. Primeiramente, deve-se fixar um valor para o número de passos temporais que serão realizados. Em cada passo temporal, são realizadas $L^{2}$ iterações do algoritmo de Metrópolis, abordado acima, sendo $L$ o lado da rede quadrada. Chamamos de Passo de Monte Carlo (MCS) cada um desses passos temporais. Em outras palavras: em cada passo de Monte Carlo, todos os sítios tem a possibilidade de realizar uma troca.

Resultados

Eficiência da dinâmica local

Como mencionado acima, a primeira dinâmica escolhida foi a que restringe os sítios a trocarem de estado para apenas aqueles que são vizinhos próximos (os 4 acima, abaixo, a direita e a esquerda). Se verificou ao analisar os dados das simulações que na verdade essa dinâmica não era apropriada, pois a aceitação do algoritmo é muito baixa e com a amostragem escolhida ($10^{6}$ MCS) o sistema ainda não estaria equilibrado. Isso pode ser percebido mais facilmente observando os histogramas das simulações. Sabe-se que a distribuição das grandezas macroscópicas de um sistema termodinâmico é a distribuição de Boltzmann , então se espera para o sistema analisado que a energia tenha uma distribuição normal de frequências. Observe a diferença entre as distribuições para a energia do sistema no caso da dinâmica local inicialmente escolhida e no caso do algoritmo tradicional de Kawasaki:

Diferença entre distribuições para a energia do sistema por dinâmica escolhida.
Distribuição da energia para a dinâmica local após $10^{6}$ MCS.	Distribuição da energia para a dinâmica do algoritmo de Kawasaki após $10^{6}$ MCS.
Distribuição da energia para a dinâmica local em alta temperatura (acima da temperatura crítica) após $10^{6}$ MCS.	Distribuição da energia para a dinâmica do algoritmo de Kawasaki em alta temperatura (acima da temperatura crítica) após $10^{6}$ MCS.

Assim podemos notar que o sistema não está equilibrado após a amostragem com a dinâmica local. No caso do algoritmo tradicional o sistema se encontra em equilíbrio. Para temperaturas mais altas, acima da temperatura crítica para cada modelo onde não há mais separação de fase, a dinâmica tem eficiência satisfatória, pois a aceitação do algoritmo de Monte Carlo é maior já que sem a formação de grandes clusters a probabilidade de uma partícula ser selecionada e poder trocar para um sítio vazio que seja seu vizinho é muito maior. Contudo, ainda assim é mais vantajoso executar as simulações na linguagem de Ising, já que a aceitação do algoritmo ainda é maior do que para a dinâmica local proposta, então a quantidade de dados descorrelacionados para a mesma amostragem de $10^{6}$ MCS é maior. Dessa maneira, para uma amostragem de mesmo tamanho, o valor estimado para uma grandeza $\langle Q\rangle$ é mais preciso no caso da dinâmica proposta pelo algoritmo de Kawasaki, e para uma mesma precisão no valor estimado seria necessária uma amostragem menor do que a da dinâmica local. A referência na qual esse trabalho mais se inspira, Monte Carlo Methods in Statistical Physics, de Newman e Barkema, não dá a devida ênfase para quão interessante computacionalmente é se manter na linguagem do modelo de Ising e adotar o algoritmo de Kawasaki. É importante notar que as simulações da dinâmica de Kawasaki foram feitas na linguagem do modelo de Ising, com spins $+1$ e $-1$, ao invés de sítios ocupados por partículas $1$ e sítios desocupados $0$. Essa mudança é importante para entender a mudança na escala das grandezas entre as simulações. Como a constante de interação $J$ do modelo de Ising é um quarto da constante $\epsilon$, a temperatura e a energia no modelo de Ising valem ${\frac {1}{4}}$ da temperatura e da energia nas simulações com a dinâmica local na linguagem do gás de rede.

Separação de fases

Animações ilustrando as duas fases acima e abaixo de $T_{c}=0.56725$.
Separação de fases do sistema do gás de rede em fases de densidade alta e baixa. Note como ocorre a formação de clusters e estes tendem a se agrupar formando clusters ainda maiores.	Fase gasosa do modelo do gás de rede com densidade constante. Note como os clusters se agrupam e se desfazem rapidamente.

Um dos aspectos interessantes do modelo do gás de rede é o fato de ser um modelo chamado de COP-Ising, ou Ising com parâmetro de ordem conservado. O fato de a densidade ser constante no modelo do gás de rede implica que o modelo de Ising análogo tem a magnetização constante (chamado de parâmetro de ordem). Abaixo da temperatura crítica $T_{c}=0.56725$ o sistema se separa em duas fases com densidades diferentes, se acomodando na forma de clusters cada vez maiores conforme a temperatura diminui. Acima de $T_{c}$ o sistema assume uma forma "gasosa", em que a formação de pequenos clusters é apenas momentânea e a densidade sobre a rede é aproximadamente homogênea.

Energia e efeito do tamanho da rede

A grandeza que foi medida nas simulações é a energia total do sistema. Para analisar os efeitos do tamanho da rede devemos comparar os valores das energias por número de sítios, já que quanto maior a rede maior o número de partículas e portanto maior o valor absoluto da energia. Abaixo temos um gráfico com os valores da energia total do sistema por número de sítios em função da temperatura, com os dados das três redes simuladas (L = 16, 32, 64).

Energia total do sistema dividida pelo número de sítios em função da temperatura para redes de tamanhos diferentes. Cada ponto no gráfico é a média de 5 amostragens com seu desvio padrão.

Aqui o comportamento da energia é como esperado para o modelo. Há uma queda no valor da energia conforme a temperatura do reservatório diminui abaixo da temperatura crítica, e conforme a temperatura aumenta acima da temperatura crítica, a energia aumenta mais devagar dado que é o estado em que o sistema está na sua fase amorfa.

Algo que se nota no gráfico acima é a diferença nos valores da energia das redes $L=32$ e $L=64$ e a rede de $L=16$. É perceptível que as duas redes maiores apresentam valores de energia muito mais próximos. Isso se dá justamente pelo tamanho da rede e o fato da densidade ser a mesma.

Ilustração dos estados de menor energia
Diagrama ilustrando o estado de menor energia para a rede de $L=16$. As cores indicam o número de sítios vizinhos ocupados.	Animação ilustrando o número de partículas internas ao cluster se tornando mais relevante que o número de partículas nas bordas conforme a rede aumenta.

Para entender como o valor da energia total do sistema por número de sítios varia com o tamanho da rede é interessante calcular essa grandeza no estado de menor energia. No caso da rede quadrada, o estado de menor energia é um cluster quadrado. Como ilustrado no diagrama acima, existem sítios com 3 vizinhos ocupados e sítios com 2 vizinhos ocupados. Estes são os da borda do cluster. Todos os outros internos ao cluster tem 4 vizinhos ocupados e portanto contribuem com $-4\epsilon$ para o Hamiltoniano.

Nesse trabalho adotamos $\rho =0.25$, o que facilita o cálculo da energia total, pois o tamanho do lado do cluster é $l={\sqrt {\frac {L^{2}}{4}}}$. Temos $(l-2)^{2}$ partículas internas ao cluster e portanto com 4 vizinhos e $4(l-2)$ partículas nas bordas com 3 vizinhos. O Hamiltoniano é portanto

${\mathcal {H}}_{GR}=(-2\epsilon )4+(-3\epsilon )4(l-2)+(-4\epsilon )(l-2)^{2}$.

O primeiro termo é a referente às interações dos 4 cantos do cluster, que tem apenas 2 vizinhos ocupados. Abrindo os termos e simplificando a equação chegamos ao valor do Hamiltoniano do estado de menor energia do gás de rede:

${\mathcal {H}}_{GR}=-(L^{2}-2L)\epsilon$.

Dividindo por 2 para encontrar o valor da energia total e então por $L^{2}$ chegamos em

$e={\frac {E}{L^{2}}}=-{\frac {\epsilon }{2}}+{\frac {\epsilon }{L}}$,

que é a grandeza apresentada no gráfico anterior. Essa elaboração foi feita para o modelo do gás de rede, porém as simulações foram feitas com o modelo de Ising, com a restrição de que a magnetização era constante (a quantidade de spins $+1$ e $-1$ é constante).

Usando a equação $(2)$ e $J={\frac {\epsilon }{4}}$ chegamos em

${\mathcal {H}}_{ising}=-J\sum _{\langle i,j\rangle }s_{i}s_{j}-8J\sum {i}s_{i}-4JL^{2}$.

Utilizamos então a equação $(3)$ e $\rho =0.25$ e percebemos que os termos do meio e da direita se cancelam, deixando somente o primeiro termo para avaliarmos em termos do estado de menor energia.

${\mathcal {H}}_{ising}=-J\sum _{\langle i,j\rangle }s_{i}s_{j}-8J(2\rho -1)L^{2}-4JL^{2}$ $=-J\sum _{\langle i,j\rangle }s_{i}s_{j}$

Agora podemos avaliar o Hamiltoniano do estado de menor energia do modelo de Ising nessas condições. Para isso podemos nos basear no diagrama usado anteriormente, com a ressalva de que os sítios imediatamente fora do cluster também interagem com ele, já que temos agora spins $-1$ ao invés de sítios vazios. Temos o mesmo número de spins internos ao cluster que de pariculas anteriormente, assim como para as bordas do cluster e os mesmos quatro cantos. Além disso, os sítios fora do cluster também tem contribuição para o Hamiltoniano. O número de spins fora do cluster com 4 vizinhos $-1$ é $L^{2}-{\frac {L^{2}}{4}}-2L$ e o número de spins fora do cluster com 1 vizinho $+1$ é $2L$. Agora podemos somar todas essas contribuições.

${\mathcal {H}}_{ising}=(-4J)(L^{2}-{\frac {L^{2}}{4}}-2L)+(-2J)(2L)+(-2J)4({\frac {L}{2}}-2)+(-4J)({\frac {L}{2}}-2)^{2}$

Da esquerda para a direita os termos são a contribuição dos spins fora do cluster com 4 vizinhos $-1$, dos spins fora do cluster mas com um vizinho no cluster, dos spins na borda do cluster e dos spins interiores ao cluster. Abrindo os termos e simplificando a equação chegamos finalmente em

${\mathcal {H}}_{ising}=-(4L^{2}-8L)J$

E como nas simulações foram utilizadas que $\epsilon =1$ para o caso do gás de rede e $J=1$ para o modelo de Ising, temos portanto que as energias medidas com o modelo de Ising são equivalentes a 4 vezes as energias no modelo do gás de rede, como se verificou ao implementar a mesma função de medida para a linguagem de spins e de sítios vazios e ocupados. E com essa equação para o Hamiltoniano do estado de menor energia do modelo de Ising nessas condições podemos perceber que o efeito do tamanho de rede é praticamente o mesmo, com a ressalva de que a energia por sítio deve ser também 4 vezes a energia por sítio do modelo do gás de rede.

Conclusão

Passamos pelo que vem a ser o modelo do gás de rede na rede quadrada e apresentamos sua relação com o modelo de Ising. O Hamiltoniano do modelo do gás de rede nos mostra que ele é idêntico ao modelo de Ising com um termo de campo magnético externo. Ao fixarmos a densidade do modelo do gás de rede, impomos a restrição de que o modelo de Ising deve ter a magnetização constante para mantermos a equivalência. Também se discutiu o efeito do tamanho das redes simuladas e vimos que quanto maior o tamanho da rede mais próximo será o valor das energias por sítio em comparação com redes de tamanhos similares e vimos que isso se deve pelo fato do número de sítios no cluster do estado de menor energia crescer mais rápido no interior do cluster do que em seu perímetro.

Programas Utilizados

Programas na linguagem C. Para utilizar os programas, abra o terminal e compile da forma

$ gcc prog.c -lm

Onde prog.c é o programa que deseja utilizar. E execute da seguinte maneira

$ ./a.out TEMP SEED

onde o segundo termo é a temperatura do banho térmico e o terceiro é a semente utilizada pelo gerador de números pseudo-aleatórios. Os programas possuem uma diretiva de compilação para visualização do sistema ao decorrer da execução. Para utilizar é necessário ter o gnuplot instalado e compilar da forma

$ gcc -DGNU prog.c -lm

e então executar da maneira

$ ./a.out TEMP SEED | gnuplot

Entretanto com a visualização no gnuplot o programa pode demorar mais para executar, então é recomendado diminuir os tempos de transiente (TRAN) e medidas (TMAX).

Gás de Rede com Densidade Constante e dinâmica local

Modelo de ising com magnetização constante e algoritmo de Kawasaki

Biblioteca com funções em C para simulações de Monte Carlo