# Demo entry 6337423

Monte Carlo simulation for 2-D Ising model

Submitted by Tianqing Zhang on Dec 09, 2016 at 09:22
Language: APL. Code size: 3.5 kB.

```import numpy as np
import random

def randomSite(L):
site=[]
for i in range(L**2):
a=random.random()
if a>0.5:
site.append(1)
else:
site.append(-1)
return site

def bondConfiguration(site,L,pfact):
listOfBond=[]
lsqu=L**2
for s in range(len(site)):
if s<lsqu-L:
rand1=random.random()
if site[s]==site[s+L] and rand1<pfact:
listOfBond.append(set([s,s+L]))
else:
rand1=random.random()
if site[s]==site[s-lsqu+L] and rand1<pfact:
listOfBond.append(set([s,s-lsqu+L]))
if (s+1)%L==0:
rand2=random.random()
if site[s]==site[s-L+1] and rand2<pfact :
listOfBond.append(set([s,s-L+1]))
else:
rand2=random.random()
if site[s]==site[s+1] and rand2<pfact:
listOfBond.append(set([s,s+1])  )
return listOfBond

def bondDic(bond):
dic={}
for b in bond:
for s in b:
if s not in dic:
dic[s]=b.difference(set([s]))
else:
dic[s]=dic[s].union(b.difference(set([s])))
return dic

def pickASite(L):
a=(L**2)*random.random()
return int(a)

def bondConnectedWithIt(s,bond):
connect=[]
for b in bond:
if s in b:
for end in b:
if end!=s:
connect.append(end)
return set(connect)

def siteToCluster(s,bond,dic):
cluster=set()
cluster=cluster.union(set([s]))
if bondConnectedWithIt(s, bond)==set([]):
Sc=1
return (Sc,cluster)
clusterNextTime=cluster
while True:
cl1=cluster
for p in clusterNextTime:
for con in dic[p]:
if con not in cluster:
cluster=cluster.union(set([con]))
cl2=cluster
clusterNextTime=cl2.difference(cl1)
if cl1==cl2:
break
Sc=len(cluster)
return (Sc,cluster)

def flipCluster(cluster,site):
for s in cluster:
site[s]=-site[s]

if __name__ == '__main__':
betaJc=0.4407
pfact=1-np.exp(-2*betaJc)
scave=[]
scstd=[]
for L in range(40,55,5):
step=L**2
numberOfSweep=25
for k in range(1):
scsum=0
sweepList=[]
for i in range(numberOfSweep):
site=randomSite(L)
Sum=0
for i in range(step):
bond=bondConfiguration(site, L, pfact)
s=pickASite(L)
dic=bondDic(bond)
(Sc,cluster)=siteToCluster(s,bond,dic)
if i>=step/2:
Sum+=Sc
flipCluster(cluster, site)
Scbar=Sum/float(step/2)
sweepList.append(Scbar)
print Scbar
scbarbar=sum(sweepList)/numberOfSweep
error=np.sqrt(sum([(v-scbarbar)**2 for v in sweepList])/(len(sweepList)**2))
scave.append(scbarbar)
scstd.append(error)
print sweepList
print "Steps in each sweep:",step,".  number of sweep:",numberOfSweep
print "The Sc average for L=",L,"is",scbarbar,"error",error
print scave
print scstd
```

This snippet took 0.00 seconds to highlight.

Back to the Entry List or Home.