Véges méretű Ising-modell fajhője¶
hasznos: https://stackoverflow.com/questions/36585517/ising-model-in-python
In [1]:
%pylab inline
import numpy as np # convert list to array
import itertools # for generating n-tuple from a list
import matplotlib.pyplot as plt
from ipywidgets import * # az interaktivitásért felelős csomag
In [2]:
def init_spin_array(rows, cols):
return np.ones((rows, cols))
def find_neighbors(spin_array, lattice, x, y):
left = (x, y - 1)
right = (x, (y + 1) % lattice)
top = (x - 1, y)
bottom = ((x + 1) % lattice, y)
return [spin_array[left[0], left[1]],
spin_array[right[0], right[1]],
spin_array[top[0], top[1]],
spin_array[bottom[0], bottom[1]]]
def energy(spin_array, lattice, x ,y):
return -0.5*spin_array[x, y] * sum(find_neighbors(spin_array, lattice, x, y))
def energy_config(config,lattice):
en_config=[]
for i in range(len(config)):
tmp=0
for x in range(lattice):
for y in range(lattice):
tmp +=energy(config[i],lattice,x,y)
en_config.append(tmp)
return array(en_config)
def Zpart(en_config,T):
return sum(exp(-en_config/T))
def en_expect(en_config,T):
tmp = sum(en_config*exp(-en_config/T))
return tmp/Zpart(en_config,T)
def cV(en_config,T):
tmp1 = sum(en_config*exp(-en_config/T))/Zpart(en_config,T)
tmp2=sum(en_config*en_config*exp(-en_config/T))/Zpart(en_config,T)
return (tmp2-tmp1**2)/T**2
def m_average(config,en_config,T):
m=0
for i in range(len(config)):
m +=sum(config[i])*exp(-en_config[i]/T)
return m/Zpart(en_config,T)
In [3]:
# Az abra kimentesehez az alabbiakat a plt.show() ele kell tenni!!!
#savefig('fig_rainbow_p2_1ray.pdf'); # Abra kimentese
#savefig('fig_rainbow_p2_1ray.eps'); # Abra kimentese
# Abra es fontmeretek
xfig_meret= 9 # 12 nagy abrahoz
yfig_meret= 6 # 12 nagy abrahoz
xyticks_meret= 15 # 20 nagy abrahoz
xylabel_meret= 20 # 30 nagy abrahoz
legend_meret= 17 # 30 nagy abrahoz
size <= 4, Figyelem: ennél nagyobb size-ra nagyon lelassul !!!¶
In [4]:
Tvec = linspace(0.1,5,981)
#Tvec = linspace(0.1,5,10)
In [5]:
size= 2 # size of the sqaure lattice
nspin = size*size # number of spin in the lattice
# all possible configurations in list
# n-Tuples in Mathematica, here n = size
config_tmp = asarray(list(itertools.product(*[(-1, 1)] * nspin)))
config = [reshape(config_tmp[i],(-1, size)) for i in range(len(config_tmp))]
#en_vec=array(energy_config(config,size))
en_config = energy_config(config,size)
#Tvec = linspace(0.1,5,2*981)
# energy as a function of T
ee=[]
[ee.append(en_expect(en_config,tt)) for tt in Tvec];
ee2=ee
# specific heat
#Tvec = linspace(0.1,10,80)
cc=[];
[cc.append(cV(en_config,tt)/size**2) for tt in Tvec];
cc2=cc
In [6]:
size= 3 # size of the sqaure lattice
nspin = size*size # number of spin in the lattice
# all possible configurations in list
# n-Tuples in Mathematica, here n = size
config_tmp = asarray(list(itertools.product(*[(-1, 1)] * nspin)))
config = [reshape(config_tmp[i],(-1, size)) for i in range(len(config_tmp))]
#en_vec=array(energy_config(config,size))
en_config = energy_config(config,size)
#Tvec = linspace(0.1,5,2*981)
# energy as a function of T
ee=[]
[ee.append(en_expect(en_config,tt)) for tt in Tvec];
ee3=ee
# specific heat
#Tvec = linspace(0.1,10,80)
cc=[];
[cc.append(cV(en_config,tt)/size**2) for tt in Tvec];
cc3=cc
In [7]:
#array([sum(config[i]) for i in range(len(config)) ]).T
In [8]:
size= 4 # size of the sqaure lattice
nspin = size*size # number of spin in the lattice
# all possible configurations in list
# n-Tuples in Mathematica, here n = size
config_tmp = asarray(list(itertools.product(*[(-1, 1)] * nspin)))
config = [reshape(config_tmp[i],(-1, size)) for i in range(len(config_tmp))]
#en_vec=array(energy_config(config,size))
en_config = energy_config(config,size)
#Tvec = linspace(0.1,5,2*981)
# energy as a function of T
ee=[]
[ee.append(en_expect(en_config,tt)) for tt in Tvec];
ee4=ee
# specific heat
#Tvec = linspace(0.1,10,80)
cc=[];
[cc.append(cV(en_config,tt)/size**2) for tt in Tvec];
cc4=cc
In [9]:
2/log(1+sqrt(2))
Out[9]:
Energy as a function $T$¶
In [10]:
### plotting E(T)
figsize(xfig_meret,yfig_meret)
plot(Tvec,ee2,'--g',lw=2,label='L=2');
plot(Tvec,ee3,'-.b',lw=2,label='L=3');
plot(Tvec,ee4,color='r',lw=2,label='L=4');
legend(loc='lower right',fontsize=legend_meret)
title("Energy, Ising model on a finite square lattice, $E(T)$",fontsize=15)
ylabel(r'$E(T)$',fontsize=xylabel_meret, rotation=0, labelpad=30)
xlabel(r'$T \, [\frac{J}{k_{\mathrm{B}}}] $',fontsize=xylabel_meret)
#annotate(r'$T_c $',
# xy=(2.5,-7.9), xytext=(2.5, -7.9),fontsize=legend_meret)
#text(4,-10,r'$T_c = \frac{2}{\ln \left(1+\sqrt{2}\right)}\, \, \frac{J}{k_{\mathrm{B}}}\approx 2.27 \, \frac{J}{k_{\mathrm{B}}} $')
#annotate(r'$T_c = \frac{2}{\ln \left(1+\sqrt{2}\right)}\, \, \frac{J}{k_{\mathrm{B}}}\approx 2.27 \, \frac{J}{k_{\mathrm{B}}} $',
# xy=(4,-5), xytext=(4, -5),fontsize=legend_meret)
fifi = r'$T_c = \frac{2}{\ln \left(1+\sqrt{2}\right)}\,' + \
r'\frac{J}{k_{\mathrm{B}}}\approx 2.27 \, \frac{J}{k_{\mathrm{B}}} $'
text(0.2, -22, fifi, fontsize=15)
text(2.4, -34.0, r'$T_c $', fontsize=15)
axvline(2.27, color='g')
#print ("Size of the square lattice = ",size)
#print ("Number of spins= ",nspin)
#print ("Number of configuration = ",2**nspin)
grid();
Specific heat as a function $T$¶
In [11]:
### plotting C_V(T)
figsize(xfig_meret,yfig_meret)
plot(Tvec,cc2,'--g',lw=2,label='L=2');
plot(Tvec,cc3,'-.b',lw=2,label='L=3');
plot(Tvec,cc4,color='r',lw=2,label='L=4');
legend(loc='upper right',fontsize=legend_meret)
title("Specific heat, Ising model on a finite square lattice, $C_V(T)$",fontsize=15)
ylabel(r'$C_V(T)$',fontsize=xylabel_meret, rotation=0, labelpad=40)
xlabel(r'$T \, [\frac{J}{k_{\mathrm{B}}}] $',fontsize=xylabel_meret)
fifi = r'$T_c = \frac{2}{\ln \left(1+\sqrt{2}\right)}\,' + \
r'\frac{J}{k_{\mathrm{B}}}\approx 2.27 \, \frac{J}{k_{\mathrm{B}}} $'
text(2.5, 0.15, fifi, fontsize=15)
#annotate(r'$T_c = \frac{2}{\ln \left(1+\sqrt{2}\right)}\, \, \frac{J}{k_{\mathrm{B}}}\approx 2.27 \, \frac{J}{k_{\mathrm{B}}} $',
# xy=(5,0.15), xytext=(5, 0.25),fontsize=legend_meret)
#annotate(r'$T_c $',
# xy=(2.5,0.01), xytext=(2.5, 0.01),fontsize=legend_meret)
text(2.1, -0.07, r'$T_c $', fontsize=15)
axvline(2.27,ls='-',color='y')
#text = 'T_c approx '
#print ("Size of the square lattice = ",size)
#print ("Number of spins= ",nspin)
print ("Number of configuration for L=4 is ",2**nspin)
print (' L=2, T_c approx = {0:.4f}'.format(Tvec[argmax(cc2)])) ### select the maximum in cc and find the corresponding T in Tvec.
print (' L=3, T_c approx = {0:.4f}'.format(Tvec[argmax(cc3)])) ### select the maximum in cc and find the corresponding T in Tvec.
print (' L=4, T_c approx = {0:.4f}'.format(Tvec[argmax(cc4)])) ### select the maximum in cc and find the corresponding T in Tvec.
grid();
#savefig('/home/cserti/Dropbox/A_Home/Ising_finite_C_V_.pdf',
# pad_inches=0.0,bbox_inches='tight'); # Abra kimentese
In [ ]: