Python: TypeError: 'int' object is not subscriptable Monte Carlo

Question

I know this is a fairly common question but looking at other answers for this has not helped me as I do not understand the problem exactly. I get the actual error itself but this problem is only occurring in one of the editions of my code where this line is the same in both that is the issue.

Here is version one:

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, )

This, though incomplete, gives me results. Then I tried to make python run through a range of T (temperature) to give me a graph of the changing heat capacity. This breaks my GenerateMove function with the error given in the title for the return i,j, S[i][j], -S[i][j].

Updated code that doesn't work is:

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()

EDIT: here is the full traceback

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

You have unnecessarily complicated your code. The extension of your first code for a fixed temperature to a range of temperatures is quite straightforward. All you had to do was to club your code performing the lattice initialization, MC equilibration and MC production runs in a for loop running over the given temperature range. I have created lists to store the heat capacity, chi, and average energy at each temperature.

Below is how it's done. I am choosing 30 temperature values instead of 100 so as to quickly give you a working solution. I am only showing the relevant code. Rest of the code with function definitions stay the same.

I am not comfortable with the behavior of your MC predicted heat capacity. It is supposed to increase with temperature. You need to cross check your implementation. You have everything now.

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*')