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

Simulação - Lennard-Jones

2023-04-21T18:53:33Z

Gustavobopsin:

!--------------------------------------------------------------------!
MODULE globais
IMPLICIT NONE
INTEGER :: Lx,Ly,N,Nf,nhis
REAL(8) :: pi,eps,sigma,sigma2,rho
END MODULE
!--------------------------------------------------------------------
PROGRAM lennard_jones
USE globais
IMPLICIT NONE
INTEGER :: i,t,k,j,m,ngr
INTEGER :: q,semente,seed(8)
INTEGER,ALLOCATABLE :: hist_g(:)
REAL(4) :: t_inicial,t_final,tempo
REAL(8) :: d,r,dt,tmax,p,delg,u
REAL(8),ALLOCATABLE :: x(:),y(:),Vx(:),Vy(:),Fx(:),Fy(:),dist(:),g(:)
REAL(8) :: e_cin,e_pot,Temp,pressao,sum_Temp,sum_pressao
!
CALL CPU_time(t_inicial)
!--------------------------------------------------------------------
! Semente para números aleatórios!
semente = 1234
DO q=0,7
seed(q+1) = semente + q*50
END DO
CALL random_seed(put = seed)
!--------------------------------------------------------------------
pi = 4.0*atan(1.0)
Lx = 25
Ly = 25
N = 144 !Número de partículas!
Nf = 2*N-2 !Nuúmero de graus de liberdade!
nhis = 200 !Número de bins!
d = 2.0 !Distância entre as partículas!
r = d*0.5 !Raio das partículas!
eps = 1.0
sigma = 1.0
sigma2 = sigma*sigma
rho = 1.0*N/(Lx*Ly) !Densidade numérica!
Temp = 1.0 !Temperatura desejada, usada para!
!reescalar as velocidades!
!--------------------------------------------------------------------
! Aloca os vetores
ALLOCATE(x(N),y(N),Vx(N),Vy(N),Fx(N),Fy(N))
ALLOCATE(g(nhis),dist(nhis),hist_g(nhis))
!--------------------------------------------------------------------
! Zera variáveis
Fx(:) = 0.0
Fy(:) = 0.0
x(:) = 0.0
y(:) = 0.0
Vx(:) = 0.0
Vy(:) = 0.0
g(:) = 0.0
hist_g(:) = 0
dist(:) = 0.0
sum_Temp = 0.0
sum_pressao = 0.0
!--------------------------------------------------------------------
! Inicializa posições, velocidades, forças e fdr
CALL inicializacao(x,y,Vx,Vy,r,Temp,e_cin)
CALL posicoes(x,y,r)
CALL forcas(x,y,Fx,Fy,e_pot,p)
CALL fdr(0,x,y,ngr,delg,dist,g)
!--------------------------------------------------------------------
! Loop para integrar as equações de movimento
dt = 0.01
tmax = 100.0/dt
OPEN(10,FILE='dados_particulas.dat',ACTION='write',STATUS='unknown')
OPEN(100,FILE='dados_termo.dat',ACTION='write',STATUS='unknown')
OPEN(200,FILE='dados_fdr.dat',ACTION='write',STATUS='unknown')
DO t=0,int(tmax)
Temp = e_cin/Nf
pressao = rho*Temp + 0.5/(Lx*Ly)*p
! Retira o valor inicial da T e p do cálculo da média, pois é muito diferente
IF (t == 0) THEN
sum_Temp = -Temp
sum_pressao = -pressao
GO TO 20
ELSE
GO TO 20
END IF
20 sum_Temp = sum_Temp + Temp
sum_pressao = sum_pressao + pressao
! Escreve os dados nos arquivos
DO i=1,N
WRITE(10,200) t,i,x(i),y(i),Vx(i),Vy(i),Fx(i),Fy(i)
END DO
WRITE(10,*) ''
WRITE(10,*) ''
200 FORMAT(i5,2x,i3,2x,f10.5,2x,f10.5,2x,f10.5,2x,f10.5,2x,f10.5,2x,f10.5)
WRITE(100,300) t,e_cin,e_pot,Temp,pressao
300 FORMAT(i5,2x,f10.5,2x,f10.5,2x,f10.5,2x,f10.5)
! Algoritmo Velocity-Verlet
DO k=1,10
CALL verlet(x,y,Vx,Vy,Fx,Fy,dt,e_cin,e_pot,p)
END DO
CALL fdr(1,x,y,ngr,delg,dist,g)
END DO
CALL fdr(2,x,y,ngr,delg,dist,g)
! Imprime a fdr
DO m=1,nhis
WRITE(200,400) dist(m),g(m)
END DO
400 FORMAT(f10.5,2x,f10.5)
CLOSE(10)
CLOSE(100)
CLOSE(200)
! Calcula a temperatura média e a pressão média
sum_Temp = sum_Temp/tmax
sum_pressao = sum_pressao/tmax
WRITE(*,600) 'Densidade:',rho
WRITE(*,700) 'Temperatura média:',sum_Temp
WRITE(*,800) 'Pressão média:',sum_pressao
600 FORMAT(a11,1x,f7.4)
700 FORMAT(a19,1x,f7.4)
800 FORMAT(a16,1x,f7.4)
!
CALL CPU_time(t_final)
tempo = t_final - t_inicial
WRITE(*,900) 'Tempo de execução:',tempo,'s'
900 FORMAT(a20,1x,f6.2,1x,a1)
!
END PROGRAM
!---------------------------------------------------------------------
! Inicialização
SUBROUTINE inicializacao(rx,ry,vx,vy,r,Temp,e_cin)
USE globais
IMPLICIT NONE
INTEGER :: num,i
REAL(8) :: d,r,x,y,Temp,e_cin,rand1,rand2
REAL(8) :: sumv2,sum_vx,sum_vy,fs
REAL(8) :: rx(N),ry(N),vx(N),vy(N)
! Inicializa as posições
num = 1
x = r
y = r
DO WHILE (y .lt. Ly-r+0.05)
DO WHILE (x .lt. Lx-r+0.05)
IF (num .le. N) THEN
rx(num) = x
ry(num) = y
num = num + 1
END IF
x = x + 2*r
END DO
x = r
y = y + 2*r
END DO
! Inicializa as velocidades
sumv2 = 0.0
sum_vx = 0.0 ; sum_vy = 0.0
DO i=1,N
CALL random_number(rand1)
CALL random_number(rand2)
vx(i) = rand1 - 0.5
vy(i) = rand2 - 0.5
sumv2 = sumv2 + vx(i)*vx(i) + vy(i)*vy(i) !energia cinética!
sum_vx = sum_vx + vx(i) !veloc x do CM!
sum_vy = sum_vy + vy(i) !veloc y do CM!
END DO
e_cin = 0.5*sumv2
sumv2 = sumv2/N
fs = sqrt(2.0*Temp/sumv2) !Reescala as velocidades!
sum_vx = sum_vx/N !veloc x do CM!
sum_vy = sum_vy/N !veloc y do CM!
DO i=1,N
vx(i) = (vx(i) - sum_vx)*fs
vy(i) = (vy(i) - sum_vy)*fs
END DO
END SUBROUTINE
!----------------------------------------------------------------------
! Implime a configuração inicial do sistema
SUBROUTINE posicoes(x,y,r)
USE globais
IMPLICIT NONE
INTEGER :: k
REAL(8) :: r
REAL(8) :: x(N),y(N)
OPEN(9,FILE='posicoes.dat',ACTION='write',STATUS='unknown')
DO k=1,N
WRITE(9,*) x(k),y(k),r
END DO
CLOSE(9)
END SUBROUTINE
!----------------------------------------------------------------------
! Calcula as forças em um instante de tempo
SUBROUTINE forcas(x,y,fx,fy,e_pot,p)
USE globais
IMPLICIT NONE
INTEGER :: i,j
REAL(8) :: a,r2,dx,dy,e_pot,p
REAL(8) :: x(N),y(N),fx(N),fy(N)
fx(:) = 0.0
fy(:) = 0.0
e_pot = 0.0
p = 0.0
DO i=1,N-1
DO j=i+1,N
! Convenção da imagem mínima
dx = x(i)-x(j)
dy = y(i)-y(j)
dx = dx - Lx*nint(dx/Lx)
dy = dy - Ly*nint(dy/Ly)
r2 = dx*dx + dy*dy
a = 48*eps/sigma2*(sigma2/r2)**4*((sigma2/r2)**3 - 0.5)
fx(i) = fx(i) + a*dx
fy(i) = fy(i) + a*dy
fx(j) = fx(j) - a*dx
fy(j) = fy(j) - a*dy
e_pot = e_pot + 4.0*eps*(sigma2/r2)**3*((sigma2/r2)**3-1.0) !Energia potencial!
p = p + a*(dx*dx + dy*dy) !Pressão!
END DO
END DO
END SUBROUTINE
!---------------------------------------------------------------------
! Implementa a condição de contorno
SUBROUTINE condicao_contorno(x,y)
USE globais
IMPLICIT NONE
INTEGER :: i
REAL(8) :: x(N),y(N)
DO i=1,N
IF (x(i) .gt. Lx) THEN
x(i) = x(i) - Lx
ELSE IF (x(i) .lt. 0) THEN
x(i) = x(i) + Lx
END IF
IF (y(i) .gt. Ly) THEN
y(i) = y(i) - Ly
ELSE IF (y(i) .lt. 0) THEN
y(i) = y(i) + Ly
END IF
END DO
END SUBROUTINE
!----------------------------------------------------------------------
! ALgoritmo Velocity-Verlet
SUBROUTINE verlet(rx,ry,vx,vy,fx,fy,deltat,e_cin,e_pot,p)
USE globais
IMPLICIT NONE
INTEGER :: i
REAL(8) :: deltat,e_cin,e_pot,p
REAL(8) :: rx(N),ry(N),vx(N),vy(N),fx(N),fy(N),g(nhis)
e_cin = 0.0
DO i=1,N
vx(i) = vx(i) + 0.5*fx(i)*deltat
vy(i) = vy(i) + 0.5*fy(i)*deltat
rx(i) = rx(i) + vx(i)*deltat
ry(i) = ry(i) + vy(i)*deltat
END DO
CALL condicao_contorno(rx,ry)
CALL forcas(rx,ry,fx,fy,e_pot,p)
DO i=1,N
vx(i) = vx(i) + 0.5*fx(i)*deltat
vy(i) = vy(i) + 0.5*fy(i)*deltat
e_cin = e_cin + vx(i)*vx(i) + vy(i)*vy(i)
END DO
e_cin = 0.5*e_cin !Energia cinética!
END SUBROUTINE
!---------------------------------------------------------------------
! Calcula a função distribuição radial
SUBROUTINE fdr(switch,x,y,ngr,delg,dist,g)
USE globais
IMPLICIT NONE
INTEGER :: i,j,k,ig,ngr,switch
! INTEGER :: hist_g(nhis)
REAL(8) :: dx,dy,r,r2,delg,vb,nid
REAL(8) :: x(N),y(N),dist(nhis),g(nhis)
IF (switch == 0) THEN !Inicializa!
ngr = 0
delg = Lx/(2.0*nhis) !Tamanho do bin!
dist(:) = 0.0
! hist_g(:) = 0
g(:) = 0.0
ELSE IF (switch == 1) THEN !amostra!
ngr = ngr + 1
DO i=1,N-1
DO j=i+1,N
! Convenção da imagem mínima
dx = x(i)-x(j)
dy = y(i)-y(j)
dx = dx - Lx*nint(dx/Lx)
dy = dy - Ly*nint(dy/Ly)
r2 = dx*dx + dy*dy
r = sqrt(r2)
ig = int(r/delg)
g(ig) = g(ig) + 2
END DO
END DO
ELSE IF (switch == 2) THEN !calcula g(r)!
DO k=1,nhis
dist(k) = delg*(k + 0.5)
vb = (2*k+1)*delg*delg
nid = pi*vb*rho
g(k) = g(k)/(ngr*N*nid) !Normaliza g(r)!
END DO
END IF
END SUBROUTINE

Simulação - Lennard-Jones

2023-04-21T18:53:19Z

Gustavobopsin:

!--------------------------------------------------------------------!
MODULE globais
IMPLICIT NONE
INTEGER :: Lx,Ly,N,Nf,nhis
REAL(8) :: pi,eps,sigma,sigma2,rho
END MODULE
!--------------------------------------------------------------------
PROGRAM lennard_jones
USE globais
IMPLICIT NONE
INTEGER :: i,t,k,j,m,ngr
INTEGER :: q,semente,seed(8)
INTEGER,ALLOCATABLE :: hist_g(:)
REAL(4) :: t_inicial,t_final,tempo
REAL(8) :: d,r,dt,tmax,p,delg,u
REAL(8),ALLOCATABLE :: x(:),y(:),Vx(:),Vy(:),Fx(:),Fy(:),dist(:),g(:)
REAL(8) :: e_cin,e_pot,Temp,pressao,sum_Temp,sum_pressao
!
CALL CPU_time(t_inicial)
!--------------------------------------------------------------------
! Semente para números aleatórios
semente = 1234
DO q=0,7
seed(q+1) = semente + q*50
END DO
CALL random_seed(put = seed)
!--------------------------------------------------------------------
pi = 4.0*atan(1.0)
Lx = 25
Ly = 25
N = 144 !Número de partículas!
Nf = 2*N-2 !Nuúmero de graus de liberdade!
nhis = 200 !Número de bins!
d = 2.0 !Distância entre as partículas!
r = d*0.5 !Raio das partículas!
eps = 1.0
sigma = 1.0
sigma2 = sigma*sigma
rho = 1.0*N/(Lx*Ly) !Densidade numérica!
Temp = 1.0 !Temperatura desejada, usada para!
!reescalar as velocidades!
!--------------------------------------------------------------------
! Aloca os vetores
ALLOCATE(x(N),y(N),Vx(N),Vy(N),Fx(N),Fy(N))
ALLOCATE(g(nhis),dist(nhis),hist_g(nhis))
!--------------------------------------------------------------------
! Zera variáveis
Fx(:) = 0.0
Fy(:) = 0.0
x(:) = 0.0
y(:) = 0.0
Vx(:) = 0.0
Vy(:) = 0.0
g(:) = 0.0
hist_g(:) = 0
dist(:) = 0.0
sum_Temp = 0.0
sum_pressao = 0.0
!--------------------------------------------------------------------
! Inicializa posições, velocidades, forças e fdr
CALL inicializacao(x,y,Vx,Vy,r,Temp,e_cin)
CALL posicoes(x,y,r)
CALL forcas(x,y,Fx,Fy,e_pot,p)
CALL fdr(0,x,y,ngr,delg,dist,g)
!--------------------------------------------------------------------
! Loop para integrar as equações de movimento
dt = 0.01
tmax = 100.0/dt
OPEN(10,FILE='dados_particulas.dat',ACTION='write',STATUS='unknown')
OPEN(100,FILE='dados_termo.dat',ACTION='write',STATUS='unknown')
OPEN(200,FILE='dados_fdr.dat',ACTION='write',STATUS='unknown')
DO t=0,int(tmax)
Temp = e_cin/Nf
pressao = rho*Temp + 0.5/(Lx*Ly)*p
! Retira o valor inicial da T e p do cálculo da média, pois é muito diferente
IF (t == 0) THEN
sum_Temp = -Temp
sum_pressao = -pressao
GO TO 20
ELSE
GO TO 20
END IF
20 sum_Temp = sum_Temp + Temp
sum_pressao = sum_pressao + pressao
! Escreve os dados nos arquivos
DO i=1,N
WRITE(10,200) t,i,x(i),y(i),Vx(i),Vy(i),Fx(i),Fy(i)
END DO
WRITE(10,*) ''
WRITE(10,*) ''
200 FORMAT(i5,2x,i3,2x,f10.5,2x,f10.5,2x,f10.5,2x,f10.5,2x,f10.5,2x,f10.5)
WRITE(100,300) t,e_cin,e_pot,Temp,pressao
300 FORMAT(i5,2x,f10.5,2x,f10.5,2x,f10.5,2x,f10.5)
! Algoritmo Velocity-Verlet
DO k=1,10
CALL verlet(x,y,Vx,Vy,Fx,Fy,dt,e_cin,e_pot,p)
END DO
CALL fdr(1,x,y,ngr,delg,dist,g)
END DO
CALL fdr(2,x,y,ngr,delg,dist,g)
! Imprime a fdr
DO m=1,nhis
WRITE(200,400) dist(m),g(m)
END DO
400 FORMAT(f10.5,2x,f10.5)
CLOSE(10)
CLOSE(100)
CLOSE(200)
! Calcula a temperatura média e a pressão média
sum_Temp = sum_Temp/tmax
sum_pressao = sum_pressao/tmax
WRITE(*,600) 'Densidade:',rho
WRITE(*,700) 'Temperatura média:',sum_Temp
WRITE(*,800) 'Pressão média:',sum_pressao
600 FORMAT(a11,1x,f7.4)
700 FORMAT(a19,1x,f7.4)
800 FORMAT(a16,1x,f7.4)
!
CALL CPU_time(t_final)
tempo = t_final - t_inicial
WRITE(*,900) 'Tempo de execução:',tempo,'s'
900 FORMAT(a20,1x,f6.2,1x,a1)
!
END PROGRAM
!---------------------------------------------------------------------
! Inicialização
SUBROUTINE inicializacao(rx,ry,vx,vy,r,Temp,e_cin)
USE globais
IMPLICIT NONE
INTEGER :: num,i
REAL(8) :: d,r,x,y,Temp,e_cin,rand1,rand2
REAL(8) :: sumv2,sum_vx,sum_vy,fs
REAL(8) :: rx(N),ry(N),vx(N),vy(N)
! Inicializa as posições
num = 1
x = r
y = r
DO WHILE (y .lt. Ly-r+0.05)
DO WHILE (x .lt. Lx-r+0.05)
IF (num .le. N) THEN
rx(num) = x
ry(num) = y
num = num + 1
END IF
x = x + 2*r
END DO
x = r
y = y + 2*r
END DO
! Inicializa as velocidades
sumv2 = 0.0
sum_vx = 0.0 ; sum_vy = 0.0
DO i=1,N
CALL random_number(rand1)
CALL random_number(rand2)
vx(i) = rand1 - 0.5
vy(i) = rand2 - 0.5
sumv2 = sumv2 + vx(i)*vx(i) + vy(i)*vy(i) !energia cinética!
sum_vx = sum_vx + vx(i) !veloc x do CM!
sum_vy = sum_vy + vy(i) !veloc y do CM!
END DO
e_cin = 0.5*sumv2
sumv2 = sumv2/N
fs = sqrt(2.0*Temp/sumv2) !Reescala as velocidades!
sum_vx = sum_vx/N !veloc x do CM!
sum_vy = sum_vy/N !veloc y do CM!
DO i=1,N
vx(i) = (vx(i) - sum_vx)*fs
vy(i) = (vy(i) - sum_vy)*fs
END DO
END SUBROUTINE
!----------------------------------------------------------------------
! Implime a configuração inicial do sistema
SUBROUTINE posicoes(x,y,r)
USE globais
IMPLICIT NONE
INTEGER :: k
REAL(8) :: r
REAL(8) :: x(N),y(N)
OPEN(9,FILE='posicoes.dat',ACTION='write',STATUS='unknown')
DO k=1,N
WRITE(9,*) x(k),y(k),r
END DO
CLOSE(9)
END SUBROUTINE
!----------------------------------------------------------------------
! Calcula as forças em um instante de tempo
SUBROUTINE forcas(x,y,fx,fy,e_pot,p)
USE globais
IMPLICIT NONE
INTEGER :: i,j
REAL(8) :: a,r2,dx,dy,e_pot,p
REAL(8) :: x(N),y(N),fx(N),fy(N)
fx(:) = 0.0
fy(:) = 0.0
e_pot = 0.0
p = 0.0
DO i=1,N-1
DO j=i+1,N
! Convenção da imagem mínima
dx = x(i)-x(j)
dy = y(i)-y(j)
dx = dx - Lx*nint(dx/Lx)
dy = dy - Ly*nint(dy/Ly)
r2 = dx*dx + dy*dy
a = 48*eps/sigma2*(sigma2/r2)**4*((sigma2/r2)**3 - 0.5)
fx(i) = fx(i) + a*dx
fy(i) = fy(i) + a*dy
fx(j) = fx(j) - a*dx
fy(j) = fy(j) - a*dy
e_pot = e_pot + 4.0*eps*(sigma2/r2)**3*((sigma2/r2)**3-1.0) !Energia potencial!
p = p + a*(dx*dx + dy*dy) !Pressão!
END DO
END DO
END SUBROUTINE
!---------------------------------------------------------------------
! Implementa a condição de contorno
SUBROUTINE condicao_contorno(x,y)
USE globais
IMPLICIT NONE
INTEGER :: i
REAL(8) :: x(N),y(N)
DO i=1,N
IF (x(i) .gt. Lx) THEN
x(i) = x(i) - Lx
ELSE IF (x(i) .lt. 0) THEN
x(i) = x(i) + Lx
END IF
IF (y(i) .gt. Ly) THEN
y(i) = y(i) - Ly
ELSE IF (y(i) .lt. 0) THEN
y(i) = y(i) + Ly
END IF
END DO
END SUBROUTINE
!----------------------------------------------------------------------
! ALgoritmo Velocity-Verlet
SUBROUTINE verlet(rx,ry,vx,vy,fx,fy,deltat,e_cin,e_pot,p)
USE globais
IMPLICIT NONE
INTEGER :: i
REAL(8) :: deltat,e_cin,e_pot,p
REAL(8) :: rx(N),ry(N),vx(N),vy(N),fx(N),fy(N),g(nhis)
e_cin = 0.0
DO i=1,N
vx(i) = vx(i) + 0.5*fx(i)*deltat
vy(i) = vy(i) + 0.5*fy(i)*deltat
rx(i) = rx(i) + vx(i)*deltat
ry(i) = ry(i) + vy(i)*deltat
END DO
CALL condicao_contorno(rx,ry)
CALL forcas(rx,ry,fx,fy,e_pot,p)
DO i=1,N
vx(i) = vx(i) + 0.5*fx(i)*deltat
vy(i) = vy(i) + 0.5*fy(i)*deltat
e_cin = e_cin + vx(i)*vx(i) + vy(i)*vy(i)
END DO
e_cin = 0.5*e_cin !Energia cinética!
END SUBROUTINE
!---------------------------------------------------------------------
! Calcula a função distribuição radial
SUBROUTINE fdr(switch,x,y,ngr,delg,dist,g)
USE globais
IMPLICIT NONE
INTEGER :: i,j,k,ig,ngr,switch
! INTEGER :: hist_g(nhis)
REAL(8) :: dx,dy,r,r2,delg,vb,nid
REAL(8) :: x(N),y(N),dist(nhis),g(nhis)
IF (switch == 0) THEN !Inicializa!
ngr = 0
delg = Lx/(2.0*nhis) !Tamanho do bin!
dist(:) = 0.0
! hist_g(:) = 0
g(:) = 0.0
ELSE IF (switch == 1) THEN !amostra!
ngr = ngr + 1
DO i=1,N-1
DO j=i+1,N
! Convenção da imagem mínima
dx = x(i)-x(j)
dy = y(i)-y(j)
dx = dx - Lx*nint(dx/Lx)
dy = dy - Ly*nint(dy/Ly)
r2 = dx*dx + dy*dy
r = sqrt(r2)
ig = int(r/delg)
g(ig) = g(ig) + 2
END DO
END DO
ELSE IF (switch == 2) THEN !calcula g(r)!
DO k=1,nhis
dist(k) = delg*(k + 0.5)
vb = (2*k+1)*delg*delg
nid = pi*vb*rho
g(k) = g(k)/(ngr*N*nid) !Normaliza g(r)!
END DO
END IF
END SUBROUTINE

Arquivo:Fdr.png

2023-04-19T01:17:41Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Fdr.png"

Arquivo:Fdr.png

2023-04-19T00:56:23Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Fdr.png"

Arquivo:Animacao2.gif

2023-04-19T00:56:04Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Animacao2.gif"

Arquivo:Animacao1.gif

2023-04-19T00:55:38Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Animacao1.gif"

Lennard-Jones e propriedades termodinâmicas

2023-04-19T00:53:09Z

Gustavobopsin: /* Simulação */

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|500px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|500px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 144, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.2304</math> e consideramos a massa das partículas como 1. A temperatura desejada que utilizamos foi de T=1.0. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os últimos 100 instantes de tempo, passo à passo e os últimos 100 instantes de tempo, de 5 em 5 passos de tempo. Quando vemos a simulação passo à passo, estamos olhando para uma visão mais microscópica (mecânica estatística), enquanto estamos vendo a simulação de 5 em 5 passo, estamos observando um comportamento mais macroscópico (termodinâmica).
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo, passo a passo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo, de 5 em 5 passos.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 300 passos da simulação.|500px]]
|}
Vemos que como o sistema já está equilibrado para uma quantidade de menos de 300 passos, plotamos à direita a série temporal contendo as 300 primeiras medidas para conseguirmos ver o começo da simulação que claramente está fora do equilíbrio.

A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 300 passos da simulação.|500px]]
|}
Podemos ver que a temperatura do sistema também já se equilibrou nos primeiros 300 passos e por isso queremos ver o comportamento do sistema fora do equilíbrio no começo da simulação. A partir do equilíbrio, poderíamos utilizar a simulação para calcular propriedades termodinâmicas.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá. Novamente, a pressão equilibra muito cedo e podemos ver que para um tempo menor que 300.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ainda calcular a função distribuição radial <math>g(\boldsymbol{r})</math> para esta simulação.
[[Arquivo:fdr.png|thumb|upright=1.0|center|Função distribuição radial.|500px]]

== Códigos ==
[[Simulação - Lennard-Jones]]

[[Figuras e animações]]

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Figuras e animações

2023-04-19T00:51:38Z

Gustavobopsin:

<source lang = 'gnuplot'>
reset
set grid
set xlabel 'x'
set ylabel 'y'
set xrange[0:25]
set yrange[0:25]
set size ratio 1

#Inicio das posições
set term png
set output 'inicializacao.png'
plot 'posicoes.dat' w circles fill solid 0.6 border 1 lc 1 t ''

#Animações
set key out
set term gif animate delay 10
set output 'animacao1.gif'
do for [i=9900:10000] {plot 'dados_particulas.dat' index i u 3:4:(1.0) w circles fill solid 0.6 border 1 t 'tempo='.i}

set term gif animate delay 50
set output 'animacao2.gif'
do for [i=9900:10000:5] {plot 'dados_particulas.dat' index i u 3:4:(1.0) w circles fill solid 0.6 border 1 t 'tempo='.i}
reset

set grid
set xlabel 't'
set ylabel 'E'
set key out
set term png size 800,500
set title 'Energia mecânica: 124.138'
set output 'energia_tempo.png'
plot 'dados_termo.dat' u 1:2 w l lc 6 t 'Energia cinética' , \
'dados_termo.dat' u 1:3 w l lc 7 t 'Energia potencial' , \
'dados_termo.dat' u 1:($2+$3) w l lc 1 t 'Energia mecânica'

unset title
set xrange [0:300]
set output 'energia_tempo_zoom.png'
plot 'dados_termo.dat' u 1:2 w l lc 6 t 'Energia cinética' , \
'dados_termo.dat' u 1:3 w l lc 7 t 'Energia potencial' , \
'dados_termo.dat' u 1:($2+$3) w l lc 1 t 'Energia mecânica'

unset xrange
set title 'Temperatura média: 0.7237'
set ylabel 'T'
set output 'temperatura_tempo.png'
plot 'dados_termo.dat' u 1:4 w l lc 7 t ''
unset title
set xrange [0:300]
set yrange [0.1:0.8]
set output 'temperatura_tempo_zoom.png'
plot 'dados_termo.dat' u 1:4 w l lc 7 t ''

unset yrange
unset xrange
set title 'Pressão média: 0.1877'
set ylabel 'P'
set output 'pressao_tempo.png'
plot 'dados_termo.dat' u 1:5 w l lc 6 t ''
unset title
set xrange [0:300]
set output 'pressao_tempo_zoom.png'
plot 'dados_termo.dat' u 1:5 w l lc 6 t ''

