Code dynamique moléculaire et tableaux

Sev07

Membre confirmé
20 Juillet 2011
10
0
Bonsoir,
Je suis étudiant en licence de physique et cette année nous avons appris des bases de C, avec cela nous devons réaliser mon binôme et moi (en suivant un algorithme donné) un programme de dynamique moléculaire pour un nombre d'atome donné dans un temps donné.

Cependant pour un grand nombre d'atomes (1000) et un grand nombre d'itération (10 000) cela est très long, on en est donc à une étape où nous devons créer une sorte de grille dans notre système pour faciliter les calculs, en fait un atome n'interagit qu'avec d'autres atomes dans un rayon donné, donc on se limite a faire les calculs de force avec les atomes dans les cases adjacentes à la sienne.

On doit donc avoir un truc du genre: (ix, iy, n), avec ix l'abscisse de la boite de l'atome x, iy l'ordonné de la boite de l'atome y et n le numéro de l'atome n°i du système dans la boite, sachant que (ix, iy, 0 renvoie le nombre d'atome dans la boite ix,iy.

A un moment de mon programme je dois actualiser des positions et définir les boites, cependant je ne sais pas du tout comment faire pour pouvoir créer un tableau de ce genre la de manière à ce que je puisse faire des opérations du genre:
(ix, iy, 0) ++; (Si un atome est dans la boite ix, iy on dit qu'on a un atome de plus dans cette boite et dans le même temps à un atome du système on lui attribue un numéro d'atome dans la boite)
(ix, iy, (ix,iy,0)) = i; (On repère l'atome i du système par sa boite et son numéro dans la boite).

Je sais créer un tableau comme ça:
type tableau[n];
J'ai essayé un truc du genre type tableau[ix][iy][n]; mais cela ne semble pas fonctionner.

Voici le programme (sans la boite) afin que vous puissiez comprendre mieux:
Bloc de code:
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>

#define nat 1024 //Nombre d'atomes
#define dt 0.005 //Pas de temps
double arrondi(double);

//////////////////////////ARRONDI//////////////////////////

double arrondi(double n) {
    
    if (n > 0) n = floor(n+0.5);
    if (n < 0) n = ceil(n-0.5);
    if (n == 0) n = 0;
    
    return n;
    
}

//////////////////////////DECLARATIONS//////////////////////////


FILE *pf1,*pf2;
double rx[nat],ry[nat],rz[nat],vx[nat],vy[nat],vz[nat]; //Positions et vitesses
double ax[nat],ay[nat],az[nat],axa[nat],aya[nat],aza[nat]; //Accelerations et accelerations au temps precedent
double fx[nat],fy[nat],fz[nat]; //Forces
double dimx=40.0, dimy=40.0, dimz=12.0; //Dimensions du syteme suivant x et y pour conditions aux limites periodiques

double dx,dy,dz,r2,ep,ec,alpha,fij,t,rxi,ryi,rzi,etot,temp; 
//Variation de position, position, energies potentielle et cinetiques, 
//coefficient de rescaling de la vitesse, forces, temps, positions, energie totale, temperature
int i,j,npas,ifixtemp,maxpas,nouta,noutb; //Variables, pas, temperature fixee ou non, maximum de pas, pas d'ecriture
double v2; //Vitesse au carre


//////////////////////////CONSTANTES ET PARAMETRES//////////////////////////

//COEFFICIENTS
double c1=2.0*dt*dt/3.0,    c2=-1*dt*dt/6.0,    c3=dt/3.0;
double c4=5.0*dt/6.0,   c5=-1*dt/6.0;
double ap=4.0,  bp=4.0,   af=24.0,bf=48.0;
double PC=2.5; //rayon de coupure de potentiel


//DEBUT DU MAIN
int main()
{
	
    //TEMPS EXECUTION
    double temps;
    clock_t t1, t2;
    
    t1 = clock();
    
    
	//LECTURE FICHIER CONTROLE
	pf1=fopen("/Users/Start/Desktop/CTRL.txt","r");
    
	fscanf(pf1,"%d %d %d \n",&maxpas,&nouta,&noutb);
	fscanf(pf1,"%d %lf \n",&ifixtemp,&temp);
    
	fclose(pf1);
	
	//LECTURE FICHIER CONFIGURATION INITIALE ET ECRITURE VISU.IN
	pf1=fopen("/Users/Start/Desktop/CONF1024.txt","r");
	pf2=fopen("/Users/Start/Desktop/VISUIN.txt","w");
    
	for (i=0;i<nat;i++) 
	{ 
		fscanf(pf1,"%lf    %lf    %lf    \n",&rx[i],&ry[i],&rz[i]);
		fscanf(pf1,"%lf    %lf    %lf    \n",&vx[i],&vy[i],&vz[i]);
		fscanf(pf1,"%lf    %lf    %lf    \n",&ax[i],&ay[i],&az[i]);
		fscanf(pf1,"%lf    %lf    %lf    \n",&axa[i],&aya[i],&aza[i]);
        
		fprintf(pf2,"%lf    %lf    %lf    \n",rx[i],ry[i],rz[i]);
	}
    
	fclose(pf1);
	fclose(pf2);
	
	
	
	//BOUCLE PRINCIPALE
	pf1=fopen("/Users/Start/Desktop/ENE.txt","w");
	pf2=fopen("/Users/Start/Desktop/TRAJ.txt","w");
	
	for (npas=1;npas<=maxpas;npas++)
	{
		
		
		//CALCUL DES NOUVELLES POSITIONS
		for (i=0;i<nat;i++)
        {   
			rx[i] = rx[i] + dt*vx[i] + c1*ax[i] + c2*axa[i];
			rx[i]=rx[i]-arrondi(rx[i]/dimx)*dimx; //Conditions aux limites periodiques (CLP)
            
			
			ry[i] = ry[i] + dt*vy[i] + c1*ay[i] + c2*aya[i];
			ry[i]=ry[i]-arrondi(ry[i]/dimy)*dimy; //CLP
			
			rz[i] = rz[i] + dt*vz[i] + c1*az[i] + c2*aza[i];
			rz[i]=rz[i]-dimz*arrondi(rz[i]/dimz); //CLP
			
			//Mise a zero des forces	 
			fx[i]=0.0;
			fy[i]=0.0;
			fz[i]=0.0;
		}
		
		
		//CALCUL DES NOUVELLES FORCES
		ep=0.0;
        
		for (i=0;i<nat-1;i++)
        {
			for (j=i+1;j<nat;j++)
            {
                
				//Distance entres les atomes
				dx=rx[i]-rx[j];
				dx=dx-dimx*arrondi(dx/dimx); //CLP
				
                dy=ry[i]-ry[j];
				dy=dy-dimy*arrondi(dy/dimy); //CLP
				
				dz=rz[i]-rz[j];
				dz=dz-dimz*arrondi(dz/dimz); //CLP
                
				//Calcul des forces (interactions entre les atomes)
				r2=dx*dx+dy*dy+dz*dz;
				
                
                //Rayon de coupure
				if(sqrt(r2) < PC) fij=bf/(r2*r2*r2*r2*r2*r2*r2)  -  af/(r2*r2*r2*r2);
				else fij = 0;
                
				
				fx[i] += fij*dx;
				fx[j] += -1*fij*dx;
				
				fy[i] += fij*dy;
				fy[j] += -1*fij*dy;
				
				fz[i] += fij*dz;
				fz[j] += -1*fij*dz;
				
				//Calcul energie potentielle
				ep += bp/(r2*r2*r2*r2*r2*r2) - ap/(r2*r2*r2);        	       
			}       
		}
		ep=ep/nat; //Calcul energie potentielle par atome
		
		
		//CALCUL DES NOUVELLES VITESSES ET CORRECTION SI ON FIXE LA TEMPERATURE
		
		ec=0.0; //Reinitialisation de l'energie cinetique pour eviter un cumul des valeurs
		for (i=0;i<nat;i++)
        {
			vx[i]=vx[i]+c3*fx[i]+c4*ax[i]+c5*axa[i];
			vy[i]=vy[i]+c3*fy[i]+c4*ay[i]+c5*aya[i];
			vz[i]=vz[i]+c3*fz[i]+c4*az[i]+c5*aza[i];
			
			//Calcul energie cinetique
			ec += vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i];
		}
        
		ec=(ec*0.5)/nat; //Calcul energie cinetique par atome
        alpha=sqrt(1.5*temp/ec); //Calcul du coefficient de rescaling
        
		
		//Rescaling de la vitesse
		if (ifixtemp == 1) 
        {
            
            for (i=0;i<nat;i++)
            {   
				vx[i]=alpha*vx[i];
				vy[i]=alpha*vy[i];
				vz[i]=alpha*vz[i];
			}
		}
		
        
        
		//ACTUALISATION DES ACCELERATIONS
		
		for (i=0;i<nat;i++)
        {
			axa[i]=ax[i];
			aya[i]=ay[i];
			aza[i]=az[i];
			ax[i]=fx[i];
			ay[i]=fy[i];
			az[i]=fz[i];
		}
		
		//Calcul du temps ecoule 
		t=(float)npas*dt;
		
		//SORTIES
        if (npas%nouta == 0)  
        {
			printf("%d / %d  \n",npas,maxpas);
            fprintf(pf1,"%lf %lf %lf \n",t,ep,ep+ec);
        }
		
        if (npas%noutb == 0)  
        {
            for (i=0;i<nat;i++) 
            {
				fprintf(pf2,"%lf %lf %lf\n",rx[i],ry[i],rz[i]);
			}
        }
		
        
	}//FIN BOUCLE PRINCIPALE
	
	//FERMETURE DES FICHIERS OUVERTS
	fclose(pf1);
	fclose(pf2);
	
	//ECRITURE FICHIER CONFIGURATION FINALE ET ECRITURE VISU.OUT
	
	pf1=fopen("/Users/Start/Desktop/CONFOUT.txt","w");
	pf2=fopen("/Users/Start/Desktop/VISUOUT.txt","w");
	for (i=0;i<nat;i++) 
	{
		fprintf(pf1,"%lf %lf %lf\n",rx[i],ry[i],rz[i]);
		fprintf(pf1,"%lf %lf %lf\n",vx[i],vy[i],vz[i]);
		fprintf(pf1,"%lf %lf %lf\n",ax[i],ay[i],az[i]);
		fprintf(pf1,"%lf %lf %lf\n",axa[i],aya[i],aza[i]);
		fprintf(pf2,"%lf    %lf    %lf    \n",rx[i],ry[i],rz[i]);
	}
	fclose(pf1);
	fclose(pf2);
	
    //TEMPS EXECUTION
	t2 = clock();
    temps = (double)(t2-t1)/CLOCKS_PER_SEC;
    printf("Il a fallu %f s pour executer ce programme, soit %lf min\n", temps, temps/60);
	
}
//FIN DU MAIN

