Física Computacional - Contribuições do usuário [pt-br]

Cálculo do tempo de colisão com aceleração

2016-07-07T16:30:16Z

Esteves: /* Cálculo do tempo de colisão com aceleração */

Tendo em mãos a condição para que ocorra uma colisão: <math>|\vec{r}_{i}(t+dt) - \vec{r}_{j}(t+dt)| = \sigma</math>, e que <math>\vec{r}_{i}(t+dt) = \vec{r}_{i} + \vec{v}_{i}dt + \vec{g}\frac{dt^2}{2}</math> e <math>\vec{r}_{j}(t+dt) = \vec{r}_{j} + \vec{v}_{j}dt + \vec{g}\frac{dt^2}{2}</math>, 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: 
<math>(\vec{r}_{i}(t+dt) - \vec{r}_{j}(t+dt))^2 = \sigma^2</math> 
<math>\vec{r}_{i}^2 + \vec{v}_{i}^2dt^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}^2dt^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</math> 
<math>\vec{r}_{i}^2 + \vec{v}_{i}^2dt^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}^2dt^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</math> 
Reajeitando os termos antes da simplificação final: 
<math>\vec{r}_{i}^2 -2\vec{r}_{i}\vec{r}_{j} + \vec{r}_{j}^2 + \vec{v}_{i}^2dt^2 -2\vec{v}_{i}\vec{v}_{j}dt^2 + \vec{v}_{j}^2dt^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</math> 
Ao simplificar obtemos a mesma equação quadrática para se calcular o <math>dt_{min}</math>, onde ficou definido <math>\Delta \vec{r} = \vec{r}_{i} - \vec{r}_{j}</math> e <math>\Delta \vec{v} = \vec{v}_{i} - \vec{v}_{j}</math> na seção anterior: 
<math>(\vec{r}_{i} - \vec{r}_{j})^2 + (\vec{v}_{i} - \vec{v}_{j})^2dt^2 + 2(\vec{r}_{i} - \vec{r}_{j})(\vec{v}_{i} - \vec{v}_{j})dt -\sigma^2 = 0</math>

Cálculo do tempo de colisão com aceleração

2016-07-07T16:30:05Z

Esteves: Criou página com '== Cálculo do tempo de colisão com aceleração == Tendo em mãos a condição para que ocorra uma colisão: <math>|\vec{r}_{i}(t+dt) - \vec{r}_{j}(t+dt)| = \sigma</math>,...'

== Cálculo do tempo de colisão com aceleração ==

Tendo em mãos a condição para que ocorra uma colisão: <math>|\vec{r}_{i}(t+dt) - \vec{r}_{j}(t+dt)| = \sigma</math>, e que <math>\vec{r}_{i}(t+dt) = \vec{r}_{i} + \vec{v}_{i}dt + \vec{g}\frac{dt^2}{2}</math> e <math>\vec{r}_{j}(t+dt) = \vec{r}_{j} + \vec{v}_{j}dt + \vec{g}\frac{dt^2}{2}</math>, 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: 
<math>(\vec{r}_{i}(t+dt) - \vec{r}_{j}(t+dt))^2 = \sigma^2</math> 
<math>\vec{r}_{i}^2 + \vec{v}_{i}^2dt^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}^2dt^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</math> 
<math>\vec{r}_{i}^2 + \vec{v}_{i}^2dt^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}^2dt^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</math> 
Reajeitando os termos antes da simplificação final: 
<math>\vec{r}_{i}^2 -2\vec{r}_{i}\vec{r}_{j} + \vec{r}_{j}^2 + \vec{v}_{i}^2dt^2 -2\vec{v}_{i}\vec{v}_{j}dt^2 + \vec{v}_{j}^2dt^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</math> 
Ao simplificar obtemos a mesma equação quadrática para se calcular o <math>dt_{min}</math>, onde ficou definido <math>\Delta \vec{r} = \vec{r}_{i} - \vec{r}_{j}</math> e <math>\Delta \vec{v} = \vec{v}_{i} - \vec{v}_{j}</math> na seção anterior: 
<math>(\vec{r}_{i} - \vec{r}_{j})^2 + (\vec{v}_{i} - \vec{v}_{j})^2dt^2 + 2(\vec{r}_{i} - \vec{r}_{j})(\vec{v}_{i} - \vec{v}_{j})dt -\sigma^2 = 0</math>

DM de potenciais descontínuos

2016-07-07T16:29:10Z

Esteves: /* Cálculo do tempo de colisão */

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 <math>dt</math>, avançar a simulação conforme as colisões forem ocorrendo. Para isso deve-se encontrar o par de partículas <math>I, J</math> que colidirá no menor intervalo de tempo entre todas as partículas, denotaremos tal intervalo por <math>dt_{min}</math>, 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.
[[File:fluxogram.png|thumb|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 <math>i, j</math> serão discos de raio <math>\sigma_i</math>, <math>\sigma_j</math>, cuja soma denotamos por <math>\sigma</math>. Portanto, segue que a condição de colisão é:
:<math>|\vec{r_i}(t + dt_{ij}) - \vec{r_j}(t + dt_{ij})| = \sigma</math>
Com <math>r_i</math> sendo o vetor posição da partícula <math>i</math> e <math>dt_{ij}</math> o tempo de colisão entre as partículas <math>i, j</math>. Tal condição nos leva a determinação de <math>dt_{ij}</math> a partir da expressão:
:<math>
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}
</math>
Onde <math> d \equiv (\Delta \vec{r} . \Delta \vec{v})^2 - (\Delta \vec{v} . \Delta \vec{v})(\Delta \vec{r} . \Delta \vec{r} - \sigma^2) </math>, <math> \Delta \vec{r} = \vec{r_i} - \vec{r_j} </math> e <math> \Delta \vec{v} = \vec{v_i} - \vec{v_j} </math>. Com isso, consegue-se determinar o valor de <math> dt_{min} </math> encontrando o menor valor de <math> dt_{ij} </math>.

===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 <math> i, j </math>, sendo impulso dado por:
:<math> \vec{J} = \frac{2m_im_j(\Delta \vec{r} . \Delta \vec{v})}{\sigma^2(m_i + m_j)}\Delta \vec{r} </math>.
Assim, a variação de velocidades pode ser determinada por:
:<math> \Delta \vec{v_i} = -\frac{\vec{J}}{m_i} </math> e <math> \Delta \vec{v_j} = \frac{\vec{J}}{m_j} </math>.

===Otimização básica===
Dado uma simulação de <math> N </math> partículas, determinar o valor de <math> dt_{min} </math> é uma operação de ordem <math> N^2 </math>, 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 <math> dt_{ij} </math> para cada partícula <math> i </math> e o índice <math> j </math> antes do loop temporal e em todo passo de tempo determinar o valor de <math> dt_{min} </math> a partir dos <math> dt_{ij} </math> armazenados. Assim, a cada passo de tempo seria necessário apenas atualizar o valor de <math> dt_{ij} </math> das partículas <math> I, J </math> e das que colidiriam com <math> I </math> ou <math> J </math>.

===Implementação computacional===
Segue a implementação computacional, na linguagem de programação C, a função utilizada para o cálculo de <math>dt_{ij}</math>. Importante notar que na formulação da funções-exemplo foram utilizadas condições de contorno periódicas no sistema.
<pre>
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;
}
}

}
</pre>
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.

