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$.

Implementação computacional

Segue a implementação computacional, na linguagem de programação C, a função utilizada para o cálculo de $dt_{ij}$.

void calc_dt(double *xx, double *yy, double *vx, double *vy, double *deltat, int *colide, int i, int j){

  double delta_t = INF;
  double dx, dy, dvx, dvy;
  double drdr, dvdv, drdv;
  double d;
  double X, Y;
  double sigma = 2*raio;
	
  X = xx[i] - xx[j];
  dx = fmod(X, Lx) - rint(fmod(X, Lx)/Lx)*Lx;
  Y = yy[i] - yy[j];
  dy = fmod(Y, Ly) - rint(fmod(Y, Ly)/Ly)*Ly;

  dvx = vx[i] - vx[j];
  dvy = vy[i] - vy[j];
  
  drdv = dx*dvx + dy*dvy;
  drdr = dx*dx + dy*dy;
  dvdv = dvx*dvx + dvy*dvy;
  
  d = pow(drdv, 2) - dvdv*(drdr - pow(sigma, 2));

  if(d > 0 && drdv < 0){
    delta_t = (-1)*(drdv + sqrt(d))/(dvdv);
      if(delta_t < deltat[i]){
	deltat[i] = delta_t;
	colide[i] = j;
      }
  }

}

Onde foram definidos INF como um número computacionalmente grande e raio como o raio das partículas, em particular consideramos os raios de todas como iguais. Para realizar a troca de velocidades do par de partículas que colidem de forma a ser coerente com uma colisão elástica, confeccionou-se a função abaixo.

void switch_veloc(double *xx, double *yy, double *vx, double *vy, int i, int j, double *Jx, double *Jy){

  double dx, dy, dvx, dvy;
  double drdr, dvdv, drdv;
  double deltavx, deltavy;
  double sigma = 2*raio, sigma2 = pow(sigma,2), d;
  double X, Y;
	
  X = xx[i] - xx[j];
  dx = fmod(X, Lx) - rint(fmod(X, Lx)/Lx)*Lx;
  Y = yy[i] - yy[j];
  dy = fmod(Y, Ly) - rint(fmod(Y, Ly)/Ly)*Ly;

  dvx = vx[i] - vx[j];
  dvy = vy[i] - vy[j];
  
  drdv = dx*dvx + dy*dvy;
  drdr = dx*dx + dy*dy;
  dvdv = dvx*dvx + dvy*dvy;

  Jx[i] = drdv*dx/sigma2;
  Jy[i] = drdv*dy/sigma2;

  deltavx = Jx[i];
  deltavy = Jy[i];

  //TROCA VELOCIDADES DA PARTICULA i
  vx[i] -= deltavx;
  vy[i] -= deltavy;

  //TROCA VELOCIDADES DA PARTICULA j
  vx[j] += deltavx;
  vy[j] += deltavy;

}

Plot exemplo de resultado de simulação. O par de partículas destacadas em vermelho são as colisoras.

Seguindo o fluxograma apresentado e utilizando as funções disponíveis chega-se, por exemplo, na animação ao lado, que mostra o decorrer da simulação conforme ocorrem as simulações.

Adição do campo gravitacional

A adição de um campo gravitacional uniforme nesta simulação é bastante simples, basta definir uma nova variável, a qual aqui será chamada de g, e atribuir um valor de sua escolha a ela. Vale lembrar que a força gravitacional aqui computada gera um campo uniforme, ou seja, é radial em relação a cada partícula do sistema, de modo que só há força atuando na direção y, no sentido negativo de y.

Como este campo só atua em y, basta modificar as linha de código acima para adicionar um termo de tempo quadrático à variação da posição (adicionar treco aqui dps da linha cima) e adicionar uma aceleração em y nesta mesma parte do programa, de modo que haverá uma atualização na velocidade em y de -g*dt. Mais acima também é possível ver uma animação de como fica o sistema com essas especificações comentadas com g = -2.