Largest imaginary element of delta(infty) e.g. of the local part of G0: 1e-10, is larger than the set parameter imag_threshold 1e-13

Dear all

I successfully followed the Nio tutorial (dft_tools) for the DMFT calculation.

I have some problem to do the calculation on my systeme which is a molecule. first I did a DFT calculation to determine the energy window of the d obital I'm interested in as you can see on the attached figure Then I decided to consider two energy windows and make the calculation for each of the windows (respectively from [Emin=-9.5; Emax=0.3] and [Emin=-5 ; Emax=0.3]) considered.

Before doing the DMFT calculation I have first of all made the calculation DFT well converged the energy with EDIFF = 1.0E-9 and I determined the double-counting potential with this formula 𝜇𝐷𝐶=U(Nd−1/2)−J/2(Nd−1) where Nd is the filling of the 3d band obtained with the DFT calculation which in my case of the Manganese molecule it is 5.4 electron. So U=4.5 , J=0.7, 𝜇𝐷𝐶=19 eV and T=160K==> 𝛽=72.5*1/eV. after having made plovasp plo.cfg and python converter.py; when I launch the dmft calculation here are the errors

for 1-) [Emin=-5; Emax=0.3] ===> EWINDOW(in plo.vasp) = -0.6 5

found 1 blocks of degenerate orbitals in shell 0 block 0 consists of orbitals: up_0 down_0 Sumk to Solver: [{('up', 0): ('up_0', 0), ('up', 1): ('up_0', 1), ('up', 2): ('up_0', 2), ('up', 3): ('up_0', 3), ('up', 4): ('up_0', 4), ('down', 0): ('down_0', 0), ('down', 1): ('down_0', 1), ('down', 2): ('down_0', 2), ('down', 3): ('down_0', 3), ('down', 4): ('down_0', 4)}] GF struct sumk: [[('up', [0, 1, 2, 3, 4]), ('down', [0, 1, 2, 3, 4])]] GF struct solver: [{'up_0': [0, 1, 2, 3, 4], 'down_0': [0, 1, 2, 3, 4]}] WARNING: gf_struct should be a list of pairs [ (str,[int,...]), ...], not a dict Greens function structure is: {'up_0': [0, 1, 2, 3, 4], 'down_0': [0, 1, 2, 3, 4]} U Matrix set to: [[0. 3.6803681 3.44110429 3.6803681 4.39815951] [3.6803681 0. 4.15889571 3.6803681 3.6803681 ] [3.44110429 4.15889571 0. 4.15889571 3.44110429] [3.6803681 3.6803681 4.15889571 0. 3.6803681 ] [4.39815951 3.6803681 3.44110429 3.6803681 0. ]] Up Matrix set to: [[5.3 4.2202454 4.0607362 4.2202454 4.69877301] [4.2202454 5.3 4.5392638 4.2202454 4.2202454 ] [4.0607362 4.5392638 5.3 4.5392638 4.0607362 ] [4.2202454 4.2202454 4.5392638 5.3 4.2202454 ] [4.69877301 4.2202454 4.0607362 4.2202454 5.3 ]] Dichotomy adjustment of Chemical Potential to obtain Total Density = 7.000000 +/- 0.010000 Warning: Imaginary part in density will be ignored (3.3587185850876708e-12) Warning: Imaginary part in density will be ignored (3.3360817220610182e-12) Chemical Potential = 0.500000
5.672753 < Total Density < 7.005509 Chemical Potential found in 6 iterations : Total Density = 7.005509;Chemical Potential = 0.409874 DC for shell 0 = 19.000000 DC energy = 133.10465457554184 10 DMFT cycles requested. Starting with iteration 0. Doing iteration: 0

╔╦╗╦═╗╦╔═╗ ╔═╗ ┌─┐┌┬┐┬ ┬┬ ┬┌┐ ║ ╠╦╝║║═╬╗╚═╗ │ │ ├─┤└┬┘├┴┐ ╩ ╩╚═╩╚═╝╚╚═╝ └─┘ ┴ ┴ ┴ ┴ └─┘