<pre>
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;

}
</pre>
Após alterar as velocidades basta avançar o sistema em <math>dt_{min}</math>. Segue um exemplo de função.
<pre>
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;
}

}
</pre>
[[File:pot_descontinuo.gif|400px|thumb|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 <math>y</math>, no sentido negativo de <math>y</math>. 
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 cálculo 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. 
De maneira geral, sem cálculos, o que deve ser feito é abrir cada vetor de deslocamento, adicionando um termo de tempo quadrádico ligado à aceleração gravitacional do campo escolhido, e elevar ao quadrado a condição de que ocorra a colisão. A partir daí tudo que se é necessário fazer é trabalhar a álgebra, que no fim, eliminará todos os termos com a aceleração. Para verificar o cálculo completo veja na página [[Cálculo do tempo de colisão com aceleração]]. 

=== Implementação ===
Como este campo só atua na componente <math>y</math>, 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:
<pre>
yy[i] = yy[i] + vy[i]*delta_t - g*delta_t*delta_t/2;
</pre>
E também adiciona-se uma aceleração em <math>y</math> nesta mesma parte do programa, de modo que haverá uma atualização na velocidade em <math>y</math> de <math>-g*dt_{min}</math>, ficando:
<pre>
vy[i] -= g*delta_t
</pre>
Mais abaixo é possível ver uma animação de como fica o sistema com essas mudanças feitas.
[[File:pot_desc_gravidade.gif|400px|thumb|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.]]

Arquivo:Pot descontinuo.gif

2016-06-27T18:30:30Z

Esteves: Esteves enviou uma nova versão de "Arquivo:Pot descontinuo.gif"

Arquivo:Pot desc gravidade.gif

2016-06-27T18:26:31Z

Esteves: Esteves enviou uma nova versão de "Arquivo:Pot desc gravidade.gif"

Arquivo:Pot desc gravidade.gif

2016-06-27T18:26:11Z

Esteves: Esteves enviou uma nova versão de "Arquivo:Pot desc gravidade.gif"

Arquivo:Pot desc gravidade.gif

2016-06-27T18:24:27Z

Esteves: Esteves enviou uma nova versão de "Arquivo:Pot desc gravidade.gif"

DM de potenciais descontínuos

2016-06-27T18:10:42Z

Esteves: /* Adição do campo gravitacional */

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 <math>dt</math>, avançar a simulação conforme as colisões forem ocorrendo. Para isso deve-se encontrar o par de partículas <math>I, J</math> que colidirá no menor intervalo de tempo entre todas as partículas, denotaremos tal intervalo por <math>dt_{min}</math>, 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.
[[File:fluxogram.png|thumb|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 <math>i, j</math> serão discos de raio <math>\sigma_i</math>, <math>\sigma_j</math>, cuja soma denotamos por <math>\sigma</math>. Portanto, segue que a condição de colisão é:
:<math>|\vec{r_i}(t + dt_{ij}) - \vec{r_j}(t + dt_{ij})| = \sigma</math>
Com <math>r_i</math> sendo o vetor posição da partícula <math>i</math> e <math>dt_{ij}</math> o tempo de colisão entre as partículas <math>i, j</math>. Tal condição nos leva a determinação de <math>dt_{ij}</math> a partir da expressão:
:<math>
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}
</math>
Onde <math> d \equiv (\Delta \vec{r} . \Delta \vec{v})^2 - (\Delta \vec{v} . \Delta \vec{v})(\Delta \vec{r} . \Delta \vec{r} - \sigma^2) </math>, <math> \Delta \vec{r} = \vec{r_i} - \vec{r_j} </math> e <math> \Delta \vec{v} = \vec{v_i} - \vec{v_j} </math>. Com isso, consegue-se determinar o valor de <math> dt_{min} </math> encontrando o menor valor de <math> dt_{ij} </math>.

===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 <math> i, j </math>, sendo impulso dado por:
:<math> \vec{J} = \frac{2m_im_j(\Delta \vec{r} . \Delta \vec{v})}{\sigma^2(m_i + m_j)}\Delta \vec{r} </math>.
Assim, a variação de velocidades pode ser determinada por:
:<math> \Delta \vec{v_i} = -\frac{\vec{J}}{m_i} </math> e <math> \Delta \vec{v_j} = \frac{\vec{J}}{m_j} </math>.

===Otimização básica===
Dado uma simulação de <math> N </math> partículas, determinar o valor de <math> dt_{min} </math> é uma operação de ordem <math> N^2 </math>, 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 <math> dt_{ij} </math> para cada partícula <math> i </math> e o índice <math> j </math> antes do loop temporal e em todo passo de tempo determinar o valor de <math> dt_{min} </math> a partir dos <math> dt_{ij} </math> armazenados. Assim, a cada passo de tempo seria necessário apenas atualizar o valor de <math> dt_{ij} </math> das partículas <math> I, J </math> e das que colidiriam com <math> I </math> ou <math> J </math>.

===Implementação computacional===
Segue a implementação computacional, na linguagem de programação C, a função utilizada para o cálculo de <math>dt_{ij}</math>.
<pre>
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;
}
}

}
</pre>
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.

