CODEWITHSUNDEEP
pythonphysicsmontecarlo
Fp3
Fp3

Reputation: 39

Ising model metropolis algorithm: lattice won't equilibrate

I have some code for the Ising model in python (2d), and the lattice won't reach an equilibrium. Here, the code prints out the number of spins that flip for each Monte Carlo sweep, and the same number are flipping for each sweep. If I'm correct, then the number that flip should decrease with each sweep, as the lattice reaches an 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]

Upvotes: 3

Views: 1289

Answers (1)

rudy
rudy

Reputation: 196

I think the loops 

    for i in range(N):
        for j in range(N):

are not good: you must choose the spin you want to test randomly.

T and other variables are integer and you will have an integer division in the exponential np.exp(-1*delta_E[i,j] / T) which is wrong.

With this code I have a saturation (note I take 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()

Upvotes: 1

Related Questions