6.4 Cálculo da Força e Energia Potencial

Como vimos no início do capítulo calcularemos a força conforme a equação 6.0.2. O potencial é o conhecido Lennard-Jones 6-12 dado pela equação 6.4.1 e a derivada pela equação 6.4.2. Em geral, como a força e a energia potencial envolve alguns parametros iguais o cálculo é realizado simultaneamente.

$$U(r) = 4\epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^{6} \right]$$ (6.4.1)

$$F(r) = - \frac{dU(r)}{dr} = 24 \frac{\sigma}{\epsilon} \left[ 2 \left( \frac{\sigma}{r} \right)^{13} - \left( \frac{\sigma}{r} \right)^{7} \right]$$ (6.4.2)

• $r$ é a distância entre o átomo $i$ e o átomo $j$;
• $\epsilon$ é o mínimo valor do gráfico $U(r)$ vs $r$;
• $\sigma$ é a distância de equilíbrio de vdW no gráfico $U(r)$ vs $r$;
• $\epsilon = 0.185$ e $\sigma = 3.868$, da referência [18].

O cálculo da energia potencial é realizada da seguinte maneira:

$$u(r) = \sum \limits_{i}^{N} \sum \limits_{i<j}^{N} U(r_{ij})$$ (6.4.3)

o mesmo deve aser aplicado para o cálculo da força. É de se observar que o somatório envolve o átomo $i$ com todos os átomos $j$, no entando não iremos calcular para todos os átomos $j$. Iremos truncar o somatório no cálculo do potencial e da força se a distância $r_{ij} > r_c$, em que $r_c$ é um raio de corte. A justificativa para essa truncagem baseia-se no fato de que as contribuições acima do raio de corte são muito pequenas, já que o potencial é de curto alcance. Utilizaremos como raio de corte um valor de $2.5\sigma$ [19].

Outro detalhe que devemos nos atentar é que as posições, velocidades e acelerações estão em coordenadas cartesianas e o potencial e a força estão em coordendas polares esféricas ($F$ e $U$ estão em função de $r$). O fato é que para atribuir os valores para as componentes da aceleração deveremos fazer uma conversão de coordenadas polares esféricas para coordenadas cartesianas. Neste ponto é o único momento em que iremos utilizar a transformação de coordenadas.

!...
! VARIAVEIS QUE DEVEM ESTA NO MODULE
! E ALOCADAS PREVIAMENTE
real(kind=8), parameter :: pi = 3.1415926535897967d0
real(kind=8), parameter :: sigma = 3.868d-10
real(kind=8), parameter :: epsilon = 0.185d0
real(kind=8), allocatable, dimension(:) :: x, y, z
real(kind=8), allocatable, dimension(:) :: fx, fy, fz, f, u
real(kind=8) :: eu_total, rc
integer :: natom
! VARIAVEIS LOCAIS
real(kind=8) :: sigma, epsilon, phi, theta
real(kind=8) :: xx, yy, zz, rij, rs, r6, r7, r12, r13
integer :: i, j 
!...
rc = sigma * 2.5d0
eu_total = 0.d0
f(1) = 0.d0; u(1) = 0.d0
!
loopi: do i = 1, natom

  loopj: do j = 1, natom
  
    ! SOMENTE PARA J DIFERENTE DE I
    condji: if ( j /= i ) then
      xx = x(i) - x(j)
      yy = y(i) - y(j)
      zz = z(i) - z(j)
      rij = dsqrt(xx*xx + yy*yy + zz*zz)
      rs  = dsqrt(xx*xx + yy*yy)
      
      ! CONDICAO DO RAIO DE CORTE
      condrij: if ( rij <= rc ) then
      
	! CALCULO DA ENERGIA POTENCIAL DO ATOMO I
	r12 = (sigma/rij)**12
	r6  = (sigma/rij)**6
	u(i) = u(i) + 4.d0 * epsilon * (r12 - r6)
	
	! CALCULO DA FORCA NO ATOMO I
	r13 = (sigma/rij)**13
	r7  = (sigma/rij)**7
	f(i) = f(i) + 24.d0 * (epsilon/sigma) * (2.d0*r13 - r7)

	! CONVERSAO DE COORDENADAS
! ESTA CONVERSAO PODE SER APLICADA PARA QUALQUER CASO ONDE
! AS COORDENADAS ESTIVEREM ESPALHADAS EM TODOS OS QUADRANTES
! DAS COORDENADAS CARTESIANAS. SEM SOMBRA DE DUVIDA QUE E MAIS
! INTELIGENTE/RAPIDO, SE ANTES DE INICIAR A SIMULACAO DESLOCAR
! O SISTEMA PARA O PRIMEIRO QUADRANTE, ASSIM DIMUNUI O NUMERO
! DE CONDICOES ABAIXO.

	! ANGULO PHI (COM RELACAO A Z)
	if ( zz == 0.0d0 ) then
	  phi = pi / 2.0d0
	else
	  phi = dacos(zz / rij))
	end if

	! ANGULO THETA (COM RELACAO A X)
	if ((xx >= 0.0d0) .and. (yy > 0.0d0)) then
	  theta = dasin(yy/rs)
	end if
	if ((xx <  0.0d0) .and. (yy > 0.0d0)) then
	  theta = pi - dasin(yy/rs)
	end if
	if ((xx <  0.0d0) .and. (yy < 0.0d0)) then
	  theta = pi + dasin(abs(yy)/rs)
	end if
	if ((xx >  0.0d0) .and. (yy < 0.0d0)) then
	  theta = (3.0d0*pi/2.0d0) + dacos(abs(yy)/rs)
	end if
	if ((yy == 0.0d0) .and. (xx > 0.0d0)) then
	  theta = 0.0d0
	end if
	if ((yy == 0.0d0) .and. (xx < 0.0d0)) then
	  theta = pi
	end if

	! FORCAS CONVERTIDAS PARA COORD CARTESIANAS
	fx(i) = fx(i) + f(i) * dsin(phi) * dcos(theta)
	fy(i) = fy(i) + f(i) * dsin(phi) * dsin(theta)
	fz(i) = fz(i) + f(i) * dcos(phi)

      end if condrij

    end if condij

  end do loopj

  eu_total = eu_total + u(i)

end do loopi
!...