unset xrange
unset yrange
set title 'Número de bins = 200'
set ylabel 'g(r)'
set xlabel 'r'
set output 'fdr.png'
plot 'dados_fdr.dat' u 1:2 w l lw 2 lc 8 t ''
reset
</source>

Figuras e animações

2023-04-19T00:50:20Z

Gustavobopsin: Criou página com 'reset set grid set xlabel 'x' set ylabel 'y' set xrange[0:25] set yrange[0:25] set size ratio 1 #Inicio das posições set term png set output 'inicializacao.png' plot 'posic...'

reset
set grid
set xlabel 'x'
set ylabel 'y'
set xrange[0:25]
set yrange[0:25]
set size ratio 1

#Inicio das posições
set term png
set output 'inicializacao.png'
plot 'posicoes.dat' w circles fill solid 0.6 border 1 lc 1 t ''

#Animações
set key out
set term gif animate delay 10
set output 'animacao1.gif'
do for [i=9900:10000] {plot 'dados_particulas.dat' index i u 3:4:(1.0) w circles fill solid 0.6 border 1 t 'tempo='.i}

set term gif animate delay 50
set output 'animacao2.gif'
do for [i=9900:10000:5] {plot 'dados_particulas.dat' index i u 3:4:(1.0) w circles fill solid 0.6 border 1 t 'tempo='.i}
reset

set grid
set xlabel 't'
set ylabel 'E'
set key out
set term png size 800,500
set title 'Energia mecânica: 124.138'
set output 'energia_tempo.png'
plot 'dados_termo.dat' u 1:2 w l lc 6 t 'Energia cinética' , \
'dados_termo.dat' u 1:3 w l lc 7 t 'Energia potencial' , \
'dados_termo.dat' u 1:($2+$3) w l lc 1 t 'Energia mecânica'

unset title
set xrange [0:300]
set output 'energia_tempo_zoom.png'
plot 'dados_termo.dat' u 1:2 w l lc 6 t 'Energia cinética' , \
'dados_termo.dat' u 1:3 w l lc 7 t 'Energia potencial' , \
'dados_termo.dat' u 1:($2+$3) w l lc 1 t 'Energia mecânica'

unset xrange
set title 'Temperatura média: 0.7237'
set ylabel 'T'
set output 'temperatura_tempo.png'
plot 'dados_termo.dat' u 1:4 w l lc 7 t ''
unset title
set xrange [0:300]
set yrange [0.1:0.8]
set output 'temperatura_tempo_zoom.png'
plot 'dados_termo.dat' u 1:4 w l lc 7 t ''

unset yrange
unset xrange
set title 'Pressão média: 0.1877'
set ylabel 'P'
set output 'pressao_tempo.png'
plot 'dados_termo.dat' u 1:5 w l lc 6 t ''
unset title
set xrange [0:300]
set output 'pressao_tempo_zoom.png'
plot 'dados_termo.dat' u 1:5 w l lc 6 t ''

unset xrange
unset yrange
set title 'Número de bins = 200'
set ylabel 'g(r)'
set xlabel 'r'
set output 'fdr.png'
plot 'dados_fdr.dat' u 1:2 w l lw 2 lc 8 t ''
reset

Arquivo:Animacao2.gif

2023-04-19T00:47:47Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Animacao2.gif"

Arquivo:Animacao1.gif

2023-04-19T00:47:19Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Animacao1.gif"

Arquivo:Fdr.png

2023-04-19T00:46:49Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Fdr.png"

Simulação - Lennard-Jones

2023-04-19T00:45:44Z

Gustavobopsin: Criou página com '!-------------------------------------------------------------------- MODULE globais IMPLICIT NONE INTEGER :: Lx,Ly,N,Nf,nhis REAL(8) :: pi,eps,sig...'

!--------------------------------------------------------------------
MODULE globais
IMPLICIT NONE
INTEGER :: Lx,Ly,N,Nf,nhis
REAL(8) :: pi,eps,sigma,sigma2,rho
END MODULE
!--------------------------------------------------------------------
PROGRAM lennard_jones
USE globais
IMPLICIT NONE
INTEGER :: i,t,k,j,m,ngr
INTEGER :: q,semente,seed(8)
INTEGER,ALLOCATABLE :: hist_g(:)
REAL(4) :: t_inicial,t_final,tempo
REAL(8) :: d,r,dt,tmax,p,delg,u
REAL(8),ALLOCATABLE :: x(:),y(:),Vx(:),Vy(:),Fx(:),Fy(:),dist(:),g(:)
REAL(8) :: e_cin,e_pot,Temp,pressao,sum_Temp,sum_pressao
!
CALL CPU_time(t_inicial)
!--------------------------------------------------------------------
! Semente para números aleatórios
semente = 1234
DO q=0,7
seed(q+1) = semente + q*50
END DO
CALL random_seed(put = seed)
!--------------------------------------------------------------------
pi = 4.0*atan(1.0)
Lx = 25
Ly = 25
N = 144 !Número de partículas!
Nf = 2*N-2 !Nuúmero de graus de liberdade!
nhis = 200 !Número de bins!
d = 2.0 !Distância entre as partículas!
r = d*0.5 !Raio das partículas!
eps = 1.0
sigma = 1.0
sigma2 = sigma*sigma
rho = 1.0*N/(Lx*Ly) !Densidade numérica!
Temp = 1.0 !Temperatura desejada, usada para!
!reescalar as velocidades!
!--------------------------------------------------------------------
! Aloca os vetores
ALLOCATE(x(N),y(N),Vx(N),Vy(N),Fx(N),Fy(N))
ALLOCATE(g(nhis),dist(nhis),hist_g(nhis))
!--------------------------------------------------------------------
! Zera variáveis
Fx(:) = 0.0
Fy(:) = 0.0
x(:) = 0.0
y(:) = 0.0
Vx(:) = 0.0
Vy(:) = 0.0
g(:) = 0.0
hist_g(:) = 0
dist(:) = 0.0
sum_Temp = 0.0
sum_pressao = 0.0
!--------------------------------------------------------------------
! Inicializa posições, velocidades, forças e fdr
CALL inicializacao(x,y,Vx,Vy,r,Temp,e_cin)
CALL posicoes(x,y,r)
CALL forcas(x,y,Fx,Fy,e_pot,p)
CALL fdr(0,x,y,ngr,delg,dist,g)
!--------------------------------------------------------------------
! Loop para integrar as equações de movimento
dt = 0.01
tmax = 100.0/dt
OPEN(10,FILE='dados_particulas.dat',ACTION='write',STATUS='unknown')
OPEN(100,FILE='dados_termo.dat',ACTION='write',STATUS='unknown')
OPEN(200,FILE='dados_fdr.dat',ACTION='write',STATUS='unknown')
DO t=0,int(tmax)
Temp = e_cin/Nf
pressao = rho*Temp + 0.5/(Lx*Ly)*p
! Retira o valor inicial da T e p do cálculo da média, pois é muito diferente
IF (t == 0) THEN
sum_Temp = -Temp
sum_pressao = -pressao
GO TO 20
ELSE
GO TO 20
END IF
20 sum_Temp = sum_Temp + Temp
sum_pressao = sum_pressao + pressao
! Escreve os dados nos arquivos
DO i=1,N
WRITE(10,200) t,i,x(i),y(i),Vx(i),Vy(i),Fx(i),Fy(i)
END DO
WRITE(10,*) ''
WRITE(10,*) ''
200 FORMAT(i5,2x,i3,2x,f10.5,2x,f10.5,2x,f10.5,2x,f10.5,2x,f10.5,2x,f10.5)
WRITE(100,300) t,e_cin,e_pot,Temp,pressao
300 FORMAT(i5,2x,f10.5,2x,f10.5,2x,f10.5,2x,f10.5)
! Algoritmo Velocity-Verlet
DO k=1,10
CALL verlet(x,y,Vx,Vy,Fx,Fy,dt,e_cin,e_pot,p)
END DO
CALL fdr(1,x,y,ngr,delg,dist,g)
END DO
CALL fdr(2,x,y,ngr,delg,dist,g)
! Imprime a fdr
DO m=1,nhis
WRITE(200,400) dist(m),g(m)
END DO
400 FORMAT(f10.5,2x,f10.5)
CLOSE(10)
CLOSE(100)
CLOSE(200)
! Calcula a temperatura média e a pressão média
sum_Temp = sum_Temp/tmax
sum_pressao = sum_pressao/tmax
WRITE(*,600) 'Densidade:',rho
WRITE(*,700) 'Temperatura média:',sum_Temp
WRITE(*,800) 'Pressão média:',sum_pressao
600 FORMAT(a11,1x,f7.4)
700 FORMAT(a19,1x,f7.4)
800 FORMAT(a16,1x,f7.4)
!
CALL CPU_time(t_final)
tempo = t_final - t_inicial
WRITE(*,900) 'Tempo de execução:',tempo,'s'
900 FORMAT(a20,1x,f6.2,1x,a1)
!
END PROGRAM
!---------------------------------------------------------------------
! Inicialização
SUBROUTINE inicializacao(rx,ry,vx,vy,r,Temp,e_cin)
USE globais
IMPLICIT NONE
INTEGER :: num,i
REAL(8) :: d,r,x,y,Temp,e_cin,rand1,rand2
REAL(8) :: sumv2,sum_vx,sum_vy,fs
REAL(8) :: rx(N),ry(N),vx(N),vy(N)
! Inicializa as posições
num = 1
x = r
y = r
DO WHILE (y .lt. Ly-r+0.05)
DO WHILE (x .lt. Lx-r+0.05)
IF (num .le. N) THEN
rx(num) = x
ry(num) = y
num = num + 1
END IF
x = x + 2*r
END DO
x = r
y = y + 2*r
END DO
! Inicializa as velocidades
sumv2 = 0.0
sum_vx = 0.0 ; sum_vy = 0.0
DO i=1,N
CALL random_number(rand1)
CALL random_number(rand2)
vx(i) = rand1 - 0.5
vy(i) = rand2 - 0.5
sumv2 = sumv2 + vx(i)*vx(i) + vy(i)*vy(i) !energia cinética!
sum_vx = sum_vx + vx(i) !veloc x do CM!
sum_vy = sum_vy + vy(i) !veloc y do CM!
END DO
e_cin = 0.5*sumv2
sumv2 = sumv2/N
fs = sqrt(2.0*Temp/sumv2) !Reescala as velocidades!
sum_vx = sum_vx/N !veloc x do CM!
sum_vy = sum_vy/N !veloc y do CM!
DO i=1,N
vx(i) = (vx(i) - sum_vx)*fs
vy(i) = (vy(i) - sum_vy)*fs
END DO
END SUBROUTINE
!----------------------------------------------------------------------
! Implime a configuração inicial do sistema
SUBROUTINE posicoes(x,y,r)
USE globais
IMPLICIT NONE
INTEGER :: k
REAL(8) :: r
REAL(8) :: x(N),y(N)
OPEN(9,FILE='posicoes.dat',ACTION='write',STATUS='unknown')
DO k=1,N
WRITE(9,*) x(k),y(k),r
END DO
CLOSE(9)
END SUBROUTINE
!----------------------------------------------------------------------
! Calcula as forças em um instante de tempo
SUBROUTINE forcas(x,y,fx,fy,e_pot,p)
USE globais
IMPLICIT NONE
INTEGER :: i,j
REAL(8) :: a,r2,dx,dy,e_pot,p
REAL(8) :: x(N),y(N),fx(N),fy(N)
fx(:) = 0.0
fy(:) = 0.0
e_pot = 0.0
p = 0.0
DO i=1,N-1
DO j=i+1,N
! Convenção da imagem mínima
dx = x(i)-x(j)
dy = y(i)-y(j)
dx = dx - Lx*nint(dx/Lx)
dy = dy - Ly*nint(dy/Ly)
r2 = dx*dx + dy*dy
a = 48*eps/sigma2*(sigma2/r2)**4*((sigma2/r2)**3 - 0.5)
fx(i) = fx(i) + a*dx
fy(i) = fy(i) + a*dy
fx(j) = fx(j) - a*dx
fy(j) = fy(j) - a*dy
e_pot = e_pot + 4.0*eps*(sigma2/r2)**3*((sigma2/r2)**3-1.0) !Energia potencial!
p = p + a*(dx*dx + dy*dy) !Pressão!
END DO
END DO
END SUBROUTINE
!---------------------------------------------------------------------
! Implementa a condição de contorno
SUBROUTINE condicao_contorno(x,y)
USE globais
IMPLICIT NONE
INTEGER :: i
REAL(8) :: x(N),y(N)
DO i=1,N
IF (x(i) .gt. Lx) THEN
x(i) = x(i) - Lx
ELSE IF (x(i) .lt. 0) THEN
x(i) = x(i) + Lx
END IF
IF (y(i) .gt. Ly) THEN
y(i) = y(i) - Ly
ELSE IF (y(i) .lt. 0) THEN
y(i) = y(i) + Ly
END IF
END DO
END SUBROUTINE
!----------------------------------------------------------------------
! ALgoritmo Velocity-Verlet
SUBROUTINE verlet(rx,ry,vx,vy,fx,fy,deltat,e_cin,e_pot,p)
USE globais
IMPLICIT NONE
INTEGER :: i
REAL(8) :: deltat,e_cin,e_pot,p
REAL(8) :: rx(N),ry(N),vx(N),vy(N),fx(N),fy(N),g(nhis)
e_cin = 0.0
DO i=1,N
vx(i) = vx(i) + 0.5*fx(i)*deltat
vy(i) = vy(i) + 0.5*fy(i)*deltat
rx(i) = rx(i) + vx(i)*deltat
ry(i) = ry(i) + vy(i)*deltat
END DO
CALL condicao_contorno(rx,ry)
CALL forcas(rx,ry,fx,fy,e_pot,p)
DO i=1,N
vx(i) = vx(i) + 0.5*fx(i)*deltat
vy(i) = vy(i) + 0.5*fy(i)*deltat
e_cin = e_cin + vx(i)*vx(i) + vy(i)*vy(i)
END DO
e_cin = 0.5*e_cin !Energia cinética!
END SUBROUTINE
!---------------------------------------------------------------------
! Calcula a função distribuição radial
SUBROUTINE fdr(switch,x,y,ngr,delg,dist,g)
USE globais
IMPLICIT NONE
INTEGER :: i,j,k,ig,ngr,switch
! INTEGER :: hist_g(nhis)
REAL(8) :: dx,dy,r,r2,delg,vb,nid
REAL(8) :: x(N),y(N),dist(nhis),g(nhis)
IF (switch == 0) THEN !Inicializa!
ngr = 0
delg = Lx/(2.0*nhis) !Tamanho do bin!
dist(:) = 0.0
! hist_g(:) = 0
g(:) = 0.0
ELSE IF (switch == 1) THEN !amostra!
ngr = ngr + 1
DO i=1,N-1
DO j=i+1,N
! Convenção da imagem mínima
dx = x(i)-x(j)
dy = y(i)-y(j)
dx = dx - Lx*nint(dx/Lx)
dy = dy - Ly*nint(dy/Ly)
r2 = dx*dx + dy*dy
r = sqrt(r2)
ig = int(r/delg)
g(ig) = g(ig) + 2
END DO
END DO
ELSE IF (switch == 2) THEN !calcula g(r)!
DO k=1,nhis
dist(k) = delg*(k + 0.5)
vb = (2*k+1)*delg*delg
nid = pi*vb*rho
g(k) = g(k)/(ngr*N*nid) !Normaliza g(r)!
END DO
END IF
END SUBROUTINE

