6.5 Integrando com Velocity Verlet

Após o cálculo das forças, podemos calcular as componentes da aceleração de cada átomos e após calcular as novas posições e velocidades.

Para calcular as componentes da aceleração utilizaremos a equação 6.5.1.

$$a_{xi} = \frac{F_{xi}}{m_i};a_{yi} = \frac{F_{yi}}{m_i};a_{zi} = \frac{F_{zi}}{m_i}$$ (6.5.1)

!...
! VARIAVEIS QUE DEVEM ESTAR NO MODULE
real(kind=8), allocatable, dimension(:) :: ax, ay, az
integer :: natom
! VARIAVEIS LOCAIS
integer :: i
!...
do i = 1, natom
  ax(i) = fx(i)/m(i)
  ay(i) = fy(i)/m(i)
  az(i) = fz(i)/m(i)
end do
!...

Com o procedimento realizado até o momento, temos as posições ($x_0$, $y_0$, $z_0$), as velocidades ($vx_0$, $vy_0$, $vz_0$), as acelerações ($ax_0$, $ay_0$, $az_0$) e as forças ($fx_0$, $fy_0$, $fz_0$), que compõe o conjunto inicial. Para obter o conjunto das posições ($x_1$, $y_1$, $z_1$), das velocidades ($vx_1$, $vy_1$, $vz_1$), das acelerações ($ax_1$, $ay_1$, $az_1$) e das forças ($fx_1$, $fy_1$, $fz_1$) devemos utilizar algum integrador das equaçõe de movimento de Newton, e a escolha aqui é o método de Velocity Verlet pela simplicidade e facilidade de implementação. Assim poderemos evoluir o sistema até o conjunto das posições ($x_n$, $y_n$, $z_n$), das velocidades ($vx_n$, $vy_n$, $vz_n$), das acelerações ($ax_n$, $ay_n$, $az_n$) e das forças ($fx_n$, $fy_n$, $fz_n$).

De acordo com as equações 5.3.14, novamente exposta abaixo, é necessário primeiro calcular as posições para o tempo ( $t+\Delta t$), na sequência calcular as componentes das forças para o tempo ( $t+\Delta t$), para termos então as componentes das velocidades e das acelerações no tempo ( $t+\Delta t$).

$$\begin{split} \vec{r}(t+\Delta t)& =\vec{r}(t)+\vec{v}(t)\Del... ...}(t+\Delta t/2) + (1/2)\vec{a}(t+\Delta t) \Delta t \end{split}$$

!...
! VARIAVEIAS QUE DEVEM ESTAR NO MODULE
real(kind=8), allocatable, dimension(:) :: x,y,z
real(kind=8), allocatable, dimension(:) :: vx,vy,vz
real(kind=8), allocatable, dimension(:) :: ax,ay,az
integer :: dt1
integer :: natom

! VARIAVEIS LOCAIS
real(kind=8), allocatable, dimension(:) :: vnx1,vny1,vnz1
real(kind=8) :: vxcm,vycm,vzcm
integer :: i

do i = 1, natom
  ! POSICAO PARA T+DELTAT
  ! NESTE MOMENTO PERDE-SE AS POSICOES NO TEMPO T
  x(i) = x(i) + vx(i)*dt1 + 0.5d0*ax(i)*dt1*dt1
  y(i) = y(i) + vy(i)*dt1 + 0.5d0*ay(i)*dt1*dt1
  z(i) = z(i) + vz(i)*dt1 + 0.5d0*az(i)*dt1*dt1
  ! VELOCIDADE PARA T+DELTAT/2
  vnx1(i) = vx(i) + 0.5d0*ax(i)*dt1
  vny1(i) = vy(i) + 0.5d0*ay(i)*dt1
  vnz1(i) = vz(i) + 0.5d0*az(i)*dt1
end do

! CALCULA FORCA BASEADA NAS NOVAS POSICOES
! NESTE MOMENTO PERDE-SE AS FORCAS NO TEMPO T
call forca

! CALCULA ACELERACAO BASEADA NAS NOVAS FORCAS
! NESTE MOMENTO PERDE-SE AS ACELERACOES NO TEMPO T
call aceleracao

! CALCULANDO NOVAS VELOCIDADES
! NESTE MOMENTO PERDE-SE AS VELOCIDADES NO TEMPO T
do i = 1, natom
  vx(i) = vnx1(i) + 0.5d0*ax(i)*dt1
  vy(i) = vny1(i) + 0.5d0*ay(i)*dt1
  vz(i) = vnz1(i) + 0.5d0*az(i)*dt1
  ! CALCULANDO AS COMPONENTES DA VEL DO CENTRO DE MASSA
  vxcm = vxcm + vx(i)
  vycm = vycm + vy(i)
  vzcm = vzcm + vz(i)
end do

! ARTIFICIO UTILIZADO PARA ZERAR O MOMENTO LINEAR
do i=1, natom
  vx(i) = vx(i) - (vxcm/natom)
  vy(i) = vy(i) - (vycm/natom)
  vz(i) = vz(i) - (vzcm/natom)
end do
!...