Traceback (most recent call last): File "/home/DMFT_MnPc/1/simple/dmft.py", line 138, in S.solve(h_int = H, p) File "/home/triqs/install/lib/python3.8/site-packages/triqs_cthyb/solver.py", line 140, in solve solve_status = SolverCore.solve(self, params_kw) RuntimeError: .. Error occurred at Wed Sep 1 18:52:12 2021

.. Error .. calling C++ overload .. solve() -> void .. in implementation of method SolverCore.solve .. C++ error was : .. Triqs runtime error at /home/triqs/cthyb.src/c++/triqs_cthyb/solver_core.cpp : 141

Largest imaginary element of delta(infty) e.g. of the local part of G0: 1.99992e-10, is larger than the set parameter imag_threshold 1e-13 .. Error occurred on node 7

Traceback (most recent call last): File "/home/DMFT/dmft.py", line 138, in S.solve(h_int = H, p) File "/home/triqs/install/lib/python3.8/site-packages/triqs_cthyb/solver.py", line 140, in solve solve_status = SolverCore.solve(self, params_kw) RuntimeError: .. Error occurred at Wed Sep 1 18:52:12 2021

for 2-) [Emin=-9.5; Emax=0.3] ===> EWINDOW(in plo.vasp) = -5 5

ound 1 blocks of degenerate orbitals in shell 0 block 0 consists of orbitals: up_0 down_0 Sumk to Solver: [{('up', 0): ('up_0', 0), ('up', 1): ('up_0', 1), ('up', 2): ('up_0', 2), ('up', 3): ('up_0', 3), ('up', 4): ('up_0', 4), ('down', 0): ('down_0', 0), ('down', 1): ('down_0', 1), ('down', 2): ('down_0', 2), ('down', 3): ('down_0', 3), ('down', 4): ('down_0', 4)}] GF struct sumk: [[('up', [0, 1, 2, 3, 4]), ('down', [0, 1, 2, 3, 4])]] GF struct solver: [{'up_0': [0, 1, 2, 3, 4], 'down_0': [0, 1, 2, 3, 4]}] WARNING: gf_struct should be a list of pairs [ (str,[int,...]), ...], not a dict Greens function structure is: {'up_0': [0, 1, 2, 3, 4], 'down_0': [0, 1, 2, 3, 4]} U Matrix set to: [[0. 3.6803681 3.44110429 3.6803681 4.39815951] [3.6803681 0. 4.15889571 3.6803681 3.6803681 ] [3.44110429 4.15889571 0. 4.15889571 3.44110429] [3.6803681 3.6803681 4.15889571 0. 3.6803681 ] [4.39815951 3.6803681 3.44110429 3.6803681 0. ]] Up Matrix set to: [[5.3 4.2202454 4.0607362 4.2202454 4.69877301] [4.2202454 5.3 4.5392638 4.2202454 4.2202454 ] [4.0607362 4.5392638 5.3 4.5392638 4.0607362 ] [4.2202454 4.2202454 4.5392638 5.3 4.2202454 ] [4.69877301 4.2202454 4.0607362 4.2202454 5.3 ]] Dichotomy adjustment of Chemical Potential to obtain Total Density = 67.000000 +/- 0.010000 Warning: Imaginary part in density will be ignored (3.318524565686638e-12) Warning: Imaginary part in density will be ignored (3.3111927824289807e-12) Chemical Potential = 0.500000
Total Density = 8.000000 Total Density = 10.137773 0.000000 < Chemical Potential < 50.500000 5.425579 < Total Density < 10.137773 FAILURE to adjust Chemical Potential to the value 67.000000 after 101 iterations. FAILURE returning (None, None) due to failure. Starting run with 24 MPI rank(s) at : 2021-09-01 19:18:30.137226 Traceback (most recent call last): File "/home/DMFT_MnPc/2/simple/dmft.py", line 119, in S.G_iw << SK.extract_G_loc()[0] File "/home/triqs/install/lib/python3.8/site-packages/triqs_dft_tools/sumk_dft.py", line 756, in extract_G_loc G_latt = self.lattice_gf( File "/home/triqs/install/lib/python3.8/site-packages/triqs_dft_tools/sumk_dft.py", line 587, in lattice_gf (idmat[ibl] mu) - (idmat[ibl] self.h_field (1 - 2 ibl)) TypeError: unsupported operand type(s) for *: 'complex' and 'NoneType'

Thanks

Hi Olivier,

for your first example everything looks good. Cthyb default threshold for imaginary parts in the local part of G0 is rather stringent. Your output says that the largest element is 1E-10 which is not too bad. In this case (if you do not want to treat spin-orbit coupling) you can safely add the following solve parameter after line 78 p['imag_threshold']=1e-08 (https://triqs.github.io/cthyb/unstable/reference/solve_parameters.html) . To get rid of those small elements you can either use the "delta interface" of cthyb and feed Hloc and the hybridization separately or increase your Vasp convergence even further (higher ENCUT and kpts). This should get you started.

For your second example I am not sure. Why does the converter believe, that you have 67 electrons in the given energy window? Your window -5 to 5 eV looks good to me. Does the electron count makes sense to you? How many electron are in the window -5 to Fermi level in DFT? Furthermore I would suggest you run the normal converter mode with:

HK = False
COMPLEMENT = False

However, this should not lead to your error. The error simply tells you it cannot find 67 electrons in your system ;-)

P.S. If you add raw text output, or source code in your comments please use the appropiate markings for markdown. Preview your message to make sure that it is readable: https://guides.github.com/features/mastering-markdown/ .

Hi, by adding p['imag_threshold']=1e-07 solve the problem for both cases.

indeed for this system NELECT=197 electrons. For the interval [-5 5] if we take the totality of charge including the other uncorrelated atoms, so it's possible to have 66 electrons. But if I use the vannier projection on the correlated atom, then for this windows I have 8 electrons in the wannier subspace(which is not the concept used in triqs i think)

Hi, that is good. So all calculations are running now?

The number of electrons is determined from the DFT bands in this window. The converter counts the electrons in the given window in the KS states up to the Fermi level. Hence, you can directly determine those from the PROCAR file of Vasp. If all your problems are solved, please consider closing the issue.

Alex

Hi Alex,

Yes, by removing proj_or_hk = 'hk' in converter.py, and setting HK = False COMPLEMENT = False the DMFT calculation turned out fine.

1-) in modified n_orbs = SK.proj_mat_csc.shape[2] to n_orbs = SK.proj_mat.shape[2] . But I still have a problem with the script local_lattice_GP.py. when I run it I get this error(please do you have an idea of how to solve it?):

