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

só um divisor

Introdução à dinâmica molecular II

Na primeira aula desta seção vimos algo das ferramentas utilizadas na dinâmica molecular necessárias para podermos realizar as simulações de esferas duras que foram vistas. Vimos também como simular sistemas onde a interação entre as partículas acontece somente quando colidem. Nesta seção analisaremos situações onde o potencial de interação tem um alcance finito maior que raio da partícula (até 3 o mais vezes o raio destas), em especial estudaremos o potencial de Lennard-Jones em algum detalhe, mas os métodos que serão desenvolvidos são aplicáveis a qualquer tipo de potencial de pares.

O potencial de Lennard-Jones

O potencial de interação inter-atômico mais famoso que se utiliza para simular sistemas atômicos, moleculares e poliméricos, entre outros, é o chamado potencia de Lennard-Jones. Esse potencial foi proposto por Sir John Edward Lennard-Jones no seu artigo intitulado On the determination of molecular fields-II. From the equation of state of a Gas. Nesse trabalho o LJ propõe que o potencial de interação inter partículas para gases deve ter a forma genérica $\lambda_nr^{-n}-\lambda_mr^{-m}$ com a exigência que $\{n,m\} \nless 4$. Num trabalho posterior intitulado The equation of state of gases and critical Phenomena, ele chega ao valor atual dos expoentes mais usados, $n=12$ e $m=6$. Na atualidade a expressão completa para o chamado potencial de Lennard-Jones 6-12 é

\[ U(r) = 4\epsilon\left[ \left( \dfrac{\sigma}{r} \right)^{12} - \left( \dfrac{\sigma}{r} \right)^6 \right] \] ou \[ U(r) = \epsilon\left[ \left( \dfrac{r_o}{r} \right)^{12} - 2\left( \dfrac{r_o}{r} \right)^6 \right] \]

Ainda que ambas expressões sejam iguais se admitimos que

\[ r_0 = 2^{1/6}\sigma \]

elas tem significados físicos diferentes, uma relaciona o potencial com o tamanho da partícula, $\sigma$, enquanto que a outra estabelece o potencial a partir do estado de equilíbrio de um dímero a temperatura zero.

Um dos exponentes do potencial de Lennard-Jones tem significado físico, pode se mostrar que a atração em sistemas que interagem através de ligações de tipo dipolo-dipolo devido a dipolo flutuante (interação de van der Waals) pode ser bem descrita por potenciais tais que $U \sim 1/r^6$. Sistemas atômicos tipo gases nobres (os quais tem todas suas camadas eletrônicas preenchidas) são muito bem modelados por está atração, igualmente algum sistemas moleculares como $CH_4$, $O_2$, $H_2$, $C_2H_4$, etc (ver JCP 104, pg. 8627, 1996). O outro expoente leva em consideração a repulsão que se obtêm quando a nuvem eletrônicas se aproximam de mais violando a regra exclusão de Pauli. A escolha da repulsão como $\sim 1/r{12}$ é arbitraria, mas por simplificar os cálculos numéricos.

Como mencionado previamente, gases nobres são bem representados pelos potencial de LJ, os parâmetros desse potencial, para o caso dos diversos gases nobre são

Elemento    $\sigma$ $[eV]$       $\epsilon$ $[10^{-10}]$   
   $He$       0.00088       2.56   
   $Ne$       0.0031       2.74   
   $Ar$       0.0104       3.40   
   $Kr$       0.014       3.65   
   $Xe$       0.0200       3.98   

A interações entre os materiais tem outras origens em geral, diferentes da de van der Waals. Nesses casos o potencial de LJ não deveria ser utilizado mas como ele é "fisicamente coerente" quando desejamos utilizar obter informação genérica dos sistemas tem sido utilizado. Uma vantagem deste potencial sobre outros potenciais mais realista é sua simplicidade, a interação de miles de átomos podem ser levada em consideração se ele é utilizado e em algumas situações o número de partículas é importante.

Outros potenciais

Um outro potencial que também é bastante utilizado é o chamado potencial de Morse

\[ U(r) = \epsilon \left[ \left( 1 - e^{-\left(r - r_o \right) / \sigma} \, \right)^2 - 1\right] \]

A forma do potencial é similar à do potencial de LJ, portanto também é um potencial razoável. Quando introduzido por P. Morse, ele foi utilizado para analisar a contribuição anharmonica da interação entre moléculas levando em consideração comportamentos quanticos (substituiu o potencial harmônico por esse potencial).

