A função radial de distribuição de pares para sistemas Lennard-Jones bidimensionais

The radial pair distribution function for two-dimensional Lennard-Jones systems

ARTIGOS GERAIS

A função radial de distribuição de pares para sistemas Lennard-Jones bidimensionais

The radial pair distribution function for two-dimensional Lennard-Jones systems

Lucas Madeira¹; Silvio A. Vitiello

Instituto de Física "Gleb Wataghin'', Universidade Estadual de Campinas, Campinas, SP, Brasil

RESUMO

Propriedades dos sistemas bidimensionais Lennard-Jones são estudadas com o objetivo de apresentarmos o método de dinâmica molecular. Em particular, damos atenção ao cálculo da função radial de distribuição de pares. A boa concordância de nossos resultados com a teoria e resultados de outra simulação mostra que esse é um método não só de fácil implementação, mas também confiável. Além disso, mostramos que o método de dinâmica molecular pode nos ajudar na interpretação dos resultados e aumentar o nosso entendimento dos mesmos.

Palavras-chave: sistema de muitos corpos, Lennard-Jones, distribuição de pares, Argônio, dinâmica molecular.

ABSTRACT

Properties of two-dimensional Lennard-Jones systems are studied with the aim of introducing the molecular dynamics method. In particular, we give attention to the calculation of the radial distribution function. The good agreement of our results with theory and results of another simulation shows that this is a method not only easy to implement, but also reliable. In addition, we have shown that the molecular dynamics method can help us with the interpretation of results and increase our understanding of them.

Keywords: many-body system, Lennard-Jones, distribution of pairs, Argon, molecular dynamics.

1. Introdução

A caracterização da matéria é importante tanto para a ciência como em aplicações tecnológicas. Um dos aspectos mais básicos que pode nos interessar em sua caracterização está em distinguir os seus três estados: o sólido, o líquido e o gasoso. Uma possível maneira de realizarmos essa tarefa é através da investigação do ordenamento espacial de seus constituintes. Os gases não possuem ordenamento espacial, líquidos possuem uma ordenação de curto alcance, e sólidos possuem uma ordenação de longo alcance.

Simulações computacionais são métodos importantes para compreender as propriedades da matéria. No método de dinâmica molecular, as leis de Newton são aplicadas a átomos ou moléculas, onde fenômenos quânticos podem ser desprezados em boa aproximação, para obter a trajetória de um grande número de partículas em um certo período de tempo, e assim calcular propriedades de interesse do sistema. Atualmente, efeitos quânticos podem ser tratados com a introdução de um tipo especial de termostato, chamado de ruído colorido [1]. Entretanto tal abordagem está fora do escopo desse trabalho. As simulações também permitem reproduzir cenários difíceis, ou até impossíveis, de serem conduzidos em laboratórios (como trabalhar com pressões e temperaturas extremas).

Uma das ferramentas que podemos utilizar para investigar o ordenamento do sistema é a função radial de pares. Neste trabalho iremos ilustrar o cálculo desta função g(r) para o Argônio bidimensional. Os aspectos essenciais deste sistema são razoavelmente bem capturados pelo modelo Lennard-Jones. A descrição de fluídos e sólidos, em algumas situações, deve ser feita em duas dimensões como, por exemplo, quando estão sendo investigados filmes finos.

Simulações bidimensionais são muito convenientes devido aos poucos recursos computacionais que exigem, bem como facilitam visualizar graficamente algumas quantidades físicas de interesse.

Na seção 2 introduzimos o potencial que descreve a interação entre as partículas. A teoria sobre a função radial de distribuição de pares é explicada na seção 3. Apresentamos os conceitos que servem de base para as simulações de dinâmica molecular na seção 4. Os resultados obtidos são discutidos na seção 5, e na seção 6 apresentamos as conclusões.

2. O potencial de interação