Arquivo:Pressao tempo zoom.png

2023-04-19T00:44:27Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Pressao tempo zoom.png"

Arquivo:Pressao tempo.png

2023-04-19T00:44:07Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Pressao tempo.png"

Arquivo:Temperatura tempo zoom.png

2023-04-19T00:43:51Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Temperatura tempo zoom.png"

Arquivo:Temperatura tempo.png

2023-04-19T00:43:31Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Temperatura tempo.png"

Arquivo:Energia tempo zoom.png

2023-04-19T00:43:12Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Energia tempo zoom.png"

Arquivo:Energia tempo zoom.png

2023-04-19T00:43:11Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Energia tempo zoom.png"

Arquivo:Energia tempo.png

2023-04-19T00:42:48Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Energia tempo.png"

Arquivo:Fdr.png

2023-04-19T00:39:33Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Fdr.png"

Lennard-Jones e propriedades termodinâmicas

2023-04-19T00:30:39Z

Gustavobopsin: /* Simulação */

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|500px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|500px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 144, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.2304</math> e consideramos a massa das partículas como 1. A temperatura desejada que utilizamos foi de T=1.0. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os últimos 100 instantes de tempo, passo à passo e os últimos 100 instantes de tempo, de 5 em 5 passos de tempo. Quando vemos a simulação passo à passo, estamos olhando para uma visão mais microscópica (mecânica estatística), enquanto estamos vendo a simulação de 5 em 5 passo, estamos observando um comportamento mais macroscópico (termodinâmica).
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo, passo a passo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo, de 5 em 5 passos.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 300 passos da simulação.|500px]]
Vemos que como o sistema já está equilibrado para uma quantidade de menos de 300 passos, plotamos à direita a série temporal contendo as 300 primeiras medidas para conseguirmos ver o começo da simulação que claramente está fora do equilíbrio.
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 300 passos da simulação.|500px]]
|}
Podemos ver que a temperatura do sistema também já se equilibrou nos primeiros 300 passos e por isso queremos ver o comportamento do sistema fora do equilíbrio no começo da simulação. A partir do equilíbrio, poderíamos utilizar a simulação para calcular propriedades termodinâmicas.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá. Novamente, a pressão equilibra muito cedo e podemos ver que para um tempo menor que 300.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ainda calcular a função distribuição radial <math>g(\boldsymbol{r})</math> para esta simulação.
[[Arquivo:fdr.png|thumb|upright=1.0|center|Função distribuição radial.|500px]]

== Códigos ==
[[Simulação - Lennard-Jones]]

[[Figuras e animações]]

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Lennard-Jones e propriedades termodinâmicas

2023-04-18T23:54:56Z

Gustavobopsin: /* Simulação */

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|500px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|500px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 144, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.724</math> e consideramos a massa das partículas como 1. A temperatura desejada que utilizamos foi de T=1.0. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os últimos 100 instantes de tempo, passo à passo e os últimos 100 instantes de tempo, de 5 em 5 passos de tempo. Quando vemos a simulação passo à passo, estamos olhando para uma visão mais microscópica (mecânica estatística), enquanto estamos vendo a simulação de 5 em 5 passo, estamos observando um comportamento mais macroscópico (termodinâmica).
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo, passo a passo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo, de 5 em 5 passos.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 300 passos da simulação.|500px]]
Vemos que como o sistema já está equilibrado para uma quantidade de menos de 300 passos, plotamos à direita a série temporal contendo as 300 primeiras medidas para conseguirmos ver o começo da simulação que claramente está fora do equilíbrio.
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 300 passos da simulação.|500px]]
|}
Podemos ver que a temperatura do sistema também já se equilibrou nos primeiros 300 passos e por isso queremos ver o comportamento do sistema fora do equilíbrio no começo da simulação. A partir do equilíbrio, poderíamos utilizar a simulação para calcular propriedades termodinâmicas.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá. Novamente, a pressão equilibra muito cedo e podemos ver que para um tempo menor que 300.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ainda calcular a função distribuição radial <math>g(\boldsymbol{r})</math> para esta simulação.
[[Arquivo:fdr.png|thumb|upright=1.0|center|Função distribuição radial.|500px]]

