Solution: c-code of Monte Carlo algorithm for 2d-ising model

#include<stdio.h>

#include<math.h>

#include<stdlib.h>

#include <string.h>

#define L 128 // ^

#define LL (L*L) // | Lattice Parameters

#define L1 (L-1) // |

#define LL1 (L*L1) // v

#define T 1.25 // T= KT/J actually

#define p1 5000 // Average magnetization calculation starts here

// i.e., Steady state starts here

#define p2 10000 // Program stops at this step

#define N 10000 // No of configurations

static long iseed=-999999999;

float ran2(long *idum); // Random number generator: See Ref.8

int main()

{

static long I,j,k,ii,jj,step,isite,a(LL),nn,nns,sp,dE,Enrg=0, spn=0, time=0, sqr_Enrg=0, sqr_spn=0, [4]={1,L,-1,-L};

static double z,r,avg_mag, avg_energy, avg_sqr_mag, avg_sqr_energy

cie[2], mgt=0, energy=0,sqr_mgt=0, sqr_energy=0, dpLL, sp_heat, mag_sus;

cie[0]=exp(-4/T);

cie[1]=exp(-8/T);

dpLL=(p2-p1)*LL;

/* Arbitraryly fillup lattice sites */

for(i=0;i<LL;i++)

{

z=ran2(&iseed);

if (z>0.4) a[i]=1;

else a[i]=-1;

spn=spn+a[i];

}

/* Calculation of initial energy */

for(i=0;i<LL;i++)

{

for(k=0;k<4;k++)

{

nn=i+iv[k];

if(k==0 && i%L==L1)nn=i-L1;

if(k==1 && i>=LL1)nn=i-LL1;

if(k==2 && i%L==0)nn=i+L1;

if(k==3 && i<L)nn=i+LL1;

Enrg+=a[i]*a[nn];

}

}

/* Algorithm starts here */

for(ii=0;ii<p2;ii++)

{

for(step=0;step<LL;step++)

{

isite=ran2(&iseed)*LL; // Choosing any arbitrary site

sp=0;

for(i=0;i<4;i++)

{

nns=isite+iv[i];

if(i==2 && isite%L==0)nns=isite+L1; // |

if(i==3 && isite<L)nns=isite+LL1; // | Finding NNs

if(i==0 && isite%L==L1)nns=isite-L1;// |

if(i==1 && isite>=LL1)nns=isite-LL1;// |

sp+=a[nns];

}

/* Metropolis Algorithm */

dE=2*a[isite]*sp;

if(dE<=0)

{

a[isite]=-a[isite]; // Flipping criterion

spn=spn+2*a[isite]; // Change in magnetization & Energy

Enrg+=dE;

}

else

{

j=(dE/4)-1;

r=cie[j];

z=ran2(&iseed); //Calling random number

if (z<=r)

{

a[isite]=-a[isite];

spn=spn+2*a[isite]; //Change in magnetization & Energy

Enrg+=dE;

}

}

} // Here we will have one state.

time++;

// After p1 steps:

if(time>=p1)

{

for(jj=0;jj<N;jj++)

{

for(step=0;step<LL;step++)

{

isite=ran2(&iseed)*LL;

sp=0;

for(i=0;i<4;i++)

{

nns=isite+iv[i];

if(i==2 && isite%L==0)nns=isite+L1;

if(i==3 && isite<L)nns=isite+LL1;

if(i==0 && isite%L==L1)nns=isite-L1;

if(i==1 && isite>=LL1)nns=isite-LL1;

sp+=a[nns];

}

dE=2*a[isite]*sp;

if(dE<=0)

{

a[isite]=-a[isite];

spn=spn+2*a[isite];

Enrg+=dE;

}

else

{

j=(dE/4)-1;

r=cie[j];

z=ran2(&iseed);

if (z<=r)

{

a[isite]=-a[isite];

spn=spn+2*a[isite];

Enrg+=dE;

}

}

}

}

mgt=mgt+(double)spn;

sqr_spn=spn*spn;

sqr_mgt=sqr_mgt+(double)sqr_spn;

energy=energy+(double)Enrg;

sqr_Enrg=Enrg*Enrg;

sqr_energy=sqr_energy+(double)sqr_Enrg;

}

// Calculating the steady state average quantity

// per spin between p1 and p2 steps

avg_mag=mgt/dpLL;

avg_sqr_mag=sqr_mgt/(dpLL*LL);

avg_energy=energy/dpLL;

avg_sqr_energy=sqr_energy/(dpLL*LL);

sp_heat=(avg_sqr_energy-avg_energy*avg_energy)/(T*T);

mag_sus=(avg_sqr_mag-avg_mag*avg_mag)/T;

// End of the Program //

return(0);

}