Metropolis Monte Carlo for the Ising Model¶
In this notebook you will do Metropolis Monte Carlo to see the properties of the 1-d Ising model, numerically, and compare to the exact results in 1d
In [1]:
#setup the notebook
%pylab inline
import numpy as np
Part 1, getting the Hamiltonian right¶
Code the energy of the 1D ising model. This will be checked in cases where we know the answer.
$H = \sum_{i=0}^{N-1} -J s_i s_{i+1} - h s_i$
It's better to start the sum at 0 and go up to $N-1$, because of how Python and many other programming languages work.
Also, use periodic boundary conditions, i.e. $s_{N} = s_0$. A good way to do this might be to use modular arithmetic. Check how the % works in python
The state of a lattice will be called configuration. It will be an array of 1's and -1's, so configuration[i]=$s_i$
We will check that the energy function is working by testing the following obvious conditions:
- The energy for all spin down is $-N J + N h = N(h-J)$
- The energy for all spin up is $-N J - N h = -N(h+J)$
- The energy for alternating spin up and half spin down is (if N is even) $+N J$
In [2]:
def energy_ising_1d(configuration,J,h):
num_spins = len(configuration)
energy = 0.0
for i in range(num_spins):
spini = configuration[i]
#set the value of spin i+1, make sure to test if i+1<num_spins, and otherwise account for periodic boundaries
#you can do this with an if statement if you have to
ip1 = (i+1)%num_spins
spinip1 = configuration[ip1]
energy = energy - J * (spini * spinip1) - h*spini
return energy
#Check that the energy is correct
test_num_spins = 10
#this should be true for any J, h
test_J = 1
test_h = 2
test_configuration_1 = -1*np.ones(test_num_spins)
test_configuration_2 = +1*np.ones(test_num_spins)
test_configuration_3 = +1*np.ones(test_num_spins)
#this sets even entries to -1
test_configuration_3[::2] = -1
print("Test Config 1:", test_configuration_1)
print("Energy Config 1:", energy_ising_1d(test_configuration_1,test_J,test_h))
print("Expected Energy Config 1:",test_num_spins*(test_h-test_J))
print()
print("Test Config 2:", test_configuration_2)
print("Energy Config 2:", energy_ising_1d(test_configuration_2,test_J,test_h))
print("Expected Energy Config 2:",-test_num_spins*(test_h+test_J))
print()
print("Test Config 3:", test_configuration_3)
print("Energy Config 3:", energy_ising_1d(test_configuration_3,test_J,test_h))
print("Expected Energy Config 3:",test_num_spins*test_J)
Part 2, simple (slow) Metropolis Monte Carlo¶
We will start with a random configuration and do the Metropolis algorithm This means you have to program the following
1) Compute the energy
While time < N steps:
2) Choose a random spin index in range [0,N-1], using numpy.random.randint(N)
3) Change that spin from 1 to -1, or -1 to 1
4) Generate a random number r in range 0 to 1, using numpy.random.random()
5) a) if r < min(1,np.exp(-beta*(delta E))):
accept move
current_energy = energy after flipping
You can use the python command `pass` to represent not doing anything
b) otherwise, reject move
change that spin back to it's old value
don't update the current energyIn [3]:
#set a seed for the random number generator here. For a given seed, all of your results should be identical
random_seed = 1
np.random.seed(random_seed)
def metropolis_mc_slow(n_steps, n_lattice_sites, beta, J, h, debug=False,save_freq=10):
# we can start with a random configuration of size n_lattice_sites by generating a random list
# of zeros and twos, then subtracting 1, the following does that, do you see why? Play around
# with this function in an empty box if you don't
configuration = 2*np.random.randint(2, size=n_lattice_sites) - 1
average_spins = []
if debug is True:
print("Starting configuration:",configuration)
current_energy = energy_ising_1d(configuration,J,h)
for i in range(n_steps):
spin_to_change = np.random.randint(n_lattice_sites)
# Change configuration[spin_to_change] to it's opposite value (1->-1, -1->1).
# There is a very simple mathematical operation that does this, regarless of it's current value
configuration[spin_to_change] *= -1
energy_flip = energy_ising_1d(configuration,J,h)
r = np.random.random()
#do metropolis test w/ this random nubmer r
if r<min(1,np.exp(-beta*(energy_flip-current_energy))):
current_energy = energy_flip
else:
#set spin back the same way you did before
configuration[spin_to_change] *= -1
#this computes the average of the spin observable
average_spin = configuration.mean()
if i%save_freq == 0:
average_spins.append(average_spin)
if debug and i%10==0:
print("%i: "%i,configuration,"Energy:",current_energy,"Spin:",average_spin)
return average_spins
#do a test high temperature simulation
print("High temperature:")
average_spins = metropolis_mc_slow(n_steps=100, n_lattice_sites=10, beta=0.1, J=1, h=2, debug=True)
#do a test on a low temperature simulation
print("Low temperature:")
average_spins = metropolis_mc_slow(n_steps=100, n_lattice_sites=10, beta=1, J=1, h=2, debug=True)
Part 3, faster metropolis monte carlo¶
Now that we have a working MC code, we can try to speed it up. We can compare the value we get for our slow code to our fast code, to make sure it is working right.
The key to speeding this monte carlo algorithm up a lot, is that the change in energy when we flip a spin only depends on the state of the 2 neighboring spins, we don't have to compute the energy of the whole system. This makes a huge difference if we want to simulate a big lattice.
The new algorithm is exactly the same, but instead of updating the whole energy, we compute the change in energy from flipping, $\Delta E$. If the move is accepted, we update the energy by adding $\Delta E$ to the total energy. In this case, we only have to flip the actual spin value in the configuration array if the move is accepted.
So, if $s_i$ goes from 1->-1, how much does that effect the energy, given values of $s_{i-1}$ and $s_{i+1}$?
If it's working, you should get identical results to the previous test simulation
In [4]:
#set a seed for the random number generator here. For a given seed, all of your results should be identical
random_seed = 1
np.random.seed(random_seed)
def energy_difference(J, h, si, sleft, sright):
#fill in the formula for the energy difference from flipping spin i,
#which has value si= 1 or -1, with spin values sleft and sright on the left and right
dE = 2*h*si + 2*J*si*(sleft+sright)
return dE
def metropolis_mc_fast(n_steps, n_lattice_sites, beta, J, h, debug=False,save_freq=10):
# we can start with a random configuration of size n_lattice_sites by generating a random list
# of zeros and twos, then subtracting 1, the following does that, do you see why? Play around
# with this function in an empty box if you don't
configuration = 2*np.random.randint(2, size=n_lattice_sites) - 1
average_spins = []
if debug is True:
print("Starting configuration:",configuration)
current_energy = energy_ising_1d(configuration,J,h)
for i in range(n_steps):
#set this like you did above:
spin_to_change = np.random.randint(n_lattice_sites)
si = configuration[spin_to_change]
#now figure out the values of the spin to the left and the spin to the right, remembering
# to take into account periodic boundary conditions for both
# you can use if statements if you have to
sright = configuration[(spin_to_change+1)%n_lattice_sites]
sleft = configuration[(spin_to_change-1)%n_lattice_sites]
dE = energy_difference(J,h,si,sleft,sright)
r = np.random.random()
# fill in the metropolis critereon, using dE
if r<min(1,np.exp(-beta*(dE))) :
#flip the spin
configuration[spin_to_change]*=-1
# update the energy
current_energy += dE
else:
pass
#this computes the average spin
average_spin = configuration.mean()
if i%save_freq == 0:
average_spins.append(average_spin)
if debug and i%10==0:
print("%i: "%i,configuration,"Energy:",current_energy,"Spin:",average_spin)
return average_spins
#do a test high temperature simulation
print("High temperature:")
average_spins = metropolis_mc_fast(n_steps=100, n_lattice_sites=10, beta=0.1, J=1, h=2, debug=True)
#do a test on a low temperature simulation
print("Low temperature:")
average_spins = metropolis_mc_fast(n_steps=100, n_lattice_sites=10, beta=1, J=1, h=2, debug=True)
Part 4, Testing for one case¶
Let's see how an mc trajectory looks for $h=0$ for one set of parameters. What is the average spin near? What do you think a histogram of the average spin should look like?
Play around with the paramters to see what happens Increase the number of lattice sites to see what happens to the histogram
In [5]:
def bin_average(bins):
return (bins[1:]+bins[:-1])/2.
test_n_lattice_sites = 100
test_beta = 0.2
test_J = 1
test_h = 0
# we usually do number of mc moves proportional to number of lattice sites.
#this is one reason we had to have fast mc code above
test_n_steps = test_n_lattice_sites*100
# First, let's one run simulation, and see how the values of spin converge
average_spins = metropolis_mc_fast(test_n_steps, test_n_lattice_sites, test_beta, test_J,test_h)
plt.plot(average_spins)
plt.ylabel("$m$")
plt.xlabel("MC Step")
plt.title("$\\beta=%.2f,$J=%.2f,$h=%.2f$"%(test_beta,test_J,test_h))
# it doesn't seem to take long for the average spin value to converage under these conditions,
# but let's compute quantities from the second half of the simulation
#let's look at the histogram of spins from this simulation
spin_hist, spin_bins = np.histogram(average_spins[len(average_spins)//2:] , bins = 20)
plt.figure()
plt.plot(bin_average(spin_bins),spin_hist)
plt.xlabel('$m$')
plt.ylabel("$P(m)$")
#for h in np.arange(-2,2,0.1):
# average_spins = metropolis_mc_fast(test_n_steps, test_n_lattice_sites, test_beta, test_J,h)
# final_average_spin = average
Out[5]:
Part 5, Scanning Parameters¶
Let's see how the average magnetization changes as we lower the temperature (increase $\beta$) or increase $h$. Is this how you expect it to look?
Fill in the formula for the exact result from Tuckerman Equation 16.6.12, so that we can compare to the exact result. Use np.sinh, np.cosh, np.exp,np.sqrt
In [6]:
def ising_spin_exact(beta,J,h):
sbh = np.sinh(beta*h)
cbh = np.cosh(beta*h)
efactor = np.exp(-4*beta*J)
numerator = sbh+sbh*cbh/np.sqrt(sbh*sbh+efactor)
denominator = cbh + np.sqrt(sbh*sbh+efactor)
return numerator/denominator
In [7]:
test_n_lattice_sites = 50
test_beta = 0.75
test_J = 1
# we usually do number of mc moves proportional to number of lattice sites.
#this is one reason we had to have fast mc code above
test_n_steps = test_n_lattice_sites*100
#scanning h
test_h_list = np.arange(-2,2.25,0.25)
spin_vs_h = []
for test_h in test_h_list:
average_spin_at_h = metropolis_mc_fast(test_n_steps, test_n_lattice_sites, test_beta, test_J,test_h)
mean_spin_from_trajectory = np.mean(average_spin_at_h[len(average_spin_at_h)//2:])
spin_vs_h.append(mean_spin_from_trajectory)
predicted_spin_v_h = ising_spin_exact(test_beta,test_J,test_h_list)
p = plt.plot(test_h_list,spin_vs_h,marker='o',label="Simulation",linestyle='')
plt.plot(test_h_list,predicted_spin_v_h,label="Theory",color=p[0].get_color())
plt.xlabel('$h$')
plt.ylabel('$\\langle m \\rangle$')
plt.axhline(0,linestyle='--',color='black')
plt.axvline(0,linestyle='--',color='black')
plt.legend(loc=0)
Out[7]:
In [8]:
test_n_lattice_sites = 50
test_n_steps = test_n_lattice_sites*100
test_J = 1
#scanning beta for a range of h's
test_h_list = np.arange(-1,1.25,0.25)
test_beta_list = np.arange(0,1.25,0.25)
for test_beta in test_beta_list:
spin_vs_h = []
for test_h in test_h_list:
average_spin_at_h = metropolis_mc_fast(test_n_steps, test_n_lattice_sites, test_beta, test_J,test_h)
mean_spin_from_trajectory = np.mean(average_spin_at_h[len(average_spin_at_h)//2:])
spin_vs_h.append(mean_spin_from_trajectory)
predicted_spin_v_h = ising_spin_exact(test_beta,test_J,test_h_list)
p = plt.plot(test_h_list,spin_vs_h,marker='o',label="$\\beta=%.2f$"%test_beta,linestyle='')
plt.plot(test_h_list,predicted_spin_v_h,label="",color=p[0].get_color())
plt.xlabel('$h$')
plt.ylabel('$\\langle m \\rangle$')
plt.axhline(0,linestyle='--',color='black')
plt.axvline(0,linestyle='--',color='black')
plt.legend(loc=0)
Out[8]:
In [9]:
#system size dependence
test_beta = 1.0
test_J = 1
#scanning beta for a range of h's
test_h_list = np.arange(-1,1.25,0.25)
n_lattice_site_list = (3,5,25,50, 100)
for test_n_lattice_sites in n_lattice_site_list:
test_n_steps = test_n_lattice_sites*500
spin_vs_h = []
for test_h in test_h_list:
average_spin_at_h = metropolis_mc_fast(test_n_steps, test_n_lattice_sites, test_beta, test_J,test_h)
mean_spin_from_trajectory = np.mean(average_spin_at_h[len(average_spin_at_h)//2:])
spin_vs_h.append(mean_spin_from_trajectory)
p = plt.plot(test_h_list,spin_vs_h,marker='o',label="$N=%i$"%test_n_lattice_sites,linestyle='')
predicted_spin_v_h = ising_spin_exact(test_beta,test_J,test_h_list)
plt.plot(test_h_list,predicted_spin_v_h,label="Theory,$N\\rightarrow\\infty$",color=p[0].get_color(),linestyle='--')
plt.xlabel('$h$')
plt.ylabel('$\\langle m \\rangle$')
plt.axhline(0,linestyle='--',color='black')
plt.axvline(0,linestyle='--',color='black')
plt.legend(loc=0)
Out[9]:
In [ ]: