Ising coupling.py

From Werner KRAUTH

import math
import random
 
L = 64; N = L * L
nbr = {i : ((i // L) * L + (i + 1) % L, (i + L) % N,
             (i // L) * L + (i - 1) % L, (i - L) % N) \
                                    for i in range(N)}
beta = 0.4  #  beta = 0.4407 is critical temperature
SLow = [-1 for site in range(N)]
SHigh = [1 for site in range(N)]
iter = 0
while True:
    iter += 1
    k = random.randint(0, N - 1)
    Upsilon = random.uniform(0.0, 1.0)
    hLow = sum(SLow[nn] for nn in nbr[k])
    hHigh = sum(SHigh[nn] for nn in nbr[k])
    SLow[k] = -1
    if Upsilon < 1.0 / (1.0 + math.exp(-2.0 * beta * hLow)):
        SLow[k] = 1
    SHigh[k] = -1
    if Upsilon < 1.0 / (1.0 + math.exp(-2.0 * beta * hHigh)):
        SHigh[k] = 1
    if SHigh == SLow:
        print(iter / N)
        print(SLow, SHigh)
        break