Grupo - Lennard Jones

O potencial devido a interação entre duas partículas separadas por uma distância $r'$ pode ser modelado pelo potencial de Lennard Jones :

$U'(r')=4\epsilon \left[\left({\frac {\sigma }{r'}}\right)^{12}-\left({\frac {\sigma }{r'}}\right)^{6}\right]$

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

$U(r)=4\left[\left({\frac {1}{r}}\right)^{12}-\left({\frac {1}{r}}\right)^{6}\right]$

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, é conveniente 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	$\sigma$	$\sigma {\sqrt {m_{p}/\epsilon }}$	$m_{p}$	$\epsilon /k_{B}$	$\epsilon$	$\epsilon /\sigma ^{3}$	$1/\sigma ^{3}$

onde $m_{p}$ é a massa da partícula e $k_{B}$ é 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, $\epsilon$ e $\sigma$ 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

Pode-se querer calcular uma integral numericamente utilizando Monte Carlo. Uma forma de fazer isso parte de que uma integral pode ser reescrita como:

$F=\int _{a}^{b}{f(x)dx}=(b-a)\langle f(x)\rangle$

Dessa forma, utiliza-se amostragem aleatória massiva para estimar $\langle f(x)\rangle$, que é a média da função no intervalo de interesse.

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, custar a estimar corretamente o valor médio da função. A fim de driblar isso, podemos utilizar uma distribuição $w(x)$ tal que a razão ${\frac {f(x)}{w(x)}}$ seja o mais constante possível. Reescrevendo a integral:

$F=\int _{a}^{b}{{\frac {f(x)}{w(x)}}w(x)dx}$

Usando isso, faz-se a mesma integral por um processo um pouco diferente; a média, agora, será dos valores da função razão ${\frac {f(x)}{w(x)}}$ com distribuição $w(x)$ 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

Dado uma amostra com $N$ partículas, a abordagem introduzida por Metropolis segue o seguinte esquema:

(1) Selecionar uma partícula aleatóriamente, e calcular sua energia $U(r)$;
(2) Dado o deslocamento $\mathbf {r_{n}} =\mathbf {r} +\mathbf {\Delta }$, calcular $U(r_{n})$;
(3) Aceitar o movimento $\mathbf {r} \rightarrow \mathbf {r_{n}}$ com probabilidade $p=min\{1;\exp[-\beta (U(r_{n})-U(r))]\}$

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). Dado isso, os sistemas são caracterizados com um densidade e uma temperatura. Com tais sistema no equilíbrio, são estimadas (média de sucessivas medidas) a energia total e a pressão, dadas respectivamente por:

$U_{T}=\sum _{i}\sum _{j>i}U(r_{ij})$ 

$P={\frac {\rho }{\beta }}+{\frac {2}{3V}}\sum _{i}\sum _{j>i}\mathbf {f(\mathbf {r_{ij}} )} \cdot \mathbf {r_{ij}}$

onde $\mathbf {r_{ij}} =\mathbf {r_{j}} -\mathbf {r_{i}}$. Além disso, é interessante a análise da capacidade térmica

$C_{v}={\frac {\langle U_{T}^{2}\rangle -\langle U_{T}\rangle ^{2}}{k_{B}T^{2}}}$

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 driblar 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 seja isotrópico.

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 $i$ a partir de uma distância de corte $r_{c}$ de forma que há uma descontinuidade no potencial dessa partícula na esfera de raio $r_{c}$. Assim o potencial simulado é:

$U_{trunc}(r)={\begin{cases}U_{LJ}(r)&,r\leq r_{c}\\0&,r>r_{c}\end{cases}}$

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 $r_{c}$. Assim, o potencial simulado é:

$U_{deslocado}(r)={\begin{cases}U_{LJ}(r)-U_{LJ}(r_{c})&,r\leq r_{c}\\0&,r>r_{c}\end{cases}}$

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á $(x_{n},y_{n},z_{n})$ tal que

 $x_{n}=x+\Delta (\alpha -0.5)$ 
 $y_{n}=y+\Delta (\alpha -0.5)$ 
 $z_{n}=z+\Delta (\alpha -0.5)$ 

sendo $\alpha$ números aleatórios uniformemente distribuídos e $\Delta$ o deslocamento máximo permitido. Se a nova posição for aceita seguindo o algorítmo de Metrópolis, a partícula assumirá essa posição.

O tamanho de $\Delta$ 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 $\eta$ entre $30\%$ e $50\%$, define-se um $\Delta _{0}$ inicial que vai se ajustando no decorrer da simulação, executando - a cada passo Monte Carlo - o seguinte algorítimo:

 Se ($\eta <0.3$): $\Delta \leftarrow \Delta \times (\eta /0.3)$  
 Se ($\eta >0.5$): $\Delta \leftarrow \Delta \times (\eta /0.5)$

Diagramas de fase

Dado um sistema tridimensional com densidade $\rho$ e temperatura $T$, os diagramas foram feitos com:

(A) $N=500$ partículas
(B) Cubo de lado $L=(N/\rho )^{1/3}$ com condições de contorno periódicas
(C) Incialização aleatória
(D) Distância de corte $r_{c}\in \{2.5\sigma ,3.0\sigma ,3.5\sigma \}$
(E) Deslocamento máximo inicial $\Delta _{0}=L/2$

T = 2.0 (acima do valor crítico)

Para $r_{c}=2.5\sigma$, temos a evolução temporal da energia total e da pressão:

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

Além disso, tem-se as diferentes configurações espaciais para as densidades $\rho =0.2$ e $\rho =0.6$, respectivamente:

Executando o mesmo esquema, são comparados os diagramas para diferentes valores de $r_{c}$:

Nessa situação (acima do valor crítico de temperatura), o gráfico pressão - densidade (esquerda) mostra sua injetividade, e por tanto, a não coexistência de fases (como esperado). O gráfico capacidade térmica - densidade (direita) é monotonicamente crescente, por consequência.

T = 0.9 (abaixo do valor crítico)

Para $r_{c}=2.5\sigma$, temos a evolução temporal da energia total e da pressão:

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

Além disso, tem-se as diferentes configurações espaciais para as densidades $\rho =0.2$ e $\rho =0.6$, respectivamente:

Executando o mesmo esquema, são comparados os diagramas para diferentes valores de $r_{c}$:

Nessa situação (abaixo do valor crítico de temperatura), o gráfico pressão - densidade (esquerda) mostra sua não-injetividade, e por tanto, a coexistência de fases líquido-vapor, apresentando metaestabilidade com pressões negativas na região. Esse efeito é devido ao tamanho pequeno da amostra, que tem um alto custo de energia livre para a criação de uma interface de vapor de líquido, por isso se torna pouco aconselhavel utilisar o ensemble NVT para diagramas de coexistência de fases.