DM: um primeiro programa

[EM CONSTRUÇÃO] Aqui mostraremos o passo a passo de como montar o primeiro programa para simulação de dinâmica molecular. A linguagem de programação utilizada será o C. Como um primeiro esboço, o programa pode ser organizado como mostra o fluxograma abaixo. Posteriormente, introduziremos partes do código comentando cada etapa. Para o leitor mais familiarizado com DM, ao final da página há o programa completo sem comentários.

A ideia básica da simulação pode ser entendida como: Leitura dos parâmetros da simulação, inicialização das posições das partículas, inicialização das velocidades das partículas, laço temporal onde se calcula a força, integra-se as equações de movimento e, de tempos em tempos, medidas são realizadas.

Abaixo introduzimos partes do código onde começamos declarando as bibliotecas standard utilizadas do C

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

Está parte é reservada para declarações de constantes que definem os parâmetros da simulação. As constantes podem ser definidas utilizando-se das diretivas de compilação do C através do comando "#define" ou podem ser lidas a partir de um arquivo externo chamado pelo programa. A vantagem do segundo método é que o programa só precisará ser compilado uma única vez para diferentes condições iniciais. Porém, por simplicidade e clareza, optaremos pelo primeiro caso.

#define Nx 16
#define Ny 16
#define Np (Nx*Ny)
#define Lx 16
#define Ly 16

#define Rc 4
#define Rc2 (Rc*Rc)
#define T 1.

#define dt 0.01
#define transient 5.
#define total_time 10.
#define measure 0.1
#define samples ((total_time-transient)/measure)

$N_{x}\;{\mbox{e}}\;N_{y}$ representam, respectivamente, o número de partículas nas direções $x\;{\mbox{e}}\;y$.

$L_{x}\;{\mbox{e}}\;L_{y}$ representam, respectivamente, o tamanho da rede nas direções $x\;{\mbox{e}}\;y$.

$Np=N_{x}N_{y}$ número total de partículas

$dt$ define o intervalo de tempo de integração na simulação.

$R_{c}$ é o raio de corte acima do qual o cálculo da força é desconsiderado.

$total\_time$ é o tempo absoluto da simulação.

$transient$ é o tempo que esperamos até o sistema supostamente equilibrar e, então, começarmos a fazer medidas.

$measure$ define o intervalo de tempo entre as medições.

$samples$ nos dá o número total de medidas que serão realizadas

Conforme o programa é construído, novas constantes podem eventualmente aparecer à medida que novas funções são implementadas. Abaixo declaramos as funções essenciais ao programa, onde cada uma delas será detalhadamente explicada no decorrer da implementação.

void initialize_positions(double *xx, double *yy ) ;
void initialize_velocities(double *vx, double *vy ) ;
void compute_forces( double *ffx, double *ffy, double *xx, double *yy, double *eenergia ) ;
void integrate_EqMotion( double *fx, double *fy, double *xx, double *yy,  double *xxp, double *yyp, double eenergia, double *energy_pp, double *vxm ) ;

A seguir introduzimos a função principal

void main(void)
{

FILE *file1,*file2;
file1=fopen("energy.dat","w");
file2=fopen("trajectories.dat","w");

unsigned i;
double *x,*y,*xo,*yo,*vx,*vy,*Fx,*Fy,Ep,Ec,Vxm=0.,t;

size_t size=N*sizeof(double);

x= (double*)malloc(size);
y= (double*)malloc(size);
xo=(double*)malloc(size);
yo=(double*)malloc(size);
vx=(double*)malloc(size);
vy=(double*)malloc(size);
Fx=(double*)malloc(size);
Fy=(double*)malloc(size);

Acima são declaradas todas as variáveis essenciais ao programa, assim como a alocação de memória para os vetores como posição, velocidade, força, etc. Com os vetores alocados, podemos agora inicializar as velocidades e posições das partículas na rede

initialize_velocities(vx,vy);
initialize_positions(x,y,xo,yo,vx,vy);

Abaixo é mostrada a parte principal do programa, o loop sobre o tempo. A cada passo de tempo, a rotina $compute\_forces$ calcula a força e a energia potencial das partículas a partir de suas posições, na sequência, a rotina $integrate\_EqMotion$ atualiza as posições das partículas. O $for$ em $i$ imprime no arquivo "trajectories.dat" as posições e velocidades de todas as partículas no instante $t$. No $if$ testamos se o tempo total é maior que o tempo de equilíbrio e se a medida deve ser feita, caso seja, no arquivo "energy.dat" imprimiremos, respectivamente, o tempo, energia cinética, energia potencial e a energia total.

for(t=0.;t<total_time;t+=dt)
{
        compute_forces(Fx,Fy,x,y,&Ep);
        integrate_EqMotion(Fx,Fy,x,y,xo,yo,&Ec,&Vxm);

        for(i=0;i<N;i++)fprintf(file2,"%lf %lf %lf %lf\n",x[i],y[i],vx[i],vy[i]);
        fflush(file2);
        rewind(file2);

        if( t>transient && (1.-fmod(t/measure,1.) )<0.00001 )
        {
                fprintf(file1,"%lf %lf %lf %lf\n",t,Ec,Ep,Ec+Ep);
                fflush(file1);
        }
}

fclose(file1);
fclose(file2);

free(x);
free(y);
free(xo);
free(yo);
free(vx);
free(vy);
free(Fx);
free(Fy);
}

Agora descremos detalhadamente todas as rotinas que aparecem na parte principal. Começamos com a inicialização das posições, a rotina recebe os ponteiros das velocidades nas direções $x$ e $y$.

void initialize_velocities(double *vx,double *vy)
{
  unsigned i;
  double Vcx=0.,Vcy=0.,V2=0.,fs;

  for(i=0;i<N;i++)
  {
      vx[i]=1.*rand()/RAND_MAX;
      vy[i]=1.*rand()/RAND_MAX;

      Vcx+=vx[i];
      Vcy+=vy[i];
      V2 +=vx[i]*vx[i] + vy[i]*vy[i];
  }

  Vcx/=N;
  Vcy/=N;
  V2 /=N;
  fs=sqrt(2.*T/V2);

  for(i=0;i<N;i++)
  {
        vx[i]=(vx[i]-Vcx)*fs;
        vy[i]=(vy[i]-Vcy)*fs;
  }

}