Dúvida Monte Carlo-Metropolis C

jeehted
(usa Outra)

Enviado em 24/07/2015 - 14:37h

Estou fazendo um programa de Monte Carlo-Metropolis com implementação do Potencial de Lennard-Jones, porém a minha energia está se mantendo constante enquanto deveria variar e não consigo achar o problema do código.
O programa se baseia em colocar as partículas em uma caixa de lado 30, fazer movimentações aleatórias nessas particulas e calcular sua energia antes e depois de cada movimentação através do Potencial de Lennard-Jones. Esta nova energia é aceita ou não através do algoritmo de Metropolis, onde aceita se a diferença entre as energias for menor que zero. Se esta é maior, gera-se um número aleatório entre 0 e 1 e compara-se este número com a energia calculada por Metropolis. Se esta for menor que o número aleatório, esta energia é aceita, se não é rejeitada.
No final do programa, ele deve me gerar um arquivo com as energias aceitas e outro com a probabilidade de cada energia ocorrer.

#include<stdio.h>

#include<stdlib.h>

#include<math.h>

#include<time.h>

#define part 60



float randomic()

{

	float f;

	f=((rand()%2001/1000.0)) - 1;

	return f;

	}



int main(){

	

	FILE *arq1;

	FILE *arq2;

	

	float raio, pos[part][3],d,dx,dy,dz,E0,sumE0,e,sigma,K,xalet,yalet,zalet,posant[part][3];

	float dxtot,dytot,dztot,dtot,E,Etot,sumE,sumEtot,dE,B,T,r,aceita,Eeq,P[60],Prob[1000];

	int lado, y,z,i,j,k,l,n,dmax,m,nn,nnn,p,q;

	

	arq1=fopen("energia.txt","w");

	arq2=fopen("probabilidade.txt","w");

		

	raio=3.41;

	lado=30;

	y=0;

	z=0;

	e=119.8;

	sigma=3.41;

	K = 1.3806503*pow(10,-3);

	dmax=5;

	T=300;

	p=0.1;

	nn=0;

	nnn=0;

	sumE0=0;

	sumE=0;

	sumEtot=0;

	E0=0;

	E=0;

	Etot=0;

	pos[k][1]=0;

	pos[k][2]=0;

	pos[k][3]=0;

	posant[k][1]=0;

	posant[k][2]=0;

	posant[k][3]=0;

	

	

	pos[1][1]=3.41; // posição inicial de 1 a partir do centro em x

	pos[1][2]=3.41; // posição inicial de 1 a partir do centro em y

	pos[1][3]=3.41; // posição inicial de 1 a partir do centro em z

	

	//colocando partículas na caixa

	for(i=2;i<=part;i++){

		pos[i][1]=pos[i-1][1] + 1.5*raio;

		pos[i][2]=3.41 + 1.5*raio*y;

		pos[i][3]=3.41 + 1.5*raio*z;

		if(pos[i][1] > lado - raio){

			y++;

			pos[i][1]=3.41;

			pos[i][2]=pos[i-1][2] + 1.5*raio;

					}

		if(pos[i][2] > lado - raio){

			y=3.41;

			z++;

			pos[i][2]=3.41;

			pos[i][3]=pos[i-1][3] + 1.5*raio;

		}

		//printf("%f %f %f  ",pos[i][1],pos[i][2],pos[i][3]);

	}

	//calculando energia inicial

	for(j=1;j<=part;j++){

		for(k=j+1;k<=part;k++){

			dx = pos[k][1] - pos[j][1];

			dy = pos[k][2] - pos[j][2];

			dz = pos[k][3] - pos[j][3];

			dx = dx - round(dx/lado)*lado;

			dy = dy - round(dy/lado)*lado;

			dz = dz - round(dz/lado)*lado;

			//printf("posicao[%d]=%f\n",k,pos[k][1]);

			//printf("posicao[%d]=%f\n",j,pos[j][1]);

			//printf("%f  ",dx);

			//printf("%f",dy);

			//printf("%f",dz);

						

			d = (dx*dx) + (dy*dy) + (dz*dz);

			d = sqrt(d);

			//printf("%f  ",d);

		    						

			E0=4*e*K*(pow((sigma/d),12)-pow((sigma/d),6));

			//printf("%f  ",E0);

			sumE0=sumE0+E0;

		}

		}

		//movimentação das partículas

	for(n=1;n<=50000;n++){

		k= 1 + ((int)rand()%(part+1)); // escolhe partícula aleatória 

		xalet=randomic();

		yalet=randomic();

		zalet=randomic();

		//printf("%f",xalet);

		

		posant[k][1]=pos[k][1];

		posant[k][2]=pos[k][2];

		posant[k][3]=pos[k][3];

		

		pos[k][1] = pos[k][1] + xalet*dmax;

		pos[k][2] = pos[k][2] + yalet*dmax;

		pos[k][3] = pos[k][3] + zalet*dmax;

		

		for(l=1;l<=3;l++){

			if(pos[k][l] > lado - raio){

				pos[k][l] = pos[k][l] + raio - (lado-raio)*(int)(pos[k][l]/(lado-raio));

			}

			if(pos[k][l] < raio){

				pos[k][l] = pos[k][l] - raio + (lado-raio)*(int)(pos[k][l]/(lado-raio));

		}

		}

	//energia depois do movimento

	for(m=1;m<=part;m++){

		if(m != k){

			

			dx = posant[k][1] - pos[m][1];

			dy = posant[k][1] - pos[m][2];

			dz = posant[k][3] - pos[m][3];

			dx = dx - round(dx/lado)*lado;

			dy = dy - round(dy/lado)*lado;

			dz = dz - round(dz/lado)*lado;

			

			d = (dx*dx) + (dy*dy) + (dz*dz);

			d = sqrt(d);

			

			E = 4*e*K*(pow((sigma/d),12) - pow((sigma/d),6));

			sumE=sumE+E;

			

			dxtot = pos[k][1] - pos[m][1];

			dytot = pos[k][2] - pos[m][2];

			dztot = pos[k][3] - pos[m][3];

			dxtot = dxtot - round(dx/lado)*lado;

			dytot = dytot - round(dy/lado)*lado;

			dztot = dztot - round(dz/lado)*lado;

			

			dtot = (dxtot*dxtot) + (dytot*dytot) + (dztot*dztot);

			dtot = sqrt(dtot);

			

			Etot = 4*e*K*(pow((sigma/dtot),12) - pow((sigma/dtot),6));

			sumEtot = sumEtot + Etot;

		} // if k

	} // if m

	//diferença das energias

	

	dE = sumEtot - sumE;

	printf("ok");

	if (dE < 0){ //aceita

		sumE0 = sumE0 + dE;

	

		}

	if(dE > 0){ // rejeita

		r = rand()%1001/1000;

		B = exp((-1.0)*(dE)/(K*T));

		if(B < r){ // rejeita

			pos[k][1] = posant[k][1];

			pos[k][2] = posant[k][2];

			pos[k][3] = posant[k][3];

		}

		if(B > r){ //aceita

			sumE0 = sumE0 + dE;

			aceita++;

		}

	}

	printf("ok1");

	// equilíbrio

	if (n > 35000){

		Eeq = (sumE0*sumE0);

		Eeq = sqrt(sumE0);

		//Eeq = Eeq/p; // energia por porção

		//P[(int)Eeq] = P[(int)Eeq] + 1;

		

		nn++;

		

		if(nn == 1){

			fprintf(arq1,"%f \n",sumE0);

			nnn++;

			if(nnn == 100){

				printf("%d %f \n",m,sumE0);

				nnn++;

				nnn=0;

			}

			nn=0;

		}

		printf("ok2");

		//probabilidades

	for(q=1;q<=1000;q++){

		Prob[q] = Prob[q]/(n-1000);

	fprintf(arq2,"%d  %f",q,Prob[q]);

	}

	

		} // n

	

}

fclose(arq1);

fclose(arq2);

return 0;

} 

Meu arquivo da energia tá voltando -7.295913 e o da probabilidade 1 e -9