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}$, cuja soma denotamos 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}$ é uma operação de ordem $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 em todo 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_ij(double *xx, double *yy, double *vx, double *vy, double *dt_ij, int *colide, int i, int j){

//dt_ij é um vetor que armazena os menores valor de dt_ij para cada i
//colide é um vetor que armazena o valor de j correspondente a partícula que i colide no menor tempo dt_ij

  double delta_t_ij = 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_ij = (-1)*(drdv + sqrt(d))/(dvdv);
      if(delta_t_ij < dt_ij[i]){
	dt_ij[i] = delta_t_ij;
  	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;

}

Após alterar as velocidades basta avançar o sistema em $dt_{min}$. Segue um exemplo de função.

void pos(double *xx, double *yy, double *vx, double *vy, double dt_min){

  int i;

  for(i = 0; i < (NP); i++){
    xx[i] = xx[i] + vx[i]*dt_min;
    yy[i] = yy[i] + vy[i]*dt_min;
  }

}

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 colisõ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 campo gravitacional aqui computado gera uma força central, 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$.
Também é importante lembrar que aqui é desconsiderado a interação gravitacional entre quaisquer duas partículas, uma vez que assume-se que massas são pequenas e desprezíveis em frente a massa que forma o campo em que estão imersas.

Cálculo do tempo de colisão

Uma preocupação que alguém pode ter com a adição de um campo gravitacional na simulação seria de que seria necessário mudar o calculo do tempo mínimo entre as colisões. No entanto, como será mostrado a seguir, isso não é necessário e o cálculo permanece o mesmo.
Tendo em mãos a condição para que ocorra uma colisão: $|{\vec {r}}_{i}(t+dt)-{\vec {r}}_{j}(t+dt)|=\sigma$, e que ${\vec {r}}_{i}(t+dt)={\vec {r}}_{i}+{\vec {v}}_{i}dt+{\vec {g}}{\frac {dt^{2}}{2}}$ e ${\vec {r}}_{j}(t+dt)={\vec {r}}_{j}+{\vec {v}}_{j}dt+{\vec {g}}{\frac {dt^{2}}{2}}$, podemos elevar os dois lados da equação da condição ao quadrado, retirar o módulo uma vez que o resultado será sempre positivo, obtendo:
$({\vec {r}}_{i}(t+dt)-{\vec {r}}_{j}(t+dt))^{2}=\sigma ^{2}$

${\vec {r}}_{i}^{2}+{\vec {v}}_{i}^{2}dt^{2}+{\vec {g}}^{2}{\frac {dt^{4}}{4}}+2{\vec {r}}_{i}{\vec {v}}_{i}dt+2{\vec {r}}_{i}{\vec {g}}{\frac {dt^{2}}{2}}+2{\vec {v}}_{i}{\frac {dt^{3}}{2}}-2({\vec {r}}_{i}+{\vec {v}}_{i}dt+{\vec {g}}{\frac {dt^{2}}{2}})({\vec {r}}_{j}+{\vec {v}}_{j}dt+{\vec {g}}{\frac {dt^{2}}{2}})+{\vec {r}}_{j}^{2}+{\vec {v}}_{j}^{2}dt^{2}+{\vec {g}}^{2}{\frac {dt^{4}}{4}}+2{\vec {r}}_{j}{\vec {v}}_{j}dt+2{\vec {r}}_{j}{\vec {g}}{\frac {dt^{2}}{2}}+2{\vec {v}}_{j}{\frac {dt^{3}}{2}}=\sigma ^{2}$

${\vec {r}}_{i}^{2}+{\vec {v}}_{i}^{2}dt^{2}+{\vec {g}}^{2}{\frac {dt^{4}}{2}}+2{\vec {r}}_{i}{\vec {v}}_{i}dt+2{\vec {r}}_{i}{\vec {g}}{\frac {dt^{2}}{2}}+2{\vec {v}}_{i}{\frac {dt^{3}}{2}}-2{\vec {r}}_{i}{\vec {r}}_{j}-2{\vec {r}}_{i}{\vec {v}}_{j}dt-2{\vec {r}}_{i}{\vec {g}}{\frac {dt^{2}}{2}}-2{\vec {r}}_{j}{\vec {v}}_{i}dt-2{\vec {v}}_{i}{\vec {v}}_{j}dt^{2}-2{\vec {v}}_{i}{\vec {g}}{\frac {dt^{3}}{2}}-2{\vec {r}}_{j}{\vec {g}}{\frac {dt^{2}}{2}}-2{\vec {v}}_{j}{\vec {g}}{\frac {dt^{3}}{2}}-{\vec {g}}^{2}{\frac {dt^{4}}{2}}+{\vec {r}}_{j}^{2}+{\vec {v}}_{j}^{2}dt^{2}+2{\vec {r}}_{j}{\vec {v}}_{j}dt+2{\vec {r}}_{j}{\vec {g}}{\frac {dt^{2}}{2}}+2{\vec {v}}_{j}{\vec {g}}{\frac {dt^{3}}{2}}=\sigma ^{2}$

Reajeitando os termos antes da simplificação final:
${\vec {r}}_{i}^{2}-2{\vec {r}}_{i}{\vec {r}}_{j}+{\vec {r}}_{j}^{2}+{\vec {v}}_{i}^{2}dt^{2}-2{\vec {v}}_{i}{\vec {v}}_{j}dt^{2}+{\vec {v}}_{j}^{2}dt^{2}+{\vec {g}}^{2}{\frac {dt^{4}}{2}}-{\vec {g}}^{2}{\frac {dt^{4}}{2}}+2{\vec {r}}_{i}{\vec {g}}{\frac {dt^{2}}{2}}-2{\vec {r}}_{i}{\vec {g}}{\frac {dt^{2}}{2}}+2{\vec {v}}_{i}{\frac {dt^{3}}{2}}-2{\vec {v}}_{i}{\vec {g}}{\frac {dt^{3}}{2}}-2{\vec {r}}_{j}{\vec {g}}{\frac {dt^{2}}{2}}+2{\vec {r}}_{j}{\vec {g}}{\frac {dt^{2}}{2}}-2{\vec {v}}_{j}{\vec {g}}{\frac {dt^{3}}{2}}+2{\vec {v}}_{j}{\vec {g}}{\frac {dt^{3}}{2}}+2{\vec {r}}_{i}({\vec {v}}_{i}-{\vec {v}}_{j})dt-2{\vec {r}}_{j}({\vec {v}}_{i}-{\vec {v}}_{j})dt=\sigma ^{2}$

Ao simplificar obtemos a mesma equação quadrática para se calcular o $dt_{min}$, onde ficou definido $\Delta {\vec {r}}={\vec {r}}_{i}-{\vec {r}}_{j}$ e $\Delta {\vec {v}}={\vec {v}}_{i}-{\vec {v}}_{j}$ na seção anterior:
$({\vec {r}}_{i}-{\vec {r}}_{j})^{2}+({\vec {v}}_{i}-{\vec {v}}_{j})^{2}dt^{2}+2({\vec {r}}_{i}-{\vec {r}}_{j})({\vec {v}}_{i}-{\vec {v}}_{j})dt-\sigma ^{2}=0$

Implementação

Como este campo só atua na componente $y$, basta modificar as linha de código da seção anterior para adicionar um termo de tempo quadrático à variação da posição assim a linha de código fica:

yy[i] = yy[i] + vy[i]*delta_t - g*delta_t*delta_t/2;

E também adiciona-se uma aceleração em $y$ nesta mesma parte do programa, de modo que haverá uma atualização na velocidade em $y$ de $-g*dt_{min}$, ficando:

vy[i] -= g*delta_t

Mais abaixo é possível ver uma animação de como fica o sistema com essas mudanças feitas.

Plot exemplo de resultado de simulação com um campo gravitacional uniforme de g=2. O par de partículas destacadas em vermelho são as que colidiram.