<pre>
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;

}
</pre>
Após alterar as velocidades basta avançar o sistema em <math>dt_{min}</math>. Segue um exemplo de função.
<pre>
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;
}

}
</pre>
[[File:pot_descontinuo.gif|400px|thumb|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 <math>y</math>, no sentido negativo de <math>y</math>. 
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: <math>|\vec{r}_{i}(t+dt) - \vec{r}_{j}(t+dt)| = \sigma</math>, e que <math>\vec{r}_{i}(t+dt) = \vec{r}_{i} + \vec{v}_{i}dt + \vec{g}\frac{dt^2}{2}</math> e <math>\vec{r}_{j}(t+dt) = \vec{r}_{j} + \vec{v}_{j}dt + \vec{g}\frac{dt^2}{2}</math>, 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: 
<math>(\vec{r}_{i}(t+dt) - \vec{r}_{j}(t+dt))^2 = \sigma^2</math> 
<math>\vec{r}_{i}^2 + \vec{v}_{i}^2dt^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}^2dt^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</math> 
<math>\vec{r}_{i}^2 + \vec{v}_{i}^2dt^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}^2dt^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</math> 
Reajeitando os termos antes da simplificação final: 
<math>\vec{r}_{i}^2 -2\vec{r}_{i}\vec{r}_{j} + \vec{r}_{j}^2 + \vec{v}_{i}^2dt^2 -2\vec{v}_{i}\vec{v}_{j}dt^2 + \vec{v}_{j}^2dt^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</math> 
Ao simplificar obtemos a mesma equação quadrática para se calcular o <math>dt_{min}</math>, onde ficou definido <math>\Delta \vec{r} = \vec{r}_{i} - \vec{r}_{j}</math> e <math>\Delta \vec{v} = \vec{v}_{i} - \vec{v}_{j}</math> na seção anterior: 
<math>(\vec{r}_{i} - \vec{r}_{j})^2 + (\vec{v}_{i} - \vec{v}_{j})^2dt^2 + 2(\vec{r}_{i} - \vec{r}_{j})(\vec{v}_{i} - \vec{v}_{j})dt -\sigma^2 = 0</math> 

=== Implementação ===
Como este campo só atua na componente <math>y</math>, 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:
<pre>
yy[i] = yy[i] + vy[i]*delta_t - g*delta_t*delta_t/2;
</pre>
E também adiciona-se uma aceleração em <math>y</math> nesta mesma parte do programa, de modo que haverá uma atualização na velocidade em <math>y</math> de <math>-g*dt_{min}</math>, ficando:
<pre>
vy[i] -= g*delta_t
</pre>
Mais abaixo é possível ver uma animação de como fica o sistema com essas mudanças feitas.
[[File:pot_desc_gravidade.gif|400px|thumb|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.]]

DM de potenciais descontínuos

2016-06-26T22:40:53Z

Esteves: /* Cálculo do tempo de colisão */

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 <math>dt</math>, avançar a simulação conforme as colisões forem ocorrendo. Para isso deve-se encontrar o par de partículas <math>I, J</math> que colidirá no menor intervalo de tempo entre todas as partículas, denotaremos tal intervalo por <math>dt_{min}</math>, 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.
[[File:fluxogram.png|thumb|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 <math>i, j</math> serão discos de raio <math>\sigma_i</math>, <math>\sigma_j</math>, cuja soma denotamos por <math>\sigma</math>. Portanto, segue que a condição de colisão é:
:<math>|\vec{r_i}(t + dt_{ij}) - \vec{r_j}(t + dt_{ij})| = \sigma</math>
Com <math>r_i</math> sendo o vetor posição da partícula <math>i</math> e <math>dt_{ij}</math> o tempo de colisão entre as partículas <math>i, j</math>. Tal condição nos leva a determinação de <math>dt_{ij}</math> a partir da expressão:
:<math>
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}
</math>
Onde <math> d \equiv (\Delta \vec{r} . \Delta \vec{v})^2 - (\Delta \vec{v} . \Delta \vec{v})(\Delta \vec{r} . \Delta \vec{r} - \sigma^2) </math>, <math> \Delta \vec{r} = \vec{r_i} - \vec{r_j} </math> e <math> \Delta \vec{v} = \vec{v_i} - \vec{v_j} </math>. Com isso, consegue-se determinar o valor de <math> dt_{min} </math> encontrando o menor valor de <math> dt_{ij} </math>.

===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 <math> i, j </math>, sendo impulso dado por:
:<math> \vec{J} = \frac{2m_im_j(\Delta \vec{r} . \Delta \vec{v})}{\sigma^2(m_i + m_j)}\Delta \vec{r} </math>.
Assim, a variação de velocidades pode ser determinada por:
:<math> \Delta \vec{v_i} = -\frac{\vec{J}}{m_i} </math> e <math> \Delta \vec{v_j} = \frac{\vec{J}}{m_j} </math>.

===Otimização básica===
Dado uma simulação de <math> N </math> partículas, determinar o valor de <math> dt_{min} </math> é uma operação de ordem <math> N^2 </math>, 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 <math> dt_{ij} </math> para cada partícula <math> i </math> e o índice <math> j </math> antes do loop temporal e em todo passo de tempo determinar o valor de <math> dt_{min} </math> a partir dos <math> dt_{ij} </math> armazenados. Assim, a cada passo de tempo seria necessário apenas atualizar o valor de <math> dt_{ij} </math> das partículas <math> I, J </math> e das que colidiriam com <math> I </math> ou <math> J </math>.

===Implementação computacional===
Segue a implementação computacional, na linguagem de programação C, a função utilizada para o cálculo de <math>dt_{ij}</math>.
<pre>
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;
}
}

}
</pre>
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.