== Códigos ==
[[Simulação - Lennard-Jones]]

[[Figuras e animações]]

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Lennard-Jones e propriedades termodinâmicas

2023-04-18T23:42:33Z

Gustavobopsin: /* Simulação */

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|500px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|500px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 144, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.724</math> e consideramos a massa das partículas como 1. A temperatura desejada que utilizamos foi de T=1.0. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os últimos 100 instantes de tempo, passo à passo e os últimos 100 instantes de tempo, de 5 em 5 passos de tempo. Quando vemos a simulação passo à passo, estamos olhando para uma visão mais microscópica (mecânica estatística), enquanto estamos vendo a simulação de 5 em 5 passo, estamos observando um comportamento mais macroscópico (termodinâmica).
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo, passo a passo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo, de 5 em 5 passos.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 1000 passos da simulação.|500px]]
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ver que a temperatura do sistema já se equilibrou nos primeiros 1000 passos. Poderíamos, a partir daí, já utilizar a simulação para calcular propriedades termodinâmicas no equilíbrio.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ainda calcular a função distribuição radial <math>g(\boldsymbol{r})</math> para esta simulação.
[[Arquivo:fdr.png|thumb|upright=1.0|center|Função distribuição radial.|500px]]

== Códigos ==
[[Simulação - Lennard-Jones]]

[[Figuras e animações]]

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Arquivo:Animacao2.gif

2023-04-18T23:37:55Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Animacao2.gif"

Arquivo:Animacao1.gif

2023-04-18T23:37:08Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Animacao1.gif"

Lennard-Jones e propriedades termodinâmicas

2023-04-18T23:02:37Z

Gustavobopsin: /* Potencial de Lennard Jones */

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|500px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|500px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 144, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.724</math> e consideramos a massa das partículas como 1. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os primeiros 100 instantes de tempo e os últimos 100 instantes de tempo. Na primeira animação o sistema ainda não está em equilíbrio, mas na segunda animação sim.
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Primeiros 100 instantes de tempo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 1000 passos da simulação.|500px]]
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ver que a temperatura do sistema já se equilibrou nos primeiros 1000 passos. Poderíamos, a partir daí, já utilizar a simulação para calcular propriedades termodinâmicas no equilíbrio.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ainda calcular a função distribuição radial <math>g(\boldsymbol{r})</math> para esta simulação.
[[Arquivo:fdr.png|thumb|upright=1.0|center|Função distribuição radial.|500px]]

== Códigos ==
[[Simulação - Lennard-Jones]]

[[Figuras e animações]]

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Arquivo:Energia tempo zoom.png

2023-04-18T22:18:44Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Energia tempo zoom.png"

Arquivo:Temperatura tempo zoom.png

2023-04-18T22:15:48Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Temperatura tempo zoom.png"

Arquivo:Pressao tempo zoom.png

2023-04-18T22:14:36Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Pressao tempo zoom.png"

Arquivo:Pressao tempo zoom.png

2023-04-18T22:14:12Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Pressao tempo zoom.png"

Arquivo:Pressao tempo zoom.png

2023-04-18T22:13:15Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Pressao tempo zoom.png"

Lennard-Jones e propriedades termodinâmicas

2023-04-18T19:29:47Z

Gustavobopsin: /* Simulação */

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|400px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|400px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 144, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.724</math> e consideramos a massa das partículas como 1. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os primeiros 100 instantes de tempo e os últimos 100 instantes de tempo. Na primeira animação o sistema ainda não está em equilíbrio, mas na segunda animação sim.
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Primeiros 100 instantes de tempo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 1000 passos da simulação.|500px]]
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ver que a temperatura do sistema já se equilibrou nos primeiros 1000 passos. Poderíamos, a partir daí, já utilizar a simulação para calcular propriedades termodinâmicas no equilíbrio.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ainda calcular a função distribuição radial <math>g(\boldsymbol{r})</math> para esta simulação.
[[Arquivo:fdr.png|thumb|upright=1.0|center|Função distribuição radial.|500px]]

== Códigos ==
[[Simulação - Lennard-Jones]]

[[Figuras e animações]]

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Lennard-Jones e propriedades termodinâmicas

2023-04-18T15:53:30Z

Gustavobopsin: /* Códigos */

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|400px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|400px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 108, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.724</math> e consideramos a massa das partículas como 1. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os primeiros 100 instantes de tempo e os últimos 100 instantes de tempo. Na primeira animação o sistema ainda não está em equilíbrio, mas na segunda animação sim.
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Primeiros 100 instantes de tempo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 1000 passos da simulação.|500px]]
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ver que a temperatura do sistema já se equilibrou nos primeiros 1000 passos. Poderíamos, a partir daí, já utilizar a simulação para calcular propriedades termodinâmicas no equilíbrio.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ainda calcular a função distribuição radial <math>g(\boldsymbol{r})</math> para esta simulação.
[[Arquivo:fdr.png|thumb|upright=1.0|center|Função distribuição radial.|500px]]

== Códigos ==
[[Simulação - Lennard-Jones]]

[[Figuras e animações]]

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Lennard-Jones e propriedades termodinâmicas

2023-04-18T15:53:09Z

Gustavobopsin: /* Códigos */

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|400px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|400px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 108, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.724</math> e consideramos a massa das partículas como 1. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os primeiros 100 instantes de tempo e os últimos 100 instantes de tempo. Na primeira animação o sistema ainda não está em equilíbrio, mas na segunda animação sim.
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Primeiros 100 instantes de tempo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 1000 passos da simulação.|500px]]
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ver que a temperatura do sistema já se equilibrou nos primeiros 1000 passos. Poderíamos, a partir daí, já utilizar a simulação para calcular propriedades termodinâmicas no equilíbrio.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ainda calcular a função distribuição radial <math>g(\boldsymbol{r})</math> para esta simulação.
[[Arquivo:fdr.png|thumb|upright=1.0|center|Função distribuição radial.|500px]]

== Códigos ==
[[Simulação]]

[[Figuras e animações]]

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Lennard-Jones e propriedades termodinâmicas

2023-04-18T15:53:00Z

Gustavobopsin: /* Códigos */

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|400px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|400px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 108, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.724</math> e consideramos a massa das partículas como 1. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os primeiros 100 instantes de tempo e os últimos 100 instantes de tempo. Na primeira animação o sistema ainda não está em equilíbrio, mas na segunda animação sim.
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Primeiros 100 instantes de tempo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 1000 passos da simulação.|500px]]
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ver que a temperatura do sistema já se equilibrou nos primeiros 1000 passos. Poderíamos, a partir daí, já utilizar a simulação para calcular propriedades termodinâmicas no equilíbrio.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ainda calcular a função distribuição radial <math>g(\boldsymbol{r})</math> para esta simulação.
[[Arquivo:fdr.png|thumb|upright=1.0|center|Função distribuição radial.|500px]]

== Códigos ==
[[Simulação]]
[[Figuras e animações]]

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Lennard-Jones e propriedades termodinâmicas

2023-04-18T15:52:41Z

Gustavobopsin:

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|400px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|400px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 108, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.724</math> e consideramos a massa das partículas como 1. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os primeiros 100 instantes de tempo e os últimos 100 instantes de tempo. Na primeira animação o sistema ainda não está em equilíbrio, mas na segunda animação sim.
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Primeiros 100 instantes de tempo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 1000 passos da simulação.|500px]]
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ver que a temperatura do sistema já se equilibrou nos primeiros 1000 passos. Poderíamos, a partir daí, já utilizar a simulação para calcular propriedades termodinâmicas no equilíbrio.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ainda calcular a função distribuição radial <math>g(\boldsymbol{r})</math> para esta simulação.
[[Arquivo:fdr.png|thumb|upright=1.0|center|Função distribuição radial.|500px]]

== Códigos ==
===[[Simulação]]===
===[[Figuras e animações]]===

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Lennard-Jones e propriedades termodinâmicas

2023-04-18T15:52:02Z

Gustavobopsin:

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|400px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|400px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 108, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.724</math> e consideramos a massa das partículas como 1. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os primeiros 100 instantes de tempo e os últimos 100 instantes de tempo. Na primeira animação o sistema ainda não está em equilíbrio, mas na segunda animação sim.
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Primeiros 100 instantes de tempo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 1000 passos da simulação.|500px]]
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ver que a temperatura do sistema já se equilibrou nos primeiros 1000 passos. Poderíamos, a partir daí, já utilizar a simulação para calcular propriedades termodinâmicas no equilíbrio.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ainda calcular a função distribuição radial <math>g(\boldsymbol{r})</math> para esta simulação.
[[Arquivo:fdr.png|thumb|upright=1.0|center|Função distribuição radial.|500px]]

== Códigos ==
=[[Simulação]]=
=[[Figuras e animações]]=

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Lennard-Jones e propriedades termodinâmicas

2023-04-18T15:51:38Z

