python - Python program for Wang Landau Algorithm halting after a point
问题描述
I wrote a python code to execute Wang Landau algorithm on 3D polymers. Hereby is attached the relevant codes.
1)The Original Code: WL3D.py
from scipy import *
import sys
import numpy as np
from pylab import *
from spinfun import SARW
from spinfun import Energy
from spinfun import Transform
from spinfun import Equality
Niter = int(input("Enter number of MC steps:"))
L=int(input("Enter number of monomers:"))
N=int(input("Enter the dimensions of lattice:"))
spin=zeros((N,N,N), dtype=int)
#1)Initialisation of lattice
#Gx,Gy,Gz=SARW(spin,N)
Gx=[300, 300, 300, 299, 299, 298, 297, 297, 298, 299, 299, 299, 299, 298, 298, 298, 297, 297, 297, 297, 298, 298, 298, 298, 298, 297, 297, 297, 297, 297, 297, 298, 298, 299, 300, 300, 300, 299, 299, 299, 300, 300, 299, 299, 299, 300, 300, 300, 300, 300, 300, 300, 300, 299, 299, 298, 297, 297, 298, 299, 299, 299, 299, 298, 298, 298, 297, 297, 297, 297, 298, 298, 298, 298, 298, 297, 297, 297, 297, 297, 297, 298, 298, 299, 300, 300, 300, 299, 299, 299, 300, 300, 299, 299, 299, 300, 300, 300, 300, 300]
Gy=[300, 301, 302, 302, 301, 301, 301, 300, 300, 300, 300, 300, 300, 300, 301, 302, 302, 302, 302, 302, 302, 302, 302, 301, 300, 300, 300, 301, 301, 301, 300, 300, 301, 301, 301, 302, 303, 303, 302, 302, 302, 301, 301, 301, 302, 302, 301, 300, 300, 300, 300, 301, 302, 302, 301, 301, 301, 300, 300, 300, 300, 300, 300, 300, 301, 302, 302, 302, 302, 302, 302, 302, 302, 301, 300, 300, 300, 301, 301, 301, 300, 300, 301, 301, 301, 302, 303, 303, 302, 302, 302, 301, 301, 301, 302, 302, 301, 300, 300, 300]
Gz=[300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 299, 298, 297, 297, 297, 297, 297, 298, 299, 300, 300, 299, 298, 298, 298, 298, 297, 297, 298, 299, 299, 299, 299, 299, 299, 299, 299, 299, 299, 298, 298, 298, 298, 297, 297, 297, 297, 297, 298, 299, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 299, 298, 297, 297, 297, 297, 297, 298, 299, 300, 300, 299, 298, 298, 298, 298, 297, 297, 298, 299, 299, 299, 299, 299, 299, 299, 299, 299, 299, 298, 298, 298, 298, 297, 297, 297, 297, 297, 298, 299]
for p in range(50):
spin[Gx[p],Gy[p],Gz[p]]=1
E1=Energy(Gx,Gy,Gz,spin)
print "Initial energy of system is",E1
f1=open("Initial_Data.txt",'w')
for i in range(len(Gx)):
f1.write("%d %d %d \n"%(Gx[i],Gy[i],Gz[i]))
f1.close()
#2) Matrices and terms required
Ener=(arange(59)-58).tolist() #List of energy values
E0=-58 #GRound state
lngE=zeros(len(Ener),dtype=float) #LogG
Hist=zeros(len(Ener),dtype=float)
lnf=1.0
#3)The WLA
for time in range(Niter):
i2=0
while i2 < L:
#E1=Energy(Gx,Gy,Gz,spin)
Gx1,Gy1,Gz1=Transform(Gx,Gy,Gz,spin) #Attempting a random move
E2=Energy(Gx1,Gy1,Gz1,spin) #Energy calculation of changed configuration
P = exp(lngE[E1-E0]-lngE[E2-E0]) #Acceptance probability
if P > uniform(0,1):
Gx,Gy,Gz=Equality(Gx1,Gy1,Gz1)
print 'Ok', Gx[0],Gy[0],Gz[0],E1,E2,time,i2
E1=E2
else:
spin[Gx1[0],Gy1[0],Gz1[0]]=0
spin[Gx[49],Gy[49],Gz[49]]=1
print 'Not Ok',Gx[0],Gy[0],Gz[0],E1,E2,time,i2
Hist[E1-E0] += 1.0
lngE[E1-E0] += lnf
i2=i2+1
if time % 1000 == 0:
Ha = sum(Hist)/(len(Hist))
Hmin = min(Hist)
if Hmin > 0.8*Ha:
print time,'Histogram is flat', Hmin, Ha, 'f=',exp(lnf),'lnf=',lnf
Hist=zeros(len(Hist))
lnf=lnf/2.0
if abs(lnf-0.0) < 0.00000001:
break
else:
print time,'Not flat',Hmin, Ha, 'f=',exp(lnf),'lnf=',lnf
f=open("Final_Data.txt",'w')
for i5 in range(len(lngE)):
f.write("%f %f \n"%(Ener[i5],lngE[i5]))
f.close()
Below is the code "spinfun.py" which has the functions used in above code.
from scipy import *
import sys
import numpy as np
from pylab import *
spin=zeros((50,50,50), dtype=int)
Gx2=[25,25,25,24,24,23,22,22,23,24,24,24,24,23,23,23,22,22,22,22,23,23,23,23,23,22,22,22,22,22,22,23,23,24,25,25,25,24,24,24,25,25,24,24,24,25,25,25,25,25]
Gy2=[25,26,27,27,26,26,26,25,25,25,25,25,25,25,26,27,27,27,27,27,27,27,27,26,25,25,25,26,26,26,25,25,26,26,26,27,28,28,27,27,27,26,26,26,27,27,26,25,25,25]
Gz2=[25,25,25,25,25,25,25,25,25,25,24,23,22,22,22,22,22,23,24,25,25,24,23,23,23,23,22,22,23,24,24,24,24,24,24,24,24,24,24,23,23,23,23,22,22,22,22,22,23,24]
#A,B,C=np.loadtxt("Self2_trials.txt",usecols=(0,1,2),unpack=True,dtype=int)
#########
def SARW(spin,N):
g2=1
g1=0
i=N/2;j=N/2;k=N/2
Gx=[N/2];Gy=[N/2];Gz=[N/2]
spin[N/2,N/2,N/2]=1
while g2 < 50:
h=uniform(0,1)
if 0.0 <= h <= 0.167:
k=k+1
elif 0.167 <= h <= 0.334:
j=j+1
elif 0.334 <= h <= 0.501:
i=i+1
elif 0.501 <= h <= 0.668:
i=i-1
elif 0.668 <= h <= 0.835:
j=j-1
else:
k=k-1
if spin[i,j,k] == 0:
print 'Ok',i,j,k,h,g2
spin[i,j,k]=1
Gx=Gx+[i]
Gy=Gy+[j]
Gz=Gz+[k]
g2=g2+1
else:
print 'Not Ok',i,j,k,h
if 0.0 <= h <= 0.167:
k=k-1
elif 0.167 <= h <= 0.334:
j=j-1
elif 0.334 <= h <= 0.501:
i=i-1
elif 0.501 <= h <= 0.668:
i=i+1
elif 0.668 <= h <= 0.835:
j=j+1
else:
k=k+1
g1=g1+1
return Gx,Gy,Gz
############
def Energy(Gx,Gy,Gz,spin):
S=0
Su=[]
for i in range(50):
i1=Gx[i]
j1=Gy[i]
k1=Gz[i]
n=spin[i1+1,j1,k1]+spin[i1-1,j1,k1]+spin[i1,j1+1,k1]+spin[i1,j1-1,k1]+spin[i1,j1,k1+1]+spin[i1,j1,k1-1]
S=S+n
Su=Su+[n]
E1=(S-98)/2
E=-E1
return E
###########
def Transform(Gx,Gy,Gz,spin):
Gx1=[];Gy1=[];Gz1=[]
m=0
while m < 50:
Gx1=Gx1+[Gx[m]]
Gy1=Gy1+[Gy[m]]
Gz1=Gz1+[Gz[m]]
m=m+1
i=Gx1[0];j=Gy1[0];k=Gz1[0]
g=0
g1=0
while g < 1:
h=uniform(0,1)
if 0.0 <= h <= 0.167:
k=k+1
elif 0.167 <= h <= 0.334:
j=j+1
elif 0.334 <= h <= 0.501:
i=i+1
elif 0.501 <= h <= 0.668:
i=i-1
elif 0.668 <= h <= 0.835:
j=j-1
else:
k=k-1
if spin[i,j,k] == 0:
spin[i,j,k]=1
spin[Gx1[49],Gy1[49],Gz1[49]]=0
for i1 in range(50):
n1=49-i1
if n1 != 0:
Gx1[n1]=Gx1[n1-1]
Gy1[n1]=Gy1[n1-1]
Gz1[n1]=Gz1[n1-1]
else:
Gx1[n1]=i;Gy1[n1]=j;Gz1[n1]=k
g=g+1
else:
if 0.0 <= h <= 0.167:
k=k-1
elif 0.167 <= h <= 0.334:
j=j-1
elif 0.334 <= h <= 0.501:
i=i-1
elif 0.501 <= h <= 0.668:
i=i+1
elif 0.668 <= h <= 0.835:
j=j+1
else:
k=k+1
g1=g1+1
return Gx1,Gy1,Gz1
########
def Equality(Gx,Gy,Gz):
Gx1=[];Gy1=[];Gz1=[]
for m in range(50):
Gx1=Gx1+[Gx[m]]
Gy1=Gy1+[Gy[m]]
Gz1=Gz1+[Gz[m]]
return Gx1,Gy1,Gz1
SO here is the problem I m facing right now. When I am executing the code, the code runs fine but it suddenly stops at a particular step. The code doesn't return any exception or error and it doesn't terminate so I have to kill the execution.
I searched for answers, and I tried the following command to trace the workings of the code.
python -m trace --trace WL3D.py
The above line obviously ran into bunch of statements printed, and after some lines it started executing the code. This is where it wanted to test the code with inputs. So I gave the necessary inputs and it is now running an infinite loop.
I am unable to trace back as to which portion of the algorithm is sending the code in an infinite loop. Is there any other possible way to work it out? Any help is much appreciated.
P.S: I want you to add the following details when you are running the code.
Enter number of MC steps: You can give any integer whichever you want.
Enter number of monomers: 50 ( I have written this code to first work for 50 monomers of the polymer)
Enter the dimensions of lattice : 600 (You can use any integer, but for the way I have sent the code, 600 is required. You can use any integer if u comment the lines 17-21 and uncomment line 16)
解决方案
这是一个相当复杂的代码，我并不真正理解不同部分应该做什么，但我很确定问题出在Transform函数中。如果你在 while 循环中添加一个计数器，你会发现有时g < 1条件会一直为真，代码永远不会退出循环。我注意到您初始化然后更新g1，但您从未使用它。也许你忘了添加一些连接g和g1？