2D Ising Model (Triangle Lattice)

Q: Using the simplest program to calculate the 2D Ising model upon triangle lattice.
A:
1. Hamiltonian:

[eq]H=J\sum\limits_{i,j}{{s_{i}}{s_{j}}}[/eq]

2. Lattice:
Fig3-Triangle-2D-Ising

Fig3-Triangle-2D-Ising

3. Summary image:

2D_Triangle_Ising_Model

2D_Triangle_Ising_Model

——CODE&SUMMARY——

*PDF file (some results shown in it):
2D_Ising_Model(Triangle)
(A little reminder: I am not really working on this kind of things. The only reason I did this is that I have to hand in a copy of my homework. So IF U R LOOKING FOR SOME REAL INSPIRING CODES, NEVER STUDY MINE!)

*CODE file (html):
CODE_2DIsingModel-Triangle

My CODE is appended here:

/*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
#include
#include

#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*/

int sum_M;
int s[L][L]; /*Size of the lattice*/
double sum_E;
double T;

struct Obsv
{
double C;
double X;
double M;
double E;
};

void Inist();
void Mag();
void Energy(double R);
void MCP(double R);
struct Obsv Cac_obsv(double R);

int main()
{
int i,j;
double R;
struct Obsv emcx[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");
}

void Inist()
{
int i,j;
for(i=0;i {
for(j=0;j {
s[i][j]=1;
}
}
}

void MCP(double R)
{
int i,j;
double sum,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;
}
}

struct Obsv Cac_obsv(double R)
{
int i;
double E_avg=0.0;
double E_as=0.0;
double M_avg=0.0;
double M_as=0.0;
struct Obsv emcx;
Mag();
Energy(R);

for(i=0; 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;

return emcx;
}

void Mag()
{
int i,j;
int dM_tmp=0;
for(i=0;i {
for(j=0;j dM_tmp+=s[i][j];
}
sum_M=dM_tmp;
}

void Energy(double R)
{
double dE_tmp=0.0;
int i,j;
for(i=0;i {
for(j=0;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;
}