Grupo - Lennard Jones: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
Sem resumo de edição
 
(177 revisões intermediárias por 2 usuários não estão sendo mostradas)
Linha 1: Linha 1:
O potencial devido a interação entre duas partículas separadas por uma distância <math>r'</math> pode ser modelado pelo potencial de Lennar Jones:
Alunos: Ben-Hur Cardoso e Vinícius Garcia


<math>U(r') = 4 \epsilon \left [ \left ( \frac{\sigma}{r'} \right )^{12} + \left ( \frac{\sigma}{r'} \right )^{6} \right ].</math>
 
O potencial devido a interação entre duas partículas separadas por uma distância <math>r'</math> pode ser modelado pelo potencial de Lennard-Jones <ref>Lennard-Jones, J. E. (1924), "On the Determination of Molecular Fields", Proc. R. Soc. Lond. A, 106 (738): 463–477</ref>:
 
<math>U'(r') = 4 \epsilon \left [ \left ( \frac{\sigma}{r'} \right )^{12} - \left ( \frac{\sigma}{r'} \right )^{6} \right ]</math>


Posto em unidades reduzidas (<math>r \equiv r'/ \sigma</math> e <math>U \equiv U' / \epsilon</math>), o potencial reduz-se a:   
Posto em unidades reduzidas (<math>r \equiv r'/ \sigma</math> e <math>U \equiv U' / \epsilon</math>), o potencial reduz-se a:   


<math>U(r) = 4  \left [ \left ( \frac{1}{r} \right )^{12} + \left ( \frac{1}{r} \right )^{6} \right ].</math>
<math>U(r) = 4  \left [ \left ( \frac{1}{r} \right )^{12} - \left ( \frac{1}{r} \right )^{6} \right ]</math>
 
[[Arquivo:LennardJones.png|500px|center]]


No caso tem-se, por conveniência, o seguintes sistema de unidades básicas:
As quantidades físicas de interesse (como pressão, tempo e temperatura) possuem ordens de grandeza muito diferentes entre si quando se usa unidades do Sistema Internacional, o que faz com que fique difícil de identificar erros. Por causa disso, é computacionalmente viável trabalhar com o seguintes sistema de unidades básicas, que mantém os valores em baixas ordens de grandeza:


{| class="wikitable"
{| class="wikitable"
|-
|-
! '''Grandeza'''
! '''Grandeza'''
! '''Comprimero'''
! '''Comprimento'''
! '''Tempo'''
! '''Tempo'''
! '''Massa'''
! '''Massa'''
Linha 30: Linha 35:
|}
|}
onde <math>m_p</math> é a massa da partícula e <math>k_B</math> é a constante de Boltzmann.
onde <math>m_p</math> é a massa da partícula e <math>k_B</math> é a constante de Boltzmann.
.
 
Note que as novas unidades são adimensionais. Assim, para comparar o modelo com um sistema físico real, deve-se converter para as unidades desejadas. É interessante notar que a simulação feita utilizando as unidades reduzidas permite visualizar que sistemas distintos serão equivalentes, já que há infinitas combinações de densidade, temperatura, <math>\epsilon </math> e <math>\sigma </math> com mesmas unidades reduzidas.
 
== Método Monte Carlo ==
== Método Monte Carlo ==


Dado uma amostra com <math>N</math> partículas, a abordagem introduzida por Metropolis segue o seguinte esquema:
Denomina-se método de Monte Carlo métodos estatísticos que se baseiam em amostragem aleatória massiva para cálculo numérico.
 
=== Amostragem simples ===
 
O cálculo numérico de uma integral utilizando Monte Carlo parte da ideia de que uma integral pode ser reescrita como:
 
<math> F = \int_a^b{f(x) dx} = (b - a)\langle f(x) \rangle </math>
 
Dessa forma, utiliza-se amostragem aleatória massiva para estimar <math> \langle f(x) \rangle </math>, que é a média da função no intervalo de interesse. O estimador para a média real é a média amostral, que por ser uma medida estatística terá sempre um erro inerente associado.
 
=== Amostragem por importância ===
 
Um problema da amostragem simples é que ela utiliza uma distribuição uniforme, que pode, para uma função que tenha regiões onde seu valor é próximo de zero por exemplo, custar a estimar corretamente o valor médio da função. A fim de contornar isso, podemos utilizar uma distribuição <math> w(x) </math> tal que a razão <math> \frac{f(x)}{w(x)} </math> seja o mais constante possível. Reescrevendo a integral:
 
<math> F = \int_a^b{\frac{f(x)}{w(x)} w(x) dx} = (b - a)\left \langle \frac{f(x)}{w(x)} \right \rangle </math>
 
 
Basicamente, é a média dos valores da função razão <math> \frac{f(x)}{w(x)} </math> com distribuição <math> w(x) </math> para sorteio dos pontos. Esse método diminui a variância amostral e faz com que seja muito mais eficiente a convergência do estimador média amostral para a média real da função.
 
=== Algoritmo de Metropolis–Hastings ===
 
Dado uma amostra com <math>N</math> partículas, a abordagem introduzida por Metropolis (em 1953, algoritmo para distribuições simétricas) e Hastings (em 1970, algoritmo generalizado) segue o seguinte esquema:
 
(1) Selecionar uma partícula aleatóriamente, e calcular sua energia <math>U(r)</math>;
(2) Dado o deslocamento <math>\mathbf{r_n} = \mathbf{r} + \mathbf{\Delta}</math>, calcular <math>U(r_n)</math>;
(3) Aceitar o movimento <math>\mathbf{r} \rightarrow \mathbf{r_n}</math> com probabilidade <math>p = min\{1; \exp[-\beta (U(r_n) - U(r))]\}</math>
 
Devido ao fato das componentes de <math>\mathbf{\Delta}</math> poderem assumir valores postivos e negativos (com a nova posição limitada por alguma condição de contorno), após ser aceito a mudança de posição de uma dada partícula, elá poderá ser selecionada (no próximo passo) com a mesma probabilidade e voltar a sua posição original. Isso garante a condição de balanço detalhado. <ref> Daan Frenkel, Berend Smit, "Understanding molecular simulation: from algorithms to applications"</ref>
 
=== Estimadores no Equilíbrio ===
Em todos os exemplos tratados aqui, será usado o ensemble NVT (com o número de partículas, volume e temperatura constantes) com <math>d</math> dimensões. Dado isso, os sistemas são caracterizados com um densidade e uma temperatura. Com tais sistemas no equilíbrio, são estimadas (média de sucessivas medidas) a energia total e a pressão, dadas respectivamente por:
<math>U_T = \sum_i \sum_{j>i}U(r_{ij}) </math>
 
<math>P = \frac{\rho}{\beta} + \frac{2\rho}{Nd}\sum_i \sum_{j>i}\mathbf{f(\mathbf{r_{ij}})} \cdot \mathbf{r_{ij}} </math>
 
onde <math>\mathbf{r_{ij}} = \mathbf{r_{j}} - \mathbf{r_{i}}</math>. Além disso, é interessante a análise da capacidade térmica
 
<math>C_v = \frac{\langle U_T^2 \rangle - \langle U_T \rangle ^2}{k_BT^2}</math>
 
A função radial de distribuição de pares, é uma função que estima o quão provável é encontrar duas partículas a uma distância <math>r</math> dentro de um sistema de várias partículas. Em um sistema de <math>N</math> partículas, o <math>g(r)</math> é definido como a média do número de partículas a uma distância <math>r</math>:
 
<math> g(r)=\frac{1}{N\rho}\langle\sum_{i=1}^N\sum_{j\neq i}^N\delta(r-|r_i-r_j|)  \rangle</math>
 
A distribuição de densidades é calculada dividindo espacialmente o sistema, no nosso caso de lado <math> L </math>, em <math> M </math> células de lado <math> \frac{L}{M} </math> e acumular em um histograma as densidades das células. Essa medida é relevante para identificação de fases no sistema.
 
== Detalhes Técnicos==
 
=== Condições de Contorno ===
 
Qualquer sistema possível de ser feito hoje com método Monte Carlo, apesar do grande poder computacional disponível, fica distante do limite termodinâmico. As condições de contorno podem ser estabelecidas de forma a tentar contornar isso. As condições utilizadas neste trabalho foram condições de contorno periódicas, que possibilitam que o sistema se comporte como se fosse muito maior do que é, desde que ele seja isotrópico.
 
=== Convenção da imagem mínima ===
 
Construindo um sistema cúbico de lado L, com um número N de partículas fixo, as condições de contorno periódicas fazem com a simulação não se limite ao cubo. Ao calcular a interação entre duas partículas, portanto, a distância entre elas não é unívoca, já que há incontáveis cópias de cada uma, repetindo-se periodicamente. Resolve-se esse problema utilizando-se da convenção da imagem mínima, em que é calculada a interação com a imagem mais próxima das outras.
 
=== Truncamento nas interações ===
 
Um problema da condição de contorno periódica é que, a princípio, cada partícula interagiria com todas as outras do sistema, que devido ao fato de ser periodicamente repetido, seriam infinitas. Como o potencial utilizado é de curto alcance, deve ser possível, de alguma forma, limitar as interações entre as partículas sem perda dos significados numéricos da simulação. Pode-se, então, truncar as interações de uma partícula <math> i </math> a partir de uma distância de corte <math> r_c </math> de forma que há uma descontinuidade no potencial dessa partícula na esfera de raio <math> r_c </math>.
 
 
[[Arquivo: mic.png|300px|center]]
 
 
Assim o potencial simulado é:
 
<math> U_{trunc}(r) = \begin{cases}U_{LJ}(r) & , r \leq r_c\\0 & , r > r_c \end{cases}  </math>
 
Esse é o chamado truncamento simples. O subíndice LJ se refere ao potencial de Lennard-Jones. O problema é que isso cria uma contribuição indesejada à pressão, já que há uma força de impulso por conta dessa descontinuidade do potencial.
 
Pode-se, então, usar o truncamento com deslocamento, que evita essa descontinuidade fazendo uma subtração em todo ponto do módulo do potencial de Lennard-Jones à distância <math> r_c </math>. Assim, o potencial simulado é:
 
<math> U_{deslocado}(r) = \begin{cases}U_{LJ}(r) - U_{LJ}(r_c) & , r \leq r_c\\0 & , r > r_c \end{cases}  </math>
 
Com esse potencial, que é contínuo, as forças serão sempre finitas, o que retira a força de impulso que alteraria medidas de pressão.
 
=== Translação ===
 
A possível nova posição de uma partícula será <math>(x_n, y_n, z_n)</math> tal que
 
  <math> x_n = x + \Delta (\alpha - 0.5)</math>
  <math> y_n = y + \Delta (\alpha - 0.5)</math>
  <math> z_n = z + \Delta (\alpha - 0.5)</math>
 
sendo <math>\alpha</math> números aleatórios uniformemente distribuídos e <math> \Delta </math> o deslocamento máximo permitido. Se a nova posição for aceita seguindo o algorítmo de Metropolis, a partícula assumirá essa posição. Caso contrário, a antiga configuração contribuirá novamente na média.  <ref> David P. Landau, Kurt Binder, "A Guide to Monte Carlo Simulations in Statistical Physics"</ref>
 
O tamanho de <math> \Delta </math> não pode ser grande ao ponto de nunca ser aceito nem pequeno ao ponto de sempre ser aceito. Dado isso, buscando uma taxa de aceitação <math> \eta </math> entre <math> 30 \%</math> e <math> 50 \% </math>, define-se um <math> \Delta_0 </math> inicial que vai se ajustando no decorrer da simulação, executando - a cada passo Monte Carlo (<math>N</math> passos de simulação) - o seguinte algorítimo:
 
  Se (<math> \eta < 0.3 </math>): <math> \Delta \leftarrow \Delta \times (\eta / 0.3)</math> 
  Se (<math> \eta > 0.5 </math>): <math> \Delta \leftarrow \Delta \times (\eta / 0.5)</math>
 
== Sistema bidimensional ==
 
Dado um sistema com densidade <math> \rho = N / L^2 </math> e temperatura <math>T</math>, os sistemas foram feitos com:
 
(A) <math> N = 256 </math> partículas
(B) Quadrado de lado <math> L =  (N / \rho)^{1/2} </math>, utilizando a convenção da imagem mínima
(C) Incialização aleatória
(D) Distância de corte <math>r_c = 2.5 \sigma</math>
(E) Deslocamento máximo inicial <math>\Delta_0 = L/2 </math>
=== Transição de fase líquido-gás ===
A partir do artigo de M Rovere et al <ref> M Rovere et al (1990) "The gas-liquid transition of the two-dimensional Lennard-Jones fluid".  J. Phys.: Condens. Matter 2 7009 </ref>, que modelou o sistema apresentado, obtemos a informação da temperatura e da densidade críticas. Os valores encontrados pelos autores do artigo são <math> T = 0.5</math> e <math> \rho = 0.3 </math>, que determinam qual intervalo de valores de temperatura e densidade apresentam coexistência de fases distintas.
 
=== Transição de fase com densidade fixa ===
 
Fixada a densidade <math> \rho = 0.3 </math>, foram realizadas simulações para as diferentes temperaturas <math> T \in \{0.45, 0.50, 0.70\}</math>.
 
Evolução temporal da Energia Total:


<math>\frac{\Psi^{n}_{j-1} - 2\Psi^{n}_{j} + \Psi^{n}_{j+1}}{\left(\Delta x \right)^2}</math>
[[Arquivo:time_evolution_U_d_03.png|500px|center]]


e as discretizações de <math>\frac{\partial \Psi}{\partial t}</math> (explícita e implícita, respectivamente):
Observa-se que o sistema se encontra em equilíbrio a <math>10^5</math> passos  Monte Carlo. Isso é corroborado pelos resultados do artigo de M Rovere et al. A partir do sistema equilibrado, são analisadas a distribuição de densidade e a função distribuição de pares para verificar a presença de uma (líquida) ou duas fases (líquida, gasosa), usando uma média de <math>10^5</math> medidas sucessivas.


<math>\frac{\partial \Psi}{\partial t} = \frac{\Psi^{n+1}_{j} - \Psi^{n}_{j}}{\Delta t}</math>  (explícita)
Para <math> T = 0.70</math> temos que:


<math>\frac{\partial \Psi}{\partial t} = \frac{\Psi^{n}_{j} - \Psi^{n-1}_{j}}{\Delta t}</math>  (implícita)
[[Arquivo: snapshot_t70_d_03.png|300px]] [[Arquivo: P_t70_d_03.png|300px]] [[Arquivo: g_t70_d_03.png|300px]]


Tanto no método explícito quanto no método implícito não é conservada a norma do estado (o que é estritamente necessário, que o estado pode ser interpretado como uma onda de probabilidade). Por esse motivo, utiliza-se o método de Crank-Nicolson, o qual tem essa propriedade.
Nesse caso observa-se apenas uma fase, com a densidade média sendo a densidade global do sistema, que é diferente do centro do pico (<math> \rho \approx 1.0 </math>). Isso caracteriza a fase gasosa, que ocupa todo o espaço disponível. A função distribuição de pares rapidamente converge para 1 com pequenas oscilações, caracterizando um gás de alta densidade.


O método de Crank-Nicolson consiste em uma média aritmética dos métodos explícito e implícito. Excetuando manipulações algébricas triviais, verifica-se que a relação de recorrência do método é dada por:
Para <math> T = 0.50</math> temos  que:


<math>a\Psi^{n+1}_{j-1} + b_{j}\Psi^{n+1}_{j} + a\Psi^{n+1}_{j+1} = a^* \Psi^{n}_{j-1} + b_{j}^{*} \Psi^{n}_{j} + a^*\Psi^{n}_{j+1},</math>
[[Arquivo: snapshot_t50_d_03.png|300px]] [[Arquivo: P_t50_d_03.png|300px]] [[Arquivo: g_t50_d_03.png|300px]]


onde
Nesse caso observa-se a coexistência das fases gasosa e líquida. Além disso, pode-se observar a densidade característica aproximada de cada fase a essa temperatura centrada nos respectivos picos (<math> \rho_g \approx 0.2 </math> e <math> \rho_l \approx 2.0 </math>). A função distribuição de pares representa a de uma fase líquida.


<math>a \equiv -\frac{i \Delta t}{4(\Delta x)^2}</math> e <math>b_{j} \equiv 1+\frac{ i\Delta t}{2} \left[\frac{1}{(\Delta x)^2} + V(j \Delta x) \right]</math>.
Para <math> T = 0.45</math> temos  que:


A integração numérica depende, portanto, do potencial em que o elétron está sujeito, bem como da sua condição inicial e suas das condições de contorno.
[[Arquivo: snapshot_t45_d_03.png|300px]] [[Arquivo: P_t45_d_03.png|300px]] [[Arquivo: g_t45_d_03.png|300px]]


Que condições podemos impor para a fronteira? Quando se trata do problema analiticamente, costuma-se considerar que a função de onda tende a zero no infinito. Numericamente, pode-se fazer uma transposição disso, criando uma condição para bordas em pontos suficientemente distantes do centro da distribuição da função de onda, igualando-as a zero. Outra forma de tratar o problema numericamente é criando condições de contorno periódicas, em que para as bordas vale <math>\Psi^{n}_{0} = \Psi^{n}_{j_{max}}</math> para todo <math>n</math> (ou, para as bordas <math>a</math> e <math>b</math> há a relação <math>\Psi (a, t) = \Psi (b, t)</math> para todo <math>t</math>).
Nesse gráfico a primeira impressão é de coexistência de fases. Apesar disso, densidade zero não corresponde a uma fase real, e sim aos espaços vazios. Dessa forma, a fase aqui é a líquida. O ''snapshot'' permite visualizar isso, que há aglomerados de densidade aparentemente constante separados por espaços praticamente sem partículas. A densidade da fase líquida a essa temperatura é aproximadamente o centro do pico da direita (<math> \rho \approx 2.3 </math>).


===Condições de contorno iguais a zero===
Para <math> T = 0.30</math> temos  que:
Para as condições de contorno <math>\Psi_0^n = \Psi_L^n = 0</math>, a iteração reduz-se à equação matricial:


:<math>
[[Arquivo: snapshot_t30_d_03.png|300px]] [[Arquivo: P_t30_d_03.png|300px]] [[Arquivo: g_t30_d_03.png|300px]]
\begin{pmatrix}
b_1 & a & 0 & \cdots & 0 & 0 & 0\\ a & b_2 & a & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & a & b_{L-2} & a \\ 0 & 0 & 0 & \cdots & 0 & a & b_{L-1}
\end{pmatrix}
\begin{pmatrix} \Psi_1^{n+1} \\ \Psi_2^{n+1} \\ \vdots \\ \Psi_{L-2}^{n+1} \\ \Psi_{L-1}^{n+1} \end{pmatrix}
=
\begin{pmatrix} b_1^* & a^* & 0 & \cdots & 0 & 0 & 0\\ a^* & b_2^* & a^* & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & a^* & b_{L-2}^* & a^* \\ 0 & 0 & 0 & \cdots & 0 & a^* & b_{L-1}^*\end{pmatrix}
\begin{pmatrix} \Psi_1^{n} \\ \Psi_2^{n} \\ \vdots \\ \Psi_{L-2}^{n} \\ \Psi_{L-1}^{n} \end{pmatrix}
</math>


Assim como no gráfico anterior, há dois picos mas apenas um representa uma fase. A principal diferença é que o pico em zero é bem maior e mais agudo que no caso <math> T = 0.45</math>, o que ilustra que há mais espaço vazio. A densidade do líquido aqui é <math> \rho \approx 2.8 </math>. Nota-se que o sistema parece estar se organizando numa fase hexagonal, apesar de ainda ser amorfa, e a densidade já é bem mais alta. A distribuição de pares tem mínimos bem menores do que para os parâmetros acima tratados, caracterizando um sistema que parece estar próximo da transição para a fase hexagonal de baixa energia.


===Condições de contorno periódicas===
De maneira geral, é possível estimar o intervalo de parâmetros para o qual ocorre a transição de fase líquido-gás. É possível, também, visualizar, a partir da distribuição de densidades, a coexistência de fases, já que cada estado tem uma densidade média característica a dada temperatura.
De maneira semelhante, a iteração do caso das condições de contorno periódicas - <math>\Psi_0^n = \Psi_L^n</math> - reduz-se à equação matricial:
:<math>
\begin{pmatrix} b_0 & a & 0 & \cdots & 0 & 0 & a\\ a & b_1 & a & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & a & b_{L-1} & a \\ a & 0 & 0 & \cdots & 0 & a & b_{L}\end{pmatrix}
\begin{pmatrix} \Psi_0^{n+1} \\ \Psi_2^{n+1} \\ \vdots \\ \Psi_{L-1}^{n+1} \\ \Psi_{L}^{n+1} \end{pmatrix}
=
\begin{pmatrix} b_0^* & a^* & 0 & \cdots & 0 & 0 & a^*\\ a^* & b_1^* & a^* & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & a^* & b_{L-1}^* & a^* \\ a^* & 0 & 0 & \cdots & 0 & a^* & b_{L}^*\end{pmatrix}
\begin{pmatrix} \Psi_0^{n} \\ \Psi_2^{n} \\ \vdots \\ \Psi_{L-1}^{n} \\ \Psi_{L}^{n} \end{pmatrix}
</math>


=== Transição de fase com temperatura fixa ===


===Condição inicial===
Dado a temperatura fixa <math> T = 0.45 </math>, é observado o sistema para as diferentes densidades <math> \rho \in \{0.10, 0.30, 0.40\}</math>.
Já a condição inicial é arbitrária, pois define o estado inicial do sistema que queremos tratar. Fazendo uma referência ao tratamento de sistemas clássicos, seria como definir posição e momento iniciais. É claro que, para ter o sentido físico de uma função de onda, deve-se ter o cuidado de criar uma condição inicial normalizada, satisfazendo 


<math>\int _{-\infty}^{\infty} \, \left| \Psi (x, 0) \right|  ^2 \, dx = 1</math>
Evolução temporal da Energia Total:


bastando, então, inseri-la no programa.


==Implementação em C==
[[Arquivo:time_evolution_U_T_45.png|500px|center]]




Observa-se que o sistema se encontra em equilíbrio a <math>10^5</math> passos  Monte Carlo. Isso é corroborado pelos resultados do artigo de M Rovere et al. A partir do sistema equilibrado, são analisadas a distribuição de densidade e a função distribuição de pares para verificar a presença de uma (líquida) ou duas fases (líquida, gasosa), usando uma média de <math>10^5</math> medidas sucessivas.


Para <math> \rho = 0.10</math> temos  que:


<br />
[[Arquivo: snapshot_t45_d_10.png|300px]] [[Arquivo: P_t45_d_10.png|300px]] [[Arquivo: g_t45_d_10.png|300px]]


<source lang="c">
Observa-se no ''snapshot'' que o sistema parece ter fase líquida e fase gasosa coexistindo. O resultado da distribuição de densidades pode ter sido afetado pelo aumento da área (consequentemente  também das células usadas para realização do histograma) já que o sistema ficou maior a fim de diminuir a densidade mantendo o número de partículas. Apesar disso, o resultado da função distribuição de pares parece razoável com a hipótese da fase líquida estar presente.
#include <stdio.h>
#include <math.h>
#include <complex.h>


#define N_steps 100000
Para <math> \rho = 0.30</math> temos  que:
#define L 1000
#define dt 1
#define dx 4.0
#define w 0.002


double V(int x_rel)
[[Arquivo: snapshot_t45_d_03.png|300px]] [[Arquivo: P_t45_d_03.png|300px]] [[Arquivo: g_t45_d_03.png|300px]]
{
return - pow(w,2) * pow(x_rel - 500,2) / 2.0; //SHO


//return 0;
Esse gráfico foi tratado anteriormente (para temperatura fixa) e está reproduzido novamente a fim de comparação.
}


double complex b(int i)
Para <math> \rho = 0.40</math> temos  que:
{
return 0.5 * dt * I * (1.0 / pow(dx,2.0) + V(i)) + 1.0;
}


double u_0(int x_rel)
[[Arquivo: snapshot_t45_d_40.png|300px]] [[Arquivo: P_t45_d_40.png|300px]] [[Arquivo: g_t45_d_40.png|300px]]
{
//return sqrt(2 / L) * sin(5*x_rel * M_PI / L); //eigenstate of infinite square well


//return (1.0/ sqrt(5 * M_PI)) * exp(I * 0.5 * x_rel) * exp(-pow(x_rel - 500, 2) / (2 * 25)); //gaussian packet
Nesse sistema há a presença de fase líquida com menos espaços vazios. O pico de densidade da fase líquida situa-se em <math> \rho \approx 1.5 </math>. Essa queda da densidade do líquido com aumento da densidade do sistema indica que uma transição de fase para gás pode estar próxima para esses parâmetros.


return (w / M_PI) * exp(-w * pow(x_rel - 500, 2) / 2.0);
Para <math> \rho = 0.70</math> temos que:
}


void CN(double complex *u, double complex *u_aux, double complex *u_next, double complex a, int bi)
[[Arquivo: snapshot_t45_d_70.png|300px]] [[Arquivo: P_t45_d_70.png|300px]] [[Arquivo: g_t45_d_70.png|300px]]
{
//bi = 1 (bounded); bi = 0 (periodic)


int i;
Observa-se que o sistema tem relativa homogeneidade e ocupa toda a área. Isso se caracteriza como um gás, como pode ser verificado na distribuição de densidade, que tem apenas um pico em  <math> \rho \approx 1.0 </math> (que não é a média global). Apesar disso, o sistema parece relativamente organizado, o que sugere uma possível transição de fase. Também sugere isso a distribuição de pares, que tem mínimos baixos que remetem à fase sólida.


for(i = bi; i <= L - bi; i++) u_aux[i] = conj(a) * (u[(i-1+L)%L] + u[(i+1+L)%L]) + conj(b(i)) * u[i];
== Sistema tridimensional ==


//thomas algorithm
Dado um sistema tridimensional com densidade <math> \rho = N / L^3</math> e temperatura <math>T</math>, os diagramas foram feitos com:


double complex c_new[L+1], d_new[L+1];
(A) <math> N = 500 </math> partículas
(B) Cubo de lado <math> L =  (N / \rho)^{1/3} </math> com convenção da imagem mínima
(C) Incialização aleatória
(D) Distância de corte <math>r_c \in \{2.5 \sigma, 3.0 \sigma, 3.5 \sigma \} </math>
(E) Deslocamento máximo inicial <math>\Delta_0 = L/2 </math>


c_new[bi] = a / b(bi);
=== Transição de fase líquido-gás ===
for(i = 1 + bi; i <= L - bi; i++) c_new[i] = a / (b(i) - c_new[i-1] * a);
Conforme o trabalho desenvolvido por Nicholas et al <ref>J.J. Nicolas, K.E. Gubbins, W.B. Streett, and D.J. Tildesley (1979). "Equation of state for
the Lennard-Jones fluid". Mol. Phys., 37</ref>, os valores críticos na transição de fase líquido-gás para um fluido de Lennard-Jones tridimensional são <math>T_{cr} =  1.35</math> e <math>\rho_{cr} = 0.35</math>.


d_new[bi] = u_aux[bi] / b(bi);
=== T = 2.0 (acima do valor crítico)===
for(i = 1 + bi; i <= L - bi; i++) d_new[i] = (u_aux[i] - d_new[i-1] * a) / (b(i) - c_new[i-1] * a);


if (bi == 1) u_next[0] = u_next[L] = 0;
Tendo em vista que a análise no equilíbrio do sistema bidimensional foi feita após <math> 10^5 </math> passos Monte Carlo, é bastante possível que os sistemas tridimensionais aqui tratados não tenham alcançado o equilíbrio. Para  <math>r_c  = 2.5\sigma</math>, temos a evolução temporal da energia total e da pressão:
u_next[L-bi] = d_new[L-bi];
for(i = L-1-bi; i >= bi; i --) u_next[i] = d_new[i] - c_new[i] * u_next[i+1];
//u = u_next


for(i = 0; i <= L; i++) u[i] = u_next[i];
[[Arquivo: t_t20_r25_U.png|500px]] [[Arquivo: t_t20_r25_P.png|500px]]
}


int main(void)
Estimando seus valores a partir do equilíbrio, determina-se o diagrama de fase:
{
int i, j, n = 0;


double complex u[L+1], u_next[L+1], u_aux[L+1], a = - 0.25 * I * dt / pow(dx,2.0);
[[Arquivo:t20_r25_Pd.png|500px|center]]  


//Initial Contition
Além disso, tem-se as diferentes configurações espaciais para as densidades <math>\rho = 0.2</math> e <math>\rho = 0.6</math>, respectivamente:
for(i = 0; i <= L; i ++) u[i] = u_0(i);


while(n < N_steps)
[[Arquivo:g_t20_d02_r25.png|500px]] [[Arquivo:g_t20_d06_r25.png|500px]]
{
CN(u, u_aux, u_next, a, 0);
printf("set title 'Time = %d'\nplot \'-' w lp pt 7 ps 0.1", n);
for(i = 0; i <= L; i ++) printf("%d\t%.10lf\n",i,pow(creal(u[i]),2) + pow(cimag(u[i]),2));
//for(i = 0; i <= L; i ++) printf("%d\t%lf\n",i,cimag(u[i]));
//for(i = 0; i <= L; i ++) printf("%d\t%lf\n",i,creal(u[i]));
printf("e\npause 0.1\n");
n ++;
}


return 0;
Executando o mesmo esquema, o diagrama é comparado para diferentes valores de  <math>r_c</math>:
}


</source><br />
[[Arquivo:t20_Pd.png|500px|center]]


Nessa situação (acima do valor crítico de temperatura), o gráfico pressão - densidade mostra sua injetividade, e portanto, a não coexistência de fases (como esperado).


==Referências==
==Referências==
*Cohen-Tannoudji C., Diu B., Laloe F. Quantum mechanics. Volume 1. Wiley, 1991.
<references />
*Numerical Resolution Of The Schrödinger Equation. Jorgensen L., Lopes Cardozo D., Thivierge E. http://web.pa.msu.edu/people/duxbury/courses/phy480/SchrodingerDynamics.pdf
*Crank, J.; Nicolson, P. (1947). "A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type". Proc. Camb. Phil. Soc. 43 (1): 50–67. doi:10.1007/BF02127704.
*Sherer, Philipp O.J., Computational Physics simulation of Classical and Quantum Systems. Springer, 2010.
*Born M., Nobel lecture: The statistical interpretation of quantum mechanics. 11 de Dezembro de 1954. https://www.nobelprize.org/nobel_prizes/physics/laureates/1954/born-lecture.pdf

Edição atual tal como às 07h48min de 23 de janeiro de 2018

Alunos: Ben-Hur Cardoso e Vinícius Garcia


O potencial devido a interação entre duas partículas separadas por uma distância pode ser modelado pelo potencial de Lennard-Jones [1]:


Posto em unidades reduzidas ( e ), o potencial reduz-se a:


LennardJones.png

As quantidades físicas de interesse (como pressão, tempo e temperatura) possuem ordens de grandeza muito diferentes entre si quando se usa unidades do Sistema Internacional, o que faz com que fique difícil de identificar erros. Por causa disso, é computacionalmente viável trabalhar com o seguintes sistema de unidades básicas, que mantém os valores em baixas ordens de grandeza:

Grandeza Comprimento Tempo Massa Temperatura Energia Pressão Densidade
Unidade

onde é a massa da partícula e é a constante de Boltzmann.

Note que as novas unidades são adimensionais. Assim, para comparar o modelo com um sistema físico real, deve-se converter para as unidades desejadas. É interessante notar que a simulação feita utilizando as unidades reduzidas permite visualizar que sistemas distintos serão equivalentes, já que há infinitas combinações de densidade, temperatura, e com mesmas unidades reduzidas.

Método Monte Carlo

Denomina-se método de Monte Carlo métodos estatísticos que se baseiam em amostragem aleatória massiva para cálculo numérico.

Amostragem simples

O cálculo numérico de uma integral utilizando Monte Carlo parte da ideia de que uma integral pode ser reescrita como:


Dessa forma, utiliza-se amostragem aleatória massiva para estimar , que é a média da função no intervalo de interesse. O estimador para a média real é a média amostral, que por ser uma medida estatística terá sempre um erro inerente associado.

Amostragem por importância

Um problema da amostragem simples é que ela utiliza uma distribuição uniforme, que pode, para uma função que tenha regiões onde seu valor é próximo de zero por exemplo, custar a estimar corretamente o valor médio da função. A fim de contornar isso, podemos utilizar uma distribuição tal que a razão seja o mais constante possível. Reescrevendo a integral:



Basicamente, é a média dos valores da função razão com distribuição para sorteio dos pontos. Esse método diminui a variância amostral e faz com que seja muito mais eficiente a convergência do estimador média amostral para a média real da função.

Algoritmo de Metropolis–Hastings

Dado uma amostra com partículas, a abordagem introduzida por Metropolis (em 1953, algoritmo para distribuições simétricas) e Hastings (em 1970, algoritmo generalizado) segue o seguinte esquema:

(1) Selecionar uma partícula aleatóriamente, e calcular sua energia ;
(2) Dado o deslocamento , calcular ;
(3) Aceitar o movimento  com probabilidade 

Devido ao fato das componentes de poderem assumir valores postivos e negativos (com a nova posição limitada por alguma condição de contorno), após ser aceito a mudança de posição de uma dada partícula, elá poderá ser selecionada (no próximo passo) com a mesma probabilidade e voltar a sua posição original. Isso garante a condição de balanço detalhado. [2]

Estimadores no Equilíbrio

Em todos os exemplos tratados aqui, será usado o ensemble NVT (com o número de partículas, volume e temperatura constantes) com dimensões. Dado isso, os sistemas são caracterizados com um densidade e uma temperatura. Com tais sistemas no equilíbrio, são estimadas (média de sucessivas medidas) a energia total e a pressão, dadas respectivamente por:

 

onde . Além disso, é interessante a análise da capacidade térmica


A função radial de distribuição de pares, é uma função que estima o quão provável é encontrar duas partículas a uma distância dentro de um sistema de várias partículas. Em um sistema de partículas, o é definido como a média do número de partículas a uma distância :


A distribuição de densidades é calculada dividindo espacialmente o sistema, no nosso caso de lado , em células de lado e acumular em um histograma as densidades das células. Essa medida é relevante para identificação de fases no sistema.

Detalhes Técnicos

Condições de Contorno

Qualquer sistema possível de ser feito hoje com método Monte Carlo, apesar do grande poder computacional disponível, fica distante do limite termodinâmico. As condições de contorno podem ser estabelecidas de forma a tentar contornar isso. As condições utilizadas neste trabalho foram condições de contorno periódicas, que possibilitam que o sistema se comporte como se fosse muito maior do que é, desde que ele seja isotrópico.

Convenção da imagem mínima

Construindo um sistema cúbico de lado L, com um número N de partículas fixo, as condições de contorno periódicas fazem com a simulação não se limite ao cubo. Ao calcular a interação entre duas partículas, portanto, a distância entre elas não é unívoca, já que há incontáveis cópias de cada uma, repetindo-se periodicamente. Resolve-se esse problema utilizando-se da convenção da imagem mínima, em que é calculada a interação com a imagem mais próxima das outras.

Truncamento nas interações

Um problema da condição de contorno periódica é que, a princípio, cada partícula interagiria com todas as outras do sistema, que devido ao fato de ser periodicamente repetido, seriam infinitas. Como o potencial utilizado é de curto alcance, deve ser possível, de alguma forma, limitar as interações entre as partículas sem perda dos significados numéricos da simulação. Pode-se, então, truncar as interações de uma partícula a partir de uma distância de corte de forma que há uma descontinuidade no potencial dessa partícula na esfera de raio .


Mic.png


Assim o potencial simulado é:


Esse é o chamado truncamento simples. O subíndice LJ se refere ao potencial de Lennard-Jones. O problema é que isso cria uma contribuição indesejada à pressão, já que há uma força de impulso por conta dessa descontinuidade do potencial.

Pode-se, então, usar o truncamento com deslocamento, que evita essa descontinuidade fazendo uma subtração em todo ponto do módulo do potencial de Lennard-Jones à distância . Assim, o potencial simulado é:


Com esse potencial, que é contínuo, as forças serão sempre finitas, o que retira a força de impulso que alteraria medidas de pressão.

Translação

A possível nova posição de uma partícula será tal que

  
  
  

sendo números aleatórios uniformemente distribuídos e o deslocamento máximo permitido. Se a nova posição for aceita seguindo o algorítmo de Metropolis, a partícula assumirá essa posição. Caso contrário, a antiga configuração contribuirá novamente na média. [3]

O tamanho de não pode ser grande ao ponto de nunca ser aceito nem pequeno ao ponto de sempre ser aceito. Dado isso, buscando uma taxa de aceitação entre e , define-se um inicial que vai se ajustando no decorrer da simulação, executando - a cada passo Monte Carlo ( passos de simulação) - o seguinte algorítimo:

 Se ():   
 Se (): 

Sistema bidimensional

Dado um sistema com densidade e temperatura , os sistemas foram feitos com:

(A)  partículas
(B) Quadrado de lado , utilizando a convenção da imagem mínima
(C) Incialização aleatória
(D) Distância de corte 
(E) Deslocamento máximo inicial 

Transição de fase líquido-gás

A partir do artigo de M Rovere et al [4], que modelou o sistema apresentado, obtemos a informação da temperatura e da densidade críticas. Os valores encontrados pelos autores do artigo são e , que determinam qual intervalo de valores de temperatura e densidade apresentam coexistência de fases distintas.

Transição de fase com densidade fixa

Fixada a densidade , foram realizadas simulações para as diferentes temperaturas .

Evolução temporal da Energia Total:

Time evolution U d 03.png

Observa-se que o sistema se encontra em equilíbrio a passos Monte Carlo. Isso é corroborado pelos resultados do artigo de M Rovere et al. A partir do sistema equilibrado, são analisadas a distribuição de densidade e a função distribuição de pares para verificar a presença de uma (líquida) ou duas fases (líquida, gasosa), usando uma média de medidas sucessivas.

Para temos que:

Snapshot t70 d 03.png P t70 d 03.png G t70 d 03.png

Nesse caso observa-se apenas uma fase, com a densidade média sendo a densidade global do sistema, que é diferente do centro do pico (). Isso caracteriza a fase gasosa, que ocupa todo o espaço disponível. A função distribuição de pares rapidamente converge para 1 com pequenas oscilações, caracterizando um gás de alta densidade.

Para temos que:

Snapshot t50 d 03.png P t50 d 03.png G t50 d 03.png

Nesse caso observa-se a coexistência das fases gasosa e líquida. Além disso, pode-se observar a densidade característica aproximada de cada fase a essa temperatura centrada nos respectivos picos ( e ). A função distribuição de pares representa a de uma fase líquida.

Para temos que:

Snapshot t45 d 03.png P t45 d 03.png G t45 d 03.png

Nesse gráfico a primeira impressão é de coexistência de fases. Apesar disso, densidade zero não corresponde a uma fase real, e sim aos espaços vazios. Dessa forma, a fase aqui é a líquida. O snapshot permite visualizar isso, já que há aglomerados de densidade aparentemente constante separados por espaços praticamente sem partículas. A densidade da fase líquida a essa temperatura é aproximadamente o centro do pico da direita ().

Para temos que:

Snapshot t30 d 03.png P t30 d 03.png G t30 d 03.png

Assim como no gráfico anterior, há dois picos mas apenas um representa uma fase. A principal diferença é que o pico em zero é bem maior e mais agudo que no caso , o que ilustra que há mais espaço vazio. A densidade do líquido aqui é . Nota-se que o sistema parece estar se organizando numa fase hexagonal, apesar de ainda ser amorfa, e a densidade já é bem mais alta. A distribuição de pares tem mínimos bem menores do que para os parâmetros acima tratados, caracterizando um sistema que parece estar próximo da transição para a fase hexagonal de baixa energia.

De maneira geral, é possível estimar o intervalo de parâmetros para o qual ocorre a transição de fase líquido-gás. É possível, também, visualizar, a partir da distribuição de densidades, a coexistência de fases, já que cada estado tem uma densidade média característica a dada temperatura.

Transição de fase com temperatura fixa

Dado a temperatura fixa , é observado o sistema para as diferentes densidades .

Evolução temporal da Energia Total:


Time evolution U T 45.png


Observa-se que o sistema se encontra em equilíbrio a passos Monte Carlo. Isso é corroborado pelos resultados do artigo de M Rovere et al. A partir do sistema equilibrado, são analisadas a distribuição de densidade e a função distribuição de pares para verificar a presença de uma (líquida) ou duas fases (líquida, gasosa), usando uma média de medidas sucessivas.

Para temos que:

Snapshot t45 d 10.png P t45 d 10.png G t45 d 10.png

Observa-se no snapshot que o sistema parece ter fase líquida e fase gasosa coexistindo. O resultado da distribuição de densidades pode ter sido afetado pelo aumento da área (consequentemente também das células usadas para realização do histograma) já que o sistema ficou maior a fim de diminuir a densidade mantendo o número de partículas. Apesar disso, o resultado da função distribuição de pares parece razoável com a hipótese da fase líquida estar presente.

Para temos que:

Snapshot t45 d 03.png P t45 d 03.png G t45 d 03.png

Esse gráfico foi tratado anteriormente (para temperatura fixa) e está reproduzido novamente a fim de comparação.

Para temos que:

Snapshot t45 d 40.png P t45 d 40.png G t45 d 40.png

Nesse sistema há a presença de fase líquida com menos espaços vazios. O pico de densidade da fase líquida situa-se em . Essa queda da densidade do líquido com aumento da densidade do sistema indica que uma transição de fase para gás pode estar próxima para esses parâmetros.

Para temos que:

Snapshot t45 d 70.png P t45 d 70.png G t45 d 70.png

Observa-se que o sistema tem relativa homogeneidade e ocupa toda a área. Isso se caracteriza como um gás, como pode ser verificado na distribuição de densidade, que tem apenas um pico em (que não é a média global). Apesar disso, o sistema parece relativamente organizado, o que sugere uma possível transição de fase. Também sugere isso a distribuição de pares, que tem mínimos baixos que remetem à fase sólida.

Sistema tridimensional

Dado um sistema tridimensional com densidade e temperatura , os diagramas foram feitos com:

(A)  partículas
(B) Cubo de lado  com convenção da imagem mínima
(C) Incialização aleatória
(D) Distância de corte 
(E) Deslocamento máximo inicial 

Transição de fase líquido-gás

Conforme o trabalho desenvolvido por Nicholas et al [5], os valores críticos na transição de fase líquido-gás para um fluido de Lennard-Jones tridimensional são e .

T = 2.0 (acima do valor crítico)

Tendo em vista que a análise no equilíbrio do sistema bidimensional foi feita após passos Monte Carlo, é bastante possível que os sistemas tridimensionais aqui tratados não tenham alcançado o equilíbrio. Para , temos a evolução temporal da energia total e da pressão:

T t20 r25 U.png T t20 r25 P.png

Estimando seus valores a partir do equilíbrio, determina-se o diagrama de fase:

T20 r25 Pd.png

Além disso, tem-se as diferentes configurações espaciais para as densidades e , respectivamente:

G t20 d02 r25.png G t20 d06 r25.png

Executando o mesmo esquema, o diagrama é comparado para diferentes valores de :

T20 Pd.png

Nessa situação (acima do valor crítico de temperatura), o gráfico pressão - densidade mostra sua injetividade, e portanto, a não coexistência de fases (como esperado).

Referências

  1. Lennard-Jones, J. E. (1924), "On the Determination of Molecular Fields", Proc. R. Soc. Lond. A, 106 (738): 463–477
  2. Daan Frenkel, Berend Smit, "Understanding molecular simulation: from algorithms to applications"
  3. David P. Landau, Kurt Binder, "A Guide to Monte Carlo Simulations in Statistical Physics"
  4. M Rovere et al (1990) "The gas-liquid transition of the two-dimensional Lennard-Jones fluid". J. Phys.: Condens. Matter 2 7009
  5. J.J. Nicolas, K.E. Gubbins, W.B. Streett, and D.J. Tildesley (1979). "Equation of state for the Lennard-Jones fluid". Mol. Phys., 37