! Ising model in two dimensions on a square lattice program ising option nolet dim spin(0,0) randomize call initialize(spin,nmax,t,h,flag,of$) clear set cursor 1,8 print "edge ="; nmax;", T =";t;"(J/k), H =";h;"(J/mu)" set cursor 1,60 print "MCS = " set cursor 1,80 print "Type p to pause/q to quit" open #1: screen 0.05,0.95,0.08,0.98 open #2: screen 0.05,0.45,0.025,0.05 open #3: screen 0.55,0.95,0.025,0.05 window #3 emax = 2+h set window -emax-0.01,emax+0.01,-0.105,0.105 set background color "white" clear set color "black" box lines -emax-0.01,emax+0.01,-0.105,0.105 window #2 set window -1.005,1.005,-0.105,0.105 set background color "white" clear set color "black" box lines -1.005,1.005,-0.105,0.105 window #1 set window 0,nmax+1,0,nmax+1 set background color "white" clear call init_spin(spin,nmax,flag,h,e,m) i=0 if of$ <> "" then open #4: name of$, organization text, create newold erase #4 print #4, using "######": i; print #4, using "<#":", "; print #4, using "#.#####^^^^": m/nmax^2; print #4, using "<#":", "; print #4, using "#.#####^^^^": e/nmax^2 end if do window #0 set cursor 1,65 print i i=i+1 window #2 set color "white" box area -1,1,-0.1,0.1 set color "red" box area 0,m/nmax^2,-0.1,0.1 window #3 set color "white" box area -emax,emax,-0.1,0.1 set color "green" box area 0,e/nmax^2,-0.1,0.1 window #1 call calculate(spin,nmax,t,h,e,m) if of$ <> "" then print #4, using "######": i; print #4, using "<#":", "; print #4, using "#.#####^^^^": m/nmax^2; print #4, using "<#":", "; print #4, using "#.#####^^^^": e/nmax^2 end if if key input then get key z c$ = chr$(z) if c$="s" or c$="p" then get key z else if c$="q" or c$="x" then exit do end if end if loop close #1 close #2 close #3 if of$ <> "" then close #4 end ! initialize variables ! nmax = side of lattice sub initialize(spin(,),nmax,t,h,flag,of$) input prompt "lattice size nmax -> ": nmax ! lattice size is nmax by nmax mat redim spin(nmax,nmax) input prompt "temperature (in units of J/k) -> ": t input prompt "field (in units of J/mu) -> ": h input prompt "initial spin config (1 up, -1 down, 2 random) -> ": flag line input prompt "output file for magnetization/energy time series -> ": of$ end sub ! initialize spin directions ! also initialize magnetization and energy sub init_spin(spin(,),nmax,flag,h,e,m) set window 0.5,nmax+0.5,0.5,nmax+0.5 set background color "white" clear for i = 1 to nmax for j = 1 to nmax if flag = 1 then spin(i,j) = 1 else if flag = -1 then spin(i,j) = -1 else ! else start with a random array if(rnd >= 0.5) then spin(i,j) = 1 else spin(i,j) = -1 end if end if if spin(i,j) = 1 then set color "black" box area i-0.1,i+0.1,j-0.1,j+0.1 end if next j next i m = 0 e = 0 for i = 1 to nmax for j = 1 to nmax m = m + spin(i,j) call neighbors(spin,nmax,i,j,sum) e = e - spin(i,j) * sum / 2 - h * spin(i,j) next j next i end sub ! do the real work here at temperature t and field h sub calculate(spin(,),nmax,t,h,e,m) for j = 1 to nmax ! sweep along each row and column for k = 1 to nmax call neighbors(spin,nmax,j,k,sum) !check each spin for possible flip call test_for_flip(spin,j,k,t,h,sum,nmax,e,m) next k next j end sub ! find the net alignment of the nearest neighbors of spin j,k ! note periodic boundary conditions sub neighbors(spin(,),nmax,j,k,sum) x1 = j y1 = k + 1 call test_for_edge(y1,nmax) ! have to deal with edges carefully x2 = j y2 = k - 1 call test_for_edge(y2,nmax) x3 = j + 1 y3 = k call test_for_edge(x3,nmax) x4 = j - 1 y4 = k call test_for_edge(x4,nmax) sum = spin(x1,y1) + spin(x2,y2) + spin(x3,y3) + spin(x4,y4) end sub sub test_for_edge(x,nmax) if x < 1 then x = nmax if x > nmax then x = 1 end sub ! apply the Monte-Carlo flipping rules here sub test_for_flip(spin(,),j,k,t,h,sum,nmax,e,m) denergy = 2 * spin(j,k) * sum + 2 * spin(j,k) * h if denergy < 0 then ! always flip if it lowers the energy spin(j,k) = - spin(j,k) if spin(j,k) = -1 then set color "blue" box area j-0.1,j+0.1,k-0.1,k+0.1 else set color "red" box area j-0.1,j+0.1,k-0.1,k+0.1 end if m = m + 2*spin(j,k) e = e + denergy else ! if energy would increase, flip with probability if exp(-denergy / t) > rnd then ! equal to Boltzmann factor spin(j,k) = - spin(j,k) if spin(j,k) = -1 then set color "blue" box area j-0.1,j+0.1,k-0.1,k+0.1 else set color "red" box area j-0.1,j+0.1,k-0.1,k+0.1 end if m = m + 2*spin(j,k) e = e + denergy end if end if end sub