DM de potenciais descontínuos

Dinâmica molecular de potenciais descontínuos é uma abordagem computacional utilizada para determinar o movimento de partículas duras que só interagem por forças de contato. Assim, fica evidente a diferença entre o potencial Lennard-Jones pois este se baseia em uma interação de curto alcance, como é mostrado em DM: um primeiro programa. Para entender como as colisões ocorrem, conhecer a forma do potencial a ser estudado é vital. Como estamos considerando corpos rígidos, ou seja, que não sofrem deformação, percebe-se que a força de contato entre as partículas será infinita e o tempo de interação zero, o que torna impossível a descrição do problema a partir de uma integração de movimento simples. O método utilizado, a ser explicitado aqui, que resolve este problema é o evento dirigido.

Evento dirigido

A ideia do método para resolver o problema do força infinita é, ao invés de avançar o sistema em pequenos passos de tempo $dt$, avançar a simulação conforme as colisões forem ocorrendo. Para isso deve-se encontrar o par de partículas $i,j$ que colidirá no menor intervalo de tempo entre todas as partículas, denotaremos tal intervalo por $dt_{min}$, e, então, avançar o sistema. Neste ponto teremos dois objetos colados, portanto aqui deve ser feita a mudança de velocidades de tal forma a respeitar uma colisão elástica.

Fluxograma de um programa simples usando o método de evento dirigido com a otimização aqui explicitada.

Determinação do tempo de colisão

Os objetos a serem usados para o cálculo do tempo de colisão entre um par de partículas $i,j$ serão discos de raio $\sigma _{i}$, $\sigma _{j}$, de distância denotada por $\sigma$. Portanto, segue que a condição de colisão é:

$|{\vec {r_{i}}}(t+dt_{ij})-{\vec {r_{j}}}(t+dt_{ij})|=\sigma$

Com $r_{i}$ sendo o vetor posição da partícula $i$ e $dt_{ij}$ o tempo de colisão entre as partículas $i,j$. Tal condição nos leva a determinação de $dt_{ij}$ a partir da expressão:

$dt_{ij}={\begin{cases}\infty &\quad {\text{se }}d<0\\\infty &\quad {\text{se }}\Delta r.\Delta v>0\\-{\frac {\Delta {\vec {r}}.\Delta {\vec {v}}+{\sqrt {d}}}{\Delta {\vec {v}}.\Delta {\vec {v}}}}&\quad {\text{nos demais casos}}\end{cases}}$

Onde $d\equiv (\Delta {\vec {r}}.\Delta {\vec {v}})^{2}-(\Delta {\vec {v}}.\Delta {\vec {v}})(\Delta {\vec {r}}.\Delta {\vec {r}}-\sigma ^{2})$, $\Delta {\vec {r}}={\vec {r_{i}}}-{\vec {r_{j}}}$ e $\Delta {\vec {v}}={\vec {v_{i}}}-{\vec {v_{j}}}$.
Com isso, consegue-se determinar o valor de $dt_{min}$ encontrando o menor valor de $dt_{ij}$.

Mudança de velocidade em uma colisão elástica

Para fazer a mudança de velocidades temos que considerar o caso de colisão elástica entre as partículas $i,j$, sendo impulso dado por:

${\vec {J}}={\frac {2m_{i}m_{j}(\Delta {\vec {r}}.\Delta {\vec {v}})}{\sigma ^{2}(m_{i}+m_{j})}}\Delta {\vec {r}}$.

Assim, a variação de velocidades pode ser determinada por:

$\Delta {\vec {v_{i}}}=-{\frac {\vec {J}}{m_{i}}}$ e $\Delta {\vec {v_{j}}}={\frac {\vec {J}}{m_{j}}}$.

Otimização básica

Dado uma simulação de $N$ partículas, determinar o valor de $dt_{min}$ é da ordem de $N^{2}$, ou seja, evitar fazer esse processo todo passo de tempo economiza grande parte do tempo computacional. Uma forma simples de fazer isso é determinar e armazenar o menor $dt_{ij}$ para cada partícula $i$ e o índice $j$ antes do loop temporal e a cada passo de tempo determinar o valor de $dt_{min}$ a partir dos $dt_{ij}$ armazenados. Assim, a cada passo de tempo seria necessário apenas atualizar o valor de $dt_{ij}$ das partículas $i,j$ e das que colidiriam com $i$ ou $j$.