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.

Delete this entry (admin only).