GW self-energy of a 2D square lattice Hubbard model

[1]:
%matplotlib inline
import sys, os
import numpy as np
from triqs.plot.mpl_interface import plt,oplot

from h5 import HDFArchive
from triqs.atom_diag import *

from triqs.gf import *
from triqs.operators import c, c_dag, n, dagger
from itertools import product


from triqs.lattice.tight_binding import TBLattice
from triqs.lattice.utils import k_space_path

from triqs_tprf.lattice import lattice_dyson_g0_wk, lattice_dyson_g_wk, lattice_dyson_g0_fk, dynamical_screened_interaction_W, lattice_dyson_g_fk
from triqs_tprf.gw import bubble_PI_wk, gw_sigma, lindhard_chi00, g0w_sigma

from triqs_Nevanlinna import Solver

from triqs_maxent import *

import seaborn as sns
import scienceplots
plt.style.use(['science','notebook'])
sns.set_palette('muted')
Warning: could not identify MPI environment!
Starting serial run at: 2024-08-08 11:54:11.515443
/mnt/sw/nix/store/29h1dijh98y9ar6n8hxv78v8zz2pqfzf-python-3.11.7-view/lib/python3.11/site-packages/numpy/core/getlimits.py:549: UserWarning: The value of the smallest subnormal for <class 'numpy.float64'> type is zero.
  setattr(self, word, getattr(machar, word).flat[0])
/mnt/sw/nix/store/29h1dijh98y9ar6n8hxv78v8zz2pqfzf-python-3.11.7-view/lib/python3.11/site-packages/numpy/core/getlimits.py:89: UserWarning: The value of the smallest subnormal for <class 'numpy.float64'> type is zero.
  return self._float_to_str(self.smallest_subnormal)
/mnt/sw/nix/store/29h1dijh98y9ar6n8hxv78v8zz2pqfzf-python-3.11.7-view/lib/python3.11/site-packages/numpy/core/getlimits.py:549: UserWarning: The value of the smallest subnormal for <class 'numpy.float32'> type is zero.
  setattr(self, word, getattr(machar, word).flat[0])
/mnt/sw/nix/store/29h1dijh98y9ar6n8hxv78v8zz2pqfzf-python-3.11.7-view/lib/python3.11/site-packages/numpy/core/getlimits.py:89: UserWarning: The value of the smallest subnormal for <class 'numpy.float32'> type is zero.
  return self._float_to_str(self.smallest_subnormal)

Setup simple two orbital 2D Hubbard model on square lattice

[2]:
n_orb = 1

loc = np.zeros((n_orb,n_orb))
for i in range(n_orb):
    for j in range(n_orb):
        if i != 0 and i==j:
            loc[i,j] = -0.3
        if j > i or j < i:
            loc[i,j] = -0.5

t = -1.0 * np.eye(n_orb)   #nearest neighbor hopping
tp = 0.1 * np.eye(n_orb)    #next nearest neighbor hopping

hop= {  (0,0)  :  loc,
        (1,0)  :  t,
        (-1,0) :  t,
        (0,1)  :  t,
        (0,-1) :  t,
        (1,1)  :  tp,
        (-1,-1):  tp,
        (1,-1) :  tp,
        (-1,1) :  tp}

TB = TBLattice(units = [(1, 0, 0) , (0, 1, 0)], hoppings = hop, orbital_positions= [(0., 0., 0.)]*n_orb)

plot dispersion along high-symmetry lines

[3]:
n_pts = 101
G = np.array([ 0.00,  0.00,  0.00])
M = np.array([0.5, 0.5, 0.0])
X = np.array([0.5, 0.0, 0.0])
R = np.array([0.5, 0.5, 0.5])
paths = [(R, G), (G, X), (X, M), (M, G)]

kvecs, k, ticks  = k_space_path(paths, num=n_pts, bz=TB.bz)

e_mat = TB.fourier(kvecs)
e_val = TB.dispersion(kvecs)
eps_k = {'k': k, 'K': ticks, 'eval': e_val, 'emat' : e_mat}
[4]:
fig, ax = plt.subplots(1,1, figsize=(8,4), dpi=150, squeeze=False)
ax = ax.reshape(-1)

for band in range(eps_k['eval'].shape[1]):
    ax[0].plot(eps_k['k'], eps_k['eval'][:,band].real, color='C0', zorder=2.5)

ax[0].axhline(y=0,zorder=2,color='gray',alpha=0.5,ls='--')
ax[0].set_xticks(eps_k['K'])
ax[0].set_xticklabels([r'R', '$\Gamma$', 'X', 'M', r'$\Gamma$'])
ax[0].set_xlim([eps_k['K'].min(), eps_k['K'].max()])
ax[0].set_ylabel(r'$\omega$ (eV)')

plt.show()
../_images/tutorials_GW_tutorial_6_0.png

GW in imaginary time

[5]:
k_grid = (30,30,1)
k_mesh = TB.get_kmesh(n_k = k_grid)
e_k = TB.fourier(k_mesh)
eps_k = TB.dispersion(k_mesh)
mu = 0.

beta = 10
n_iw = 400
iw_mesh = MeshImFreq(beta=beta, S='Fermion', n_max=n_iw)
G0_iwk = lattice_dyson_g0_wk(mu=mu, e_k=e_k, mesh=iw_mesh)

setup bare interaction

[6]:
def construct_U_kan(n_orb, U, J, Up=None, Jc=None):

    orb_range = range(0, n_orb)
    U_kan = np.zeros((n_orb, n_orb, n_orb, n_orb))

    if not Up:
        Up = U-2*J
    if not Jc:
        Jc = J

    for i,j,k,l in product(orb_range, orb_range, orb_range, orb_range):
        if i == j == k == l: # Uiiii
            U_kan[i, j, k, l] = U
        elif i == k and j == l: # Uijij
            U_kan[i, j, k, l] = Up
        elif i == l and j == k: # Uijji
            U_kan[i, j, k, l] = J
        elif i == j and k ==l: # Uiijj
            U_kan[i, j, k, l] = Jc
    return U_kan
[7]:
U=1
Up=0.2
J=0.4
V_k = Gf(mesh=k_mesh, target_shape=[n_orb*1]*4)

V_bare = np.zeros((n_orb,n_orb,n_orb,n_orb))
# simple onsite intra orbital Coulomb repulsion
for i in range(n_orb):
    for j in range(n_orb):
        if i == j:
            V_bare[i,i,j,j] = U
        else:
            V_bare[i,i,j,j] = Up

V_bare = construct_U_kan(n_orb=n_orb,U=U,J=J,Up=Up)


V_k.data[:] = V_bare

Run one GW loop

[8]:
print('--> pi_bubble')
PI_iwk = bubble_PI_wk(G0_iwk)

print('--> screened_interaction_W')
Wr_iwk = dynamical_screened_interaction_W(PI_iwk, V_k)

print('--> gw_self_energy')
Σ_iwk = gw_sigma(Wr_iwk, G0_iwk)

print('--> lattice_dyson_g_wk')
G_wk = lattice_dyson_g_wk(mu, e_k, Σ_iwk)

Σ_Γ_iw = Σ_iwk[:, Idx(0,0,0)]
Σ_X_iw = Σ_iwk[:, Idx(k_grid[1]-1,0,0)]
--> pi_bubble
--> screened_interaction_W
--> gw_self_energy
--> lattice_dyson_g_wk

Plot results

[9]:
fig, ax = plt.subplots(n_orb, n_orb, figsize=(8*n_orb,4*n_orb), dpi=150, squeeze=False,sharex=True)
ax = ax.reshape(-1)

ax[0].oplot(Σ_Γ_iw[i,j].imag, label=f'[{i},{j}]@$\Gamma$')
ax[0].oplot(Σ_X_iw[i,j].imag, label=f'[{i},{j}]@X')

ax[0].set_xlim(0,50)

ax[0].set_xlabel(r'$i\omega_n$')
ax[0].set_ylabel(r'Im $\Sigma (i\omega_n)$ (eV)')

plt.tight_layout(pad=0.4, w_pad=0.1, h_pad=0.4)
plt.show()
../_images/tutorials_GW_tutorial_15_0.png

GW on real frequency axis

[10]:
# make sure 0 is not in the mesh! Divergence for q=[0,0,0]
# no large freq mesh is needed. kmesh critical for convergence here
n_w = 100
delta = 0.1
GW_window = (-15, 15)
w_mesh = MeshReFreq(window=GW_window, n_w=n_w)


G0_wk = lattice_dyson_g0_fk(mu=mu, e_k=e_k, mesh=w_mesh, delta=delta)
[11]:
print('--> pi_bubble')
PI_wk = lindhard_chi00(e_k=e_k, mesh=w_mesh, beta=beta, mu=mu, delta=delta)

print('--> screened_interaction_W')
Wr_wk = dynamical_screened_interaction_W(PI_wk, V_k)

print('--> gw_self_energy')
Σ_wk = g0w_sigma(mu=mu, beta=beta, e_k=e_k, W_fk=Wr_wk, v_k=V_k, delta=delta)

print('--> lattice_dyson_g_wk')
g_fk = lattice_dyson_g_fk(mu, e_k, Σ_wk, delta)

Σ_Γ_w = Σ_wk[:, Idx(0,0,0)]
Σ_X_w = Σ_wk[:, Idx(k_grid[1]-1,0,0)]
--> pi_bubble
--> screened_interaction_W
--> gw_self_energy
--> lattice_dyson_g_wk

Analytic continuation

Nevanlinna

[12]:
# setup Nevanlinna kernel solver
solver = Solver(kernel='CARATHEODORY')

# For the caratheodory formalism we have to subtract the Hartree shift
solver.solve(Σ_Γ_iw-Σ_Γ_iw(Σ_Γ_iw.mesh.last_index()))

Σ_w_mesh = MeshReFreq(window=GW_window, n_w=n_w*2)
Σ_Γ_w_CN = solver.evaluate(Σ_w_mesh, delta)
[13]:
# setup Nevanlinna kernel solver
solver = Solver(kernel='NEVANLINNA')

solver.solve(Σ_Γ_iw)

Σ_w_mesh = MeshReFreq(window=GW_window, n_w=n_w*2)
Σ_Γ_w_N = solver.evaluate(Σ_w_mesh, delta)
This is Nevanlinna analytical continuation. All off-diagonal elements will be ignored.

Pade

[14]:
Σ_Γ_w_P = Σ_Γ_w_N.copy()

Σ_Γ_w_P.set_from_pade(Σ_Γ_iw, n_points=n_iw, freq_offset=delta)

MaxEnt

[15]:
# Initialize SigmaContinuator
isc = DirectSigmaContinuator(Σ_Γ_iw)
[16]:
tm = TauMaxEnt()
tm.set_G_iw(isc.Gaux_iw)
tm.set_error(1.e-4)
tm.omega = HyperbolicOmegaMesh(omega_min=GW_window[0], omega_max=GW_window[1], n_points=300)
tm.alpha_mesh = LogAlphaMesh(alpha_min=1e-2, alpha_max=5e2, n_points=60)
result = tm.run()
print(('LineFit: ', result.analyzer_results['LineFitAnalyzer']['alpha_index']))
2024-08-08 11:56:43.802678
MaxEnt run
TRIQS application maxent
Copyright(C) 2018 Gernot J. Kraberger
Copyright(C) 2018 Simons Foundation
Authors: Gernot J. Kraberger and Manuel Zingl
This program comes with ABSOLUTELY NO WARRANTY.
This is free software, and you are welcome to redistributeit under certain conditions; see file LICENSE.
Please cite this code and the appropriate original papers (see documentation).

Minimal chi2: 3.1012476779564424
scaling alpha by a factor 2401 (number of data points)
alpha[ 0] =   1.20050000e+06, chi2 =   1.86098288e+05, n_iter=      22
alpha[ 1] =   9.99352274e+05, chi2 =   1.47692182e+05, n_iter=       4
alpha[ 2] =   8.31907511e+05, chi2 =   1.17033657e+05, n_iter=       4
alpha[ 3] =   6.92518670e+05, chi2 =   9.25372229e+04, n_iter=       4
alpha[ 4] =   5.76484888e+05, chi2 =   7.29614635e+04, n_iter=       4
alpha[ 5] =   4.79892947e+05, chi2 =   5.73366260e+04, n_iter=       4
alpha[ 6] =   3.99485305e+05, chi2 =   4.48992103e+04, n_iter=       4
alpha[ 7] =   3.32550227e+05, chi2 =   3.50392785e+04, n_iter=       4
alpha[ 8] =   2.76830342e+05, chi2 =   2.72618500e+04, n_iter=       4
alpha[ 9] =   2.30446507e+05, chi2 =   2.11604692e+04, n_iter=       4
alpha[10] =   1.91834436e+05, chi2 =   1.63995950e+04, n_iter=       4
alpha[11] =   1.59691945e+05, chi2 =   1.27025480e+04, n_iter=       4
alpha[12] =   1.32935034e+05, chi2 =   9.84259726e+03, n_iter=       4
alpha[13] =   1.10661332e+05, chi2 =   7.63572388e+03, n_iter=       4
alpha[14] =   9.21196614e+04, chi2 =   5.93432809e+03, n_iter=       4
alpha[15] =   7.66847089e+04, chi2 =   4.62157115e+03, n_iter=       4
alpha[16] =   6.38359336e+04, chi2 =   3.60624445e+03, n_iter=       4
alpha[17] =   5.31400128e+04, chi2 =   2.81813047e+03, n_iter=       4
alpha[18] =   4.42362288e+04, chi2 =   2.20384924e+03, n_iter=       4
alpha[19] =   3.68243030e+04, chi2 =   1.72320363e+03, n_iter=       4
alpha[20] =   3.06542699e+04, chi2 =   1.34605177e+03, n_iter=       4
alpha[21] =   2.55180461e+04, chi2 =   1.04973178e+03, n_iter=       4
alpha[22] =   2.12424135e+04, chi2 =   8.17037936e+02, n_iter=       4
alpha[23] =   1.76831772e+04, chi2 =   6.34707633e+02, n_iter=       8
alpha[24] =   1.47203026e+04, chi2 =   4.92342230e+02, n_iter=       8
alpha[25] =   1.22538675e+04, chi2 =   3.81665499e+02, n_iter=       8
alpha[26] =   1.02006917e+04, chi2 =   2.96023747e+02, n_iter=       9
alpha[27] =   8.49153220e+03, chi2 =   2.30046206e+02, n_iter=       9
alpha[28] =   7.06874803e+03, chi2 =   1.79404990e+02, n_iter=       9
alpha[29] =   5.88435603e+03, chi2 =   1.40634231e+02, n_iter=       9
alpha[30] =   4.89841281e+03, chi2 =   1.10984279e+02, n_iter=       9
alpha[31] =   4.07766763e+03, chi2 =   8.82981058e+01, n_iter=      10
alpha[32] =   3.39444099e+03, chi2 =   7.09038474e+01, n_iter=      10
alpha[33] =   2.82569123e+03, chi2 =   5.75209639e+01, n_iter=      10
alpha[34] =   2.35223737e+03, chi2 =   4.71790629e+01, n_iter=      11
alpha[35] =   1.95811226e+03, chi2 =   3.91489155e+01, n_iter=      11
alpha[36] =   1.63002410e+03, chi2 =   3.28851983e+01, n_iter=      12
alpha[37] =   1.35690820e+03, chi2 =   2.79802848e+01, n_iter=      12
alpha[38] =   1.12955376e+03, chi2 =   2.41281078e+01, n_iter=      13
alpha[39] =   9.40293314e+02, chi2 =   2.10968267e+01, n_iter=      13
alpha[40] =   7.82744074e+02, chi2 =   1.87088091e+01, n_iter=      14
alpha[41] =   6.51592729e+02, chi2 =   1.68263575e+01, n_iter=      15
alpha[42] =   5.42416222e+02, chi2 =   1.53416908e+01, n_iter=      15
alpha[43] =   4.51532599e+02, chi2 =   1.41699005e+01, n_iter=      16
alpha[44] =   3.75876826e+02, chi2 =   1.32438974e+01, n_iter=      18
alpha[45] =   3.12897427e+02, chi2 =   1.25106612e+01, n_iter=      18
alpha[46] =   2.60470433e+02, chi2 =   1.19283579e+01, n_iter=      18
alpha[47] =   2.16827755e+02, chi2 =   1.14640690e+01, n_iter=      19
alpha[48] =   1.80497551e+02, chi2 =   1.10919807e+01, n_iter=      21
alpha[49] =   1.50254592e+02, chi2 =   1.07919330e+01, n_iter=      22
alpha[50] =   1.25078941e+02, chi2 =   1.05482528e+01, n_iter=      23
alpha[51] =   1.04121553e+02, chi2 =   1.03488096e+01, n_iter=      25
alpha[52] =   8.66756437e+01, chi2 =   1.01842461e+01, n_iter=      27
alpha[53] =   7.21528544e+01, chi2 =   1.00473518e+01, n_iter=      28
alpha[54] =   6.00634061e+01, chi2 =   9.93255836e+00, n_iter=      30
alpha[55] =   4.99995848e+01, chi2 =   9.83554439e+00, n_iter=      32
alpha[56] =   4.16219898e+01, chi2 =   9.75293415e+00, n_iter=      35
alpha[57] =   3.46480884e+01, chi2 =   9.68207603e+00, n_iter=      38
alpha[58] =   2.88426872e+01, chi2 =   9.62088069e+00, n_iter=      40
alpha[59] =   2.40100000e+01, chi2 =   9.56769988e+00, n_iter=      42
MaxEnt loop finished in 0:00:41.606576
('LineFit: ', 40)
[17]:
isc.set_Gaux_w_from_Aaux_w(result.analyzer_results['LineFitAnalyzer']['A_out'], result.omega, np_interp_A=n_w*2,
                           np_omega=n_w, w_min=GW_window[0], w_max=GW_window[1])

