Para gerar as co rdenadas construiremos um programa que leia do teclado a quantidade de átomos desejados, o tamanho dos vértices da caixas para garantir que todos átomos estarão dentro da caixa. Como comentado anteriormente utilizaremos o gerador de números aleatórios para atribuir as coodenadas iniciais.
No programa a seguir, é utilizado o tamanho do raio de van der Walls (vdw) para
fazer a distribuição aleatória ou ordenada dos átomos de Argônio (Ar) dentro de
uma caixa. Na distribuição aleatória é respeitada a distância entre átomos de Ar
de no mínimo Åpois o raio de vdw para o átomo de Ar é de
Å. Na distribuições de posições ordenadas fixa-se as
coordenadas
e
e distribui-se ao longo da coordenada
. Quando
é
completamente preenchido adiciona-se é acrescentado um delta em
(
Å) e depois distribui-se ao longo de
novamente, primeiro é preenchido as
posições nas coordenadas em
, depois em
e por último em
. O número de
átomos total é baseado nas dimensões da caixa (x,y,z) e no diâmetro do átomo de
Ar ditado pelo raio de vdw. Supondo que as dimensões da caixa seja
,
,
então teremos
,
e
com o
total de
átomos, que é o número máximo que caberia
dentro da caixa.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! VARIAVEL N SEMPRE SERA COMPARTILHADA ENTRE AS SUBROTINAS
! QUE GERAM OS NUMEROS ALEATORIOS.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module aleatorio
implicit none
integer, save :: n=6677
end module aleatorio
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! PROGRAMA QUE GERAM POSICOES BASEADO NO TAMANHO DA CAIXA
! PARA ENTRADA DO PROGRAMA DE DM
! RAIO DE vdW PARA Ar R_vdw(Ar)= 1.88 Angstron
! http://www.webelements.com/periodicity/van_der_waals_radius/
! VOLUME DA ESFERA = (4.d0/3.d0)*pi*r**3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
program gera_xyz
implicit none
! VARIAVEIS DO GERADOR DE NUMEROS ALEATORIOS
real(kind=8) :: nran ! NUMERO ALEATORIO
integer :: isemente ! SEMENTE INICIAl
! VARIAVEIS DA GERACAO DAS POSICOES ATOMICAS
real(kind=8) :: ax,ay,az ! TAMANHO DA CAIXA
real(kind=8) :: rx,ry,rz ! X(I+1)-X(I), Y(I+1)-Y(I), Z(I+1)-Z(I)
real(kind=8) :: rij ! DISTANCIA ENTRE DOIS ATOMOS
real(kind=8) :: rmin ! DISTANCIA MINIMA ENTRE ATOMOS
real(kind=8), allocatable, dimension(:) :: x,y,z ! POSICAO DOS ATOMOS
real(kind=8), parameter :: vdw_ar = 1.88d0
integer :: i, j, k, l ! UTILIZADOS NOS DO'S
integer :: natom,natomx,natomy,natomz ! NUMERO DE ATOMOS
integer :: refaca ! PARA NO SORTEIO DE POSICOES
! VERIFICANDO QUANTIDADE DE ATOMOS DENTRO DA CAIXA
real(kind=8) :: diametro
integer :: lx,ly,lz
integer :: natomos
integer :: tipo
! ARQUIVO DE SAIDA
open(unit=20, file='coord.xyz', status='replace', action='write')
! DIMENSOES DA CAIXA
write(*,*) 'Entre com o tamanho de x, y e z (em Angstrons):'
write(*,*) '(Dimensoes da caixa)'
read(*,*) ax, ay, az
! INDICANDO O NUMERO MAXIMO DE ATOMOS QUE CAIBAM DENTRO DA CAIXA
! 1 ATOMO DE Ar OCUPA UM QUADRADO DE LADO 2.d0*1.88d0 Angstrons
diametro = vdw_ar * 2.d0
lx = int(ax/diametro)
ly = int(ay/diametro)
lz = int(az/diametro)
natomos = lx*ly*lz
write(*,*)' '
write(*,*)'O numero maximo de atomos que cabem dentro da caixa e:',natomos
write(*,*)' '
! NUMEROS DE ATOMOS DESEJADO
write(*,*) 'Entre com o numero de atomos (deve ser inteiro):'
write(*,*) 'OBS: Caso seja aleatorio nao colocar mais do que',int(0.9d0*dfloat(natomos)),'atomos (90% do total)'
read(*,*) natom
! CONDICAO VOLUME DO NUMERO DE ATOMOS DESEJADO NAO PODE SER MAIOR
! QUE O VOLUME DA CAIXA DESEJADA
if (natom > natomos) then
write(*,*) 'ATENCAO >>> O Volume de atomos desejado eh muito'
write(*,*) '>>> maior que o volume da caixa desejado !!!'
stop
end if
! CONDICAO ENTRE POSICAO ORDENADA OU ALEATORIA
write(*,*)' '
write(*,*)'Sorteio das posicoes sera aleatorio (1) ou ordenado (2) ?'
read(*,*) tipo
write(*,*)' '
! ALOCANDO AS POSICOES ATOMICAS
allocate(x(natom))
allocate(y(natom))
allocate(z(natom))
! DISTANCIA MINIMA ENTRE ATOMOS
rmin=2.d0*vdw_ar
! ESCREVENDO CABECALHO DO ARQUIVO DE SAIDA
write(20,*) natom
write(20,500) ax,ay,az
500 format('dimensoes da caixa',f10.6,f10.6,f10.6)
aleat: if ( tipo == 1 ) then
write(*,*) 'Entre com a semente: '
read(*,*) isemente
! CHAMANDO A SUBROTINA SEED
call semente(isemente)
loop1: do i = 1, natom
! CHAMADA DO GERADOR DE NUM.ALEAT.
call num_aleatorio(nran)
x(i) = nran * ax
call num_aleatorio(nran)
y(i) = nran * ay
call num_aleatorio(nran)
z(i) = nran * az
! VERIFICANDO POSICAO ATOMICA
j = i
if ( j > 1 ) then
loop2: do
refaca = 0
! CALCULA RIJ
! SE MENOR QUE RMIN SORTEIA NOVA POSICAO
do k=1,j-1
rx = x(j) - x(k)
ry = y(j) - y(k)
rz = z(j) - z(k)
rij = dsqrt( rx*rx + ry*ry + rz*rz )
if (rij < rmin) then
refaca=1
end if
end do
! SORTEIO DA NOVA POSICAO
if ( refaca == 1 ) then
call num_aleatorio(nran)
x(j) = nran * ax
call num_aleatorio(nran)
y(j) = nran * ay
call num_aleatorio(nran)
z(j) = nran * az
end if
! CONDICAO DE SAIDA
if ( refaca == 0 ) exit
end do loop2
end if
write(20,501) x(i),y(i),z(i)
end do loop1
endif aleat
! POSICOES NO MODO ORDENADO
! NESTE TIPO DE SORTEIO EH PURAMENTE PARA FINS DIDATICOS
! NAO SEI SE TEM SENTIDO FISICO
orden: if ( tipo == 2 ) then
natomx = int( ax / rmin )
natomy = int( ay / rmin )
natomz = int( az / rmin )
l = 0
do i = 1, natomz
! CONDICAO DE PARADA: SENAO ERRO CRITICO
if ( l == natom ) exit
do j = 1, natomy
! CONDICAO DE PARADA: SENAO ERRO CRITICO
if ( l == natom ) exit
do k = 1, natomx
l = l + 1
if ( k == 1 ) then
x(l) = 0.5d0 * rmin
else
x(l) = x(1) + dfloat(k-1) * rmin
end if
if ( j == 1 ) then
y(l) = 0.5d0 * rmin
else
y(l) = y(1) + dfloat(j-1) * rmin
end if
if ( i == 1 ) then
z(l) = 0.5d0 * rmin
else
z(l) = z(1) + dfloat(i-1) * rmin
end if
! CONDICAO DE PARADA: SENAO ERRO CRITICO
if ( l == natom ) exit
end do
end do
end do
do i = 1, natom
write(20,501) x(i),y(i),z(i)
end do
end if orden
501 format('Ar',2x,f16.6,f16.6,f16.6)
stop
end program gera_xyz
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine num_aleatorio(nran)
use aleatorio
implicit none
real(kind=8), intent(out) :: nran
n = mod( 8127*n+28417 , 134453 )
nran = dfloat(n) / 134453.d0
end subroutine num_aleatorio
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine semente(isemente)
use aleatorio
implicit none
integer, intent(in) :: isemente
n=abs(isemente)
end subroutine semente