Part 1¶
a. Write a routine that creates a 2d array and fills it with a random distribution of $\pm 1$.
b. Write a routine to calculate the total energy from a 2d array of spin values.
c. Write a routine that calculates the difference in energy if a single spin is flipped. This can be computed locally.
d. Check that when you use your energy routine to calculate the energy difference that you get the same answer as when you use the energy difference routine.
e. Write a routine that calculates the magnetization.
import numpy as np
#Solution
#create a 2dimensional Lx times Ly spin array and fill with random +-1 values
#does not need to be done efficiently with vectors because only called at beginning
#creates and returns created spin array
def fill_spin(Lx,Ly):
spinarray = np.zeros((Lx,Ly))
for i in range(0,Lx):
for j in range(0,Ly):
x = np.random.random() #[0,1)
if (x >= 0.5):
spinarray[i,j] = 1
else:
spinarray[i,j] = -1
return spinarray
# create an initial condition that is nearly but not completely
# magnetized at the level of frac in [0,1)
# if frac is 0.1 then 90 percent of spins should be 1
def fill_spin_f(Lx,Ly,frac):
spinarray = np.zeros((Lx,Ly))
for i in range(0,Lx):
for j in range(0,Ly):
x = np.random.random() #[0,1)
if (x >= frac):
spinarray[i,j] = 1
else:
spinarray[i,j] = -1
return spinarray
# compute total energy of spin array
# using a Hamiltonian H = sum_nn (- J s_i s_j)
# nn = nearest neighbors (up down and right left)
# where spins are s_i, s_j
# inputs: spinarray, a 2d array of spins
# J is energy of spin interaction
# spins are 1,-1
# spinarray is a previously allocated 2x2 array
# nearest neighbors are 4 adjacent on 2d grid
# using periodic boundary conditions for 2d grid
# returns total summed energy
def energy(spinarray,J):
center = spinarray
up = np.roll(spinarray,-1,axis=0) # Axis shift
right= np.roll(spinarray,-1,axis=1) # Axis shift
sum2 = up + right
sumt = np.sum(sum2*center) # each pair is covered once!
#print(sumt)
E = -J*sumt
return E
# compute change in energy if you flip a single spin
# if you flip the spin at index i,j, compute change in energy
# return new_energy - old_energy
# J is energy of spin interaction
# spin is not actually flipped in spinarray
# i.e. spinarray is not changed
def delta_energy(spinarray,i,j,J):
Lx = spinarray.shape[0]
Ly = spinarray.shape[1]
iL = (i-1)%Lx # periodic boundary conditions
iR = (i+1)%Lx
jL = (j-1)%Ly
jR = (j+1)%Ly
ssum = spinarray[iL,j] + spinarray[iR,j] + spinarray[i,jL] + spinarray[i,jR]
#print('ssum',ssum)
dE = 2*ssum*spinarray[i,j]*J # I have checked signs
# factor of two comes from change of old to new energy
return dE
# return magnitization
def magnetization(spinarray):
n = np.size(spinarray)
return np.sum(spinarray)/n
# we could speed this up by keeping track of the change in magnetization
# rather than than the magnetization itself
# for testing energy, delta energy, everything looks good!
# solution!
J=1
sa = fill_spin(3,4)
#print("M=", magnetization(sa))
print(sa)
E1 = energy(sa,J)
print("E1 ",E1)
i0 = 1
j0 = 2
dE = delta_energy(sa,i0,j0,J)
sa[i0,j0] = -sa[i0,j0]
print(sa)
E2 = energy(sa,J)
print("E2 ",E2)
print("E2-E1 ",E2-E1)
print("dE ", dE) # this should be the same as the above line
print('dE should be the same as E2-E1')
print("M=", magnetization(sa))
print(np.size(sa))
The Metropolis-Hastings algorithm¶
The Metropolis algorithm is the following Monte Carlo method.
a. Randomly select a grid position $i,j$ in two-dimensions.
b. Compute the difference in energy if the spin of the atom at $i,j$ is flipped. If $E_1$ is the original energy and $E_2$ is the energy with that atom's spin flipped, then compute $$dE = E_2 - E_1$$
c. If $dE < 0$ flip the spin. The new configuration has lower energy.
d. If $dE>0$ flip the spin using a probability $$ P = e^{-\beta dE}$$ with $\beta = 1/k_BT$ depending on temperature. The new configuration has higher energy. To do this you can randomly generate a number between [0,1] and flip the spin if your randomly generated number is below probability $P$.
Repeat!
If the temperature is low, the system should become magnetized (all spins are in the same direction) and if the temperature is high the system should stay with magnetization near zero (spins are randomly oriented).
Part 2¶
a. Write routines to implement the Metropolis algorithm.
Plot magnetization as a function of time including steps where the spin flip was not accepted.
b. Show that if $T<1$ (setting $J=1$ and $k_B=1$) the system becomes magnetized.
c. Show that if $T>1$ the system stays with magnetization near zero.
d. Show that two simulations begun at different magnetization converge to the same final magnetization.
# Solution:
# Compute a single Metropolis step in Markov Chain
# randomly choose a grid element i,j
# compute the change in energy caused by a flip in spin
# if change in energy is negative then flip the spin
# if change in energy is positive then flip spin with
# probability P=e^(-beta dE) where dE is change in energy
# and beta is kT
# inputs:
# J is energy of spin interaction
# spinarray is a 2d array
# beta = 1/(k_B T)
# pp if you want printouts
# calls above routine delta_energy
# returns: nothing (but array spinarray could have 1 element flipped)
def metropolis_step(spinarray,J,beta,pp):
Lx = spinarray.shape[0]
Ly = spinarray.shape[1]
# randomly choose a grid position
i = np.random.randint(0,Lx) # return integers in [0,Lx-1] inclusive
j = np.random.randint(0,Ly)
#print(i,j)
dE = delta_energy(spinarray,i,j,J) #compute energy change
isflip=0 # was there a flip
if (dE < 0): # new state has lower energy, accept change
spinarray[i,j] *= -1.0 # flip the spin!
isflip=1
else:
#beta = 1.0/(kB*Temperature)
expfac = np.exp(-1.0*beta*dE) # we need to accept pased on probability of this
x = np.random.random()
if (x < expfac): #accept
spinarray[i,j] *= -1.0 # flip the spin!
isflip=1
if (pp==1):
print(i,j,dE,isflip)
# do N metropolis steps and return a vector
# of magnetizations for each iteration step
kB = 1.0 # Boltzmann constant
def Nmetropolis(spinarray,J,Temperature,N,pp):
beta = 1.0/(kB*Temperature)
magarr=np.zeros(0)
for i in range(0,N):
metropolis_step(spinarray,J,beta,pp)
mag = magnetization(spinarray)
magarr=np.append(magarr,[mag])
return magarr
# short test
sa = np.zeros((5,5)) + 1.0 # magnetized initial condition
mag = magnetization(sa)
print("initial magnetization ", mag)
J=1.0
N=10
Temperature=5 #high temperature
#metropolis_step(spinarray,J,beta,1)
magarr=Nmetropolis(sa,J,Temperature,N,1)
# short test
sa = fill_spin(5,5) # not magnetized initial condition
mag = magnetization(sa)
print("initial magnetization ", mag)
J=1.0
N=10
Temperature=5 #high temperature
#metropolis_step(spinarray,J,Temperature,1)
magarr=Nmetropolis(sa,J,Temperature,N,1)
# short test
sa = fill_spin(5,5) # not magnetized initial condition
mag = magnetization(sa)
print("initial magnetization ", mag)
J=1.0
N=10
Temperature=0.6 #low temperature
#metropolis_step(spinarray,J,Temperature,1)
magarr=Nmetropolis(sa,J,Temperature,N,1)
# short test
sa = np.zeros((5,5)) + 1.0 # magnetized initial condition
mag = magnetization(sa)
print("initial magnetization ", mag)
J=1.0
N=10
Temperature=0.6 #low temperature
#metropolis_step(spinarray,J,beta,1)
magarr=Nmetropolis(sa,J,Temperature,N,1)
print(np.exp(-8/0.6)) # very unlikely to change state
# so this is correct but it will take it forever to get
# out of a completely magnetized state.
# short test
sa = fill_spin_f(5,5,0.1) # nearly magnetized initial condition
mag = magnetization(sa)
print("initial magnetization ", mag)
J=1.0
N=10
Temperature=0.6 #low temperature
#metropolis_step(spinarray,J,beta,1)
magarr=Nmetropolis(sa,J,Temperature,N,1)
%matplotlib inline
from matplotlib import pyplot as plt
# solution
sa = fill_spin(50,50)
mag = magnetization(sa)
print("initial mag ", mag)
J=1
N=100000
Temperature=5.0
magarr=Nmetropolis(sa,J,Temperature,N,0)
sa = np.zeros((50,50)) + 1.0 # a different initial condition
print("initial mag ", mag)
# fully magnetized initiall condition
magarr2=Nmetropolis(sa,J,Temperature,N,0)
plt.plot(magarr)
plt.plot(magarr2)
plt.ylabel("M");
With $T>1$ magnetization stays near zero.
# solution
sa = fill_spin(50,50) # near zero mag initial condition
mag = magnetization(sa)
print("initial mag ", mag)
J=1
N=500000
Temperature=0.6
magarr=Nmetropolis(sa,J,Temperature,N,0)
#print(magarr)
sa = fill_spin_f(50,50,0.1) # magnetized initial condition
print("initial mag ", mag)
magarr2=Nmetropolis(sa,J,Temperature,N,0)
# make sure we plot both with magnetization >0 at the end
if (magarr[-1] < 0):
magarr*= -1
plt.plot(magarr)
plt.plot(magarr2)
plt.ylabel("M");
With $T<1$ magnetization increases! Spins align.