6.3 Atribuição das Velocidades

Inicialmente não possuimos velocidades para os átomos e precisamos de alguma forma atribuir respeitando algumas condições iniciais. Uma delas o tipo de distribuição (distribuição tipo Boltzman) e a outra que o momento linear seja zero. Utilizaremos o gerador de números aleatórios e então atribuiremos valores entre $-1$ e $+1$, em seguida zeramos o momento linear.

! SOMENTE NESTA SUBROTINA QUE SERA UTILIZADO 
! A ROTINA DE NUMEROS RANDOMICOS
! VARIAVEIS RECOMENDADO ESTAR NO MODULE
real(kind=8), allocatable, dimension(:) :: vx, vy, vz
real(kind=8) :: svx, svy, svz
! VARIAVEL LOCAL
integer :: i
! SORTEANDO VELOCIDADES INICIAIS, VALORES NO INTERVALOR [-1,1]
svx = 0.d0; svy = 0.d0; svz = 0.d0
do i = 1, natom
  call num_aleatorio(nran)
  vx(i) = (-1.d0)**nint(nran) * nran
  svx = svx + vx(i)
  call num_aleatorio(nran)
  vy(i) = (-1.d0)**nint(nran) * nran
  svy = svy + vy(i)
  call num_aleatorio(nran)
  vz(i) = (-1.d0)**nint(nran) * nran  
  svz = svz + vz(i)
end do
! ZERANDO O MOMENTO LINEAR
do i = 1, natom
  vx(i) = vx(i) - (svx/natom)
  vy(i) = vy(i) - (svy/natom)
  vz(i) = vz(i) - (svz/natom)
end do
! ESCALANDO AS VELOCIDADES PARA A TEMPERATURA ALVO
call vel_scal

Para escalar a velocidade para a temperatura alvo, precisamos inicialmente calcular a temperatura atual. Para isso chamamos a rotina que calcula a temperatura atual. O calculo da temperatura instantânea é baseado na equação 6.3.1 e 6.3.2.

$$\sum \limits_{i=1}^{N} \frac{1}{2} m_i \vec{v}_i^2 = \frac{3}{2}N k_b T$$ (6.3.1)

$$T = \frac{2 \sum \limits_{i=1}^{N} \frac{1}{2} m_i \vec{v}_i^2}{3 N k_b}$$ (6.3.2)

Talvez seja desejável criar uma subrotina que calcule a energia cinética, pois iremos utilizar esta informação todas as vezes para calcular a temperatura instantânea e para calcular a energia mecânica total. Rotulamos a energia cinética total como $K$, assim teremos:

$$T = \frac{2 K}{3 N k_b}$$ (6.3.3)

e a velocidade escalada para a temperatura desejada é dada pela equação 6.3.4

$$v_{xi}^{nova} = v_{xi}^{antiga} \sqrt{\frac{T_{alvo}}{T_{inst}}}$$ (6.3.4)

calculando a energia cinética

!..
ek = 0.d0
do i = 1, natom
  v = dsqrt(vx(i)*vx(i) + vy(i)*vy(i) + vz(i)*vz(i))
  ek = ek + m(i)*v
end do
ek = ek*0.5d0
!...

escalando as velocidades

!...
! kb = 1.3800D-23 J/K = 8.6200D-5 eV/K
real(kind=8), parameter :: kb = 1.3800D-23

t_inst = (2.d0 * ek) / (natom * kb * 3.d0)
fator = dsqrt{temperatura/t_inst}
do i = 1, natom
  vx(i) = vx(i) * fator
  vy(i) = vy(i) * fator
  vz(i) = vz(i) * fator
end do
!...