# Wegner 2d Exact.py

```import math, cmath, random
random.seed('Werner')

L = 4
pp = 2.0 * math.pi / L

dim = 2
N = L ** dim # take even, dim = 2 is assumed
beta = math.sqrt(2.0)

knumbers = [(k, l) for l in range(L) for k in range(L)]
kminusnumbers = {}
for (k, l) in knumbers:
kminusnumbers[(k, l)] = ((L - k) % L, (L - l) % L)
rvalues = knumbers[:]
PhiR = {}
neighbours = {}
for rvec in rvalues:
PhiR[rvec] = (random.uniform(-1.0, 1.0)) # this is an example
neighbours[rvec] = [((rvec[0] + 1) % L, rvec[1]), (rvec[0], (rvec[1] + 1) % L),
((rvec[0] - 1) % L, rvec[1]), (rvec[0], (rvec[1] - 1) % L)]
epskvalues = {}
for (kx, ky) in knumbers:
epskvalues[(kx, ky)] = 4.0 * (math.sin(kx * pp / 2.0) ** 2 +
math.sin(ky * pp / 2.0) ** 2)
#
# Energy in PhiR
#
Energy_real = - 2.0 * dim * N
for rvec in rvalues:
for rvecp  in neighbours[rvec]:
Energy_real += (PhiR[rvec] - PhiR[rvecp]) ** 2 / 2.0
print Energy_real, 'Energy in PhiR'
#
# Fourier modes PhiK
#
PhiK = {}
for kvec in knumbers:
PhiK[kvec] = 0.0
for rvec in rvalues:
PhiK[kvec] += cmath.exp(- 1j * pp * (kvec[0] * rvec[0] + kvec[1] *
rvec[1])) * PhiR[rvec] / math.sqrt(N)
#
# Energy in PhiK
#
Energy_fourier = - 2.0 * dim * N
for kvec in knumbers:
kvecp = kminusnumbers[kvec]
Energy_fourier += epskvalues[kvec] * PhiK[kvec] * PhiK[kvecp]
print Energy_fourier, ' Energy in PhiK'
#
# Definition of positive Brillouin zone, and of symmetric points
#
kdummy = knumbers[:]
kdummy.remove((0, 0))
bplus = []
bsymm = []
while kdummy != []:
kvec = kdummy.pop()
if kvec ==  kminusnumbers[kvec]:
bsymm.append(kvec)
else:
kdummy.remove(kminusnumbers[kvec])
bplus.append(kvec)
print bplus, 'bplus'
print bsymm, 'bsymm'
PsiKpos = {}
PsiKneg = {}
#
# Fourier modes PsiK (real valued)
#
for kvec in bplus:
kvecp = kminusnumbers[kvec]
PsiKpos[kvec] = ((PhiK[kvec] + PhiK[kvecp]) / math.sqrt(2.0)).real
PsiKneg[kvec] = ((PhiK[kvec] - PhiK[kvecp]) / (math.sqrt(2.0) * 1j)).real
PsiKsymm = {}
for kvec in bsymm:
PsiKsymm[kvec] = PhiK[kvec].real
#
# Energy in PsiK
#
Energy_PsiK = - 2.0 * dim * N
for kvec in bplus:
Energy_PsiK += (PsiKpos[kvec] ** 2 + PsiKneg[kvec] ** 2) * epskvalues[kvec]
for kvec in bsymm:
Energy_PsiK += PsiKsymm[kvec] ** 2 * epskvalues[kvec]
print Energy_PsiK, ' Energy in PsiK'
#
# Correlation between site 0 and site k
#
print 'correlations using PhiR and PsiK'
for rvec in rvalues:
Corr_PsiK = 0.0
for kvec in bplus:
kr = (kvec[0] * rvec[0] + kvec[1] * rvec[1]) * pp
Corr_PsiK += \
(PsiKpos[kvec] * (1.0 - math.cos(kr)) + PsiKneg[kvec] * math.sin(kr)) \
* math.sqrt(2.0 / N)
for kvec in bsymm:
kr = (kvec[0] * rvec[0] + kvec[1] * rvec[1]) * pp
Corr_PsiK += \
PsiKsymm[kvec] * (1.0 - math.cos(kr)) * math.sqrt(1.0 / N)
```
```    print rvec, Corr_PsiK, PhiR[(0, 0)] - PhiR[rvec], ' r, Fourier term, direct term'
#
# Expectation value of correlation
#
xdata = []
ydata = []
for rvec  in rvalues:
summe = 0.0
for kvec in bplus:
kr = (kvec[0] * rvec[0] + kvec[1] * rvec[1]) * pp
summe += math.sin(kr / 2) ** 2 / epskvalues[kvec]
for kvec in bsymm:
kr = (kvec[0] * rvec[0] + kvec[1] * rvec[1]) * pp
summe += (1 - math.cos(kr)) ** 2 / (8.0 *epskvalues[kvec])
print rvec, math.exp(- 2.0 * summe / N / beta)
ydata.append( math.exp(- 2.0 * summe / N / beta))
xdata.append(rvec)
```