Trabalharemos com um sistema de N partículas e vamos supor que a dinâmica dele pode ser tratada classicamente, ou seja, as equações de Newton descrevem o comportamento das partículas. Também vamos supor que a interação entre um par de partículas depende apenas da distância entre elas. O potencial de Lennard-Jones descreve essas interações em razoável aproximação para diversos sistemas reais, como os gases nobres, e será aqui adotado

onde r é a distância entre as partículas, σ é um parâmetro de comprimento e ε é um parâmetro de energia. Na Fig. 1 apresentamos esse potencial parametrizado para o Argônio, onde podemos ver a parte repulsiva (r → 0) e a parte atrativa. A repulsão é devido ao princípio de exclusão de Pauli, o qual produz um efeito de curto alcance, consequência de orbitais eletrônicos que não se sobrepõem. A força de atração é de van der Waals e descreve uma fraca atração de longo alcance.

Para facilitar as simulações, utilizamos como unidade de comprimento σ, de energia e ε de massa m, a massa do átomo. A Tabela 1 ilustra as unidades empregadas. Uma das vantagens de se trabalhar com as chamadas unidades reduzidas é que uma única simulação pode ser utilizada para obter resultados para elementos diferentes, basta utilizarmos os valores numéricos dos parâmetros do elemento desejado.

3. Função radial de distribuição de pares

A função radial de distribuição de pares é uma medida da correlação entre as partículas de um sistema de muitos corpos. Escolhendo-se arbitrariamente uma das partículas como a origem, Fig. 2, temos que o número médio de partículas entre uma distância r e r + dr é de ρg(r)dr, onde g(r) é a função radial de distribuição de pares.

Na prática, ela funciona como uma "densidade local" , ou seja, quando g(r) = 1, a densidade é aquela média do sistema. Valores mais altos da função indicam que a densidade local é maior que a média do sistema, enquanto valores mais baixos mostram que a densidade é menor.

A condição de normalização adotada é dada por

onde o limite superior de integração é metade do lado da célula cúbica de simulação, fato que será melhor discutido na seção 4.2. Essencialmente, a Eq. (2) nos diz que, posicionando uma partícula na origem e contando todas as outras restantes, obtemos N - 1 partículas. Como o valor de N é tipicamente grande em simulações de dinâmica molecular, poderíamos ter feito a aproximação N - 1 ≈ N. Para o potencial de Lennard-Jones, esperamos que g(r) → 0 quando r → 0, pois as partículas não podem penetrar umas nas outras.

A ordenação espacial de um sistema pode ser estudada através da função g(r) e, consequentemente, temos informações sobre a fase em que o sistema se encontra. Consideramos um sistema bidimensional de empacotamento hexagonal. Se r₀ é o parâmetro da rede, teremos seis primeiros vizinhos a uma distância de r₀, seis segundos vizinhos a uma distância de , seis terceiros vizinhos a 2 r₀, e assim sucessivamente (, 3 r₀,...) [4]. A Fig. 3 ilustra geometricamente os três primeiros vizinhos. Neste caso a função radial de distribuição de pares, g(r), seria dada por picos nas distâncias dos primeiros vizinhos, ver Fig. 8. No caso de um gás ideal, em contraste, g(r) seria uma constante, conveniente normalizada para 1.

Para um sólido devemos esperar que a função g(r) apresente picos periódicos correspondentes a essas posições razoavelmente bem definidas, ou seja, os cinco primeiros vizinhos devem aparecer com máximos relativos em distâncias proporcionais a 1 : : 2 : : 3.

4. Dinâmica molecular

4.1. O algoritmo numérico

Como já vimos na seção 2, precisamos computar as trajetórias dos N átomos do sistema para realizar a simulação de dinâmica molecular. Para isso, vamos precisar de um algoritmo com as seguintes características: ele deve ser capaz de resolver numericamente as equações de Newton, de calcular as posições e velocidades das partículas em um instante t e em um instante t + Δt, de permitir o uso de um Δt relativamente grande. Deve ainda, nos intervalos de tempo envolvidos na simulação, calcular as trajetórias clássicas com precisão, conservar energia e momentum, ser reversível no tempo e computacionalmente não dispendioso.

O método empregado no programa para computar a trajetória de cada partícula é o chamado Velocity Verlet. Esse algoritmo é reversível no tempo e é preciso para iterações pequenas de tempo. A posição x_n+1 e a velocidade v_n+1, de uma componente do movimento da partícula, são dadas por [2]

sendo que o índice n denota as quantidades no instante atual e n + 1 naquele seguinte, após uma iteração de tempo Δt. Essas equações são obtidas através da expansão de x_n e x_n+1 em séries de Taylor [5].

As acelerações a_n são obtidas através da Segunda Lei de Newton e, como no nosso sistema de unidades a massa é unitária, a aceleração será numericamente igual à força. Essa, por sua vez, é obtida através do gradiente do potencial

4.2. Condições periódicas de contorno e convenção da mínima imagem

Um sistema físico macroscópico possui um número de partículas da ordem de 10²³, que seguramente não está ao alcance de nenhuma simulação. Através da dinâmica molecular, grandes simulações envolvem um número de partículas da ordem de 10⁵. Nessa situação, a fração de partículas nas bordas do sistema é significativa, logo, não podemos ignorar os efeitos de superfície, os quais são desprezíveis em sistemas reais. Uma maneira de evitar essa dificuldade é adotar as condições periódicas de contorno, juntamente com a convenção da mínima imagem. As condições periódicas de contorno consistem na replicação da caixa de simulação, a qual contém N partículas, em todas as direções, Fig. 4. Assim, cada vez que uma partícula deixa a caixa de simulação na qual ela está, uma outra de uma caixa adjacente entra nessa mesma caixa.

Na convenção da mínima imagem, somente partículas, e suas imagens definidas pela replicação da caixa central, dentro de um raio igual à metade do lado L da caixa de simulação irão interagir com a partícula central, ou seja, o potencial de Lennard-Jonnes (Eq. (1)) só é calculado para distâncias menores do que L/2. A contribuição das demais partículas, com uma distância de separação maior que L/2, é calculada de forma aproximada.

4.3. Implementação do programa

O sistema de N partículas, e densidade ρ = , é posicionado na célula quadrada de lado L e as posições iniciais das N partículas obedecem uma rede triangular, por simples comodidade dos cálculos.

Desejamos controlar a temperatura inicial do sistema, para isso devemos analisar a forma como é calculada a temperatura. Sabemos, do Teorema da Equipartição da Energia, que para cada grau de liberdade do sistema temos

kT de energia associada, onde k é a constante de Boltzmann. Logo, a temperatura como função do tempo é dada por [2]

No entanto, devemos impor a condição de que o centro de massa do sistema esteja fixo. Num experimento real, as paredes do recipiente garantiriam o momento do centro de massa ser nulo. Na simulação impomos que o centro de massa permaneça estático, subtraindo da velocidade de cada partícula aquela do centro de massa. Essa imposição faz com que não tenhamos mais 2N graus de liberdade para o sistema, mas sim 2(N - 1)

onde denota a média temporal.

4.4. Os ensembles NVE e NPT

Os processos descritos até agora formam um sistema com um número de partículas (N) fixo, bem como a área (análoga bidimensional do volume V) e a enegia total E (energia cinética mais a potencial). Por possuir essas quantidades fixas, dizemos que esse sistema é um ensemble NVE. Esse ensemble pode ser entendido como um conjunto de sistemas onde todos os estados consistentes com os valores de NVE estão presentes.

No entanto, experimentos físicos são geralmente conduzidos em condições de temperatura (T) e pressão (P) constantes. Logo, desejamos passar do ensemble NVE para o NPT. Isso é feito através de dois mecanismos de controle, um termostato para a temperatura e um barostato para a pressão.

O termostato consiste em reescalarmos a velocidade de todas as partículas, multiplicando-as por um fator λ. Como já discutimos na seção 4.3, mudando as velocidades do sistema, teremos uma temperatura diferente

