User Guide
This is the user guide for the omegamaxent_interface. The information presented here is sufficient for many cases, however, for a more advanced use, see the \(\Omega MaxEnt\) user guide.
In the following, we refer to the function to be analytically continued as the “Green’s function”, but it can be also a two-particle correlation function, a self-energy, or any function having a representation of the form
where \(\omega_n\) is a fermionic or bosonic Matsubara frequency, or
where the + and - signs apply to the fermionic and bosonic cases, respectively, and \(\mathbf{A}(\omega)\) is a positive semi-definite matrix in the fermionic case or \(\mathbf{A}(\omega)/\omega\) is positive semi-definite in the bosonic case.
Note
The code works internally with the imaginary frequency Green’s function. Therefore, if you provide \(\mathbf{G}(\tau)\), it will be Fourier-transformed by the code before the analytic continuation begins.
Obtain the retarded Green’s function
Whether you have a scalar-, a matrix-valued or a block Green’s function, the real frequency Green’s function is obtained with the function compute_GfReFreq(), which takes a TRIQS Matsubara Green’s function as input, and returns the corresponding retarded Green’s function. The function signature is:
GR=compute_GfReFreq(G, **kwa)
Input parameters
- G:
Only mandatory parameter of type Gf, GfImFreq, GfImTime, or BlockGf containing objects of one of those types.
Input Matsubara Green’s function.
The following parameters, to be passed as keyword arguments, are the most common ones. The more advanced parameters are defined in the dictionaries OmegaMaxEnt_input_params and OmegaMaxEnt_other_params (defined in OmegaMaxEnt_parameters.py) and are described in details in the \(\Omega MaxEnt\) user guide:
- ERR:
Optional real or complex numpy array.
Standard deviation if G is scalar or has a single element.
ERR must have the same shape as G.data.
See section Setting errors for non-diagonal covariance matrix.
- output_grid_params:
Optional list of the form \([\omega_{min}, \Delta\omega, \omega_{max}]\).
Defines the minimum frequency, the frequency step, and the maximum frequency, respectively, of the frequency grid on which the output Green’s function will be defined.
If not provided, the output grid is set by \(\Omega MaxEnt\). See section Frequency grids for more details.
Note: this is not the grid used in the calculation, which can be set instead with comp_grid_params.
- name:
Optional string.
Name parameter of the returned GfReFreq object.
- interactive_mode:
Optional boolean. Default: True.
See section Interactive mode for more details.
- save_figures_data:
Optional boolean. Default: True (scalar G only).
If G is scalar, \(\Omega MaxEnt\) saves files that allow to display the output figures after execution with the function display_figures(). See section Display figures for more details. Always False for non-scalar Green’s functions.
- save_G:
Optional boolean. Default: True.
By default, the result is saved in hdf5 format as ‘G’ in the file G_Re_Freq.h5.
- comp_grid_params:
Optional list of the form [\(\Delta\omega\)] or [\(\Delta\omega\), SW] or [\(\Delta\omega\), SW, SC].
Grid parameters used in the computation.
\(\Delta\omega\) is the frequency step used in the main spectral region, namely, the part of the grid where most of the spectral weight is located. It can also be set separately with parameter freq_step.
SW is the width of the main spectral region (should be between 2 and 4 standard deviations typically). It can also be set separately with parameter spectrum_width.
SC is the center of the main spectral region. It can also be set separately with parameter spectrum_center.
See section Frequency grids for more details.
- non_uniform_grid:
Optional boolean. Default: False.
Tells \(\Omega MaxEnt\) to use a non-uniform grid in the main spectral region for the computation.
This accelerates the calculation if the spectrum has a peak at zero frequency that is very narrow compared to the total width of the spectrum. See section Frequency grids for more details.
- inv_sym:
Optional boolean. Default: False.
If G is a matrix or a BlockGf, set inv_sym=True if \(G_{ji}=G_{ij}\). The calculation time for the off diagonal elements will then be halved.
- mu, nu:
Optional floats. Default: mu=1, nu=1.
Parameters involved in the calculation of off-diagonal elements of matrix-valued Green functions. See section Matrix-valued functions for more details.
- inv_sym_time:
Optional boolean. Default: False.
For bosonic Green’s function: if \(G(i\omega_n)=G(-i\omega_n)\), i.e. \(G(i\omega_n)\) is purely real, i.e. \(G(\tau)=G(-\tau)\), i.e. \(G(\beta-\tau)=G(\tau)\), set inv_sym_time=True. The MaxEnt calculation will then be performed over positive real frequencies only.
Return parameter
GR:
Interactive mode
If interactive_mode =True, \(\Omega MaxEnt\) displays figures during the execution. Also, if displ_preproc_figs=True, figures are displayed during the preprocessing stage. Otherwise, figures are displayed at the end of the calculation, showing the resulting Green’s function, along with different quantities used as diagnostic tools. Using those tools is very useful to assess, first, if the result is valid and, second, if it is the best result possible given the data. Therefore, when processing a set of data for the first time, it is strongly advised to use the interactive mode. Details about how to interpret the diagnostic quantities are given in the \(\Omega MaxEnt\) user guide.
You can also force the calculation to pause and display the results at different points by setting the minimal value of alpha with parameter alpha_min, or by setting the maximum number of values of alpha to be computed with parameter n_alpha_values and look at the results at different stages (i.e. different values of alpha). If you use the latter option, you can resume the calculation at the point of interuption after closing the figures.
Note that you can change parameters during a pause by modifying the file OmegaMaxEnt_input_params.dat and the changes will be applied when execution is resumed at the point of interruption. For example, if you have set alpha_min to a certain value (that parameter will appear as “minimum value of alpha:” in the file), you can modify that value (or remove the corresponding line completely) during the pause occuring when alpha_min is reached, and resume the computation at the next value of \(\alpha\) after having closed the figures. On the other hand, if you would want to add a new parameter during a pause, say new_param, in OmegaMaxEnt_input_params.dat you have to define it using the string OmegaMaxEnt_input_params[‘new_param’]. The parameter names are only understood by the python interface.
When the calculation is over and you are satisfied with the result displayed, you can exit the execution by closing all the figures and entering any character other than ‘y’ in the terminal. This will resume the execution of the python function compute_GfReFreq().
If interactive_mode=False, \(\Omega MaxEnt\) will not display any figure and compute_GfReFreq() will resume as soon as the calculation is over.
If you find that there are too many figures, instead of completely disable the interactive mode to eliminate all of them, you can reduce the number of figures displayed by setting to False one of the parameters: displ_alpha_opt_figs, displ_alpha_min_figs and displ_alpha_curves.
Note
For the continuation of matrix-valued Green’s functions, \(\Omega MaxEnt\) is called the same number of times as there are elements in the matrix (or in the upper part if inv_sym=True). If you are in interactive mode, figures showing the result will appear each time and, once you have closed them, you have to tell the program not to continue execution, to let the analytic continuation of the matrix or block function continue.
Setting errors
In the current version, you can provide errors only for a scalar-valued Green’s function. If the covariance matrix is diagonal, you can use parameter ERR to provide the standard deviation as a real or complex numpy array having the same shape as G.data. For a non-diagonal covariance, you can provide the name of the files containing the covariance matrix with parameter cov_tau for imaginary time data or cov_re_re, cov_im_im and cov_re_im for imaginary frequency data. The file type must be one of the valid armadillo types.
For matrix-valued Green’s function, the error is assumed to be constant. The value of that constant is not relevant since it has no effect on the results.
Imaginary time data
If your data is a scalar GfImTime and you do not have an estimate of the error, or the error is constant, do not provide errors. Otherwise, because \(\Omega MaxEnt\) works internally in Matsubara frequency, it will Fourier transform the covariance matrix, which takes time, but is not useful in that case because the result will also be a constant diagonal covariance in frequency, while the result does not depend on the absolute value of the error.
On the other hand, if the error depends on \(\tau\) and you do provide errors, either with parameter ERR or cov_tau (file name for a covariance matrix), note that the Fourier transform of the Green function is saved by default as a GfImFreq object called ‘G’ in file G_im_freq.h5 and the Fourier transform of the covariance matrix is saved in files covar_ReRe.dat, covar_ImIm.dat and covar_ReIm.dat in directory Fourier_transformed_data. This can be useful if you want to perform the continuation again on the same data, without having to wait during the Fourier transform of the covariance matrix, which takes some time if there are many \(\tau\) points. To do so, you pass to compute_GfReFreq() the saved GfImFreq object and the paths to the covariance files with parameters cov_re_re, cov_im_im and cov_re_im instead of the original GfImTime object and the error on \(G(\tau)\) with ERR.
Display figures
For a scalar Green’s function, if save_figures_data =True, regardless of the value of interactive_mode, you can display the same figures that are displayed in interactive mode by calling the function display_figures() after the execution of compute_GfReFreq(). For the matrix case, save_figures_data is always False. Details about the output figures are given in the \(\Omega MaxEnt\) user guide.
Frequency grids
There are two different real frequency grids: the output grid and the computational grid.
The output grid is the grid on which the output Green’s function is defined. You can control it with parameter output_grid_params. This frequency grid has a uniform density and is defined between \(\omega_{min}\) and \(\omega_{max}\) with a step \(\Delta\omega\). This is an optional parameter. If not provided, \(\Omega MaxEnt\) generates an output frequency grid that is usually well adapted to the spectrum.
For computational efficiency reasons, the real frequency grid used in the calculation is different from the output grid. In many cases the default computational grid generated by \(\Omega MaxEnt\) is well suited to the spectrum and there is no need to modify it. It can however happen that the calculation fails (no optimal value of alpha is found) because the default grid is not appropriate. Even when the calculation terminates successfully, the result might not always be as disired. For those cases you can use the input parameter comp_grid_params to control the computational grid in the region where most of the spectral weight is located. Outside that region, a particular non-uniform grid is always used by \(\Omega MaxEnt\). More advanced parameters are also available to control the computational grid in the dictionary OmegaMaxEnt_input_params (section FREQUENCY GRID PARAMETERS in OmegaMaxEnt_parameters.py). See the \(\Omega MaxEnt\) user guide for more details on those parameters.
For a spectrum having a peak centered at zero frequency that is very narrow compared to the total width of the spectrum, a simple way to optimize the computational grid is to set non_uniform_grid =True. \(\Omega MaxEnt\) will then use a grid with a density that is high in a narrow region around \(\omega=0\) and decreases as \(|\omega|\) increases. The detailed definition of this grid are given in the \(\Omega MaxEnt\) user guide.
Maximum entropy method and choice of entropy weight \(\alpha\)
The maximum entropy method consists in minimizing the functional
with
where \(G\) is the input data vector, \(A\) is a vector containing the spectrum values on a discretized frequency grid, \(\mathbf{K}\) is a matrix such that \(\mathbf{K}A\) performs the integral (1) (in \(\Omega MaxEnt\), but also (2) or another representation in general), \(C\) is the data’s covariance matrix, \(\alpha\) is an adjustable parameter, and \(S\) is a differential entropy, defined as
where \(D(\omega)\) is called the default model and is the solution that minimizes \(Q\) if \(\chi^2\) is negligible, i.e., at very large \(\alpha\). \(D(\omega)\) is defined in a way to include information known in advanced about the spectrum. By default, in \(\Omega MaxEnt\), \(D(\omega)\) is a gaussian with the same norm and first two moments associated with \(A(\omega)\), which are extracted from the input data.
The weight \(\alpha\) of the entropy term can be chosen in different ways. In \(\Omega MaxEnt\), the spectra are computed for a large range of \(\alpha\), starting at large \(\alpha\), and the optimal value is chosen where the curvature of \(log(\chi^2)\) as a function of \(\gamma log(\alpha)\) is maximal [1]. Here \(\gamma<1\) (parameter name: gamma) reduces the probability of a wrong value of \(\alpha\) to be chosen (default value: \(\gamma=0.2\)). Despite the use of \(\gamma\) and some smoothing of the curve \(log(\chi^2)\) vs \(\gamma log(\alpha)\) in the computation of the curvature, there is still a chance that a wrong value of \(\alpha\) will be selected because of some irregularities in \(log(\chi^2)\) vs \(\gamma log(\alpha)\) that produce parasitic peaks in the curvature. This is one of the reasons why the diagnostic tools are useful.
Matrix-valued functions
If the Green’s function is matrix-valued, the calculation is done using the auxiliary Green’s function approach described in [2] and also in appendix C of the \(\Omega MaxEnt\) user guide. In that calculation, the off-diagonal elements of the retarded function are obtained indirectly from the spectral functions of the diagonal elements and auxiliary functions that are linear combinations of diagonal and non-diagonal elements. Those auxiliary functions are constructed to have positive semi-definite spectral functions, so that they can be computed with the standard maximum entropy approach. As for the diagonal elements (of the form \(\langle T_{\tau} c_i(\tau) c_i^\dagger\rangle\)), they always have positive semi-definite spectral functions. Therefore, in that calculation, the retarded functions corresponding to the diagonal and the auxiliary Matsubara functions are first computed with \(\Omega MaxEnt\), and are then combined at the python level to obtain the retarded off-diagonal elements.
mu and nu
For Green or correlation functions of the form \(\langle T_{\tau} o_i(\tau) o_j^\dagger\rangle\), where \(o_i\) and \(o_j\) are operators corresponding to the same type of excitations, e.g. electronic excitations, the parameters mu and nu should be left equal to 1. On the other hand, for a correlation function of the form \(\langle T_{\tau} p(\tau) q^\dagger\rangle\), where \(p\) and \(q\) correspond to different types of excitations, different values of mu and nu should be tried to find a stable result. See Ref. [2] for more details.
If your matrix Green’s function has the symmetry \(G_{ji}=G_{ij}\), set inv_sym =True. Then, only the upper part of the matrix will actually be computed, which reduces the required computational demand by a factor of two.
Simple example of usage
Suppose you have saved a Matsubara Green’s function as a TRIQS object ‘G’ in a hdf5 file “G.h5”. The quickest way to obtain the corresponding real frequency Green’s function is:
from h5 import HDFArchive as HA
import OmegaMaxEnt_TRIQS as OT
#load the Matsubara Green's function
with HA("G.h5",'r') as A:
G=A['G']
#obtain the retarded Green's function
GR=OT.compute_GfReFreq(G)
Additionally, if you know that the spectrum has a sharp peak at \(\omega=0\) and that the spectrum is mostly located between \(\omega=-2\) and \(\omega=2\), you can use:
Dw=0.001
SW=4
GR=OT.compute_GfReFreq(G, comp_grid_params=[Dw,SW], non_uniform_grid=True)
For more advanced examples, see the tutorials.