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

DM de potenciais descontínuos

2016-10-13T00:44:28Z

Amorim: /* Determinação 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>\vec{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.]]

DM de potenciais descontínuos

2016-07-10T01:06:23Z

Amorim: /* Mudança de velocidade em uma colisão elástica */

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

DM de potenciais descontínuos

2016-07-01T12:44:34Z

Amorim: /* Implementação computacional */

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 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-07-01T12:42:44Z

Amorim: /* Implementação computacional */

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

<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-21T18:41:48Z

Amorim: /* Evento dirigido */

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 a força gravitacional aqui computada gera um campo uniforme, ou seja, é radial em relação a cada partícula do sistema, de modo que só há força atuando na direção <math>y</math>, no sentido negativo de <math>y</math>.

Como este campo só atua em <math>y</math>, basta modificar as linha de código acima 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 -2. O par de partículas destacadas em vermelho são as que colidiram.]]

DM de potenciais descontínuos

2016-06-21T18:41:27Z

Amorim: /* Implementação computacional */

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 a força gravitacional aqui computada gera um campo uniforme, ou seja, é radial em relação a cada partícula do sistema, de modo que só há força atuando na direção <math>y</math>, no sentido negativo de <math>y</math>.

Como este campo só atua em <math>y</math>, basta modificar as linha de código acima 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 -2. O par de partículas destacadas em vermelho são as que colidiram.]]

DM de potenciais descontínuos

2016-06-21T14:38:41Z

Amorim: /* Implementação computacional */

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 a força gravitacional aqui computada gera um campo uniforme, ou seja, é radial em relação a cada partícula do sistema, de modo que só há força atuando na direção <math>y</math>, no sentido negativo de <math>y</math>.

Como este campo só atua em <math>y</math>, basta modificar as linha de código acima 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 -2. O par de partículas destacadas em vermelho são as que colidiram.]]

DM de potenciais descontínuos

2016-06-21T14:38:09Z

Amorim: /* Implementação computacional */

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 a força gravitacional aqui computada gera um campo uniforme, ou seja, é radial em relação a cada partícula do sistema, de modo que só há força atuando na direção <math>y</math>, no sentido negativo de <math>y</math>.

Como este campo só atua em <math>y</math>, basta modificar as linha de código acima 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 -2. O par de partículas destacadas em vermelho são as que colidiram.]]

DM de potenciais descontínuos

2016-06-21T14:28:22Z

Amorim: /* Determinação do tempo de colisão */

DM de potenciais descontínuos

2016-06-21T13:51:00Z

Amorim: /* Implementação computacional */

DM de potenciais descontínuos

2016-06-21T13:50:25Z

Amorim: /* Implementação computacional */

DM de potenciais descontínuos

2016-06-21T13:50:13Z

Amorim: /* Implementação computacional */

DM de potenciais descontínuos

2016-06-21T13:49:56Z

Amorim: /* Implementação computacional */

DM de potenciais descontínuos

2016-06-21T13:49:30Z

Amorim: /* Implementação computacional */

DM de potenciais descontínuos

2016-06-21T13:48:18Z

Amorim: /* Otimização básica */

DM de potenciais descontínuos

2016-06-21T13:46:04Z

Amorim: /* Evento dirigido */

Arquivo:Fluxogram.png

2016-06-21T13:45:10Z

Amorim:

Arquivo:Flux.png

2016-06-21T13:43:26Z

Amorim: Amorim enviou uma nova versão de "Arquivo:Flux.png"

Arquivo:Flux.png

2016-06-21T13:41:57Z

Amorim: Amorim enviou uma nova versão de "Arquivo:Flux.png"

DM de potenciais descontínuos

2016-06-21T13:41:29Z

Amorim: /* Evento dirigido */

DM de potenciais descontínuos

2016-06-20T18:32:59Z

Amorim: /* Implementação computacional */

DM de potenciais descontínuos

2016-06-20T18:31:53Z

Amorim: /* Implementação computacional */

DM de potenciais descontínuos

2016-06-20T14:09:02Z

Amorim: /* Implementação computacional */

DM de potenciais descontínuos

2016-06-20T14:08:46Z

Amorim: /* Implementação computacional */

DM de potenciais descontínuos

2016-06-20T14:08:09Z

Amorim: /* Implementação computacional */

DM de potenciais descontínuos

2016-06-20T14:06:43Z

Amorim: /* Implementação computacional */

Arquivo:Pot descontinuo.gif

2016-06-20T14:06:36Z

Amorim:

DM de potenciais descontínuos

2016-06-20T13:55:06Z

Amorim: /* Implementação computacional */

DM de potenciais descontínuos

2016-06-20T13:32:49Z

Amorim: /* Implementação computacional */

DM de potenciais descontínuos

2016-06-20T13:02:38Z

Amorim:

DM de potenciais descontínuos

2016-06-19T01:40:26Z

Amorim: /* Evento dirigido */

DM de potenciais descontínuos

2016-06-19T01:39:20Z

Amorim: /* Otimização básica */

Arquivo:Flux.png

2016-06-19T01:37:33Z

Amorim:

DM de potenciais descontínuos

2016-06-19T00:41:17Z

Amorim: /* Otimização básica */

DM de potenciais descontínuos

2016-06-19T00:25:04Z

Amorim: /* Mudança de velocidade em uma colisão elástica */

DM de potenciais descontínuos

2016-06-19T00:11:42Z

Amorim: /* Mudança de velocidade em uma colisão elástica */

DM de potenciais descontínuos

2016-06-19T00:11:18Z

Amorim: /* Evento dirigido */

DM de potenciais descontínuos

2016-06-19T00:10:04Z

Amorim: /* Mudança de velocidade em uma colisão elástica */

DM de potenciais descontínuos

2016-06-19T00:04:25Z

Amorim: /* Otimização */

DM de potenciais descontínuos

2016-06-18T23:53:06Z

Amorim: /* Determinação do tempo de colisão */

DM de potenciais descontínuos

2016-06-18T23:52:51Z

Amorim: /* Determinação do tempo de colisão */

DM de potenciais descontínuos

2016-06-18T23:52:18Z

Amorim: /* Determinação do tempo de colisão */

DM de potenciais descontínuos

2016-06-18T23:52:07Z

Amorim: /* Determinação do tempo de colisão */

DM de potenciais descontínuos

2016-06-18T23:48:25Z

Amorim: /* Evento dirigido */

DM de potenciais descontínuos

2016-06-18T23:46:42Z

Amorim: /* Determinação do tempo de colisão */

DM de potenciais descontínuos

2016-06-18T23:44:25Z

Amorim: /* Determinação do tempo de colisão */

DM de potenciais descontínuos

2016-06-18T23:43:33Z

Amorim: /* Determinação do tempo de colisão */

DM de potenciais descontínuos

2016-06-18T23:42:44Z

Amorim: /* Determinação do tempo de colisão */

DM de potenciais descontínuos

2016-06-18T23:38:54Z

Amorim: /* Determinação do tempo de colisão */

DM de potenciais descontínuos

2016-06-18T22:56:02Z

Amorim: /* Determinação do tempo de colisão */