Gustavobopsin: /* Códigos */

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|400px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|400px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 108, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.724</math> e consideramos a massa das partículas como 1. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os primeiros 100 instantes de tempo e os últimos 100 instantes de tempo. Na primeira animação o sistema ainda não está em equilíbrio, mas na segunda animação sim.
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Primeiros 100 instantes de tempo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 1000 passos da simulação.|500px]]
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ver que a temperatura do sistema já se equilibrou nos primeiros 1000 passos. Poderíamos, a partir daí, já utilizar a simulação para calcular propriedades termodinâmicas no equilíbrio.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ainda calcular a função distribuição radial <math>g(\boldsymbol{r})</math> para esta simulação.
[[Arquivo:fdr.png|thumb|upright=1.0|center|Função distribuição radial.|500px]]

== Códigos ==
==[[Simulação]]==
==[[Figuras e animações]]==

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Lennard-Jones e propriedades termodinâmicas

2023-04-18T15:50:37Z

Gustavobopsin:

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|400px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|400px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 108, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.724</math> e consideramos a massa das partículas como 1. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os primeiros 100 instantes de tempo e os últimos 100 instantes de tempo. Na primeira animação o sistema ainda não está em equilíbrio, mas na segunda animação sim.
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Primeiros 100 instantes de tempo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 1000 passos da simulação.|500px]]
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ver que a temperatura do sistema já se equilibrou nos primeiros 1000 passos. Poderíamos, a partir daí, já utilizar a simulação para calcular propriedades termodinâmicas no equilíbrio.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ainda calcular a função distribuição radial <math>g(\boldsymbol{r})</math> para esta simulação.
[[Arquivo:fdr.png|thumb|upright=1.0|center|Função distribuição radial.|500px]]

== Códigos ==

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Arquivo:Energia tempo.png

2023-04-18T15:45:05Z

Gustavobopsin: Gustavobopsin enviou uma nova versão de "Arquivo:Energia tempo.png"

Lennard-Jones e propriedades termodinâmicas

2023-04-18T15:34:16Z

Gustavobopsin: /* Simulação */

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|400px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|400px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 108, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.724</math> e consideramos a massa das partículas como 1. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os primeiros 100 instantes de tempo e os últimos 100 instantes de tempo. Na primeira animação o sistema ainda não está em equilíbrio, mas na segunda animação sim.
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Primeiros 100 instantes de tempo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 1000 passos da simulação.|500px]]
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ver que a temperatura do sistema já se equilibrou nos primeiros 1000 passos. Poderíamos, a partir daí, já utilizar a simulação para calcular propriedades termodinâmicas no equilíbrio.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ainda calcular a função distribuição radial <math>g(\boldsymbol{r})</math> para esta simulação.
[[Arquivo:fdr.png|thumb|upright=1.0|center|Função distribuição radial.|500px]]

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Arquivo:Fdr.png

2023-04-18T15:33:46Z

Gustavobopsin:

Lennard-Jones e propriedades termodinâmicas

2023-04-18T15:33:28Z

Gustavobopsin: /* Simulação */

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|400px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|400px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 108, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.724</math> e consideramos a massa das partículas como 1. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os primeiros 100 instantes de tempo e os últimos 100 instantes de tempo. Na primeira animação o sistema ainda não está em equilíbrio, mas na segunda animação sim.
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Primeiros 100 instantes de tempo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 1000 passos da simulação.|500px]]
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ver que a temperatura do sistema já se equilibrou nos primeiros 1000 passos. Poderíamos, a partir daí, já utilizar a simulação para calcular propriedades termodinâmicas no equilíbrio.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ainda calcular a função distribuição radial <math>g(\boldsymbol{r})</math> para esta simulação.
|[[Arquivo:fdr.png|thumb|upright=1.0|center|Função distribuição radial.|500px]]

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Lennard-Jones e propriedades termodinâmicas

2023-04-18T15:25:06Z

Gustavobopsin: /* Simulação */

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|400px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|400px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 108, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.724</math> e consideramos a massa das partículas como 1. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os primeiros 100 instantes de tempo e os últimos 100 instantes de tempo. Na primeira animação o sistema ainda não está em equilíbrio, mas na segunda animação sim.
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Primeiros 100 instantes de tempo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 1000 passos da simulação.|500px]]
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 1000 passos da simulação.|500px]]
|}
Podemos ver que a temperatura do sistema já se equilibrou nos primeiros 1000 passos. Poderíamos, a partir daí, já utilizar a simulação para calcular propriedades termodinâmicas no equilíbrio.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|500px]]
|}

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Lennard-Jones e propriedades termodinâmicas

2023-04-18T15:24:43Z

Gustavobopsin: /* Simulação */

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|400px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|400px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 108, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.724</math> e consideramos a massa das partículas como 1. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os primeiros 100 instantes de tempo e os últimos 100 instantes de tempo. Na primeira animação o sistema ainda não está em equilíbrio, mas na segunda animação sim.
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Primeiros 100 instantes de tempo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|500px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 1000 passos da simulação.|500px]]
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|400px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 1000 passos da simulação.|400px]]
|}
Podemos ver que a temperatura do sistema já se equilibrou nos primeiros 1000 passos. Poderíamos, a partir daí, já utilizar a simulação para calcular propriedades termodinâmicas no equilíbrio.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|400px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|400px]]
|}

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Lennard-Jones e propriedades termodinâmicas

2023-04-18T15:23:07Z

Gustavobopsin: /* Simulação */

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|400px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|400px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 108, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.724</math> e consideramos a massa das partículas como 1. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os primeiros 100 instantes de tempo e os últimos 100 instantes de tempo. Na primeira animação o sistema ainda não está em equilíbrio, mas na segunda animação sim.
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Primeiros 100 instantes de tempo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo.|340px]]
|}
A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|400px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 1000 passos da simulação.|400px]]
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|400px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 1000 passos da simulação.|400px]]
|}
Podemos ver que a temperatura do sistema já se equilibrou nos primeiros 1000 passos. Poderíamos, a partir daí, já utilizar a simulação para calcular propriedades termodinâmicas no equilíbrio.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|400px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|400px]]
|}

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Lennard-Jones e propriedades termodinâmicas

2023-04-18T15:20:42Z

Gustavobopsin:

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

== Potencial de Lennard Jones ==

O potencial de Lennard Jones representa a interação entre os àtomos ou moléculas. Possui dois termos distintos: um repulsivo devido às forças de troca do princípio de exclusão de Pauli (origem quântica) e um atrativo devido às forças de van der Walls.

<center><math> V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]. </math></center>

O termo proporcional à <math>\frac{1}{r^{6}}</math> é o termo atrativo, enquanto o termo proporcional à <math>\frac{1}{r^{12}}</math> é o termo repulsivo, muito mais intenso à curtas distâncias. A constante <math>\epsilon</math> descreve a profundidade do poço deste potencial, enquanto <math>\sigma</math> descreve a posição do ponto de equilíbrio estável. Nas imagens abaixo, podemos ver o comportamento do potencial para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.
{|
|[[Arquivo:figura_pot2.png|thumb|upright=0.0|center|Potencial de Lennard-Jones.|400px]]
|[[Arquivo:figura_pot.gif|thumb|upright=0.0|center|Comportamento do potencial de Lennard-Jones para diferentes valores de <math>\epsilon</math> e <math>\sigma</math>.|400px]]
|}
A força que uma partícula sofre, quando sujeita a este potencial é dada por
<center><math>\boldsymbol{F}(\boldsymbol{r}) = -\nabla V(\boldsymbol{r}),</math></center>
<center><math>\boldsymbol{F}(\boldsymbol{r}) = 4\epsilon\left[12\frac{\sigma^{12}}{r^{13}} - 6\frac{\sigma^6}{r^7}\right]\hat{\boldsymbol{r}}.</math></center>

== Dinâmica Molecular ==

A abordagem da dinâmica molecular se baseia em integrar as equações de movimento newtonianas para todas as partículas do sistema. Abaixo, podemos sintetizar um resumo dos passos para realizar uma simulação de dinâmica molecular.

'''1.''' Inicializar a posição e velocidade das partículas do sistema.

'''2.''' Calcular as forças de interação entre todas as partículas do sistema.

'''3.''' Integrar as equações de movimento utilizando algum algoritmo integrador.

'''4.''' Medir grandezas físicas, como energia cinética, energia potencial, temperatura e pressão.

A seguir, vamos explicitar cada um dos passos.

=== Inicialização ===
O primeiro passo é inicializar a posição inicial das partículas. Podemos iniciar a simulação com as partículas distribuídas em qualquer forma que desejarmos, entretanto a forma mais comum de arranjar as partículas é em um formato cúbico ou em um formato triangular. O formato triangular é mais denso do que o formato cúbico. Além disso, a densidade calculada em uma simulação é a densidade numérica, ou seja, quanto de área, ou volume se for em três dimensões, da caixa está ocupada. Esta densidade é calculada como a soma da área ocupada por cada partícula dividida pela area total da caixa. É importante mencionar que as partículas são pontuais, entretanto como o Potencial de Lennard-Jones tende à infinito quanto a distância entre as partículas diminui, cada partícula é como se tivesse um raio efetivo, na qual outra partícula não consegue atravessar. É como se partículas fossem esferas que sofrem colisões umas com as outras.
Além disso, podemos inicializar a velocidade das partículas de diferentes formas.

'''a)''' Velocidade nula:

Se dermos velocidade nula para todas as partículas, a simulação vai demorar mais, uma vez que elas começarão o movimento somente devido à distância de interação entre elas. Se elas estiverem muito afastadas, a força entre elas será muito pequena e demorará muito para que o sistema comece a evoluir e equilibrar, que é o nosso objetivo.

'''b)''' Velocidade aleatória:

A opção mais utilizada. Ao iniciarmos as velocidades de forma aleatória, as partículas não demoram tanto para interagir e conseguimos identificar a parte não-equilibrada e a equilibrada do sistema.

