6.2 Início do Programa de DM

Em princípio, os programas para realizar algum tipo de simulação necessita de dados de entrada, alguns somente de controle, outros somente coordenadas iniciais ou alguns com informações combinadas. É importante pensar em um programa que tenha um caráter genérico de entrada de dados, dentro dessa generalidade é importante que o programa leia os dados de entrada de um arquivo. O porque desse esquema se baseia no fato de que as simulações podem levar um longo tempo para serem finalizadas ou então curtos intervalos de tempo, mas várias simulações com vários intervalos de tempos. No ambiente unix-like (linux) temos uma ferramenta importante e poderosa que nos permite automatizar uma série de tarefas a serem realizadas (programa a serem rodados), mas para isso precisamos de um programa (feito em fortran ou em c/c++) que permite essa automatização.

No caso específico de DM necessitamos das coordendas iniciais e também de um arquivo de controle. Para leitura do arquivo de coordendas utilizaremos o formato xyz que é um formato utilizado por vários programas que também realizam DM e que visualizam a trajetória (uma sequência de arquivos xyz) como os programas xmakemol [22], vmd [23], gopenmol [24], etc. O arquivo de controle da DM faremos de acordo com as nossas necessidades.

O arquivo xyz conforme o programa anterior que gera as coordenadas é composta por três partes. A primeira é o número total de átomos do arquivo que está localizado na primeira linha do arquivo, a segunda é uma linha de comentário localizado na segunda linha do arquivo e a terceira e última é a sequência de posições atômicas em que cada linha possui quatro colunas contendo o rótulo do átomo e as coordenadas x, y, e z. Um exemplo pode ser visto abaixo:

 310
linha de comentario
Ar   7.167486  26.499595  28.546332
Ar  12.382171  16.242033   5.338966
Ar  16.117528  13.491108  28.578685
Ar   5.317546  22.036176  24.344641
Ar   5.236997  27.417759  20.470945

Para ler o arquivo de coordendas, é recomendado que utilizemos o recurso de alocação dinâmica de memória, assim não precisamos compilar o programa todas as vezes que trocarmos a quantidade de átomos que desejarmos estudar. Parte do programa poderia ser feito da seguinte forma:

! IMPORTANTE COLOCAR ESSAS VARIAVEIS DENTRO DE UM MODULE
integer :: natom
character(40) :: comentario
character(2), allocatable, dimension(:) :: rotulo
real(kind=8), allocatable, dimension(:) :: x, y, z, m
!
integer :: i
!...
open(unit=20,file='coord.xyz,status='old',action='read')
!...
! LEITURA DO ARQUIVO DAS COORDENADAS
read(20,*)natom
read(20,*)comentario
allocate(x(natom))
allocate(y(natom))
allocate(z(natom))
allocate(rotulo(natom))
!
! Massa Ar (kg) = 39.948 uma x 1.6605402d-27
!
do i = 1, natom
  read(20,*)rotulo(i),x(i),y(i),z(i)
  x(i) = x(i) * 1.d-10 ! CONVERTENDO PARA METROS
  y(i) = y(i) * 1.d-10 ! CONVERTENDO PARA METROS
  z(i) = z(i) * 1.d-10 ! CONVERTENDO PARA METROS
  m(i) = 6.633526d-27 ! PESO EM KG
end do
!...

O programa será simples e os controles tambéms serão simples, mas com uma base bastante sólida para um programa mais complexo e destinado à pesquisa científica. Teremos como dados de controle as dimensões da caixa (em angstrons), temperatura desejada (em Kelvin), número de passos total (número inteiro), intervalo de tempo de integração (real), o intervalo que serão escritos os dados de energia e a trajetória (inteiro) e o intervalo em que serão reescaladas as velocidades (inteiro).

O arquivo de controle poderia ser da seguinte maneira:

30.d0 30.d0 30.d0
300.d0
1000
1.d0
10
20
13

[1] Dimensoes da caixa lx, ly, lz
[2] Temperatura em Kelvin
[3] Numero de passos da simulacao
[4] dt - intervalo de tempo de integracao em femto segundos
[5] Intervalo para escrita dos dados e trajetoria
[6] Intervalo para reescalar as velocidades
[7] Semente para distribuicao das velocidades

Parte do código para leitura desses dados poderia ser da seguinte maneira:

! IMPORTANTE COLOCAR ESSAS VARIAVEIS DENTRO DE UM MODULE
! DIMENSOES DA CAIXA
real(kind=8) :: lx, ly, lz
! TEMPERATURA DESEJADA
real(kind=8) :: temperatura
! NUMERO DE PASSOS TOTAL
integer :: ntotal
! INTERVALO DE TEMPO PARA INTEGRACAO DAS EQ. DE MOVIMENTO
real(kind=8) :: dt1
! INTERVALO DE TMEPO PARA ESCRITA DOS DADOS E BANHO TERMICO
integer :: nw, nb
! SEMENTE PARA GERACAO DE NUMEROS ALEATORIOS
integer :: isemente
!...
open(unit=21,file='control.inp,status='old',action='read')
!...
read(21,*) lx, ly, lz
lx = lx * 1.d-10 ! CONVETENDO PARA METROS
ly = ly * 1.d-10 ! CONVETENDO PARA METROS
lz = lz * 1.d-10 ! CONVETENDO PARA METROS
read(21,*) temperatura
read(21,*) ntotal
read(21,*) dt1
dt1 = dt1 * 1.d-15 ! CONVETENDO PARA METROS
read(21,*) nw
read(21,*) nb
read(21,*) isemente
!...

Nestes primeiros passos da construção do programa de DM, vimos como gerar e ler os arquivos contendo as coordendas e também como ler o arquivo de controle. Na figura 6.1 é apresentado um esquema básico do funcionamento do programa de DM. É praticamente impossível desenvolver algo desta dimensão se não tivermos pelo menos uma idéia de como irá funcionar o programa. É muito interessante que haja um questionamento para cada pedaço do fluxograma e como poderia ser transformar em um código computacional.

Figura 6.1: Fluxograma básico do programa de Dinêmica Molecular.

Com base no fluxograma a idéia do programa principal de MD poderia assumir a forma abaixo, contendo várias subrotinas, utilizando o compartilhamento de variáveis module e a alocação dinâmica de memória:

!...
! LENDO CONTRLE DA SIMULACAO
call r_control
! LENDO COORDENADAS INICIAIS
call r_xyz
! INICIANDO AS VELOCIDADES
call inicia_vel
! CALCULANDO AS FORCAS
call forca
! CALCULANDO A ACELERACAO
call acel
! CALCULANDO ENERGIA POTENCIAL
call e_potencial
! CALCULANDO ENERGIA CINETICA
call e_cinetica
! CALCULANDO TEMPERATURA
call temper
! ESCREVENDO POSICOES INICIAIS NA TRAJETORIA
call w_trj
! ESCREVENDO ENERGIA, TEMPERATURA, PASSO, ETC
call w_data
!
banho_termico = 0
escreve = 0
!
do i = 2, ntotal
    ! CALCULA NOVA POSICAO
    call evolucao
    ! CALCULA FORCA, ACELERACAO, ENERGIAS, TEMPERATURA
    call forca
    call aceleracao
    call e_potencial
    call e_cinetica
    call temper
    
    ! CONDICAO DE ESCRITA
    if ( escreve == nw ) then
      call w_trj
      escreve = 0
    else
      escreve = escreve + 1
    end if

    ! CONDICAO DO BANHO TERMICO
    if ( banho_termico == nb ) then
      call bt
      banho_termico = 0
    else
      banho_termico = banho_termico + 1
    end if
end do