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
Populating the interactive namespace from numpy and matplotlib

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)
Test Config 1: [-1. -1. -1. -1. -1. -1. -1. -1. -1. -1.]
Energy Config 1: 10.0
Expected Energy Config 1: 10

Test Config 2: [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
Energy Config 2: -30.0
Expected Energy Config 2: -30

Test Config 3: [-1.  1. -1.  1. -1.  1. -1.  1. -1.  1.]
Energy Config 3: 10.0
Expected Energy Config 3: 10

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 energy
In [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)
High temperature:
Starting configuration: [ 1  1 -1 -1  1  1  1  1  1 -1]
0:  [-1  1 -1 -1  1  1  1  1  1 -1] Energy: -6.0 Spin: 0.2
10:  [-1  1  1 -1 -1  1  1  1  1  1] Energy: -10.0 Spin: 0.4
20:  [ 1 -1  1  1  1 -1  1 -1  1  1] Energy: -6.0 Spin: 0.4
30:  [ 1 -1 -1  1  1  1 -1 -1  1  1] Energy: -6.0 Spin: 0.2
40:  [-1 -1 -1 -1  1  1  1 -1 -1  1] Energy: 2.0 Spin: -0.2
50:  [-1  1 -1  1 -1  1 -1 -1 -1 -1] Energy: 10.0 Spin: -0.4
60:  [-1  1  1  1  1  1  1  1  1 -1] Energy: -18.0 Spin: 0.6
70:  [-1 -1 -1 -1 -1 -1  1 -1  1 -1] Energy: 10.0 Spin: -0.6
80:  [ 1  1  1  1 -1 -1  1  1 -1 -1] Energy: -6.0 Spin: 0.2
90:  [-1  1  1 -1  1 -1 -1  1  1 -1] Energy: 2.0 Spin: 0.0
Low temperature:
Starting configuration: [-1 -1  1 -1 -1 -1  1 -1  1  1]
0:  [ 1 -1  1 -1 -1 -1  1 -1  1  1] Energy: 2.0 Spin: 0.0
10:  [ 1 -1  1  1  1  1  1  1  1  1] Energy: -22.0 Spin: 0.8
20:  [1 1 1 1 1 1 1 1 1 1] Energy: -30.0 Spin: 1.0
30:  [1 1 1 1 1 1 1 1 1 1] Energy: -30.0 Spin: 1.0
40:  [1 1 1 1 1 1 1 1 1 1] Energy: -30.0 Spin: 1.0
50:  [1 1 1 1 1 1 1 1 1 1] Energy: -30.0 Spin: 1.0
60:  [1 1 1 1 1 1 1 1 1 1] Energy: -30.0 Spin: 1.0
70:  [1 1 1 1 1 1 1 1 1 1] Energy: -30.0 Spin: 1.0
80:  [1 1 1 1 1 1 1 1 1 1] Energy: -30.0 Spin: 1.0
90:  [1 1 1 1 1 1 1 1 1 1] Energy: -30.0 Spin: 1.0

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)
High temperature:
Starting configuration: [ 1  1 -1 -1  1  1  1  1  1 -1]
0:  [-1  1 -1 -1  1  1  1  1  1 -1] Energy: -6.0 Spin: 0.2
10:  [-1  1  1 -1 -1  1  1  1  1  1] Energy: -10.0 Spin: 0.4
20:  [ 1 -1  1  1  1 -1  1 -1  1  1] Energy: -6.0 Spin: 0.4
30:  [ 1 -1 -1  1  1  1 -1 -1  1  1] Energy: -6.0 Spin: 0.2
40:  [-1 -1 -1 -1  1  1  1 -1 -1  1] Energy: 2.0 Spin: -0.2
50:  [-1  1 -1  1 -1  1 -1 -1 -1 -1] Energy: 10.0 Spin: -0.4
60:  [-1  1  1  1  1  1  1  1  1 -1] Energy: -18.0 Spin: 0.6
70:  [-1 -1 -1 -1 -1 -1  1 -1  1 -1] Energy: 10.0 Spin: -0.6
80:  [ 1  1  1  1 -1 -1  1  1 -1 -1] Energy: -6.0 Spin: 0.2
90:  [-1  1  1 -1  1 -1 -1  1  1 -1] Energy: 2.0 Spin: 0.0
Low temperature:
Starting configuration: [-1 -1  1 -1 -1 -1  1 -1  1  1]
0:  [ 1 -1  1 -1 -1 -1  1 -1  1  1] Energy: 2.0 Spin: 0.0
10:  [ 1 -1  1  1  1  1  1  1  1  1] Energy: -22.0 Spin: 0.8
20:  [1 1 1 1 1 1 1 1 1 1] Energy: -30.0 Spin: 1.0
30:  [1 1 1 1 1 1 1 1 1 1] Energy: -30.0 Spin: 1.0
40:  [1 1 1 1 1 1 1 1 1 1] Energy: -30.0 Spin: 1.0
50:  [1 1 1 1 1 1 1 1 1 1] Energy: -30.0 Spin: 1.0
60:  [1 1 1 1 1 1 1 1 1 1] Energy: -30.0 Spin: 1.0
70:  [1 1 1 1 1 1 1 1 1 1] Energy: -30.0 Spin: 1.0
80:  [1 1 1 1 1 1 1 1 1 1] Energy: -30.0 Spin: 1.0
90:  [1 1 1 1 1 1 1 1 1 1] Energy: -30.0 Spin: 1.0

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]:
Text(0,0.5,'$P(m)$')

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]:
<matplotlib.legend.Legend at 0x10df959e8>
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]:
<matplotlib.legend.Legend at 0x10e0f6828>
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]:
<matplotlib.legend.Legend at 0x112bbd198>
In [ ]: