DM: um primeiro programa

[EM CONSTRUÇÃO] A seguir 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, podemos organizar o programa da seguinte maneira, onde cada etapa será explicada detalhadamente.

Aqui declaramos quais são 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 e condições iniciais utilizadas na 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 Lx 16
#define Ly 16
#define dt 0.01
<\pre>

<math> N_x \mbox{e} N_y <\math> representam, respectivamente, o número de partículas na direção <math> x \mbox{e} y <\math>. 

/***********      declaring functions       ************/
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 *ffx, double *ffy, double *xx, double *yy,  double *xxp, double *yyp, double eenergia, double *energy_pp, double *vxm ) ;


void main()
{
  /***********  declaring variables    ************/
  // declara as variaveis necessarias
  // ex: double *fx ;

  /**********   allocating memory     ************/
  // aloca a memoria
  // ex. de alocacao DINAMICA de memoria:
  // fx  = (double *)malloc( NP * sizeof(double) );  
  
  initialize_positions( px, py); 
  initialize_velocities( vx, vy);
  for(i = 0; i < NP; i++)
    { ppx[i] = px[i] - vx[i]*deltat;   ppy[i] = py[i] - vy[i]*deltat; } // positions in the initial time 
  
  tt=0;  iframe=0;
  // faz um loop no tempo 
  //       chama a funcao que calcula forcas
  //       chama a funcao que integra as Eqs de movimento
  //	   incrementa o tempo de simulacao
  //       a cada intervalo de tempo definido, fazer medidas


}


void initialize_positions( double *xx, double *yy )
{
  int      i ;
  double   deltax, deltay; 
  
  deltax = (double)Lx/Nx;   deltay = (double)Ly/Ny;
  for(i = 0; i < NP; i++)
    {
      xx[i] = (i%Nx)*deltax ; 
      yy[i] = (i/Ny)*deltay ;
      // printf("%lf %lf   %d \n", xx[i], yy[i], i);
    }
}

void initialize_velocities(double *vx, double *vy )
{ 
  int      i ;
  double   vcmx, vcmy, v2, fs;

  vcmx=0.0;  vcmy=0.0;   v2=0.0;
  for(i = 0; i < NP; i++)
    {
      vx[i] = (double)rand()/RAND_MAX ;
      vy[i] = (double)rand()/RAND_MAX ;
      vcmx += vx[i];    vcmy += vy[i];
      v2 += vx[i]*vx[i] + vy[i]*vy[i] ;  // kinetic energy
    }
  vcmx /= NP ;  vcmy /= NP ;  v2 /= NP;
  fs=sqrt(2*TEMP/v2) ;  // fator de escala -- escolhe uma dada temp T e computa isto : SE escolher uma dist maxwelliana, isto nao eh necessario (me parece)
  
  // set the velocity center of mass to zero (linear momentum of the system is zero)
  // and set its value in such a way that \sum m_i v^2 = Kb T for each degree of freedom
  for(i = 0; i < NP; i++)
    {  vx[i] = (vx[i] - vcmx)*fs ;         vy[i] = (vy[i]-vcmy)*fs ;  } 
 }

// fx = -dV/dx = -(x/r) (dV/dr)
// fx= 48x/r² (1/r^12 - 1/r^6) --  LJ in reduced Units
void compute_forces( double *ffx, double *ffy, double *xx, double *yy, double *eenergia ) 
{
  int     i, j;
  double  x1, y1, x2, y2,dxx, dyy, r2, r2i, r6, ff, energia=0.0 ;
  
  for( i=0; i < NP; i++ )
    {  ffx[i] = 0.0;    ffy[i] = 0.0; }
	    
  for( i=0; i < NP; i++ )
    {
      x1  = xx[i];    y1 = yy[i]; 
     for( j=i+1; j < NP; j++ )
	{
	  x2  = xx[j];    y2 = yy[j];
	  
	  dxx = x1 - x2;          // compute the distance between particles taking into account the periodic boundary conditions   
	  dxx = fmod(dxx, Lx);  
	  dxx=dxx-rint(dxx/Lx)*Lx ;   //note: rint(arg) : rounds its argument to the nearest whole number
		  
	  dyy = y1 - y2;           
	  dyy = fmod(dyy, Ly);
	  dyy=dyy-rint(dyy/Ly)*Ly ;
	 	  
	  r2 = dxx*dxx + dyy*dyy ;
	  if(r2<RC2) // rc2: critical value or r
	    {
	      r2i = 1/r2 ;         
	      r6 = r2i*r2i*r2i;
	      ff = 48*r2i*r6*(r6-0.5);                         // LJ potential integrated                         
	      ffx[i] += ff* dxx ;   ffy[i] += ff* dyy   ;      // update forces
	      ffx[j] -= ff* dxx ;   ffy[j] -= ff* dyy   ;
	      energia += 4*r6*(r6-1) ;                         // compute energy
	    }
	}
    }
  *eenergia=energia;
}


void integrate_EqMotion( double *ffx, double *ffy, double *xx, double *yy,  double *xxp, double *yyp, double eenergia, double *energy_pp, double *vxm ) 
{
  int     i; 
  double  xp, yp, vxi, vyi,  vcmx, vcmy, v2, temp  ;

  vcmx=0.0;  vcmy=0.0;   v2=0.0;
  for(i = 0; i < NP; i++)
    {
      xp = 2*xx[i] - xxp[i] + deltat2 * ffx[i] ;
      yp = 2*yy[i] - yyp[i] + deltat2 * ffy[i] ;
      vxi = (xp - xx[i]) / (2*deltat);
      vyi = (yp - yy[i]) / (2*deltat);
      
      vcmx += vxi;    vcmy += vyi;         // to monitor the average velocity
      v2 += vxi*vxi + vyi*vyi ;            // kinetic energy
      xxp[i] = xx[i] ;  yyp[i] = yy[i] ;   // update the "previous position"
      yyp[i] = yy[i] ; 
      xx[i] = xp;                          // update the positions in a current time
      yy[i] = yp;
   } 
  temp=v2/(2*NP) ;                         // kbT = < m*v^2 > per degree of freedom
  *energy_pp=(eenergia + 0.5*v2)/NP;        // energia total = ( U + K ) -- energy per particle =energia/NP
  *vxm = vcmx /NP ;
}