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:
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.
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 MAINMa 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:
