Python программа для алгоритма Ван Ландау, остановка после точки - PullRequest
Новая
       83

Python программа для алгоритма Ван Ландау, остановка после точки

1 голос
/ 28 апреля 2020

Я написал python код для выполнения алгоритма Ван Ландау на трехмерных полимерах. Настоящим прилагаются соответствующие коды.

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

Вышеприведенная строка, очевидно, натолкнулась на несколько напечатанных операторов, и через несколько строк она начала выполнять код. Здесь он хотел протестировать код с входными данными. Поэтому я дал необходимые данные, и теперь он работает бесконечно l oop.

Я не могу отследить, какая часть алгоритма отправляет код в бесконечном l oop. Есть ли другой возможный способ решить это? Любая помощь очень ценится.

PS: я хочу, чтобы вы добавили следующие детали, когда вы запускаете код.

Введите число шагов M C: вы можете дать любое целое число, в зависимости от того, что Вы хотите.

Введите число мономеров: 50 (Я написал этот код для первой работы для 50 мономеров полимера)

Введите размеры решетки: 600 (Вы можете использовать любое целое число , но для способа, которым я отправил код, требуется 600. Вы можете использовать любое целое число, если вы прокомментируете строки 17-21 и раскомментируйте строку 16)

1 Ответ

0 голосов
/ 28 апреля 2020

Это довольно сложный код, и я не совсем понимаю, что должны делать разные части, но я почти уверен, что проблема в функции Transform. Если вы добавите счетчик в while l oop, вы увидите, что иногда условие g < 1 будет оставаться истинным, и код никогда не выйдет из l oop. Я заметил, что вы инициализируете, а затем обновляете g1, но вы никогда не используете его. Может быть, вы забыли добавить что-то, что соединяет g и g1?

...