/*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 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<L;i++)
{
for(j=0;j<L;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<N_t; 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<L;i++)
{
for(j=0;j<L;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<L;i++)
{
for(j=0;j<L;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;
}