Wegner 1d Direct.py

From Werner KRAUTH

Revision as of 23:16, 1 December 2019 Werner (Talk | contribs) ← Previous diff		Current revision Werner (Talk | contribs)
Line 1:		Line 1:
	import math, cmath, random		import math, cmath, random
-		+
	N = 8 # take even		N = 8 # take even
	beta = math.sqrt(2.0)		beta = math.sqrt(2.0)
	NIter = 500000000		NIter = 500000000
	print NIter, 'NIter'		print NIter, 'NIter'
-		+
	rvalues = range(N)		rvalues = range(N)
	kvalues = [2.0 * math.pi / N * float(l) for l in range(N)]		kvalues = [2.0 * math.pi / N * float(l) for l in range(N)]
Line 22:		Line 22:
	PhiK = [0] * N		PhiK = [0] * N
	PhiK[0] = 0.0		PhiK[0] = 0.0
-		+
	sigma = 1.0 / math.sqrt(8.0 * beta)		sigma = 1.0 / math.sqrt(8.0 * beta)
	PhiK[N / 2] = random.gauss(0.0, sigma)		PhiK[N / 2] = random.gauss(0.0, sigma)
Line 28:		Line 28:
	PhiK[k] = (PsiKpos[k] + 1j * PsiKneg[k]) / math.sqrt(2.0)		PhiK[k] = (PsiKpos[k] + 1j * PsiKneg[k]) / math.sqrt(2.0)
	PhiK[N - k] = (PsiKpos[k] - 1j * PsiKneg[k]) / math.sqrt(2.0)		PhiK[N - k] = (PsiKpos[k] - 1j * PsiKneg[k]) / math.sqrt(2.0)
-		+
	PhiR = []		PhiR = []
	for r in rvalues:		for r in rvalues:

---

Current revision

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'