IndexError: index 5 is out of bounds for axis 1 with size 5

[quantumable]

I did some DMFT calculations using dft_tools but I'm a bit stuck on the maxent analysis.

indeed i made the calculation for two correlated atoms which i want to plot each of the orbitals d(dxy,dxy,dz2,dxz,dx2-y2) for atom1 and 2. i used the model of the maxent script present here MAXENT.

i attach here the example of my maxent file. and the error message please do you have an idea for these two atoms? and when the calculation is spin polarized?

ERROR:

Traceback (most recent call last):
  File "max.py", line 40, in <module>
    gf = gf + G_latt['up'][iO,iO]
  File /TRIQS/install/lib/python3.8/site-packages/triqs/gf/gf.py", line 359, in __getitem__
    dat = self._data[ self._rank*(slice(0,None),) + key_tpl ]
IndexError: index 5 is out of bounds for axis 1 with size 5

``

[script.zip](https://github.com/TRIQS/maxent/files/7720110/script.zip)

[the-hampel]

Hi, this is not an issue in MaxEnt. your Green function object has a shape (5,5) so you can access indices 0,1,2,3,4. But in your script you are at some point looping over something with one element being d1xy_orbs=[5] . Then you loop also over this list for iO in orb: where orb=[5] . So the script will try to access G_latt['up'][5,5] which does not exist. So this cannot work. What are the list entries

d1xy_orbs = [5]
d1yz_orbs = [6]
d1z2_orbs = [7]
d1xz_orbs = [8]
d1x2y2_orbs = [9]

for? Possible solution: Why are you not directly looping over range(5)? I.e. change line 28 to

orbs = list(range(5))

and then remove the additional loop for iO in orb and directly access G_latt['up'][orb, orb]?

You can also directly continue the whole matrix values Gf G_latt['up'] : https://triqs.github.io/maxent/latest/guide/elementwise.html .

[quantumable]

indeed I have two correlated atoms(see plo.cfg file) and I thought that Glatt was written like this(d1xy,d1yz,d1z2,d1xz,d1x2-y2, d2xy,d2yz,d2z2,d2xz,d2x2-y2) where the index 1 and 2 are those of the atoms 1 and 2.

I want to plot for each of these atoms the different d orbitals. following your instructions I can't do it.

or I didn't understand your instructions!!

and for the polarized spin calculation how to do it? what order in the GF?

from triqs.gf import *
from h5 import *
from triqs_maxent import *

filename = 'MnPc'

ar = HDFArchive(filename+'.h5','a')
if 'iteration_count' in ar['DMFT_results']:
    iteration_offset = ar['DMFT_results']['iteration_count']+1
    G_latt = ar['DMFT_results']['Iterations']['G_latt_orb_it'+str(iteration_offset-1)]

tm = TauMaxEnt(cost_function='bryan', probability='normal')
print((G_latt['up'][0,0]))
dxy_orbs = [0]
dyz_orbs = [1]
dz2_orbs = [2]
dxz_orbs = [3]
dx2y2_orbs = [4]  

d1xy_orbs = [5]
d1yz_orbs = [6]
d1z2_orbs = [7]
d1xz_orbs = [8]
d1x2y2_orbs = [9]

#orbs = [dxy_orbs, dyz_orbs, dz2_orbs, dxz_orbs, dx2y2_orbs, d1xy_orbs, d1yz_orbs, d1z2_orbs, d1xz_orbs, d1x2y2_orbs]
orbs = list(range(5))
#orbs = [t2g_orbs]

for orb in orbs:

    print('\n'+str(orb[0])+'\n')

    gf = 0*G_latt['up'][0,0]
   # for iO in orb:
     #   gf = gf + G_latt['up'][iO,iO]
    tm.set_G_iw(gf)
    tm.omega =LinearOmegaMesh(omega_min=-10, omega_max=10, n_points=301)
    tm.alpha_mesh = LogAlphaMesh(alpha_min=0.01, alpha_max=20000, n_points=100)

    tm.set_error(1.e-3)
    result=tm.run()
    result.get_A_out('LineFitAnalyzer')

    if 'iteration_count' in ar['DMFT_results']:
        iteration_offset = ar['DMFT_results']['iteration_count']+1
        for oo in orb:
            ar['DMFT_results']['Iterations']['G_latt_orb_w_o'+str(oo)+'_it'+str(iteration_offset-1)] = result.analyzer_results['LineFitAnalyzer']['A_out']
        ar['DMFT_results']['Iterations']['w_it'+str(iteration_offset-1)] = result.omega
del ar

[the-hampel]

Yes I see that in your plo.cfg file you have two shells. I just tried to explain to you what the error means. I don't want to do your whole calculation. I can only tell you that your Green functions has only 5 orbitals.

I also see that in your dmft script you never loops over the two shells. You have to solve two impurity problems. One for each shell. So you have to add another for loop over the shells in your dmft script. So for example

S.G_iw << SK.extract_G_loc()[0]

should be change to

S[icrsh].G_iw << SK.extract_G_loc()[icrsh]

Otherwise you only solve one impurity problem for your first correlated shell. Also the script Mn_local_lattice_GF.py is only written to use the first correlated shell. Once you have that you can also calculated the lattice GF for each shell. The scripts will never append the orbitals from one shell to the other in the way you wrote the script. For example in your Mn_local_lattice_GF.py script the line

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

will extract the 0 shell downfolded lattice Gf at kpoint ik. the second argument of the function is the shell index.

I hope that helps.

[quantumable]

I understand the philosophy of your code. But for someone new, it's very very difficult to use it! And to be honest, I don't necessarily understand what you say every time.

by changing S.G_iw << SK.extract_G_loc()[0] by S[icrsh]. G_iw << SK.extract_G_loc()[icrsh]

NameError : le nom 'icrsh' n'est pas défini
Retraçage (dernier appel le plus récent) :
Fichier "dmft.py", ligne 117, dans <module>
S[icrsh]. G_iw << SK.extract_G_loc()[icrsh]
NameError : le nom 'icrsh' n'est pas défini

it is the same type of modification that it is necessary to bring for the calculations SPIN polarized? this example #80 talks about a similar case as mine of two correlated atoms and initializing the magnetic moment but it doesn't use vasp. you had suggested me to make first a paramagnetic calculation and then restart the calculation with:SK.deg_shells = [[]]and initializing the magnetic moment as follows: Sigma_iw['up'] << 1.2Sigma_iw['up'] Sigma_iw['down'] << 0.8Sigma_iw['down'] . for the suggestions that you had made to me, the calculation turns just with:cSK.deg_shells = [[]] ; by adding Sigma_iw['up'] << 1.2Sigma_iw['up'] Sigma_iw['down'] << 0.8Sigma_iw['down'] I have this error initial moment. the fact of not initializing the magnetic moment and of putting just SK.deg_shells = [[]] would not change the results of the calculation?

2-) and for the script Mn_local_lattice_GF.py I did not understand what I should change.

@cmtcl18 I saw that you wanted to do a spin-polarized self consistent calculation with several correlated atoms. did you succeed? can you help me with the script? i use VASP

[quantumable]

Dear Alex,

Please could you give me feedback on my message so that I know if I should wait for the resolution or move on and do something else. thank you!

best regards.

[the-hampel]

Hi, I am not quite sure how to help you. I think my answer was quite detailed, and targets you to the right problems. You also cannot demand any service for something that you do not pay for in the way you did the last messages. I am happy to help you working with triqs, but many of the things you ask for have nothing to do with triqs. These are very fundamental questions for DMFT or python, and I simply to not have the time to answer all of it.

If you have troubles writing your own DMFT loop you can always use ready to use DMFT programs. These also exist for triqs, e.g. solid_dmft (https://flatironinstitute.github.io/solid_dmft) or dcore (https://issp-center-dev.github.io/DCore). They work for most problems and utilize a simple to use input file.

To answer your question above regarding spin-polarized: You can also skip setting the initial moment and hope the solution gets magnetic by itself.

1) Again: My change suggestions are only part of the solution. You have to also implement a for loop in python of the shells. You inserted now in the code icrsh. How should python know what this variable means? You have to define it, by making a for loop over the shells, with icrsh the loop variable.

2) Mn_local_lattice_GF.py: I think I was quite clear that the function you called explicitly has the shell as an argument. In your script you past 0 as variable. So you only obtain the local lattice Gf for the 0 shell. Not for any other.

[quantumable]

Dear Alex,

Thank you for your message, I am a native French speaker, perhaps the translation can cause problems in understanding and I apologize.

I was saying that for the case of two correlated atoms I made the modification you suggested: S.G_iw << SK.extract_G_loc()[0] by. S[icrsh]. G_iw << SK.extract_G_loc()[icrsh] without success and that I did not know how to solve it.

Best

[the-hampel]

Yes because you need to write a python for loop that loops over the shells. Something like

for icrsh in range(SK.n_inequiv_shells):
    S[icrsh].G_iw << SK.extract_G_loc()[icrsh]
    ... do more

[quantumable]

Dear Alex,

I have made the changes and here is the script I am shooting. Please can you look at it and tell me if you have any comments or changes?

One-shot