Ma question est donc: Comment créer un tel tableau ?

Je précise que mon système est un cristal, une plaque d'atome espacé de 1.2 de longueur 14.4 suivant X et Y et de longueur 6 suivant Z (on ne s'occupe donc pas de faire une grille 3D mais juste une 2D pour X et Y).
J'ai essayé de commenter un maximum et j'espère être clair, je précise que j'utilise Xcode et que je suis sous 10.7.2. J'aurais voulu mettre le code en couleur pour plus de clarté mais je ne sais pas comment faire :confuses: .

Dans l'incrémentation des positions en fait j'ai essayé de mettre, en ayant déclaré int grd[1][1][1];:
ix = (int)(rx/2); (boites de largeur 2)
ix = (int)(rx/2);
grd[ix][iy][0]++;
grd[ix][iy][grd[ix][iy][0]] = i;

Il y a compilation mais lors de l'exécution j'ai un: EXC_BAD_ACCES, cela signifie quoi ?


Merci d'avance pour votre aide et la peine que vous prenez à me lire. :)
 
Dernière édition:

ntx

Membre vénérable
Club MacG
15 Octobre 2004
12 130
377
92
En C un tableau est uniquement défini par l'adresse de son élément "0" et la taille de cet élément. L'allocation mémoire des éléments du tableau se fait de manière consécutive : n éléments occuperont n cases mémoires consécutives, toutes de tailles identique.

