#include #include #include #include // Parameters #define N 100 #define nmcs 10000 #define border 0 #define scale 0.4 #define Pi 4.0*atan(1.0) double ransu() { return rand()/(RAND_MAX+1.0); //return genrand_real2();//Mersenne-Twister } void initialize1(double alpha[][N],double s[][N][2],int *nspin,int *P,double theta) { int i,j; double r; for(i = 0; i < N ; i++){ for(j = 0; j < N ; j++){ alpha[i][j] = 0.0; s[i][j][0] = 0.0; s[i][j][1] = 0.0; } } for(i = 0; i < *nspin ; i++){ for(j = 0; j < *nspin ; j++){ r = ransu(); //r = genrand_real2();//Mersenne-Twister alpha[i][j] = ((int)(r*(*P)))*theta; } } } double delta_E(int i,int j,double alpha[][N],double H[],double check,double *J,int *nspin) { int iprev,inext,jprev,jnext; double topS, bottomS,leftS,rightS; double H_part1,H_part2; //境界条件 iprev = (i == 0 ? *nspin-1 : i-1); inext = (i == *nspin-1 ? 0 : i+1); jprev = (j == 0 ? *nspin-1 : j-1); jnext = (j == *nspin-1 ? 0 : j+1); topS = cos(alpha[i][j]-alpha[iprev][j])-cos(check-alpha[iprev][j]); bottomS = cos(alpha[i][j]-alpha[inext][j])-cos(check-alpha[inext][j]); leftS = cos(alpha[i][j]-alpha[i][jprev])-cos(check-alpha[i][jprev]); rightS = cos(alpha[i][j]-alpha[i][jnext])-cos(check-alpha[i][jnext]); H_part1 = H[0]*cos(alpha[i][j])+H[1]*sin(alpha[i][j]); H_part2 = H[0]*cos(check)+H[1]*sin(check); return ((*J)*(leftS + rightS + topS + bottomS)-H_part1+H_part2); } double total_magnetization_x(double alpha[][N],int *nspin) { int i,j; double Mx; Mx = 0.0; for( i = 0 ; i < *nspin ; i++){ for( j = 0 ; j < *nspin ; j++){ Mx += cos(alpha[i][j]); } } return Mx; } double total_magnetization_y(double alpha[][N],int *nspin) { int i,j; double My; My = 0.0; for( i = 0 ; i < *nspin ; i++){ for( j = 0 ; j < *nspin ; j++){ My += sin(alpha[i][j]); } } return My; } double total_energy(double alpha[][N],double H[],double *J,int *nspin) { int i,j; double E; int iprev,inext,jprev,jnext; double topS, bottomS,leftS,rightS; double H_part; E = 0.0; for( i = 0 ; i < *nspin ; i++){ for( j = 0 ; j < *nspin ; j++){ //境界条件 //三項演算子 iprev = (i == 0 ? *nspin-1 : i-1); // if(i==0) {iprev=nspin-1} else { iprev=i-1} inext = (i == *nspin-1 ? 0 : i+1); // jprev = (j == 0 ? *nspin-1 : j-1); jnext = (j == *nspin-1 ? 0 : j+1); topS = cos(alpha[i][j]-alpha[iprev][j]); bottomS = cos(alpha[i][j]-alpha[inext][j]); leftS = cos(alpha[i][j]-alpha[i][jprev]); rightS = cos(alpha[i][j]-alpha[i][jnext]); H_part = H[0]*cos(alpha[i][j])+H[1]*sin(alpha[i][j]); E += -1.0/2.0*(*J)*(leftS + rightS + topS + bottomS)-H_part; } } return E; } int main() { FILE *fp0,*fp1,*fp2; double alpha[N][N],s[N][N][2]; int i,j,istep,mcs; int ispin,jspin; double hanten,dE,check,theta; double H[2],T; double magnetization[2],magnetization_sum[2]; // double magnetization_2[2],magnetization_sum2[2]; double energy,energy_sum; // double energy_2,energy_sum2; // double magnetic_susceptibility[2],specific_heat; int nspin,P; double J,Hx_max,Hy_max,T_max; int start; fp0 = fopen("input_step_xy.dat","r"); fp1 = fopen("xy_data.dat","w"); fp2 = fopen("xy_spin.dat","w"); if(fp0 == NULL || fp1 == NULL || fp2 == NULL ){ printf("failure\n"); return 0; } fscanf(fp0,"%d",&nspin); fscanf(fp0,"%d",&P); fscanf(fp0,"%lf",&J); fscanf(fp0,"%lf",&T_max); fscanf(fp0,"%lf",&Hx_max); fscanf(fp0,"%lf",&Hy_max); fclose(fp0); theta = 2.0*Pi/P; printf(" initialize \n spin = %d*%d P=%d J=%6.3f \n T=%6.3f (Hx,Hy)=(%5.3f,%5.3f) \n input number start=1 or exit=2\n",nspin,nspin,P,J,T_max,Hx_max,Hy_max); scanf("%d",&start); if(start != 1){ printf("exit\n"); return 0; } else{ printf("calculation start! \n"); } initialize1(alpha,s,&nspin,&P,theta); H[0] = Hx_max; H[1] = Hy_max; T = T_max; mcs = 1; magnetization_sum[2] = 0.0; energy_sum = 0.0; magnetization[2] = 0.0; energy = 0.0; //magnetization_sum2 = 0.0; //energy_sum2 = 0.0; //magnetization_2 = 0.0; //energy_2 = 0.0; //magnetic_susceptibility = 0.0; //specific_heat = 0.0; fprintf(fp1,"# ising modelb step\n"); fprintf(fp1,"# nspin=%d,J=%g,T=%g,(Hx,Hy)=(%g,%g),border=%d \n",nspin,J,T,H[0],H[1],border); fprintf(fp1,"# step T E[i] Mx My C[i] x[i] \n"); fprintf(fp2,"# ising modelb step\n"); fprintf(fp2,"# nspin=%d,J=%g,T=%g,(Hx,Hy)=(%g,%g),border=%d \n",nspin,J,T,H[0],H[1],border); fprintf(fp2,"# x y x + dx y + dy \n"); for( istep = 1 ; istep <= nmcs ; istep++ ){ for( i = 0 ; i < nspin ; i++ ){ for( j = 0 ; j < nspin ; j++ ){ //ランダムにスピンを選ぶ ispin = (int)(nspin * ransu()); jspin = (int)(nspin * ransu()); hanten = ransu(); check = ((int)(hanten*P))*theta; while(check == alpha[ispin][jspin]){ hanten = ransu(); check = ((int)(hanten*P))*theta; } dE = delta_E(ispin,jspin,alpha,H,check,&J,&nspin); if( dE < 0.0){ alpha[ispin][jspin] = check; } else{ if( hanten < exp(-dE/T)){ alpha[ispin][jspin] = check; } } }// j_loop }// i_loop if(mcs > border){ energy_sum += total_energy(alpha,H,&J,&nspin); //magnetization_sum += fabs(total_magnetization(s));//|M| magnetization_sum[0] += total_magnetization_x(alpha,&nspin);// M magnetization_sum[1] += total_magnetization_y(alpha,&nspin);// M magnetization[0] = magnetization_sum[0]/(mcs-border);// <|M|> or magnetization[1] = magnetization_sum[1]/(mcs-border);// <|M|> or energy = energy_sum/(mcs-border);// magnetization[0] /= (nspin*nspin);//=/spin magnetization[1] /= (nspin*nspin);//=/spin energy /= (nspin*nspin); fprintf(fp1,"%5d %5.2f %8.6f %8.6f %8.6f \n",mcs,T,energy,magnetization[0],magnetization[1]); for( i = 0 ; i < nspin ; i++ ){ for( j = 0 ; j < nspin ; j++ ){ s[i][j][0] = cos(alpha[i][j]); s[i][j][1] = sin(alpha[i][j]); fprintf(fp2,"%5.2f %5.2f %4.2f %4.2f\n",i-scale*s[i][j][0]/2.0,j-scale*s[i][j][1]/2.0,s[i][j][0]*scale,s[i][j][1]*scale); } } fprintf(fp2,"\n\n"); mcs += 1; } } fclose(fp1); fclose(fp2); return 0; }