Monte Carlo Investigation of Frustrated Heisenberg Model

· by · in · [ Co-Reading ] ·

Draft ( PDF file, click to view or download):
Frustrated-Honeycomb Heisenberg-Model

We did not do much on this because of the limitation of time.
One should be careful of the possible mistakes in this draft.

Here we have some slides to show the main idea of this work.
2D_Honeycomb_Frustrated_Heisenberg_Model_Monte Carlo

Codes:
原始的：Honeycomb-2DIsingModel
大致看了下类似刚度之类的性质：Rigid-Honeycomb-2DIsingModel
OPENGL的可视化：Visual-Honeycomb
但是这个可视化不太好，跟slide里面那个可视化不太一样。这里使用颜色深度表示方向的。我在slides里面其实只是把某时刻的状态记录到一个文件中，然后用sigmaplot画的……囧~这个的code在这里：Final-pirntspin-Honeycomb
三角格子的（因为在某极限下本模型趋向于三角格子，所以就单独算了一下三角格子的作为对照）：Rvaries-Triangle-2DIsingModel

============= Divider of Main Text================
A simple vector model in 2 dimensions.
Its Hamiltonian can be written as

[eq]H={{J}{1}}\sum\limits{i,j}{{{s}{i}}{{s}{j}}}+{{J}{2}}\sum\limits{i’,j’}{{{s}{i’}}{{s}{j’}}}[/eq]

,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

Scheme of a honeycomb lattice

To make it clear for programming, we reshaped the lattice under homeomorphic：

Transformation of honeycomb lattice

Transformation of honeycomb lattice for MC program

I 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

include

include

include

#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;
}

长程相关:

  1. 2D Ising Model (Triangle Lattice)
  2. 关于Grandcanonical ensemble的巨配分函数的一些讨论
  3. A Brief Summary of Statistical Physics
  4. 相对论心理学1
Tags: , ,