python ising magnetisation and energy

1import numpy as np
2import random
3
4#creating the initial array
5def init_spin_array(rows, cols):
6    return np.ones((rows, cols))
7
8#calcuating the nearest neighbours
9def find_neighbors(spin_array, lattice, x, y):
10    left   = (x, y - 1)
11    right  = (x, (y + 1) % lattice)
12    top    = (x - 1, y)
13    bottom = ((x + 1) % lattice, y)
14
15    return [spin_array[left[0], left[1]],
16            spin_array[right[0], right[1]],
17            spin_array[top[0], top[1]],
18            spin_array[bottom[0], bottom[1]]]
19
20#calculating the energy of the configuration
21def energy(spin_array, lattice, x ,y):
22    return 2 * spin_array[x, y] * sum(find_neighbors(spin_array, lattice, x, y))
23
24
25#main code
26def main10():
27    #defining the number of initial sweeps, the lattice size, and number of monte carlo sweeps
28    RELAX_SWEEPS = 50
29    lattice = 10
30    sweeps = 1000
31    e1= e0 = 0
32    for temperature in np.arange(0.1, 4.0, 0.2):
33        #setting up initial variables
34        spin_array = init_spin_array(lattice, lattice)
35        mag = np.zeros(sweeps + RELAX_SWEEPS)
36        spec = np.zeros(sweeps + RELAX_SWEEPS)
37        Energy = np.zeros(sweeps + RELAX_SWEEPS)
38        # the Monte Carlo
39        for sweep in range(sweeps + RELAX_SWEEPS):
40            for i in range(lattice):
41                for j in range(lattice):
42                    e = energy(spin_array, lattice, i, j)
43                    if e <= 0:
44                        spin_array[i, j] *= -1
45                    elif np.exp((-1.0 * e)/temperature) > random.random():
46                        spin_array[i, j] *= -1
47
48            #Thermodynamic Variables 
49
50            #Magnetization
51            mag[sweep] = abs(sum(sum(spin_array))) / (lattice ** 2)
52
53            #Energy
54            Energy[sweep] = energy(spin_array,lattice,i,j)/ (lattice ** 2)
55
56            #Specific Heat
57            e0 = e0 + energy(spin_array,lattice,i,j)               
58            e1 = e1 + energy(spin_array,lattice,i,j) *energy(spin_array,lattice,i,j)
59            spec[sweep]=((e1/(sweeps*lattice) - e0*e0/(sweeps*sweeps*lattice*lattice)) / (temperature * temperature))
60
61        #Printing the thermodynamic variables    
62
63        print(temperature,sum(Energy[RELAX_SWEEPS:]) / sweeps, sum(mag[RELAX_SWEEPS:]) / sweeps,sum(spec[RELAX_SWEEPS:]) / sweeps)
64
65
66
67main10()
68