from itertools import *
 import numpy as np
 import triqs.utility.mpi as mpi
 from h5 import *
 from triqs.gf.tools import delta, descriptors,gf_fnt
 from triqs.gf import *
 import sys, triqs.version as triqs_version
 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 triqs_cthyb.version as cthyb_version
 import triqs_dft_tools.version as dft_tools_version
 import warnings
 warnings.filterwarnings("ignore", category=FutureWarning)
 filename = 'Mn'
 SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False)
 beta = 223.16  
 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)
 S = []
 h_int = [None] * SK.n_inequiv_shells
 H = [None] * SK.n_inequiv_shells
 for icrsh in range(SK.n_inequiv_shells):
     n_orb = SK.corr_shells[icrsh]['dim']
     orb_names = [i for i in range(icrsh,n_orb)]
     l = SK.corr_shells[icrsh]['l']
     spin_names = ["up", "down"]
     gf_struct = [(block, indices) for block, indices in SK.gf_struct_solver[icrsh].items()] 
     S.append(Solver(beta=beta, gf_struct=gf_struct))
 mpi.report('Sumk to Solver: %s'%SK.sumk_to_solver)
 mpi.report('GF struct sumk: %s'%SK.gf_struct_sumk)
 mpi.report('GF struct solver: %s'%SK.gf_struct_solver)
 H = Operator() 
 U = 4.5
 J = 0.7
 Umat, Upmat = U_matrix_kanamori(n_orb=n_orb, U_int=[U][icrsh], J_hund=[J][icrsh]) 
 h_int[icrsh] = h_int_density(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[icrsh], U=Umat, Uprime=Upmat) 
 mpi.report('Greens function structure is: %s '%gf_struct)
 mpi.report('U Matrix set to:\n%s'%Umat)
 mpi.report('Up Matrix set to:\n%s'%Upmat)
 p = {}
 p["max_time"] = -1
 p["random_name"] = ""
 p["random_seed"] = 123 * mpi.rank + 567
 p["length_cycle"] = 100
 p["n_warmup_cycles"] = 8000
 p["n_cycles"] = 200000
 p["fit_max_moment"] = 4
 p["fit_min_n"] = 30
 p['imag_threshold']=4e-10
 p["fit_max_n"] = 50
 p["perform_tail_fit"] = True
 dc_type = 0
 dc_value = 42.92
 n_iterations = 15
 iteration_offset = 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 not 'DMFT_input' in ar: ar.create_group('DMFT_input')
     if not 'Iterations' in ar['DMFT_input']: ar['DMFT_input'].create_group('Iterations')
     if not 'code_versions' in ar['DMFT_input']: ar['DMFT_input'].create_group('code_versio\
 ns')
     ar['DMFT_input']['code_versions']["triqs_version"] = triqs_version.version
     ar['DMFT_input']['code_versions']["triqs_git"] = triqs_version.git_hash
     ar['DMFT_input']['code_versions']["cthyb_version"] = cthyb_version.version
     ar['DMFT_input']['code_versions']["cthyb_git"] = cthyb_version.triqs_cthyb_hash
     ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.version
     ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.triqs_dft_tools_hash
     if 'iteration_count' in ar['DMFT_results']:
         iteration_offset = ar['DMFT_results']['iteration_count']+1
         S[icrsh].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
     ar['DMFT_input']["dmft_script_it"+str(iteration_offset)] = open(sys.argv[0]).read()
 iteration_offset = mpi.bcast(iteration_offset)
 S[icrsh].Sigma_iw =  mpi.bcast(S[icrsh].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.symm_deg_gf(S[icrsh].Sigma_iw,ish=icrsh)
 SK.put_Sigma(Sigma_imp = [S[icrsh].Sigma_iw])
 SK.calc_mu(precision=0.01)
 S[icrsh].G_iw << SK.extract_G_loc()[icrsh]
 SK.symm_deg_gf(S[icrsh].G_iw, ish=icrsh)
 if iteration_offset == 0:
     dm = S[icrsh].G_iw.density()
     SK.calc_dc(dm, U_interact=U, J_hund=J, orb=icrsh, use_dc_formula=dc_type,use_dc_value=dc_value)
     S[icrsh].Sigma_iw << SK.dc_imp[icrsh]['up'][0,0]
  mpi.report('%s DMFT cycles requested. Starting with iteration %s.'%(n_iterations,iteration_offset))
  for it in range(iteration_offset, iteration_offset + n_iterations):    
     mpi.report('Doing iteration: %s'%it)
     S[icrsh].G0_iw << inverse(S[icrsh].Sigma_iw + inverse(S[icrsh].G_iw))
     S[icrsh].solve(h_int = h_int[icrsh], **p)
     if mpi.is_master_node(): 
         ar['DMFT_input']['Iterations']['solver_dict_it'+str(it)] = p
         ar['DMFT_results']['Iterations']['G_0_it'+str(it)] = S[icrsh].G0_iw
         ar['DMFT_results']['Iterations']['Gimp_it'+str(it)] = S[icrsh].G_iw
         ar['DMFT_results']['Iterations']['Gtau_it'+str(it)] = S[icrsh].G_tau
         ar['DMFT_results']['Iterations']['Sigma_uns_it'+str(it)] = S[icrsh].Sigma_iw
     dm = S[icrsh].G_iw.density()
     SK.calc_dc(dm, U_interact=U, J_hund=J, orb=icrsh, use_dc_formula=dc_type,use_dc_value=dc_value)
     SK.symm_deg_gf(S[icrsh].Sigma_iw,ish=icrsh)
     SK.put_Sigma(Sigma_imp=[S[icrsh].Sigma_iw])
     SK.calc_mu(precision=0.01)
     S[icrsh].G_iw << SK.extract_G_loc()[icrsh]
     for sig,gf in S[icrsh].G_iw:
         mpi.report("Orbital %s density: %.6f"%(sig,dm[sig][0,0].real))
     mpi.report('Total charge of Gloc : %.6f'%S[icrsh].G_iw.total_density().real)
     if mpi.is_master_node():
         ar['DMFT_results']['iteration_count'] = it
         ar['DMFT_results']['Iterations']['Sigma_it'+str(it)] = S[icrsh].Sigma_iw
         ar['DMFT_results']['Iterations']['Gloc_it'+str(it)] = S[icrsh].G_iw
         ar['DMFT_results']['Iterations']['G0loc_it'+str(it)] = S[icrsh].G0_iw
         ar['DMFT_results']['Iterations']['dc_imp'+str(it)] = SK.dc_imp
         ar['DMFT_results']['Iterations']['dc_energ'+str(it)] = SK.dc_energ
         ar['DMFT_results']['Iterations']['chemical_potential'+str(it)] = SK.chemical_potential
 if mpi.is_master_node(): del ar

and the. Mn_local_lattice_GF.py file

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 = '2MnPc'
SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False)

beta = 223.16 

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-4)

for icrsh in range(SK.n_inequiv_shells):
    n_orb = SK.corr_shells[icrsh]['dim']
    orb_names = [i for i in range(icrsh,n_orb)]
    l = SK.corr_shells[icrsh]['l']
    spin_names = ["up", "down"]

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)))

n_orbs = SK.corr_shells[icrsh]['dim']
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, icrsh, bname, G_latt_KS[bname], gf, shells='corr', ir=None)
        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-hampel]

Hi, I had a quick look and your changes to the dmft.py script are very hard to read, because it is not formatted. Why is the local_lattice_GF.py script correctly formatted and the other not? Anyway, I noticed that you added the loop over shells not correct. You have to loop over the impurities within the DMFT loop:

G_loc_all = SK.extract_G_loc()
for it in range(iteration_offset, iteration_offset + n_iterations):
    for icrsh in range(SK.n_inequiv_shells):
        S[icrsh].G_iw << G_loc_all[icrsh]
        ...

If you sent the script correctly formatted I can help you better. But does this make sense to you that the loop over impurities is in the DMFT iteration loop?

Also, if you have the feeling that this is all too cumbersome I would like to emphasize that we also have a DFT+DMFT code ready to be used, for various scenarios (the tutorials are only mend for specific cases): flatironinstitute.github.io/solid_dmft . This is a triqs application is working similar to large DFT codes, with input file where you set all relevant parameters. No coding required.

Best, Alex

[quantumable]

Dear Alex,

best wishes for the year 2022!!

I have now changed the syntax so you can now check the script.

Indeed I tried to install soliDMFT it seems to me that this is not compatible with the version of triqs 3.0.0 present because I made a make and I had this error:

WARNING: autodoc: failed to import module 'csc_flow' from module 'solid_dmft'; the following exception was raised:
cannot import name 'gf_struct_flatten' from 'triqs_dft_tools.block_structure' (triqs/install/lib/python3.8/site-packages/triqs_dft_tools/block_structure.py)
WARNING: autodoc: failed to import module 'dmft_cycle' from module 'solid_dmft'; the following exception was raised:
**cannot import name 'gf_struct_flatten'** from 'triqs_dft_tools.block_structure' (/triqs/install/lib/python3.8/site-packages/triqs_dft_tools/block_structure.py)
WARNING: autodoc: failed to import module 'solver' from module 'solid_dmft.dmft_tools'; the following exception was raised:
cannot import name 'gf_struct_flatten' from 'triqs_dft_tools.block_structure' (//triqs/install/lib/python3.8/site-packages/triqs_dft_tools/block_structure.py)
WARNING: autodoc: failed to import module 'postproc_plotbands' from module 'solid_dmft.postprocessing'; the following exception was raised:
No module named 'TB_functions'

[the-hampel]

Hi, I believe you are trying to install solid_dmft against triqs/3.0.x? Then you have to use also this version of solid_dmft. You are trying to install a newer version. To do so check out version 3.0.x of solid_dmft before building, i.e.

git clone git@github.com:flatironinstitute/solid_dmft.git solid_dmft.src
git checkout 3.0.x

then follow the normal installation steps.

Yes now the script is readable and indeed the problem is the one I mentioned above. please try to incorporate my suggestion and make the use of icrsh coherent.

best, Alex

[quantumable]

Dear Alex,

I have made the following changes as suggested however I have some doubts about line 81 and line 110(S[icrsh].Sigma_iw << SK.dc_imp[icrsh]['up'][0,0]) I don't know if it should be S[icrsh].Sigma_iw << SK.dc_imp[icrsh]['up'][icrsh,icrsh] and line 143 (mpi.report("Orbital %s density: %.6f"%(sig,dm[sig][0,0].real))

can you please check it?

from itertools import *
import numpy as np
import triqs.utility.mpi as mpi
from h5 import *
# from triqs.archive import *
from triqs.gf import *
import sys, triqs.version as triqs_version
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 triqs_cthyb.version as cthyb_version
import triqs_dft_tools.version as dft_tools_version

import warnings
warnings.filterwarnings("ignore", category=FutureWarning)

filename = 'Mn'

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

beta = 101.79
p = {}
p["max_time"] = -1
p["random_name"] = ""
p["random_seed"] = 123 * mpi.rank + 567
p["length_cycle"] = 100
p["n_warmup_cycles"] = 10000
p["n_cycles"] = 100000
p["fit_max_moment"] = 4
p["fit_min_n"] = 30
p['imag_threshold']=50e-10
p["fit_max_n"] = 50
p["perform_tail_fit"] = True
DC_type = 0
DC_value = 20.51
n_iterations = 15
U = 4.5
J = 0.7

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-4)
h_int = [None] * SK.n_inequiv_shells
S = []
for icrsh in range(SK.n_inequiv_shells):
    n_orb = SK.corr_shells[icrsh]['dim']
    spin_names = ['up','down']
    orb_names = [i for i in range(icrsh,n_orb)]
    gf_struct = [(block, indices) for block, indices in SK.gf_struct_solver[icrsh].items()]
    mpi.report('Sumk to Solver: %s'%SK.sumk_to_solver)
    mpi.report('GF struct sumk: %s'%SK.gf_struct_sumk)
    mpi.report('GF struct solver: %s'%SK.gf_struct_solver)
    S.append(Solver(beta=beta, gf_struct=gf_struct))
    U_sph = U_matrix(l=2, U_int=U, J_hund=J)
    U_cubic = transform_U_matrix(U_sph, spherical_to_cubic(l=2, convention=''))
    Umat, Upmat = reduce_4index_to_2index(U_cubic)
    h_int[icrsh] = h_int_density(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U=Umat, Uprime=Upmat)
    mpi.report('Greens function structure is: %s '%gf_struct)
    mpi.report('U Matrix set to:\n%s'%Umat)
    mpi.report('Up Matrix set to:\n%s'%Upmat)
iteration_offset = 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 not 'DMFT_input' in ar: ar.create_group('DMFT_input')
    if not 'Iterations' in ar['DMFT_input']: ar['DMFT_input'].create_group('Iterations')
    if not 'code_versions' in ar['DMFT_input']: ar['DMFT_input'].create_group('code_versio\
ns')
    ar['DMFT_input']['code_versions']["triqs_version"] = triqs_version.version
    ar['DMFT_input']['code_versions']["triqs_git"] = triqs_version.git_hash
    ar['DMFT_input']['code_versions']["cthyb_version"] = cthyb_version.version
    ar['DMFT_input']['code_versions']["cthyb_git"] = cthyb_version.triqs_cthyb_hash
    ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.version
    ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.triqs_dft_tools_hash
    if 'iteration_count' in ar['DMFT_results']:
        iteration_offset = ar['DMFT_results']['iteration_count']+1
        for icrsh in range(SK.n_inequiv_shells):
            S[icrsh].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
    ar['DMFT_input']["dmft_script_it"+str(iteration_offset)] = open(sys.argv[0]).read()

for icrsh in range(SK.n_inequiv_shells):
    iteration_offset = mpi.bcast(iteration_offset)
    S[icrsh].Sigma_iw = mpi.bcast(S[icrsh].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)
    # Calc the first G0
    SK.symm_deg_gf(S[icrsh].Sigma_iw,ish=icrsh)
    SK.put_Sigma(Sigma_imp = [S[icrsh].Sigma_iw])
    SK.calc_mu(precision=0.01)
    S[icrsh].G_iw << SK.extract_G_loc()[icrsh]
    SK.symm_deg_gf(S[icrsh].G_iw, ish=icrsh)
#Init the DC term and the self-energy if no previous iteration was found

#Init the DC term and the self-energy if no previous iteration was found

for icrsh in range(SK.n_inequiv_shells):
    if iteration_offset == 0:
        dm = S[icrsh].G_iw.density()
        SK.calc_dc(dm, U_interact=U, J_hund=J, orb=icrsh, use_dc_formula=DC_type,use_dc_value=DC_value)
        S[icrsh].Sigma_iw << SK.dc_imp[icrsh]['up'][0,0]

mpi.report('%s DMFT cycles requested. Starting with iteration %s.'%(n_iterations,iteration_offset))

mpi.report('%s DMFT cycles requested. Starting with iteration %s.'%(n_iterations,iteration_offset))

# The infamous DMFT self consistency cycle
for it in range(iteration_offset, iteration_offset + n_iterations):

    mpi.report('Doing iteration: %s'%it)
    for icrsh in range(SK.n_inequiv_shells):
    # Get G0
        S[icrsh].G0_iw << inverse(S[icrsh].Sigma_iw + inverse(S[icrsh].G_iw))
    # Solve the impurity problem
        S[icrsh].solve(h_int = h_int[icrsh], **p)
        if mpi.is_master_node():
            ar['DMFT_input']['Iterations']['solver_dict_it'+str(it)] = p
            ar['DMFT_results']['Iterations']['G_0_it'+str(it)] = S[icrsh].G0_iw
            ar['DMFT_results']['Iterations']['Gimp_it'+str(it)] = S[icrsh].G_iw
            ar['DMFT_results']['Iterations']['Gtau_it'+str(it)] = S[icrsh].G_tau
            ar['DMFT_results']['Iterations']['Sigma_uns_it'+str(it)] = S[icrsh].Sigma_iw
    # Calculate double counting
        dm = S[icrsh].G_iw.density()
        SK.calc_dc(dm, U_interact=U, J_hund=J, orb=icrsh, use_dc_formula=DC_type,use_dc_value=DC_value)
    # Get new G
        SK.symm_deg_gf(S[icrsh].Sigma_iw,ish=icrsh)
        SK.put_Sigma(Sigma_imp=[S[icrsh].Sigma_iw])
        SK.calc_mu(precision=0.01)
        S[icrsh].G_iw << SK.extract_G_loc()[icrsh]

    # print densities
        for sig,gf in S[icrsh].G_iw:
            mpi.report("Orbital %s density: %.6f"%(sig,dm[sig][0,0].real))
        mpi.report('Total charge of Gloc : %.6f'%S[icrsh].G_iw.total_density().real)

        if mpi.is_master_node():
            ar['DMFT_results']['iteration_count'] = it
            ar['DMFT_results']['Iterations']['Sigma_it'+str(it)] = S[icrsh].Sigma_iw
            ar['DMFT_results']['Iterations']['Gloc_it'+str(it)] = S[icrsh].G_iw
            ar['DMFT_results']['Iterations']['G0loc_it'+str(it)] = S[icrsh].G0_iw
            ar['DMFT_results']['Iterations']['dc_imp'+str(it)] = SK.dc_imp
            ar['DMFT_results']['Iterations']['dc_energ'+str(it)] = SK.dc_energ
            ar['DMFT_results']['Iterations']['chemical_potential'+str(it)] = SK.chemical_potential
if mpi.is_master_node(): del ar

[the-hampel]

Hi, that looks already quite good. I think there is some confusion at the bottom on how to set the new self energies into the lattice green function. I made a correction:

        if mpi.is_master_node():
            ar['DMFT_input']['Iterations']['solver_dict_it'+str(it)+'_'+str(icrsh)] = p
            ar['DMFT_results']['Iterations']['G_0_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G0_iw
            ar['DMFT_results']['Iterations']['Gimp_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G_iw
            ar['DMFT_results']['Iterations']['Gtau_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G_tau
            ar['DMFT_results']['Iterations']['Sigma_uns_it'+str(it)+'_'+str(icrsh)] = S[icrsh].Sigma_iw
    # Calculate double counting
        dm = S[icrsh].G_iw.density()
        SK.calc_dc(dm, U_interact=U, J_hund=J, orb=icrsh, use_dc_formula=DC_type,use_dc_value=DC_value)
    # Get new G 
        SK.symm_deg_gf(S[icrsh].Sigma_iw,ish=icrsh)

    # print densities
        for sig,gf in S[icrsh].G_iw:
            mpi.report("Orbital %s density: %.6f"%(sig,dm[sig][0,0].real))
        mpi.report('Total charge of Gloc : %.6f'%S[icrsh].G_iw.total_density().real)

        if mpi.is_master_node():
            #here an impurity counter is missing on the left side
            ar['DMFT_results']['iteration_count'] = it
            ar['DMFT_results']['Iterations']['Sigma_it'+str(it)+'_'+str(icrsh)] = S[icrsh].Sigma_iw
            ar['DMFT_results']['Iterations']['Gloc_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G_iw
            ar['DMFT_results']['Iterations']['G0loc_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G0_iw
            ar['DMFT_results']['Iterations']['dc_imp'+str(it)+'_'+str(icrsh)] = SK.dc_imp
            ar['DMFT_results']['Iterations']['dc_energ'+str(it)+'_'+str(icrsh)] = SK.dc_energ
            ar['DMFT_results']['Iterations']['chemical_potential'+str(it)+'_'+str(icrsh)] = SK.chemical_potential

    sigma_list = [S[icrsh].Sigma_iw for icrsh in range(SK.n_inequiv_shells)]
    SK.put_Sigma(Sigma_imp=sigma_list)
    SK.calc_mu(precision=0.01)
    gloc_all = SK.extract_G_loc()
    for icrsh in range(SK.n_inequiv_shells):
        S[icrsh].G_iw << gloc_all[icrsh]

if mpi.is_master_node(): del ar

which replaces all lines below:

S[icrsh].solve(h_int = h_int[icrsh], **p)

please check carefully what I changed and let me know if this makes sense. No guarantee, that this works already but you have to implement some changes like this. put_sigma has to be called with a list of self-energies for all icrsh at once.

Best, Alex

[quantumable]

dear Alex , I have made the changes as you suggested the calculation seems to finish well here is the outpout file and .h5. making a h5ls there is. that the _0 files I thought we should have also the files of type _1 for the respective atoms 1 and 2. ? outpout.zip Mn.h5

Following your last remark, did you omit some points? I wonder if I should change map_operator_structure=SK.sumk_to_solver[0] to map_operator_structure=SK.sumk_to_solver[icrsh] S[icrsh].Sigma_iw << SK.dc_imp[icrsh]['up'][0,0] to S[icrsh].Sigma_iw << SK.dc_imp[icrsh]['up'][icrsh,icrsh]

dmft script.

from itertools import *
import numpy as np
import triqs.utility.mpi as mpi
from h5 import *
# from triqs.archive import *
from triqs.gf import *
import sys, triqs.version as triqs_version
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 triqs_cthyb.version as cthyb_version
import triqs_dft_tools.version as dft_tools_version

import warnings
warnings.filterwarnings("ignore", category=FutureWarning)

filename = 'Mn'

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

beta = 223.16
p = {}
p["max_time"] = -1
p["random_name"] = ""
p["random_seed"] = 123 * mpi.rank + 567
p["length_cycle"] = 100
p["n_warmup_cycles"] = 20000
p["n_cycles"] = 200000
p["fit_max_moment"] = 4
p["fit_min_n"] = 30
p['imag_threshold']=5e-6
p["fit_max_n"] = 60
p["perform_tail_fit"] = True
DC_type = 0
DC_value = 42.92
n_iterations = 20
U = 4.5
J = 0.7

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-4)
h_int = [None] * SK.n_inequiv_shells
S = []
for icrsh in range(SK.n_inequiv_shells):
    n_orb = SK.corr_shells[icrsh]['dim']
    spin_names = ['up','down']
    orb_names = [i for i in range(icrsh,n_orb)]
    gf_struct = [(block, indices) for block, indices in SK.gf_struct_solver[icrsh].items()]
    mpi.report('Sumk to Solver: %s'%SK.sumk_to_solver)
    mpi.report('GF struct sumk: %s'%SK.gf_struct_sumk)
    mpi.report('GF struct solver: %s'%SK.gf_struct_solver)
    S.append(Solver(beta=beta, gf_struct=gf_struct))
    U_sph = U_matrix(l=2, U_int=U, J_hund=J)
    U_cubic = transform_U_matrix(U_sph, spherical_to_cubic(l=2, convention=''))
    Umat, Upmat = reduce_4index_to_2index(U_cubic)
    h_int[icrsh] = h_int_density(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U=Umat, Uprime=Upmat)
    mpi.report('Greens function structure is: %s '%gf_struct)
    mpi.report('U Matrix set to:\n%s'%Umat)
    mpi.report('Up Matrix set to:\n%s'%Upmat)
iteration_offset = 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 not 'DMFT_input' in ar: ar.create_group('DMFT_input')
    if not 'Iterations' in ar['DMFT_input']: ar['DMFT_input'].create_group('Iterations')
    if not 'code_versions' in ar['DMFT_input']: ar['DMFT_input'].create_group('code_versio\
ns')
    ar['DMFT_input']['code_versions']["triqs_version"] = triqs_version.version
    ar['DMFT_input']['code_versions']["triqs_git"] = triqs_version.git_hash
    ar['DMFT_input']['code_versions']["cthyb_version"] = cthyb_version.version
    ar['DMFT_input']['code_versions']["cthyb_git"] = cthyb_version.triqs_cthyb_hash
    ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.version
    ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.triqs_dft_tools_hash
    if 'iteration_count' in ar['DMFT_results']:
        iteration_offset = ar['DMFT_results']['iteration_count']+1
        for icrsh in range(SK.n_inequiv_shells):
            S[icrsh].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
    ar['DMFT_input']["dmft_script_it"+str(iteration_offset)] = open(sys.argv[0]).read()

for icrsh in range(SK.n_inequiv_shells):
    iteration_offset = mpi.bcast(iteration_offset)
    S[icrsh].Sigma_iw = mpi.bcast(S[icrsh].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)
    # Calc the first G0
    SK.symm_deg_gf(S[icrsh].Sigma_iw,ish=icrsh)
    SK.put_Sigma(Sigma_imp = [S[icrsh].Sigma_iw])
    SK.calc_mu(precision=0.01)
    S[icrsh].G_iw << SK.extract_G_loc()[icrsh]
    SK.symm_deg_gf(S[icrsh].G_iw, ish=icrsh)
#Init the DC term and the self-energy if no previous iteration was found

#Init the DC term and the self-energy if no previous iteration was found

for icrsh in range(SK.n_inequiv_shells):
    if iteration_offset == 0:
        dm = S[icrsh].G_iw.density()
        SK.calc_dc(dm, U_interact=U, J_hund=J, orb=icrsh, use_dc_formula=DC_type,use_dc_value=DC_value)
        S[icrsh].Sigma_iw << SK.dc_imp[icrsh]['up'][0,0]

mpi.report('%s DMFT cycles requested. Starting with iteration %s.'%(n_iterations,iteration_offset))

mpi.report('%s DMFT cycles requested. Starting with iteration %s.'%(n_iterations,iteration_offset))

# The infamous DMFT self consistency cycle
for it in range(iteration_offset, iteration_offset + n_iterations):

    mpi.report('Doing iteration: %s'%it)
    for icrsh in range(SK.n_inequiv_shells):
    # Get G0
        S[icrsh].G0_iw << inverse(S[icrsh].Sigma_iw + inverse(S[icrsh].G_iw))
    # Solve the impurity problem
        S[icrsh].solve(h_int = h_int[icrsh], **p)

        if mpi.is_master_node():
            ar['DMFT_input']['Iterations']['solver_dict_it'+str(it)+'_'+str(icrsh)] = p
            ar['DMFT_results']['Iterations']['G_0_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G0_iw
            ar['DMFT_results']['Iterations']['Gimp_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G_iw
            ar['DMFT_results']['Iterations']['Gtau_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G_tau
            ar['DMFT_results']['Iterations']['Sigma_uns_it'+str(it)+'_'+str(icrsh)] = S[icrsh].Sigma_iw
    # Calculate double counting
        dm = S[icrsh].G_iw.density()
        SK.calc_dc(dm, U_interact=U, J_hund=J, orb=icrsh, use_dc_formula=DC_type,use_dc_value=DC_value)
    # Get new G 
        SK.symm_deg_gf(S[icrsh].Sigma_iw,ish=icrsh)

    # print densities
        for sig,gf in S[icrsh].G_iw:
            mpi.report("Orbital %s density: %.6f"%(sig,dm[sig][0,0].real))
        mpi.report('Total charge of Gloc : %.6f'%S[icrsh].G_iw.total_density().real)

        if mpi.is_master_node():
            #here an impurity counter is missing on the left side
            ar['DMFT_results']['iteration_count'] = it
            ar['DMFT_results']['Iterations']['Sigma_it'+str(it)+'_'+str(icrsh)] = S[icrsh].Sigma_iw
            ar['DMFT_results']['Iterations']['Gloc_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G_iw
            ar['DMFT_results']['Iterations']['G0loc_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G0_iw
            ar['DMFT_results']['Iterations']['dc_imp'+str(it)+'_'+str(icrsh)] = SK.dc_imp
            ar['DMFT_results']['Iterations']['dc_energ'+str(it)+'_'+str(icrsh)] = SK.dc_energ
            ar['DMFT_results']['Iterations']['chemical_potential'+str(it)+'_'+str(icrsh)] = SK.chemical_potential

    sigma_list = [S[icrsh].Sigma_iw for icrsh in range(SK.n_inequiv_shells)]
    SK.put_Sigma(Sigma_imp=sigma_list)
    SK.calc_mu(precision=0.01)
    gloc_all = SK.extract_G_loc()
    for icrsh in range(SK.n_inequiv_shells):
        S[icrsh].G_iw << gloc_all[icrsh]

if mpi.is_master_node(): del ar

local_lattice_GF.py

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 = '2Mn'
SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False)

beta = 222.16

# 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-4)

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

# 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)+'_'+str(icrsh)]
        SK.dc_imp = ar['DMFT_results']['Iterations']['dc_imp'+str(iteration_offset-1)+'_'+str(icrsh)]
        SK.dc_energ = ar['DMFT_results']['Iterations']['dc_energ'+str(iteration_offset-1)+'_'+str(icrsh)]
        SK.chemical_potential = ar['DMFT_results']['Iterations']['chemical_potential'+str(iteration_offset-1)+'_'+str(icrsh)].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:
for icrsh in range(SK.n_inequiv_shells):
    n_orbs = SK.corr_shells[icrsh]['dim']
    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, icrsh, bname, G_latt_KS[bname], gf, shells='corr', ir=None)
        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)+'_'+str(icrsh)] = G_latt_orb
if mpi.is_master_node(): del ar

[the-hampel]

Hi, I had a look at your output files. This h5 file contains in dft_input only one in-equivalent correlated shell, hence SK.n_inequiv_shells is only 1 and the script does exactly what it is supposed to do. There are two correlated shells in the h5 but they are marked as equivalent. I cannot tell you why, because this is how you run the converter. Please check the output of plovasp for that.
If you want to run for multiple impurities please check after running the converter that you have two in equivalent impurities. For example right in the beginning of your output because for example the line:

mpi.report('Sumk to Solver: %s'%SK.sumk_to_solver)

is only printed once, not twice. But the script is correct. Your correction to the script regarding map operator structure is indeed correct.

Best, Alex

[quantumable]

please this is the files

plo.cfg file

[General]
BASENAME = 2Mn
HK = False
COMPLEMENT = False
[Group 1]
SHELLS = 1 2
EWINDOW = -5.4 4
NORMION = False
NORMALIZE = True
[Shell 1]
IONS = 1
LSHELL = 2
CORR = True
TRANSFORM = 1.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  0.0
            0.0  0.0  0.0  1.0  0.0
            0.0  0.0  0.0  0.0  1.0
[Shell 2]
LSHELL = 2
IONS = 2
CORR = True
TRANSFORM = 1.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  0.0
            0.0  0.0  0.0  1.0  0.0
            0.0  0.0  0.0  0.0  1.0

converter

from triqs_dft_tools.converters.vasp import *
Converter = VaspConverter(filename = '2MnPc')
Converter.convert_dft_input()

INCAR

SYSTEM = MnPc
NPAR = 12
SIGMA = 0.02
EMIN = -10
EMAX = 3
LORBIT=14
LOCPROJ = 1 : d : Pr 1
LOCPROJ = 2 : d : Pr 1
LMAXMIX=4
ISYM=-1
ENCUT   = 450
EDIFF   = 1.0E-9
PREC = Accurate
ISMEAR = 0
NELM = 150

[the-hampel]

This looks okay to me. What does the converter tell you when you run it, please include such outputs? I cannot identify any problems from inputs alone. And also include a clear statement what you think the problem is. In general the converter will write to the output how many correlated shells there are, how they are sorted and if they are treated as equivalent. You see something like:

  Shell         : 1
  Orbital l     : 2
  Number of ions: 4
  Dimension     : 3
  Correlated    : True
  Ion sort      : [1, 1, 1, 1]

You should be able to identify the problem from the output of the converter. If this still gives you two equivalent shells we can have a closer look.

[quantumable]

Dear Alex,

thank you for your answer. yes indeed this is the input file with which I got the results I sent two days ago.

by doing plovasp plo.cfg a get this:

  Generating 2 shells...

    Shell         : 1
    Orbital l     : 2
    Number of ions: 1
    Dimension     : 5
    Correlated    : True
    Ion sort      : [1]

    Shell         : 2
    Orbital l     : 2
    Number of ions: 1
    Dimension     : 5
    Correlated    : True
    Ion sort      : [1]
Density matrix:
  Shell 1
    Site 1
     0.5035866     0.0003440     0.0000289    -0.0006677    -0.0004275
     0.0003440     1.2096755    -0.0000488    -0.0000021    -0.0001146
     0.0000289    -0.0000488     1.0608837     0.0000377     0.0005269
    -0.0006677    -0.0000021     0.0000377     1.2126918    -0.0000984
    -0.0004275    -0.0001146     0.0005269    -0.0000984     1.3818024
      trace:  5.3686399683957235
  Shell 2
    Site 1
     1.3818681    -0.0002922     0.0001148     0.0001472     0.0015300
    -0.0002922     1.2126044    -0.0002974    -0.0000715     0.0000435
     0.0001148    -0.0002974     1.0609893     0.0002330    -0.0001436
     0.0001472    -0.0000715     0.0002330     1.2102155    -0.0002432
     0.0015300     0.0000435    -0.0001436    -0.0002432     0.5035406
      trace:  5.369217830947333

[the-hampel]

But this should not mark the two sites as equivalent. This looks correct. If you open the h5 archive that comes out, what is the value of dft_input/n_inequiv_shells ? This should be 2. Otherwise something is not working correctly.

As a side note, your impurities seem to be very similar. Only that the dxy and dx2y2 orbitals are switched. You could make them maybe equivalent.

[quantumable]

Yes indeed it is equivalent this represents the system of the two atoms in vacuum. but if I put the system on a substrate I don't have the same thing so I need to treat them differently.

in substrat i have:

  Generating 2 shells...

    Shell         : 1
    Orbital l     : 2
    Number of ions: 1
    Dimension     : 5
    Correlated    : True
    Ion sort      : [1]

    Shell         : 2
    Orbital l     : 2
    Number of ions: 1
    Dimension     : 5
    Correlated    : True
    Ion sort      : [1]
Density matrix:
  Shell 1
    Site 1
     0.0757358     0.0006750    -0.0002519    -0.0011429    -0.0241535
     0.0006750     1.2409396     0.0017361     0.0000721     0.0009295
    -0.0002519     0.0017361     0.9714247     0.0010677     0.0086713
    -0.0011429     0.0000721     0.0010677     1.2581879    -0.0004599
    -0.0241535     0.0009295     0.0086713    -0.0004599     1.3450732
      trace:  4.891361148064956
  Shell 2
    Site 1
     1.5794863     0.0005720     0.0006804     0.0001488     0.0164668
     0.0005720     1.1283297    -0.0113357     0.0008812     0.0076627
     0.0006804    -0.0113357     1.0264248    -0.0005380    -0.0047712
     0.0001488     0.0008812    -0.0005380     1.1707086    -0.0000503
     0.0164668     0.0076627    -0.0047712    -0.0000503     0.0094296
      trace:  4.9143789646406

[the-hampel]

Before I investigate further can you please check: If you open the h5 archive that comes out from the converter, what is the value of dft_input/n_inequiv_shells ?

[quantumable]

I send you the .h5 file. please how can I see this value?

[the-hampel]

you can do:

h5ls -d 2Mn.h5/dft_input/n_inequiv_shells

[quantumable]

1-) this what I get:

h5ls -d 2Mn.h5/dft_input/n_inequiv_shells
n_inequiv_shells         Dataset {SCALAR}
    Data:
        (0) 1

for python converter.py I have this:

Starting run with 1 MPI rank(s) at : 2022-01-28 08:53:19.091456
Reading input from 2MnPc.ctrl...
{
    "ngroups": 1,
    "nk": 1,
    "ns": 1,
    "kvec1": [
        0.043263997869955176,
        0.0,
        0.0
    ],
    "kvec2": [
        -0.02497848081643802,
        0.04995696163287604,
        0.0
    ],
    "kvec3": [
        0.0,
        0.0,
        0.02702702702702703
    ],
    "nc_flag": 0
}

  No. of inequivalent shells: 1

2-) I don't know if this has to do with the SORT parameter: I didn't understand the use of this parameter. I tried to use SORT = 1 for the [Shell 1] and SORT = 2f or [Shell 2] I get this:

 h5ls  -d  2MnPc.h5/dft_input/n_inequiv_shells
n_inequiv_shells         Dataset {SCALAR}
    Data:
        (0) 2

but by running this: mpirun -np 24 python dmft.py

this is the error:

Starting run with 24 MPI rank(s) at : 2022-01-29 11:00:28.598788
Traceback (most recent call last):
  File "dmft.py", line 43, in <module>
    SK.put_Sigma([Sigma])
  File "/install/lib/python3.8/site-packages/triqs_dft_tools/sumk_dft.py", line 620, in put_Sigma
    Sigma_imp = self.transform_to_sumk_blocks(Sigma_imp)
  File "/install/lib/python3.8/site-packages/triqs_dft_tools/sumk_dft.py", line 675, in transform_to_sumk_blocks
    assert len(Sigma_imp) == self.n_inequiv_shells,\
AssertionError: transform_to_sumk_blocks: give exactly one Sigma for each inequivalent corr. shell!
Traceback (most recent call last):

the line 43 is

Sigma = SK.block_structure.create_gf(beta=beta)
SK.put_Sigma([Sigma])

[the-hampel]

Hi,

yes indeed that seems to be the problem. The way you enforced to become two inequivalent shells with manually specifying SORT=1 and SORT=2 gives the correct results. This must be then actually a bug in the vasp converter, as it should per default make the shells inequivalent. I will check this on Monday.

The error you then receive is of course reasonable. Because you did not specify in put_sigma the correct list of self energies. You have to give a list of self energies for each shell. So you have to change the line to give on self energy for each impurity problem solved.

[quantumable]

Dear Alex,

thank you for you response.

please how to specified a list of self energies?

[the-hampel]

I think this is pretty straightforward. You have the lines right there:

Sigma = SK.block_structure.create_gf(beta=beta)
SK.put_Sigma([Sigma, Sigma])

this would be a valid input for two impurities. Later in the script you have to of course make sure to always set the solutions of both impurities in the loop.

[quantumable]

Dear Alex,

I'm stuck on the script. I have the same type of error on this line:

S[icrsh].Sigma_iw = mpi.bcast(S[icrsh].Sigma_iw)
SK.put_Sigma(Sigma_imp = [S[icrsh].Sigma_iw]

error: AssertionError: transform_to_sumk_blocks: give exactly one Sigma for each inequivalent corr. shell!

I tried this:

Sigma_imp = ([S[icrsh].Sigma_iw,S[icrsh].Sigma_iw])
SK.put_Sigma(Sigma_imp)

[the-hampel]

I am having a hard time explaining this to you. I honestly don't know anymore how to explain. What you wrote will work, but it will not solve all your problems. This whole block is not correct for multiple impurities:

for icrsh in range(SK.n_inequiv_shells):
    iteration_offset = mpi.bcast(iteration_offset)
    S[icrsh].Sigma_iw = mpi.bcast(S[icrsh].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)
    # Calc the first G0
    SK.symm_deg_gf(S[icrsh].Sigma_iw,ish=icrsh)
    SK.put_Sigma(Sigma_imp = [S[icrsh].Sigma_iw])
    SK.calc_mu(precision=0.01)
    S[icrsh].G_iw << SK.extract_G_loc()[icrsh]
    SK.symm_deg_gf(S[icrsh].G_iw, ish=icrsh)

you need change this. The things you need to do are: 1) do first the mpi.bcast lines 2) Take the lines out of the loop that do not depend on icrsh 3) you can also here remove the lines that do symm_deg_gf completely for now to make life simpler 4) do `SK.put_Sigma outside any loop with the lines below 5) call extract_G_loc only once and store the result as a list

In the bottom of the script this is done correctly, please take a look:

sigma_list = [S[icrsh].Sigma_iw for icrsh in range(SK.n_inequiv_shells)]
SK.put_Sigma(Sigma_imp=sigma_list)
SK.calc_mu(precision=0.01)
gloc_all = SK.extract_G_loc()
for icrsh in range(SK.n_inequiv_shells):
    S[icrsh].G_iw << gloc_all[icrsh]

It seems that the concept is unclear to you, let me know which part is not understandable.

[quantumable]

dear Alex,

iteration_offset = mpi.bcast(iteration_offset)
SK.dc_imp = mpi.bcast(SK.dc_imp)
SK.chemical_potential = mpi.bcast(SK.chemical_potential)
SK.dc_energ = mpi.bcast(SK.dc_energ)
for icrsh in range(SK.n_inequiv_shells):
    S[icrsh].Sigma_iw = mpi.bcast(S[icrsh].Sigma_iw)
    SK.symm_deg_gf(S[icrsh].Sigma_iw,ish=icrsh)
    SK.put_Sigma(Sigma_imp = [S[icrsh].Sigma_iw])
    SK.calc_mu(precision=0.01)
    S[icrsh].G_iw << SK.extract_G_loc()[icrsh]
    SK.symm_deg_gf(S[icrsh].G_iw, ish=icrsh)

Yes, it's a problem because I don't understand your instructions. The problem is related to the fact that the line SK.put_Sigma(Sigma_imp = [S[icrsh].Sigma_iw]) does not take into account the two sigma. and following your instructions I can't do it.

Is there a concrete example that I can use as inspiration to finish?

iteration_offset = mpi.bcast(iteration_offset)
SK.dc_imp = mpi.bcast(SK.dc_imp)
SK.chemical_potential = mpi.bcast(SK.chemical_potential)
SK.dc_energ = mpi.bcast(SK.dc_energ)
for icrsh in range(SK.n_inequiv_shells):
    S[icrsh].Sigma_iw = mpi.bcast(S[icrsh].Sigma_iw) 
    SK.symm_deg_gf(S[icrsh].Sigma_iw,ish=icrsh)
    SK.put_Sigma(Sigma_imp = [S[icrsh].Sigma_iw])
    SK.calc_mu(precision=0.01)
    S[icrsh].G_iw << SK.extract_G_loc()[icrsh]
    SK.symm_deg_gf(S[icrsh].G_iw, ish=icrsh)
#Init the DC term and the self-energy if no previous iteration was found

#Init the DC term and the self-energy if no previous iteration was found

for icrsh in range(SK.n_inequiv_shells):
    if iteration_offset == 0:
        dm = S[icrsh].G_iw.density()
        SK.calc_dc(dm, U_interact=U, J_hund=J, orb=icrsh, use_dc_formula=DC_type,use_dc_value=DC_value)
        S[icrsh].Sigma_iw << SK.dc_imp[icrsh]['up'][0,0]

mpi.report('%s DMFT cycles requested. Starting with iteration %s.'%(n_iterations,iteration_offset))

mpi.report('%s DMFT cycles requested. Starting with iteration %s.'%(n_iterations,iteration_offset))

# The infamous DMFT self consistency cycle
for it in range(iteration_offset, iteration_offset + n_iterations):

    mpi.report('Doing iteration: %s'%it)
    for icrsh in range(SK.n_inequiv_shells):
    # Get G0
        S[icrsh].G0_iw << inverse(S[icrsh].Sigma_iw + inverse(S[icrsh].G_iw))
    # Solve the impurity problem
        S[icrsh].solve(h_int = h_int[icrsh], **p)

        if mpi.is_master_node():
            ar['DMFT_input']['Iterations']['solver_dict_it'+str(it)+'_'+str(icrsh)] = p
            ar['DMFT_results']['Iterations']['G_0_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G0_iw
            ar['DMFT_results']['Iterations']['Gimp_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G_iw
            ar['DMFT_results']['Iterations']['Gtau_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G_tau
            ar['DMFT_results']['Iterations']['Sigma_uns_it'+str(it)+'_'+str(icrsh)] = S[icrsh].Sigma_iw
    # Calculate double counting
        dm = S[icrsh].G_iw.density()
        SK.calc_dc(dm, U_interact=U, J_hund=J, orb=icrsh, use_dc_formula=DC_type,use_dc_value=DC_value)
    # Get new G 
        SK.symm_deg_gf(S[icrsh].Sigma_iw,ish=icrsh)

    # print densities
        for sig,gf in S[icrsh].G_iw:
            mpi.report("Orbital %s density: %.6f"%(sig,dm[sig][0,0].real))
        mpi.report('Total charge of Gloc : %.6f'%S[icrsh].G_iw.total_density().real)

        if mpi.is_master_node():
            #here an impurity counter is missing on the left side
            ar['DMFT_results']['iteration_count'] = it
            ar['DMFT_results']['Iterations']['Sigma_it'+str(it)+'_'+str(icrsh)] = S[icrsh].Sigma_iw
            ar['DMFT_results']['Iterations']['Gloc_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G_iw
            ar['DMFT_results']['Iterations']['G0loc_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G0_iw
            ar['DMFT_results']['Iterations']['dc_imp'+str(it)+'_'+str(icrsh)] = SK.dc_imp
            ar['DMFT_results']['Iterations']['dc_energ'+str(it)+'_'+str(icrsh)] = SK.dc_energ
            ar['DMFT_results']['Iterations']['chemical_potential'+str(it)+'_'+str(icrsh)] = SK.chemical_potential

    sigma_list = [S[icrsh].Sigma_iw for icrsh in range(SK.n_inequiv_shells)]
    SK.put_Sigma(Sigma_imp=sigma_list)
    SK.calc_mu(precision=0.01)
    gloc_all = SK.extract_G_loc()
    for icrsh in range(SK.n_inequiv_shells):
        S[icrsh].G_iw << gloc_all[icrsh]

if mpi.is_master_node(): del ar

[the-hampel]

I wrote you a specific example:

sigma_list = [S[icrsh].Sigma_iw for icrsh in range(SK.n_inequiv_shells)]
SK.put_Sigma(Sigma_imp=sigma_list)

This is exactly the line you need. You just cannot do that in the for loop. But exactly those two lines you can use to set the self energy as a list from the solver self energies.

Best, Alex

[quantumable]

Dear Alex,

I replace this SK.put_Sigma(Sigma_imp = [S[icrsh].Sigma_iw]) by this SK.put_Sigma(Sigma_imp = [S[icrsh].Sigma_iw for icrsh in range(SK.n_inequiv_shells)]) and I got errors later for this line:

S[icrsh].G0_iw << inverse(S[icrsh].Sigma_iw + inverse(S[icrsh].G_iw))

the error:

Traceback (most recent call last):
  File "d.py", line 121, in <module>
    S[icrsh].G0_iw << inverse(S[icrsh].Sigma_iw + inverse(S[icrsh].G_iw))
  File "nstall/lib/python3.8/site-packages/triqs/gf/tools.py", line 37, in inverse
    return x.inverse()
  File "/install/lib/python3.8/site-packages/triqs/gf/block_gf.py", line 417, in inverse
    return self.__class__( name_block_generator = [ (n, g.inverse()) for n,g in self], make_copies=False)
  File "/install/lib/python3.8/site-packages/triqs/gf/block_gf.py", line 417, in <listcomp>
    return self.__class__( name_block_generator = [ (n, g.inverse()) for n,g in self], make_copies=False)
  File "/install/lib/python3.8/site-packages/triqs/gf/gf.py", line 600, in inverse
    r.invert()
  File "/install/lib/python3.8/site-packages/triqs/gf/gf.py", line 586, in invert
    wrapped_aux._gf_invert_data_in_place(d)
RuntimeError: .. Error occurred at Tue Feb  8 17:40:32 2022

I was talking about an example, maybe you could help me to rectify the end of the script so that it is adapted to this type of calculation. I have no experience in this type of calculation.

[the-hampel]

I don't see how this connects to the changes. Can you attach or point to the d.py script? There must be a problem somewhere else in the script.

[quantumable]

from itertools import *
import numpy as np
import triqs.utility.mpi as mpi
from h5 import *
# from triqs.archive import *
from triqs.gf import *
import sys, triqs.version as triqs_version
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 triqs_cthyb.version as cthyb_version
import triqs_dft_tools.version as dft_tools_version

import warnings
warnings.filterwarnings("ignore", category=FutureWarning)

filename = '2Mn'

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

beta = 223.16
p = {}
p["max_time"] = -1
p["random_name"] = ""
p["random_seed"] = 123 * mpi.rank + 567
p["length_cycle"] = 100
p["n_warmup_cycles"] = 20000
p["n_cycles"] = 200000
p["fit_max_moment"] = 4
p["fit_min_n"] = 30
p['imag_threshold']=5e-6
p["fit_max_n"] = 60
p["perform_tail_fit"] = True
DC_type = 0
DC_value = 42.92
n_iterations = 1
U = 4.5
J = 0.7

Sigma = SK.block_structure.create_gf(beta=beta)
SK.put_Sigma([Sigma,Sigma])
G = SK.extract_G_loc()
SK.analyse_block_structure_from_gf(G, threshold = 1e-4)
h_int = [None] * SK.n_inequiv_shells
S = []
for icrsh in range(SK.n_inequiv_shells):
    l = SK.corr_shells[icrsh]['l']
    n_orb = SK.corr_shells[icrsh]['dim']
    spin_names = ['up','down']
    orb_names = [i for i in range(icrsh,n_orb)]
    gf_struct = [(block, indices) for block, indices in SK.gf_struct_solver[icrsh].items()]
    mpi.report('Sumk to Solver: %s'%SK.sumk_to_solver)
    mpi.report('GF struct sumk: %s'%SK.gf_struct_sumk)
    mpi.report('GF struct solver: %s'%SK.gf_struct_solver)
    S.append(Solver(beta=beta, gf_struct=gf_struct))
    U_sph = U_matrix(l=2, U_int=U, J_hund=J)
    U_cubic = transform_U_matrix(U_sph, spherical_to_cubic(l=2, convention=''))
    Umat, Upmat = reduce_4index_to_2index(U_cubic)
    h_int[icrsh] = h_int_density(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[icrsh], U=Umat, Uprime=Upmat)
    mpi.report('Greens function structure is: %s '%gf_struct)
    mpi.report('U Matrix set to:\n%s'%Umat)
    mpi.report('Up Matrix set to:\n%s'%Upmat)
iteration_offset = 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 not 'DMFT_input' in ar: ar.create_group('DMFT_input')
    if not 'Iterations' in ar['DMFT_input']: ar['DMFT_input'].create_group('Iterations')
    if not 'code_versions' in ar['DMFT_input']: ar['DMFT_input'].create_group('code_versio\
ns')
    ar['DMFT_input']['code_versions']["triqs_version"] = triqs_version.version
    ar['DMFT_input']['code_versions']["triqs_git"] = triqs_version.git_hash
    ar['DMFT_input']['code_versions']["cthyb_version"] = cthyb_version.version
    ar['DMFT_input']['code_versions']["cthyb_git"] = cthyb_version.triqs_cthyb_hash
    ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.version
    ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.triqs_dft_tools_hash
    if 'iteration_count' in ar['DMFT_results']:
        iteration_offset = ar['DMFT_results']['iteration_count']+1
        for icrsh in range(SK.n_inequiv_shells):
            S[icrsh].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
    ar['DMFT_input']["dmft_script_it"+str(iteration_offset)] = open(sys.argv[0]).read()
iteration_offset = mpi.bcast(iteration_offset)
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.calc_mu(precision=0.01)

for icrsh in range(SK.n_inequiv_shells):
    S[icrsh].Sigma_iw = mpi.bcast(S[icrsh].Sigma_iw)
    SK.symm_deg_gf(S[icrsh].Sigma_iw,ish=icrsh)
    SK.put_Sigma(Sigma_imp = [S[icrsh].Sigma_iw for icrsh in range(SK.n_inequiv_shells)])
    S[icrsh].G_iw << SK.extract_G_loc()[icrsh]
    SK.symm_deg_gf(S[icrsh].G_iw, ish=icrsh)
#Init the DC term and the self-energy if no previous iteration was found

#Init the DC term and the self-energy if no previous iteration was found

if iteration_offset == 0:
    for icrsh in range(SK.n_inequiv_shells):

        dm = S[icrsh].G_iw.density()
        SK.calc_dc(dm, U_interact=U, J_hund=J, orb=icrsh, use_dc_formula=DC_type,use_dc_value=DC_value)
        S[icrsh].Sigma_iw << SK.dc_imp[icrsh]['up'][0,0]

        mpi.report('%s DMFT cycles requested. Starting with iteration %s.'%(n_iterations,iteration_offset))

        mpi.report('%s DMFT cycles requested. Starting with iteration %s.'%(n_iterations,iteration_offset))

# The infamous DMFT self consistency cycle
        for it in range(iteration_offset, iteration_offset + n_iterations):

            mpi.report('Doing iteration: %s'%it)
            S[icrsh].G0_iw << inverse(S[icrsh].Sigma_iw + inverse(S[icrsh].G_iw))
    # Solve the impurity problem
            S[icrsh].solve(h_int = h_int[icrsh], **p)

    if mpi.is_master_node():
        ar['DMFT_input']['Iterations']['solver_dict_it'+str(it)+'_'+str(icrsh)] = p
        ar['DMFT_results']['Iterations']['G_0_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G0_iw
        ar['DMFT_results']['Iterations']['Gimp_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G_iw
        ar['DMFT_results']['Iterations']['Gtau_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G_tau
        ar['DMFT_results']['Iterations']['Sigma_uns_it'+str(it)+'_'+str(icrsh)] = S[icrsh].Sigma_iw
    # Calculate double counting
        dm = S[icrsh].G_iw.density()
        SK.calc_dc(dm, U_interact=U, J_hund=J, orb=icrsh, use_dc_formula=DC_type,use_dc_value=DC_value)
    # Get new G 
        SK.symm_deg_gf(S[icrsh].Sigma_iw,ish=icrsh)

    # print densities
        for sig,gf in S[icrsh].G_iw:
            mpi.report("Orbital %s density: %.6f"%(sig,dm[sig][0,0].real))
        mpi.report('Total charge of Gloc : %.6f'%S[icrsh].G_iw.total_density().real)

        if mpi.is_master_node():
            #here an impurity counter is missing on the left side
            ar['DMFT_results']['iteration_count'] = it
            ar['DMFT_results']['Iterations']['Sigma_it'+str(it)+'_'+str(icrsh)] = S[icrsh].Sigma_iw
            ar['DMFT_results']['Iterations']['Gloc_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G_iw
            ar['DMFT_results']['Iterations']['G0loc_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G0_iw
            ar['DMFT_results']['Iterations']['dc_imp'+str(it)+'_'+str(icrsh)] = SK.dc_imp
            ar['DMFT_results']['Iterations']['dc_energ'+str(it)+'_'+str(icrsh)] = SK.dc_energ
            ar['DMFT_results']['Iterations']['chemical_potential'+str(it)+'_'+str(icrsh)] = SK.chemical_potential

    sigma_list = [S[icrsh].Sigma_iw for icrsh in range(SK.n_inequiv_shells)]
    SK.put_Sigma(Sigma_imp=sigma_list)
    SK.calc_mu(precision=0.01)
    gloc_all = SK.extract_G_loc()
    for icrsh in range(SK.n_inequiv_shells):
        S[icrsh].G_iw << gloc_all[icrsh]

if mpi.is_master_node(): del ar

[the-hampel]

I think I am not sure how to help you. It is again very chaotic in the script. I think this is not what I asked for. You need to take one step at a time and make sure that you understand every line of the script. Let's start with this bit:

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.calc_mu(precision=0.01)

for icrsh in range(SK.n_inequiv_shells):
    S[icrsh].Sigma_iw = mpi.bcast(S[icrsh].Sigma_iw)
    SK.symm_deg_gf(S[icrsh].Sigma_iw,ish=icrsh)
    SK.put_Sigma(Sigma_imp = [S[icrsh].Sigma_iw for icrsh in range(SK.n_inequiv_shells)])
    S[icrsh].G_iw << SK.extract_G_loc()[icrsh]
    SK.symm_deg_gf(S[icrsh].G_iw, ish=icrsh)
#Init the DC term and the self-energy if no previous iteration was found

#Init the DC term and the self-energy if no previous iteration was found

this should look like this:

# first broadcast loaded variables
SK.dc_imp = mpi.bcast(SK.dc_imp)
SK.dc_energ = mpi.bcast(SK.dc_energ)
SK.chemical_potential = mpi.bcast(SK.chemical_potential)
for icrsh in range(SK.n_inequiv_shells):
    S[icrsh].Sigma_iw = mpi.bcast(S[icrsh].Sigma_iw)
    SK.symm_deg_gf(S[icrsh].Sigma_iw,ish=icrsh)

# now put sigma into the lattice Green function sumk
SK.put_Sigma(Sigma_imp = [S[icrsh].Sigma_iw for icrsh in range(SK.n_inequiv_shells)])

# recalc mu with the new self_energy
SK.calc_mu(precision=0.01)

# now extract the local Green function
Gloc_all = SK.extract_G_loc()

# loop over imp to set the local Gf into the solver
for icrsh in range(SK.n_inequiv_shells):
    S[icrsh].G_iw << Gloc_all[icrsh]
    SK.symm_deg_gf(S[icrsh].G_iw, ish=icrsh)

Can you explain to me what was wrong in your version? It really helps to make little comments so that I can better understand your code.

If this makes sense to you we can go on. The next part of your script afterwards is also not correct.

[quantumable]

Thank you for your message and correction. However the error message the same type of as before in this line:

S[icrsh].G0_iw << inverse(S[icrsh].Sigma_iw + inverse(S[icrsh].G_iw))

the error:

Traceback (most recent call last):
  File "d.py", line 121, in <module>
    S[icrsh].G0_iw << inverse(S[icrsh].Sigma_iw + inverse(S[icrsh].G_iw))
  File "nstall/lib/python3.8/site-packages/triqs/gf/tools.py", line 37, in inverse
    return x.inverse()
  File "/install/lib/python3.8/site-packages/triqs/gf/block_gf.py", line 417, in inverse
    return self.__class__( name_block_generator = [ (n, g.inverse()) for n,g in self], make_copies=False)
  File "/install/lib/python3.8/site-packages/triqs/gf/block_gf.py", line 417, in <listcomp>
    return self.__class__( name_block_generator = [ (n, g.inverse()) for n,g in self], make_copies=False)
  File "/install/lib/python3.8/site-packages/triqs/gf/gf.py", line 600, in inverse
    r.invert()
  File "/install/lib/python3.8/site-packages/triqs/gf/gf.py", line 586, in invert
    wrapped_aux._gf_invert_data_in_place(d)
RuntimeError: .. Error occurred at Tue Feb  

The second part of the script:

# The infamous DMFT self consistency cycle
        for it in range(iteration_offset, iteration_offset + n_iterations):

            mpi.report('Doing iteration: %s'%it)
            S[icrsh].G0_iw << inverse(S[icrsh].Sigma_iw + inverse(S[icrsh].G_iw))
    # Solve the impurity problem
            S[icrsh].solve(h_int = h_int[icrsh], **p)

    if mpi.is_master_node():
        ar['DMFT_input']['Iterations']['solver_dict_it'+str(it)+'_'+str(icrsh)] = p
        ar['DMFT_results']['Iterations']['G_0_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G0_iw
        ar['DMFT_results']['Iterations']['Gimp_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G_iw
        ar['DMFT_results']['Iterations']['Gtau_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G_tau
        ar['DMFT_results']['Iterations']['Sigma_uns_it'+str(it)+'_'+str(icrsh)] = S[icrsh].Sigma_iw
    # Calculate double counting
        dm = S[icrsh].G_iw.density()
        SK.calc_dc(dm, U_interact=U, J_hund=J, orb=icrsh, use_dc_formula=DC_type,use_dc_value=DC_value)
    # Get new G 
        SK.symm_deg_gf(S[icrsh].Sigma_iw,ish=icrsh)

    # print densities
        for sig,gf in S[icrsh].G_iw:
            mpi.report("Orbital %s density: %.6f"%(sig,dm[sig][0,0].real))
        mpi.report('Total charge of Gloc : %.6f'%S[icrsh].G_iw.total_density().real)

        if mpi.is_master_node():
            #here an impurity counter is missing on the left side
            ar['DMFT_results']['iteration_count'] = it
            ar['DMFT_results']['Iterations']['Sigma_it'+str(it)+'_'+str(icrsh)] = S[icrsh].Sigma_iw
            ar['DMFT_results']['Iterations']['Gloc_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G_iw
            ar['DMFT_results']['Iterations']['G0loc_it'+str(it)+'_'+str(icrsh)] = S[icrsh].G0_iw
            ar['DMFT_results']['Iterations']['dc_imp'+str(it)+'_'+str(icrsh)] = SK.dc_imp
            ar['DMFT_results']['Iterations']['dc_energ'+str(it)+'_'+str(icrsh)] = SK.dc_energ
            ar['DMFT_results']['Iterations']['chemical_potential'+str(it)+'_'+str(icrsh)] = SK.chemical_potential

    sigma_list = [S[icrsh].Sigma_iw for icrsh in range(SK.n_inequiv_shells)]
    SK.put_Sigma(Sigma_imp=sigma_list)
    SK.calc_mu(precision=0.01)
    gloc_all = SK.extract_G_loc()
    for icrsh in range(SK.n_inequiv_shells):
        S[icrsh].G_iw << gloc_all[icrsh]

if mpi.is_master_node(): del ar

[the-hampel]

Can you sent again the whole script? Maybe even with the input h5? Then I can try the script. From what you showed I can not see a problem.

You can also upload them to github in your own repo and I can always clone the most recent version.

Best, Alex

[quantumable]

Dear Alex,

thank you for your message. here is the complet script. 1)+'_'+str(icrsh)] SK.dc_imp = ar['DMFT_results']['Iterations']['dc_imp'+str(iterationoffset-1)+''+str(icrsh)] SK.dc_energ = ar['DMFT_results']['Iterations']['dc_energ'+str(iterationoffset-1)+''+str(icrsh)] SK.chemical_potential = ar['DMFT_results']['Iterations']['chemical_potential'+str(iterationoffset-1)+''+str(icrsh)].real

[quantumable]

Dear Alex,

please have a look at my files and script?

thank you

[the-hampel]

Hi, I think I told that many times before, please never share vasp executables or pseudopotentials. Please remove the tar archive above!

I checked your script and it seems more or less fine. I changed a few things to make it run better. But it is correct. The problem seems to be your input. The projectors are not correct for the second shell. Here is the updated script: dmft.zip

If you run it on a fresh h5 archive, i.e. after running the converter you can see the following output:

ccqlin027>cp orig_2Mn.h5 2Mn.h5 && mpirun -n 12 python dmft.py                            
cp: overwrite ‘2Mn.h5’? y
hwloc/linux: Ignoring PCI device with non-16bit domain.
Pass --enable-32bits-pci-domain to configure to support such devices                            
(warning: it would break the library ABI, don't enable unless really needed).                   
Starting run with 12 MPI rank(s) at : 2022-02-16 09:05:03.888718    
Dichotomy adjustment of Chemical Potential to obtain Total Density = 142.000000 +/- 0.000100    
Warning: Imaginary part in density will be ignored (6.285629850611754e-12)                      
    -0.054273 < Chemical Potential < -0.054273                                                  
    142.000034 < Total Density < 142.000034                                                     
    Chemical Potential found in 0 iterations :                                                  
    Total Density = 142.000034;Chemical Potential = -0.054273
non interacting density matrix:
up                                        
[[ 0.81095 -0.00044  0.00013  0.0002   0.00097]
 [-0.00044  0.7802  -0.00048  0.00002  0.00005]
 [ 0.00013 -0.00048  0.50475  0.00025 -0.00002]
 [ 0.0002   0.00002  0.00025  0.78122 -0.00017]
 [ 0.00097  0.00005 -0.00002 -0.00017  0.25178]]
--------------                            
down                                      
[[ 0.81095 -0.00044  0.00013  0.0002   0.00097]
 [-0.00044  0.7802  -0.00048  0.00002  0.00005]
 [ 0.00013 -0.00048  0.50475  0.00025 -0.00002]
 [ 0.0002   0.00002  0.00025  0.78122 -0.00017]
 [ 0.00097  0.00005 -0.00002 -0.00017  0.25178]]                                                
--------------                                                                                  
total occupation 6.2578
non interacting density matrix:
up
[[0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]]
--------------
down
[[0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]]
--------------
total occupation 0.0000

I modified the script so that it prints first thing the local occupations per impurity per spin channel. Which are somehow 0 for the second shell, whereas the first occupations look okayish! This will not work. I don't have time to debug this further. Please try to rerun Vasp and see what goes wrong. When running plovasp it will print the occupation matrix per site. Can you show me the output?

[quantumable]

Dear Alex, Thank you for your help. Indeed I have redone the correctly vasp calculation with high convergence: EDIFF =1.0E-10 and in "plovasp pro.cfg" I have indeed the two respective occupations of 5.322245776207522 and 5.322812115924464.

plovasp plo.cfg

Density matrix:
  Shell 1
    Site 1
     0.5028892     0.0003465     0.0000181    -0.0006829    -0.0004421
     0.0003465     1.1986041    -0.0000567    -0.0000107    -0.0001673
     0.0000181    -0.0000567     1.0225813     0.0000274     0.0001125
    -0.0006829    -0.0000107     0.0000274     1.1997720    -0.0000980
    -0.0004421    -0.0001673     0.0001125    -0.0000980     1.3983992
      trace:  5.322245776207522
  Shell 2
    Site 1
     1.3985505    -0.0001587     0.0001543     0.0000931     0.0015584
    -0.0001587     1.1998179    -0.0003013    -0.0000229     0.0000721
     0.0001543    -0.0003013     1.0227893     0.0002040    -0.0000495
     0.0000931    -0.0000229     0.0002040     1.1988454    -0.0002391
     0.0015584     0.0000721    -0.0000495    -0.0002391     0.5028090
      trace:  5.322812115924464

  Impurity density: 10.645057892131986

but the same message you got with 0. occupation for the second non interacting density matrix:

mpirun -np 24 python dmft.py

Dichotomy adjustment of Chemical Potential to obtain Total Density = 142.000000 +/- 0.000100
Warning: Imaginary part in density will be ignored (1.0454899506743113e-11)
    -0.054273 < Chemical Potential < -0.054273
    142.000034 < Total Density < 142.000034
    Chemical Potential found in 0 iterations : 
    Total Density = 142.000034;Chemical Potential = -0.054273
non interacting density matrix:
up
[[ 0.81095 -0.00044  0.00013  0.0002   0.00097]
 [-0.00044  0.7801  -0.00048  0.00002  0.00005]
 [ 0.00013 -0.00048  0.50466  0.00025 -0.00002]
 [ 0.0002   0.00002  0.00025  0.78113 -0.00017]
 [ 0.00097  0.00005 -0.00002 -0.00017  0.25139]]
--------------
down
[[ 0.81095 -0.00044  0.00013  0.0002   0.00097]
 [-0.00044  0.7801  -0.00048  0.00002  0.00005]
 [ 0.00013 -0.00048  0.50466  0.00025 -0.00002]
 [ 0.0002   0.00002  0.00025  0.78113 -0.00017]
 [ 0.00097  0.00005 -0.00002 -0.00017  0.25139]]
--------------
total occupation 6.2565
non interacting density matrix:
up
[[0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]]
--------------
down
[[0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]]
--------------
total occupation 0.0000
impurity #0

and the modifications made do not change the original error message. if we let the code run for a few minutes we get the same error message sur la ligne. S[icrsh].G0_iw << inverse(S[icrsh].Sigma_iw + inverse(S[icrsh].G_iw))

impurity #1
Traceback (most recent call last):
  File "dmft.py", line 154, in <module>
    S[icrsh].G0_iw << inverse(S[icrsh].Sigma_iw + inverse(S[icrsh].G_iw))
  File "/install/lib/python3.8/site-packages/triqs/gf/tools.py", line 37, in inverse
    return x.inverse()
  File "/install/lib/python3.8/site-packages/triqs/gf/block_gf.py", line 417, in inverse
    return self.__class__( name_block_generator = [ (n, g.inverse()) for n,g in self], make_copies=False)
  File "/install/lib/python3.8/site-packages/triqs/gf/block_gf.py", line 417, in <listcomp>
    return self.__class__( name_block_generator = [ (n, g.inverse()) for n,g in self], make_copies=False)
  File "/install/lib/python3.8/site-packages/triqs/gf/gf.py", line 600, in inverse
    r.invert()
  File "/install/lib/python3.8/site-packages/triqs/gf/gf.py", line 586, in invert
    wrapped_aux._gf_invert_data_in_place(d)
RuntimeError: .. Error occurred at Mon Feb 21 00:48:30 2022

[the-hampel]

Hm, I slowly start to wonder if something is wrong in the Vasp converter. Could you please try the following plo.cfg file instead:

[General]
BASENAME = 2Mn
EFERMI = -4.1790368
[Group 1]
SHELLS = 1
EWINDOW = -5.4 4
#NORMION = False
NORMALIZE = True
[Shell 1]
IONS = 1 2
LSHELL = 2

I cannot check it right now, but it should do exactly the same thing. So just replace everything in your plo.cfg with this file and try again.

[quantumable]

Dear Alex, thanks for your help. using the pro.cfg you suggested this looks good.

i will now test the self-consistent calculation.

Best

[the-hampel]

Okay. I will close this issue for now. Please open a new one if necessary. I will take a look why your plo.cfg failed.