title: How to parallelize Ising Model in Python (multiprocessing package)?
body: <p>I have a Python script that does Monte Carlo simulations of Ising model of a 2D lattice. MC simulations are embarrassingly parallel, sampling for every temperature can be distributed to different threads. I want to use multiprocessing module for this, but new to this package. How can I do this?</p>
  
  <pre><code>from __future__ import division
  import numpy as np
  from numpy.random import rand
  import matplotlib.pyplot as plt
  import multiprocessing as mp
  
  nt    = 5
  N     = 16
  
  Energy = np.zeros(nt)
  Magnetization = np.zeros(nt)
  
  T  = np.linspace(1, 5, nt) 
  
  ## monte carlo moves
  def mcmove(config, beta):
      for i in range(N):
          for j in range(N):
                  a = np.random.randint(0, N)
                  b = np.random.randint(0, N)
                  s =  config[a, b]
                  nb = config[(a+1)%N,b] + config[a,(b+1)%N] + config[(a-1)%N,b] + config[a,(b-1)%N]
                  cost = 2*s*nb
                  if cost &lt; 0:
                      s *= -1
                  elif rand() &lt; np.exp(-cost*beta):
                      s *= -1
                  config[a, b] = s
      return config
  
  # calculate thermodynamic variables
  def calcEnergy(config):
      energy = 0
      for i in range(len(config)):
          for j in range(len(config)):
              S = config[i,j]
              nb = config[(i+1)%N, j] + config[i,(j+1)%N] + config[(i-1)%N, j] + config[i,(j-1)%N]
              energy += -nb*S
      return energy/4.
  
  def calcMag(config):
      mag = np.sum(config)
      return mag
  
  def simulate():
  # --&gt;: parallelize T  
      for m in range(len(T)):
          E1 = M1 = 0
          def initialstate(N):  
              state = 2*np.random.randint(2, size=(N,N))-1
              return state
  
          eqSteps = 2000
          config = initialstate(N)
          for i in range(eqSteps):
              mcmove(config, 1.0/T[m])
  
          mcSteps = 2000
          for i in range(mcSteps):
              mcmove(config, 1.0/T[m])   
              Ene = calcEnergy(config)        
              Mag = calcMag(config)           
  
              E1 = E1 + Ene
              M1 = M1 + Mag
  
              Energy[m]         = E1/(mcSteps*N*N)
              Magnetization[m]  = M1/(mcSteps*N*N)
  
      return Magnetization, Energy
  
  Magnetization, Energy = simulate() 
  print Magnetization, Energy
  </code></pre>