Prog c

bonjour:
voila audessous le programme qui prend le dynamique d'un reseau d'atome au cours de temps(nouvelles positions,forces...).mais moi je voulais prendre l'exemple d'une poutre = resau d'atomes et encastrée=ecouches où les atomes sont fixes ainsi comment introduire une force appliquee sur l' extrimité.
merci
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define nat 121
#define dt 0.005

// CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
// DECLARATIONS
//

FILE *pf1,*pf2,*pf3,pf4;
double rx[nat],ry[nat],vx[nat],vy[nat];
double ax[nat],ay[nat],axa[nat],aya[nat];//accel avant et apres
double fx[nat],fy[nat];
double dx,dy,r2,ep,ec,alpha,fij,t,rxi,ryi,etot;
int i,j,npas,ifixtemp;
double v2;

// CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
// CONSTANTES ET PARAMETRES
//

double c1=2.0*dt*dt/3.0,c2=-dt*dt/6.0,c3=dt/3.0;
double c4=5.0*dt/6.0,c5=-dt/6.0;

double ap=4.0,bp=4.0,af=24.0,bf=48.0;


int ifixtemp=1;
double temp=0.1;
int maxpas=10000;

main()
{

// CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
// LECTURE FICHIER CONFIGURATION INITIALE
// ET ECRITURE VISU.IN
//
pf1=fopen("CONF.IN","r");
pf2=fopen("VISU.IN","w");
for (i=0;i<nat;i++)
{
fscanf(pf1,"%lf %lf \n",&rx[i],&ry[i]);
fscanf(pf1,"%lf %lf \n",&vx[i],&vy[i]);
fscanf(pf1,"%lf %lf \n",&ax[i],&ay[i]);
fscanf(pf1,"%lf %lf \n",&axa[i],&aya[i]);
fprintf(pf2,"%lf %lf %d \n",rx[i],ry[i],i);
}
close(pf1);
close(pf2);



// CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
// BOUCLE PRINCIPALE
//
pf1=fopen("ENE","w");
pf2=fopen("TRAJ","w");
pf3=fopen("VIT","w");
for (npas=1;npas<=maxpas;npas++)
{


//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
// calcul des nlles positions
//
for (i=0;i<nat;i++){
rx[i] = rx[i] + dt*vx[i] + c1*ax[i] + c2*axa[i];
ry[i] = ry[i] + dt*vy[i] + c1*ay[i] + c2*aya[i];
fx[i]=0.0;
fy[i]=0.0;
}


//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
// calcul des nouvelles forces
//
ep=0.0;
for (i=0;i<nat-1;i++){
for (j=i+1;j<nat;j++){
dx=rx[i]-rx[j];
dy=ry[i]-ry[j];
r2=dx*dx+dy*dy;
fij=bf/(r2*r2*r2*r2*r2*r2*r2) - af/(r2*r2*r2*r2);
fx[i]+=fij*dx;
fy[i]+=fij*dy;
fx[j]+=-fij*dx;
fy[j]+=-fij*dy;
ep+= bp/(r2*r2*r2*r2*r2*r2) -ap/(r2*r2*r2);
}
}
ep=ep/nat;


//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
// calcul des nouvelles vitesses
// + corection pout fixer T
//

ec=0.0;
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];
ec+=vx[i]*vx[i]+vy[i]*vy[i];
}
ec=ec*0.5/nat;
if (ifixtemp==1) {
alpha=sqrt(temp/ec);
for (i=0;i<nat;i++){
vx[i]=alpha*vx[i];
vy[i]=alpha*vy[i];

}
}

//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
// actualisation des tableaux
//
for (i=0;i<nat;i++){
axa[i]=ax[i];
aya[i]=ay[i];
ax[i]=fx[i];
ay[i]=fy[i];
}



//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
// sorties
//
t=(float)npas*dt;
if (fmod(npas,1000)==0) {
printf("%d / %d \n",npas,maxpas);
}

if (fmod(npas,10)==0) {
fprintf(pf1,"%d %lf %lf \n",npas,ep,ep+ec);
}

if (fmod(npas,500)==0) {
for (i=0;i<nat;i++) { fprintf(pf2,"%lf %lf \n",rx[i],ry[i]);
}
}
if (fmod(npas,100)==0) {
for (i=0;i<nat;i++) { fprintf(pf3," %lf \n",sqrt(vx[i]*vx[i]+vy[i]*vy[i]));
}
}


// FIN BOUCLE PRINCIPALE
}
close(pf1);
close(pf2);
close(pf3);

// CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
// ECRITURE FICHIER CONFIGURATION FINALE
// ET ECRITURE VISU.OUT
//
//
pf1=fopen("CONF.OUT","w");
pf2=fopen("VISU.OUT","w");
for (i=0;i<nat;i++)
{
fprintf(pf1,"%lf %lf \n",rx[i],ry[i]);
fprintf(pf1,"%lf %lf \n",vx[i],vy[i]);
fprintf(pf1,"%lf %lf \n",ax[i],ay[i]);
fprintf(pf1,"%lf %lf \n",axa[i],aya[i]);
fprintf(pf2,"%lg %lg %d \n",rx[i],ry[i],i);
}
close(pf1);
close(pf2);


// CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
// FIN DU MAIN
}