Ising model of the metropolis: the grid will not balance
I have some code for an Ising model in python (2d) and the lattice will not reach equilibrium. Here, the code prints the number of spins that are flipped for each sweep in Monte Carlo, and the same number is flipped for each sweep. If I'm right, then the number that flips should decrease with each sweep as the grid reaches equilibrium. Can anyone see any errors with the code?
import numpy as np
from numpy import random as rn
N=20
a=[1,-1]
b=[N,N]
#first make an array
init_lattice=rn.choice(a,size=(N,N))
#specify how many Monte Carlo sweeps
number_of_sweeps=1000
#main code to flip spins and record how many are flipped
def new_lattice(lattice,T):
delta_E=np.zeros(b)
switch1=np.zeros(number_of_sweeps)
switch=0
for sweep in range(number_of_sweeps):
for i in range(N):
for j in range(N):
Si=lattice[i,j]
sum_Sj=lattice[i,(j+1)%N]+lattice[(i+1)%N,j]+lattice[i,(j-1)%N]+lattice[(i-1)%N,j]
delta_E[i,j]=2*Si*sum_Sj
if delta_E[i,j]<0:
lattice[i,j]*=-1
switch+=1
elif np.exp(-1*delta_E[i,j]/(T))>rn.random():
lattice[i,j]*=-1
switch+=1
switch1[sweep]=switch
return lattice,switch1
#print how many spins have flipped
switch= new_lattice(init_lattice,5)[1]
print switch
for i in range(len(switch)):
print switch[i+1]-switch[i-1]
+3
source to share
1 answer
I think the loops
for i in range(N):
for j in range(N):
not good: you have to choose a spin that you want to test at random.
T and other variables are integers and you will have integer division in exponential np.exp (-1 * delta_E [i, j] / T), which is not correct.
With this code, I have saturation (note that I accept T = 1. instead of 5. and sweeps = 10000 instead of 1000):
import numpy as np
from numpy import random as rn
import random
import matplotlib.pyplot as plt
latticeDimensions = [20, 20]
nx = latticeDimensions[0]
ny = latticeDimensions[1]
nodesNumber = nx*ny
spinValues = [-1,1]
sweeps = 10000
T = 1.
lattice = init_lattice = rn.choice(spinValues,size=(nx,ny))
#lattice = [ random.choice([-1,1]) for i in range(nodesNumber) ]
t = 0
switch = 0
res = []
while t < sweeps:
selectedNode = random.choice(nodeIdx)
i = random.choice(range(nx))
j = random.choice(range(ny))
Si = lattice[i,j]
sum_Sj=lattice[i,(j+1)%ny] + lattice[(i+1)%nx,j] + lattice[i,(j-1)%ny] + lattice[(i-1)%nx,j]
delta = 2 * Si * sum_Sj
probaTemp = np.exp(- delta / T)
if delta > 0. and rn.random() < probaTemp:
lattice[i,j] *= -1
switch += 1
elif delta <= 0.:
lattice[i,j] *= -1
switch += 1
t += 1
print t
res.append(switch)
plt.plot(res)
plt.show()
+1
source to share