Monte Carlo Investigation of Frustrated Heisenberg Model
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 ;
}