Termostato de Andersen

INTRODUÇÃO

A dinâmica molecular é uma técnica que naturalmente simula sistemas clássicos compostos por N partículas interagindo dentro de um volume V. Nesse contexto, as posições das partículas são atualizadas com base no potencial de interação escolhido. Sob a suposição de ergodicidade — ou seja, que as médias temporais equivalem às médias de ensemble —, as simulações resultam em amostragens do ensemble microcanônico (NVE). Nesse ensemble, o número de partículas N, o volume V, e a energia total E permanecem constantes (aproximadamente).

Ao colocar um sistema em contato com um reservatório térmico a uma temperatura T, mudamos do ensemble microcanônico (NVE), onde a energia é mantida constante, para o ensemble canônico (NVT), no qual a temperatura do sistema é constante. Nesse novo ensemble, a distribuição de probabilidade das velocidades das partículas segue a forma da distribuição de Maxwell-Boltzmann, uma distribuição gaussiana associada à temperatura T. Um dos métodos mais simples para realizar uma amostragem correta do ensemble canônico é o termostato de Andersen. Neste estudo, focaremos na análise desse termostato, explorando sua implementação, características e aplicação na simulação de sistemas termodinâmicos.

FUNDAMENTO TEÓRICO E DETALHES DA IMPLEMENTAÇÃO

Condições de Contorno Periódicas

As condições de contorno periódicas são fundamentais para garantir que simulações computacionais de sistemas físicos representem com precisão o comportamento de sistemas grandes e infinitos. Elas ajudam a evitar efeitos de borda que não são representativos do comportamento real de partículas em um espaço muito grande, permitindo simulações mais realistas e precisas.

É utilizada a convenção da imagem mínima, que calcula a menor distância entre as partículas, sendo essa sempre menor ou igual a $L/2$.

Integração de Velocity-Verlet

Escolhida principalmente por não ser de $O(N^{2})$ como o RK2, por exemplo. Precisamos calcular apenas 1 vez a força em cada passo do algoritmo. Assim, é o melhor "custo-benefício" dentre os possíveis algoritmos, possibilitando precisão suficiente.

Condições Iniciais

Em todas as simulações, as partículas foram inicialmente posicionadas em uma grade cúbica uniformemente espaçada. Essa abordagem foi adotada para evitar que as partículas começassem muito próximas umas das outras, o que poderia gerar erros devido ao rápido aumento do potencial de Lennard-Jones quando as partículas estão muito próximas.

As velocidades foram escolhidas a partir de uma distribuição uniforme no intervalo $[-0.5,0.5]$, sendo reescaladas para a temperatura desejada.

Potencial de Lennard-Jones

O potencial de Lennard-Jones é amplamente utilizado para descrever as interações entre partículas em simulações de dinâmica molecular. Ele captura de forma simples e eficiente o comportamento atrativo de longo alcance e o repulsivo de curto alcance, que são características comuns em sistemas interatômicos e intermoleculares.

Matematicamente, o potencial de Lennard-Jones é dado por:

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

onde:

• $r$ é a distância entre dois centros de partículas,
• $\varepsilon$ é a profundidade do poço de potencial, representando a intensidade da interação atrativa,
• $\sigma$ é a distância na qual o potencial é igual a zero.

Nas simulações deste estudo foi utilizado o potencial de Lennard-Jones truncado para otimização, onde $r_{c}=2.5\sigma$ foi escolhido como o raio de corte, a partir do qual o potencial é considerado zero. Além disso, para evitar descontinuidades, o potencial também é deslocado, temos então: $U_{\text{truncado e deslocado}}(r)={\begin{cases}4\varepsilon \left[\left({\frac {\sigma }{r}}\right)^{12}-\left({\frac {\sigma }{r}}\right)^{6}\right]-U(r_{c}),&{\text{se }}r\leq r_{c},\\0,&{\text{se }}r>r_{c}.\end{cases}}$

