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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
import random
from numpy import around, amax, argmax, array, zeros, linspace, arange, where, abs, zeros, ones, full
from scipy.integrate import simps
from scipy.optimize import curve_fit
import warnings
import os
import time
import numpy as np
from numpy import array, around, linspace, exp, sqrt, log, cosh
import matplotlib.pyplot as plt
from scipy.ndimage.interpolation import shift
from scipy.optimize import curve_fit
from math import *
import numpy as np
import matplotlib.pyplot as plt
'''
series of functions which will extract all major order parameters from the
Ising simulation
'''
def magnetization(lattice):
N_s = np.prod(lattice.shape)
return abs((1/N_s) * np.sum(lattice))
def magnetization_2(grid):
''' calculate <m^2>'''
N_s = np.prod(grid.shape)
# if we do grid*grid, this is some constant number...
return ((1/N_s))*np.sum(grid*grid)
def energy(grid, J=1):
'''
this function must be able to generallize to arbitrary dimensions
:param grid:
:param J:
:return:
'''
N = grid.shape
dimension = len(N) # dimension of system being simulated
NN = np.copy(grid)
E = 0
neighbors = 0
for i in range(dimension):
for j in [-1, 1]:
neighbors += np.roll(NN, shift=j, axis=i)
E += J*np.sum(grid*np.roll(NN, shift=j, axis=i))
DeltaE = J * (grid * neighbors)/(np.prod(N))
# -(E/np.prod(N))/2; #return is avg energy per site
return -np.sum(DeltaE)/2
def energy_2(grid, J=1):
''' measurement of energy^2 per site, which should be constant'''
N_s = np.prod(grid.shape)
N = grid.shape
dimension = len(N) # dimension of system being simulated
NN = np.copy(grid)
E = 0
for i in range(dimension):
for j in [-1, 1]:
E += J*(grid*np.roll(NN, shift=j, axis=i))
return np.sum(E*E)/N_s/2
# compute the two-point spin-spin correlator of the 2D Ising model
def s0sr(lattices, moment=None):
N = lattices[0].shape
points = []
for lattice in lattices:
s0 = lattice[0, 0] # this could be any point on the lattice
for i in range(N[0]):
x = min([i, N[0]-i]) # because of the periodic boundary conditions
for j in range(N[1]):
y = min([j, N[1]-j]) # because of the periodic boundary conditions
r = np.around(sqrt(x**2 + y**2), 3)
if moment is None: # the correlator <s(0)s(r)>
points.append([r, s0*lattice[i, j]])
else: # some moment like <r^2s(0)s(r)>
points.append([r, r**moment*s0*lattice[i, j]])
# sort the list by increasing values of the distance r
points = sorted(points, key=lambda x: x[0])
new_points = []
k = 0
while k < len(points):
counter = 1
new_point = points[k]
if k < len(points)-1:
next_point = points[k+1]
# they have the same value of r, so average them together
while new_point[0] == next_point[0]:
new_point[1] += next_point[1]
counter += 1
if k+counter < len(points):
next_point = points[k+counter]
else:
next_point = [100000, 100000] # arbitrary point that breaks loop
new_points.append([new_point[0], new_point[1]/counter])
k += counter
r_vals = []
s0sr_vals = []
for point in new_points:
r_vals.append(point[0])
s0sr_vals.append(point[1])
return r_vals, s0sr_vals
def susceptibility(grid, beta):
'''
formula chi = 1/T(<m^2> - <m>^2)
:param grid:
:param beta:
:return:
'''
m = magnetization(grid)
m_2 = magnetization_2(grid)
chi = beta*(m_2 - m**2)
return chi
def heat_capacity(grid, K):
E = energy(grid, K)
E2 = energy_2(grid, K)
cv = K**2*(E2 - E**2)
return cv
def spin_spin_correlator(x, cl, a, eta):
return exp(-x/cl)*(a/(a+x))**(eta)
def correlation_length(x, c, nu, Tc):
return c*(x-Tc)**(-nu)
def fit_spin_correlator_data(beta_vals, beta_r_vals, beta_s0sr_vals, rmin, rmax):
cl_vals = []
cl_std_vals = []
beta_s0sr_fit_vals = []
truncated_beta_r_vals = []
truncated_beta_s0sr_vals = []
# truncate r values over the range rmin to rmax
for k in range(len(beta_vals)):
r_vals = []
s0sr_vals = []
for j in range(len(beta_r_vals[k])):
r = beta_r_vals[k][j]
if r >= rmin and r <= rmax:
r_vals.append(beta_r_vals[k][j])
s0sr_vals.append(beta_s0sr_vals[k][j])
truncated_beta_r_vals.append(r_vals)
truncated_beta_s0sr_vals.append(s0sr_vals)
# Now do the fit:
for k in range(len(beta_vals)):
cl, cl_cov = curve_fit(spin_spin_correlator, truncated_beta_r_vals[k], truncated_beta_s0sr_vals[k], bounds=([1, 0.2, 0.24999999], [64, 0.20000001, 0.25]))
print("beta = ",beta_vals[k])
print("from fit cl = ",around(cl[0],3))
print("from fit a = ",around(cl[1],3))
print("from fit eta = ",around(cl[2],3))
cl_vals.append(around(cl[0],3))
cl_std_vals.append(around(sqrt(np.diag(cl_cov)),3))
s0sr_fit_vals = []
for r in beta_r_vals[k]:
s0sr_fit_vals.append(spin_spin_correlator(r,cl[0], cl[1], cl[2]))
beta_s0sr_fit_vals.append(s0sr_fit_vals)
return beta_r_vals, beta_s0sr_vals, cl_vals, cl_std_vals, beta_s0sr_fit_vals
def fit_temp_cl_data(T_vals, cl_vals):
N64_T_vals = T_vals
N64_cl_vals = cl_vals
print("N64_T_vals =", np.around(N64_T_vals, 3))
print("N64_cl_vals =", N64_cl_vals)
# Do the fit for the temperature dependence of the correlation length:
popt, pcov = curve_fit(correlation_length, T_vals,
cl_vals, bounds=([0, 0, 1.5], [10, 3, 3]))
print("\n Best fit to cl = c*(T-Tc)^(-nu) is c= ", np.around(popt[0], 3), " nu=", np.around(
popt[1], 3), " Tc=", np.around(popt[2], 4), " beta_c=", np.around(1/popt[2], 4), "\n")
cl_fit_vals = []
cl_exact_vals = []
cl_finite_size_vals = []
for T in T_vals:
cl_fit_vals.append(correlation_length(T, popt[0], popt[1], popt[2]))
# Here we use the large volume limit value of Tc
cl_exact_vals.append(correlation_length(T, 1.0, 1, 2.269))
# Here was use the extracted value of Tc
cl_finite_size_vals.append(correlation_length(T, 1.0, 1, np.around(popt[2], 4)))
return cl_fit_vals, cl_exact_vals, cl_finite_size_vals