Wegner 2d Exact.py

From Werner KRAUTH

(Difference between revisions)
Jump to: navigation, search
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)
Personal tools