Don't like ads? PRO users don't see any ads ;-)
Public Pastes
Guest

monte-carlo

By: a guest on Apr 30th, 2012  |  syntax: Python  |  size: 6.26 KB  |  hits: 21  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. #!/usr/bin/env python
  2. import sys
  3. from numpy import array, zeros, exp, arange, sum, ravel
  4. from random import randint, random
  5.  
  6. class Lattice:
  7.     def __init__(self, J=-1.0, S=1.0, L=10, T=1.0):
  8.         self.J = J
  9.         self.S = S
  10.         self.dim = L
  11.         self.lattice = array(map(lambda x: x+1, zeros((self.dim, self.dim))))
  12.         self.kT = 1.0*T
  13.         self.eTot = 0.0
  14.  
  15.         # The possible values for the change in energy when
  16.         # flipping a spin.
  17.         # All up, the middle one is shifted to down: dE = +8.
  18.         # One do, the middle one is shifted to down: dE = +4.
  19.         # [...]
  20.         # All do, the middle one is shifted to down: dE = -8.
  21.         self.dE = [+8, +4, 0, -4, -8]
  22.         #self.dE = [-8, -4, 0, +4, +8]
  23.  
  24.         self.dEB = array(self.dE)
  25.         # Precalculate the corresponding Boltzmann factor
  26.         # for the given temperature.
  27.         self.dEB = exp(-self.dEB/self.kT)
  28.  
  29.         # Convenience lists for the properties to be sampled.
  30.         self.E = []
  31.         self.M = []
  32.  
  33.     # Convenience methods.
  34.     def getLattice(self):
  35.         return self.lattice
  36.  
  37.     def flip(self, i, j):
  38.         self.lattice[i][j] *= -1
  39.  
  40.     def getNeighbours(self, r, c):
  41.         # r: row, c: column
  42.         # (r, c) is the position of the spin of which
  43.         # the positions of the neighbours are returned.
  44.         # Periodic boundaries:
  45.         # if m == 0, then the neighbour is the element
  46.         # in the last row, i.e. the lattice is
  47.         # on a torus.
  48.         up = [(r-1)%self.dim, c]
  49.         le = [r, (c-1)%self.dim]
  50.         ri = [r, (c+1)%self.dim]
  51.         do = [(r+1)%self.dim, c]
  52.         return [up, le, ri, do]
  53.  
  54.     # Called upon initialization of the simulation.
  55.     def calculateInitialTotalEnergy(self):
  56.         eIni = 0.0
  57.         # Run over the lattice and over neighbours of spin at
  58.         # site 'i,j'.
  59.         for i in range(self.dim):
  60.             for j in range(self.dim):
  61.                 for n in self.getNeighbours(i,j):
  62.                      eIni += self.J*self.lattice[i][j]*self.lattice[n[0]][n[1]]
  63.  
  64.         # Devide by 2 to eliminate double counting.
  65.         self.eTot = -eIni/2.0
  66.  
  67.     def calculateEnergyChangeOfFlip(self,i,j):
  68.         # Calculating dE depending on number of neighbours which are
  69.         # already down.
  70.         # By convention, initially all spins up and J>0.
  71.         # 'counter' counts how many spins are of like orientation.
  72.         counter = 0
  73.         for n in self.getNeighbours(i,j):
  74.             # We flipped the i,j spin. Now we check its environment.
  75.             # E.g. for all neighbour spins up, and the i,j spin down,
  76.             # 'if' never evaluates to 'True'. We end up with '+8' in
  77.             # dE.
  78.             if self.getLattice()[n[0]][n[1]] == self.getLattice()[i][j]:
  79.                 counter += 1
  80.         # Returns the 'counter' too so the call to the 'index' function
  81.         # can be avoided.
  82.         return counter, self.dE[counter]
  83.  
  84.     def run(self, step):
  85.         # Define random site to flip.
  86.         k, l = randint(0,self.dim-1), randint(0,self.dim-1)
  87.         self.flip(k, l)
  88.         # Obtain the index of the energy value in the prestored
  89.         # energy array, and the change in energy for the
  90.         # prospective flip.
  91.         ind, dE = self.calculateEnergyChangeOfFlip(k,l)
  92.         if dE <= 0.0:
  93.             if step % 10000 == 0:
  94.                 print "%05i: Accepted. dE = %0.2f" % (step, dE)
  95.             # Only add change in energy to the total energy.
  96.             self.eTot += dE
  97.         else:
  98.             # Uniform distributed random number.
  99.             x = random()
  100.             # Obtain the corresponding, precalculated Boltzmann
  101.             # factor, avoiding the call to the 'index' builtin
  102.             # function.
  103.             if self.dEB[ind] >= x:
  104.                 if step % 10000 == 0:
  105.                     print "%05i: " % step +\
  106.                             "kT: %2.1f, " % self.kT +\
  107.                             "Accepted  Bf: %3.2f, " % self.dEB[ind] +\
  108.                             "%1.4f" % x
  109.                 # Adding the change in energy.
  110.                 self.eTot += dE
  111.             else:
  112.                 # The step is being undone.
  113.                 if step % 10000 == 0:
  114.                     print "%05i: " % step +\
  115.                             "kT: %2.1f, " % self.kT +\
  116.                             "Not Acce. Bf: %3.2f, " % self.dEB[ind] +\
  117.                             "%1.4f" % x
  118.                 # Set the spin state back.
  119.                 self.flip(k,l)
  120.         # Append the total energy of the temperature.
  121.         self.E.append(self.eTot)
  122.         # Append the total magnetization to 'self.M'.  
  123.         self.M.append(sum(self.lattice))
  124.  
  125.     def heat(self, steps):
  126.         E = array(self.E)
  127.         E2 = E*E
  128.         return 1.0/self.dim*1.0/self.dim*\
  129.                 (self.J/self.kT)*(self.J/self.kT)*\
  130.                 (1.0/steps*sum(E2) - (1.0/steps*sum(E))*(1.0/steps*sum(E)))
  131.  
  132.     def magnet(self, steps):
  133.         M = array(self.M)
  134.         M /= (steps*self.dim*self.dim)
  135.         return M
  136.  
  137.     def susz(self, steps):
  138.         M = array(self.M)
  139.         M /= (self.dim*self.dim)
  140.         #M /= M[0]
  141.         M2 = M*M
  142.         return 1.0/self.dim*1.0/self.dim*\
  143.                 (self.J/self.kT)*\
  144.                 (1.0/steps*sum(M2) - (1.0/steps*sum(M))*(1.0/steps*sum(M)))
  145.  
  146. if __name__ == '__main__':
  147.     J = 1.0; S = 1.0; L = 20
  148.     nsteps=20000
  149.     heat = []
  150.     magn = []
  151.     susz = []
  152.     tempRange = arange(1.5, 4.0, 0.1)
  153.     for t in tempRange:
  154.         print "*"*50, t
  155.         lat = Lattice(J, S, L, t)
  156.         lat.calculateInitialTotalEnergy()
  157.         for step in range(nsteps):
  158.             lat.run(step)
  159.         heat.append(lat.heat(nsteps))
  160.         #magn.append(sum(lat.magn(nsteps)))
  161.         magn.append(sum(lat.magnet(nsteps)))
  162.         susz.append(lat.susz(nsteps))
  163.     print heat
  164.  
  165.     import pylab
  166.     print "Simulation done."
  167.     pylab.plot(tempRange, heat, 'bo')
  168.     pylab.savefig('heat.png')
  169.     print "Save heat done."
  170.     pylab.clf()
  171.     pylab.plot(tempRange, magn, 'bo')
  172.     pylab.savefig('magn.png')
  173.     print "Save magn done."
  174.     pylab.clf()
  175.     pylab.plot(tempRange, susz, 'bo')
  176.     pylab.savefig('susz.png')
  177.     print "Save susz done."
  178.     """
  179.    """
create a new version of this paste RAW Paste Data