'''c)''' Distribuição de Maxwell-Boltzmann de velocidades:

Já esperamos que quando o sistema for equilibrar, a distribuição de velocidades das partículas seja a distribuição de Maxwell-Boltzman. Iniciar as velocidades desta forma acelera a simulação de forma significativa, entretanto se você tem algum interesse no sistema antes dele entrar em equilíbrio, você perderá informações sobre.

Conforme veremos a seguir, gostaríamos que o nosso sistema tivesse uma determinada temperatura e como ela é calculada utilizando a energia cinética, podemos reescalar as velocidades para obtermos a temperatura desejada. Isto não é um problema uma vez que a temperaturá vai variar durante a simulação.

=== Cálculo das forças ===
Esta etapa é a mais custosa do ponto de vista computacional, uma vez que temos que calcular a força sofrida por cada partícula do sistema devido à todas as outras N-1 partículas, e com isso precisamos computar N(N-1) forças. Entretanto, se a interação exercida em i por j for a mesma que a exercida em j por i, então este número cai pela metade. Ainda é possível reduzir o tempo e o gasto computacional se considerarmos apenas as partículas dentro de um "raio crítico" <math>r_c</math>. Este raio crítico é geralmente tomado como metade do tamanho da caixa. Em nossa simulação, contabilizamos todas as interações pois o número de partículas não é muito grande.

=== Integração das equações de movimento ===
Uma vez que nós obtemos a força sofrida por cada partícula, estamos preparados para integrar as equações de movimento. Para isso, precisamos escolher algum integrador que atenda à nossas necessidades. O integrador mais adequado aos nossos propósitos é um algoritmo do tipo Verlet. Em particular, utilizamos em nossa simulação o algoritmo de Velocity-Verlet.

<center><math>\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)=\boldsymbol{v}(t)+\frac{1}{2}\boldsymbol{f}(t)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right)=\boldsymbol{r}(t)+\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\Delta t</math></center>

<center><math>\boldsymbol{r}\left(t+\Delta t\right) \to \boldsymbol{f}(t+\Delta t)</math></center>

<center><math>\boldsymbol{v}\left(t+\Delta t\right)=\boldsymbol{v}\left(t+\frac{\Delta t}{2}\right)\frac{1}{2} \boldsymbol{f}\left(t+\Delta t\right)\Delta t</math></center>

As vantagens desse algoritmo são que ele preserva a área do espaço de fase, é reversível temporalmente, como são as euqações de Newton, conserva a energia para alcances curtos, é de fácil implementação e é rápido. A desvantagem é que ele não é um dos métodos mais precisos. Devemos nos lembrar que a precisão está diretamente relacionada ao tempo e gasto computacional, além é claro, da escolha adequada do algoritmo.

=== Medidas ===
Nesta etapa, calculamos as grandezas termodinâmicas de interesse. A primeira grandeza relevante é a energia mecânica do sistema, uma vez que gostaríamos de saber se ela está sendo conservada. Para isso, precisamos calcular a energia cinética e a energia potencial do sistema. Naturalmente, a energia mecânica deste sistema será
<center><math>E = \frac{1}{2}m(v_x^2+v_y^2) + \sum_{i<J} 4\epsilon \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^6 \right].</math></center>

Outra grandeza importante é a temperatura. Podemos calcular a temperatura através da energia cinética utilizando o teorema da equipartição da energia da seguinte forma
<center><math>\left\langle \frac{1}{2}mv_\alpha^2\right\rangle=k_BT,</math></center>
onde estamos considerando nosso sistema com 2 dimensões. Como a energia cinética flutua durante a simulação, a temperatura também o fará.A temperatura será
<center><math>T(t)=\sum_{i=1}^N\frac{1}{2}\frac{m_iv_i^2(t)}{k_BN_f},</math></center>
onde <math>N_f=2N-2</math> é o número de graus de liberdade para um número fixo de partículas com momento total fixo.

Por último, gostaríamos de medir a pressão do nosso sistema, que será dada por
<center><math>p=\rho k_BT+\frac{1}{dA}\left\langle \sum_{i<j}\boldsymbol{f}(\boldsymbol{r}_{ij})\cdot\boldsymbol{r}_{ij} \right\rangle,</math></center>
onde <math>\rho</math> é a densidade numérica, <math>k_B</math> é a constante de Boltzmann, <math>T</math> é a temperatura, <math>d</math> é a dimensionalidade do sistema (em nosso caso 2) e <math>A</math> é a aŕea da caixa. Devemos nos lembrar que estamos considerando em nossa simulação <math>k_B</math> como sendo 1, pois como ela é uma constante muito pequena, poderia acumular erros numéricos. Setamos <math>k_B</math> como 1 e todas as unidades das grandezas ficam em termos de <math>k_B</math>.

Por último, devemos nos lembrar que todas estas grandezas estão flutuando no tempo. Para obter, por exemplo, a temperatura média do sistema, deveríamos fazer uma média sobre todas essas flutuações.

=== Convenção da imagem mímina e PBC ===
[[Arquivo:imagem_minima.png|thumb|upright=1.1|right|Convenção da imagem mínima.|300px]]
Existem algumas características que devem ser levadas em consideração quando fazemos uma simulação de dinâmica molecular. Como integrar as equações de movimento requer um poder computacional elevado para muitas partículas, fazemos uma simulação de um "pedaço" do sistema de interesse. Construímos este "pedaço" como uma caixa de lado L e confinamos nossas partículas neste espaço. O número de partículas permanecesse constante e quando uma partícula atravessa a fronteira da caixa, ela retorna do outro lado. Isso simula um sistema infinito. Esta condição de contorno é chamada de condição de contorno periódica (Periodic Boundary Condition - PBC).

Outra característica importante é a convenção da imagem mínima. Esta característica diz respeito à interação entre as partículas. Como cada partícula interage com todas as outras, se a distância entre as partículas for maior do que a metade da caixa, ou seja, L/2, então devemos imaginar que ela vai interagir com a partícula correspondente na caixa "ao lado". Isso é uma extensão do conceito de sistema infinito discutido anteriormente. Ao lado temos uma figura que esquematiza este conceito.

=== Função Distribuição Radial ===
Outra grandeza importante é a função distribuição radial (FDR), que nos diz como é a organização local ao redor de um átomo. A FDR é um histograma e pode ser calculada como
<center><math>g(\boldsymbol{r})=\frac{2A}{N^2}\sum_{i<j}\delta(\boldsymbol{r}-\boldsymbol{r}_{ij}),</math></center>
onde A é a área da caixa e N é o número de partículas.

== Simulação ==

Nesta simulação, utilizamos L = 25 e N = 108, onde N é o número de partículas. Neste contexto, obtemos <math>\rho = 0.724</math> e consideramos a massa das partículas como 1. Inicializamos as partículas em uma rede quadrada e utilizamos um tempo igual à 10000 para que o sistema equilibrasse. Abaixo, mostramos a configuração inicial, os primeiros 100 instantes de tempo e os últimos 100 instantes de tempo. Na primeira animação o sistema ainda não está em equilíbrio, mas na segunda animação sim.
{|
|[[Arquivo:inicializacao.png|thumb|upright=1.0|center|Inicialização do sistema.|340px]]
|[[Arquivo:animacao1.gif|thumb|upright=1.0|center|Primeiros 100 instantes de tempo.|340px]]
|[[Arquivo:animacao2.gif|thumb|upright=1.0|center|Últimos 100 instantes de tempo.|340px]]
|}
Podemos ver que nos instantes iniciais, as partículas tem mais velocidade

A primeira grandeza importante que devemos calcular utilizando nossa simulação é a energia mecânica e verificar se ela está sendo conservada. Para isso, vamos precisar também calcular a energia cinética e a energia potencial do sistema.
{|
|[[Arquivo:energia_tempo.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica durante todos os 10000 passos da simulação.|400px]]
|[[Arquivo:energia_tempo_zoom.png|thumb|upright=1.0|center|Evolução da energia cinética, potencial e mecânica nos primeiros 1000 passos da simulação.|400px]]
|}
A seguir mostramos como a temperatura, calculada na seção 2.4, evolui com tempo.
{|
|[[Arquivo:temperatura_tempo.png|thumb|upright=1.0|center|Evolução da temperatura durante todos os 10000 passos da simulação.|400px]]
|[[Arquivo:temperatura_tempo_zoom.png|thumb|upright=1.0|center|Evolução da temperatura nos primeiros 1000 passos da simulação.|400px]]
|}
Podemos ver que a temperatura do sistema já se equilibrou nos primeiros 1000 passos. Poderíamos, a partir daí, já utilizar a simulação para calcular propriedades termodinâmicas no equilíbrio.

Da mesma forma, podemos calcular a pressão do sistema utilizando a expressão da seção 2.4. Esta expressão depende da temperatura, e como a temperatura flutua no tempo, então a pressão também irá.
{|
|[[Arquivo:pressao_tempo.png|thumb|upright=1.0|center|Evolução da pressão durante todos os 10000 passos da simulação.|400px]]
|[[Arquivo:pressao_tempo_zoom.png|thumb|upright=1.0|center|Evolução da pressão nos primeiros 1000 passos da simulação.|400px]]
|}

== Bibliografia==

D. Frenkel, B. Smit. '''Understanding Molecular Simulation''': From Algorithms to Applications. Academic Press. 2002.

Arquivo:Energia tempo zoom.png

2023-04-18T15:19:17Z

Gustavobopsin:

Arquivo:Energia tempo.png

2023-04-18T15:18:50Z

Gustavobopsin: