-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathIsingModel.py
146 lines (115 loc) · 4.44 KB
/
IsingModel.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
import random
import matplotlib.pyplot as plt
import math
# ---Some constants---
N = 25 # Size of grid
J = 1 # Coupling strength
kBT = 2.5 * J # Critical temperature occurs at approximately 2.25~2.27 * J without diagonal contributions and at approximately 5.5 * J with diagonal contributions
# ---2D case---
# ---Initiate a grid of zeros---
grid = []
for i in range(N):
row = []
for j in range(N):
row.append(0)
grid.append(row)
# ---Fill the grid with spins of -1 and +1 chosen randomly---
for i in range(N):
for j in range(N):
spin = random.choice([-1, 1])
grid[i][j] = spin
print(grid)
# ---Actual Metropolis algorithm implementation---
for i in range (N * 10000): # Iteration count
# Pick a random spin
x = random.randrange(0, N)
y = random.randrange(0, N)
# Neighbouring spins
x_left = x - 1
x_right = x + 1
y_down = y - 1
y_up = y + 1
# Periodic boundary conditions
if x_left < 0:
x_left = x_left + N
if x_right >= N:
x_right = x_right - N
if y_down < 0:
y_down = y_down + N
if y_up >= N:
y_up = y_up - N
'''
If the spin at (x, y) flips, the only affected terms in the sum of the total energy
will be the terms that contain the flipped spin.
Instead of calculating the total initial and final energies of the system
if the flip occurs, we only need to calculate the initial and final energies
of the neighbouring pairs to the flipped spin.
'''
# No diagonal contributions
Ei = -J * ((grid[x][y]) * (grid[x_left][y] + grid[x_right][y] + grid[x][y_down] + grid[x][y_up]))
Ef = -J * ((-grid[x][y]) * (grid[x_left][y] + grid[x_right][y] + grid[x][y_down] + grid[x][y_up]))
# With diagonal contributions (uncomment to include diagonal neighbours)
# Ei = -J * ((grid[x][y]) * (grid[x_left][y] + grid[x_right][y] + grid[x][y_down] + grid[x][y_up] + grid[x_left][y_up] + grid[x_right][y_up] + grid[x_left][y_down] + grid[x_right][y_down]))
# Ef = -J * ((-grid[x][y]) * (grid[x_left][y] + grid[x_right][y] + grid[x][y_down] + grid[x][y_up] + grid[x_left][y_up] + grid[x_right][y_up] + grid[x_left][y_down] + grid[x_right][y_down]))
# Calculate the change in energy of the system if the flip occurs
DeltaE = Ef - Ei
# Flip condition
if DeltaE <= 0:
grid[x][y] = -grid[x][y]
else:
P = math.exp(-DeltaE / kBT) # Probability
p = random.uniform(0, 1) <= P
if p == True: # Flip
grid[x][y] = -grid[x][y]
else: # No flip
grid[x][y] = grid[x][y]
# ---Plotting (I like blue)---
fig, ax = plt.subplots()
ax.matshow(grid, cmap='Blues_r') # _r reverses the colormap spectrum
plt.show()
'''
# ---1D case---
# ---Initiate an array of zeros---
array = []
for i in range(N):
array.append(0)
# ---Fill the array with spins of -1 and +1 chosen randomly---
for i in range(N):
spin = random.choice([-1, 1])
array[i] = spin
# ---Actual Metropolis algorithm implementation---
for i in range (N * 10000): # Iteration count
# Pick a random spin
x = random.randrange(0, N)
# Neighbouring spins
x_left = x - 1
x_right = x + 1
# Periodic boundary conditions
if x_left < 0:
x_left = x_left + N
if x_right >= N:
x_right = x_right - N
# If the spin at (x, y) flips, the only affected terms in the sum of the total energy will be the terms that contain the flipped spin.
# Instead of calculating the total initial and final energies of the system if the flip occurs, we only need to calculate the initial and final energies of the neighbouring pairs to the flipped spin.
Ei = -J * ((array[x]) * (array[x_left] + array[x_right]))
Ef = -J * ((-array[x]) * (array[x_left] + array[x_right]))
# Calculate the change in energy of the system if the flip occurs
DeltaE = Ef - Ei
# Flip condition
if DeltaE <= 0:
array[x] = -array[x]
else:
P = math.exp(-DeltaE / kBT) # Probability
p = random.uniform(0, 1) <= P
if p == True: # Flip
array[x]= -array[x]
else: # No flip
array[x] = array[x]
# ---Plotting (I like blue)---
grid = []
grid.append(array)
fig, ax = plt.subplots()
ax.get_yaxis().set_visible(False)
ax.matshow(grid, cmap='Blues_r') # _r reverses the colormap spectrum
plt.show()
'''