File contents

//2D ising model,b_c=2.269
//L=10^1,tmax=10^4,ens=10^1 ==> 46"


# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# include "ran2.c"

# define L          10
# define tmax       10000
# define ens        10

int main()
{

int    i,j,k,e,t,m,n,p,p1;
double y[L][L]={},eng1,eng2,eng,q,b,E;
long   sd =-99676;

FILE *fp ;

fp=fopen("file1","w");

for(b=1;b<5;b+=0.1){ E=0;

for(e=0;e<ens;e++){

//initialization
for(i=0;i<L;i++)for(j=0;j<L;j++){y[i][j]=0;}
for(i=0;i<L;i++)for(j=0;j<L;j++){y[i][j]=-1;}
//for(i=0;i<L;i++)for(j=0;j<L;j++){if(ran2(&sd)<0.5)y[i][j]=1;}
k=0; while (k<L*L*0.5){ i=ran2(&sd)*L; j=ran2(&sd)*L; if(y[i][j]==-1){y[i][j]=1; k++ ;} }

//print matrix
//for(i=0;i<L;i++){for(j=0;j<L;j++){printf("   %d   ",(int)y[i][j]);}printf("\n");}

//dynamics starts
for(t=0;t<tmax;t++){

for(i=0;i<L;i++)for(j=0;j<L;j++){
m=ran2(&sd)*L; n=ran2(&sd)*L;
p = y[(m-1+L)%L][n] + y[m][(n+1)%L] + y[m][(n-1+L)%L] + y[(m+1)%L][n];
//printf("%d %d %d\n",m,n,p);
eng1 = 1.0*p*y[m][n];
eng2 = -1.0*p*y[m][n];
eng = eng2 - eng1;

if(eng<0){y[m][n]=-y[m][n];} 
else if(eng>0){q=exp(-eng/b);if(ran2(&sd)<q){y[m][n]=-y[m][n];} }
                                }

//calculate energy
for(i=0;i<L;i++)for(j=0;j<L;j++){
p1 = y[(i-1+L)%L][j] + y[i][(j+1)%L] + y[i][(j-1+L)%L] + y[(i+1)%L][j];
E += 1.0*p1*y[i][j];
                                }

                    }//tmax ends
            }//ens ends

fprintf(fp,"%lf  %lf\n",b,E*1.0/tmax/ens/L/L);

    }//b loop ends

fclose(fp);
}//main ends