## Source Code for Program "Ising"

For use with An Introduction to Thermal Physics by Daniel V. Schroeder.

A few comments: Unlike the pseudocode in the text, this version of the program runs indefinitely until you press the mouse button. To modify it to run for a fixed number of iterations, remove the "do . . . loop" lines and the "get mouse" instruction, then change the value 100 in the "for" instruction to whatever number you like. Also note that this version prompts you for the lattice size and temperature, so you don't have to change the code every time you run the program. Feel free to change "blue" and "yellow" in the colorsquare subroutine to your two favorite colors. Enjoy!

```program ising          ! Monte Carlo simulation of a 2D Ising model
! using the Metropolis algorithm (importance sampling)

dim s(50,50)                        ! declare array of dipoles
input prompt "Size of lattice (max 50): ": size
input prompt "Temperature (units of epsilon/k): ": T
set window 0,(size+2)*1.55,0,(size+2)*1.1   ! coordinates for graphics
randomize                           ! different random numbers for each run
call initialize(s,size)
do                                  ! main iteration loop
for iteration = 1 to 100
let i = int(rnd*size+1)           ! choose a random row number
let j = int(rnd*size+1)           ! and a random column number
call deltaU(s,i,j,size,Ediff)     ! compute hypothetical delta-U of flip
if Ediff <= 0 then                ! if flipping reduces the energy...
let s(i,j) = -s(i,j)            ! then flip it!
call colorsquare(i,j,s(i,j))
else
if rnd < exp(-Ediff/T) then     ! otherwise the Boltzmann factor
let s(i,j) = -s(i,j)          ! gives the probability of flipping
call colorsquare(i,j,s(i,j))
end if
end if
next iteration
get mouse x,y,mousestate           ! check mouse every 100 iterations
loop until mousestate = 1
end

sub deltaU(s(,),i,j,size,Ediff)       ! compute delta-U of flipping a dipole
! (complicated because of pbc)
if i=1 then let top = s(size,j) else let top = s(i-1,j)
if i=size then let bottom = s(1,j) else let bottom = s(i+1,j)
if j=1 then let left = s(i,size) else let left = s(i,j-1)
if j=size then let right = s(i,1) else let right = s(i,j+1)
let Ediff = 2*s(i,j)*(top+bottom+left+right)
end sub

sub initialize(s(,),size)             ! initialize to random array
for i = 1 to size
for j = 1 to size
if rnd < .5 then let s(i,j) = 1 else let s(i,j) = -1
call colorsquare(i,j,s(i,j))
next j
next i
end sub

sub colorsquare(i,j,s)                ! color square i,j according to s
if s = 1 then set color "blue" else set color "yellow"
box area i,i+1,j,j+1
end sub
```