hw1.c.tar.gz_Ising model_ising

#include<stdio.h> #include<stdlib.h> #include<math.h> //#include "mysubroutine_u.h" //#include "mysubroutine_u.c" #include "ran2.c" //H ( \sigma )= -J \sum_<i j> \sigma_i \sigma_j, J =-1 double ham(int *sig,int m) { int j1,j2,j1x,j1y,j2x,j2y, n =m*m,flagx,flagy; double h=0.0; for(j1=0;j1<n;j1++){ j1x=j1/m; j1y=j1%m; for(j2=j1+1;j2<n;j2++){ j2x=j2/m; j2y=j2%m; if(j1x==j2x) flagx=0; else if( j2x ==j1x+1 || (j1x==0 && j2x== m-1)) flagx=1; else flagx=-1; if(j1y==j2y) flagy=0; else if( j2y ==j1y+1 || j2y==j1y-1 || (j1y==0 && j2y== m-1) || (j2y==0 && j1y== m-1)) flagy=1; else flagy=-1; if((flagx==0 && flagy==1)||(flagy==0 && flagx==1)) h=h+(-1.0)*sig[j1]*sig[j2]; } } return h; } void compu_nn(int *sig,int *sig1,int *nn,int m) { int j,j1,j2,j3,j4,n=m*m,near_j; for(j=0;j<10;j++) nn[j]=0; for(j=0;j<n;j++){ j1=j-1; if(j1<0) j1=j1+n; j2=(j+1)%n; j3=j-m; if(j3<0) j3=j3+n; j4=(j+m)%n; near_j=(sig[j1]+sig[j2]+sig[j3]+sig[j4])/2; if(sig[j]>0){ sig1[j]=2-near_j; nn[2-near_j]=nn[2-near_j]+1; } else{ sig1[j]=7-near_j; nn[7-near_j]=nn[7-near_j]+1; } } return; } int flip(int *sig1,int m,int cla_i,int r1) { int j,num=-1; for(j=0;j<m*m && num<r1;j++){ if(sig1[j]==cla_i) num++; } return(j-1); } double dham(int *sig,int m,int j) { int n=m*m, j1,j2,j3,j4, s,s1,s2,s3,s4; double dh; s=sig[j]; j1=j-1; if(j1<0) j1=j1+n; j2=(j+1)%n; j3=j-m; if(j3<0) j3=j3+n; j4=(j+m)%n; s1=sig[j1];s2=sig[j2];s3=sig[j3];s4=sig[j4]; dh=-1.0*(s1+s2+s3+s4)*2*s; return dh; } void dnn(int *sig1,int *nn,int m,int j,int cla_i) { int j1,state=0,n=m*m; nn[cla_i]=nn[cla_i]-1; nn[(cla_i+5)%10]=nn[(cla_i+5)%10]+1; sig1[j]=(cla_i+5)%10; if(sig1[j]>4) state=1; j1=j-1; if(j1<0) j1=j1+n; if(state==0){ nn[sig1[j1]]=nn[sig1[j1]]-1; sig1[j1]=sig1[j1]-1; nn[sig1[j1]]=nn[sig1[j1]]+1; } else{ nn[sig1[j1]]=nn[sig1[j1]]-1; sig1[j1]=sig1[j1]+1; nn[sig1[j1]]=nn[sig1[j1]]+1; } j1=(j+1)%n; if(state==0){ nn[sig1[j1]]=nn[sig1[j1]]-1; sig1[j1]=sig1[j1]-1; nn[sig1[j1]]=nn[sig1[j1]]+1; } else{ nn[sig1[j1]]=nn[sig1[j1]]-1; sig1[j1]=sig1[j1]+1; nn[sig1[j1]]=nn[sig1[j1]]+1; } j1=j-m; if(j1<0) j1=j1+n; if(state==0){ nn[sig1[j1]]=nn[sig1[j1]]-1; sig1[j1]=sig1[j1]-1; nn[sig1[j1]]=nn[sig1[j1]]+1; } else{ nn[sig1[j1]]=nn[sig1[j1]]-1; sig1[j1]=sig1[j1]+1; nn[sig1[j1]]=nn[sig1[j1]]+1; } j1=(j+m)%n; if(state==0){ nn[sig1[j1]]=nn[sig1[j1]]-1; sig1[j1]=sig1[j1]-1; nn[sig1[j1]]=nn[sig1[j1]]+1; } else{ nn[sig1[j1]]=nn[sig1[j1]]-1; sig1[j1]=sig1[j1]+1; nn[sig1[j1]]=nn[sig1[j1]]+1; } return; } double phi(int *sig,int m) { double p=0.0; int n=m*m,j; for(j=0;j<n;j++) p=p+sig[j]; return p; } void inisig(int *sig,int m) { int j,n=m*m; for(j=0;j<n;j++) sig[j]=1; } void outphi(double *q,int n) { FILE *fp; int k; if( (fp=fopen("phi.txt","wb"))==NULL) {printf("open file wrong!\n");getchar();exit(0);} for(k=0;k<n;k=5+k){ fprintf(fp,"%e",q[k]); fputc(13,fp);fputc(10,fp); } fclose(fp); } /************************************************************* npre //取样之前的迭代次数为 Npre* Nintv nstep //总取样次数 nintv //每演化Nintv次取样一次 m //M*M grid **************************************************************/ void KMC(FILE *fp,double beta,int npre,int nstep,int nintv,int m) { int n=m*m,*sig,*sig1,k,kk,j,i,cla_i,rr,*nn; double r,h,dh, *P_,*Q,q1_,q2_,q3_,qq,qq1, u,c,ph; long idum; //printf("A"); idum=-10; //random number seed sig = (int*)malloc(sizeof(int)*n); sig1= (int*)malloc(sizeof(int)*n); //printf("C"); //q1=(double*)malloc(nstep*sizeof(double)); //q2=(double*)malloc(nstep*sizeof(double)); //q3=(double*)malloc(nstep*sizeof(double)); //qq=(double*)malloc(nstep*sizeof(double)); P_=(double*)malloc(10*sizeof(double)); Q =(double*)malloc(10*sizeof(double)); nn = (int*)malloc(10*sizeof(int)); //printf("B"); //initial state inisig(sig,m); h=ham(sig,m); P_[0]=exp(-8*beta);P_[1]=exp(-4*beta);P_[2]=1;P_[3]=1;P_[4]=1; P_[5]=1;P_[6]=1;P_[7]=1;P_[8]=exp(-4*beta);P_[9]=exp(-8*beta); //begin sample for(k=0;k<npre;k++){ if(k%10000==1) idum=idum-1; //evlve Nintv times for(kk=0;kk<nintv;kk++){ //sig ' r=ran2(&idum); j=(int)(r*n); sig[j]=-sig[j]; //accept or not? (metro strategy) dh=dham(sig,m,j); if(dh > 0){ r=ran2(&idum); if (r> exp(-beta*dh)) sig[j]=-sig[j];//not accept else h=h+dh; } else h=h+dh; //accept } } compu_nn(sig,sig1,nn,m); q1_=0.0;q2_=0.0;q3_=0.0;qq1=0.0; for(k=0;k<nstep;k++){ //printf("8"); if(k%3000==1) idum=idum-1; for(kk=0;kk<nintv;kk++){ // //printf("0"); Q[0]=nn[0]*P_[0]; for(i=1;i<10;i++) Q[i]=Q[i-1]+nn[i]*P_[i]; //printf("1"); r=ran2(&idum); r=r*Q[9]; for(i=8;i>=0 && r<Q[i];i--) ; cla_i=i+1; //printf("2"); r=ran2(&idum); rr=(int)(r*nn[cla_i]); //printf("3"); j=flip(sig1,m,cla_i,rr); //printf("4"); sig[j]=-sig[j]; dnn(sig1,nn,m,j,cla_i); //printf("5"); h=h+dham(sig,m,j); //if(cla_i<5) h=h+(4-cla_i)*4-8; //else h=h+(cla_i-5)*4-8; r=ran2(&idum); qq=-log(r)/Q[9]; q1_=q1_+h*qq; q2_=q2_+h*h*qq; qq1=qq1+qq; q3_=q3_+qq*phi(sig,m); //printf("%e",qq); //printf("7\n"); // printf("%d\t%e\n",k,h,q3[k-npre]); } } // if(k%100==1) // printf("%d\n",k); //outphi(q3,nstep); //printf("1"); // for(k=0;k<nstep;k++){ // q1_=q1_+q1[k]*qq[k];q2_=q2_+q2[k]*qq[k];q3_=q3_+q3[k]*qq[k]; // } q1_=q1_/qq1; q2_=q2_/qq1; q3_=q3_/qq1;//(double)(nstep); u=1/(double)(n)*q1_; c=beta*beta/(double)(n)*(q2_-q1_*q1_); ph=1/(double)(n)*q3_; // printf("%e",qq1); fprintf(fp,"%e\t%e\t%e\t%e\n",beta,u,c,ph); // fprintf(fp,"\n"); // fputc(13,fp);fputc(10,fp); //printf("3"); free(sig);free(sig1);free(P_);free(Q);free(nn); //printf("4"); } void main() { FILE *fp; double beta; if( (fp=fopen("result1.txt","wb"))==NULL) {printf("open file wrong!\n");getchar();exit(0);} for(beta=0.01;beta<=1;beta=beta+0.01){ printf("\n%e",beta); KMC(fp,beta,10000,30000,1,30); } fclose(fp); }