python - Python:TypeError:'int'对象不可下标蒙特卡洛
问题描述
我知道这是一个相当普遍的问题,但查看其他答案并没有帮助我,因为我不完全理解这个问题。我自己得到了实际的错误,但是这个问题只发生在我的代码的一个版本中,其中这行在两个版本中都是相同的,这就是问题所在。
这是第一版:
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
解决方案
您的代码不必要地复杂化了。将固定温度的第一个代码扩展到温度范围非常简单。您所要做的就是将执行晶格初始化、MC 平衡和 MC 生产的代码在一个for
循环中运行,并在给定的温度范围内运行。我创建了列表来存储每个温度下的热容量、气和平均能量。
下面是它是如何完成的。我选择 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*')
推荐阅读
- json - 无法解析通用 Codable
- python - 如何让 PyCharm 分析器只显示我的源代码的时间,而不是任何库?
- javascript - 节点错误“无法读取未定义的属性 'insertId'”
- laravel - 如何从 Vue 组件访问值?
- c# - 无法模拟已经存在的电子邮件功能 - Selenium
- javascript - 编写一个程序来查找由其他单词组成的最长单词,该单词也出现在 javascript 的数组中
- python - 递归计算对
- r - R Shiny:修改数据框中的输出值
- javascript - 如何在组件中使用多个订阅者调用而不覆盖所有订阅者调用的数据
- python - 有在线运行心理脚本的地方吗?