the first error message was this:

$ mpirun -np 24 python st.py ('offset', 20) Starting run with 24 MPI rank(s) at : 2021-09-11 01:35:35.235327 Traceback (most recent call last): File "st.py", line 80, in add_g_ik << SK.downfold(ik, 0, bname, G_latt_KS[bname], gf, shells='csc', ir=None) File "/home/triqs/install/lib/python3.8/site-packages/triqs_dft_tools/sumk_dft.py", line 342, in downfold gf_downfolded.from_L_G_R( File "/home/triqs/install/lib/python3.8/site-packages/triqs/gf/gf.py", line 673, in from_L_G_R assert L.shape[1] == G.target_shape[0], "Dimension mismatch between L and G" AssertionError: Dimension mismatch between L and G

so I modified the following two lines:

ikarray = numpy.array(list(range(SK.n_k))) to ikarray = numpy.array(range(SK.n_k)) block_structure = [list(range(n_orbs)) for sp in spn] to block_structure = [range(n_orbs) for sp in spn]

an the a got this

mpirun -np 24 python st.py ('offset', 20) Starting run with 24 MPI rank(s) at : 2021-09-11 01:48:24.617819 Traceback (most recent call last): File "st.py", line 71, in block_list=glist(), make_copies=False) File "st.py", line 67, in glist = lambda: [GfImFreq(indices=inner, mesh=mesh) File "st.py", line 67, in glist = lambda: [GfImFreq(indices=inner, mesh=mesh) File "/home/triqs/install/lib/python3.8/site-packages/triqs/gf/backwd_compat/gf_imfreq.py", line 73, in init delegate(self, kw) File "home/triqs/install/lib/python3.8/site-packages/triqs/gf/backwd_compat/gf_imfreq.py", line 66, in delegate super(GfImFreq, self).init( File "/home/triqs/install/lib/python3.8/site-packages/triqs/gf/gf.py", line 198, in init delegate(self, kw) File "/home/triqs/install/lib/python3.8/site-packages/triqs/gf/gf.py", line 141, in delegate assert isinstance(indices, (type(None), list, GfIndices)), "Type of indices incorrect : should be None, Gfindices, list of str, or list of list of str" AssertionError: Type of indices incorrect : should be None, Gfindices, list of str, or list of list of str

Or this

('offset', 7) Starting run with 144 MPI rank(s) at : 2021-09-10 10:47:01.006994 Traceback (most recent call last): File "/home/DMFT/g.py", line 79, in add_g_ik << SK.downfold(ik, 0, bname, G_latt_KS[bname], gf, shells='csc', ir=None) File "/home/triqs/install/lib/python3.8/site-packages/triqs_dft_tools/sumk_dft.py", line 340, in downfold projmat = self.proj_mat_csc[ik, isp, :, 0:n_orb] TypeError: 'NoneType' object is not subscriptable

Thanks

Hi,

1) If you change the projector to the normal mode (Hk=False) than you have to change the line:

add_g_ik << SK.downfold(ik, 0, bname, G_latt_KS[bname], gf, shells='csc', ir=None)

add_g_ik << SK.downfold(ik, 0, bname, G_latt_KS[bname], gf, shells='corr')

With shells='csc' it will fail, because it expects different shapes of the projectors. Try that instead your other changes.

2) I think there is a misunderstanding. You choose to use the Vasp projector interface to construct your low energy Hamiltonian to perform DMFT calculations with triqs. But you can also directly use your wannier90 output to do the same. Especially if you already used wannier90 to construct such Hamiltonian. The interface needs only the wannier90_hr.dat file to construct a triqs readable h5 input file: https://triqs.github.io/dft_tools/unstable/guide/conv_W90.html . If you do so, everything works more or less the same. This uses triqs in exactly the same way., plus you have access to all three things you ask for in both interfaces.

-- Alex

Dear Alex,

I have done what you suggested but it doesn't always work and I have this error now :

Starting run with 24 MPI rank(s) at : 2021-09-14 20:30:52.372143
Traceback (most recent call last):
  File "gl.py", line 80, in <module>
    add_g_ik << SK.downfold(ik, 0, bname, G_latt_KS[bname], gf, shells='corr')
  File "/home/triqs/install/lib/python3.8/site-packages/triqs_dft_tools/sumk_dft.py", line 342, in downfold
    gf_downfolded.from_L_G_R(
  File "/home/triqs/install/lib/python3.8/site-packages/triqs/gf/gf.py", line 676, in from_L_G_R
    assert L.shape[0] == self.target_shape[0], "Dimension mismatch between L and self"
AssertionError: Dimension mismatch between L and self

MODIFICATIONS:

1-HK = False ; COMPLEMENT = False 2- n_orbs = SK.proj_mat.shape[2] 3-) add_g_ik << SK.downfold(ik, 0, bname, G_latt_KS[bname], gf, shells='corr') (if i remove this line the analysis will still correct?)

Thanks

Hi,

1) I don't have a clear idea what goes wrong now. Can you put your exact script somewhere, so I can check it? In principle this just tells us, that the projectors in SK.proj_mat have the wrong shape. Not sure why. The changes look good to me. And I cannot say if this line is necessary or not. For that I need to see the whole script.

2) If you want to analyze a 5x5 Gf object you should use matrix valued continuation. You can do that a bit simpler for BlockGf object:

results = {}
for name, giw in G_iw:
     tm = TauMaxEnt()
     tm.set_G_iw(giw)
     tm.set_error(1.e-3)
     results[name] = tm.run().data

see also: https://triqs.github.io/maxent/latest/guide/blockgf.html . In general the documentation of triqs/maxent should help you a lot here.

The Matsubara self energies are written to the h5 archive in your calculation I guess? Then you can load them afterwards:

with HDFArchive('your_results.h5','r') as h5:
    Sigma_iw = ar['DMFT_results']['Iterations']['Sigma_itX']

The dynamic crystal field is not something that we provide. You have to calculate that by yourself from Sigma_iw.

the local_lattice_GF.py is

from itertools import *
import numpy as np
import triqs.utility.mpi as mpi
from h5 import *
from triqs.gf import *
from triqs_dft_tools.sumk_dft import *
from triqs_dft_tools.sumk_dft_tools import *
from triqs.operators.util.hamiltonians import *
from triqs.operators.util.U_matrix import *
from triqs_cthyb import *
import warnings
warnings.filterwarnings("ignore", category=FutureWarning)

filename = 'Mn'
SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False)

beta = 101.79

