Skip to content

Commit

Permalink
Fix test wick
Browse files Browse the repository at this point in the history
  • Loading branch information
@j-otsuki
j-otsuki committed Feb 27, 2018
1 parent 5bf763c commit deee9d6
Showing 1 changed file with 52 additions and 74 deletions.
126 changes: 52 additions & 74 deletions test/python/wick.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,21 +18,23 @@
h = 0.15 # Magnetic field

# Levels of the bath sites
epsilon = [-1.0, 1.0]
# epsilon = [-1.0, 1.0]
epsilon = [-1.0, 1.2]

# Hopping amplitudes
V = [0.5, 0.5]
# V = [0.5, 0.5]
V = [0.5, 0.6]

spin_names = ("up", "dn")

# Number of Matsubara frequencies for GF calculation
n_iw = 200

# Number of bosonic Matsubara frequencies for G^2 calculations

g2_n_iw = 5
g2_n_wb = 5

# Number of fermionic Matsubara frequencies for G^2 calculations
g2_n_inu = 5
g2_n_wf = 5

# GF structure
gf_struct = {'up' : [0], 'dn' : [0]}
Expand Down Expand Up @@ -73,76 +75,52 @@
# G^{(2)} #
###########

common_g2_params = {'gf_struct' : gf_struct,
common_g2_params = {'channel' : "PH",
'gf_struct' : gf_struct,
'beta' : beta,
'n_iw' : g2_n_iw}

# Compute G^{(2),ph}(i\omega;i\nu,i\nu'), AABB block order
G2_iw_inu_inup_ph_AABB = ed.G2_iw_inu_inup(channel = "PH",
block_order = "AABB",
n_inu = g2_n_inu,
**common_g2_params)

# Compute G^{(2),ph}(i\omega;i\nu,i\nu'), ABBA block order
G2_iw_inu_inup_ph_ABBA = ed.G2_iw_inu_inup(channel = "PH",
block_order = "ABBA",
n_inu = g2_n_inu,
**common_g2_params)

# Compute G^{(2),pp}(i\omega;i\nu,i\nu'), AABB block order
G2_iw_inu_inup_pp_AABB = ed.G2_iw_inu_inup(channel = "PP",
block_order = "AABB",
n_inu = g2_n_inu,
**common_g2_params)

# Compute G^{(2),pp}(i\omega;i\nu,i\nu'), ABBA block order
G2_iw_inu_inup_pp_ABBA = ed.G2_iw_inu_inup(channel = "PP",
block_order = "ABBA",
n_inu = g2_n_inu,
**common_g2_params)

b_mesh = MeshImFreq(beta = beta, S = 'Boson', n_max = g2_n_iw)
f_mesh = MeshImFreq(beta = beta, S = 'Fermion', n_max = g2_n_inu)
g2_mesh = MeshProduct(b_mesh, f_mesh, f_mesh)

# Check meshes of G2 containers
for bn1, bn2 in product(spin_names, spin_names):
assert G2_iw_inu_inup_ph_AABB[bn1, bn2].mesh == g2_mesh
assert G2_iw_inu_inup_ph_ABBA[bn1, bn2].mesh == g2_mesh
assert G2_iw_inu_inup_pp_AABB[bn1, bn2].mesh == g2_mesh
assert G2_iw_inu_inup_pp_ABBA[bn1, bn2].mesh == g2_mesh

G2_iw_inu_inup_ph_AABB_wick = G2_iw_inu_inup_ph_AABB.copy()
G2_iw_inu_inup_ph_ABBA_wick = G2_iw_inu_inup_ph_AABB.copy()
G2_iw_inu_inup_pp_AABB_wick = G2_iw_inu_inup_ph_AABB.copy()
G2_iw_inu_inup_pp_ABBA_wick = G2_iw_inu_inup_ph_AABB.copy()
'n_f' : g2_n_wf,
'n_b' : g2_n_wb, }

G2_ph_uuuu = ed.G2_iw( index1=('up',0), index2=('up',0), index3=('up',0), index4=('up',0), **common_g2_params )
G2_ph_dddd = ed.G2_iw( index1=('dn',0), index2=('dn',0), index3=('dn',0), index4=('dn',0), **common_g2_params )
G2_ph_uudd = ed.G2_iw( index1=('up',0), index2=('up',0), index3=('dn',0), index4=('dn',0), **common_g2_params )
G2_ph_dduu = ed.G2_iw( index1=('dn',0), index2=('dn',0), index3=('up',0), index4=('up',0), **common_g2_params )
G2_ph_uddu = ed.G2_iw( index1=('up',0), index2=('dn',0), index3=('dn',0), index4=('up',0), **common_g2_params )
G2_ph_duud = ed.G2_iw( index1=('dn',0), index2=('up',0), index3=('up',0), index4=('dn',0), **common_g2_params )

G2_ph_uuuu_wick = G2_ph_uuuu.copy()
G2_ph_dddd_wick = G2_ph_dddd.copy()
G2_ph_uudd_wick = G2_ph_uudd.copy()
G2_ph_dduu_wick = G2_ph_dduu.copy()
G2_ph_uddu_wick = G2_ph_uddu.copy()
G2_ph_duud_wick = G2_ph_duud.copy()

G = lambda s, i: G_iw[s].data[i + n_iw, 0, 0]

g2_mesh_enum = product(enumerate(b_mesh), enumerate(f_mesh), enumerate(f_mesh))
for (m, w), (i, nu), (ip, nup) in g2_mesh_enum:
m -= g2_n_iw - 1
i -= g2_n_inu
ip -= g2_n_inu
d_w_0 = int(m == 0)
d_w_nu_nup = int(m == i + ip + 1)
d_nu_nup = int(i == ip)


for s1, s2 in product(spin_names, spin_names):
d_s1_s2 = int(s1 == s2)

G2_iw_inu_inup_ph_AABB_wick[s1, s2].data[m + (g2_n_iw - 1), i + g2_n_inu, ip + g2_n_inu] = \
beta * d_w_0 * G(s1, i) * G(s2, ip) - beta * d_nu_nup * d_s1_s2 * G(s1, m + i) * G(s2, i)
G2_iw_inu_inup_ph_ABBA_wick[s1, s2].data[m + (g2_n_iw - 1), i + g2_n_inu, ip + g2_n_inu] = \
beta * d_w_0 * d_s1_s2 * G(s1, i) * G(s2, ip) - beta * d_nu_nup * G(s2, m + i) * G(s1, i)
G2_iw_inu_inup_pp_AABB_wick[s1, s2].data[m + (g2_n_iw - 1), i + g2_n_inu, ip + g2_n_inu] = \
beta * d_w_nu_nup * G(s1, i) * G(s2, ip) - beta * d_nu_nup * d_s1_s2 * G(s1, m - i - 1) * G(s2, i)
G2_iw_inu_inup_pp_ABBA_wick[s1, s2].data[m + (g2_n_iw - 1), i + g2_n_inu, ip + g2_n_inu] = \
beta * d_w_nu_nup * d_s1_s2 * G(s1, i) * G(s2, ip) - beta * d_nu_nup * G(s2, m - i - 1) * G(s1, i)

for s1, s2 in product(spin_names, spin_names):
assert_gfs_are_close(G2_iw_inu_inup_ph_AABB_wick[s1, s2], G2_iw_inu_inup_ph_AABB[s1, s2])
assert_gfs_are_close(G2_iw_inu_inup_ph_ABBA_wick[s1, s2], G2_iw_inu_inup_ph_ABBA[s1, s2])
assert_gfs_are_close(G2_iw_inu_inup_pp_AABB_wick[s1, s2], G2_iw_inu_inup_pp_AABB[s1, s2])
assert_gfs_are_close(G2_iw_inu_inup_pp_ABBA_wick[s1, s2], G2_iw_inu_inup_pp_ABBA[s1, s2])
for (i, j1, j2) in product(range(g2_n_wb),range(2*g2_n_wf),range(2*g2_n_wf)):
wb = i
wf1 = j1 - g2_n_wf
wf2 = j2 - g2_n_wf
d_w_0 = int(wb == 0)
# d_w_nu_nup = int(m == i + ip + 1)
d_nu_nup = int(j1 == j2)

G2_ph_uudd_wick[i,j1,j2] = beta * d_w_0 * G("up", wf1) * G("dn", wf2)
G2_ph_dduu_wick[i,j1,j2] = beta * d_w_0 * G("dn", wf1) * G("up", wf2)

G2_ph_uuuu_wick[i,j1,j2] = beta * d_w_0 * G("up", wf1) * G("up", wf2) - beta * d_nu_nup * G("up", wf1+wb) * G("up", wf1)
G2_ph_dddd_wick[i,j1,j2] = beta * d_w_0 * G("dn", wf1) * G("dn", wf2) - beta * d_nu_nup * G("dn", wf1+wb) * G("dn", wf1)

G2_ph_uddu_wick[i,j1,j2] = - beta * d_nu_nup * G("dn", wf1+wb) * G("up", wf1)
G2_ph_duud_wick[i,j1,j2] = - beta * d_nu_nup * G("up", wf1+wb) * G("dn", wf1)


# assert_gfs_are_close(G2_ph_uddu_wick, G2_ph_uddu)
gfs_are_close = lambda g1, g2: np.sum(g1-g2) < 1e-7

assert( gfs_are_close(G2_ph_uuuu_wick, G2_ph_uuuu) )
assert( gfs_are_close(G2_ph_dddd_wick, G2_ph_dddd) )
assert( gfs_are_close(G2_ph_uudd_wick, G2_ph_uudd) )
assert( gfs_are_close(G2_ph_dduu_wick, G2_ph_dduu) )
assert( gfs_are_close(G2_ph_uddu_wick, G2_ph_uddu) )
assert( gfs_are_close(G2_ph_duud_wick, G2_ph_duud) )

0 comments on commit deee9d6

Please sign in to comment.