Analytically continue using SOM

[abdelghany0]

I tried to follow this tutorial to perform analytically continue using SOM. I used this Python script "somAC.py":

from triqs.gf import *
from triqs.operators import *
from h5 import HDFArchive
import triqs.utility.mpi as mpi
from triqs_cthyb import Solver
from math import pi
from som import Som

def cut_coefficients(g_l, n_remaining_coeffs):
    g_cut = GfLegendre(indices = [0], beta = g_l.beta, n_points = n_remaining_coeffs)
    g_cut.data[:,:,:] = g_l.data[:n_remaining_coeffs, :, :]
    g_cut.enforce_discontinuity(np.identity(g_cut.N1))
    return g_cut

n_tau = 500                 # Number of tau-slices for the input GF
n_l = 30                    # Number of Legendre polynomials for the input GF
n_w = 801                   # Number of energy slices for the solution
energy_window = (-6.0, 6.0)  # Energy window to search the solution in

# Parameters for Som.run()
run_params = {'energy_window' : energy_window}
# Verbosity level
run_params['verbosity'] = 2
# Number of particular solutions to accumulate
run_params['l'] = 2000
# Number of global updates
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'] = False

def SOM(input_filename="dmft_results.h5", output_filename="som_results.h5"):

# Read G(\tau) from archive
# Could be G(i\omega_n) or G_l as well.
    with HDFArchive(input_filename,'r') as QMC:
        G_tau = QMC["G_tau"]
        G_iw = QMC["G_iw"]
        G_l = QMC["G_l"]

    # Paramagnetic case: average spin up and spin down GF
    g_tau = (G_tau['up'] + G_tau['down']) / 2.
    g_iw = (G_iw['up'] + G_iw['down']) / 2.
    g_l = (G_l['up'] + G_l['down']) / 2.

    # Prepare input data: reduce the number of \tau-slices from 10001 to n_tau
    # reduce the number of Legendre coefficients to n_l
    g_tau_rebinned = rebinning_tau(g_tau, n_tau)
    g_l_cut = cut_coefficients(g_l, n_l)

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

    S_l = g_l_cut.copy()
    S_l.data[:] = 1.0

    # Construct a SOM object
    #cont = Som(g_tau_rebinned, S_tau, kind = "FermionGf")
    cont = Som(g_l_cut, S_l, kind = "FermionGf")

    # Run!
    cont.run(**run_params)

    # Create a real frequency GF obtained with SOM
    g_w = GfReFreq(window = energy_window, n_points = n_w, indices = [0])
    g_w << cont

    # G(\tau) reconstructed from the SOM solution
    g_rec_tau = g_tau_rebinned.copy()
    g_rec_tau << cont

    # On master node, save results to an archive
    if mpi.is_master_node():
        with HDFArchive(output_filename,'w') as Results:
            Results['g_rec_tau'] = g_rec_tau
            Results['g_w'] = g_w

if __name__ == "__main__":
    SOM()

but I got the following error:

Starting run with 4 MPI rank(s) at : 2023-05-31 12:25:26.161903
Traceback (most recent call last):
  File "/home/ragab/test/test2AC/somAC.py", line 86, in <module>
    SOM()
  File "/home/ragab/test/test2AC/somAC.py", line 53, in SOM
    g_l_cut = cut_coefficients(g_l, n_l)
  File "/home/ragab/test/test2AC/somAC.py", line 11, in cut_coefficients
    g_cut = GfLegendre(indices = [0], beta = g_l.beta, n_points = n_remaining_coeffs)
AttributeError: 'Gf' object has no attribute 'beta'

[the-hampel]

Dear @abdelghany0, the script is for a quite old version of triqs. With the rework of the meshes that happened a few years ago the beta of a Gf is accessed within the mesh of a Gf as Gf.mesh.beta:

mesh = MeshLegendre(S= 'Fermion', beta=40, max_n=20)
Gl = Gf(mesh = mesh, target_shape=[])
print(Gl.mesh.beta)
>>> 40

Hence, you have to change the script slightly to g_l.mesh.beta to work again. The rest of the script looks okay to me, but let me know if there are further problems.

Best, Alex

[abdelghany0]

Hi Alex, After modifying the script I got: AttributeError: 'GfLegendre' object has no attribute 'N1'

best,

Ragab

[the-hampel]

Ah I think this is supposed to be the dimension of the Gf. Changing g_cut.N1 to g_cut.target_shape[0] should be the new correct way of doing this.

[abdelghany0]

Does it usually take long time to finish the analytic continuation? One hour so far, and still running!!

[the-hampel]

I do not have any experience with the SOM analytic continuation. Typically I use MaxEnt or Pade which takes less long. But it is definitely time consuming.

[Wentzell]

You are using https://github.com/krivenko/som, correct? Maybe @krivenko can comment on the typical run-times?

[krivenko]

@abdelghany0 Let us close this issue and continue discussing in https://github.com/krivenko/som/issues/7

Also, please refrain from opening identical issues in multiple repositories. This repository is devoted to the TRIQS library. Issues pertaining to the applications, such as pomerol2triqs and SOM, should be reported in their respective repositories.

[abdelghany0]

@krivenko I realized that after reporting the issue here. However, I reported the issue in SOM repository, where you graciously resolved it.

Plese close the issue here.