- Forum
- General C++ Programming
- Ising 2D with Metropolis Algorithm
Ising 2D with Metropolis Algorithm
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;
}
]
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;
}
]
closed account (48T7M4Gy)
Astoundingly 'your' program runs with the following output. What's wrong with that, looks pretty good to me.
Temperatur: 0.1 mean magnetisation -1.9941 mean energy: -1.99713 mean squared magnetisation :3.98064 mean squared energy: 3.98977 specific heat: 0.12358 analytical Cv:8.24461e-07 susceptibility : 16.8616 analytical sus: 4.85165e+09 Temperatur: 0.2 mean magnetisation -2 mean energy: -2 mean squared magnetisation :4 mean squared energy: 4 specific heat: 0 analytical Cv:0.00453958 susceptibility : 0 analytical sus: 110132 Temperatur: 0.3 mean magnetisation -2 mean energy: -2 mean squared magnetisation :4 mean squared energy: 4 specific heat: 0 analytical Cv:0.0564178 susceptibility : 0 analytical sus: 2619.24 Temperatur: 0.4 mean magnetisation -2 mean energy: -2 mean squared magnetisation :4 mean squared energy: 4 specific heat: 0 analytical Cv:0.166201 susceptibility : 0 analytical sus: 371.033 Temperatur: 0.5 mean magnetisation -2 mean energy: -2 mean squared magnetisation :4 mean squared energy: 4 specific heat: 0 analytical Cv:0.282603 susceptibility : 0 analytical sus: 109.196 Temperatur: 0.6 mean magnetisation -1.99999 mean energy: -1.99999 mean squared magnetisation :3.99998 mean squared energy: 3.99996 specific heat: 5.55278e-07 analytical Cv:0.369541 susceptibility : 3.33167e-05 analytical sus: 46.7194 Temperatur: 0.7 mean magnetisation -1.99996 mean energy: -1.99991 mean squared magnetisation :3.99983 mean squared energy: 3.99966 specific heat: 3.45464e-06 analytical Cv:0.419293 susceptibility : 0.000241825 analytical sus: 24.8739 Temperatur: 0.8 mean magnetisation -1.99985 mean energy: -1.9997 mean squared magnetisation :3.99941 mean squared energy: 3.99883 specific heat: 9.39527e-06 analytical Cv:0.438148 susceptibility : 0.000751622 analytical sus: 15.2281 Temperatur: 0.9 mean magnetisation -1.99946 mean energy: -1.99893 mean squared magnetisation :3.99783 mean squared energy: 3.99576 specific heat: 2.73652e-05 analytical Cv:0.43562 susceptibility : 0.00266026 analytical sus: 10.2531 Temperatur: 1 mean magnetisation -1.99833 mean energy: -1.99673 mean squared magnetisation :3.99334 mean squared energy: 3.98701 specific heat: 6.84484e-05 analytical Cv:0.419974 susceptibility : 0.00738527 analytical sus: 7.38906 Temperatur: 1.1 mean magnetisation -1.99717 mean energy: -1.99454 mean squared magnetisation :3.98871 mean squared energy: 3.97832 specific heat: 0.000100447 analytical Cv:0.397188 susceptibility : 0.0125325 analytical sus: 5.60059 Temperatur: 1.2 mean magnetisation -1.99417 mean energy: -1.98873 mean squared magnetisation :3.97679 mean squared energy: 3.95527 specific heat: 0.000162504 analytical Cv:0.371194 susceptibility : 0.0220648 analytical sus: 4.41208 Temperatur: 1.3 mean magnetisation -1.98948 mean energy: -1.98004 mean squared magnetisation :3.95817 mean squared energy: 3.92101 specific heat: 0.000255354 analytical Cv:0.344415 susceptibility : 0.0408614 analytical sus: 3.58263 Temperatur: 1.4 |
Last edited on
> My measurement values are very far away from the analytical results.
Can you please show us "the analytical results"?
Can you please show us "the analytical results"?
the analytical results are behind the measurements in the program.
the analytical result for the specific heat is in my program called "analytical Cv".
and the analytical result for the susceptibility is called "analytical sus".
Sorry that i posted this in both forums. But I did´t know where to post it.
Thank you for your help.
the analytical result for the specific heat is in my program called "analytical Cv".
and the analytical result for the susceptibility is called "analytical sus".
Sorry that i posted this in both forums. But I did´t know where to post it.
Thank you for your help.
Topic archived. No new replies allowed.