Σ_Γ_w_ME = isc.S_w

Comparison: Pade, Nevanlinna, Pade, and Maxent

[18]:
fig, ax = plt.subplots(n_orb, n_orb, figsize=(8*n_orb,4*n_orb), dpi=150, squeeze=False,sharex=True)
ax = ax.reshape(-1)

ax[0].oplot(Σ_Γ_w[i,j].imag, '-o', c='gray', label=f'ref')
ax[0].oplot(Σ_Γ_w_N[i,j].imag, label=f'Nev', zorder=10)
ax[0].oplot(Σ_Γ_w_CN[i,j].imag, label=f'Car', zorder=6)
ax[0].oplot(Σ_Γ_w_P[i,j].imag, label=f'Pade')
ax[0].oplot(Σ_Γ_w_ME[i,j].imag, label=f'MaxEnt')


ax[0].set_xlim(GW_window)
ax[0].set_ylim(0,0.15)

ax[0].set_xlabel(r'$\omega$')
ax[0].set_ylabel(r'Im $\Sigma (\omega)$ (eV)')


plt.tight_layout(pad=0.4, w_pad=0.1, h_pad=0.4)
plt.show()

fig, ax = plt.subplots(n_orb, n_orb, figsize=(8*n_orb,4*n_orb), dpi=150, squeeze=False,sharex=True)
ax = ax.reshape(-1)

# plotting results
ax[0].oplot(Σ_Γ_w[i,j].real, '-o', c='gray', label=f'ref')
ax[0].oplot(Σ_Γ_w_N[i,j].real, label=f'Nev', zorder=10)
ax[0].oplot((Σ_Γ_w_CN[i,j]+Σ_Γ_iw(Σ_Γ_iw.mesh.last_index())[i,j]).real, label=f'Car', zorder=6)
ax[0].oplot(Σ_Γ_w_P[i,j].real, label=f'Pade')
ax[0].oplot(Σ_Γ_w_ME[i,j].real, label=f'MaxEnt')

ax[0].set_xlim(GW_window)
ax[0].set_xlabel(r'$\omega$')
ax[0].set_ylabel(r'Re $\Sigma (\omega)$ (eV)')

plt.tight_layout(pad=0.4, w_pad=0.1, h_pad=0.4)
plt.show()
../_images/tutorials_GW_tutorial_30_0.png
../_images/tutorials_GW_tutorial_30_1.png