%matplotlib inline
%config InlineBackend.figure_format = 'svg'
from triqs.plot.mpl_interface import plt, oplot

import numpy as np
from triqs.gf import *

The Two-Particle Self Consistency (TPSC) method

The Two-Particle Self Consistent (TPSC) method is an approach to solving the single-band Hubbard Model first introduced by Vilk and Tremblay (https://arxiv.org/abs/cond-mat/9702188).

Theory

Among the simplest approximations to self-energy and susceptibilities are the Hartree-Fock (HF) and Random-Phase-Approximation (RPA), which, in the single-band Hubbard model are

\[\Sigma^\text{HF}_\sigma(i\omega_n, \mathbf{k}) = U \braket{n_{-\sigma}}\]

\[\chi^\text{RPA}_\text{sp}(i\omega_n, \mathbf{k}) = \frac{\chi_0(i\omega_n, \mathbf{k})}{1 - \frac{U}{2}\chi_0(i\omega_n, \mathbf{k})}, \qquad \chi^\text{RPA}_\text{ch}(i\omega_n, \mathbf{k}) = \frac{\chi_0(i\omega_n, \mathbf{k})}{1 + \frac{U}{2}\chi_0(i\omega_n, \mathbf{k})}\]

The HF-self-energy has the immediate disadvantage of being constant and thus absorbed in the chemical potential.

For the susceptibilities, there exist two sum rules

\[\begin{split}\begin{align} \frac{T}{N}\sum_{i\omega_n, \mathbf{k}}{\chi_\text{sp}(i\omega_n, \mathbf{k})} &= n - 2\braket{n_\uparrow n_\downarrow} \\ \frac{T}{N}\sum_{i\omega_n, \mathbf{k}}{\chi_\text{ch}(i\omega_n, \mathbf{k})} &= n + 2\braket{n_\uparrow n_\downarrow} - n^2 \end{align}\end{split}\]

which are violated by the RPA-approximations. Also, since the denominator of \(\chi_\text{sp}\) can vanish, RPA predicts a magnetically ordered phase, which, in two dimensions is prohibited by the Mermin-Wagner theorem.

To fix this, TPSC introduces RPA-like susceptibilities with screened interaction vertices \(U_\text{ch}\) and \(U_\text{sp}\)

\[\chi^\text{TPSC}_\text{sp}(i\omega_n, \mathbf{k}) = \frac{\chi_0(i\omega_n, \mathbf{k})}{1 - \frac{U_\text{sp}}{2}\chi_0(i\omega_n, \mathbf{k})}, \qquad \chi^\text{TPSC}_\text{ch}(i\omega_n, \mathbf{k}) = \frac{\chi_0(i\omega_n, \mathbf{k})}{1 + \frac{U_\text{ch}}{2}\chi_0(i\omega_n, \mathbf{k})}\]

\(U_\text{ch}\) and \(U_\text{sp}\) are chosen such that the sum rules are fulfilled; for this to work, we need to now the double occupancy \(\braket{n_\uparrow n_\downarrow}\) for which TPSC makes the Ansatz

\[U_\text{sp} \braket{n_\uparrow}\braket{n_\downarrow} = U \braket{n_\uparrow n_\downarrow}\]

Having obtained the renormalized vertices and susceptibilities, TPSC then approximates the self-energy as

\[\Sigma^\text{TPSC}_\sigma(i\omega_n, \mathbf{k}) = U \braket{n_{-\sigma}} + \frac{U}{8}\frac{T}{N}\sum_{i\nu_m, \mathbf{q}}{\left[ 3 U_\text{sp}\chi_\text{sp}(i\nu_m, \mathbf{q}) + U_\text{ch}\chi_\text{ch}(i\nu_m, \mathbf{q}) \right] G_{0\sigma}(i\nu_m+i\omega_n, \mathbf{q}+\mathbf{k})}\]

Triangular lattice

In the following, we apply the TPSC method to a triangular two-dimensional Hubbard model. The triangular lattice in 2 dimensions is spanned by the primitive lattice vectors

\[\mathbf{a}_1 = a\hat{x}, \qquad \mathbf{a}_2 = \frac{a}{2} \hat{x} + \frac{\sqrt{3}a}{2}\hat{y}\]

For simplicity, we set \(a = 1\); in TRIQS we can create such a lattice as follows:

from triqs.lattice.tight_binding import TBLattice

t = -1.0
mu = -2.0

units = [(1, 0, 0), (1/2, np.sqrt(3)/2, 0)]

hoppings = {(+1,+0) : [[t]],
            (-1,+0) : [[t]],
            (+0,+1) : [[t]],
            (+0,-1) : [[t]],
            (+1,-1) : [[t]],
            (-1,+1) : [[t]],
            (+0,+0) : [[mu]]}

H = TBLattice(units=units, hoppings=hoppings)

TRIQS also allows us to calculate the dispersion relation and store it in an object e_k. We can plot this dispersion along high symmetry lines:

kmesh = H.get_kmesh(n_k=(64, 64, 1))

e_k = H.fourier(kmesh)

# define the high-symmetry points of the triangular lattice
G = [0.0, 0.0, 0.0]
M = [1/2, 0.0, 0.0]
K = [1/3, -1/3, 0.0]

paths = [(G, M), (M, K), (K, G)]

from triqs_tprf.lattice_utils import k_space_path
k_vecs, k_plot, k_ticks = k_space_path(paths, bz=H.bz, relative_coordinates=False)
k_vecs_rel, k_plot, k_ticks = k_space_path(paths, bz=H.bz, relative_coordinates=True)

e_k_interp_ref = [ H.tb.fourier(k)[0,0].real for k in k_vecs_rel]
e_k_interp = np.vectorize(lambda k : e_k(k)[0, 0].real, signature='(n)->()')

plt.plot(k_plot, e_k_interp(k_vecs))
plt.plot(k_plot, e_k_interp_ref, '--')
plt.xticks(k_ticks, [r'$\Gamma$',r'$M$',r'$K$',r'$\Gamma$'])
plt.ylabel(r'$\epsilon(\mathbf{k})$')
plt.grid(True)

Since it is a triangular lattice, the dispersion is not particle-hole symmetric. As we can see, this dispersion has a saddle point at the Fermi-level. In 2D, this will lead to a logarithmic divergence in the density of states, called a van Hove-singularity. The density of states can be calculated with the dos function:

from triqs.lattice.tight_binding import dos

DOS = dos(tight_binding=H.tb, n_kpts=2**10, n_eps=200, name='dos')[0]
plt.plot(DOS.eps, DOS.rho)
plt.xlabel('$\\varepsilon$')
plt.ylabel('$\\rho(\\varepsilon)$')

Text(0, 0.5, '$\\rho(\\varepsilon)$')

To get a better feeling for the dispersion, we can plot it as a function of \(\mathbf{k}\) within the first Brillouin-Zone. This also allows us to visualize the Fermi-Surface, which in this case is a hexagon (dashed lines):

e_k = H.fourier(H.get_kmesh(n_k=512))
k = np.linspace(-1.0, 1.0, num=200) * 2. * np.pi
Kx, Ky = np.meshgrid(k, k)

e_k_interp = np.vectorize(lambda kx, ky : e_k([kx, ky, 0])[0, 0].real)
e_k_interp = e_k_interp(Kx, Ky)

plt.figure()

g_abs = H.bz.units.T @ G
k_abs = H.bz.units.T @ K
m_abs = H.bz.units.T @ M

abs_path = [ g_abs, k_abs, m_abs, g_abs ]

K_X = [ p[0] for p in abs_path ]
K_Y = [ p[1] for p in abs_path ]

plt.plot(K_X, K_Y, '-', color='white')

plt.plot(g_abs[0], g_abs[1], 'co', label='G')
plt.plot(k_abs[0], k_abs[1], 'ro', label='K')
plt.plot(m_abs[0], m_abs[1], 'mo', label='M')

extent = (k.min(), k.max(), k.min(), k.max())
plt.imshow(e_k_interp, cmap=plt.get_cmap('RdBu'),
           extent=extent, origin='lower')
plt.colorbar()
plt.contour(Kx, Ky, e_k_interp, levels=[0], linestyles='dashed')
plt.xlabel(r'$k_x$'); plt.ylabel(r'$k_y$');
plt.legend(loc='upper left');

Non-interacting susceptibility \(\chi_0\)

TRIQS TPRF offers functions to calculate non-interacting Green’s functions and susceptibilities from a known dispersion e_k. Here, we calculate the static susceptibility \(\chi_0(i\omega_n=0, \mathbf{k})\):

from triqs_tprf.lattice import lattice_dyson_g0_wk
from triqs_tprf.lattice_utils import imtime_bubble_chi0_wk

beta = 100.

wmesh = MeshDLRImFreq(
    beta=beta, statistic='Fermion',
    eps=1e-6, # Tune to desired accuracy
    w_max=10.0, # Tune to energy range in problem
    )

g0_wk = lattice_dyson_g0_wk(mu=0., e_k=e_k, mesh=wmesh)
# get the density
g0_w = Gf(mesh=wmesh, target_shape=())
g0_w.data[:] = 1/len(e_k.mesh)*np.squeeze(np.sum(g0_wk.data, axis=1))
target_density = density(g0_w).real
chi0_wk = imtime_bubble_chi0_wk(g0_wk, nw=1, verbose=True)
iw0 = chi0_wk.mesh.components[0][0] # get zero Matsubara point to evaluate later

╔╦╗╦═╗╦╔═╗ ╔═╗  ┌┬┐┌─┐┬─┐┌─┐
 ║ ╠╦╝║║═╬╗╚═╗   │ ├─┘├┬┘├┤
 ╩ ╩╚═╩╚═╝╚╚═╝   ┴ ┴  ┴└─└
Two-Particle Response Function tool-box

beta  = 100.0
nk    = 262144
nw    = 32
norb  = 1

Approx. Memory Utilization: 0.38 GB

--> fourier_wk_to_wr
--> fourier_wr_to_tr
--> chi0_w0r_from_grt_PH (bubble in tau & r)
--> chi_wk_from_chi_wr (r->k)

We can now plot the static susceptibility:

chi0_k_interp = np.vectorize(lambda kx, ky : chi0_wk(Idx(0), [kx, ky, 0])[0, 0].real)
chi0_k_interp = chi0_k_interp(Kx, Ky)

plt.figure()

extent = (k.min(), k.max(), k.min(), k.max())
plt.imshow(chi0_k_interp, cmap=plt.get_cmap('RdBu'),
           extent=extent, origin='lower')
plt.colorbar()
plt.contour(Kx, Ky, e_k_interp, levels=[0], linestyles='dashed')
plt.xlabel(r'$k_x$'); plt.ylabel(r'$k_y$');

chi0_stat_k = np.array([chi0_wk(Idx(0), k)[0,0,0,0].real for k in k_vecs_rel])

fig, axes = plt.subplots(nrows=1, ncols=1)

axes.plot(k_plot, chi0_stat_k, label='$\\chi_0$')
axes.set_title('Absolute comparison')
axes.legend(bbox_to_anchor=(1,1))
axes.set_xticks(k_ticks, [r'$\Gamma$',r'$M$',r'$K$',r'$\Gamma$'])
axes.set_xlabel('$\\mathbf{k}$')
axes.grid(True)
axes.set_box_aspect(1)

The analytic form of this quantity is known:

\[\chi_0^\text{stat}(\mathbf{k}) = -\frac{2}{N} \sum_{\mathbf{q}}^{\text{1.BZ}} { \frac{f_\mu(\varepsilon(\mathbf{q})) - f_\mu(\varepsilon(\mathbf{q} + \mathbf{k}))} {\varepsilon(\mathbf{q}) - \varepsilon(\mathbf{q} + \mathbf{k})}},\]

where \(f_\mu\) is the Fermi-Dirac-statistic. From this expression, we can see that the peaks of \(\chi_0\) correspond to all possible particle-hole excitations. Whenever there exists a wave vector \(\mathbf{k}\) which maps a large portion of the Fermi-surface to another part of the Fermi-surface, then many terms in this sum will contribute to this peak; this phenomenon is known as Fermi-surface nesting and can be observed in the image above. In this case, it is visible at the M-point.

Interacting charge (\(\chi_{ch}\)) and spin response (\(\chi_{sp}\))

TRIQS TPRF has a tpsc_solver-class that can be used to do TPSC-calculations for one-band Hubbard models. The solver has to be supplied with the target density, the strength of the Hubbard interaction, an imaginary frequency mesh and a dispersion relation:

from triqs_tprf.tpsc_solver import tpsc_solver

U = 2.0
S = tpsc_solver(n=target_density, U=U, wmesh=wmesh, e_k=e_k)

S.solve(calc_sigma=False, calc_g=False, check_self_consistency=False)

╔╦╗╦═╗╦╔═╗ ╔═╗  ┌┬┐┌─┐┌─┐┌─┐
 ║ ╠╦╝║║═╬╗╚═╗   │ ├─┘└─┐│
 ╩ ╩╚═╩╚═╝╚╚═╝   ┴ ┴  └─┘└─┘
Two-Particle Self-Consistent method

  Using the TPSC-Ansatz for the double occupancy.

len(wmesh) = 21
nk = 262144
beta = 10.0
n = 0.7507420673321673
U = 2.0

Calculating non-interacting quantities...

Calculating first level of approximation...

   Calculating Usp...

   Calculating Uch...

   Calculating chisp_wk and chich_wk...

   Calculating double occupancy...

Summary first level approximation:
    Usp = 1.5896434157730261, Uch = 2.388267598779043
    Double Occupation <n_up*n_down> = 0.11199309130057132

 DONE!

FIXME need custom function to make chi DLR (target_shape!)

def make_chi_dlr(chi_wk):
    temp = Gf(mesh=chi_wk.mesh, target_shape=(1,1))
    temp.data[:,:,0,0] = chi_wk.data[:,:,0,0,0,0]
    return make_gf_dlr(temp)

We now plot the non-interacting, and the TPSC-charge- and spin-susceptibilities along high-symmetry paths in the 1. BZ. This allows us to easily compare them:

chi0_stat_k = make_chi_dlr(S.chi0_wk)(iw0, all)
chi0_stat_k = np.array([chi0_stat_k(k)[0,0].real for k in k_vecs_rel])

chich_stat_k = make_chi_dlr(S.chich_wk)(iw0, all)
chich_stat_k = np.array([chich_stat_k(k)[0,0].real for k in k_vecs_rel])

chisp_stat_k = make_chi_dlr(S.chisp_wk)(iw0, all)
chisp_stat_k = np.array([chisp_stat_k(k)[0,0].real for k in k_vecs_rel])

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 10))

