Python：TypeError：'int'对象不可下标蒙特卡洛

Question

我知道这是一个相当常见的问题，但查看其他答案并没有帮助我，因为我不完全理解这个问题。我自己得到了实际的错误，但这个问题只发生在我的代码的一个版本中，这行在这两个版本中都是相同的。

这是第一版：

import numpy as np
import numpy.random as rnd
import matplotlib.pyplot as plt

#S[i, j] = 1 (spin up)
#S[i, j] = -1 (spin down)

def InitSpins(S, N):
    for i in range(N):
        for j in range(N):
            S[i, j] = 2*rnd.randint(2)-1
    return S

def GenerateMove(S, N):
    i = rnd.randint(N)
    j = rnd.randint(N)
    return i, j, S[i][j], - S[i][j]

def ComputeEnergy(S, N):
    E = 0.0
    for i in range(N):
        for j in range(N):
            E += -J*S[i,j]*S[(i+1)%N,j] - J*S[i,j]*S[i,(j+1)%N] - H*S[i,j]
    return E

def ComputeDiffEnergy(S, i, j, old, new, N):
    Eold = ComputeEnergy(S, N)
    S[i,j] = new
    Enew = ComputeEnergy(S, N)
    S[i,j] = old
    return Enew - Eold

def MCStep(S, N, E, nacc, T):
    i,j,old,new = GenerateMove(S, N)
    DE = ComputeDiffEnergy(S, i, j, old, new, N)
    if np.any(DE <= 0) or np.any(rnd.random() < np.exp(-DE/T)):
        AcceptMove(S, i, j, old, new, N)
        E += DE
        nacc += 1
#    elif 
#        AcceptMove(S, i, j, old, new, N)
#        E += DE
#        nacc += 1
    else:
        RejectMove(S, i, j, old, new, N)
    return E, nacc

def AcceptMove(S, i, j, old, new, N):
    S[i, j] = new

def RejectMove(S, i, j, old, new, N):
    S[i, j] = old

def ComputeX(S, N):
   X = 0.0
   for i in range(N):
       for j in range(N):
           X += S[i,j]
   X += 1/(N**2)*np.sum(X)
   return X

N = 10 #dimension of lattice
NIter = 10000 #iterations for production run
NEquil = NIter//10 #iterations in actual calculation
NT = 100 #number of time steps
T = 2.4
H = 0.0 #set outside magnetisation
J = 1.0 #set internal magnetisation
S = np.empty([N, N]) #set initlal spin array

S = InitSpins(S, N)
print(S)
print('energy',ComputeEnergy(S, N))

# Equilibration:
nacc = 0
E = ComputeEnergy(S, N)
for i in range(NEquil):
    E, nacc = MCStep(S, N, E, nacc, T)

# Production run
nacc = 0
sum_E = 0.0
sum_E2 = 0.0
E = ComputeEnergy(S, N)
for i in range(NIter):
    E, nacc = MCStep(S, N, E, nacc, T)
    sum_E += E
    sum_E2 += E**2

X = ComputeX(S, N)

def SumX(X, N):
    sum_X = 0.0
    sum_X2 = 0.0
    for i in range(NIter):
        sum_X += X
        sum_X2 += X**2
    return sum_X, sum_X2

sum_X, sum_X2 = SumX(X, N)

# Calculate averages
av_E = sum_E/float(NIter)
av_E2 = sum_E2/float(NIter)
av_X = sum_X/float(NIter)
av_X2 = sum_X2/float(NIter)
CV = 1/(1*(T**2))*(av_E2-av_E**2)
chi = 1/1*(T**2)*(av_X2-av_X**2)
#M = 

# print results
print("acceptance ratio",nacc/float(NIter))
print("average energy",av_E)
print("heat capacity",CV)
print("magnetic susceptibility", chi)
#print("magnetisation, )

这虽然不完整，但给了我结果。 然后我试图让 python 在T（温度）的范围内运行，给我一张热容量变化的图表。这会破坏我的GenerateMove 函数，并出现return i,j, S[i][j], -S[i][j] 标题中给出的错误。

不起作用的更新代码是：

import numpy as np
import numpy.random as rnd
import matplotlib.pyplot as plt

#S[i, j] = 1 (spin up)
#S[i, j] = -1 (spin down)

def InitSpins(S, N):
    for i in range(N):
        for j in range(N):
            S[i, j] = 2*rnd.randint(2)-1
    return S

def GenerateMove(S, N):
    i = rnd.randint(N)
    j = rnd.randint(N)
    return i, j, S[i][j], - S[i][j]


def ComputeEnergy(S, N):
    E = 0.0
    for i in range(N):
        for j in range(N):
            E += -J*S[i,j]*S[(i+1)%N,j] - J*S[i,j]*S[i,(j+1)%N] - H*S[i,j]
    return E

def ComputeDiffEnergy(S, i, j, old, new, N):
    Eold = ComputeEnergy(S, N)
    S[i,j] = new
    Enew = ComputeEnergy(S, N)
    S[i,j] = old
    return Enew - Eold

def MCStep(S, N, E, nacc, T):
    i,j,old,new = GenerateMove(S, N)
    DE = ComputeDiffEnergy(S, i, j, old, new, N)
    if np.any(DE <= 0) or np.any(rnd.random() < np.exp(-DE/T)):
        AcceptMove(S, i, j, old, new, N)
        E += DE
        nacc += 1
#    elif 
#        AcceptMove(S, i, j, old, new, N)
#        E += DE
#        nacc += 1
    else:
        RejectMove(S, i, j, old, new, N)
    return E, nacc

