Duvida Monte Carlo-Metropolis C

[jeeh.ted]

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 (sumE0) e outro com a probabilidade de cada energia ocorrer (q,P[q]).

Segue abaixo o código

#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 de energia tá voltando -7.295913 e o da probabilidade 1 e -9

Editado Julho 24, 2015 por jeeh.ted