<pre>
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;

}
</pre>
Após alterar as velocidades basta avançar o sistema em <math>dt_{min}</math>. Segue um exemplo de função.
<pre>
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;
}

}
</pre>
[[File:pot_descontinuo.gif|400px|thumb|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 <math>y</math>, no sentido negativo de <math>y</math>. 
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: <math>|\vec{r}_{i}(t+dt) - \vec{r}_{j}(t+dt)| = \sigma</math>, e que <math>\vec{r}_{i}(t+dt) = \vec{r}_{i} + \vec{v}_{i}dt + \vec{g}\frac{dt^2}{2}</math> e <math>\vec{r}_{j}(t+dt) = \vec{r}_{j} + \vec{v}_{j}dt + \vec{g}\frac{dt^2}{2}</math>, 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: 
<math>(\vec{r}_{i}(t+dt) - \vec{r}_{j}(t+dt))^2 = \sigma^2</math> 
<math>r_{i}^2 + v_{i}^2dt^2 + g^2\frac{dt^4}{4} + 2r_{i}v_{i}dt + 2r_{i}g\frac{dt^2}{2} + 2v_{i}\frac{dt^3}{2} -2(r_{i} + v_{i}dt + g\frac{dt^2}{2})(r_{j} + v_{j}dt + g\frac{dt^2}{2}) + r_{j}^2 + v_{j}^2dt^2 + g^2\frac{dt^4}{4} + 2r_{j}v_{j}dt + 2r_{j}g\frac{dt^2}{2} + 2v_{j}\frac{dt^3}{2}= \sigma^2</math> 
<math>r_{i}^2 + v_{i}^2dt^2 + g^2\frac{dt^4}{2} + 2r_{i}v_{i}dt + 2r_{i}g\frac{dt^2}{2} + 2v_{i}\frac{dt^3}{2} -2r_{i}r_{j} -2r_{i}v_{j}dt -2r_{i}g\frac{dt^2}{2} -2r_{j}v_{i}dt -2v_{i}v_{j}dt^2 -2v_{i}g\frac{dt^3}{2} -2r_{j}g\frac{dt^2}{2} -2v_{j}g\frac{dt^3}{2} -g^2\frac{dt^4}{2} + r_{j}^2 + v_{j}^2dt^2 + 2r_{j}v_{j}dt + 2r_{j}g\frac{dt^2}{2} + 2v_{j}\frac{dt^3}{2}= \sigma^2</math> 
Reajeitando os termos antes da simplificação final: 
<math>r_{i}^2 -2r_{i}r_{j} + r_{j}^2 + v_{i}^2dt^2 -2v_{i}v_{j}dt^2 + v_{j}^2dt^2 + g^2\frac{dt^4}{2} -g^2\frac{dt^4}{2} + 2r_{i}g\frac{dt^2}{2} -2r_{i}g\frac{dt^2}{2} + 2v_{i}\frac{dt^3}{2} -2v_{i}\frac{dt^3}{2} -2r_{j}g\frac{dt^2}{2} + 2r_{j}g\frac{dt^2}{2} -2v_{j}\frac{dt^3}{2} + 2v_{j}\frac{dt^3}{2} + 2r_{i}(v_{i} - v_{j})dt -2r_{j}(v_{i} - v_{j})dt = \sigma^2</math> 
Ao simplificar obtemos a mesma equação quadrática para se calcular o <math>dt_{min}</math>, onde ficou definido <math>\Delta \vec{r} = \vec{r}_{i} - \vec{r}_{j}</math> e <math>\Delta \vec{v} = \vec{v}_{i} - \vec{v}_{j}</math> na seção anterior: 
<math>(\vec{r}_{i} - \vec{r}_{j})^2 + (\vec{v}_{i} - \vec{v}_{j})^2dt^2 + 2(\vec{r}_{i} - \vec{r}_{j})(\vec{v}_{i} - \vec{v}_{j})dt -\sigma^2 = 0</math> 

=== Implementação ===
Como este campo só atua na componente <math>y</math>, 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:
<pre>
yy[i] = yy[i] + vy[i]*delta_t - g*delta_t*delta_t/2;
</pre>
E também adiciona-se uma aceleração em <math>y</math> nesta mesma parte do programa, de modo que haverá uma atualização na velocidade em <math>y</math> de <math>-g*dt_{min}</math>, ficando:
<pre>
vy[i] -= g*delta_t
</pre>
Mais abaixo é possível ver uma animação de como fica o sistema com essas mudanças feitas.
[[File:pot_desc_gravidade.gif|400px|thumb|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.]]

DM de potenciais descontínuos

2016-06-26T22:38:37Z

Esteves: /* Cálculo do tempo de colisão */

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 <math>dt</math>, avançar a simulação conforme as colisões forem ocorrendo. Para isso deve-se encontrar o par de partículas <math>I, J</math> que colidirá no menor intervalo de tempo entre todas as partículas, denotaremos tal intervalo por <math>dt_{min}</math>, 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.
[[File:fluxogram.png|thumb|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 <math>i, j</math> serão discos de raio <math>\sigma_i</math>, <math>\sigma_j</math>, cuja soma denotamos por <math>\sigma</math>. Portanto, segue que a condição de colisão é:
:<math>|\vec{r_i}(t + dt_{ij}) - \vec{r_j}(t + dt_{ij})| = \sigma</math>
Com <math>r_i</math> sendo o vetor posição da partícula <math>i</math> e <math>dt_{ij}</math> o tempo de colisão entre as partículas <math>i, j</math>. Tal condição nos leva a determinação de <math>dt_{ij}</math> a partir da expressão:
:<math>
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}
</math>
Onde <math> d \equiv (\Delta \vec{r} . \Delta \vec{v})^2 - (\Delta \vec{v} . \Delta \vec{v})(\Delta \vec{r} . \Delta \vec{r} - \sigma^2) </math>, <math> \Delta \vec{r} = \vec{r_i} - \vec{r_j} </math> e <math> \Delta \vec{v} = \vec{v_i} - \vec{v_j} </math>. Com isso, consegue-se determinar o valor de <math> dt_{min} </math> encontrando o menor valor de <math> dt_{ij} </math>.

===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 <math> i, j </math>, sendo impulso dado por:
:<math> \vec{J} = \frac{2m_im_j(\Delta \vec{r} . \Delta \vec{v})}{\sigma^2(m_i + m_j)}\Delta \vec{r} </math>.
Assim, a variação de velocidades pode ser determinada por:
:<math> \Delta \vec{v_i} = -\frac{\vec{J}}{m_i} </math> e <math> \Delta \vec{v_j} = \frac{\vec{J}}{m_j} </math>.

===Otimização básica===
Dado uma simulação de <math> N </math> partículas, determinar o valor de <math> dt_{min} </math> é uma operação de ordem <math> N^2 </math>, 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 <math> dt_{ij} </math> para cada partícula <math> i </math> e o índice <math> j </math> antes do loop temporal e em todo passo de tempo determinar o valor de <math> dt_{min} </math> a partir dos <math> dt_{ij} </math> armazenados. Assim, a cada passo de tempo seria necessário apenas atualizar o valor de <math> dt_{ij} </math> das partículas <math> I, J </math> e das que colidiriam com <math> I </math> ou <math> J </math>.

===Implementação computacional===
Segue a implementação computacional, na linguagem de programação C, a função utilizada para o cálculo de <math>dt_{ij}</math>.
<pre>
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;
}
}

}
</pre>
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.

