# monte-carlo

By: a guest on Apr 30th, 2012  |  syntax: Python  |  size: 6.26 KB  |  hits: 21  |  expires: Never
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
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.    """