On Jan 6, 2:17 pm, Jan <jan....@gmail.com> wrote:
> I was wondering if anyone had an implementation of the conventional 2D
> Ising Model implemented in IDL that they would share?
>
> Jan
This is a fun model. Try the program below (sorry for the lengthy
listing here, folks). You'll need Craig Markwardt's routine
PLOTIMAGE.PRO to display the evolution of the spin lattice as the
temperature is reduced. It's not a particularly pretty program but it
does some basic things like plot the magnetization and energy as a
function of temperature.
;+
; NAME:
; ISING2D
;
; PURPOSE:
;
; Simple implementation of the 2D Ising model with no applied magnetic
field.
; As the temperature drops, the lattice of spins is plotted. After the
simulation
; has ended, the magnetization and the energy as a function of
temperature is
; plotted. Periodic boundary conditions are implemented.
;
;
; AUTHOR:
;
; Robert M. Dimeo, Ph.D.
; NIST Center for Neutron Research
; 100 Bureau Drive
; Gaithersburg, MD 20899
; Phone: (301) 975-8135
; E-mail: robert.dimeo@nist.gov
; http://www.ncnr.nist.gov/staff/dimeo
;
; CATEGORY:
;
; Physics, cellular automata
;
; CALLING SEQUENCE:
;
; IDL> ising2d
;
; PARAMETERS
;
; The physical parameters of the model can be set in the procedure
named ISING2D below.
; These parameters include the number of temperature steps,
;
; COMMON BLOCKS:
;
; None
;
; REQUIREMENTS:
;
; Uses PLOTIMAGE.PRO from Craig Marwardt (http://www.physics.wisc.edu/
~craigm/idl/)
;
; DISCLAIMER
;
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this disclaimer
; is included unchanged.
;
; MODIFICATION HISTORY:
;
; Written January 29, 2012
; ************************************ ;
function ising2d_linspace,xlo,xhi,nx
compile_opt idl2,hidden
dx = (xhi-xlo)/(nx-1.0)
return,xlo+dx*findgen(nx)
end
; ************************************ ;
function ising2d_lattice,n,random = random
compile_opt idl2,hidden
; Generates a random Ising lattice of size n x n if the keyword RANDOM
is set. Otherwise
; the lattice is a checkerboard of alternating spin-up and spin-down
if keyword_set(random) then lattice = byte(round(randomu(s,n,n)))
if ~keyword_set(random) then begin
lattice = bytarr(n,n)
for j = 0,n-1 do begin
for i = 0,n-1 do begin
if ((i mod 2) eq 1) xor ((j mod 2) eq 1) then lattice[i,j] = 1B
endfor
endfor
endif
return,lattice
end
; ************************************ ;
function ising2d_energy,jcouple,lattice
compile_opt idl2,hidden
; Calculates the energy of the 2D Ising lattice
; inflate the lattice by two elements in each direction
nx = (size(lattice,/dim))[0] & nnx = nx + 2
ll = bytarr(nnx,nnx)
ll[1:nx,1:nx] = lattice
ll[0,1:nx] = lattice[nx-1,0:nx-1]
ll[nnx-1,1:nx] = lattice[0,0:nx-1]
ll[1:nx,0] = lattice[0:nx-1,nx-1]
ll[1:nx,nnx-1] = lattice[0:nx-1,0]
; Muster array mojo to calculate the energy of the lattice in one line
rather
; than a double loop
return,-jcouple*total(((shift(ll,0,1) + shift(ll,0,-1) + shift(ll,1,0)
+ shift(ll,-1,0))*ll)[1:nx,1:nx])
end
; ************************************ ;
function ising2d_magnetism,lattice
compile_opt idl2,hidden
nx = (size(lattice,/dim))[0]
return,total(lattice)/(nx*nx)
end
; ************************************ ;
function ising2d_plot_lattice,lattice
compile_opt idl2,hidden
nx = (size(lattice,/dim))[0]
xr = [0,nx-1]
plotimage,bytscl(lattice),xrange = xr,yrange = xr,imgxrange = xr, $
imgyrange = xr,xstyle = 5,ystyle = 5,xmargin = [0,0],ymargin
= [0,0]
return,1B
end
; ************************************ ;
pro ising2d
; Main program driver and loop through decreasing temperatures
device,decomposed = 0
loadct,1,/silent
; ************ SIMULATION PARAMETERS *******************************
nx = 35 ; side of the ising lattice
jcouple = 1.0 ; spin coupling
nt = 100 ; number of temperature steps
tlo = 1.5 & thi = 5.0 ; lower and upper temperature
nmc = 2 ; number of monte-carlo sweeps/site/temperature
random = 1B ; set to 1B for a random initial array
; set to 0B for alternating spin-up and spin-down
; ************ SIMULATION PARAMETERS *******************************
lattice = ising2d_lattice(nx,random = random) ; creates the spin array
winvis = 0
xsize = (ysize = 400)
window,winvis,xsize = xsize,ysize = ysize, title = '2D Ising Model'
window,/free,xsize = xsize,ysize = ysize,/pixmap
winpix = !d.window
temperature = ising2d_linspace(tlo,thi,nt)
mag = fltarr(nt)
e = fltarr(nt)
for k = 0,nt-1 do begin
t = temperature[nt-1-k] ; decrement the temperature
for ss = 0,nmc-1 do begin
energy = ising2d_energy(jcouple,lattice) ; we need the energy for
comparison in the Metropolis algorithm
; Sweep over the lattice nmc times and flip a spin
testlattice = lattice ; this is the lattice which we'll calculate
the energy with the flipped spins
sto_lattice = lattice ; this is the lattice in which we'll store
the final spin array at each temperature
for j = 0,nx-1 do begin
for i = 0,nx-1 do begin
testlattice = lattice
if lattice[i,j] eq 0B then testlattice[i,j] = 1B else
testlattice[i,j] = 0B
newenergy = ising2d_energy(jcouple,testlattice)
de = newenergy - energy
if (de le 0) then sto_lattice[i,j] = testlattice[i,j]
if ((de gt 0) and (randomu(seed,1) lt exp(-de/t))) then
sto_lattice[i,j] = testlattice[i,j]
endfor
endfor
lattice = sto_lattice
endfor
e[nt-1-k] = ising2d_energy(jcouple,sto_lattice)
mag[nt-1-k] = ising2d_magnetism(sto_lattice)
lattice = sto_lattice
wset,winpix
ret = ising2d_plot_lattice(lattice)
wset,winvis
device,copy = [0,0,!d.x_size,!d.y_size,0,0,winpix]
endfor
;
window,/free,xsize = 500,ysize = 500
plot,temperature,e,psym = -4,ytitle = 'energy',xtitle = 'T'
window,/free,xsize = 500, ysize = 500
plot,temperature,mag,psym = -4,ytitle = 'magnetization',xtitle = 'T'
wdelete,winpix
end
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
![Jan[2] is currently offline](theme/default/images/xoffline.png.pagespeed.ic.XRkd1fkXye.png "Jan[2] is currently offline")
