A simple vector model in 2 dimensions. Its Hamiltonian can be written as
$$ H={ {J} { 1} }\sum\limits{i,j}{ { {s}{i} }{ {s}{j} } }+{ {J}{2} }\sum\limits{i’,j’}{ { { s}{i’} }{ {s}{j’} } }, $$
in which i,j means the nearest neighbor interaction when i’, j’ means the next nearest neighbor interaction. We designed the program under Honeycomb lattice, i.e.,
Scheme of a honeycomb lattice
To make it clear for programming, we reshaped the lattice under homeomorphic：
Transformation of honeycomb lattice for MC program
We wrote the Metropolis CODE with C:
/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 100
#define NGL 1000000 /*Samples Neglected*/
#define N_t 51000000 /*Total samples*/
#define N (N_t-NGL)
#define T_i 0.1
#define T_f 1.0
#define Tstep 0.1
#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 50.0 /*Ration of J2/J1. Initial value.*/
#define R_step 0.1 /*step*/
#define R_f 50.0 /*End of R*/
#define Pi 3.1415926536
double sum_M[2]; /*sum_M[0] denotes the cos part; sum_M[1] denotes the sin part*/
double sum_Mag;
double s[L][L]; /*Size of the lattice*/
double sum_E;
double T;
struct Obsv
{
double C;
double X;
double M;
double E;
};
/********************************/
// FILE *fchk; //chk_sum_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");
/**********************************/
//fchk=fopen("Chk_dE.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(fr);
fclose(ft);
fclose(fx);
fclose(fc);
fclose(fm);
fclose(fe);
/*****************************/
//fclose(fchk);
/*****************************/
printf("The End!\n");
}
void Inist() /*Initialized lattice. All oriented to one direction*/
{
int i,j;
for(i=0;i {
for(j=0;j {
s[i][j]=0;
}
}
}
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);
/********************************/
//printf("chk.Cac_obsv.sum_Mag=%f\n",sum_Mag);
//printf("chk.Cac_obsv.sum_E=%f\n",sum_E);
//fprintf(fchk,"%f\n",sum_E);
/********************************/
if(i>=NGL)
{
E_avg+=(double)sum_E/N;
E_as+=(double)sum_E*sum_E/N;
M_avg+=(double)sum_Mag/N;
M_as+=(double)sum_Mag*sum_Mag/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;
}
/*Monte Carlo process.*/
void MCP(double R)
{
int i,j;
double agl_rdm;
double sum,sum_tmp,dE; /*sum_tmp is used to cal the energy difference*/
i=rand()%L;
j=rand()%L;
agl_rdm=(double)rand()/RAND_MAX*2*Pi; /*What if I just make it (double)rand()*/
/*Cal the energy before we add a random angle to s[i][j]*/
if(j%2==0)
{
if(i%2==0)
{
sum_tmp=cos(s[i][j]-s[i][(j-1+L)%L])+cos(s[i][j]-s[(i+1)%L][j])+cos(s[i][j]-s[(i-1+L)%L][j]);
sum_tmp+=R*(cos(s[i][j]-s[(i+1)%L][(j-1+L)%L])+cos(s[i][j]-s[(i+1)%L][(j+1)%L])+cos(s[i][j]-s[(i-1+L)%L][(j+1)%L])+cos(s[i][j]-s[(i-1+L)%L][(j-1+L)%L])+cos(s[i][j]-s[(i+2)%L][j])+cos(s[i][j]-s[(i-2+L)%L][j]));
}
else if(i%2==1)
{
sum_tmp=cos(s[i][j]-s[i][(j+1)%L])+cos(s[i][j]-s[(i+1)%L][j])+cos(s[i][j]-s[(i-1+L)%L][j]);
sum_tmp+=R*(cos(s[i][j]-s[(i+1)%L][(j-1+L)%L])+cos(s[i][j]-s[(i+1)%L][(j+1)%L])+cos(s[i][j]-s[(i-1+L)%L][(j+1)%L])+cos(s[i][j]-s[(i-1+L)%L][(j-1+L)%L])+cos(s[i][j]-s[(i+2)%L][j])+cos(s[i][j]-s[(i-2+L)%L][j]));
}
}
else if(j%2==1)
{
if(i%2==1)
{
sum_tmp=cos(s[i][j]-s[i][(j-1+L)%L])+cos(s[i][j]-s[(i+1)%L][j])+cos(s[i][j]-s[(i-1+L)%L][j]);
sum_tmp+=R*(cos(s[i][j]-s[(i+1)%L][(j-1+L)%L])+cos(s[i][j]-s[(i+1)%L][(j+1)%L])+cos(s[i][j]-s[(i-1+L)%L][(j+1)%L])+cos(s[i][j]-s[(i-1+L)%L][(j-1+L)%L])+cos(s[i][j]-s[(i+2)%L][j])+cos(s[i][j]-s[(i-2+L)%L][j]));
}
else if(i%2==0)
{
sum_tmp=cos(s[i][j]-s[i][(j+1)%L])+cos(s[i][j]-s[(i+1)%L][j])+cos(s[i][j]-s[(i-1+L)%L][j]);
sum_tmp+=R*(cos(s[i][j]-s[(i+1)%L][(j-1+L)%L])+cos(s[i][j]-s[(i+1)%L][(j+1)%L])+cos(s[i][j]-s[(i-1+L)%L][(j+1)%L])+cos(s[i][j]-s[(i-1+L)%L][(j-1+L)%L])+cos(s[i][j]-s[(i+2)%L][j])+cos(s[i][j]-s[(i-2+L)%L][j]));
}
}
else
{
printf("Checkpoint 1! j is neither odd nor even! Attention!");
}
/****************************/
//printf("chk.MCP.sum_tmp=%f\n",sum_tmp);
/****************************/
/*Add a random angle to s[i][j]*/
s[i][j]+=agl_rdm; /*Exceeding 2*Pi makes not sense in cos or sin func*/
/*Cal the energy after we add a random to s[i][j]*/
if(j%2==0)
{
if(i%2==0)
{
sum=cos(s[i][j]-s[i][(j-1+L)%L])+cos(s[i][j]-s[(i+1)%L][j])+cos(s[i][j]-s[(i-1+L)%L][j]);
sum+=R*(cos(s[i][j]-s[(i+1)%L][(j-1+L)%L])+cos(s[i][j]-s[(i+1)%L][(j+1)%L])+cos(s[i][j]-s[(i-1+L)%L][(j+1)%L])+cos(s[i][j]-s[(i-1+L)%L][(j-1+L)%L])+cos(s[i][j]-s[(i+2)%L][j])+cos(s[i][j]-s[(i-2+L)%L][j]));
}
else if(i%2==1)
{
sum=cos(s[i][j]-s[i][(j+1)%L])+cos(s[i][j]-s[(i+1)%L][j])+cos(s[i][j]-s[(i-1+L)%L][j]);
sum+=R*(cos(s[i][j]-s[(i+1)%L][(j-1+L)%L])+cos(s[i][j]-s[(i+1)%L][(j+1)%L])+cos(s[i][j]-s[(i-1+L)%L][(j+1)%L])+cos(s[i][j]-s[(i-1+L)%L][(j-1+L)%L])+cos(s[i][j]-s[(i+2)%L][j])+cos(s[i][j]-s[(i-2+L)%L][j]));
}
}
else if(j%2==1)
{
if(i%2==1)
{
sum=cos(s[i][j]-s[i][(j-1+L)%L])+cos(s[i][j]-s[(i+1)%L][j])+cos(s[i][j]-s[(i-1+L)%L][j]);
sum+=R*(cos(s[i][j]-s[(i+1)%L][(j-1+L)%L])+cos(s[i][j]-s[(i+1)%L][(j+1)%L])+cos(s[i][j]-s[(i-1+L)%L][(j+1)%L])+cos(s[i][j]-s[(i-1+L)%L][(j-1+L)%L])+cos(s[i][j]-s[(i+2)%L][j])+cos(s[i][j]-s[(i-2+L)%L][j]));
}
else if(i%2==0)
{
sum=cos(s[i][j]-s[i][(j+1)%L])+cos(s[i][j]-s[(i+1)%L][j])+cos(s[i][j]-s[(i-1+L)%L][j]);
sum+=R*(cos(s[i][j]-s[(i+1)%L][(j-1+L)%L])+cos(s[i][j]-s[(i+1)%L][(j+1)%L])+cos(s[i][j]-s[(i-1+L)%L][(j+1)%L])+cos(s[i][j]-s[(i-1+L)%L][(j-1+L)%L])+cos(s[i][j]-s[(i+2)%L][j])+cos(s[i][j]-s[(i-2+L)%L][j]));
}
}
else
{
printf("Checkpoint 2! j is neither odd nor even! Attention!");
}
/********************************/
//printf("chk.MCP.sum=%f\n",sum);
/********************************/
dE=sum-sum_tmp; /*Attention!*/
/********************************/
//printf("chk.MCP.dE=%f\n",dE);
/********************************/
if(dE< =0)
{
sum_M[0]+=cos(s[i][j])-cos(s[i][j]-agl_rdm); /*Accept&Cal*/
sum_M[1]+=sin(s[i][j])-sin(s[i][j]-agl_rdm);
sum_Mag=sqrt(sum_M[0]*sum_M[0]+sum_M[1]*sum_M[1]);
sum_E+=dE;
}
else if(((double)rand()/RAND_MAX) {
sum_M[0]+=cos(s[i][j])-cos(s[i][j]-agl_rdm); /*Accept&Cal*/
sum_M[1]+=sin(s[i][j])-sin(s[i][j]-agl_rdm);
sum_Mag=sqrt(sum_M[0]*sum_M[0]+sum_M[1]*sum_M[1]);
sum_E+=dE;
}
else
{
s[i][j]=s[i][j]-agl_rdm; /*Reject the change*/
}
/********************************/
//printf("chk.MCP.sum_M[0]=%f\n",sum_M[0]);
//printf("chk.MCP.sum_M[1]=%f\n",sum_M[1]);
//printf("chk.MCP.cos(s[i][j])-cos(s[i][j]-agl_rdm)=%f\n",cos(s[i][j])-cos(s[i][j]-agl_rdm));
/********************************/
}
void Mag()
{
// int i,j;
// double dM_tmp[2]={0};
// for(i=0;i // {
// for(j=0;j // dM_tmp[0]+=cos(s[i][j]);
// dM_tmp[1]+=sin(s[i][j]);
// }
// sum_M[0]=dM_tmp[0];
// sum_M[1]=dM_tmp[1];
sum_M[0]=1.0*L*L;
sum_M[1]=0.0;
/********************************/
//printf("chk.Mag.sum_M[0]=%f\n",sum_M[0]);
//printf("chk.Mag.sum_M[1]=%f\n",sum_M[1]);
/********************************/
}
void Energy(double R)
{
// double dE_tmp;
// dE_tmp=(3.0+R*6.0)*L*L/2.0;
sum_E=(3.0+R*6.0)*L*L/2.0;
}