Unidades Adimensionais

Para simplificar e possibilitar a comparação entre diferentes sistemas, todas as grandezas foram convertidas para unidades adimensionais, normalizando com base na massa e nos parâmetros do potencial:

• $m$
• $\sigma$
• $\epsilon$

Para facilitar os cálculos, em todas as simulações foram adotados $m=1$, $\sigma =1$ e $\epsilon =1$.

Distribuição de Maxwell-Boltzmann

A distribuição de Maxwell-Boltzmann descreve como as velocidades das partículas em um sistema no ensemble canônico estão distribuídas em equilíbrio térmico. Para um sistema de partículas com temperatura $T$, a função de distribuição de Maxwell-Boltzmann $f(v)$ para a velocidade $v$ de uma partícula é dada pela seguinte equação:

$f(v)=4\pi \left({\frac {m}{2\pi k_{B}T}}\right)^{3/2}v^{2}e^{-{\frac {mv^{2}}{2k_{B}T}}}$

Temperatura

A definição de temperatura em um sistema clássico em função das velocidades das partículas pode ser obtida pela equipartição de energia:

$\left\langle {\frac {1}{2}}mv_{\alpha }^{2}\right\rangle ={\frac {1}{2}}k_{B}T$

Temos então:

$T(t)={\frac {1}{k_{B}N_{f}}}\sum _{i=1}^{N_{f}}m_{i}v_{i}^{2}$

onde $N_{f}$ é o número de graus de liberdade do sistema. No ensemble microcanônico, como o momento é conservado, temos $N_{f}=3N-3$, mas no ensemble canônico, e portanto no termostato de Andersen, o momento não é conservado e utilizamos $N_{f}=3N$.

TERMOSTATO DE ANDERSEN

Para simular o contato do sistema com um reservatório térmico, a cada passo de tempo Δt é realizado um procedimento de Monte Carlo. Nesse processo, algumas partículas (nesse estudo foi escolhido $N$, ou seja, cada partícula é, em media, sorteada uma vez a cada passo) são selecionadas aleatoriamente, uma de cada vez, e suas velocidades são atualizadas. Essas novas velocidades são sorteadas a partir da distribuição gaussiana centrada em zero com $\sigma ={\sqrt {T}}$. Essa abordagem introduz um elemento de estocasticidade ao modelo, permitindo que a temperatura média do sistema flutue em torno da temperatura T do reservatório.

Na prática, a interação entre o sistema e o reservatório é controlada pela frequência de colisão $\nu$. Para cada partícula, é gerado um número aleatório $r$ no intervalo [0,1]. Caso $r<\nu *\Delta t$ a velocidade da partícula é atualizada conforme descrito.

Algoritmo

Podemos descrever o algoritmo em 4 passos:

1. Inicia-se com um conjunto de posições e velocidades.

2. Integra-se as equações do movimento para um passo $\Delta t$.

3. N partículas são selecionadas para colidir com o reservatório térmico.

4. Para cada partícula selecionada, definir nova velocidade a partir da distribuição de Maxwell-Boltzmann correspondente à temperatura T do reservatório.

RESULTADOS

Para a análise, foram escolhidas 3 temperaturas diferentes. Para cada uma delas, foi simulado $N=108$ partículas, com densidade $p=0.8442$, $\nu =0.5$, $\Delta t=0.001$ e número de passos $n=50.000$. Para checar de forma mais consistente se o sistema no equilíbrio realmente corresponde ao ensemble canônico, foi comparada a curva teórica da distribuição de Maxwell-Boltzmann para a temperatura correspondente, com um histograma das velocidades em todas as dimensões obtidas nos últimos 40.000 passos (os primeiros 10.000 passos são o transiente para o equilíbrio). Na figura abaixo, vemos que houve uma aproximação muito boa da distribuição de velocidades da simulação com os valores esperados pela distribuição de Maxwell-Boltzmann.