Outros potenciais são bastante complexos se comparados como os dois prévios. Para modelarmos cristais ionicos podemos utilizar o potencial de Born-Mayer:

\[ U = \frac{1}{2}\sum_{i,j\neq i}^N \frac{Z_iZ_je}{r_{ij}} + \sum_{i,j\neq i}^{viz}A_{ij}\left( 1 + \frac{Z_i}{z_i} + \frac{Z_j}{z_j} \right) \exp\left[ \left( \alpha_i + \alpha_j - r_{ij}\right) / \rho_{ij} \right] \]

nessa expressão $Z$ é o número atômico, $z$ é a quantidade de elétrons na camada externa, $\sigma$ e $\rho$ são parâmetros a serem ajustados dependente do tipo de íons que participam da interação.

Existem outras situações onde a ligação depende da existência de outras ligações, por exemplo, ao formar um cristal colocamos primeiro 2 átomos, quando colocamos o terceiro átomo a interação dos 2 primeiros átomos se debilita e a energia global da interação das 3 partículas é menor do que quando tínhamos só 2 partículas, a este comportamento se lhe conhece como ordem de ligação, dessa forma uma sistema pode preferir um arranjo bcc em lugar de um fcc. O potencial de Tersoff simula este comportamento: \[ U = \frac{1}{2}\sum_{i,j\neq i}^N\left[ \phi_R(r_{ij})+b(z_{ij})\phi_A(r_{ij}) \right] \] onde \[ \phi_R(r) = Ae^{-\lambda_1 r}\;\;\;\;\;\;\;\;\;\phi_A(r) = -Be^{-\lambda_2 r} \] \[ b(z) = \left( 1 + \delta^n z^n \right)^{-1/2n}\;\;\;\;\;\;\;\;\; z_{ij} = \sum_{k,k\neq i,j} f(r_{ij}g\left(\theta_{ijk} \right) e^{\lambda_3^3\left( r_{ik} - r_{jk} \right)} \] onde \[ f(r) = \left\{\begin{array}{l l} 1 & r < R-D\\ \frac{1}{2} - \frac{1}{2}\sin\left( \frac{\pi}{2}\frac{r-R}{D} \right) & R-D < r < R+D\\ 0 & r > R+D \end{array}\right. \] e \[ g\left(\theta \right) = 1 + \frac{c^2}{d^2}- \frac{c^2}{d^2 + \left( h - \cos \theta \right)} \] sendo que $\theta_{ijk}$ é o ângulo formado entre as ligações $i-j$ e $i-k$. Note que este potencial possui $12$ parâmetros $A,\,B,\,\lambda_1,\,\lambda_2,\,\lambda_3, \, \delta, \, n, \, R, \, D, \, c, \, d, \, h$

Unidades reduzidas

Potenciais simples como o de LJ e Morse permitem que sejam definidas unidades de medida convenientes a serem utilizadas durante a simulação. A vantagem desta abordagem é que como este tipo de potencial é genérico, utilizando unidades reduzidas fica ainda mais claro esse carater pois os resultados se referem a todo uma família de sistema e não a um em particular. Outra vantagem é que muitos dos parâmetros da simulação assumem valores unitários o que tem um impacto positivo nos cálculos, fora o fato de expressar as equações de forma mais compacta.

Para o caso do potencial de Lennard-Jones as unidades reduzidas mais usuais são: escolher $m\equiv 1$, $\sigma \equiv 1 $ e $\epsilon \equiv 1$ dentro do código o que implica que:

grandeza unidade
comprimento       $r^* = \dfrac{r}{\sigma}$
energia    $E^* = \dfrac{E}{\epsilon}$
tempo    $t^* = \left( \dfrac{\epsilon}{m\sigma^2} \right)^{1/2}t$
temperatura    $T^* = \dfrac{k_BT}{\epsilon}$
força    $f^* = f\dfrac{\sigma}{\epsilon}$
pressão $P^* = P \dfrac{\sigma^3}{\epsilon}$
densidade    $\rho^* = \rho \dfrac{\sigma^3}{m}$

Força para o caso de partículas que interagem via LJ

Da mecânica clássica sabemos que a força que age sobre a partícula $i$ se obtém a partir de \[ \vec{F}_i = -\vec{\nabla}_{r_i}U\left( \vec{r}_1, \vec{r}_2, \vec{r}_3, \ldots, \vec{r}_n \right) \nonumber \] sendo $U$ o potencial de interação dado por \[ U\left( \vec{r}_1, \vec{r}_2, \vec{r}_3, \ldots, \vec{r}_n \right) = \sum_i \sum_{j > i} U\left( r_{ij} \right) \] onde \[ r_{ij} = \left| \vec{r}_i - \vec{r}_j \right| = \sqrt{\left(x_i - x_j\right) + \left(y_i - y_j\right)+\left(z_i - z_j\right)} \nonumber \] de forma que \[ \begin{aligned} \vec{F}_i =& -\vec{\nabla}_{r_i} U\left( r_{ij} \right)\\ =& -\sum_{j \neq i} \left( \dfrac{\vec{r}_i - \vec{r}_j}{\left| \vec{r}_i - \vec{r}_j \right|} \right) \dfrac{dU\left( r_{ij} \right)}{dr_{ij}}\nonumber\\ =& -\sum_{j \neq i} f_{ij} \end{aligned} \] assim, substituindo \[ f_{ij} = 4\epsilon \left[ -12 \dfrac{\sigma^{12}}{r_ij^{13}} - 6\dfrac{\sigma^6}{r_ij^7} \right] \nonumber \] de forma que \[ \vec{F}_i = -24\sum_{j \neq i} \dfrac{\epsilon \sigma^6 \left( \vec{r}_i - \vec{r}_j \right) }{r_{ij}^8}\left[ 1 - 2 \left( \dfrac{\sigma}{r_{ij}}\right)^6\right] \]

Raio de Corte

O potencial de LJ, similar a vários outros, tem alcance infinito, consequentemente se formos considera todas as interações que acontecem entre as partículas do sistema este número cresceria como $N^2/2$, só como exemplo se nosso sistema é formado por $N=1024$ partículas teríamos que avaliar $524.288$ interações, mas se nos truncarmos a interação digamos até $8\overset{\circ}{A}$ o número de vizinhos próximos vai para $50$ o que significa que o total de interações cai para $51200$. Um outro problema com o alcance infinito da interação é que devido às condições periódicas a partícula $i$ poderia interagir com um de suas imagens contudo, há métodos de tratar estas situações.

Dessa forma vemos que seria útil limitar o raio de ação da interação mas isso não pode ser feito para todos os potenciais, se a limitação da interação resulta em um erro enorme no calculo da energia devemos considerar o alcance infinito desta, no caso do potencial de LJ isto não acontece.

Algoritmos de integração

Nas aulas anteriores analisamos alguns potenciais de interação que se mostraram adequados para integrar as equações de movimento. O problema com estes algoritmos é que no caso do melhor deles (algoritmo de Runge-Kutta) o custo computacional resulta alto. A fim de melhorar os problemas encontrados com o RK introduziremos um algoritmo com a precisão similar à do RK que analisamos mas com um custo computacional menor, o algoritmo de que falamos pertence ao tipo de algoritmo denominado de Verlet.

Algoritmo clássico de Verlet

Este algoritmo foi primeiramente utilizado por Delambre em 1791 e redescoberto por Verlet em 1960. Também foi utilizado por Cowell e Crommelin en 1909 para calcular a orbita do cometa Halley e por Stömer em 1907 para o calculo do movimento de partículas num campo magnético.

Para o calculo do algoritmo utilizamos uma expansão em serie de Taylor até IV ordem da posição de uma partícula \[ \begin{aligned} r(t + \Delta t) =& r(t) + v(t) \Delta t + \frac{1}{2}a(t)\Delta T^2 + \frac{1}{6}a(t)\Delta T^3 + O(\Delta t^4) \nonumber\\ r(t - \Delta t) =& r(t) - v(t) \Delta t + \frac{1}{2}a(t)\Delta T^2 - \frac{1}{6}a(t)\Delta T^3 + O(\Delta t^4) \nonumber \end{aligned} \] somando ambas expressões \[ r(t + \Delta t) = 2r(t) - r(t - \Delta t) + a(t)\Delta T^2 + O(\Delta t^4) \] Para o calculo da velocidade observamos que \[ \begin{aligned} r(t + \Delta t) =& r(t) + v(t) \Delta t + O(\Delta t^2) \nonumber\\ r(t - \Delta t) =& r(t) - v(t) \Delta t + O(\Delta t^2) \nonumber \end{aligned} \] de onde \[ v(t) = \frac{r(t + \Delta t) - r(t - \Delta t)}{2\Delta t} \]

Velocity Verlet

Tarefa

Comparar o diagrama de fases do potencial de Lennard-Jones utilizando o potencial de Morse com $\alpha=6,\,4,\,8$

Código do LJ

principal.py

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

      import numpy as np
      import matplotlib.pyplot as plt

      import confInicial as inicial
      import filme as grava
      import forca
      import dinamica


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

      print L

      R  = inicial.fcc(Ncelas, L)

      print 'Maximo', L, np.amax(R)

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

      A  = np.zeros_like(R)
      Ap = np.zeros_like(R)

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

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

      Epot, A = forca.aceleracao(L, R, A)

      nPasos   = int(tempEq / dt)
      invPasos = int(0.1 / dt)

      nome = 'filme.xyz'

      for i in xrange(nPasos):
        R, V = dinamica.velocityVerlet(L, dt, R, V, A, Ap)
        if (i%invPasos == 0.0):
          modV2   = np.sum(V*V, axis=1)
          ecinSum = np.sum(modV2)
          print '%f %f'%(i*dt, ecinSum/(3.0*N))
          grava.xyz(nome, R)

      print Epot      
      

confInicial.py

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

      import numpy as np

      def fcc(Ncelas, L):

        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 = L / 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.5*L):
              R[i,j] = R[i,j] - L
            if (R[i,j] <= -0.5*L):
              R[i,j] = R[i,j] + L

        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.zeros(3, dtype=np.float)
        vcm[:] = np.sum(V, axis=0)
        vcm[:] = vcm[:] / N
        for i in range(3):
          V[:,i] = V[:,i] - vcm[i]

        return V
      

filme.py

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

      import numpy as np

      def xyz(nome, R):
        N = R.shape[0]
        arq = open(nome, 'a')
        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()
      

forca.py

      def aceleracao(L, R, A):
        N      = R.shape[0]
        A[:,:] = 0.0
        Epot   = 0.0
        
        for i in xrange(N-1):
          dr = np.zeros(3, dtype=np.float)
          for j in xrange(i+1, N):
              dr[:] = R[i,:] - R[j,:]
              for k in xrange(3):
                if (dr[k] >= 0.5*L):
                  dr[k] = dr[k] - L
                if (dr[k] <= -0.5*L):
                  dr[k] = dr[k] + L

              dr2    = np.sum(dr*dr)
              invdr2 = 1.0 / dr2
              invdr6 = invdr2*invdr2*invdr2
              
              Epot   = Epot + 4.0 * invdr6 * (invdr6 - 1.0)

              fij    = 48.0 * invdr6*invdr2 * (invdr6 - 0.5)
              
              A[i,:] = A[i,:] + fij * dr[:]
              A[j,:] = A[j,:] - fij * dr[:]

        return Epot, A
      

dinamica.py

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

      import numpy as np
      import forca

      def velocityVerlet(L, dt, R, V, A, Ap):

        R[:,:]  = R[:,:] + dt * V[:,:] + 0.5 * A[:,:] * dt*dt
        Ap[:,:] = A[:,:]
        
        #condicoes periodicas
        for i in xrange(R.shape[0]):
          for j in xrange(3):
            if (R[i,j] >= 0.5*L):
              R[i,j] = R[i,j] - L
            if (R[i,j] <= -0.5*L):
              R[i,j] = R[i,j] + L
        
        Epot, A = forca.aceleracao(L, R, A)
        
        V[:,:]  = V[:,:] + 0.5 * (A[:,:] + Ap[:,:]) * dt

        return R, V
      

Referencias

  1. Computer "Experiments" on Classical Fluids. I Thermodynamical properties of Lennard-Jones molecules. Verlet L. Phys. Rev. vol 159, pag. 98, 1967.
  2. Phase diagrams of Lennard‐Jones fluids. Smit B. J. Chem. Phys. 96, 8639, 1992.
  3. The Structure and Stability of Atomic Liquids: From Clusters to Bulk. Doye J. P. K. e Wales D. J. Science 271, 484, 1996.
  4. Application of the Morse potential function to cubic metals. Girifalco L. A. e Weizer V. G. Phys. Rev. 114, 687, 1959.
  5. Application of Morse potential in the molecular-metallic-framework. Sharma K. S. e Kachhava.
  6. Solution of Morse Potential for Face Centre Cube Using Embedded Atom Method. Abajingin, D. D. Advances in Physics Theories and Applications 8, 36, 2012