<pre>
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;

}
</pre>
Após alterar as velocidades basta avançar o sistema em <math>dt_{min}</math>. Segue um exemplo de função.
<pre>
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;
}

}
</pre>
[[File:pot_descontinuo.gif|400px|thumb|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 <math>y</math>, no sentido negativo de <math>y</math>. 
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: <math>|\vec{r}_{i}(t+dt) - \vec{r}_{j}(t+dt)| = \sigma</math>, e que <math>\vec{r}_{i}(t+dt) = \vec{r}_{i} + \vec{v}_{i}dt + \vec{g}\frac{dt^2}{2}</math> e <math>\vec{r}_{j}(t+dt) = \vec{r}_{j} + \vec{v}_{j}dt + \vec{g}\frac{dt^2}{2}</math>, 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: 
<math>(\vec{r}_{i}(t+dt) - \vec{r}_{j}(t+dt))^2 = \sigma^2</math> 
<math>r_{i}^2 + v_{i}^2dt^2 + g^2\frac{dt^4}{4} + 2r_{i}v_{i}dt + 2r_{i}g\frac{dt^2}{2} + 2v_{i}\frac{dt^3}{2} -2(r_{i} + v_{i}dt + g\frac{dt^2}{2})(r_{j} + v_{j}dt + g\frac{dt^2}{2}) + r_{j}^2 + v_{j}^2dt^2 + g^2\frac{dt^4}{4} + 2r_{j}v_{j}dt + 2r_{j}g\frac{dt^2}{2} + 2v_{j}\frac{dt^3}{2}= \sigma^2</math>
<math>r_{i}^2 + v_{i}^2dt^2 + g^2\frac{dt^4}{2} + 2r_{i}v_{i}dt + 2r_{i}g\frac{dt^2}{2} + 2v_{i}\frac{dt^3}{2} -2r_{i}r_{j} -2r_{i}v_{j}dt -2r_{i}g\frac{dt^2}{2} -2r_{j}v_{i}dt -2v_{i}v_{j}dt^2 -2v_{i}g\frac{dt^3}{2} -2r_{j}g\frac{dt^2}{2} -2v_{j}g\frac{dt^3}{2} -g^2\frac{dt^4}{2} + r_{j}^2 + v_{j}^2dt^2 + 2r_{j}v_{j}dt + 2r_{j}g\frac{dt^2}{2} + 2v_{j}\frac{dt^3}{2}= \sigma^2</math> 
Reajeitando os termos antes da simplificação final: 
<math>r_{i}^2 -2r_{i}r_{j} + r_{j}^2 + v_{i}^2dt^2 -2v_{i}v_{j}dt^2 + v_{j}^2dt^2 + g^2\frac{dt^4}{2} -g^2\frac{dt^4}{2} + 2r_{i}g\frac{dt^2}{2} -2r_{i}g\frac{dt^2}{2} + 2v_{i}\frac{dt^3}{2} -2v_{i}\frac{dt^3}{2} -2r_{j}g\frac{dt^2}{2} + 2r_{j}g\frac{dt^2}{2} -2v_{j}\frac{dt^3}{2} + 2v_{j}\frac{dt^3}{2} + 2r_{i}(v_{i} - v_{j})dt -2r_{j}(v_{i} - v_{j})dt = \sigma^2</math> 
Ao simplificar obtemos a mesma equação quadrática para se calcular o <math>dt_{min}</math>, onde ficou definido <math>\Delta \vec{r} = \vec{r}_{i} - \vec{r}_{j}</math> e <math>\Delta \vec{v} = \vec{v}_{i} - \vec{v}_{j}</math> na seção anterior: 
<math>(\vec{r}_{i} - \vec{r}_{j})^2 + (\vec{v}_{i} - \vec{v}_{j})^2dt^2 + 2(\vec{r}_{i} - \vec{r}_{j})(\vec{v}_{i} - \vec{v}_{j})dt -\sigma^2 = 0</math> 

=== Implementação ===
Como este campo só atua na componente <math>y</math>, 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:
<pre>
yy[i] = yy[i] + vy[i]*delta_t - g*delta_t*delta_t/2;
</pre>
E também adiciona-se uma aceleração em <math>y</math> nesta mesma parte do programa, de modo que haverá uma atualização na velocidade em <math>y</math> de <math>-g*dt_{min}</math>, ficando:
<pre>
vy[i] -= g*delta_t
</pre>
Mais abaixo é possível ver uma animação de como fica o sistema com essas mudanças feitas.
[[File:pot_desc_gravidade.gif|400px|thumb|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.]]

DM de potenciais descontínuos

2016-06-26T22:22:08Z

Esteves: /* Adição do campo gravitacional */

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 <math>dt</math>, avançar a simulação conforme as colisões forem ocorrendo. Para isso deve-se encontrar o par de partículas <math>I, J</math> que colidirá no menor intervalo de tempo entre todas as partículas, denotaremos tal intervalo por <math>dt_{min}</math>, 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.
[[File:fluxogram.png|thumb|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 <math>i, j</math> serão discos de raio <math>\sigma_i</math>, <math>\sigma_j</math>, cuja soma denotamos por <math>\sigma</math>. Portanto, segue que a condição de colisão é:
:<math>|\vec{r_i}(t + dt_{ij}) - \vec{r_j}(t + dt_{ij})| = \sigma</math>
Com <math>r_i</math> sendo o vetor posição da partícula <math>i</math> e <math>dt_{ij}</math> o tempo de colisão entre as partículas <math>i, j</math>. Tal condição nos leva a determinação de <math>dt_{ij}</math> a partir da expressão:
:<math>
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}
</math>
Onde <math> d \equiv (\Delta \vec{r} . \Delta \vec{v})^2 - (\Delta \vec{v} . \Delta \vec{v})(\Delta \vec{r} . \Delta \vec{r} - \sigma^2) </math>, <math> \Delta \vec{r} = \vec{r_i} - \vec{r_j} </math> e <math> \Delta \vec{v} = \vec{v_i} - \vec{v_j} </math>. Com isso, consegue-se determinar o valor de <math> dt_{min} </math> encontrando o menor valor de <math> dt_{ij} </math>.

===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 <math> i, j </math>, sendo impulso dado por:
:<math> \vec{J} = \frac{2m_im_j(\Delta \vec{r} . \Delta \vec{v})}{\sigma^2(m_i + m_j)}\Delta \vec{r} </math>.
Assim, a variação de velocidades pode ser determinada por:
:<math> \Delta \vec{v_i} = -\frac{\vec{J}}{m_i} </math> e <math> \Delta \vec{v_j} = \frac{\vec{J}}{m_j} </math>.

===Otimização básica===
Dado uma simulação de <math> N </math> partículas, determinar o valor de <math> dt_{min} </math> é uma operação de ordem <math> N^2 </math>, 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 <math> dt_{ij} </math> para cada partícula <math> i </math> e o índice <math> j </math> antes do loop temporal e em todo passo de tempo determinar o valor de <math> dt_{min} </math> a partir dos <math> dt_{ij} </math> armazenados. Assim, a cada passo de tempo seria necessário apenas atualizar o valor de <math> dt_{ij} </math> das partículas <math> I, J </math> e das que colidiriam com <math> I </math> ou <math> J </math>.

===Implementação computacional===
Segue a implementação computacional, na linguagem de programação C, a função utilizada para o cálculo de <math>dt_{ij}</math>.
<pre>
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;
}
}

}
</pre>
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.

