Introdução à dinâmica molecular clássica

só um divisor

Esferas duras

Dizemos que a física dos problemas em DM está contida, principalmente, na forma do potencial de interação. Nesta seção vamos analisar o potencial mais simples que podemos pensar \[ u(r) = \left\{\begin{array}{lccr} \infty & r &\leq & \sigma\\ 0 \;\;\;\;\;\; & r & > & \sigma \end{array}\right. \nonumber \] A partir dessa equação vemos que somente existira força quando as partículas tiverem contato direto antes disso são partículas livres e por tant o viajam em linha reta, determinada pela sua velocidade. A colisão descrita por essa equação é uma colisão elástica, não há intercambio de energia entre as partículas e estas não sofrem deformação. A partir das leis de conservação da energia e do momento seremos capazes de predizer quando duas esferas colidem e como está informação poder gerar a dinâmica do sistema.

O algoritmo que governa a colisão de esferas duras é

  1. Procurar quando uma colisão entre duas esferas vai acontecer.
  2. Mover todas as partículas (linha reta) até o momento da colisão.
  3. Implementar as mudanças que deveram acontecer entre as partículas que colidem.
  4. Calcular as propriedades de interesse e retornar a (1).
Agora que sabemos o algoritmo devemos vemos que nosso problema é o calculo do tempo de colisão e como se modificam o estado das partículas que participam da colisão.

Cinemática das colisões

Inicialmente vamos considerara o movimento em uma dimensão, nessas circunstancias teremos:

Seja $m_1$, $\vec{v}_{1i}$ e $\vec{v}_{1f}$ a massa, velocidade antes da colisão e velocidade depois da colisão do corpo $1$, respectivamente; e $m_2$, $\vec{v}_{2i}$ e $\vec{v}_{2f}$ a massa, velocidade antes da colisão e velocidade depois da colisão do corpo $2$, respectivamente, as lei de conservação nos dizem que: \[ \begin{align} m_1v_{1i} + m_2v_{2i} &= m_1v_{1f} + m_2v_{2f} \label{p1}\\ \frac{1}{2}m_1v_{1i}^2 + \frac{1}{2}m_2v_{2i}^2 &= \frac{1}{2}m_1v_{1f}^2 + \frac{1}{2}m_2v_{2f}^2\label{e1} \end{align} \] cancelando os $1/2$ da expressão $\ref{e1}$ e agrupando em ambas equações os termos com massa iguais: \[ \begin{align} m_1\left( v_{1i} - v_{1f} \right) &= m_2\left( v_{2f} - v_{2i} \right) \label{p2}\\ m_1\left( v^2_{1i} - v^2_{1f} \right) &= m_2\left( v^2_{2f} - v^2_{2i} \right) \label{e2} \end{align} \] observe que podemos reescrever a equação $\ref{e2}$ da seguinte forma \[ m_1\left( v_{1i} - v_{1f} \right)\left( v_{1i} + v_{1f} \right) = m_2\left( v_{2f} - v_{2i} \right)\left( v_{2f} + v_{2i} \right) \nonumber \] dividindo esta última equação pela equação $\ref{p2}$ \[ \frac{m_1\bcancel{\left( v_{1i} - v_{1f} \right)}\left( v_{1i} + v_{1f} \right)} {\bcancel{m_1\left( v_{1i} - v_{1f} \right)}}= \frac{\bcancel{m_2\left( v_{2f} - v_{2i} \right)}\left( v_{2f} + v_{2i} \right)} {\bcancel{m_2\left( v_{2f} - v_{2i} \right)}} \nonumber \] de onde \[ \begin{align} v_{1i} + v_{1f} =& v_{2f} + v_{2i} \nonumber\\ v_{1f} =& v_{2f} + v_{2i} - v_{1i} \label{pe1} \end{align} \] colocando em $\ref{p2}$ \[ \begin{align} m_1\left( v_{1i} - v_{2f} - v_{2i} + v_{1i} \right) &= m_2\left( v_{2f} - v_{2i} \right) \nonumber\\ m_1\left( 2v_{1i} - v_{2f} - v_{2i} \right) &= m_2\left( v_{2f} - v_{2i} \right) \nonumber\\ v_{2f} \left( m_2 + m_1 \right) &= v_{2i} \left( m_2 - m_1 \right) + 2m_1v_{1i}\nonumber\\ v_{2f} &= \frac{m_2 - m_1}{m_2 + m_1} v_{2i} + \frac{2m_1}{m_2 + m_1} v_{1i} \label{eq_v2fa} \end{align} \] substituindo em $\ref{pe1}$: \[ \begin{align} v_{1f} =& \frac{m_2 - m_1}{m_2 + m_1} v_{2i} + \frac{2m_1}{m_2 + m_1} v_{1i} + v_{2i} - v_{1i}\nonumber\\ v_{1f} =& \frac{m_1 - m_2}{m_2 + m_1} v_{1i} + \frac{2m_2}{m_2 + m_1} v_{2i} \label{eq_v1fa} \end{align} \] Aqui o calculo foi feito para situação mais genêrica, mas no caso de nossas simulações trataremos com partículas idênticas (mesma massa, $m$, e raio, $\sigma$), assim \[ \begin{align} v_{2f} =& v_{1i}\\ v_{1f} =& v_{2i} \end{align} \] Como esse resultado somente é válido no caso de colisões em uma dimensão devemos estender ele para um número maior de dimensões. A ideia é perceber que quando acontece a colisão só afeta a componente da velocidade que está ao longo da reta que une os centros das partículas no momento da colisão, $\vec{r}_{12} = \vec{r}_2 - \vec{r}_1$. As fim de realizar a analise definimos 3 vetores ortogonais $\left\{\hat{r}_{12}, \hat{j}, \hat{k}\right\}$ e expressamos a velocidade das partículas em termos desses vetores unitários ($1$ e $2$ são os índices das partículas que colidem): \[ \vec{v}_{io} = \left( \vec{v}_{io} \cdot \hat{r}_{12} \right)\hat{r}_{12} + \left( \vec{v}_{io} \cdot \hat{j} \right) \hat{j} + \left( \vec{v}_{io} \cdot \hat{k} \right)\hat{k} \label{eq_vio1} \] onde $o$ indica que se trata da situação antes da colisão, depois da colisão teremos \[ \vec{v}_{if} = \left( \vec{v}_{if} \cdot \hat{r}_{12} \right)\hat{r}_{12} + \left( \vec{v}_{if} \cdot \hat{j} \right) \hat{j} + \left( \vec{v}_{if} \cdot \hat{k} \right)\hat{k} \label{eq_vif1} \] do resultado unidimensional podemos escrever \[ \begin{align*} \vec{v}_{1f} \cdot \hat{r}_{12} &= \vec{v}_{2o} \cdot \hat{r}_{12}\\ \vec{v}_{2f} \cdot \hat{r}_{12} &= \vec{v}_{1o} \cdot \hat{r}_{12}\\ \end{align*} \] enquanto que as componente perpendiculares não se modificam: \[ \begin{align*} \vec{v}_{1f} \cdot \hat{j} &= \vec{v}_{1o} \cdot \hat{j}\\ \vec{v}_{1f} \cdot \hat{k} &= \vec{v}_{1o} \cdot \hat{k}\\ \vec{v}_{2f} \cdot \hat{j} &= \vec{v}_{2o} \cdot \hat{j}\\ \vec{v}_{2f} \cdot \hat{k} &= \vec{v}_{2o} \cdot \hat{k} \end{align*} \] substituindo essas 6 equações em $\ref{eq_vif1}$: \[ \begin{align*} \vec{v}_{1f} =& \left( \vec{v}_{2o} \cdot \hat{r}_{12} \right)\hat{r}_{12} + \left( \vec{v}_{1o} \cdot \hat{j} \right) \hat{j} + \left( \vec{v}_{1o} \cdot \hat{k} \right)\hat{k}\\ =& \left( \vec{v}_{1o} \cdot \hat{r}_{12} \right)\hat{r}_{12} + \left( \vec{v}_{1o} \cdot \hat{j} \right) \hat{j} + \left( \vec{v}_{1o} \cdot \hat{k} \right)\hat{k} + \left[ \left( \vec{v}_{2o} - \vec{v}_{1o} \right) \right] \cdot \hat{r}_{12}\\ \end{align*} \] de onde \[ \vec{v}_{1f} = \vec{v}_{1o} - \left[ \left( \vec{v}_{1o} - \vec{v}_{2o} \right) \right] \cdot \hat{r}_{12} \label{eq_v1fb} \] da mesma forma \[ \vec{v}_{2f} = \vec{v}_{2o} - \left[ \left( \vec{v}_{2o} - \vec{v}_{1o} \right) \right] \cdot \hat{r}_{12} \label{eq_v2fb} \]

Tempo de colisão

A condição para colisão é (eleva ao quadrado) \[ \left[ \vec{r}_1\left( t_c \right) - \left( t_c \right)\right]^2 = \sigma ^2 \nonumber \] onde $t_c$ é o instante em que acontece a colisão. A posição de qualquer partícula está dada por \[ \vec{r}_i\left( t \right) = \vec{r}_i\left( t_o \right) + \left( t - t_o \right) \vec{v}_i\left( t_o \right) \nonumber \] substituindo na equação anterior \[ \begin{align*} \left[ \vec{r}_1\left( t_o \right) + \left( t_c - t_o \right) \vec{v}_1\left( t_o \right) - \vec{r}_2\left( t_o \right) - \left( t_c - t_o \right) \vec{v}_2\left( t_o \right) \right]^2 &= \sigma ^2 \\ \left[ \vec{r}_{12} + \left( t_c - t_o \right)\vec{v}_{12}\right]\,^2 = \sigma ^2 \end{align*} \] onde \[ \begin{align*} \vec{r}_{12} =& \vec{r}_1\left( t_o \right) - \vec{r}_2\left( t_o \right)\\ \vec{v}_{12} =& \vec{v}_1\left( t_o \right) - \vec{v}_2\left( t_o \right)\\ \end{align*} \] Abrindo a equação quadrática obtemos (note a necessidade de passar para escalar a equação, por isso quadrática) \[ \vec{r}_{12}\cdot\vec{r}_{12} - 2\left( t_c - t_o \right)\vec{r}_{12}\cdot\vec{v}_{12} + \left( t_c - t_o \right)^2\vec{v}_{12}\cdot\vec{v}_{12} = 0 \nonumber \] de onde \[ t_c = t_o + \frac{-\vec{r}_{12}\cdot\vec{v}_{12}\pm \sqrt{\left( \vec{r}_{12}\cdot\vec{v}_{12} \right)^2 - v_{12}\left( r_{12}^2 - \sigma^2 \right)}}{r_{12}^2} \label{tc1} \] Contudo está solução deve ser interpretada fisicamente. Observe que se toda colisão deve verificar \[ \vec{r}_{12}\cdot\vec{v}_{12} < 0 \] isto é, devemos verificar que a partícula $1$ e $2$ estão indo ao encontro uma da outra (note que a solução $t_c$ pode dar um valor negativo que diz que no passado houve colisão). Mas está condição só é necessária mas não suficiente, uma situação que não interessa é quando o conjunto de partículas se movem de tal forma que o contato é tal que não se altera a velocidade das partículas depois da colisão, por exemplo, quando a partícula 2 está em repouso e a partícula 1 se move em sua direção mas passa a uma distancia de exatamente $\sigma$, como se mostra no desenho:

Do desenho, essa condição se expressa matematicamente como \[ r_{12}^2 - v_{12}^2 = \sigma^2 \label{ras} \] Aplicando a lei dos cosenos podemos escrever \[ \sigma^2 = r_{12}^2 + v_{12}^2 - 2r_{12}v_{12}\cos \theta \nonumber \] A colisão acontecer efetivamente se o ângulo de colisão $\theta$ for menor que um valor crítico $\theta _L$ ou, equivalentemente \[ 2r_{12}v_{12}\cos \theta \geq 2r_{12}v_{12}\cos \theta _L \nonumber \] que do desenho \[ \cos \theta _L = \frac{v_{12}}{r_{12}} \nonumber \] substituindo na equação previa \[ \begin{align*} r_{12}v_{12}\cos \theta &\geq r_{12}\left( \frac{v_{12}^2}{r_{12}} \right)\\ \left( r_{12}v_{12}\cos \theta\right)^2 &\geq v_{12}^4\\ \left(\vec{r}_{12} \cdot \vec{v}_{12}\right)^2 &\geq v_{12}^4 \end{align*} \] de $\ref{ras}$ \[ v_{12}^2 = r_{12}^2 - \sigma^2 \nonumber \] reescrevemos a equação previa \[ \left(\vec{r}_{12} \cdot \vec{v}_{12}\right)^2 - v_{12}^2\left( r_{12}^2 - \sigma^2 \right)\geq 0 \]

Programa em python sobre dinâmica de esferas duras

No primeiro conjunto de programa temos duas sub rutinas destinadas a definir a configuração inicial das partículas:

      #!/usr/bin/env python
      #-*- coding: utf-8 -*- 

      import numpy as np
      
      def fcc(Ncelas):

        N     = 4 * Ncelas**3
        Runit = np.zeros((4, 3), dtype=np.float)
        R     = np.zeros((N, 3), dtype=np.float)

        #tamanho da cela unitaria
        a = 1.0 / Ncelas

        #posição das partículas na cela unitária
        Runit[0,0] = 0.0
        Runit[0,1] = 0.0
        Runit[0,2] = 0.0
        Runit[1,0] = 0.0
        Runit[1,1] = 0.5 * a
        Runit[1,2] = 0.5 * a
        Runit[2,0] = 0.5 * a
        Runit[2,1] = 0.0
        Runit[2,2] = 0.5 * a
        Runit[3,0] = 0.5 * a
        Runit[3,1] = 0.5 * a
        Runit[3,2] = 0.0

        #coloca todas as partículas
        cont = 0
        for i in xrange(Ncelas):
          for j in xrange(Ncelas):
            for k in xrange(Ncelas):
              for ij in xrange(4):
                if (cont <= N-4):
                  R[cont+ij, 0] = Runit[ij,0] + a*k
                  R[cont+ij, 1] = Runit[ij,1] + a*j
                  R[cont+ij, 2] = Runit[ij,2] + a*i
              cont = cont + 4

        #condicoes periodicas
          for i in xrange(R.shape[0]):
            for j in xrange(3):
              if (R[i,j] >= 0.50):
                R[i,j] = R[i,j] - 1.0
              if (R[i,j] <= -0.50):
                R[i,j] = R[i,j] + 1.0
              
        return R

      #------------------------------------------------------------

      def velIni(N, semente, temperatura):
        sigma = np.sqrt(temperatura)
        V = np.zeros((N, 3), dtype=float)
        conjAleaNum = np.random.RandomState(seed=semente)
        for i in range(3):
          V[:,i] = conjAleaNum.normal(loc=0.0, scale=sigma, size=N)
        vcm = np.sum(V, axis=1)
        for i in range(3):
          V[:,i] = V[:,i] - vcm[i]
        
        return V
      

No segundo conjunto de programas temos 3 sub rutinas destinadas ao calculo das colisões das partículas

      #!/usr/bin/env python
      #-*- coding: utf-8 -*- 

      import numpy as np

      def tempo(sigSq, ri, rj, vi, vj):
        tmin  = 1.0e300
        rij = ri - rj
        for i in xrange(3):
          if (rij[i] >= 0.5):
            rij[i] = rij[i] - 1.0
          if (rij[i] <= -0.5):
            rij[i] = rij[i] + 1.0
        
        vij = vi - vj
        bij = np.sum(rij*vij)  
        
        if (bij < 0.0):
          rij2 = np.sum(rij*rij)
          vij2 = np.sum(vij*vij)
          disc = bij**2 - vij2*(rij2 - sigSq)
          
          if (disc > 0.0):
            tempoCol = (-bij - np.sqrt(disc)) / vij2
            if (tempoCol < tmin):
              tmin = tempoCol
              
        return(tmin)
        
      #------------------------------------------------------------

      def tabela(sigSq, R, V):
        N = R.shape[0]
        tabelaColisao    = np.ones(N, dtype=np.float)
        colideCom        = np.zeros(N, dtype=np.int)
        tabelaColisao[:] = 1.0E300
        
        ri = np.zeros(3, dtype=np.float)
        rj = np.zeros(3, dtype=np.float)
        vi = np.zeros(3, dtype=np.float)
        vj = np.zeros(3, dtype=np.float)
        
        for i in range(N-1):
          ri[:] = R[i,:]
          vi[:] = V[i,:]
          for j in range(i+1,N):
            rj[:] = R[j,:]
            vj[:] = V[j,:]
            
            tempCol = tempo(sigSq, ri, rj, vi, vj)
            if (tempCol < tabelaColisao[i]):
              tabelaColisao[i] = tempCol
              colideCom[i]     = j
            if (tempCol < tabelaColisao[j]):
              tabelaColisao[j] = tempCol
              colideCom[j]     = i
        return(colideCom, tabelaColisao)

      #------------------------------------------------------------

      def atualizaTabela(sigSq, test_i, test_j, colideCom, tabelaColisao, R, V):
        
        ri = np.zeros(3, dtype=np.float)
        rj = np.zeros(3, dtype=np.float)
        vi = np.zeros(3, dtype=np.float)
        vj = np.zeros(3, dtype=np.float)
        
        N = R.shape[0]
        for i in range(N):
          checa = False
          if (i == test_i or colideCom[i] == test_i):
            checa = True
          if (i == test_j or colideCom[i] == test_j):
            checa = True
        
          if (checa):
            ri[:] = R[i,:]
            vi[:] = V[i,:]
            tabelaColisao[i] = 1.0E300
            for j in range(N):
              if (i != j):
                rj[:] = R[j,:]
                vj[:] = V[j,:]

                tempoCol = tempo(sigSq, ri, rj, vi, vj)
                if (tempoCol < tabelaColisao[i]):
                  tabelaColisao[i] = tempoCol
                  colideCom[i]    = j
                if (tempoCol < tabelaColisao[j]):
                  tabelaColisao[j] = tempoCol
                  colideCom[j]    = i
        return(colideCom, tabelaColisao)

      #------------------------------------------------------------

      def proxCol(colideCom, tabelaColisao):
        tempoCol  = tabelaColisao.min()
        colDaPart = tabelaColisao.argmin()

        return (colDaPart, colideCom[colDaPart], tempoCol)      
      

No terceiro conjunto de programas temos uma sub rutinas destinadas ao gravar a configuração das partículas a fim de poder serem visualizadas em programas externos tipo vmd ou jmol

      #!/usr/bin/env python
      #-*- coding: utf-8 -*- 

      import numpy as np

      def xyz(nome, R):
        N = R.shape[0]
        R = 5.0 * R
        arq = open(nome, 'w')
        atomo = 'H'
        print >> arq, N
        print >> arq, 'atoms'
        for i in range(N):
          line = '%s %13.6e %13.6e %13.6e' % (atomo, R[i,0], R[i,1], R[i,2])
          print >> arq, line
        arq.close()      
      

No quarto conjunto de programas temos a sub rutina que modifica a posição e velocidade das partículas

   
      #!/usr/bin/env python
      #-*- coding: utf-8 -*- 

      import numpy as np
      import sys

      def movePos(T, tabelaColisao, R, V):
        tabelaColisao = tabelaColisao - T
        
        teste         = tabelaColisao < 0.0
        if (np.any(teste)):
          sys.exit('Tempo de colisao negativo')    
        
        R  = R + V*T

        #condicoes periodicas
        for i in xrange(R.shape[0]):
          for j in xrange(3):
            if (R[i,j] >= 0.5):
              R[i,j] = R[i,j] - 1.0
            if (R[i,j] <= -0.5):
              R[i,j] = R[i,j] + 1.0
        
        return(tabelaColisao, R)
        
      #------------------------------------------------------------
        
      def modifVel(eCin, sigSq, ri, rj, vi, vj):
        
        dr = ri - rj
        
        for i in xrange(3):
          if (dr[i] >= 0.5):
            dr[i] = dr[i] - 1.0
          if (dr[i] <= -0.5):
            dr[i] = dr[i] + 1.0
            
        dv  = vi - vj
        
        ddv = dr * ( np.sum(dr*dv) / np.sum(dr*dr) )
        
        eCin = eCin - sum(vi*vi) - sum(vj*vj)

        vi = vi - ddv
        vj = vj + ddv
        
        eCin = eCin + sum(vi*vi) + sum(vj*vj)
        
        return(eCin, vi, vj)      
      

Finalmente temos o programa principal

      #!/usr/bin/env python
      #-*- coding: utf-8 -*-
      
      import numpy as np
      import matplotlib.pyplot as plt

      import confInicial as inicial
      import gravaXYZ as grava
      import colisao, dinamica


      Ncelas  = 3
      eta     = 0.3
      rho     = 6.0 * eta / np.pi
      N       = 4 * Ncelas**3
      L       = (float(N) / rho)**(1.0/3.0)
      sigma   = 1.0 / L
      sigSq   = sigma*sigma
      temp    = 0.5
      semente = 39284
      nColEq  = 200
      dt       = 0.01
      tempEq   = 10.0

      R = inicial.fcc(Ncelas)
      N = R.shape[0]
      V  = inicial.velIni(N, semente, temp)

      #Leitura do estado inicial, arq. em binario
      arq = open('estadoInicial.dat', 'rb')
      B = np.fromfile(file=arq, dtype=np.float).reshape(N,6)
      R = B[:,0:3]
      V = B[:,3:N]

      ri = np.zeros(3, dtype=np.float)
      rj = np.zeros(3, dtype=np.float)
      vi = np.zeros(3, dtype=np.float)
      vj = np.zeros(3, dtype=np.float)

      nome = 'ini.xyz'
      grava.xyz(nome, R)

      V2      = V*V
      modV2   = np.sum(V2, axis=1)
      ecinSum = np.sum(modV2)

      colideCom, tabelaColisao = colisao.tabela(sigSq, R, V)

      nome = 'filme.xyz'
      tempoTot = 0.0
      tempoCap = 0.0
      continua = True
      while (continua):
        i, j, temProxCol = colisao.proxCol(colideCom, tabelaColisao)
        #print '%3d %3d %20.15f'%(i, j, temProxCol)
        tempoTot = tempoTot + temProxCol
        
        ri[:] = R[i,:]; rj[:] = R[j,:]
        vi[:] = V[i,:]; vj[:] = V[j,:]
        tabelaColisao, R = dinamica.movePos(temProxCol, tabelaColisao, R, V)
        ecinSum, V[i,:], V[j,:] = dinamica.modifVel(ecinSum, sigSq, ri, rj, vi, vj)
        
        colideCom, tabelaColisao = colisao.atualizaTabela(sigSq, i, j, colideCom,\
                                                          tabelaColisao, R, V)
        
        if (tempoTot >= tempoCap):
          grava.xyz(nome, R)
          tempoCap = tempoCap + dt
          print tempoTot, ecinSum / (3.0*N)
          
        if (tempoTot >= tempEq):
          continua = False   
      

Versão em C utilizando vetor de ponteiros

Versão em Fortran 95

Tarefas

  1. Utilizar a sub rotina que inicializa a velocidade e mostrar que as velocidades respeitam a distribuição de Maxwell-Boltzman, teste se a temperatura obtida a partir das velocidades está em acordo com a temperatura imposta ao sistema
  2. Verificar o diagrama de fase para as esferas duras apresentado no livro Molecular Dynamics Simulation: Elementary Methods de J. M. Haile (página 129 - fig. 3.11). Como referencia pode dar uma olhada também nos artigos originais do tema:
    1. PHASE TRANSITION FOR A HARD SPHERE SYSTEM
    2. STUDIES IN MOLECULAR DYNAMICS .1. GENERAL METHOD
    3. STUDIES IN MOLECULAR DYNAMICS .2. BEHAVIOR OF A SMALL NUMBER OF ELASTIC SPHERES
    4. PHASE TRANSITION IN ELASTIC DISKS
  3. Calcule as funções de distribuição radial utilizando como fator de empacotamento 0.45, 0.55 e 0.65.