Wegner 2d Exact.py
From Werner KRAUTH
(Difference between revisions)
Revision as of 23:20, 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 | ||
random.seed('Werner') | random.seed('Werner') | ||
- | + | ||
L = 4 | L = 4 | ||
pp = 2.0 * math.pi / L | pp = 2.0 * math.pi / L | ||
- | + | ||
dim = 2 | dim = 2 | ||
N = L ** dim # take even, dim = 2 is assumed | N = L ** dim # take even, dim = 2 is assumed | ||
beta = math.sqrt(2.0) | beta = math.sqrt(2.0) | ||
- | + | ||
knumbers = [(k, l) for l in range(L) for k in range(L)] | knumbers = [(k, l) for l in range(L) for k in range(L)] | ||
kminusnumbers = {} | kminusnumbers = {} | ||
Line 119: | Line 119: | ||
ydata.append( math.exp(- 2.0 * summe / N / beta)) | ydata.append( math.exp(- 2.0 * summe / N / beta)) | ||
xdata.append(rvec) | xdata.append(rvec) | ||
- | ~ |
Current revision
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)