def AcceptMove(S, i, j, old, new, N):
    S[i, j] = new

def RejectMove(S, i, j, old, new, N):
    S[i, j] = old

def ComputeX(S, N):
   X = 0.0
   for i in range(N):
       for j in range(N):
           X += S[i,j]
   X += 1/(N**2)*np.sum(X)
   return X

N = 10 #dimension of lattice
NIter = 10000 #iterations for production run
NEquil = NIter//10 #iterations in actual calculation
NT = 100 #number of time steps
T = np.linspace(0.1,5.0,NT) #set temperaure range
C = np.zeros(NT) #inital heat capacity
H = 0.0 #set outside magnetisation
J = 1.0 #set internal magnetisation
S = np.empty([N, N]) #set initlal spin array

S = InitSpins(S, N)
print(S)
print('energy',ComputeEnergy(S, N))

## Equilibration:
#nacc = 0
#E = ComputeEnergy(S, N)
#for i in range(NEquil):
#    E, nacc = MCStep(S, N, E, nacc, T)
#
## Production run
#nacc = 0
#sum_E = 0.0
#sum_E2 = 0.0
#E = ComputeEnergy(S, N)
#for i in range(NIter):
#    E, nacc = MCStep(S, N, E, nacc, T)
#    sum_E += E
#    sum_E2 += E**2

X = ComputeX(S, N)

def SumX(X, N):
    sum_X = 0.0
    sum_X2 = 0.0
    for i in range(NIter):
        sum_X += X
        sum_X2 += X**2
    return sum_X, sum_X2

sum_X, sum_X2 = SumX(X, N)

for t in range(NT):
    nacc = 0
    sum_E  = 0.0
    sum_E2 = 0.0
    S = InitSpins(S, N)
    E = ComputeEnergy(S, N)
    for i in range(NIter):
        S = MCStep(S, N, E, nacc, T[t])
        for i in range(NEquil):
            S = MCStep(S, N, E, nacc, T[t])
            E, nacc = MCStep(S, N, E, nacc, T)
            sum_E += E
            sum_E2 += E**2
        av_E = sum_E/float(NIter)
        av_E2 = sum_E2/float(NIter)
        C[t] = (sum_E/(float(NIter*N*N)) - sum_E*sum_E/(NIter*NIter*N*N))/T[t]**2

# Calculate averages
av_X = sum_X/float(NIter)
av_X2 = sum_X2/float(NIter)
CV = 1/(1*(T**2))*(av_E2-av_E**2)
chi = 1/1*(T**2)*(av_X2-av_X**2)
#M = 

# print results
print("acceptance ratio",nacc/float(NIter))
print("average energy",av_E)
print("heat capacity",CV)
print("magnetic susceptibility", chi)
#print("magnetisation, )

#plotting resultts
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(T, C)
plt.show()

编辑：这是完整的回溯

Traceback (most recent call last):

  File "<ipython-input-3-52258536ecbb>", line 1, in <module>
    runfile('F:/adv. numerical project graph.py', wdir='F:')

  File "C:\Users\User\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 704, in runfile
    execfile(filename, namespace)

  File "C:\Users\User\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 108, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "F:/adv. numerical project graph.py", line 121, in <module>
    S = MCStep(S, N, E, nacc, T[t])

  File "F:/adv. numerical project graph.py", line 42, in MCStep
    i,j,old,new = GenerateMove(S, N)

  File "F:/adv. numerical project graph.py", line 24, in GenerateMove
    return i, j, S[i][j], - S[i][j]

IndexError: tuple index out of range

Answer 1

您的代码不必要地复杂化了。将固定温度的第一个代码扩展到温度范围非常简单。您所要做的就是在给定温度范围内运行for 循环中执行晶格初始化、MC 平衡和 MC 生产运行的代码。我创建了列表来存储每个温度下的热容量、气和平均能量。

下面是它是如何完成的。我选择 30 个温度值而不是 100 个，以便快速为您提供有效的解决方案。我只显示相关代码。具有函数定义的其余代码保持不变。

我对你的 MC 预测热容量的行为感到不舒服。它应该随着温度而增加。您需要交叉检查您的实施。你现在拥有一切。

NT = 30 # number of temperature steps
T_list = np.linspace(0.1, 5.0, NT) # set temperaure range
av_E_list, CV_list, chi_list = [], [], []

for T in T_list: # <--- Loop over different temperatures
    print ("Performing MC for T = %.2f" %T)
    S = np.empty([N, N]) #set initlal spin array
    S = InitSpins(S, N)

    # Equilibration:
    nacc = 0
    E = ComputeEnergy(S, N)
    for i in range(NEquil):
        E, nacc = MCStep(S, N, E, nacc, T)

    # Production run
    nacc = 0
    sum_E = 0.0
    sum_E2 = 0.0
    E = ComputeEnergy(S, N)
    for i in range(NIter):
        E, nacc = MCStep(S, N, E, nacc, T)
        sum_E += E
        sum_E2 += E**2

    X = ComputeX(S, N)
    sum_X, sum_X2 = SumX(X, N)

    # Calculate averages
    av_E = sum_E/float(NIter)
    av_E2 = sum_E2/float(NIter)
    av_X = sum_X/float(NIter)
    av_X2 = sum_X2/float(NIter)
    CV = 1/(1*(T**2))*(av_E2-av_E**2)
    chi = 1/1*(T**2)*(av_X2-av_X**2)
    CV_list.append(CV)
    chi_list.append(chi)
    av_E_list.append(av_E)

# Plotting 
plt.plot(T_list, CV_list, '-b*')