DMFT lattice susceptibility
In this guide we will compute the uniform magnetic susceptibility \(\chi = \chi(\mathbf{Q} = \mathbf{0})\) of the single band Hubbard model on the square lattice with nearest neighbour hopping using dynamical mean-field theory (DMFT).
We will do this in two very different ways from
self consistent DMFT calculations in applied fields, and from
the lattice Bethe-Salpeter Equation using the DMFT local vertex.
Since DMFT is thermodynamically consistent [Hafermann et al., PRB 90, 235105 (2014)] these two approaches gives the same susceptibility within error-bars.
Self consistent calculations in applied field
First we compute the uniform magnetic susceptibility \(\chi\) by appyling an external magnetic field \(B\) and measuring the resulting magnetization \(M\). The susceptibility \(\chi\) is then given as the derivative at zero field, i.e.
In practice this is done by performing a series of self-consistent DMFT calculations at finite fields \(B \ne 0\). To do this we will use a lightweight framework for doing self consistent DMFT calculations taylored to the Hubbard model on the square lattice, for details see DMFT self consistent framework.
The framework python module common.py is available here: common.py, and is based on the TPRF helper class ParameterCollection that makes it possible to store the results of a calculation in a single object that, in turn, can be stored and passed as argument.
To determine \(\chi\) we perform a series of self consistent calculations for several values of the magnetic field \(B\). Using the DMFT framework this looks like
from common import *
p = ParameterCollection(
t=1., B=0., U=10., mu=0., n_k=16, n_iter=10, G_l_tol=2e-5,
solve = ParameterCollection(
length_cycle = 10, n_warmup_cycles = 10000, n_cycles = int(2.5e6),
move_double = False, measure_G_l = True
),
init = ParameterCollection(
beta = 1., n_l = 10, n_iw = 400, n_tau = 4000,
gf_struct = [('up',[0]), ('do',[0])]
),
)
p = setup_dmft_calculation(p)
for B in [0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.4, 0.6, 0.8, 1.0]:
p.B = B
ps = solve_self_consistent_dmft(p)
p = ps[-1]
if mpi.is_master_node():
with HDFArchive('data_B_{:f}.h5'.format(p.B), 'w') as a:
a['ps'] = ParameterCollections(ps)
The result is a magnetization \(M\) vs. applied magnetic field \(B\) curve, to which we fit a low order polynomial, and compute its zero field derivative.
The visualization script is available here: plot_field.py.
The resulting homogenous static spin susceptbility has the value
This concludes the first method for computing static susceptibilities within DMFT.
Lattice susceptibility from the Bethe-Salpeter Equation
Instead of multiple calculations in applied field the susceptibility \(\chi\) can be obtained from a direct calculation of the two-particle susceptbility using the DMFT local vertex and the lattice Bethe-Salpeter equation.
To this end one has to perform the steps
compute the DMFT impurity two-particle Green’s function \(G^{(2)}\),
compute the DMFT impurity two-particle vertex \(\Gamma\), and
solve the lattice Beht-Salpeter Equation for the lattice susceptibility \(\chi(\mathbf{Q})\).
DMFT local vertex
To obtain the DMFT impurity magnetic vertex \(\Gamma_m\) one first computes the DMFT impurity two-particle Green’s function in the particle-hole channel \(G^{(2,PH)}_{abcd}(i\omega, i\nu, i\nu')\), where \(a,b,c,d\) are spin-orbital indices. The two-particle Green’s function can be sampled using triqs_cthyb and to compute the static susceptibility it is sufficient to only keep the zero Bosonic frequency \(\omega = 0\).
For the single band model the magnetic susceptbility \(\chi_m\) is directly related to \(G^{(2,PH)}\) as
and the magnetic bubble susceptibility \(\chi^{(0)}_m\) is given by \(\chi^{(0)}_m(i\omega. i\nu, i\nu') = - \beta \delta_{\nu, \nu'} G(i\nu) G(i\omega + i\nu)\). The two susceptibilities are related by the impurity Bethe-Salpeter Equation giving the corresponding magetic vertex function \(\Gamma_m\) as
Starting from the self consistent DMFT solution at zero-field we run triqs_cthyb to sample the two-particle Green’s function and then compute \(\chi_m\), \(\chi_m^{(0)}\), and \(\Gamma_m\) see below
from common import *
if mpi.is_master_node():
with HDFArchive('data_B_0.000000.h5', 'r') as a:
p = a['ps'].objects[-1]
else: p = None
p = mpi.bcast(p)
# -- Sample G2
p.solve.n_cycles = int(1e9 / 40.)
p.solve.measure_G_l = False
p.solve.measure_G_tau = False
p.solve.measure_G2_iw_ph = True
p.solve.measure_G2_blocks = set([('up','up'), ('up','do')])
p.solve.measure_G2_n_bosonic = 1
p.solve.measure_G2_n_fermionic = 20
cthyb = triqs_cthyb.Solver(**p.init.dict())
cthyb.G0_iw << p.G0_w
cthyb.solve(**p.solve.dict())
p.G2_iw_ph = cthyb.G2_iw_ph.copy()
# -- Compute DMFT impurity vertex
from triqs_tprf.linalg import inverse_PH
from triqs_tprf.chi_from_gg2 import chi0_from_gg2_PH
p.chi_m = p.G2_iw_ph[('up','up')] - p.G2_iw_ph[('up','do')]
p.chi0_m = chi0_from_gg2_PH(p.G_w['up'], p.chi_m)
p.gamma_m = inverse_PH(p.chi0_m) - inverse_PH(p.chi_m)
del p.solve.measure_G2_blocks
if mpi.is_master_node():
with HDFArchive('data_g2.h5', 'w') as a: a['p'] = p
The resulting response functions are plotted below
The visualization script is available here: plot_g2.py.
Lattice Bethe-Salpeter Equation
Equipped with the DMFT local vertex \(\Gamma_m\) it is possible to compute the DMFT lattice susceptibility \(\chi(\mathbf{Q})\) from the lattice Bethe-Salpeter Equation (BSE)
TPRF comes with an OpenMP amd MPI parallelized BSE solver triqs_tprf.bse.solve_lattice_bse. However, the calculation is done with a fixed number of frequencies \(n_\nu\) in the fermionic frequencies \(\nu\) and \(\nu'\), and the solution converges only linearly with the size of the frequency window. Therefore we solve the BSE for a range of window sizes to enable extrapolation \(N_\nu \rightarrow \infty\).
from common import *
from triqs_tprf.linalg import inverse_PH
from triqs_tprf.chi_from_gg2 import chi0_from_gg2_PH
from triqs_tprf.utilities import G2_loc_fixed_fermionic_window_python
from triqs_tprf.lattice import lattice_dyson_g_wk
from triqs_tprf.bse import solve_lattice_bse
# -- Solve the lattice BSE for several fermionic window sizes
for nwf in [8, 10, 12, 20]:
with HDFArchive('data_g2.h5', 'r') as a: p = a['p']
p.nwf = nwf
# -- DMFT impurity vertex
p.chi_m = p.G2_iw_ph[('up','up')] - p.G2_iw_ph[('up','do')]
p.chi_m = G2_loc_fixed_fermionic_window_python(p.chi_m, nwf=p.nwf)
p.chi0_m = chi0_from_gg2_PH(p.G_w['up'], p.chi_m)
p.gamma_m = inverse_PH(p.chi0_m) - inverse_PH(p.chi_m)
# -- Lattice BSE
g_wk = lattice_dyson_g_wk(mu=p.mu, e_k=p.e_k, sigma_w=p.sigma_w)[0:1, 0:1]
p.chi_kw, p.chi0_kw = solve_lattice_bse(g_wk, p.gamma_m)
with HDFArchive('data_bse_nwf{:03d}.h5'.format(nwf), 'w') as a: a['p'] = p
The resuls along the high symmetry path of the Brillouin zone is shown below for fixed \(N_\nu\) (left panel) and the extrapolation for the \(\Gamma\)-point is also shown (right panel).
The visulaization script is available here: plot_bse.py.
The result for the homogeneous magnetic susceptbilitiy \(\chi(\mathbf{0})\) from the BSE is
in quantitative agreement with the applied field value.
Summary
Now we can compare the two results for the homogeneous static magnetic susceptibility, from
the applied field calculation, and
the BSE calculation.
The results are in quantitative agreement
and the accuracy is limited by the stochastic Monte Carlo noise in the vertex. This can be improved further by increasing the number of samples in the triqs_cthyb two-particle Green’s function calculation.
Note, that the BSE approach is much more general than the applied field approach. The BSE calculation gives the whole momentum dependent lattice susceptibility \(\chi(\mathbf{Q})\) and provides dynamical information when using finite Bosonic freqiencies \(|\omega| > 0\).