Ir ao conteúdo
  • Cadastre-se

Dúvida Monte Carlo-Metropolis C


Posts recomendados

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 60float 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;} 

O arquivo da energia retorna -7.295913 e da probabilidade 1 e -9

Link para o comentário
Compartilhar em outros sites

Minhas matemáticas vão bem enferrujadas, porém si você explicar em que parte do seu código você acha que esta o erro posso dar uma olhada e tentar achar o erro. Si não teríamos que averiguar a traves do depurador linha a linha para ver o que não esta saindo bem. É pouco porém é o que eu posso oferecer =(

Link para o comentário
Compartilhar em outros sites

Visitante
Este tópico está impedido de receber novas respostas.

Sobre o Clube do Hardware

No ar desde 1996, o Clube do Hardware é uma das maiores, mais antigas e mais respeitadas comunidades sobre tecnologia do Brasil. Leia mais

Direitos autorais

Não permitimos a cópia ou reprodução do conteúdo do nosso site, fórum, newsletters e redes sociais, mesmo citando-se a fonte. Leia mais

×
×
  • Criar novo...