From Werner KRAUTH
import math, cmath, random
N = 8 # take even
beta = math.sqrt(2.0)
NIter = 500000000
print NIter, 'NIter'
rvalues = range(N)
kvalues = [2.0 * math.pi / N * float(l) for l in range(N)]
epskvalues = [4.0 * math.sin( k / 2.0) **2 for k in kvalues]
correlation = [0.0] * N
for iter in range(NIter):
#
# real-valued Fourier coefficients, don't forget to count PsiK[N/2] = PhiK[N/2]
#
PsiKpos = [0.0]
PsiKneg = [0.0]
for k in range(1, N / 2):
sigma = 1.0 / math.sqrt(2.0 * beta * epskvalues[k])
PsiKpos.append(random.gauss(0.0, sigma))
PsiKneg.append(random.gauss(0.0, sigma))
PhiK = [0] * N
PhiK[0] = 0.0
sigma = 1.0 / math.sqrt(8.0 * beta)
PhiK[N / 2] = random.gauss(0.0, sigma)
for k in range(1, N / 2):
PhiK[k] = (PsiKpos[k] + 1j * PsiKneg[k]) / math.sqrt(2.0)
PhiK[N - k] = (PsiKpos[k] - 1j * PsiKneg[k]) / math.sqrt(2.0)
PhiR = []
for r in rvalues:
dummy = 0.0
for k in range(N):
term = cmath.exp(+ 1j * r * kvalues[k]) * PhiK[k] / math.sqrt(N)
dummy += term
PhiR.append(dummy.real)
#
# Compute correlation between site 0 and site k
#
for i in range(N):
correlation[i] += math.cos(PhiR[0] - PhiR[i])
for i in range(N):
print i, correlation[i] / NIter, 'correlation with i'