O Potencial de Lennard-Jones

Introdução

O modelo para o Potencial de Lennard-Jones foi proposto por Sir John Edward Lennard-Jones, por volta de 1924, e tem como objetivo descrever a energia potencial de interação entre duas partículas não covalentes, que são átomos ou moléculas que não realizam ligações químicas com compartilhamento de elétrons , mas que exercem forças atrativas (dipolo-dipolo, dipolo-dipolo induzido e interações de London) e forças repulsivas uns sobre os outros .

A equação do Potencial Lennard-Jones considera tanto forças atrativas quanto repulsivas, sendo descrito como:

$V_{LJ}(r)=4ϵ[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6]$

Nesta equação, os parâmetros possuem os seguintes significados físicos:

• $ϵ$ (Épsilon): Representa a profundidade do poço de potencial, ou seja, a força da atração entre as partículas. Corresponde à energia mínima do sistema.
• $\sigma$ (Sigma): É a distância finita na qual o potencial interpartículas é nulo. Fornece uma estimativa do diâmetro efetivo do átomo ou molécula.
• O termo ${\left(\frac{\sigma }{r}\right)}^{12}$: Modela a repulsão de curto alcance, originada pelo Princípio de Exclusão de Pauli, que impede a sobreposição das nuvens eletrônicas.
• O termo ${\left(\frac{\sigma }{r}\right)}^6$: Modela a atração de longo alcance (forças de dispersão de London ou de van der Waals), originada pela interação entre dipolos instantâneos induzidos.

As aplicações do Potencial Lennard-Jones são várias, mas a mais conhecida delas é para moléculas diatômicas, como o oxigênio e o hidrogênio. Se a força resultante sobre dois átomos inicialmente separados (ionizados ou neutros) impede a separação deles, pode-se formar uma molécula simples. Considerando o caso em que dois átomos neutros estão à uma distância muito maior que o raio deles, há uma Força de Van de Waals induz um momento dipolar não nulo devido ao movimento dos elétrons ao redor do núcleo, atraindo os dois átomos. Ao se aproximarem, a força de repulsão eletrostática devida às cargas de mesmo sinal começa a aumentar sua intensidade. Ao final da ação destas duas forças, uma molécula é formada (desprezamos deformações da nuvem de elétrons, pela aproximação de Bohr-Oppenheimer).

Como estas duas forças são conservativas, pois só dependem da distância, podemos associar um potencial $V_{LJ}(r)$ pela relação:

$F(r)=-\mathrm{\nabla }{V}_{LJ}(r)$

Um potencial que satisfaz essas condições foi proposto por Lennard-Jones e é conhecido como o Potencial Lennard-Jones.

Objetivos

Com isso, agora buscamos estudar alguns parâmetros de um sistema simular um gás sob efeito das forças de interação entre átomos utilizando uma rede bidimensional com N partículas. Inicialmente, pretendemos implementar essa simulação para calcular a energia e o parâmetro de ordem g(r), a fim de simular esse sistema para diferentes temperaturas e identificar as fases do sistema. Em seguida, vamos encontrar o tempo de correlação pela função de autocorrelação e depois analisar a série temporal da energia para calcular as barras de erros. Por fim, vamos aplicar o Finite-Size Scaling (FSS) para compreender o comportamento real do sistema no limite termodinâmico.

Metodologia

Para realizar as simulações do Potencial de Lennard-Jones, simulamos uma rede bidimensional quadrada com N partículas que se deslocam na rede de forma aleatória, com condições de contorno periódicas e vizinhança esférica. Para isso, utilizamos métodos de Monte Carlo, em que a cada passo sorteamos N partículas (uma de cada vez) e calculamos a sua energia local. Então, escolhemos aleatoriamente outra posição para a partícula e, se a nova energia local for menor, trocamos a sua posição para a nova.

Todo o código foi realizado em Python3, no ambiente do Google Colab, com o uso das bibliotecas Numpy, Matplotlib, Numba e pode ser acessado no link a seguir:

https://colab.research.google.com/drive/1PeGA2pAOo4MMLUyam2wevbuhHsaLrJAm?usp=sharing

Resultados e Discussão

Na seção a seguir, discutimos os resultados de uma única amostra, com semente aleatória.

Evolução Temporal e Termalização da Energia

O monitoramento da energia potencial total do sistema em função dos passos de Monte Carlo (MCS) é fundamental para garantir que o sistema atingiu o equilíbrio termodinâmico antes da coleta de dados.

Evolução da energia potencial total do sistema em unidades reduzidas. A linha vermelha tracejada demarca o fim do regime de termalização e o início da fase de produção.

O gráfico de evolução da energia apresenta duas regiões físicas distintas:

• Regime Transiente (Termalização): Nos passos iniciais da simulação (antes da linha tracejada vermelha marcando 5000 MCS), observa-se uma queda abrupta e colossal da energia. Isso ocorre porque as posições iniciais das partículas foram geradas aleatoriamente no espaço contínuo da caixa bidimensional. Consequentemente, algumas partículas "nascem" sobrepostas ou demasiadamente próximas umas das outras. Devido ao forte termo repulsivo $r^{-12}$ do potencial de Lennard-Jones, essa proximidade artificial gera uma energia positiva tendendo ao infinito. O algoritmo de Metropolis, buscando minimizar a energia, atua rejeitando as sobreposições e afastando essas partículas violentamente, o que causa o despencar imediato da curva de energia.
• Regime de Produção (Equilíbrio): Após o período de termalização, o sistema relaxa. A energia estabiliza-se e passa a variar em torno de um valor médio negativo. Estas variações visíveis não são ruído numérico ou erro do código, mas sim as flutuações térmicas naturais inerentes a um sistema no ensemble canônico. A magnitude dessas flutuações está intimamente ligada a propriedades macroscópicas do fluido, fornecendo a base para o cálculo da capacidade térmica a volume constante ($C_v$). Apenas os dados situados nesta região de equilíbrio são utilizados para o cálculo de médias e propriedades estruturais.

Estrutura do Fluido: Função de Distribuição Radial $g(r)$

A Função de Distribuição Radial, $g(r)$, é um dos observáveis microscópicos mais importantes obtidos na simulação, pois revela a estrutura espacial local do sistema, calculando a probabilidade de encontrar partículas vizinhas a uma dada distância da partícula de referência em relação a um gás ideal.

Função de Distribuição Radial $g(r)$ do fluido simulado. O pico pronunciado denota a primeira camada de coordenação típica de uma fase líquida.

O gráfico de $g(r)$ obtido para a densidade ${\rho }^{*}=0.5$ e temperatura $T^{*}=0.5$ demonstra claramente um estado fortemente estruturado, análogo a uma fase líquida, caracterizada por três regiões específicas:

• Caroço Repulsivo (Volume Excluído): Para distâncias curtas ($r^{*}<0.9$), observa-se que $g(r)=0$. Isso comprova fisicamente a inacessibilidade dessa região: é a manifestação da repulsão de curto alcance, evidenciando o diâmetro efetivo das partículas onde a sobreposição das nuvens eletrônicas é proibida.
• Ordem de Curto Alcance (Camadas de Solvatação): O sistema apresenta um pico primário extremamente pronunciado em torno de $r^{*}\approx 1.12$. Este valor coincide analiticamente com o mínimo do poço atrativo do potencial de Lennard-Jones ($2^{1/6}\sigma$). Este pico representa a primeira camada de coordenação, ou seja, a "casca" de vizinhos mais próximos que circundam uma partícula. Os picos subsequentes (perto de $r^{*}=2.2$ e $3.2$), com amplitudes sucessivamente menores, representam a segunda e terceira camadas de vizinhos. A presença dessas variações indica uma forte ordem de curto alcance estrutural, uma assinatura de fluidos líquidos.
• Limite de Longo Alcance: Para grandes distâncias, as interações intermoleculares decaem e as correlações espaciais desaparecem. A função suaviza-se e tende assintoticamente a $g(r)=1$ (linha tracejada cinza), que é o comportamento estatístico homogêneo esperado para distâncias macroscopicamente grandes. As pequenas flutuações observadas nessa região se devem ao tamanho finito do sistema simulado e das estatísticas de amostragem finitas das partículas.

Determinação do Tempo de Correlação

O tempo de correlação de uma simulação é o tempo necessário entre duas configurações diferentes do nosso sistema para que sejam independentes. Para fazer esse cálculo, usamos a Função de Autocorrelação C(t) que mede o quanto a configuração atual do sistema se assemelha com a configuração há t passos. Por definição, temos que no instante t=0 a função é C(0) = 1, onde sua correlação é máxima, enquanto que o instante t em que C(t) = 1/e (aproximadamente 0.36) é o valor de t para o tempo de correlação que procuramos.

Observando o gráfico, vemos que o valor de $t\approx 102MCS$ corresponde a o C(t) = 1/e, encontrando o tempo de correlação da simulação. Com esse valor, sabemos que a cada pelo menos 103 passos para o sistema perder a maior parte da sua correlação temporal.

Contudo, deve-se esclarecer que para obter mais precisão ainda nas nossas medidas devemos escolher amostras com $2t$ de "distância".

Análise de Blocos

A análise de blocos é um método utilizado para encontrarmos o erro real da nossa simulação através da divisão das amostras em blocos. Dessa forma, podemos calcular a média da energia de cada bloco de amostras e utilizar como uma nova amostra, agora independente das outras. Assim, podemos calcular a variância dessas novas amostras e encontrar a variância real do sistema.

Vemos que o valor da variância real é de ${\sigma }_{Real}^2\approx 25000$. No código, calculamos que o erro real é $\approx 1.589$, visto que lidamos com 10000 passos salvos e também $err{o}_{Real}=({\sigma }_{Real}^2/N_{dados}{)}^{1/2}$. Enquanto que o erro ingênuo (o erro calculado pelo desvio padrão) foi de $\approx 0.101$, o que deixa claro que as amostra são dependentes entre si e que não podemos usar esse método para calcular a variância e o erro.

Finite-Size Scaling (FSS)

O Escalonamento de Tamanho Finito é uma técnica que utilizamos para compreender o comportamento de um sistema finito no limite termodinâmico. Para isso, repetimos a simulação para diferentes números de partículas (N = 100, 256, 500, 1000) e calculamos a energia por partícula (E/N) em cada simulação. Então, geramos um gráfico em função de 1/N, visto que quando N tende ao infinito, 1/N tende a zero, e então onde a curva corta o eixo Y (X=0) temos o valor da energia no limite termodinâmico.

Do nosso gráfico, obtivemos que a energia no limite termodinâmico é $E\approx -1.755$.

Conclusão

Vimos que a implementação do método de Monte Carlo e do algoritmo de Metropolis foram bem sucedidos na tentativa de simular um fluido sob o potencial de Lennard-Jones, de modo que conseguimos extrair propriedades termodinâmicas do sistema. O parâmetro de ordem g(r) foi a ferramenta crucial que possibilitou a identificação das diferentes fases do sistema. Além disso, percebemos que a dinâmica de Metropolis gera amostras altamente correlacionadas, o que nos levou a identificar o tempo de correlação necessário para obter amostras independentes e a utilizar o método de análise de blocos para encontrar o erro real das amostras, provando que o erro ingênuo subestima drasticamente as incertezas das amostras. Por fim, com o FSS conseguimos encontrar o valor da energia do sistema no limite termodinâmico.