Example: Matsubara correlator at zero temperature¶

Observable kind: ZeroTemp.

Formally speaking, the imaginary time segment $$\tau\in[0;\beta)$$ turns into an infinite interval as $$\beta\to\infty$$. Similarly, spacing between Matsubara frequencies goes to 0 in that limit, and the difference between fermionic and bosonic Matsubaras disappears.

One can still define a correlation function on a finite time mesh $$\tau_i\in[0;\tau_{max}]$$, and assume the function is zero for $$\tau>\tau_{max}$$. In the frequency representation this corresponds to fictitious Matsubara spacing $$2\pi/\tau_{max}$$.

The spectral function is defined only on the positive half-axis of energy, since $$(1\pm e^{-\beta\epsilon})^{-1}$$ vanishes for negative $$\epsilon$$ in the zero temperature limit.

Run analytical continuation¶

# Import some TRIQS modules and NumPy
from pytriqs.gf.local import *
from pytriqs.archive import HDFArchive
import pytriqs.utility.mpi as mpi
import numpy

# Import main SOM class
from pytriqs.applications.analytical_continuation.som import Som

n_w = 1001               # Number of energy slices for the solution
energy_window = (0,10.0) # Energy window to search the solution in
# The window must entirely lie on the positive half-axis

# Parameters for Som.run()
run_params = {'energy_window' : energy_window}
# Verbosity level
run_params['verbosity'] = 3
# Number of particular solutions to accumulate
run_params['l'] = 5000
run_params['f'] = 100
# Number of local updates per global update
run_params['t'] = 50
# Accumulate histogram of the objective function values
run_params['make_histograms'] = True

# Could be g(i\omega_n) or g_l as well.
g_tau = HDFArchive('example.h5', 'r')['g_tau']

# Set the weight function S to a constant (all points of g_iw are equally important)
S = g_tau.copy()
S.data[:] = 1.0

# Construct a SOM object
# Norm of spectral function is known to be 0.5
cont = Som(g_tau, S, kind = "ZeroTemp", norms = numpy.array([0.5]))

# Run!
# Takes 1-2 minutes on 16 cores ...
cont.run(**run_params)

# Evaluate the solution on an energy mesh
# NB: we can use *any* energy window at this point, not necessarily that from run_params
g_w = GfReFreq(window = (0,10.0), n_points = n_w, indices = [0])
g_w << cont

# G(i\omega_n) reconstructed from the solution
g_rec_tau = g_tau.copy()
g_rec_tau << cont

# On master node, save results to an archive
if mpi.is_master_node():
with HDFArchive("results.h5",'w') as ar:
ar['g_tau'] = g_tau
ar['g_rec_tau'] = g_rec_tau
ar['g_w'] = g_w
ar['histograms'] = cont.histograms


Download input file example.h5.

Plot input and reconstructed imaginary-time correlators¶

from pytriqs.gf.local import *
from pytriqs.archive import HDFArchive
from matplotlib import pyplot as plt
from pytriqs.plot.mpl_interface import oplot

ar = HDFArchive('results.h5', 'r')

# Plot input and reconstructed correlator
oplot(ar['g_tau'],     mode='R', linewidth=0.8, label="$g(\\tau)$")
oplot(ar['g_rec_tau'], mode='R', linewidth=0.8, label="$g_\mathrm{rec}(\\tau)$")

plt.xlim((0, 10))
plt.ylabel("$g(\\tau)$")
plt.legend(loc="upper right")


Plot the spectral function¶

from pytriqs.gf.local import *
from pytriqs.archive import HDFArchive
from matplotlib import pyplot as plt
from pytriqs.plot.mpl_interface import oplot

oplot(ar['g_w'], mode='S', linewidth=0.8, label="$A(\\omega)$")
plt.ylabel("$A(\\omega)$")