Documentation

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 the Legendre basis coefficients). Nonetheless, a typical script will the following basic parts.

  • Import TRIQS and SOM Python modules.

    # Green's function containers
    from pytriqs.gf.local import *
    # HDFArchive interface to .h5 files
    from pytriqs.archive import HDFArchive
    # HDF5 archives must be modified only by one MPI rank.
    # is_master_node() checks we are on rank 0.
    from pytriqs.utility.mpi import is_master_node
    
    # Main SOM class
    from pytriqs.applications.analytical_continuation.som import Som
    
  • Optional: Load input data from an HDF5 archive.

    # Open an HDF5 file in read-only mode
    arch = HDFArchive('input.h5', 'h')
    
    # Read input function
    inp = arch['input']
    # Read importance function S
    S = arch['S']
    

    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 Som object. Here you must specify the kind of the physical observable in question, or equivalently, the integral kernel in the analytical continuation equation.

    # Create Som object
    cont = Som(inp, S, kind = "FermionGf", norms = norms)
    

    inp and S must be Green’s function containers of the same type (one of GfImTime, GfImFreq or GfLegendre), defined on the same mesh and having the same target shape.

    Currently supported observable kinds are FermionGf (fermionic Green’s function or self-energy), BosonCorr (correlator of a boson-like operator with its Hermitian conjugate), BosonAutoCorr (correlator of a Hermitian operator with itself), and ZeroTemp (correlator computed in Matsubara formalism but at zero temperature). Refer to ‘Supported integral kernels’ for a detailed description of the implemented observable kinds.

    If target shape of inp is bigger than 1x1, SOM will only construct analytical continuation of the diagonal matrix elements.

    norms is an optional one-dimensional NumPy array containing 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, which is correct for the Green’s functions of the fermions. However, adjustments would normally be needed for self-energies and bosonic correlators/autocorelators.

  • Set simulation parameters.

    # Dictionary contaning all run() parameters
    run_params = {}
    
    # SOM will construct a spectral function within this energy window
    run_params['energy_window'] = (-5.0,5.0)
    
    # Number of particular solutions to calculate. It controls smoothness of
    # the final solution, and should in practice be chosen from 100-10000 range
    # (computation time scales linearly with it).
    run_params['l'] = 1000
    

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

    # Set verbosity level to the maximum on MPI rank 0, and silence all other ranks
    run_params['verbosity'] = 3 if is_master_node() else 0
    
    # Do not auto-adjust the number of particular solutions to accumulate
    # (default and recommended)
    run_params['adjust_l'] = False
    
    # Do not auto-adjust the number of global updates (default and recommended)
    run_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')
    run_params['f'] = 200
    
    # Number of local updates per global update
    # Has a similar effect on the algorithm as 'f', and scales linearly too.
    run_params['t'] = 500
    
    # Accumulate histogram of the objective function values,
    # useful to analyse quality of solution
    run_params['make_histograms'] = True
    

    See ‘run() parameters’ for a full list of accepted parameters.

  • Run simulation.

    cont.run(**run_params)
    
  • Extract solution and reconstructed input.

    SOM internally represents spectral functions as sums of rectangles. This representation is accessible as cont.solutions. In most cases one is interested in the resulting observable defined on a real frequency mesh.

    sol = cont.solutions
    # 'sol' is now a list of spectral functions, len(sol) == len(inp.indices)
    # Each element of 'sol' represents a sum of rectangles; its type is
    # list((rect_center,rect_width,rect_height))
    
    # Evaluate the solution on a frequency mesh.
    # The tail coefficients will be written into f_w too.
    #
    # NB: we can use *any* energy window at this point, not necessarily that from run_params
    f_w = GfReFreq(window = (-5.0, 5.0), n_points = 1000, indices = inp.indices)
    f_w << cont
    
    # Imaginary time/frequency/Legendre data reconstructed from the solution
    rec = inp.copy()
    rec << cont
    

    It is neccessary (but not enough) to have rec close to inp to ensure correctness of the 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['f_w'] = f_w
            ar['rec'] = rec
            # Save histograms for further analysis
            ar['histograms'] = cont.histograms
    

There are a few examples explaining how to run SOM for specific observable kinds.