Grupo - Conservação do Parâmetro de Ordem

Introdução

O modelo de Ising possui características universais que permitem aplicá-lo a situações diversas sendo tão versátil a ponto de descrever desde ferromagnetos até interações sociais. Dentro dessa gama de possibilidades existe o modelo de Conservação do Parâmetro de Ordem (CPO) em que, como o nome indica, mantém-se o parâmetro de ordem constante. No caso de um ferromagneto o parâmetro de ordem é a magnetização, portanto, um modelo de Ising sujeito CPO a grandeza análoga à magnetização se manteria constante a cada passo da simulação.

Apesar da estrutura matemática muito similar ao modelo de Ising, o modelo de CPO com sua simples condição de conservação do parâmetro de ordem aliado a condições de contorno permite que se modele sistemas marcadamente diferentes do tradicional sistema de ferromagneto tais como o gás de rede onde é possível estudar o comportamento de interfaces vapor-sólido ou vapor-líquido em condições de equilíbrio como por exemplo o equilíbrio entre água líquida e seu vapor ou entre gelo e vapor d'água.

O gás de rede é um modelo simplificado de um gás real onde se associa a cada ponto da rede uma partícula (átomo) ou sua ausência (vacância). Ao contrário do gás real a coordenada do movimento não é contínua, pois as partículas se movem de maneira discreta somente pelos vértices da rede. Pode-se refinar o modelo de diversas formas:

• Atribuindo inércia às partículas
• Alterando a forma da rede (quadrada, hexagonal, fcc, bcc, cúbica)
• Incluindo partículas de tipos diferentes com interações comum a seu respectivo tipo
• Presença e/ou tipos de colisões

No entanto, uma versão simplificada (e simples de simular) desse modelo é suficiente para reproduzir qualitativamente o comportamento de interfaces.

Teoria

No modelo simplificado do gás de rede as partículas (sem inércia), movem-se de forma aleatória sob excitação térmica e satisfazem as seguintes condições:

1. O número total de partículas é fixo: nenhuma partícula deixa ou entra no sistema, portanto, caso desapareça a partícula deve reaparecer em outro ponto da rede no mesmo passo de simulação.
2. Um ponto da rede pode ser ocupado por uma única partícula ou permanecer vazio (não ocupado). Essa é uma maneira grosseira de assimilar o caráter físico de repulsão do gás real onde partículas não podem interpenetrar-se devido a exclusão de Pauli.
3. Se duas partículas são primeiras vizinhas uma da outra elas sentem uma atração $\epsilon$ que é a mesma para qualquer par de partículas. Essa condição modela o efeito de atração entre partículas de um gás real.

As forças de atração e repulsão num gás real não possuem alcance de mesma ordem. A repulsão é de curto alcance enquanto a atração é de longo-alcance. Embora o presente modelo trate as partículas como se o alcance de repulsão e atração fossem da mesma ordem, ainda é possível extrair propriedades físicas que tem paralelo com o gás real tais como transições de fase e formato de interfaces.

A cada ponto da rede associamos o valor $+1$ se houver uma partícula nesse ponto ou $0$ caso contrário. Representamos essa variável por $\sigma _{i}$, ou seja, no iésimo ponto da rede a variável $\sigma _{i}$ pode assumir apenas os valores $+1$ ou $0$, ou resumidamente:

${\begin{cases}0,{\text{ponto da rede vazio}}\\1,{\text{ponto da rede ocupado}}\\\end{cases}}$

A conservação do número de partículas exige que se tenha:

$\sum _{i}\sigma _{i}=\rho N$

Onde $\rho$ é a densidade de partículas e $N$ é o número total de partículas, sendo, portanto, $\rho N$ o número de pontos ocupados da rede.

O hamiltoniano do sistema é modelado a partir da condição 2 exposta acima em que é especificado que um par de primeiros vizinhos na rede contribui para a diminuição da energia do sistema por uma quantidade $\epsilon$:

$H=-\epsilon \sum _{<ij>}\sigma _{i}\sigma _{j}$

Onde $\sum _{<ij>}$ denota soma sobre todos os pares de primeiros vizinhos da rede.

Equivalência ao modelo de Ising

Para mostrar a equivalência com o modelo de Ising definimos a seguinte variável:

$s_{i}=2\sigma _{i}-1$

Essa nova variável é nada mais do que o spin no modelo de Ising para um ferromagneto assumindo os valores:

• $+1$ quando $\sigma _{i}=+1$, ou seja, posição ocupada por partícula; ou
• $-1$ quando $\sigma _{i}=0$, ou seja, posição não ocupada

Em termos da variável de spin $\sigma _{i}$ é dada por:

$\sigma _{i}={\frac {1}{2}}(s_{i}+1)$

Substituindo no Hamiltoniano tem-se:

$H=-{\frac {1}{4}}\epsilon \sum _{<ij>}(s_{i}+1)(s_{j}+1)$

$H=-{\frac {1}{4}}\epsilon \sum _{<ij>}s_{i}s_{j}-{\frac {1}{4}}\epsilon \sum _{<ij>}s_{i}-{\frac {1}{4}}\epsilon \sum _{<ij>}s_{j}-{\frac {1}{4}}\epsilon \sum _{<ij>}1$

Seja $z$ o número de coordenação da rede, ou seja, o número de primeiros vizinhos ($z=4$ para rede quadrada e $z=6$ para rede cúbica simples). Para uma dada rede existem <math\frac{1}{2}z N</math> possíveis pares distintos de primeiros vizinhos

Pode-se simplificar esssa expressão com base nas seguintes observações:

• Os somatórios em $s_{i}$ e $s_{j}$ são idênticos exceto pelo índice.
• A soma $\sum _{<ij>}s_{i}$sobre pares de vizinhos é equivalente a somar $z$ vezes sobre o número de pontos da rede:

$\sum _{<ij>}s_{i}=2z\sum _{i}s_{i}$

• $\sum _{i}s_{i}$ pode ser escrito em termos das constantes $\rho$ e $N$ assim como ocorre com $\sigma _{i}$

$\sum _{i}\sigma _{i}=\rho N$

${\frac {1}{2}}\sum _{i}(s_{i}+1)=\rho N$

${\frac {1}{2}}\sum _{i}s_{i}=\rho N-{\frac {1}{2}}\sum _{i}1$

${\frac {1}{2}}\sum _{i}s_{i}=\rho N-{\frac {1}{2}}N$

$\sum _{i}s_{i}=N(2\rho -1)$

Dessa forma o Hamiltoniano se reduz a:

$H=-{\frac {1}{4}}\epsilon \sum _{<ij>}s_{i}s_{j}-{\frac {1}{2}}z\epsilon \sum _{i}s_{i}-{\frac {1}{4}}\epsilon (2zN)$

$H=-{\frac {1}{4}}\epsilon \sum _{<ij>}s_{i}s_{j}-{\frac {1}{2}}z\epsilon (N(2\rho -1))-{\frac {1}{4}}\epsilon (2zN)$

$H=-{\frac {1}{4}}\epsilon \sum _{<ij>}s_{i}s_{j}-z\epsilon N\rho +{\frac {1}{2}}\epsilon zN-{\frac {1}{2}}\epsilon zN$

$H=-{\frac {1}{4}}\epsilon \sum _{<ij>}s_{i}s_{j}-z\epsilon N\rho$

Seja J = ${\frac {1}{4}}\epsilon$ e observando que $-z\epsilon N\rho$ é uma constante pois todos seus termos são constantes, chegamos na equivalência com o Hamiltoniano do modelo de Ising na ausência de campo magnético:

$H=-J\sum _{<ij>}s_{i}s_{j}+ctc=H_{B=0}^{Ising}+ctc$

O valor esperado $<Q>$ de qualquer quantidade física $Q$ não é alterado pela adição de uma constante ao hamiltoniano:

