Nicholas Miller's ising model code, results, and report.

mill1825/python/grid.py

@@ -0,0 +1,99 @@
+import math
+from math import exp
+import random
+import numpy as np
+
+random.seed()
+
+def pause():
+    garbage = raw_input("[PAUSE]")
+
+class ising_grid(object):    
+    def __init__(self, x_dim=10, y_dim=10):
+        self.grid = []
+        self.x_max = x_dim
+        self.y_max = y_dim
+        self.J = 1.0
+        self.B = 1.0
+        self.T = 1.0
+        self.beta = 1/self.T   
+
+    def initialize(self):
+        self.grid = np.ones([self.y_max, self.x_max])
+
+    def compute_boltzmann_factor(self, S):
+        return (exp(2*(self.J*S + self.B)*self.beta), exp(-2*(self.J*S + self.B)*self.beta))
+
+    def metropolis(self):
+
+        self.initialize()
+
+        boltzmann_factor = {}
+        for S in [-4, -2, 0, 2, 4]: 
+            boltzmann_factor[S] = self.compute_boltzmann_factor(S)
+
+        M = []
+
+        for k in range(self.xM*self.y_max*self.x_max):
+
+            # ------------------------
+            # Metropolis algorithm 
+            # ------------------------
+
+            i = random.randint(0,self.x_max-1)
+            j = random.randint(0,self.y_max-1)
+            val = self.spin(i,j)
+
+            # if spin i is down, choose first column in precomputed boltzmann factors (corresponds to flipping the spin)
+            if(val == -1):
+                w = boltzmann_factor[self.get_S(i, j)][0]
+            # if spin i is down, choose second column in precomputed boltzmann factors (corresponds to flipping the spin)
+            elif(val == 1):
+                w = boltzmann_factor[self.get_S(i, j)][1]
+            if(w > 1):
+                self.flip(i, j)
+            else:
+                if(w > random.random()):
+                    self.flip(i, j)
+
+            M.append(self.get_M())
+
+        M = np.array(M)
+
+        return np.mean(M), self.std_block(M)
+
+    def std_block(self, data):
+        std = []
+        blocks = 50
+        split = len(data)/blocks 
+        for i in range(blocks):
+            std.append(np.std(data[(i*split):((i+1)*split)]))
+        return np.std(np.array(std))
+
+    def set_T(self,T):
+        self.T = float(T)
+        self.beta = 1./self.T
+
+    def spin(self, i, j):
+        return self.grid[j][i]    
+
+    def get_M(self):
+        return np.sum(self.grid)/self.x_max/self.y_max
+
+    def flip(self, i, j):
+        self.grid[j][i] *= -1
+
+    def get_E(self):
+        E = 0.0
+        for j in range(self.y_max):
+            for i in range(self.x_max):
+                E += self.get_S(i, j) 
+        return E
+
+    def get_S(self, i_in, j_in, p=False):
+        i, j = i_in%self.x_max, j_in%self.y_max
+        # This is where the periodic BCs are implemented
+        ip1, jp1 = (i_in+1)%self.x_max, (j_in+1)%self.y_max
+        im1, jm1 = (i_in-1)%self.x_max, (j_in-1)%self.y_max
+        return (self.grid[jm1][i] + self.grid[jp1][i] + self.grid[j][ip1] + self.grid[j][im1])
+

mill1825/python/gui.py

@@ -0,0 +1,90 @@
+from Tkinter import *
+import colorsys
+
+class ising_gui(object):
+    MIN_RED = MIN_GREEN = MIN_BLUE = 0x0
+    MAX_RED = MAX_GREEN = MAX_BLUE = 0xFF    
+    PIXEL_WIDTH = 10
+
+    def __init__(self, grid_x, grid_y):
+        self.grid_x = grid_x
+        self.grid_y = grid_y
+        self.pixels = [(0, 0, 0) for i in range(self.grid_x * self.grid_y)]
+
+        self.root = Tk()
+        self.root.title("Ising Simulation %d x %d" % (self.grid_x, self.grid_y))
+        self.root.resizable(0, 0)
+        self.frame = Frame(self.root, bd=5, relief=SUNKEN)
+        self.frame.pack(side = TOP)
+
+        self.canvas = Canvas(self.frame,
+                             width=self.PIXEL_WIDTH * self.grid_x,
+                             height=self.PIXEL_WIDTH * self.grid_y,
+                             bd=0, highlightthickness=0)
+        self.canvas.pack()
+
+
+        self.TString = StringVar()
+        self.EString = StringVar()
+        self.GString = StringVar()    
+
+        self.bframe = Frame(self.root)
+        self.TLabel = Label(self.bframe, anchor=N, textvariable = self.TString)
+        self.ELabel = Label(self.bframe, anchor=S, textvariable = self.EString)
+        self.GLabel = Label(self.bframe, anchor=CENTER, textvariable = self.GString)   
+        self.bframe.pack(side = BOTTOM)          
+        self.root.update()
+
+    def draw_grid(self, grid):
+        self.set_T(grid.T)
+
+        for j in range(grid.y_max):
+            for i in range(grid.x_max):
+                spin = grid.spin(i, j)
+                if(spin == 1):
+                    self.set_pixel(i, j, (0, 0, 0x00))
+                elif(spin == -1):
+                    self.set_pixel(i, j, (0, 0, 0xFF))
+                else:
+                    self.set_pixel(i, j, (0x00, 0xFF, 0xFF))
+
+        self.draw()
+
+    def set_pixel(self, x, y, hsv):
+        self.pixels[self.grid_x * y + x] = hsv
+
+    def get_pixel(self, x, y):
+        return self.pixels[self.grid_x * y + x]
+
+    def set_T(self, temp):
+        self.TString.set("T = %f" % temp)
+        self.TLabel.pack()
+
+    def draw(self):
+        self.canvas.delete(ALL)
+        for x in range(len(self.pixels)):
+            x_0 = (x % self.grid_x) * self.PIXEL_WIDTH
+            y_0 = (x / self.grid_x) * self.PIXEL_WIDTH
+            x_1 = x_0 + self.PIXEL_WIDTH
+            y_1 = y_0 + self.PIXEL_WIDTH
+            hue = "#%02x%02x%02x" % self._get_rgb(self.pixels[x])
+            self.canvas.create_rectangle(x_0, y_0, x_1, y_1, fill=hue)
+        self.canvas.update()
+
+    def clear(self):
+        for i in range(self.width * self.height):
+            self.pixels[i] = (0, 0, 0)
+
+    def _hsv_to_rgb(self, hsv):
+        rgb = colorsys.hsv_to_rgb(*hsv)
+        red = self.MAX_RED * rgb[0]
+        green = self.MAX_GREEN * rgb[1]
+        blue = self.MAX_BLUE * rgb[2]
+        return (red, green, blue)
+
+    def _get_rgb(self, hsv):
+        red, green, blue = self._hsv_to_rgb(hsv)
+        red = int(float(red) / (self.MAX_RED - self.MIN_RED) * 0xFF)
+        green = int(float(green) / (self.MAX_GREEN - self.MIN_GREEN) * 0xFF)
+        blue = int(float(blue) / (self.MAX_BLUE - self.MIN_BLUE) * 0xFF)
+        return (red, green, blue)

mill1825/python/ising.py

@@ -0,0 +1,47 @@
+def pause():
+    garbage = raw_input("[PAUSE]")
+
+from grid import *
+from gui import ising_gui
+import numpy as np
+import matplotlib.pyplot as plt
+
+gui = True
+
+# ---------------------------------------
+# Change lattice dimensions here
+Ni = 20
+Nj = 20
+# ---------------------------------------
+
+grid = ising_grid(Ni, Nj)
+if gui: display = ising_gui(Ni, Nj)    
+
+M, M_err = [], []
+Ch, ChB = [], []
+T = np.linspace(1.0, 4.0, 30)
+
+grid.B = 0.00
+grid.J = 1.00
+
+# ---------------------------------------
+# Change Metropolis iteration scale here
+grid.xM = 50
+# ---------------------------------------
+
+for T_i in T:
+    grid.set_T(T_i)
+    M_i, M_std = grid.metropolis()
+    if gui: display.draw_grid(grid)
+    print "T = %.3e, M = %.3e"%(T_i, M_i)
+    M.append(M_i)
+    M_err.append(M_std)
+
+M = np.array(M)
+fig = plt.figure()
+ax_1 = fig.add_subplot(111)
+ax_1.errorbar(T, M, yerr=M_err, fmt="b.")
+plt.xlabel("Temperature (K)")
+plt.ylabel("Magnetization")
+
+plt.show()

mill1825/report/ising.aux

@@ -0,0 +1,16 @@
+\relax 
+\@writefile{toc}{\contentsline {section}{\numberline {1}Theory}{1}}
+\@writefile{toc}{\contentsline {subsection}{\numberline {1.1}Ising Model}{1}}
+\newlabel{eq:spins}{{1}{1}}
+\@writefile{toc}{\contentsline {subsection}{\numberline {1.2}Monte Carlo and the Metropolis Algorithm}{1}}
+\newlabel{eq:mag_avg}{{3}{1}}
+\newlabel{eq:magnetization}{{4}{2}}
+\@writefile{toc}{\contentsline {section}{\numberline {2}Some Coding Tricks}{2}}
+\@writefile{toc}{\contentsline {subsection}{\numberline {2.1}Precomputing Boltzmann Factors}{2}}
+\@writefile{toc}{\contentsline {subsection}{\numberline {2.2}Periodic Boundary Conditions}{2}}
+\@writefile{toc}{\contentsline {section}{\numberline {3}Results}{2}}
+\@writefile{lof}{\contentsline {figure}{\numberline {1}{\ignorespaces Magnetization over temperature of a $10\times 10$ periodic lattice with a $N_x=10,000$ scale to the number of Metropolis iterations.\relax }}{3}}
+\providecommand*\caption@xref[2]{\@setref\relax\@undefined{#1}}
+\newlabel{fig:mag10x10}{{1}{3}}
+\@writefile{lof}{\contentsline {figure}{\numberline {2}{\ignorespaces Final configuration for three different temperatures.\relax }}{3}}
+\newlabel{fig:mag50x50}{{2}{3}}