For an animation of the run see: http://imgur.com/a/zpTFJ

Ising2d.py

'''
2d Implementation of the Lenz-Ising Model

Created on May 13, 2014
Updated on Nov 02, 2016

@author: Andrew Bennett

The Ising model can be used to calculate the ground state energy of a spin lattice.  You can think 
of a spin lattice as a grid of particles with a property called spin, which can be plus or minus 1.  
Energy can be stored between two particles, a and b, and is defined by 

E(a,b) = J if spin(a)==spin(b),  
         0 otherwise.

We "solve" the Ising model with the Metropolis algorithm with selection probabilities related to
thermal equilibrium.  First, let dE = E(~a,b) - E(a,b) where ~a refers to a, but with the spin flipped.  
Then, the probability of keeping a random flip is given by 

P(switch) = 1 if dE < 0,
            exp(-b*dE) otherwise

Given this information the algorithm is as follows:
    1. Choose a random lattice site
    2. Calculate the energy associated with it
    3. Calculate the energy of the site if the spin were flipped
    4. If the energy is lowered, keep the flip, otherwise keep the flip with P(switch)
    5. Goto 1
    
Note that this does _not_ simulate the lattice in time, but is rather an optimization to find the ground
state energy of the system.  

See https://en.wikipedia.org/wiki/Ising_model#Monte_Carlo_methods_for_numerical_simulation 
for more information about the implementation.  
'''

# Imports
import math 
import random
import numpy as np
import matplotlib
import matplotlib.pyplot as plot

def neighbors(a,b):
    """
    Nearest neighbors (up, down, left, right)
    """
    return [[a+1,b],[a-1,b],[a,b+1],[a,b-1]]


def E1(x,y):
    """
    Calculates the energy at a site
    """
    energy=0
    for n in neighbors(x,y):
        # If spins are the same energy is added (ie opposites attract)
        energy += J if lat[x][y] == lat[n[0]][n[1]] else 0
    return energy


def E2(x,y):
    """
    Calculates the energy at a site, given a flipped spin
    """
    energy=0
    for n in neighbors(x,y):
        # If spins are the same energy is added (ie opposites attract)
        energy += J if -lat[x][y] == lat[n[0]][n[1]] else 0
    return energy


def base_energy():
    """
    Calculates the energy of the initial lattice
    """
    energy=0
    # Loop over all cells
    for j in range(N-2):
        for k in range(N-2):
            # At each cell grab 4 nearest neighbors
            for n in neighbors(j+1,k+1):
                # Interaction energy defined in E1
                energy+=E1(lat[j+1][k+1],lat[n[0],n[1]])
    return energy


# Variables 
N_draws =  400    # Number of images to create
N_it    =  50     # Number of iterations before drawing
N       =  52     # reassigning from parsed args may be 
J       = -1      # Strength of the interaction
H       =  0      # Background field
beta    =  1e-20  # kT (very small unless ridiculous temperature)

# Initialize the lattice with random spins
lat=2*np.random.randint(2,size=[N,N])-1 

# Calculate the initial energy
energy = np.zeros(N_draws)
energy[1] = base_energy()

# Run the simulation
for i in range(N_draws-1):
    fig = plot.figure()

    # Compute N_it flips
    net_dE = 0  # Reset change in energy
    for t in range(N_it):
        # Choose a random site
        x  = random.randint(0,N-4)+2  # Need to pad the boundaries to account for 
        y  = random.randint(0,N-4)+2  #   different neighbor arrangements
        
        # Compute the energy difference of flipped - current 
        dE = E2(x,y)-E1(x,y)

        # Selection criteria:
        #   - If energy was lowered by flipping, flip
        #   - Otherwise, still flip with probability exp(-b*dE)
        if dE < 0 or np.random.rand() > math.exp(-beta*dE):
            lat[x][y] =- lat[x][y]   # Flip the spin
            net_dE += dE             # Keep track of energy shift

    # Record the energy 
    energy[i+1] = energy[i] + net_dE

    # Plotting and saving
    ax1 = fig.add_subplot(211)    
    ax1.plot(energy)
    ax2 = fig.add_subplot(212)  
    ax2.imshow(lat, interpolation='nearest', cmap='bone')
    plot.draw()
    plot.savefig("out_" +  str(i).zfill(5) + ".png")