<pre>
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;

}
</pre>
Após alterar as velocidades basta avançar o sistema em <math>dt_{min}</math>. Segue um exemplo de função.
<pre>
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;
}

}
</pre>
[[File:pot_descontinuo.gif|400px|thumb|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 <math>y</math>, no sentido negativo de <math>y</math>. 
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: <math>|r_{i}(t+dt) - r_{j}(t+dt)| = \sigma</math>, e que <math>r_{i}(t+dt) = r_{i} + v_{i}dt + g\frac{dt^2}{2}</math> e <math>r_{j}(t+dt) = r_{j} + v_{j}dt + g\frac{dt^2}{2}</math>, 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: 
<math>(r_{i}(t+dt) - r_{j}(t+dt))^2 = \sigma^2</math> 
<math>r_{i}^2 + v_{i}^2dt^2 + g^2\frac{dt^4}{4} + 2r_{i}v_{i}dt + 2r_{i}g\frac{dt^2}{2} + 2v_{i}\frac{dt^3}{2} -2(r_{i} + v_{i}dt + g\frac{dt^2}{2})(r_{j} + v_{j}dt + g\frac{dt^2}{2}) + r_{j}^2 + v_{j}^2dt^2 + g^2\frac{dt^4}{4} + 2r_{j}v_{j}dt + 2r_{j}g\frac{dt^2}{2} + 2v_{j}\frac{dt^3}{2}= \sigma^2</math>
<math>r_{i}^2 + v_{i}^2dt^2 + g^2\frac{dt^4}{2} + 2r_{i}v_{i}dt + 2r_{i}g\frac{dt^2}{2} + 2v_{i}\frac{dt^3}{2} -2r_{i}r_{j} -2r_{i}v_{j}dt -2r_{i}g\frac{dt^2}{2} -2r_{j}v_{i}dt -2v_{i}v_{j}dt^2 -2v_{i}g\frac{dt^3}{2} -2r_{j}g\frac{dt^2}{2} -2v_{j}g\frac{dt^3}{2} -g^2\frac{dt^4}{2} + r_{j}^2 + v_{j}^2dt^2 + 2r_{j}v_{j}dt + 2r_{j}g\frac{dt^2}{2} + 2v_{j}\frac{dt^3}{2}= \sigma^2</math> 
Reajeitando os termos antes da simplificação final: 
<math>r_{i}^2 -2r_{i}r_{j} + r_{j}^2 + v_{i}^2dt^2 -2v_{i}v_{j}dt^2 + v_{j}^2dt^2 + g^2\frac{dt^4}{2} -g^2\frac{dt^4}{2} + 2r_{i}g\frac{dt^2}{2} -2r_{i}g\frac{dt^2}{2} + 2v_{i}\frac{dt^3}{2} -2v_{i}\frac{dt^3}{2} -2r_{j}g\frac{dt^2}{2} + 2r_{j}g\frac{dt^2}{2} -2v_{j}\frac{dt^3}{2} + 2v_{j}\frac{dt^3}{2} + 2r_{i}(v_{i} - v_{j})dt -2r_{j}(v_{i} - v_{j})dt = \sigma^2</math> 
Ao simplificar obtemos a mesma equação quadrática para se calcular o <math>dt_{min}</math>, onde ficou definido <math>\Delta r = r_{i} - r_{j}</math> e <math>\Delta v = v_{i} - v_{j}</math>: 
<math>(r_{i} - r_{j})^2 + (v_{i} - v_{j})^2dt^2 + 2(r_{i} - r_{j})(v_{i} - v_{j})dt -\sigma^2 = 0</math> 
=== Implementação ===
Como este campo só atua na componente <math>y</math>, 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:
<pre>
yy[i] = yy[i] + vy[i]*delta_t - g*delta_t*delta_t/2;
</pre>
E também adiciona-se uma aceleração em <math>y</math> nesta mesma parte do programa, de modo que haverá uma atualização na velocidade em <math>y</math> de <math>-g*dt_{min}</math>, ficando:
<pre>
vy[i] -= g*delta_t
</pre>
Mais abaixo é possível ver uma animação de como fica o sistema com essas mudanças feitas.
[[File:pot_desc_gravidade.gif|400px|thumb|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.]]

DM de potenciais descontínuos

2016-06-26T22:11:12Z

Esteves: /* Adição do campo gravitacional */

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 <math>dt</math>, avançar a simulação conforme as colisões forem ocorrendo. Para isso deve-se encontrar o par de partículas <math>I, J</math> que colidirá no menor intervalo de tempo entre todas as partículas, denotaremos tal intervalo por <math>dt_{min}</math>, 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.
[[File:fluxogram.png|thumb|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 <math>i, j</math> serão discos de raio <math>\sigma_i</math>, <math>\sigma_j</math>, cuja soma denotamos por <math>\sigma</math>. Portanto, segue que a condição de colisão é:
:<math>|\vec{r_i}(t + dt_{ij}) - \vec{r_j}(t + dt_{ij})| = \sigma</math>
Com <math>r_i</math> sendo o vetor posição da partícula <math>i</math> e <math>dt_{ij}</math> o tempo de colisão entre as partículas <math>i, j</math>. Tal condição nos leva a determinação de <math>dt_{ij}</math> a partir da expressão:
:<math>
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}
</math>
Onde <math> d \equiv (\Delta \vec{r} . \Delta \vec{v})^2 - (\Delta \vec{v} . \Delta \vec{v})(\Delta \vec{r} . \Delta \vec{r} - \sigma^2) </math>, <math> \Delta \vec{r} = \vec{r_i} - \vec{r_j} </math> e <math> \Delta \vec{v} = \vec{v_i} - \vec{v_j} </math>. Com isso, consegue-se determinar o valor de <math> dt_{min} </math> encontrando o menor valor de <math> dt_{ij} </math>.

===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 <math> i, j </math>, sendo impulso dado por:
:<math> \vec{J} = \frac{2m_im_j(\Delta \vec{r} . \Delta \vec{v})}{\sigma^2(m_i + m_j)}\Delta \vec{r} </math>.
Assim, a variação de velocidades pode ser determinada por:
:<math> \Delta \vec{v_i} = -\frac{\vec{J}}{m_i} </math> e <math> \Delta \vec{v_j} = \frac{\vec{J}}{m_j} </math>.

===Otimização básica===
Dado uma simulação de <math> N </math> partículas, determinar o valor de <math> dt_{min} </math> é uma operação de ordem <math> N^2 </math>, 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 <math> dt_{ij} </math> para cada partícula <math> i </math> e o índice <math> j </math> antes do loop temporal e em todo passo de tempo determinar o valor de <math> dt_{min} </math> a partir dos <math> dt_{ij} </math> armazenados. Assim, a cada passo de tempo seria necessário apenas atualizar o valor de <math> dt_{ij} </math> das partículas <math> I, J </math> e das que colidiriam com <math> I </math> ou <math> J </math>.

===Implementação computacional===
Segue a implementação computacional, na linguagem de programação C, a função utilizada para o cálculo de <math>dt_{ij}</math>.
<pre>
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;
}
}

}
</pre>
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.

<pre>
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;

}
</pre>
Após alterar as velocidades basta avançar o sistema em <math>dt_{min}</math>. Segue um exemplo de função.
<pre>
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;
}

}
</pre>
[[File:pot_descontinuo.gif|400px|thumb|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 <math>y</math>, no sentido negativo de <math>y</math>. 
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. 
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 de interação entre as colisões. No entanto, como será mostrado a seguir, isso não é necessário e cálculo permanece o mesmo. 
Tendo em mãos a condição para que ocorra uma colisão: <math>|r_{i}(t+dt) - r_{j}(t+dt)| = \sigma</math>, e que <math>r_{i}(t+dt) = r_{i} + v_{i}dt + g\frac{dt^2}{2}</math> e <math>r_{j}(t+dt) = r_{j} + v_{j}dt + g\frac{dt^2}{2}</math>, 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: 
<math>(r_{i}(t+dt) - r_{j}(t+dt))^2 = \sigma^2</math> 
<math>r_{i}^2 + v_{i}^2dt^2 + g^2\frac{dt^4}{4} + 2r_{i}v_{i}dt + 2r_{i}g\frac{dt^2}{2} + 2v_{i}\frac{dt^3}{2} -2(r_{i} + v_{i}dt + g\frac{dt^2}{2})(r_{j} + v_{j}dt + g\frac{dt^2}{2}) + r_{j}^2 + v_{j}^2dt^2 + g^2\frac{dt^4}{4} + 2r_{j}v_{j}dt + 2r_{j}g\frac{dt^2}{2} + 2v_{j}\frac{dt^3}{2}= \sigma^2</math>
<math>r_{i}^2 + v_{i}^2dt^2 + g^2\frac{dt^4}{2} + 2r_{i}v_{i}dt + 2r_{i}g\frac{dt^2}{2} + 2v_{i}\frac{dt^3}{2} -2r_{i}r_{j} -2r_{i}v_{j}dt -2r_{i}g\frac{dt^2}{2} -2r_{j}v_{i}dt -2v_{i}v_{j}dt^2 -2v_{i}g\frac{dt^3}{2} -2r_{j}g\frac{dt^2}{2} -2v_{j}g\frac{dt^3}{2} -g^2\frac{dt^4}{2} + r_{j}^2 + v_{j}^2dt^2 + 2r_{j}v_{j}dt + 2r_{j}g\frac{dt^2}{2} + 2v_{j}\frac{dt^3}{2}= \sigma^2</math> 
Reajeitando os termos antes da simplificação final: 
<math>r_{i}^2 -2r_{i}r_{j} + r_{j}^2 + v_{i}^2dt^2 -2v_{i}v_{j}dt^2 + v_{j}^2dt^2 + g^2\frac{dt^4}{2} -g^2\frac{dt^4}{2} + 2r_{i}g\frac{dt^2}{2} -2r_{i}g\frac{dt^2}{2} + 2v_{i}\frac{dt^3}{2} -2v_{i}\frac{dt^3}{2} -2r_{j}g\frac{dt^2}{2} + 2r_{j}g\frac{dt^2}{2} -2v_{j}\frac{dt^3}{2} + 2v_{j}\frac{dt^3}{2} + 2r_{i}(v_{i} - v_{j})dt -2r_{j}(v_{i} - v_{j})dt = \sigma^2</math> 
Ao simplificar obtemos a mesma equação quadrática para se calcular o <math>dt_{min}</math>: 
<math>(r_{i} - r_{j})^2 + (v_{i} - v_{j})^2dt^2 + 2(r_{i} - r_{j})(v_{i} - v_{j})dt -\sigma^2 = 0</math> 
Como este campo só atua na componente <math>y</math>, basta modificar as linha de código da sessão anterior para adicionar um termo de tempo quadrático à variação da posição assim a linha de código fica:
<pre>
yy[i] = yy[i] + vy[i]*delta_t - g*delta_t*delta_t/2;
</pre>
E também adiciona-se uma aceleração em <math>y</math> nesta mesma parte do programa, de modo que haverá uma atualização na velocidade em <math>y</math> de <math>-g*dt_{min}</math>, ficando:
<pre>
vy[i] -= g*delta_t
</pre>
Mais abaixo é possível ver uma animação de como fica o sistema com essas mudanças feitas.
[[File:pot_desc_gravidade.gif|400px|thumb|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.]]

DM de potenciais descontínuos

2016-06-26T22:10:48Z

Esteves: /* Adição do campo gravitacional */

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 <math>dt</math>, avançar a simulação conforme as colisões forem ocorrendo. Para isso deve-se encontrar o par de partículas <math>I, J</math> que colidirá no menor intervalo de tempo entre todas as partículas, denotaremos tal intervalo por <math>dt_{min}</math>, 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.
[[File:fluxogram.png|thumb|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 <math>i, j</math> serão discos de raio <math>\sigma_i</math>, <math>\sigma_j</math>, cuja soma denotamos por <math>\sigma</math>. Portanto, segue que a condição de colisão é:
:<math>|\vec{r_i}(t + dt_{ij}) - \vec{r_j}(t + dt_{ij})| = \sigma</math>
Com <math>r_i</math> sendo o vetor posição da partícula <math>i</math> e <math>dt_{ij}</math> o tempo de colisão entre as partículas <math>i, j</math>. Tal condição nos leva a determinação de <math>dt_{ij}</math> a partir da expressão:
:<math>
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}
</math>
Onde <math> d \equiv (\Delta \vec{r} . \Delta \vec{v})^2 - (\Delta \vec{v} . \Delta \vec{v})(\Delta \vec{r} . \Delta \vec{r} - \sigma^2) </math>, <math> \Delta \vec{r} = \vec{r_i} - \vec{r_j} </math> e <math> \Delta \vec{v} = \vec{v_i} - \vec{v_j} </math>. Com isso, consegue-se determinar o valor de <math> dt_{min} </math> encontrando o menor valor de <math> dt_{ij} </math>.

===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 <math> i, j </math>, sendo impulso dado por:
:<math> \vec{J} = \frac{2m_im_j(\Delta \vec{r} . \Delta \vec{v})}{\sigma^2(m_i + m_j)}\Delta \vec{r} </math>.
Assim, a variação de velocidades pode ser determinada por:
:<math> \Delta \vec{v_i} = -\frac{\vec{J}}{m_i} </math> e <math> \Delta \vec{v_j} = \frac{\vec{J}}{m_j} </math>.

===Otimização básica===
Dado uma simulação de <math> N </math> partículas, determinar o valor de <math> dt_{min} </math> é uma operação de ordem <math> N^2 </math>, 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 <math> dt_{ij} </math> para cada partícula <math> i </math> e o índice <math> j </math> antes do loop temporal e em todo passo de tempo determinar o valor de <math> dt_{min} </math> a partir dos <math> dt_{ij} </math> armazenados. Assim, a cada passo de tempo seria necessário apenas atualizar o valor de <math> dt_{ij} </math> das partículas <math> I, J </math> e das que colidiriam com <math> I </math> ou <math> J </math>.

===Implementação computacional===
Segue a implementação computacional, na linguagem de programação C, a função utilizada para o cálculo de <math>dt_{ij}</math>.
<pre>
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;
}
}

}
</pre>
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.

<pre>
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;

}
</pre>
Após alterar as velocidades basta avançar o sistema em <math>dt_{min}</math>. Segue um exemplo de função.
<pre>
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;
}

}
</pre>
[[File:pot_descontinuo.gif|400px|thumb|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 <math>y</math>, no sentido negativo de <math>y</math>. 
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. 
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 de interação entre as colisões. No entanto, como será mostrado a seguir, isso não é necessário e cálculo permanece o mesmo. 
Tendo em mãos a condição para que ocorra uma colisão: <math>|r_{i}(t+dt) - r_{j}(t+dt)| = \sigma</math>, e que <math>r_{i}(t+dt) = r_{i} + v_{i}dt + g\frac{dt^2}{2}</math> e <math>r_{j}(t+dt) = r_{j} + v_{j}dt + g\frac{dt^2}{2}</math>, 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: 
<math>(r_{i}(t+dt) - r_{j}(t+dt))^2 = \sigma^2</math> 
<math>r_{i}^2 + v_{i}^2dt^2 + g^2\frac{dt^4}{4} + 2r_{i}v_{i}dt + 2r_{i}g\frac{dt^2}{2} + 2v_{i}\frac{dt^3}{2} -2(r_{i} + v_{i}dt + g\frac{dt^2}{2})(r_{j} + v_{j}dt + g\frac{dt^2}{2}) + r_{j}^2 + v_{j}^2dt^2 + g^2\frac{dt^4}{4} + 2r_{j}v_{j}dt + 2r_{j}g\frac{dt^2}{2} + 2v_{j}\frac{dt^3}{2}= \sigma^2</math>
<math>r_{i}^2 + v_{i}^2dt^2 + g^2\frac{dt^4}{2} + 2r_{i}v_{i}dt + 2r_{i}g\frac{dt^2}{2} + 2v_{i}\frac{dt^3}{2} -2r_{i}r_{j} -2r_{i}v_{j}dt -2r_{i}g\frac{dt^2}{2} -2r_{j}v_{i}dt -2v_{i}v_{j}dt^2 -2v_{i}g\frac{dt^3}{2} -2r_{j}g\frac{dt^2}{2} -2v_{j}g\frac{dt^3}{2} -g^2\frac{dt^4}{2} + r_{j}^2 + v_{j}^2dt^2 + 2r_{j}v_{j}dt + 2r_{j}g\frac{dt^2}{2} + 2v_{j}\frac{dt^3}{2}= \sigma^2</math>
Reajeitando os termos antes da simplificação final: 
<math>r_{i}^2 -2r_{i}r_{j} + r_{j}^2 + v_{i}^2dt^2 -2v_{i}v_{j}dt^2 + v_{j}^2dt^2 + g^2\frac{dt^4}{2} -g^2\frac{dt^4}{2} + 2r_{i}g\frac{dt^2}{2} -2r_{i}g\frac{dt^2}{2} + 2v_{i}\frac{dt^3}{2} -2v_{i}\frac{dt^3}{2} -2r_{j}g\frac{dt^2}{2} + 2r_{j}g\frac{dt^2}{2} -2v_{j}\frac{dt^3}{2} + 2v_{j}\frac{dt^3}{2} + 2r_{i}(v_{i} - v_{j})dt -2r_{j}(v_{i} - v_{j})dt = \sigma^2</math> 
Ao simplificar obtemos a mesma equação quadrática para se calcular o <math>dt_{min}</math>: 
<math>(r_{i} - r_{j})^2 + (v_{i} - v_{j})^2dt^2 + 2(r_{i} - r_{j})(v_{i} - v_{j})dt -\sigma^2 = 0</math> 
Como este campo só atua na componente <math>y</math>, basta modificar as linha de código da sessão anterior para adicionar um termo de tempo quadrático à variação da posição assim a linha de código fica:
<pre>
yy[i] = yy[i] + vy[i]*delta_t - g*delta_t*delta_t/2;
</pre>
E também adiciona-se uma aceleração em <math>y</math> nesta mesma parte do programa, de modo que haverá uma atualização na velocidade em <math>y</math> de <math>-g*dt_{min}</math>, ficando:
<pre>
vy[i] -= g*delta_t
</pre>
Mais abaixo é possível ver uma animação de como fica o sistema com essas mudanças feitas.
[[File:pot_desc_gravidade.gif|400px|thumb|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.]]

DM de potenciais descontínuos

2016-06-20T18:42:42Z

Esteves: /* Adição do campo gravitacional */

DM de potenciais descontínuos

2016-06-20T18:42:24Z

Esteves: /* Adição do campo gravitacional */

DM de potenciais descontínuos

2016-06-20T18:41:35Z

Esteves: /* Adição do campo gravitacional */

DM de potenciais descontínuos

2016-06-20T18:35:52Z

Esteves: /* Adição do campo gravitacional */

DM de potenciais descontínuos

2016-06-20T17:10:47Z

Esteves: /* Adição do campo gravitacional */

Arquivo:Pot desc gravidade.gif

2016-06-20T14:39:34Z

Esteves: Esteves enviou uma nova versão de "Arquivo:Pot desc gravidade.gif"

DM de potenciais descontínuos

2016-06-20T14:36:15Z

Esteves: /* Adição do campo gravitacional */

Arquivo:Pot desc gravidade.gif

2016-06-20T14:35:01Z

Esteves:

DM de potenciais descontínuos

2016-06-20T14:26:34Z

Esteves: /* Adição do campo gravitacional */

DM de potenciais descontínuos

2016-06-20T14:26:15Z

Esteves: /* Adição do campo gravitacional */

DM de potenciais descontínuos

2016-06-20T14:23:11Z

Esteves: /* Adição do campo gravitacional */