axes[0].plot(k_plot, chi0_stat_k, label='$\\chi_0$')
axes[0].plot(k_plot, chich_stat_k, label='$\\chi_{ch}$')
axes[0].plot(k_plot, chisp_stat_k, label='$\\chi_{sp}$')
axes[0].set_title('Absolute comparison')
axes[0].legend(bbox_to_anchor=(1,1))
axes[0].set_xticks(k_ticks, [r'$\Gamma$',r'$M$',r'$K$',r'$\Gamma$'])
axes[0].set_xlabel('$\\mathbf{k}$')
axes[0].grid(True)
axes[0].set_box_aspect(1)

axes[1].plot(k_plot, chi0_stat_k - chi0_stat_k[0], label='$\\chi_0$')
axes[1].plot(k_plot, chich_stat_k - chich_stat_k[0], label='$\\chi_{ch}$')
axes[1].plot(k_plot, chisp_stat_k - chisp_stat_k[0], label='$\\chi_{sp}$')
axes[1].set_title('Relative comparison')
axes[1].legend(bbox_to_anchor=(1,1))
axes[1].set_xticks(k_ticks, [r'$\Gamma$',r'$M$',r'$K$',r'$\Gamma$'])
axes[1].set_xlabel('$\\mathbf{k}$')
axes[1].grid(True)
axes[1].set_box_aspect(1)

plt.tight_layout()

In the left image, we have plotted the absolute values of the susceptibilities; this shows that, when compared to the bare bubble, the spin fluctuations are increased in the interacting system, whereas charge fluctuations are reduced. At the same time, the pronounced peak at the M-point due to Fermi-surface nesting is still visible in all three susceptibilities.