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. Importante notar também que na formulação da função exemplo acima foram utilizadas condições de contorno periódicas no sistema. 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.