Tutorial

Running SOM to analytically continue input data requires writing a simple Python script. Details of the script will vary depending on the physical quantity to be continued, and its representation (function of imaginary time, Matsubara frequencies, or a list of Legendre basis coefficients). Nonetheless, a typical script will have the following basic parts.

  • Import TRIQS and SOM Python modules.

    # Green's function containers
from triqs.gf import *
# HDFArchive interface to .h5 files
from h5 import HDFArchive
# HDF5 archives must be modified only by one MPI rank.
# is_master_node() checks whether we are on rank 0.
from triqs.utility.mpi import is_master_node

# Main SOM class
from som import Som

  • Optional: Load input data from an HDF5 archive .

    # Open an HDF5 file in read-only mode
arch = HDFArchive('input.h5', 'r')

# Read input function
inp = arch['input']

# Read estimated error bars or a full covariance matrix of
# the input data
error_bars = arch['error_bars']
cov_matrices = arch['cov_matrices']

    This step can be omitted if the input data comes from a different source. For instance, it could be loaded from text files or generated in the very same script by a quantum Monte-Carlo impurity solver.

  • Construct the Som object. Here, you must specify the kind of the physical observable in question.

    # Create Som object using the estimated error bars
cont = Som(inp, error_bars, kind = "FermionGf", norms = norms)

    # Create Som object using the covariance matrices
cont = Som(inp,
           cov_matrices,
           kind = "FermionGf",
           norms = norms,
           filtering_levels = 1e-5)

    inp and error_bars must be Green’s function containers of the same type triqs.gf.gf.Gf, defined on the same mesh (triqs.gf.meshes.MeshImTime, triqs.gf.meshes.MeshImFreq or triqs.gf.meshes.MeshLegendre) and having the same square target shape.

    Currently supported observable kinds are

    • FermionGf - fermionic Green’s function or self-energy;

    • FermionGfSymm - fermionic Green’s function or self-energy with enforced particle-hole symmetry.

    • BosonCorr - bosonic Green’s function, transverse magnetic susceptibility, general dynamical correlator of the form \(\chi_{OO^\dagger}(\tau) = \langle \mathbb{T}_\tau \hat O(\tau) \hat O^\dagger(0)\rangle\), where \(\hat O\) is a boson-like (fermion-number-conserving) operator.

    • BosonAutoCorr - charge susceptibility, longitudinal magnetic susceptibility, optical conductivity, general dynamical correlator of the form \(\chi_{OO}(\tau) = \langle \mathbb{T}_\tau \hat O(\tau) \hat O(0)\rangle\).

    • ZeroTemp - dynamical correlator computed in the Matsubara formalism but at zero temperature.

    Refer to ‘Supported observables’ for a detailed description of the implemented observable kinds.

    If target shape of inp is \(M{\times}M\) with \(M>1\), SOM will solve an independent analytic continuation problem for each of its diagonal matrix elements.

    norms is an optional one-dimensional NumPy array containing \(M\) expected norms of the spectral functions to be found, one positive real number per one diagonal element of inp. Instead of setting all elements of norms to the same constant x, one may simply pass norms = x. By default, all norms are set to 1.0 for the FermionGf, FermionGfSymm and ZeroTemp observables, which is the correct choice for the fermionic Green’s functions. However, adjustments would normally be needed for the fermionic self-energies. BosonCorr and BosonAutoCorr always require the norms to be explicitly specified. If those normalization constants are not known from some analytical considerations, they can be estimated from the input data itself using estimate_boson_corr_spectrum_norms().

    The covariance matrices \(\hat\Sigma\) must be provided as a container of type Gf(mesh=MeshProduct(inp.mesh, inp.mesh), target_shape=[M]). Each target element of this container holds one covariance matrix. It is also advised to specify positive filtering levels via the filtering_levels argument.

  • Set parameters for accumulation of particular solutions.

    # Dictionary containing all accumulate() parameters
acc_params = {}

# SOM will construct a spectral function within this energy window
acc_params['energy_window'] = (-5.0, 5.0)

# Number of particular solutions to accumulate. It controls smoothness
# of the final solution, and should in practice be chosen from the
# 100-10000 range (computation time scales linearly with it).
acc_params['l'] = 1000

    These two parameters are required for any simulation. Other useful parameters include verbosity and settings of the underlying Markov chain algorithm. You can find a detailed discussion of the latter in Sections 3 and 4 of A. S. Mishchenko’s lecture.

    # Set verbosity level to the maximum on MPI rank 0, and silence all
# other ranks.
acc_params['verbosity'] = 3 if is_master_node() else 0

# Do not auto-adjust the number of particular solutions to accumulate
# (default and recommended).
acc_params['adjust_l'] = False

# Do not auto-adjust the number of global updates (default and
# recommended).
acc_params['adjust_f'] = False

# Number of global updates per particular solution.
# Bigger values improve ergodicity of the Markov chain algorithm
# (computation time scales linearly with 'f').
acc_params['f'] = 200

# Number of local updates per global update.
# Has a similar effect on the algorithm as 'f', and scales linearly too.
acc_params['t'] = 500

# Enable the Consistent Constraints updates to speed-up accumulation.
acc_params['cc_update'] = True

# Accumulate histogram of the objective function values \chi,
# useful to analyse quality of solution.
acc_params['make_histograms'] = True

    See the documentation of the method Som.accumulate() for a full list of accepted parameters.

  • Accumulate particular solutions.

    cont.accumulate(**acc_params)

    Normally, this is the most time consuming step. Calling Som.accumulate() multiple times will incrementally extend the pool of the accumulated solutions.

  • Construct the final solution out of the accumulated particular solutions. One can choose to use either the standard SOM procedure and average selected particular solutions with equal weights,

    cont.compute_final_solution()

    or the more advanced SOCC procedure introduced in SOM 2.0.

    cont.compute_final_solution_cc()

    Behavior of Som.compute_final_solution_cc() can be fine-tuned via the many keyword arguments it accepts.

  • Extract the final solution, recover the real-frequency version of inp, its tail expansion coefficients and reconstruct the input.

    SOM internally represents spectral functions (configurations) as sums of rectangles. The final solutions are accessible as a list of Configuration objects via an attribute Som.solutions.

    sol = cont.solutions
# 'sol' is now a list of spectral functions,
# len(sol) == len(inp.indices)

# Recover the real-frequency version of 'inp'
#
# NB: we can use *any* energy window at this point, not necessarily that
# from 'acc_params'
w_mesh = MeshReFreq(window=(-5.0, 5.0), n_w=1000)
f_w = Gf(mesh=w_mesh, indices=inp.indices)
fill_refreq(f_w, cont)

# High frequency expansion (tail) coefficients up to order 'max_order'
# can be computed by function 'compute_tail()' and returned as a
# 3-dimensional complex NumPy array.
# The first index of the array is the zero-based expansion order.
tail = compute_tail(5, cont) # max_order = 5

# Imaginary time/frequency/Legendre data reconstructed from the solution
rec = inp.copy()
reconstruct(rec, cont)

    It is necessary (but not sufficient) to have rec close to inp to ensure good quality of the final solution.

  • Optional: Save results to an HDF5 archive.

    # On master node, save results to an archive
if mpi.is_master_node():
    with HDFArchive("results.h5", 'w') as ar:
        ar['inp'] = inp
        ar['sol'] = sol
        ar['f_w'] = f_w
        ar['tail'] = tail
        ar['rec'] = rec
        # Save histograms for further analysis
        ar['histograms'] = cont.histograms

  • Optional: Study statistical properties of the accumulated particular solutions.

    from som.spectral_stats import spectral_avg, spectral_disp, spectral_corr

# Compute statistical characteristics of the accumulated solution ensemble
# for the 0-th diagonal matrix element of 'inp' on the real-frequency
# mesh 'w_mesh' and using the rectangular resolution function.

# Averages over the ensemble
avg = spectral_avg(cont, 0, f_w.mesh, "rectangle")
# Dispersions of the ensemble
disp = spectral_disp(cont, 0, f_w.mesh, avg, "rectangle")
# Two-point correlations of the ensemble
corr = spectral_corr(cont, 0, f_w.mesh, avg, "rectangle")

There are a few examples explaining how to run SOM for various observables and how to use its more advanced features.