Ir ao conteúdo
  • Cadastre-se

Jéssica Tedesco

Membro Júnior
  • Posts

    1
  • Cadastrado em

  • Última visita

Reputação

0
  1. 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

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