solver result depends on how the blocks of the GF are named

[richtma93]

I just recognized, that the result of the ctint solver is effected by which block names are chosen
for the Green's function.
Its seems to me like if the plausible ('correct') solution is just achieved if the block_names are either 'up','dn' or 'down'. If one or more blocks are called something like 'ud', 'up_0' , dn_0' ... or seemingly any thing else but the previously mentioned something goes wrong.

And the amount how wrong it is, depends on how the interaction U_Hubbard and the frequency range is chosen (and the more blocks are called 'wrongly').

Here is a script that should reproduce the Problem in it's extremest case (large U, large frequency range and both block names are neither 'up','dn' nor 'down'):
In this case the outcomming Sigma is basically just a constant.

from pytriqs.lattice import BravaisLattice, BrillouinZone
from pytriqs.lattice.tight_binding import TightBinding, dos, energies_on_bz_path
import trilex.tools as tools
import trilex.pytools as pytools
from pytriqs.gf import *
from pytriqs.operators import n
from triqs_ctint import Solver
%matplotlib notebook
from ipywidgets import *
import matplotlib.pyplot as plt 
from pytriqs.plot.mpl_interface import oplot, plt

plt.figure()
for block_names in [['up','dn'],['up_0','dn_0']]:
    gf_struct = [[bn,[0]] for bn in block_names]


    BL = BravaisLattice(units = [(1, 1) , (-1, 1) ], orbital_positions= [ (0, 0)] )

    t = -1.0               # nearest neighbor hopping
    n_orb = BL.n_orbitals
    hop = { (1,0)  :  [[t]],
            (-1,0) :  [[t]],
            (0,1)  : [[t]],
            (0,-1) :[[t]]
            }
    TB = TightBinding(BL,hop)

    U_Hubbard = 4
    filling = 0.5
    beta = 30
    n_k = 5
    n_iw = 180 #18
    iw_mesh = MeshImFreq(beta,'Fermion',n_iw)
    k_mesh = MeshBrillouinZone(BrillouinZone(BL),n_k)
    k_iw_mesh = MeshProduct(k_mesh,iw_mesh)


    G = BlockGf(mesh = k_iw_mesh, gf_struct = gf_struct)
    G0 = tools.make_local_gf(pytools.make_G0(G, TB))


    h_int = U_Hubbard*n(block_names[0],0)*n(block_names[1],0)

    S = Solver(beta = beta, 
               gf_struct = gf_struct, 
               n_iw = n_iw, 
               n_tau = 10001)

    solve_params_dmft = { 
        'h_int' : h_int,
        'n_cycles' : 54000,
        'length_cycle' : 100, 
        'n_warmup_cycles' : 1000,
        'random_seed' : 1234,
        'delta' : 0.1
    }


    for bn in block_names:
        S.Sigma_iw[bn] = U_Hubbard*filling
        S.G0_iw[bn] = G0[bn]

    S.solve(**solve_params_dmft)

    plt.subplot(2,1,1)
    oplot(S.Sigma_iw[block_names[0]].imag)
    plt.subplot(2,1,2)
    oplot(S.Sigma_iw[block_names[0]].real)

• Put an X between the brackets on this line if you have done all of the following:
  • Checked that your issue isn't already filed: https://github.com/issues?utf8=%E2%9C%93&q=is%3Aissue+user%3Atriqs

Description

[Description of the issue]

Steps to Reproduce

1. [First Step]
2. [Second Step]
3. [and so on...]

Expected behavior: [What you expect to happen]

Actual behavior: [What actually happens]

Reproduces how often: [What percentage of the time does it reproduce?]

Versions

You can get this information from copy and pasting the output of

python -c "from app4triqs.version import *; show_version(); show_git_hash();"

from the command line. Also, please include the OS and what version of the OS you're running.

Additional Information

Any additional information, configuration or data that might be necessary to reproduce the issue.

[richtma93]

I just saw that this probably belongs to: #5

[richtma93]

this is from solver.py

if bl == 'up':
              alpha.append( np.array([[n_o + delta] for n_o in dens_HF ]) )
elif bl == 'dn' or bl == 'down' or bl == 'do':
              alpha.append( np.array([[n_o - delta] for n_o in dens_HF ]) )
else:
              alpha.append( np.array([[n_o] for n_o in dens_HF ]) )

it seems like delta is only added/substracted if the names are chosen right

[Wentzell]

Dear @richtma93,

Thank you for pointing this out!

Ctint has an internal degree of freedom which can be used to improve both the sign problem
as well as the average perturbation order. This parameter is deduced automatically
unless explicitly provided. As you point out, the choice of this value is currently dependent
on the choice of the block names for your impurity model.

The final converged solution however should not depend on the choice of this alpha.
Certain choices of alpha may however take much longer to converge.

I have been working on an improvement to the ctint solver that will circumvent this explicit
dependency on the block names. I will keep this issue open until this functionality is merged.

[Wentzell]

In fact my statement is not fully correct. In some cases, like in the example given in #5, the automatic
choice of alpha introduces ergodicity problems in the absence of a double vertex insertion move. Due to the lack of rank 3/4 matrix updates in the det_manip class of triqs however
those updates are currently not possible.

[richtma93]

Dear @Wentzell,
Thank you for clarification!

I see that setting alpha right automatically is not a trivial thing. For my purpose it should be doable to work this around by choosing alpha explicitly. Thanks!