Figura 1: Comparação da distribuição de Maxwell-Boltzmann (linha sólida) com a distribuição de velocidades obtida nos últimos 40.000 passos da simulação (círculos) para diferentes temperaturas.

Além disso, são mostradas abaixo a série temporal para cada uma das temperaturas, com a média e desvio padrão da temperatura dos últimos $40.000$ passos.

T=1

$\mu =1.011$ $\sigma =0.067$

Figura 1: Evolução da temperatura durante a simulação com temperatura do reservatório $T=1$.

T=2

$\mu =2.025$ $\sigma =0.166$

Figura 3: Evolução da temperatura durante a simulação com temperatura do reservatório $T=2$.

T=3

$\mu =3.027$ $\sigma =0.247$

Figura 5: Evolução da temperatura durante a simulação com temperatura do reservatório $T=3$.

Análise do parâmetro da frequência de colisão

A frequência de colisão $\nu$ controla o nível de interação entre o sistema e o reservatório térmico. Quando $\nu$ é muito baixa, o número de colisões entre o sistema e o reservatório em um intervalo de tempo é pequeno, exigindo um tempo maior de simulação para equilibrar a temperatura. Nesse caso, o sistema permanece por períodos mais longos em estados que não refletem o equilíbrio térmico desejado, o que prejudica a convergência para o ensemble canônico. Na figura abaixo vemos que a flutuação parece estar se aproximando lentamente de $T=1.0$, com mais tempo de simulação parece que essa temperatura será alcançada.

Figura 6: Evolução temporal da temperatura para $\nu =0.01$ e $T=1$.

Por outro lado, valores muito altos de $\nu$ podem introduzir um excesso de estocasticidade no modelo, fazendo com que a interação entre as partículas percam importância. Esse regime pode ser indesejado em simulações onde as interações interpartículas são críticas para o comportamento do sistema. Além disso, o termostato de Andersen introduz descontinuidades nas velocidades e, por consequência, no momento, sendo assim, não é possível medir propriedades de transporte, como difusão e viscosidade.

Na figura abaixo é visto claramente que a convergência demorou mais passos do que no caso de $\nu =0.5$.

Figura 6: Evolução temporal da temperatura para $\nu =0.1$ e $T=1$.

CONCLUSÃO

Este estudo implementou e analisou o termostato de Andersen como ferramenta para simular sistemas no ensemble canônico (NVT). A estocasticidade introduzida pela atualização aleatória das velocidades das partículas, baseada na distribuição de Maxwell-Boltzmann, mostrou-se eficiente para reproduzir a temperatura do reservatório térmico e as propriedades estatísticas de equilíbrio termodinâmico. Os resultados demonstraram a acurácia do método. A precisão diminui com o aumento da temperatura, como esperado para essa distribuição.

No entanto, foi observado um aumento das flutuações em temperaturas mais altas, como evidenciado pelos desvios padrão maiores. Isso é consistente com o comportamento esperado para sistemas térmicos e reflete o impacto da energia térmica sobre as velocidades das partículas.

Ao longo das simulações realizadas, foi necessário utilizar um valor de frequência de colisão $\nu =0.5$. Para valores menores de $\nu$, observou-se que o sistema não convergia adequadamente para a temperatura do reservatório térmico, apresentando desvios significativos em relação ao valor esperado. Esse comportamento pode ser explicado pelo $\Delta t$ escolhido. Se fosse menor, apesar de as colisões demorarem mais para acontecer, a dinâmica também seria mais lenta, e talvez o valor de $\nu$ pudesse ser menor.

Em suma, o termostato de Andersen se mostrou uma ferramenta valiosa e relativamente simples para a simulação de sistemas interagindo com um reservatório térmico. A abordagem permite a obtenção de propriedades relevantes no ensemble canônico e abre caminho para estudos mais aprofundados, como medições de outras propriedades termodinâmicas, a comparação com outros termostatos ou a aplicação a sistemas mais complexos.

código utilizado