Metropolis algorithm for 2D Ising model

[Vena]

Hi

at the moment i am writing my Bachelor theses about Montecarlo simulation for the 2D Ising Model with the Metropolis algorithm.
But my code does´t work. and i can´t finde my mistakes.
my measurement values are very far away from the analytical results.
Please help me to find my mistakes.
I am thankful for any advice and help.




[
//Ising_2D mit Metropolisalgorithmus


#include <iostream>
#include <cmath>
#include <ctime>
#include <random>
#include <fstream>
#include <cstdlib>
#include <stdio.h>

using namespace std;


const int Nx = 20; //Spins in direction x
const int Ny = 20;
const int N2 = 400; // all Spins
int MCSteps = 10000; // number of Metropolis Steps
const int k_B = 1; // Boltzmann constant
const int B = 0; //magnetic field
const int J = +1; // ferromagnetic coupling
const int mu = 1; //permeability
double beta; // invers temperatur
double T; // temperatur
int s[Nx][Ny]; //spin
int li[N2] , lj[N2] , ri[N2] , rj[N2]; //spin neightbours
int i; //spin position in direction x
int j; // spin position in direction y
double dE; // energy difference (E_new - E_old)
double E; // energy
double Cv; //specific heat
double a_Cv; // analytical results specific heat
double sus; // susceptibility
double a_sus; // analytical results susceptibility
double tprob; // transmission probability
double rprob; // random probability betwenn 0 and 1
double mag; // magnetisation
double med_mag; // mean magnetisation
double med_mag2; // mean squared magnetisation
double med_E; // mean energy
double med_E2; //mean squared energy
int count1; //counter
int count2; //counter



//generates an inital spin
//configuration randomly spin up (+1) spin down (-1)
void Initialisierung () //inizialisation


// random start +1 oder -1

{ random_device seed;
mt19937 engine(seed());
uniform_real_distribution<double> values(0.0, 1.0 );
for(int i=0; i<Nx; ++i)
for(int j=0; j<Ny; ++j)
{
if(values(engine)<0.5) { s[i][j] = 1;}
else { s[i][j] = -1;}
}

//cout<<"starting position:"<<s[i][j]<<'\n';

}

// measurement of the magnetisation
double Magnetisierung() // magnetisation
{
int Magnetisierung = 0;
for(int i =0; i<Nx; i++)
for(int j=0; j<Ny; j++)
{
Magnetisierung += s[i][j];
}
return Magnetisierung;
}

//measurement of the energy
double Energie() // energy
{

int sumN =0;

for(int i=0; i<Nx; i++)
for(int j=0; j<Ny; j++)
{
mag += s[i][j];
int ri = i == Nx-1 ? 0 : i+1;
int rj = j == Ny-1 ? 0 : j+1;
sumN += s[i][j] * (s[ri][j] + s[i][rj]);


}
return -(J*sumN + B*mu*mag);
}
//One of N2 spins is choosen randomly
//defines next neightbours
//implement periodic boundarys
//makes one metropolis step
void OneMetropolisStep()
{
random_device seed;
mt19937 engine(seed());
uniform_int_distribution<int> value1(0, Nx);
int i = value1(engine);
uniform_int_distribution<int> value2(0, Ny);
//cout <<"i:"<<i<<'\n';
int j = value2(engine);
//cout<<"j:"<<j<<'\n';

// nearest neightbours
int li;
if(i==0)
li = Nx-1;
else li = i-1;

int ri;
if(i == Nx-1)
ri = 0;
else ri = i+1;

int lj;
if(j==0)
lj= Ny-1;
else lj = j-1;

int rj;
if(j == Ny-1)
rj = 0;
else rj = j+1;

mag = Magnetisierung();
E = Energie();

//mag= N2*s[i][j];
//E=-J*(s[i][ri]*s[j][rj]) - B*mu*mag;

beta=1/(k_B*T);

dE = 2*J*s[i][j]*(s[ri][j] + s[li][j] + s[i][rj] + s[i][lj]) - 2*mu*B*s[i][j];

if(dE<=0)

{
s[i][j]=-s[i][j];

mag= mag+2*s[i][j];
E=E+dE;
}
else
{ random_device seed;
mt19937 engine(seed());
uniform_real_distribution<double> values(0.0, 1.0 );
tprob=exp(-beta*dE);
rprob= values(engine);
if (rprob<=tprob)
{
s[i][j]=-s[i][j];

mag=mag+2*s[i][j];
E=E+dE;
}
else
{
s[i][j]=s[i][j];
mag = mag;
E=E;

}

}


}


// N2 Metropolis Steps
// magnetisation per spin
// energy per spin
void oneMetropolisStepPerSpin()
{
for (int i =0; i<N2; i++)
{
OneMetropolisStep();

}
mag = mag / N2;
E = E / N2;

}



int main(int argc, const char * argv[])
{
for(count1=0; count1<MCSteps; count1++)
{Initialisierung();}
//T=0.1;
for (T=0.1; T<=5; T+=0.1)
{cout << " Temperatur: " << T << '\n';
beta = 1/(k_B*T);

med_mag =0;
med_E =0;
med_mag2 =0;
med_E2 =0;

int MM = int(0.4* MCSteps); //perform thermalization steps to let the system come to equilibrium
for(int count1=0; count1< MM; count1++)
{ oneMetropolisStepPerSpin();

med_mag += mag;
med_mag2 += mag*mag;
med_E+= E;
med_E2 += E*E;

}

med_mag=med_mag/(MM);
med_mag2=med_mag2/(MM);
med_E=med_E/(MM);
med_E2=med_E2/(MM);

Cv=(med_E2-med_E*med_E)/(T*T); //specific heat
sus=N2*(med_mag2-med_mag*med_mag)/(T); // susceptibility

// analytical solution:
a_Cv = (k_B* pow(beta, 2)* pow(J,2))/(pow(cosh(beta * J),2));
a_sus = beta*exp(2*beta*J);


cout <<"mean magnetisation "<<med_mag<<'\n';
cout <<"mean energy: "<<med_E<<'\n';
cout <<"mean squared magnetisation :"<<med_mag2<<'\n';
cout <<"mean squared energy: "<<med_E2<<'\n';
cout <<"specific heat: "<<Cv<<'\n';
cout <<"analytical Cv:"<<a_Cv<<'\n';
cout <<"susceptibility : "<<sus<<'\n';
cout <<"analytical sus: "<<a_sus<<'\n';




}



return 0;

}






]

[48T7M4Gy]

http://www.cplusplus.com/forum/general/195277/

Please don't double post.

[Naughty Albatross]

TheIdeasMan wrote:

Just a note for the future :+)

There also shouldn't be 2 topics about essentially the same subject. There is a risk that the same things will be said in both, possibly making it a time waster for those who reply.