Isospin modelinmg via metropolis algprithm

import numpy
import math
import random


def init_grid(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 step_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__":
    # initalize grid of spins in a nxm grid
    U = init_grid(10, 10)
    b = -1.0
    h = 0
    for _ in range(100*100):
        U = step_metropolis(U, b, h)
        print(U)