Scripts | Codes

/***********************************/
/*  Find more codes and scripts on
/*   codes-n-scripts.blogspot.com
/*  ~Ising model~
/***********************************/

#include 
#include 
#include 
#include 
#include 

#define N 30
#define M 30
#define nmax  100000
#define nzero  40000

double rg_1() ;
int rg_2() ;
int rg_3() ;

int main()
{

 FILE *output_file;
 output_file = fopen("output.dat", "w");

int i,j,n;
int i1,j1,i2,i3,i4,j4,i5,j5,i6,j6,i7,i8,i9,i10; 
int s1,s2,randsigma ;
int sigma[N][M] ;
int H1,H2,H2a,H2b,H2c,H2d,H3,H3a,H3b,H3c,H3d,H4;
double H[nmax] ;
double Ma[nmax] ;
double r,z ;
double deltaE ;
double J,B,beta,T;
double E, sigmaE, Mag, sigmaMag;

/* parameters */

J = 1.0 ;
B = 1.0 ;
T = 3.5 ;

beta =1.0/T  ;

/* initial configuration */
 
for(i = 0; i < N; i++){  
   for(j = 0; j < M; j++){
     sigma[i][j] = rg_2()*2-1  ; 
  }
}

printf("\n" );
printf("Lattice in:\n" );
printf("\n" );

for(i5 = 0; i5 < N; i5++){  
   for(j5 = 0; j5 < M; j5++){
     if(sigma[i5][j5]==-1){
     printf("  ");
     }else{
     printf("* ");
    }
}
printf("\n" ); 
}

printf("\n" ); 

for(n = 1; n < nmax+1; n++){


 if(n>2){

 s1 = rg_3() ;
 s2 = rg_3() ;
 
 randsigma = sigma[s1][s2] ;

 if(randsigma == 1 ){
 sigma[s1][s2] = -1 ;
 }
 else{
 sigma[s1][s2] = 1 ;
 }

 }

   /* calculation of the energy */
 
 H1=0;

for(i1 = 1; i1 < N-1 ; i1++){
  for(j1 = 1; j1 < M-1; j1++){ 
   H1 = sigma[i1][j1]*sigma[i1][j1+1]+sigma[i1][j1]*sigma[i1][j1-1]+sigma[i1][j1]*sigma[i1-1][j1]+sigma[i1][j1]*sigma[i1+1][j1]+ H1;
  }
}

H2=0; H2a=0; H2b=0; H2c=0; H2d=0;

 for(i2 = 1; i2 < N-1 ; i2++){
   H2a= sigma[i2][0]*sigma[i2][1]+sigma[i2][0]*sigma[i2][M-1]+sigma[i2][0]*sigma[i2-1][0]+sigma[i2][0]*sigma[i2+1][0]+ H2a;
   H2b= sigma[0][i2]*sigma[0][i2-1]+sigma[0][i2]*sigma[0][i2+1]+sigma[0][i2]*sigma[1][i2]+sigma[0][i2]*sigma[N-1][i2]+ H2b;
   H2c= sigma[i2][N-1]*sigma[i2-1][N-1]+sigma[i2][N-1]*sigma[i2+1][N-1]+sigma[i2][N-1]*sigma[i2][0]+sigma[i2][N-1]*sigma[i2][M-2]+ H2c;
   H2d= sigma[N-1][i2]*sigma[N-1][i2-1]+sigma[N-1][i2]*sigma[N-1][i2+1]+sigma[N-1][i2]*sigma[N-2][i2]+sigma[N-1][i2]*sigma[0][i2]+ H2d ;
 }

 H2=H2a+H2b+H2c+H2d; 
 
 H3=0; H3a=0; H3b=0; H3c=0; H3d=0;

 H3a=sigma[0][0]*sigma[0][1]+sigma[0][0]*sigma[0][N-1]+sigma[0][0]*sigma[M-1][0]+sigma[0][0]*sigma[1][0] ;
 H3b=sigma[N-1][N-1]*sigma[N-1][M-2]+sigma[N-1][N-1]*sigma[N-1][0]+sigma[N-1][N-1]*sigma[N-2][M-1]+sigma[N-1][N-1]*sigma[0][M-1] ;
 H3c=sigma[0][N-1]*sigma[0][N-2]+sigma[0][N-1]*sigma[0][0]+sigma[0][N-1]*sigma[1][N-1]+sigma[0][N-1]*sigma[M-1][N-1] ;
 H3d=sigma[N-1][0]*sigma[N-1][1]+sigma[N-1][0]*sigma[N-1][M-1]+sigma[N-1][0]*sigma[0][0]+sigma[N-1][0]*sigma[N-2][0] ;

 H3=H3a+H3b+H3c+H3d; 

 H4=0;
 for(i4 = 0; i4 < N; i4++){  
   for(j4 = 0; j4 < M; j4++){
     H4 = sigma[i4][j4]+H4 ; 
  }
}

 Ma[n-1] =( (double) H4 )/((double) N*M) ;  /* magnetisation */

 H[n-1] = -J*((double) H1+H2+H3 )-B*( (double) H4) ; /* energy */

 fprintf( output_file , " %d %f %f \n",n-1,H[n-1],Ma[n-1] );  /* print to the file */


/* flip the spin */

 if(n==1){

/* choose the spin randomly */

 s1 = rg_3() ;   
 s2 = rg_3() ;

 randsigma = sigma[s1][s2] ; 
  

 if(randsigma == 1 ){
 sigma[s1][s2] = -1 ;
 }
 else{
 sigma[s1][s2] = 1 ;
 }

 }
 
 if(n>1){

/* Metropolis  algorithm */
 
 deltaE = 0 ; 

 deltaE = H[n-1]- H[n-2] ;    

 r = exp( (double) -beta*deltaE ) ; 

 z = rg_1() ;

 /* printf("r=%f, z=%f, Energy[%d]=%f, Energy[%d]=%f\n", r,z,n-1,H[n-1], n-2,H[n-2]  ) ; */

if( H[n-1]> H[n-2]  && z>r ){

if(sigma[s1][s2]==1){
 sigma[s1][s2] = -1 ;
 }
 else{
 sigma[s1][s2] = 1 ;
}

 H[n-1]= H[n-2];

}

}

}

printf("\n" );
printf("Lattice out :\n" );
printf("\n" );

for(i6 = 0; i6 < N; i6++){  
   for(j6 = 0; j6 < M; j6++){
   if(sigma[i6][j6]==1){
     printf("* ");
     }else{ 
     printf("  ");
    }
 }
printf("\n" ); 
}

printf("\n" ); 


/*  energy  */

E=0.0 ; sigmaE=0.0 ;


for(i7 = nzero; i7 < nmax; i7++){
  E=H[i7]+E ;
  }

E = E/((double) nmax-nzero  ) ;

/*  energy dispersion */

for(i8 = nzero; i8 < nmax; i8++){
sigmaE =pow((double) E-H[i8],2.0) + sigmaE ; 
  }

sigmaE = sqrt( (double) sigmaE/((double) nmax-nzero-1)) ;

/*  magnetisation   */

Mag=0.0 ; sigmaMag=0.0 ;

for(i9 = nzero; i9 < nmax; i9++){
  Mag=Ma[i9]+Mag ;
  }

 Mag = Mag/((double) nmax-nzero  ) ;

/*  magnetisation  dispersion  */

for(i10 = nzero; i10 < nmax; i10++){
sigmaMag =pow((double) Mag-Ma[i10],2.0) + sigmaMag ; 
  }

sigmaMag = sqrt( (double) sigmaMag/((double) nmax-nzero-1)) ;

/*  print mean values */

printf("beta=%f, J=%f, B=%f \n",beta,J,B) ;

printf("\n" ); 

printf(" Energy = %f ± %f \n",E ,sigmaE );

printf("\n" ); 

printf(" Magnetisation = %f ± %f \n",Mag ,sigmaMag );

printf("\n" ); 

fclose(output_file);

return(0) ; 

}


/*  pseudo-random number generators  */


/*   0..1 real numbers generator   */

double rg_1()
{
  double ps;
  static int flag = 1 ;
      time_t date ;
      if(flag)
 { srand(time(&date) );
 flag=0;
 }
  ps =(double) rand()/(RAND_MAX+1.0) ;
  return(ps);
}

/*   O,1 integers generator    */

int rg_2()
{
  int ps;
  static int flag = 1 ;
      time_t date ;
      if(flag)
 { srand(time(&date) );
 flag=0;
 }
      ps =(int) rand()/(((unsigned)RAND_MAX+1.0)/2) ;
  return(ps);
}


/*   O,1,2...N-1 integers generator    */

int rg_3()
{
  int ps;
  static int flag = 1 ;
      time_t date ;
      if(flag)
 { srand(time(&date) );
 flag=0;
 }
      ps =(int) rand()/(((unsigned)RAND_MAX+1.0)/N) ;
  return(ps);
}