ising2d.py

#!/usr/bin/python3
import numpy as np
from numpy import random

N = 10
size = [N, N]
theta = 0.42
n_steps = 1000000

x = 2 * random.random_integers(0, 1, size=size) - 1

print('initial x =\n', x)

flip_count = 0
f1_sum = 0.0
f2_sum = 0.0
for step in range(n_steps):
    i, j = random.random_integers(0, N-1), random.random_integers(0, N-1)
    i1 = i - 1 if i > 0 else N - 1
    i2 = i + 1 if i < N - 1 else 0
    j1 = j - 1 if j > 0 else N - 1
    j2 = j + 1 if j < N - 1 else 0
    sum_neighbors = x[i1, j] + x[i2, j] + x[i, j1] + x[i, j2]
    r = np.exp(-2.0*theta*x[i, j] * sum_neighbors)
    R = random.uniform()
    # print('r = %f, R = %f' % (r, R))
    if R < r:
        x[i, j] = -x[i, j]
        flip_count += 1
    f1_sum += x[0, 0]
    f2_sum += x[0, 0]*x[N/2, N/2]

print('flips = %d, flip ratio =%f' % (flip_count, flip_count/n_steps))
print('average f1 = %f' % (f1_sum / n_steps))
print('average f2 = %f' % (f2_sum / n_steps))