Quand tu essaies d'atteindre l'élément x, le programme se contente d'incrémenter l'adresse de base de x cases de la taille de l'élément "0" pour accéder à ce nouvel élément. Pour les adresse, cela donne quelque chose du genre pour un tableau de type "elem" :

adresse de tab[x] = adresse de tab[0] + x * sizeof(elem)

Le programme ne contrôle pas du tout ce que tu fais. Donc si tu lui donnes un mauvais indice d'élément, il va aller taper dans une case mémoire de ton programme qui n'a rien à voir avec ton tableau : en langage de développeur, tu jardines :D:D:D Pour l'OS, EXC_BAD_ACCES. :zen:
 

Sev07

Membre confirmé
20 Juillet 2011
10
0
Merci pour ton explication ntx.
Cependant, comment donc créer un tel tableau ?
Tel qu'il contiennent trois variables ix et iy (qui varient entre 0 et 19) et n qui varie entre 0 et 199 , de plus (par exemple) (1,3,5) et (6,14,9) doivent renvoyer des valeurs différentes et bien précises correspondant au numéro de l'atome du système qui est dans la boite (ix,iy) et qui est l'atome n de cette boite (par exemple l'atome 3 de la boite peut être l'atome 984 du système) ? Si je fais tableau[3] = {ix, iy, n} je ne peux pas manipuler les trois variables à la fois, et si je fais tableau[1][1][1] ça me renvoie un EXC_BAD_ACCES.
 

ntx

Membre vénérable
Club MacG
15 Octobre 2004
12 130
377
92
Un tableau à deux dimensions se déclare de la façon suivante :

int tab[n][m];

Il est considéré par le compilo comme un tableau de tableau. Et là tu pourras accéder à un élément de ton tableau avec un appel du style tab[x][y].