# 2d Ising Monte Carlo -- Wolff units: J=1, T<=> kT/J from __future__ import division, print_function from visual import * from visual.graph import * from random import random N=64 #inputs T=2.3 maxsteps=200 N2=N*N #setup print('maxsteps=',maxsteps,'T=',T) s=ones( (N,N), dtype=int ) acount=0 M=N2 print("N=", N, "M(0) =", M) p=1-exp(-2/T) spins=[] #graphics for ij in arange(N2): i,j=divmod(ij,N) c=color.yellow spins=spins+[box(pos=(i-N/2+0.5,j-N/2+0.5,0), length=0.8, height=0.8, width=0.1, color=c)] pm=gdots(size=10) Mave=0. cnave=0. for steps in arange(maxsteps): # computation starts rate(10000) clist=[] #start a new cluster nlist=[] nn=0 cn=0 m=int(N2*random()) #pick a spin i,j=divmod(m,N) s0=s[i,j] #sign of cluster nlist.append([i,j]) nn=nn+1 s[i,j]=-s0 #flip first spin while nn>0: v=nlist.pop() nn=nn-1 i=v[0] j=v[1] clist.append([i,j]) cn=cn+1 #s[i,j]=-s0 ip=mod(i+1,N) #periodic boundary conditions im=mod(i-1,N) jp=mod(j+1,N) jm=mod(j-1,N) if s[ip,j]==s0: u=random() if u maxsteps/4: acount=acount+1 cnave=cnave+cn Mave=Mave+M #end of computation Mave=Mave/acount cnave=cnave/acount print('T M ') print(T,' ',Mave,' ',cnave)