/*This program calculates the energy, heat capacity, magnetization, magnetic permitivity.*//*Boundary condition: periodic boundary condition.*//*In this program, J&k are put into T. To fill it to the expressions, just consider the dimensions.*/#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#define L 32
#define NGL 100000 /*Samples Neglected*/
#define N_t 5100000 /*Total samples*/
#define N (N_t-NGL)
#define T_i 0.5
#define T_f 10.5
#define Tstep 0.05
#define SIZE_T 500 /*The size of emcx in T axis. Check this before compiling!*/
#define SIZE_R 60 /*The size of emcx in R axis. Check this before compiling!*/
#define R_i -6 /*Ration of J2/J1. Initial value.*/
#define R_step 0.1 /*step*/
#define R_f 0 /*End of R*/intsum_M;ints[L][L];/*Size of the lattice*/doublesum_E;doubleT;structObsv{doubleC;doubleX;doubleM;doubleE;};voidInist();voidMag();voidEnergy(doubleR);voidMCP(doubleR);structObsvCac_obsv(doubleR);intmain(){inti,j;doubleR;structObsvemcx[SIZE_T][SIZE_R];FILE*ft,*fm,*fx,*fc,*fe,*fr;fr=fopen("R.txt","w");ft=fopen("T.txt","w");fx=fopen("X.txt","w");fc=fopen("C.txt","w");fm=fopen("M.txt","w");fe=fopen("E.txt","w");R=R_i;srand(time(NULL));for(j=0;R<=R_f;j++){T=T_i;for(i=0;T<=T_f;i++){Inist();emcx[i][j]=Cac_obsv(R);fprintf(fr,"%f\n",R);fprintf(ft,"%f\n",T);fprintf(fx,"%f\n",emcx[i][j].X);fprintf(fc,"%f\n",emcx[i][j].C);fprintf(fm,"%f\n",emcx[i][j].M);fprintf(fe,"%f\n",emcx[i][j].E);printf("-------------\n");printf("R=%f\n",R);printf("T=%f\n",T);printf("X=%f\n",emcx[i][j].X);printf("C=%f\n",emcx[i][j].C);printf("M=%f\n",emcx[i][j].M);printf("E=%f\n",emcx[i][j].E);T=T+Tstep;}R+=R_step;}fclose(ft);fclose(fx);fclose(fc);fclose(fm);fclose(fe);printf("The End!\n");}voidInist(){inti,j;for(i=0;i<L;i++){for(j=0;j<L;j++){s[i][j]=1;}}}voidMCP(doubleR){inti,j;doublesum,dE;i=rand()%L;j=rand()%L;sum=s[i][j]*(s[(L+i-1)%L][j]+s[(i+1)%L][j]+s[i][(j-1+L)%L]+s[i][(j+1)%L]+R*(s[(i+1)%L][(j-1+L)%L]+s[(i-1+L)%L][(j+1)%L]));//printf("chk2.sum=%d\n",sum);//switch(sum)//{// case 4: dE= 8; break;// case 2: dE= 4; break;// case 0: dE= 0; break;// case -2: dE= -4; break;// case -4: dE= -8; break;// default: printf("sum=%d. Attention!\n",sum); break;//}dE=2*sum;if(dE>0){if(((double)rand()/RAND_MAX)<=exp((-1.0)*dE/T)){s[i][j]=-s[i][j];sum_M+=2*s[i][j];sum_E+=dE;}}else{s[i][j]=-s[i][j];sum_M+=2*s[i][j];sum_E+=dE;}}structObsvCac_obsv(doubleR){inti;doubleE_avg=0.0;doubleE_as=0.0;doubleM_avg=0.0;doubleM_as=0.0;structObsvemcx;Mag();Energy(R);for(i=0;i<N_t;i++){MCP(R);if(i>=NGL){E_avg+=(double)sum_E/N;E_as+=(double)sum_E*sum_E/N;M_avg+=(double)sum_M/N;M_as+=(double)sum_M*sum_M/N;}}emcx.E=E_avg;emcx.M=M_avg;emcx.C=(E_as-E_avg*E_avg)/T/T;emcx.X=(M_as-M_avg*M_avg)/T;returnemcx;}voidMag(){inti,j;intdM_tmp=0;for(i=0;i<L;i++){for(j=0;j<L;j++)dM_tmp+=s[i][j];}sum_M=dM_tmp;}voidEnergy(doubleR){doubledE_tmp=0.0;inti,j;for(i=0;i<L;i++){for(j=0;j<L;j++){dE_tmp+=0.5*(s[i][j]*(s[(i-1+L)%L][j]+s[(i+1)%L][j]+s[i][(j-1+L)%L]+s[i][(j+1)%L])+R*s[i][j]*(s[(i+1)%L][(j-1+L)%L]+s[(i-1+L)%L][(j+1)%L]));// printf("chk1.dE_tmp=%f\n",dE_tmp);}}sum_E=dE_tmp;}