Ising Model - code¶
This code is done using C++.
Here we present the code to simulate the ising model saving the state configuration so the thermodynamic quantities are calculated after the simulation using this notebook. The code is constructed in such a way that is compatible with the way the analysis is going to be done, but there are other ways to approach this problem.
First, This code is constructed to simulate a certain amount of monte carlos steps an Ising model at a given temperature and field recieved from the terminal, this would mean that in principle this same code can be used to simulate consatnt temperature conditions but that particular case present an interesting behavior when the configurations are not independent,i.e. the final state of a simulation is the initial of the followng, because a hysteresis cycle can be found there, this code and analysis can be found at link
Code¶
First, let us import the libraries we are going to use,
#include<iostream> // To Print
#include<cstdlib> // To Convert Data Types
#include<vector> // To Construct the Matrix
#include<cmath> // To Calculate Mathematical Expresions
In this part, we define the size of the grid we are about to use and also how many montecarlo steps we are going to use. $J$ is the constant interaction that we have select to have arbitrary units.
const int Nx=100;
const int Ny=100;
const int N=100;
const double J=1.;
We define a matrix structure as a vector which each compound is also a vector. So here is where our system will be,
typedef std::vector<std::vector<int> > matrix;
We need to have some functions which are going to be useful in our code, first, we need to initialize the matrix, as this particular code is done in such a way that the matrix has to me allocated in memory, the function that do both, the allocation and initialization is called get_mem.
Then we may need to print the matrix on the terminal.
And finally the function Step is not what is known as a MonteCarlo step, but just the metropolis algorithm for a single cell.
void get_mem(matrix & m);
void print(const matrix & m);
void Step(matrix & m, double kbT,double h);
The main function basically creates the matrix, initializes it, then use the values of temperature and field gotten from the terminal and applies $N$ MonteCarlo steps, so it prints the final matrix and finishes the program
int main(int argc, char *argv[]){
srand(1); //seed
//Creation of the matrix
matrix ising;
//Initial conditions, randomly
get_mem(ising);
int Number=Nx*Ny;
double kbT=atof(argv[1])/100;
double h=atof(argv[2])/100;
//Monte Carlo stabilization
for (int t=0;t<N*Number;t++){
Step(ising,kbT,h);
}
//Print the final matrix
print(ising);
return 0;
}
Resize and initialization of the matrix
void get_mem(matrix & m){
m.resize(Nx);
for(int i=0;i<Nx;i++){
m[i].resize(Ny);
for(int j=0;j<Ny;j++){
m[i][j]=2*(rand() % 2)-1;
}
}
}
The function to print the matrix
void print(const matrix & m){
for(int i=0;i<Nx;i++){
for(int j=0;j<Ny;j++){
std::cout<<m[i][j]<<' ';
}
std::cout<<std::endl;
}
}
The metropolis algorithm implemented for this particular case,
void Step(matrix & m,double kbT,double h){
//Auxiliar variables
double x=rand()%Nx; double y=rand()%Ny;
int xp,xl,yp,yl;
double E,Eold,Enew,dE,r;
//Neighbors calculation and
//periodic Boundary conditions applied
xp=x+1;xl=x-1; yp=y+1;yl=y-1;
if(x==0){xl=Nx-1;}
if(x==Nx-1){xp=0;}
if(y==0){yl=Ny-1;}
if(y==Ny-1){yp=0;}
//Calculation of the Energy
E=-J*(m[xl][y]+m[xp][y]+m[x][yl]+m[x][yp])+h;
Eold=m[x][y]*E;
Enew=-m[x][y]*E;
dE=Enew-Eold;
//Metropolis acceptance
if(dE<0){
m[x][y]=-m[x][y];
}else{
r=drand48();
if(r<exp(-dE/kbT)){
m[x][y]=-m[x][y];
}
}
}
The code itself is very simple, but from there we can make great conclusions, here there is a first one of them, which basically shows the transiton, but some things have to me taken into account before making conclusions, this results should be averaged on ensemble averages which we are leaving outside for a while just to see if our code works just as it is. So let us make some remarks,
- Averaging: These quantities are supposed to be averaged on ensemble or if we assume ergodicity, in time.
- Ergodicity: To guarantee ergodicity condition is reached, the system should be on equilibrium, so we are going to do the ensemble averaging.
- Finite size effects: Using periodic boundary conditions do not guarantee the system to be infinite, so there are some size effects that may smooth some behavior we expect to be strong from the theory.
import numpy as np
import os
import matplotlib.pylab as plt
import matplotlib as mpl
mpl.rcParams['mathtext.fontset'] = 'cm'
%config InlineBackend.figure_format = 'retina'
First, let us consider the case of $h=0$, that is, without external magnetic field,
As in the code we are using the temperature as 1/100 the parameter we are about to use, we are going to go from 1 up to 500, so, in our units is from 1 to 5 (We use this range because beforehand we knew that the transition is in the middle).
%%bash
rm *.dat
g++ ising2d.cpp
start=`date +%s`
for t in `seq 1 1 500`
do
./a.out $t 0 > $t.dat
done
end=`date +%s`
echo Time execution: $((end-start)) seg
To make our analysis general, let us select the files we just created, assuming that they will be those the only ones that finishes on .dat, and print the first 50.
for _,_,filenames in os.walk("."):
break
files=[name[:-4] for name in filenames if name.endswith('dat')]
print(files[:50])
Now, we get into each of the files and calculate the average magnetization and the fluctuations which correspond to the magnetic susceptibillity,
M=[]
chi=[]
T=[]
for i in files:
data=np.genfromtxt(i+'.dat')
M.append(np.mean(data))
chi.append((np.var(data))/((float(i)/100)**2))
T.append(float(i)/100)
If we plot the magnetization, we may get something as
plt.plot(T,np.abs(M),'.')
plt.ylabel('Magnetization $|<M>|$')
plt.xlabel('Temperature $k_BT$')
plt.show()
But this do not look anythin as the results we expected, it looks more like a noisi signal. Let us see the susceptibillity,
plt.plot(T,chi,'.')
plt.ylabel('Susceptibillity $\chi$')
plt.xlabel('Temperature $k_BT$')
plt.show()
The first points have a small temperature and therefore, the amplitud increases rapidly, let us take a different rangem
T, chi = (list(l) for l in zip(*sorted(zip(T, chi))))
plt.plot(T[100:],chi[100:],'.')
plt.ylabel('Susceptibillity $\chi$')
plt.xlabel('Temperature $k_BT$')
plt.show()
But this is nothing similar to what we expected to have, let us explore why,
At small temperatures, the system should be all going in the same direction, so let us plot the first temperature we have
data=np.genfromtxt('1.dat')
plt.imshow(data)
And there we see the problem, and it is that clusteres are created, as there is no prefered direction, it is possible to have an stable configuration (which correspond to a local minima) and the montecarlo method may take longer to get to the global minima or simply remain there.
To solve this, let us apply a small perturbation adding an external magnetic field,
%%bash
rm *.dat
g++ ising2d.cpp
start=`date +%s`
for t in `seq 1 1 500`
do
./a.out $t 10 > $t.dat
done
end=`date +%s`
echo Time execution: $((end-start)) seg
We repeat the exact same procedure as before
for _,_,filenames in os.walk("."):
break
files=[name[:-4] for name in filenames if name.endswith('dat')]
print(files[:50])
M=[]
chi=[]
T=[]
for i in files:
data=np.genfromtxt(i+'.dat')
M.append(np.mean(data))
chi.append((np.var(data))/((float(i)/100)**2))
T.append(float(i)/100)
T, chi, M = (list(l) for l in zip(*sorted(zip(T, chi, M))))
So now when we plot the magnetization, we get
plt.plot(T,np.abs(M))
plt.ylabel('Magnetization $|<M>|$')
plt.xlabel('Temperature $k_BT$')
plt.show()
And this is the result we were expecting to have
plt.plot(T,chi)
plt.ylabel('Susceptibillity $\chi$')
plt.xlabel('Temperature $k_BT$')
plt.show()
But we were expecting to have a strong decay and no magnetization after the critical point, as we had in the first plot. So the addition of the external field broke the symetry of the problem but also created an effect we did not want after the simulation, so what we are about to do next is to turn off the external field after the half of the time has passed, so we break the symetry, but afterwards we let the system equilibrate without external field. This will require an addition to one single line in our code,
%%bash
rm *.dat
g++ ising2d_field_off.cpp
start=`date +%s`
for t in `seq 1 1 500`
do
./a.out $t 10 > $t.dat
done
end=`date +%s`
echo Time execution: $((end-start)) seg
The great thing about this, is that the analysis part remain exactly the same, so that is why we decided to do things this way, just with our execution script and the analysis code,
for _,_,filenames in os.walk("."):
break
files=[name[:-4] for name in filenames if name.endswith('dat')]
print(files[:50])
M=[]
chi=[]
T=[]
for i in files:
data=np.genfromtxt(i+'.dat')
M.append(np.mean(data))
chi.append((np.var(data))/((float(i)/100)**2))
T.append(float(i)/100)
T, chi, M = (list(l) for l in zip(*sorted(zip(T, chi, M))))
So when we plot the final results, we get for the magnetization
plt.plot(T,np.abs(M))
plt.ylabel('Magnetization $|<M>|$')
plt.xlabel('Temperature $k_BT$')
plt.show()
And for the susceptibllity,
plt.plot(T,chi)
plt.ylabel('Susceptibillity $\chi$')
plt.xlabel('Temperature $k_BT$')
plt.show()
Which now looks way better to reduyce the noise the system had to be left evolve more time steps, but remmember that we still have no averaged the results, which may solve all these problems.
Fow now, two things are left, which are the averaging and the finite size effects.
As our objective is leave the code as simple as possible, to do the averaging we are going to create a python script to calculate it without affecting our C++, at the end we can construct a .makefile to wrap all these different simulations we have done in this notebook so that everything is always clean and clear for the usage.
So, let us run the code 10 times for each temperature (Of course this can be very optimized reducing the writing part and calculating everything inside the C++ code, but we are just doing a first approach to the problem).
Of course 10 is not enough to be called and ensemble, but just to see how can we do it is enough
%%bash
rm *.dat
g++ ising2d_average.cpp
start=`date +%s`
for t in `seq 1 1 500`
do
for rep in `seq 0 1 9`
do
./a.out $t 10 $rep > $t-$rep.dat
done
done
end=`date +%s`
echo Time execution: $((end-start)) seg
So, the averaging for every temperatures lead to,
%%time
import numpy as np
import os
import matplotlib.pylab as plt
for _,_,filenames in os.walk("."):
break
files=[int(name[:-6]) for name in filenames if (name.endswith('dat') and '-' in (name))]
M=[]
chi=[]
T=[]
for i in np.unique(files):
mag=0
sus=0
temp=0
for j in range(10):
data=np.genfromtxt(str(i)+'-'+str(j)+'.dat')
mag+=(np.abs(np.mean(data)))
sus+=((np.var(data))/((i/100)**2))
M.append(mag/10)
chi.append(sus/10)
T.append(i/100)
T, chi, M = (list(l) for l in zip(*sorted(zip(T, chi, M))))
plt.plot(T,np.abs(M))
plt.ylabel('Magnetization $|<M>|$')
plt.xlabel('Temperature $k_BT$')
plt.show()
plt.plot(T,chi)
plt.ylabel('Susceptibillity $\chi$')
plt.xlabel('Temperature $k_BT$')
plt.show()
Where it is clear that the noise generated is reduced and we can see more clearly the tendencies.