$<Q>=\sum _{\mu }Q_{\mu }p_{\mu }={\frac {\sum _{\mu }{Q_{\mu }}e^{-\beta (E_{\mu }+ctc)}}{\sum _{\mu }e^{-\beta (E_{\mu }+ctc)}}}={\frac {e^{-\beta ctc}\sum _{\mu }{Q_{\mu }}e^{-\beta E_{\mu }}}{e^{-\beta ctc}\sum _{\mu }e^{-\beta E_{\mu }}}}={\frac {1}{Z}}\sum _{\mu }{Q_{\mu }}e^{-\beta E_{\mu }}=<Q>$

Conservação do parâmetro de ordem

A magnetização do sistema é nada mais do que a soma de spins que já calculamos acima:

$M=\sum _{i}s_{i}=N(2\rho -1)$

No entanto, $\rho$ e $N$ devem permanecer constantes durante toda a simução, isso implica que a magnetização também é sempre constante, ou seja, a magnetização é o parâmetro de ordem conservado nesse sistema fato que dá nome ao método.

É vantajoso tratar o modelo de gás de rede sob a perspectiva de um modelo de Ising pois todo o arcabouço de técnicas amplamente conhecidas e extensivamente estudadas para o modelo de Ising podem ser aplicadas.

Apesar das similaridades, o gás de rede, como definido, possui muito menos estados válidos (${\frac {1}{2}}zN$) pois não é permitido alterar a magnetização do sistema enquanto no modelo de Ising qualquer spin individual pode ser invertido sem restrições pois a magnetização não precisa se manter constante ($2^{N}$ estados possíveis).

Transição de fase

Aproveitando a equivalência estabelecida entre gás de rede e o modelo de Ising sabe-se que o sistema possui uma transição de fase que ocorre a uma temperatura crítica $T_{c}$. Rearranjando a densidade de pontos (equivalente agora à densidade de spins up) tem-se:

$2\rho -1={\frac {M}{N}}$

$\rho ={\frac {1}{2}}(1+m)$

No modelo de Ising sabe-se também que abaixo da temperatura crítica existem dois valores de equilíbrio para a magnetização que são $+|m|$ e $-|m|$, portanto, para favorecer a coexistência de fases tem-se que:

${\frac {1}{2}}(1-|m|)=\rho _{-}\leq \rho \leq \rho _{+}={\frac {1}{2}}(1+|m|)$

Para valores de $\rho$ fora do intervalo $\rho _{-}\leq \rho \leq \rho _{+}$ ainda é possível que uma região do sistema favoreça uma das duas densidades preferenciais. Suponha que se tenha $\rho <\rho _{-}$. Nesse caso o sistema possui menos partículas do que precisa pra atingir o a densidade $\rho _{+}$. Ainda que localmente seja possível o sistema atingir a densidade $\rho _{+}$ isso leva a uma falta ainda maior de partículas em outras regiões do sistema sendo, portanto, energeticamente custoso. A opção energeticamente mais favorável adotada pelo sistema é distribuir as poucas partículas homegeneamente pela rede. Esse comportamento é observado na simulação.

Dessa forma, no caso de $J>0$ o sistema possui duas fases:

• Uma em que $\rho \in [\rho _{-},\rho _{+}]$ se dividindo em dois domínios cada qual favorecendo uma das duas densidades $\pm \rho$
• E outra em que $\rho \not \in [\rho _{-},\rho _{+}]$ tendo densidade homogênea

Com $\rho$ sujeito ao intervalo ${\frac {1}{2}}(1-|m|)\leq \rho \leq {\frac {1}{2}}(1+|m|)$ conclui-se que $\rho$ pode assumir um intervalo menor de valores a medida que $|m|$ diminui. A magnetização diminui sob o aumento da temperatura. Acima da temperatura crítica a $|m|=0$ e portanto o intervalo ${\frac {1}{2}}(1-|m|)\leq \rho \leq {\frac {1}{2}}(1+|m|)$ reduz-se a zero evidenciando que não existe mais um valor de $\rho$ que evite a homogeinização da rede.

A discussão acima pode ser apresentada resumidamente pelo diagrama de fases:

Diagrama de fases do modelo CPO. Fase homogênea para temperaturas além da temperatura crítica e fase coexistente abaixo com densidades preferenciais $\rho _{+}$ e $\rho _{-}$

