In [1]:
%pylab inline
From:
In [3]:
import numpy as np
SIZE=5
STEPS=10
#----------------------------------------------------------------------#
# Check periodic boundary conditions
#----------------------------------------------------------------------#
def bc(i):
return i%SIZE
#----------------------------------------------------------------------#
# Calculate internal energy
#----------------------------------------------------------------------#
def energy(system, N, M):
return -1 * system[N,M] * (system[bc(N-1), M]
+ system[bc(N+1), M]
+ system[N, bc(M-1)]
+ system[N, bc(M+1)])
#----------------------------------------------------------------------#
# Build the system
#----------------------------------------------------------------------#
def build_system():
system = np.random.random_integers(0,1,(SIZE,SIZE))
system[system==0] =- 1
return system
#----------------------------------------------------------------------#
# The Main monte carlo loop
#----------------------------------------------------------------------#
def main(T):
system = build_system()
for step, x in enumerate(range(STEPS)): #list(enumerate(range(STEPS)))
M = np.random.randint(0,SIZE)
N = np.random.randint(0,SIZE)
E = -2. * energy(system, N, M)
if E <= 0.:
system[N,M] *= -1
elif np.exp(-1./T*E) > np.random.rand():
system[N,M] *= -1
print system
#----------------------------------------------------------------------#
# Run the menu for the monte carlo simulation
#----------------------------------------------------------------------#
def run():
print '='*30
print '\tMonte Carlo Statistics for an ising model with'
print '\t\tperiodic boundary conditions'
print '='*30
print "Choose the temperature for your run (0.1-100)"
T = float(raw_input())
main(T)
In [4]:
run()
In [3]:
#!/usr/bin/env python
"""
Monte Carlo simulation of the 2D Ising model
"""
from scipy import *
from scipy import weave
from pylab import *
Nitt = 1000000 # total number of Monte Carlo steps
N = 10 # linear dimension of the lattice, lattice-size= N x N
warm = 1000 # Number of warmup steps
measure=100 # How often to take a measurement
def CEnergy(latt):
"Energy of a 2D Ising lattice at particular configuration"
Ene = 0
for i in range(len(latt)):
for j in range(len(latt)):
S = latt[i,j]
WF = latt[(i+1)%N, j] + latt[i,(j+1)%N] + latt[(i-1)%N,j] + latt[i,(j-1)%N]
Ene += -WF*S # Each neighbor gives energy 1.0
return Ene/2. # Each par counted twice
def RandomL(N):
"Radom lattice, corresponding to infinite temerature"
latt = zeros((N,N), dtype=int)
for i in range(N):
for j in range(N):
latt[i,j] = sign(2*rand()-1)
return latt
def SamplePython(Nitt, latt, PW):
"Monte Carlo sampling for the Ising model in Pythons"
Ene = CEnergy(latt) # Starting energy
Mn=sum(latt) # Starting magnetization
Naver=0 # Measurements
Eaver=0.0
Maver=0.0
N2 = N*N
for itt in range(Nitt):
t = int(rand()*N2)
(i,j) = (t % N, t/N)
S = latt[i,j]
WF = latt[(i+1)%N, j] + latt[i,(j+1)%N] + latt[(i-1)%N,j] + latt[i,(j-1)%N]
P = PW[4+S*WF]
if P>rand(): # flip the spin
latt[i,j] = -S
Ene += 2*S*WF
Mn -= 2*S
if itt>warm and itt%measure==0:
Naver += 1
Eaver += Ene
Maver += Mn
return (Maver/Naver, Eaver/Naver)
def SampleCPP(Nitt, latt, PW, T):
"The same Monte Carlo sampling in C++"
Ene = float(CEnergy(latt)) # Starting energy
Mn = float(sum(latt)) # Starting magnetization
# Measurements
aver = zeros(5,dtype=float) # contains: [Naver, Eaver, Maver]
code="""
using namespace std;
int N2 = N*N;
for (int itt=0; itt<Nitt; itt++){
int t = static_cast<int>(drand48()*N2);
int i = t % N;
int j = t / N;
int S = latt(i,j);
int WF = latt((i+1)%N, j) + latt(i,(j+1)%N) + latt((i-1+N)%N,j) + latt(i,(j-1+N)%N);
double P = PW(4+S*WF);
if (P > drand48()){ // flip the spin
latt(i,j) = -S;
Ene += 2*S*WF;
Mn -= 2*S;
}
if (itt>warm && itt%measure==0){
aver(0) += 1;
aver(1) += Ene;
aver(2) += Mn;
aver(3) += Ene*Ene;
aver(4) += Mn*Mn;
}
}
"""
weave.inline(code, ['Nitt','latt','N','PW','Ene','Mn','warm', 'measure', 'aver'],
type_converters=weave.converters.blitz, compiler = 'gcc')
aE = aver[1]/aver[0]
aM = aver[2]/aver[0]
cv = (aver[3]/aver[0]-(aver[1]/aver[0])**2)/T**2
chi = (aver[4]/aver[0]-(aver[2]/aver[0])**2)/T
return (aM, aE, cv, chi)
if __name__ == '__main__':
latt = RandomL(N)
PW = zeros(9, dtype=float)
wT = linspace(4,0.5,100)
wMag=[]
wEne=[]
wCv=[]
wChi=[]
for T in wT:
# Precomputed exponents
PW[4+4] = exp(-4.*2/T)
PW[4+2] = exp(-2.*2/T)
PW[4+0] = exp(0.*2/T)
PW[4-2] = exp( 2.*2/T)
PW[4-4] = exp( 4.*2/T)
#(maver, eaver) = SamplePython(Nitt, latt, PW)
(aM, aE, cv, chi) = SampleCPP(Nitt, latt, PW, T)
wMag.append( aM/(N*N) )
wEne.append( aE/(N*N) )
wCv.append( cv/(N*N) )
wChi.append( chi/(N*N) )
print T, aM/(N*N), aE/(N*N), cv/(N*N), chi/(N*N)
plot(wT, wEne, label='E(T)')
plot(wT, wCv, label='cv(T)')
plot(wT, wMag, label='M(T)')
xlabel('T')
legend(loc='best')
show()
plot(wT, wChi, label='chi(T)')
xlabel('T')
legend(loc='best')
show()
In [11]:
STEPS=5
list(enumerate(range(STEPS)))
Out[11]:
In [10]:
a.next()
Out[10]:
In [ ]: