我编写了一个python代码来在3D聚合物上执行Wang Landau算法。现附上相关代码
1)原始代码: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()
下面是代码“spinfun.py”,它具有上述代码中使用的函数
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
这就是我现在面临的问题。当我执行代码时,代码运行良好,但在特定步骤突然停止。代码不会返回任何异常或错误,也不会终止,因此我必须终止执行
我搜索答案,并尝试使用以下命令跟踪代码的工作方式
python -m trace --trace WL3D.py
上面这一行显然遇到了一堆打印出来的语句,在几行之后,它开始执行代码。这就是它想要用输入测试代码的地方。所以我给出了必要的输入,它现在运行一个无限循环
我无法追溯到算法的哪一部分在无限循环中发送代码。有没有其他可能的方法来解决这个问题?非常感谢您的帮助
附言:我希望您在运行代码时添加以下详细信息
输入MC步数:您可以根据需要给出任意整数
输入单体的数量:50(我已经写了这个代码来第一次为50个聚合物单体工作)
输入lattice:600的维度(您可以使用任何整数,但对于我发送代码的方式,600是必需的。如果您注释第17-21行并取消注释第16行,则可以使用任何整数)
这是一段相当复杂的代码,我不太明白不同的部分应该做什么,但我很确定问题出在
Transform
函数中。如果向while循环添加一个计数器,您将看到有时g < 1
条件将保持为true,代码将永远不会退出循环。我注意到您初始化并更新了g1
,但从未使用过它。也许你忘了添加连接g
和g1
的东西相关问题 更多 >
编程相关推荐