Esse comportamento é observado quando se diminui a temperatura de vapor d'agua que passa a formar gotas líquidas que coexistem com o vapor para um intervalo de temperaturas. A fase condensada do gás de rede, no entanto, é mais adequadamente interpretada como um sólido devido a posição fixa das partículas (análogas a moléculas ou átomos) na rede, dessa forma, falamos de interface vapor/sólido ao invés de vapor/líquido.

Teoria

Sistemas físicos em equilíbrio com muitos graus de liberdade e no limite termodinâmico comportam-se de tal forma que ao flutuarem de um estado $\mu$ para um estado $\nu$ tem-se que $\nu$ difere pouco de $\mu$. Outra maneira de dizer isso é que as flutuações dessa tipo de sistema físico são muito pequenas em relação ao número de configurações possíveis e que portanto o sistema passa a maior parte do tempo alternando entre um pequeno conjunto de configurações. A consequência disso é que pode-se escolher uma estratégia de visitar com maior probabilidade apenas a fração de estados do sistema, as quais mais contribuem para atingir o equilíbrio ao invés de se visitar todos os estados indistintamente. No modelo de ferromagneto, por exemplo, com uma rede $10\times 10\times 10$, há $2^{1000}\simeq 10^{300}$ configurações possíveis sendo que mesmo com um supercomputador seria impraticável realizar essa simulação. O método de Monte Carlo consiste em visitar eficientemente uma pequena fração desses estados e atingir rapidamente o equilíbrio em poucos passos e o peso que define como visitar o estado seguinte é dado pela distribuição de Boltzmann $e^{\beta (E_{\nu }-E_{\mu })}$ onde fica claro que quanto mais diferente $\nu$ for de $\mu$ menor a change de fazer a transição $\mu \to \mu$

Dessa forma impõe-se que no equilíbrio o sistema obedeça a distribuição de Boltzmann, portanto a condição de balanço detalhado dá liberdade na escolha de $P(\mu \to \nu )$ e $P(\nu \to \mu )$ desde que seja satisfeita:

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

Uma possível escolha para $A(\mu \to \nu )$ seria:

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

Tomando $g(\mu )=g(\nu )$ e como $A_{0}$ é cancelada na razão entre probabilidades de aceitação temos a liberdade na sua escolha desde que mantenha a probabilidade menor ou igual a um. No modelo de ferromagneto, por exemplo, a maior diferença de energia que se pode obter entre estados é $\Delta E=E_{\nu }-E_{\mu }=\pm 2zJ$ o que significa que o maior valor de $e^{-{\frac {1}{2}}\beta (E_{\nu }-E_{\mu })}$ é justamente $e^{\beta zJ}$. Assim, para garantir que a probabilidade seja menor ou igual a 1 deve-se escolher $A_{0}\leq e^{-\beta zJ}$

Para que o algoritmo seja eficiente deseja-se que a probabilidade de aceitação seja a maior possível, pois do contrário estaríamos utilizando tempo computacional apenas para rejeitar trocas de estado. Portanto queremos que $A_{0}$ assuma o maior valor possível $A_{0}=e^{-\beta zJ}$, maximizando $A(\mu \to \nu )$:

$A(\mu \to \nu )=e^{-{\frac {1}{2}}\beta (E_{\nu }-E_{\mu }+2\beta zJ)}$

Devido a condição de balanço detalhado, essa escolha implica:

$A(\nu \to \mu )=e^{-\beta zJ}$

Metropolis percebeu que desde que a condição de balanço detalhado seja satisfeita tem-se liberdade na escolha das probabilidades de aceitação. Então ele decidiu atribuir o maior valor possível para a probabilidade de aceitação que tem o maior valor entre as duas, no caso $A(\nu \to \mu )$, ou seja:

$A(\nu \to \mu )=1$

O que implica:

$A(\mu \to \nu )=e^{-\beta (E_{\nu }-E_{\mu })}$

