title: Ising Model 2D Python
body: <p>I've coded out the Ising Model simulation based on Monte Carlo simulation and the Metropolis Algorithm but I am having some trouble, namely the <code>ValueError: setting an array element with a sequence</code>. Any help to solve this would be appreciated.</p>
  
  <p>The following shows where the error lies in the code:</p>
  
  <p><img src="https://i.stack.imgur.com/jCLHO.png" alt="The image embedded shows where the error lies in the code">.</p>
  
  <p>If anyone spots any other errors that may inhibit the code from producing the graphs then please let me know.</p>
  
  <pre><code># Code Mark 1.4: Ising Model in 2D
  
  # Importing all the necessary packages I will use:
  import numpy as np
  import scipy as sp
  import matplotlib.pyplot as plt
  import numpy.random as rnd
  
  # Before I do anything I need to make sure the specifics of the size of the lattice as a 12x12 sized, number of steps, number of temperature points used, number of points in order to reach an equilibrium state are all accounted for
  num=16 
  steps_monte=100
  temp_points=50 
  steps_equil=100 
  num_a=1/(steps_monte*(num**2))
  num_b=1/((steps_monte**2)*(num**2))
  
  # Next I must create an intitial spin state for my molecule: 
  def startspin(num):   
      spin=rnd.randint(2,size=(num,num))-1
      return spin
  
  # Now I need to define the energy function. This will give the energy of any given spin state:
  def Energy(Q):
      starting_energy=0
      for i in range(len(Q)):
          for j in range(len(Q)):
              g=Q[i,j]
              n_y=Q[(i+1)%num,j]+Q[i,(j+1)%num]+Q[(i-1)%num,j]+Q[i,(j-1)%num]
              starting_energy+=g*-n_y
      return starting_energy/4
  
  # Now I need to code the Monte Carlo Steps with Metro ALgorithm for the ising model:
  def montestep(Q,P):
      for i in range(num):
          for j in range(num):
                  x=rnd.randint(0,num)
                  y=rnd.randint(0,num)
                  z=Q[x,y]
                  n_y=Q[(x+1)%num,y]+Q[x,(y+1)%num]+Q[(x-1)%num,y]+Q[x,(y-1)%num]
                  l=2*z*n_y
                  if l&lt;0:
                      z*=-1
                  elif rnd.rand()&lt;np.exp(P*-l):
                      z*=-1
                  Q[x,y]=z
      return Q
  
  #Now I need to define the magnetization function. This gives the Magnetization of any given spin state:
  def Magnetization(Q):
      magnetization=np.sum(Q)
      return magnetization
  
  # Here I am saying that the temperature points need to have a midpoint where we expect the disorder to begin and to give us random points between this value as 1&lt;T&lt;5
  temp_mid = 2.25;    
  Temp=rnd.normal(temp_mid,0.5,temp_points)
  Temp=Temp[(Temp&gt;1)&amp;(Temp&lt;5)]   
  temp_points=np.size(Temp)
  
  # I need to define 0 points for each of the 4 physical quantities I wish to find:
  Magnetization=np.zeros(temp_points)
  SpecificHeat=np.zeros(temp_points)  
  Energy=np.zeros(temp_points)
  Susceptibility=np.zeros(temp_points)
  
  # Now I need to implement each function above by creating an Ising Function with them:
  for j in range(len(Temp)):
      E_a=M_a=0
      E_b=M_b=0
      Q=startspin(num)
      init_Temp_a=1/Temp[j] 
      init_Temp_b=init_Temp_a**2
  
      for i in range(steps_monte):
          montestep(Q,init_Temp_a)            
          Energy_calc=Energy[Q]  
          Mag_calc=Magnetization[Q]       
          E_a=E_a+(Energy_calc)
          M_a=M_a+(Mag_calc)
          M_b=M_b+(Mag_calc**2)
          E_b=E_b+(Energy_calc**2)
          Energy[j]=num_a*E_a[i][j]
          Magnetization[j]=num_a*M_a[i][j]
          SpecificHeat[j]=(num_a*E_b[i][j]-num_b*(E_a[i][j]**2))*init_Temp_b
          Susceptibility[j]=(num_a*M_b[i][j]-num_b*(M_a[i][j]**2))*init_Temp_a
  
      for k in range(steps_equil):        
          montestep(Q,init_Temp_a)   
  
  #Finally I can plot the 4 graphs I need to show the monte carlo steps of the ising model and can then conduct my research thereafter.
  
  plt.plot(T, Energy,'d', color="#8A2BE2")
  plt.xlabel("Temperature (Kelvin)")
  plt.ylabel("Energy (Arbitrary Units)")
  plt.show()
  
  plt.plot(T, abs(Magnetization),'x', color="#7FFF00")
  plt.xlabel("Temperature (Kelvin)")
  plt.ylabel("Magnetization (Arbitrary Units)")
  plt.show()
  
  plt.plot(T, SpecificHeat, 'd', color="#00FFFF")
  plt.xlabel("Temperature (Kelvin)")
  plt.ylabel("Specific Heat (Arbitrary Units)")
  plt.show()
  
  plt.plot(T, Susceptibility, 'x', color="#0000FF")
  plt.xlabel("Temperature (Kelvin)")
  plt.ylabel("Susceptibility (Arbitrary Units)")
  plt.show()
  </code></pre>