Já, o barostato funciona de maneira similar ao termostato, porém neste caso são reescalonadas as distâncias, de tal forma que com mudanças de área da célula de simulação a pressão seja mantida constante. Como a pressão é definida como força por área (no caso bidimensional é força por unidade de comprimento), multiplicamos todas as distâncias envolvidas na simulação por um fator µ, tal que

onde τ é um parâmetro que determina a velocidade com que a pressão desejada (P_t) será alcançada, e P₀ é a pressão calculada antes de reescalar as distâncias.

4.5. Detalhes técnicos da simulação

O programa utilizado nas simulações foi escrito em FORTRAN 95 com base em um código-modelo [2], e implementado em um computador com arquitetura intel i3. Para as simulações com 450 partículas o tempo gasto foi de aproximadamente 1 minuto, enquanto para a simulação de 16200, foi de aproximadamente 13 horas. Os cálculos dependem dos seguintes parâmetros de entrada: número de partículas (N), densidade (ρ), valor da iteração de tempo (Δt), tempo de simulação (unidades reduzidas), temperatura (T) e pressão (P). Além da função g(r), podemos calcular algumas grandezas termodinâmicas, bem como produzir vídeos da posição das partículas. As posições iniciais obedecem uma rede triangular, enquanto as velocidades iniciais são provenientes de um gerador aleatório.

Para os resultados da seção 5 consideramos um ensemble NPT com 450 partículas (N = 450) e o intervalo de tempo entre as interações de 0,01 unidades (Δt = 0,01 ). Os resultados foram obtidos em 10000 passos após o equilíbrio ter sido atingido, o qual teve duração de 1000 passos. É igualmente importante fazermos um tratamento estatístico dos resultados. Estes procedimentos são bem estabelecidos e estão discutidos nas Referências [6, 7].

5. Resultados

A Fig. 5 mostra a função g(r) para duas temperaturas diferentes, T = 0,8 e T = 1,8.

Analisando apenas a curva para T = 0,8, observamos três picos da função g(r), nas distâncias 1,1; 2,1 e 3,1 σ. Sabemos que o empacotamento para o sistema em duas dimensões é hexagonal, logo os picos devem ocorrer na proporção 1 : : 2 : : 3. Os segundos e terceiros vizinhos contribuem para o segundo pico observado na figura, pois a diferença entre e 2 é relativamente pequena. Da mesma maneira, os quartos e quintos vizinhos contribuem para o terceiro pico da figura.

Comparando as duas curvas obtidas através da simulação, T = 0,8 e T = 1,8, notamos que a proporção entre os máximos é preservada. A diferença está no alargamento dos picos e na diminuição no valor absoluto dos máximos. Ambos resultados são explicados pelo maior grau de desordem do sistema, decorrente do aumento da energia.

Para uma densidade correspondente a um líquido, a função g(r) mostra ainda uma certa ordenação do sistema. Analisando a função radial do líquido, Fig. 6, vemos que o terceiro máximo desaparece e o segundo apresenta um alargamento considerável. Já na Fig. 7, com densidade correspondente a gás, podemos notar que o segundo pico desaparece. Esses resultados são consequência da desordenação do sistema, causada pelo espaço livre disponível [8].

Se desejamos observar os cinco primeiros picos separados da função g(r) para o empacotamento hexagonal, 1 : : 2 : : 3, devemos simular uma densidade mais alta do que as anteriores [8]. A Fig. 8 mostra o gráfico de g(r) para ρ = 1. Podemos observar que os máximos estão na proporção 1 : 1,76 : 2 : 2,70 : 3,05, em concordância com os valores esperados pela teoria (≈ 1,73 e ≈ 2,65).

Para fins de validação do nosso método de simulação, podemos comparar os resultados obtidos com outras simulações de dinâmica molecular da literatura. O estudo de Argônio bidimensional, ou seja, filmes, geralmente é feito através da deposição de átomos de Argônio em um substrato. Escolhemos um artigo em que a simulação consiste em observar as três primeiras camadas de Argônio no substrato de grafite [9]. A função g(r), para cada camada, pode ser vista na Fig. 9. Apesar dos parâmetros das simulações serem diferentes, notamos que a forma da função g(r) é igual nas Figs. 8 e 9.

Em um sólido perfeito, esperaríamos que a função de distribuição radial de pares não tendesse a um, como em um líquido ou um gás. Em um sólido ideal, as posições das partículas são altamente correlacionadas, logo devemos ter os máximos da função g(r), dados pela estrutura hexagonal do sólido, em todo o intervalo r considerado.

Em nosso trabalho realizamos uma simulação onde o comportamento real dos átomos é aproximadamente reproduzido. Nele, os átomos vibram em torno das posições da rede que caracteriza o cristal perfeito, comprometendo desta forma a ordem de longo alcance do sistema. Logo, a função g(r) deverá evidenciar esse comportamento. Para verificar esse fato, fizemos uma simulação com os mesmos parâmetros da anterior (T = 0,8 e ρ = 1), porém com 16200 partículas. Isso garantiu que as distâncias r fossem suficientemente grandes para testar a hipótese considerada. Na Fig. 10 podemos ver o resultado. Observamos os máximos correspondentes aos vizinhos no empacotamento hexagonal, porém o valor absoluto deles vai diminuindo com o aumento da distância r, até que g(r) → 1 para distâncias muito grandes.

6. Conclusão

O estudo da função radial de distribuição de pares é uma ótima oportunidade de se ter um primeiro contato com simulações de dinâmica molecular. Com o intuito de calcular g(r), passamos por vários aspectos importantes desse tipo de simulação: o potencial de interação entre as partículas, as unidades utilizadas, o algoritmo integrador, as condições periódicas de contorno e a convenção da mínima imagem, condições iniciais do sistema e propriedades termodinâmicas.

Analisando apenas a função g(r), vemos que ela é uma ótima ferramenta para melhorarmos a nossa compreensão dos estados da matéria. Mudanças na estrutura da função correspondem à mudança no ordenamento das partículas e, consequentemente, das fases do sistema. Os resultados que obtivemos através das simulações de dinâmica molecular mostram boa concordância com a literatura.

7. Agradecimentos

Os autores agradecem a Fundação de Amparo à Pesquisa do Estado de São Paulo e o Conselho Nacional de Desenvolvimento Científico e Tecnológico pelo apoio financeiro.

Recebido em 13/10/2011; Aceito em 30/4/2012; Publicado em 7/12/2012

• [1] M. Ceriotti, G. Bussi and M. Parrinello, Phys. Rev. Lett. 103, 030603 (2009).
• [2] H. Gould, An Introduction to Computer Simulation Methods: Applications to Physical Systems (Addison-Wesley Longman Publishing Co., Boston, 2006), third ed.
• [3] http://www.bss.phy.cam.ac.uk/~amd3/teaching/A_Donald/Amorphous_solids.htm, acessado em 30/11/2012.
  » link
• [4] M. Bishop and C. Bruin, American Journal of Physics 52, 1106-1108 (1984).
• [5] M.P. Allen and D.J. Tildesley, Computer Simulation of Liquids (Clarendon Press, Oxford, 1989).
• [6] R. Canabrava and S.A. Vitiello, Revista Brasileira de Ensino de Física 25, 35 (2003).
• [7] M.A. Dos Reis and S.A. Vitiello, Revista Brasileira de Ensino de Física 28, 45 (2006).
• [8] R.M. Sperandeo-Mineo and G. Tripi, European Journal of Physics 8, 117 (1987).
• [9] Ai Lan Cheng and William A. Steele, Langmuir 5, 60 (1989).

Datas de Publicação

• Publicação nesta coleção
  31 Jan 2013
• Data do Fascículo
  Dez 2012

Histórico

• Recebido
  13 Out 2011
• Aceito
  30 Abr 2012