Dessa forma a transição de estados sempre ocorre se $E_{\nu }<E_{\mu }$ ou seja $\Delta E<0$ mas pode ou não ocorrer caso seja $\Delta E>0$ com uma probabilidade dada por $e^{-\beta (E_{\nu }-E_{\mu })}$. Em suma:

$A(\mu \to \nu )={\begin{cases}e^{-\beta (E_{\nu }-E_{\mu })}\ {\mbox{ se }}\ E_{\nu }-E_{\mu }>0\\1\ {\mbox{ caso contrario}}\\\end{cases}}$

Implementação

Para obedecer a condição de conservação da magnetização não é permitido alterar um spin individualmente (ou um número ímpar de spins). Uma maneira de tratar a dinâmica desse sistema foi proposta por Kawasaki e consiste em simplesmente alternar o estado de spin de um par de partículas que tenham estados de spin oposto, ou seja:

${\begin{cases}\uparrow \uparrow ou\downarrow \downarrow \quad \to \quad {\text{nada a fazer}}\\\uparrow \downarrow \quad \to \quad \downarrow \uparrow \\\downarrow \uparrow \quad \to \quad \uparrow \downarrow \\\end{cases}}$

É evidente que nesse caso a mudança conserva a magnetização, pois a troca de spins resulta em variação de magnetização nula.

Cada ponto da rede possui $z$ vizinhos e portanto a cada passo de iteração deve-se sortear com qual dos vizinhos será feita uma tentativa de troca de spins. Essa escolha é feita aleatoriamente (uniforme). Uma vez escolhido um vizinho deve-se decidir se a troca deve ser feita ou não. Essa decisão é tomada com base no método de Monte Carlo, em particular, com a probabilidade de aceitação de Metropolis exatamente como exposto na seção acima.

A ergodicidade é satisfeita pelo sistema pois um passo de Monte Carlo corresponde a uma troca entre vizinhos que numa rede finita pode ser efetuada a partir de outro estado qualquer em número finito de passos

Como já foi mencionado a rede possui $N$ pontos e número de coordenação $z$ o que resulta em ${\frac {1}{2}}zN$ pares de primeiros vizinhos, portanto, a probabilidade de selecionar um par qualquer é dada por:

$g(\mu \to \nu )={\frac {1}{{\frac {1}{2}}zN}}={\frac {2}{zN}}$

A probabilidade de seleção $g(\nu \to \mu )$ é a mesma fazendo com que esses termos se cortem na condição de balanço detalhado e permitindo que se aplique a escolha de Metropolis discutida acima sem alterações.

Para efetivamente tomar a decisão sobre a troca entre vizinhos onde $\Delta E$ é necessário especificar como é feito seu cálculo. $\Delta E$ é dado pela seguinte expressão:

$\Delta E=E_{\nu }-E_{\mu }=2J\left[s_{k}^{\mu }\sum _{i\not \in \{k',k\}}s_{i}^{\mu }+s_{k'}^{\mu }\sum _{j\not \in \{k,k'\}}s_{j}^{\mu }\right]$

Seja um ponto da rede $k$ e $k'$ primeiro vizinho de $k$. Deseja-se calcular a energia de interação entre esse par. A expressão acima apenas diz que deve-somar os produtos do spin de $k$, $(s_{k})$, com seus primeiros vizinhos $s_{i}$ excluindo-se da soma tanto $k$ quanto $k'$. Faz-se o mesmo procedimento para $k'$, ou seja, soma-se os produtos do spin $k'$, $(s_{k'})$, com todos os seus primeiros vizinhos $s_{j}$ exceto ele mesmo e $k$. A soma dessas duas quantidades multiplicadas por $2J$ é igual a diferença de energia entre a configuração $\mu$ e a $\nu$

A simulação é extremamente simplificada pelo fato de o estado final não figurar na expressão para o cálculo da diferença de energia da transição, pois pode-se verificar de antemão se a transição deve ser feita sem que seja necessário efetivamente colocar o sistema no estado final $\nu$.

Simulação

Foram simulados três sistemas diferentes os quais são discutidos a seguir. O objetivo das simulações em determinar a tendência do formato de cada interface após vários passos de monte carlo sem verificar rigorosamente se o equilíbrio foi atingido. Uma análise das condições de equilíbrio é discutida na última seção. Todas as simulações demandam a geração de muitos números aleatórios por isso optou-se por usar Mersenne Twister - um algoritmo veloz para geração de números aleatórios.

Interface linear

Esse sistema consiste de uma rede quadrada de aresta $l$. A rede tem inicialmente sua metade inferior completamente populada por partículas (spins up) e o restante da rede possui vacâncias (spins down).

Cada ponto da rede é visitado sequencialmente e um dos seus 4 vizinhos é sorteado para que se faça uma tentativa de troca de posição entre o par. Caso a energia do sistema diminua, a troca é feita com probabilidade 1. Caso contrário a troca ocorre com probabilidade dada pelo peso de Boltzmann (algoritmo de Metropolis).

A primeira condição de contorno refere-se às paredes superior e inferior. A ultima linha de pontos da rede possui spins apontando pra baixo permanentemente assim como a primeira linha de pontos da rede tem-se spins apontando pra cima. Para que essa configuração seja mantida ao longo da simulação, evita-se qualquer tentativa de troca entre pares envolvendo esses pontos. Essa condição de contorno garante que a interface se mantenha fixa, do contrário, ela poderia transitar pela rede. Nas paredes laterais aplica-se condição de contorno periódica onde por exemplo um ponto da rede na coluna $l-1$ pode trocar de lugar com um primeiro vizinho da coluna $0$ e vice-versa.

Ao iterar pela rede é comum deparar-se com pares de vizinhos que possuem mesma orientação de spin e portanto são ignorados imediatamente pois em nada contribuem para a dinâmica do sistema.

Interface linear entre sólido e vapor. Cada frame corresponde a 10 passos de Monte Carlo de um total de 500 passos. Primeira simulação com alta temperatura $T>T_{C}$. Segunda simulação com temperatura intermediária $T<T_{C}$. Terceira simulação com baixa temperatura $T\ll T_{C}$

Como observado, a alta temperatura recai-se no regime homogêneo em que a alternativa mais adequada para que o sistema minize sua energia é distribuir os spins up (partículas) uniformemente pela rede. A uma temperatura abaixo da crítica percebe-se a formação de uma interface persistente entre a região inferior com densidade preferencial $\rho _{+}$ e a superior com densidade preferencial $\rho _{-}$ em coexistência. Ao diminuir ainda mais a temperatura o efeito fica mais acentuado, ou seja, a interface é menos deformada em relação ao caso anterior.

É possível escrever um algoritmo mais eficiente. Ao invés de percorrer toda a rede e testar ponto a ponto pela possibilidade de uma transição, armazena-se em memória duas listas, uma de spins up e outra de spins down e sorteia-se dois elementos de cada lista para que se faça a transição. Dessa forma, não se perde tempo com rejeições. No entanto, para sistemas pequenos com alguns milhares de passos de Monte Carlo e relativamente poucos pontos como o sistema estudado não foi necessário esse ganho de performance.

Interface circular

Nesse sistema exclui-se a condição de contorno que fixa os spins das paredes superior e inferior e atribui-se a elas as mesmas condições das paredes laterais, ou seja, condições de contorno periódicas onde, por exemplo, uma partícula na linha $l-1$ pode trocar de lugar com sua primeira vizinha da linha $0$ Como condição inicial posiciona-se as partículas num formato quadrado no centro da rede e observa-se como o formato da interface varia ao longo da simulação. A densidade de partículas é menor que no exemplo anterior da interface linear.

Interface circular entre sólido e vapor. Cada frame corresponde a 10 passos de Monte Carlo de um total de 500 passos. Primeira simulação com alta temperatura $T>T_{C}$. Segunda simução com temperatura intermediária $T<T_{C}$. Terceira simulação com baixa temperatura $T\ll T_{C}$

A interface é energeticamente custosa pois para cada par de spins antialinhados o sistema aumenta de energia por um fator 2J e como na rede quadrada um spin da interface possui 3 vizinhos antialinhados, sistema aumenta de energia por um fator 6J. Portanto fisicamente espera-se que o sistema evolua de tal forma a minimizar a extensão da sua interface, minimizando sua tensão superficial. No caso simulado espera-se que um domínio circular seja gerado pois o círculo é a forma geométrica que possui menor perímetro. No entanto, como a simulação demonstra, mesmo pra baixas temperaturas a forma nunca é perfeitamente circular e isso se deve ao tamanho finito da rede o faz com que seu formato (da rede) influencie o formato da interface.

A animação abaixo ilustra o mesmo processo mas com menos frames por segundo permitindo acompanhar detalhes da dinâmica do sistema.

Animação com 100 passos de Monte Carlo. Primeira simulação com alta temperatura $T>T_{C}$. Segunda simução com temperatura intermediária $T<T_{C}$. Terceira simulação com baixa temperatura $T\ll T_{C}$

Interface esférica

A simulação da interface esférica é uma extensão direita da simulação da interface circular apenas adicionando mais uma dimensão. Observa-se os mesmos efeitos de redução de tensão superficial pela deformação do cubo em uma região aproximadamente esférica quando a temperatura é menor que a temperatura crítica. Acima da a temperatura crítica a densidade fica homogênea como esperado.

Cada ponto da rede agora possui 6 vizinhos ao invés de 4 e isso faz com que as haja valores adicionais de diferenças de energia entre estados (na rede quadrada era possível obter $\Delta E\in \{0,\pm 4,\pm 8,\pm 12\}$ e agora na rede cúbica $\Delta E\in \{0,\pm 4,\pm 8,\pm 12,\pm 16,\pm 20\}$

Visualização da vizinhança de um ponto de uma rede cúbica. Download for free at http://cnx.org/contents/85abf193-2bd2-4908-8563-90b8a7ac8df6@9.311

Interface esférica entre sólido e vapor. Cada frame corresponde a 10 passos de Monte Carlo de um total de 500 passos. Primeira simulação com alta temperatura $T>T_{C}$. Segunda simução com temperatura intermediária $T<T_{C}$. Terceira simulação com baixa temperatura $T\ll T_{C}$

A mesma simulação com menos partículas, vista mais distante e com uma pequena diferença na quantidade de passos.

Interface esférica entre sólido e vapor. Cada frame corresponde a 5 passos de Monte Carlo de um total de 250 passos. Primeira simulação com alta temperatura $T>T_{C}$. Segunda simução com temperatura intermediária $T<T_{C}$. Terceira simulação com baixa temperatura $T\ll T_{C}$

Introduzindo interações entre segundos vizinhos é possível reproduzir formatos de cristais cúbicos como por exemplo o cristal de face centrada ou de corpo centrado.

Equilíbrio

Para medir qualquer grandeza de um sistema simulado pelo método de Carlo é necessária que a medida seja feita sob regime de equilíbrio. Torna-se então importante saber quando o sistema atinge o equilíbrio. No caso de um ferromagneto no modelo de Ising pode-se monitorar a magnetização ou calor específico até que a grandeza atinja um valor estacionário. No caso do gás de rede podemos monitorar o formato da interface, ou mais simplesmente, a densidade média de spins up (partículas) em cada coluna da rede quadrada. Acompanhando a mudança percentual nessa densidade ao passar de um passo de Monte Carlo para o passo seguinte é possível ter uma idéia de quanto passos são necessários para atingir o equilíbrio.

Densidade percentual média em função de número de passos de Monte Carlo (até 20000 passos)

Olhando mais de perto o início da curva percebe-se que em torno de 1500 passos o equilíbrio ja foi seguramente atingido.

Densidade percentual média em função de número de passos de Monte Carlo (até 2000 passos)

As simulações exibidas acima estão, portanto, na região fora do equilíbrio mas como observado acima o objetivo das simulações era determinar a tendência do formato das interfaces e não o seu formato no equilíbrio.

Códigos

Conservação de parâmetro de ordem - Interface linear

Conservação de parâmetro de ordem - Interface circular

Conservação de parâmetro de ordem - Interface esférica

Conservação de parâmetro de ordem - Equilíbrio