ISEG isospin model Nov 20, 2018 Snippet: import numpy import math import random def init_isig(Nx, Ny, random=True): if random: return 2*numpy.random.randint(2, size=(Nx, Ny)) -1 else: return numpy.ones(shape=(Nx,Ny)) def delta_energy(U, x, y, b, h): (Nx, Ny) = U.shape x_l = (x + Nx -1) % Nx x_h = (x + Nx +1) % Nx y_l = (y + Ny -1) % Ny y_h = (y + Ny +1) % Ny E = 2 * U[x,y] * (h + b * (U[x_l,y]+U[x_h,y]+U[x,y_l]+U[x,y_h])) return E def iseg_metropolis(U, b, h): (Nx, Ny) = U.shape x = random.randint(0,Nx-1) y = random.randint(0,Ny-1) dE = delta_energy(U, x, y, b, h) if random.random() < math.exp(-dE): U[x,y] = -U[x,y] return U if __name__ == "__main__": U = init_isig(10, 10) b = -1.0 h = 0 for _ in range(100*100): U = iseg_metropolis(U, b, h) print(U)