Lennard-Jones e propriedades termodinâmicas

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

Potencial de Lennard Jones

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

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

O termo proporcional à ${\frac {1}{r^{6}}}$ é o termo atrativo, enquanto o termo proporcional à ${\frac {1}{r^{12}}}$ é o termo repulsivo, muito mais intenso à curtas distâncias. A constante $\epsilon$ descreve a profundidade do poço deste potencial, enquanto $\sigma$ descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de $\epsilon$ e $\sigma$.

Potencial de Lennard-Jones.

Comportamento do potencial de Lennard-Jones para diferentes valores de $\epsilon$ e $\sigma$.

A força que uma partícula sofre, quando sujeita a este potencial é dada por

${\boldsymbol {F}}({\boldsymbol {r}})=-\nabla V({\boldsymbol {r}}),$ ${\boldsymbol {F}}({\boldsymbol {r}})=4\epsilon \left[12{\frac {\sigma ^{12}}{r^{13}}}-6{\frac {\sigma ^{6}}{r^{7}}}\right]{\hat {\boldsymbol {r}}}.$

Dinâmica Molecular

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

1. Inicializar a posição e velocidade das partículas do sistema.

2. Calcular as forças de interação entre todas as partículas do sistema.

3. Integrar as equações de movimento utilizando algum algoritmo integrador.

4. Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

Inicialização

O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras. Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

a) Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

b) Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

c) Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

Cálculo das forças

Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" $r_{c}$. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

Integração das equações de movimento

Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

${\boldsymbol {v}}\left(t+{\frac {\Delta t}{2}}\right)={\boldsymbol {v}}(t)+{\frac {1}{2}}{\boldsymbol {f}}(t)\Delta t$

${\boldsymbol {r}}\left(t+\Delta t\right)={\boldsymbol {r}}(t)+{\boldsymbol {v}}\left(t+{\frac {\Delta t}{2}}\right)\Delta t$

${\boldsymbol {r}}\left(t+\Delta t\right)\to {\boldsymbol {f}}(t+\Delta t)$

${\boldsymbol {v}}\left(t+\Delta t\right)={\boldsymbol {v}}\left(t+{\frac {\Delta t}{2}}\right){\frac {1}{2}}{\boldsymbol {f}}\left(t+\Delta t\right)\Delta t$

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

Medidas

Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será

$E={\frac {1}{2}}m(v_{x}^{2}+v_{y}^{2})+\sum _{i<J}4\epsilon \left[\left({\frac {\sigma }{r_{ij}}}\right)^{12}-\left({\frac {\sigma }{r_{ij}}}\right)^{6}\right].$

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma

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

onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será

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

onde $N_{f}=2N-2$ é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por

$p=\rho k_{B}T+{\frac {1}{dA}}\left\langle \sum _{i<j}{\boldsymbol {f}}({\boldsymbol {r}}_{ij})\cdot {\boldsymbol {r}}_{ij}\right\rangle ,$

onde $\rho$ é a densidade numérica, $k_{B}$ é a constante de Boltzmann, $T$ é a temperatura, $d$ é a dimensionalidade do sistema (em nosso caso 2) e $A$ é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação $k_{B}$ como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos $k_{B}$ como 1 e todas as unidades das grandezas ficam em termos de $k_{B}$.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

Convenção da imagem mímina e PBC

Convenção da imagem mínima.

Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

Função Distribuição Radial

Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como

$g({\boldsymbol {r}})={\frac {2A}{N^{2}}}\sum _{i<j}\delta ({\boldsymbol {r}}-{\boldsymbol {r}}_{ij}),$

onde A é a área da caixa e N é o número de partículas.

Simulação

Nesta simulação, utilizamos L = 25 e N = 144, onde N é o número de partículas. Neste contexto, obtemos $\rho =0.2304$ e consideramos a massa das partículas como 1. A temperatura desejada que utilizamos foi de T=1.0. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os últimos 100 instantes de tempo, passo à passo e os últimos 100 instantes de tempo, de 5 em 5 passos de tempo. Quando vemos a simulação passo à passo, estamos olhando para uma visão mais microscópica (mecânica estatística), enquanto estamos vendo a simulação de 5 em 5 passo, estamos observando um comportamento mais macroscópico (termodinâmica).

Inicialização do sistema.

Últimos 100 instantes de tempo, passo a passo.

Últimos 100 instantes de tempo, de 5 em 5 passos.

A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.

Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.

Evolução da energia cinética, potencial e mecânica nos primeiros 300 passos da simulação.

Vemos que como o sistema já está equilibrado para uma quantidade de menos de 300 passos, plotamos à direita a série temporal contendo as 300 primeiras medidas para conseguirmos ver o começo da simulação que claramente está fora do equilíbrio.

A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.

Evolução da temperatura durante todos os 10000 passos da simulação.

Evolução da temperatura nos primeiros 300 passos da simulação.

Podemos ver que a temperatura do sistema também já se equilibrou nos primeiros 300 passos e por isso queremos ver o comportamento do sistema fora do equilíbrio no começo da simulação. A partir do equilíbrio, poderíamos utilizar a simulação para calcular propriedades termodinâmicas.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá. Novamente, a pressão equilibra muito cedo e podemos ver que para um tempo menor que 300.

Evolução da pressão durante todos os 10000 passos da simulação.

Evolução da pressão nos primeiros 1000 passos da simulação.

Podemos ainda calcular a função distribuição radial $g({\boldsymbol {r}})$ para esta simulação.

Função distribuição radial.

Códigos

Simulação - Lennard-Jones

Figuras e animações

Bibliografia

D. Frenkel, B. Smit. Understanding Molecular Simulation: From Algorithms to Applications. Academic Press. 2002.