# We analyze the block structure of the Hamiltonian
Sigma = SK.block_structure.create_gf(beta=beta)

SK.put_Sigma([Sigma])
G = SK.extract_G_loc()
SK.analyse_block_structure_from_gf(G, threshold = 1e-5)

# Setup CTQMC Solver
n_orb = SK.corr_shells[0]['dim']
spin_names = ['up','down']
orb_names = [i for i in range(0,n_orb)]

# Print some information on the master node
iteration_offset = 0
Sigma_iw = 0
if mpi.is_master_node():
    ar = HDFArchive(filename+'.h5','a')
    if not 'DMFT_results' in ar: ar.create_group('DMFT_results')
    if not 'Iterations' in ar['DMFT_results']: ar['DMFT_results'].create_group('Iterations')
    if 'iteration_count' in ar['DMFT_results']: 
        iteration_offset = ar['DMFT_results']['iteration_count']+1
        print(('offset',iteration_offset))
        Sigma_iw = ar['DMFT_results']['Iterations']['Sigma_it'+str(iteration_offset-1)]
        SK.dc_imp = ar['DMFT_results']['Iterations']['dc_imp'+str(iteration_offset-1)]
        SK.dc_energ = ar['DMFT_results']['Iterations']['dc_energ'+str(iteration_offset-1)]
        SK.chemical_potential = ar['DMFT_results']['Iterations']['chemical_potential'+str(iteration_offset-1)].real

iteration_offset = mpi.bcast(iteration_offset)
Sigma_iw = mpi.bcast(Sigma_iw)
SK.dc_imp = mpi.bcast(SK.dc_imp)
SK.dc_energ = mpi.bcast(SK.dc_energ)
SK.chemical_potential = mpi.bcast(SK.chemical_potential)

SK.put_Sigma(Sigma_imp = [Sigma_iw])

ikarray = numpy.array(list(range(SK.n_k)))

# set up the orbitally resolved local lattice greens function:
n_orbs = SK.proj_mat.shape[2]
spn = SK.spin_block_names[SK.SO]
mesh = Sigma_iw.mesh
block_structure = [list(range(n_orbs)) for sp in spn]
gf_struct = [(spn[isp], block_structure[isp])
         for isp in range(SK.n_spin_blocks[SK.SO])]
block_ind_list = [block for block, inner in gf_struct]
glist = lambda: [GfImFreq(indices=inner, mesh=mesh)
    for block, inner in gf_struct]

G_latt_orb = BlockGf(name_list=block_ind_list,
             block_list=glist(), make_copies=False)
G_latt_orb.zero()

for ik in mpi.slice_array(ikarray):
    G_latt_KS = SK.lattice_gf(ik=ik, beta=beta)
    G_latt_KS *= SK.bz_weights[ik]
    for bname, gf in G_latt_orb:
        add_g_ik = gf.copy()
        add_g_ik.zero()
        add_g_ik << SK.downfold(ik, 0, bname, G_latt_KS[bname], gf, shells='corr')
        gf << gf + add_g_ik

G_latt_orb << mpi.all_reduce(
                mpi.world, G_latt_orb, lambda x, y: x + y)

mpi.barrier()

if mpi.is_master_node(): 
    ar['DMFT_results']['Iterations']['G_latt_orb_it'+str(iteration_offset-1)] = G_latt_orb
if mpi.is_master_node(): del ar

The error message:

Starting run with 24 MPI rank(s) at : 2021-09-18 21:02:04.119798
Traceback (most recent call last):
  File "G_local.py", line 80, in <module>
    add_g_ik << SK.downfold(ik, 0, bname, G_latt_KS[bname], gf)
  File "/home/triqs/install/lib/python3.8/site-packages/triqs_dft_tools/sumk_dft.py", line 342, in downfold
    gf_downfolded.from_L_G_R(
  File "/home/triqs/install/lib/python3.8/site-packages/triqs/gf/gf.py", line 676, in from_L_G_R
    assert L.shape[0] == self.target_shape[0], "Dimension mismatch between L and self"
AssertionError: Dimension mismatch between L and self

Hi, I think I found the problem in your script G_local.py, and could reproduce the problem. This line is the problem:

n_orbs = SK.proj_mat.shape[2]

proj_mat.shape[2] is not the number of orbitals, it is the number of sites. You can check this in the documentation: https://triqs.github.io/dft_tools/unstable/reference/h5structure.html#groups-and-their-formats which says:

numpy.array.complex, dim [n_k,SP+1-SO,n_corr_shells,max(corr_shell[‘dim’]),max(n_orbitals)]

Index 2 will point to n_corr_shells, not number of orbitals. To obtain the number of orbitals in the correlated shell the safest option is to replace the line by:

n_orbs = SK.corr_shells[0]['dim']

where 0 marks the site you want to target. Changing this line will make your script work. I changed the line and tried it with some of my files and then it worked fine.

Sorry but such kind of problems take a long time to find. Please try to provide a minimal example in the future and only provide the necessary script. I first tried to look at your first script dmft.py you posted, until I realized the problem is only in the second script.

Dear Alex,

i) Thanks for your suggestion it solved the problem. in a previous post you suggested to do the analysis like this

results = {}
for name, giw in G_iw:
     tm = TauMaxEnt()
     tm.set_G_iw(giw)
     tm.set_error(1.e-3)

i didn't understand well and how to insert it in the script maxent.py of the example NiO i simply modified this script like this

print((G_latt['up'][0,0]))
dxy_orbs = [0]
dyz_orbs = [1]
dz2_orbs = [2]
dxz_orbs = [3]
dx2y2_orbs = [4]
orbs = [dxy_orbs, dyz_orbs, dz2_orbs, dxz_orbs, dx2y2_orbs]

did i make a mistake?

ii) regarding the second point you mentioned: the representation of imaginary part of Matsubara self energies, and that it was only necessary to load

with HDFArchive('your_results.h5','r') as h5:
    Sigma_iw = ar['DMFT_results']['Iterations']['Sigma_itX']

indeed I have the Sigma_itX files(where X=0;...;19) when I try to visualize the last file Sigma_it19 with HDFView for instance, but when i try to visualize the last file (which has 6 columns i suppose that the 5 columns are the one of each of the 5 orbitals) it looks more like a density of state curve than the one of a Matsubara self energies (because i googled to see what this quantity could look like) and moreover i don't know if it's the imaginary part or not. please can you enlighten me on this point?(or suggest me a script to get this quantity?)

Indeed the second part is more of an analysis. If I have to post this question on maxent, tell me I'm closing this ticket because the initial problem is solved!

Best regards,

Hi, Although quite complicated your setup for the analytic continuation looks okay, but you will drop all off-diagonal elements in the analytic continuation! Check also the convergence of your G_imp / cthyb run with respect to the number of cycles to make sure you obtain a reasonable result.

Regarding your second point you can maybe just open a discussion here: https://github.com/TRIQS/dft_tools/discussions. This is not really an issue and maybe also someone else is able to help you with that. In principle I would encourage you to use the python oplot interface in triqs (https://triqs.github.io/triqs/unstable/documentation/manual/triqs/plotting_protocols/plotting/plotting.html?highlight=oplot) . Also there is an example how the self-energy should look like for SrVO3 here: https://triqs.github.io/dft_tools/unstable/tutorials/srvo3.html#the-dmft-calculation .

I will close this issue now, as the original problem is solved. Feel free to open a new discussion / issue if necessary.

Best, Alex

TRIQS / dft_tools

Largest imaginary element of delta(infty) e.g. of the local part of G0: 1e-10, is larger than the set parameter imag_threshold 1e-13 #178