Issue with Gloc DOS version

[lcrippa]

Hi,

a bit of context. I was looking with @alexkowalski at the imaginary part of the delta function of a given model, and he realized the tails were slightly decreasing. This is of course wrong.

That specific model was a bethe lattice, so calculations were done with the DOS. I have now written a comparative driver for a square lattice which can calculate quantities both with sum over k or with the dos. What I notice is the following: clearly the self-energy is identical (see here the last Matsubara frequency of 10000 for beta=80)

But the Gloc is actually slightly different (again, last frequency)

And so is Gloc ^-1

So that finally, when Delta is calculated (the routine is the same in both cases) we have the tail issue

Since this issue is persistent across drivers, and also in the bethe case, I doubt it's a consequence of my routine calculating the DOS of the square lattice from H(k).

There are indeed some issues I have in the GLOC calculating routine:

For example, the code block

if(present(diagonal))then dos_diag=diagonal else dos_diag = .not.(size(Dbands,1) < size(Ebands,1)) endif

is missing in all the overloaded functions bar rank2, so I don't know if the dos is being considered diagonal or not. Secondly, I don't understand precisely this condition...

I'm looking into it but I would appreciate help on this :)

[lcrippa]

Update: actually, I might have been mistaken. I have tried on the Bethe lattice to compare the Delta function coming out from dmft_self_consistency with the one from the analytic formula. They look quite different. However, the Gloc obviously now is the same. Perhaps there is something in dmft_self_consistency? The issue on the "diagonal" flag in the previous message remains valid though

[SamueleGiuli]

Hi Lorenzo,

The line about the diagonal flag is to compute Gloc on multiorbital systems with different bandwidths. I used it for example for calculations on 2 orbitals Bethe lattice where the orbitals had different "hoppings" since the DOS is really the same but Eband is orbital dependent (e.g. Eband[1,:] = alpha* Eband[2,:] ) otherwise you can only treat SU(2N) band dispersion without considering different Ebands. Whenever you are passing arrays with "size(Dbands,1) < size(Ebands,1)" it means you are in this situation of different Ebands for different spin-orbitals. At the time DMFT_GLOC had only the rank2 routines if I remember well (and the GLOC_MATSUBARA.f90 + many other files structure) and I guess it was not extended to the other ranks when the shift to various ranks was made.

On Thu, Oct 31, 2024 at 12:13 PM Lorenzo Crippa @.***> wrote:

Update: actually, I might have been mistaken. I have tried on the Bethe lattice to compare the Delta function coming out from dmft_self_consistency with the one from the analytic formula. They look quite different. However, the Gloc obviously now is the same. Perhaps there is something in dmft_self_consistency? The issue on the "diagonal" flag in the previous message remains valid though

— Reply to this email directly, view it on GitHub https://github.com/aamaricci/DMFTtools/issues/6#issuecomment-2450276098, or unsubscribe https://github.com/notifications/unsubscribe-auth/ASVA45B77SWTT2A4AT3ZLJTZ6JJJZAVCNFSM6AAAAABQ6SDX76VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDINJQGI3TMMBZHA . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[aamaricci]

mmmmhhh....

1. I entirely rewrote the GK, GLOC and SELF_CONSISTENCY modules recently (left the legacy somewhere available though) to fix a huge potential memory problem.

2. Not sure to follow completely the issue without any math expression. dmft_self_consistency is, by default, returning the Weiss field (calG0 obtained from calG0 = [Gloc^-1 + Sigma]^-1 ) To get Delta once need to pass a specific variable. In that case \Delta is obtained essentially from \calG0. In the Bethe lattice this can be skipped as we know that \Delta = t^2 G where G can be either G_loc or G_imp depending on how you implement the code.

As Samuele noticed, I might have not extended the full support to all possible ranks for the DOS case. For sure there should be rank-2 [Nso,Nso] and rank-4 [Nspin,Nspin,Norb,Norb] (using them right now). The multi-orbital support for the DOS was extended by Samuele as he explained.

Let's fix some points: what are we comparing here?

On Thu, Oct 31, 2024 at 5:13 PM Lorenzo Crippa @.***> wrote:

Update: actually, I might have been mistaken. I have tried on the Bethe lattice to compare the Delta function coming out from dmft_self_consistency with the one from the analytic formula. They look quite different. However, the Gloc obviously now is the same. Perhaps there is something in dmft_self_consistency? The issue on the "diagonal" flag in the previous message remains valid though

— Reply to this email directly, view it on GitHub https://github.com/aamaricci/DMFTtools/issues/6#issuecomment-2450276098, or unsubscribe https://github.com/notifications/unsubscribe-auth/AARRDWCXF73E6SS4PFQ2HKTZ6JJJZAVCNFSM6AAAAABQ6SDX76VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDINJQGI3TMMBZHA . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[lcrippa]

@ Samuele,

hm, but why would the situation when Ebands are different (i.e., if I understand correctly, the dos have a different width per orbital) be off-diagonal?

@ Adriano,

I have written a minimal version of the driver I was using, to show the difference in the Delta. However, with this one, when I run it it fails with a segfault at the stage of calculating Gloc. I can't see what is the problem (probably something stupid) and unfortunately debug flags don't help. If you want to take a look, i attach it.

Once it's running, the only thing would be to compare bethesc=True with bethesc=False. Sorry for the zip file, github doesn't let you share txt in linux

test_delta.zip

[SamueleGiuli]

yeah the name may be misleading, it is the error of a younger me. I also think at a certain point a developed a non-diagonal version since I needed inter-orbital hoppings but then I never ended up pushing it to master.

On the delta problem: I would not expect the sigmas and Gloc to be the same for DoS and k-sum but just very similar, as they are.Although as you notice the delta should not behave as such but delta really depends only on the parameters of the impurity problem and the routine is the same right? So the parameters are different but still you should not see this effect since Delta should go as ~ 1 / ( i \omega ). Are you sure you are using the same routines for delta? One check could be to just taking the converged parameters for the DoS solution and make only one loop with the k-point calculation and viceversa and look if the problem is still there

On Thu, Oct 31, 2024 at 1:31 PM Lorenzo Crippa @.***> wrote:

@Samuele https://github.com/Samuele,

hm, but why would the situation when Ebands are different (i.e., if I understand correctly, the dos have a divverent width per orbital) be off-diagonal?

@adriano https://github.com/adriano,

I have written a minimal version of the driver I was using, to show the difference in the Delta. However, with this one, when I run it it fails with a segfault at the stage of calculating Gloc. I can't see what is the problem (probably something stupid) and unfortunately debug flags don't help. If you want to take a loo, i attach it.

Once it's running, the only thing would be to compare bethesc=True with bethesc=False. Sorry for the zip file, github doesn't let you share txt in linux

test_delta.zip https://github.com/user-attachments/files/17591195/test_delta.zip

— Reply to this email directly, view it on GitHub https://github.com/aamaricci/DMFTtools/issues/6#issuecomment-2450445226, or unsubscribe https://github.com/notifications/unsubscribe-auth/ASVA45G3ZEQZBB4SI52G6RDZ6JSQ7AVCNFSM6AAAAABQ6SDX76VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDINJQGQ2DKMRSGY . You are receiving this because you commented.Message ID: @.***>

[lcrippa]

Hi Samuele,

no problem for the naming, i was just curious about the origin. I'm pretty sure I'm using the same routine... i start from the same bath, do only one loop, stop after calculating delta. Yes, I also think the Gloc may not be the problem (see second message). If there indeed is a bug, it could be in the function that calculates delta from the dos. In fact, if I calculate the Weiss field, the does and the k-sum are much more reasonable

[lcrippa]

Ok, some more debugging:

test_delta.zip

find attached a zip folder. It contains a driver (test_delta.f90), a script (run.sh) and some files. Compile the driver and run the script, it will produce the delta and weiss files (.data extension) to be compared. Or, use the ones already there. Weiss are on top, ImDelta from the generic routine is decreasing for ω→∞.

Also, now the driver doesn't segfault anymore. It works with the rank-2 routine for obtaining gloc, but not with the rank-4. That one sometimes segfaults (in this driver it does, but not in another one...)

Thanks for the help Lorenzo

[SamueleGiuli]

Even more misterious... But maybe I spotted a critically dangerous typo: In get_glocal you wrote:

call dmft_get_gloc(Ebands,Dbands,H0,Smats,Smats,axis='mats',diagonal=.true.)

Therefore passing Smats twice! One of those should be a Gmats and maybe that has to be with the errors you are finding. Your Gmats at that point is really zero and Smats will definitely be changed inside dmft_get_gloc in two different ways for DOS e ksum. I can't compile and run now but I may find time later or tomorrow. I hope that helps!

On Thu, Oct 31, 2024 at 3:16 PM Lorenzo Crippa @.***> wrote:

Ok, some more debugging:

test_delta.zip https://github.com/user-attachments/files/17592242/test_delta.zip

find attached a zip folder. It contains a driver (test_delta.f90), a script (run.sh) and some files. Compile the driver and run the script, it will produce the delta and weiss files (.data extension) to be compared. Or, use the ones already there. Weiss are on top, delta from the generic routine is decreasing.

Also, now the driver doesn't segfault anymore. It works with the rank-2 routine for obtaining gloc, but not with the rank-4. That one sometimes segfaults (in this driver it does, but not in another one...)

Thanks for the help Lorenzo

— Reply to this email directly, view it on GitHub https://github.com/aamaricci/DMFTtools/issues/6#issuecomment-2450643733, or unsubscribe https://github.com/notifications/unsubscribe-auth/ASVA45G42F7XCMSKWI3YIYTZ6J6Z7AVCNFSM6AAAAABQ6SDX76VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDINJQGY2DGNZTGM . You are receiving this because you commented.Message ID: @.***>

[lcrippa]

Hehe no... that was a test I made to see that I allocated stuff with the correct dimension, maybe Gmats was wrong... sorry for the confusion, the first version of the drive I sent had that bug on purpose. The one in the archive (last attachment, the one with many files) is correct. You can spot the difference because there I convert to rank-2

:)

[aamaricci]

The seg fault happens in the rank-4 because clearly dos_diag is undefined in the get_gloc_dos_normal_rank4

I forgot to include the if;else statement which defines it but then this variable is used in the code. I will fix it.

Now I'll look into the delta issue.

[lcrippa]

Hm, I added the if;else yesterday in my local copy to see if that was the problem, but it kept segfaulting... did it work for you?

[aamaricci]

I see. I'll check and let you know...

[aamaricci]

Concerning the seg fault: still may bad... (compiled DMFT_TOOLS with BUILD_TYPE=debug, loaded as module, recompiled your test_delta code) At line 607 of file /Users/amaricci/projects_src/dmft_tools/src/DMFT_GF/GF_GLOC.f90 Fortran runtime error: Index '1' of dimension 1 of array 'wfreq' above upper bound of 0

Fix1: added at line 596:

 if(present(diagonal))then
       dos_diag=diagonal
    else
       dos_diag = .not.(size(Dbands,1) < size(Ebands,1))
    endif

Fix2: added at line 601: call build_frequency_array(axis)

Code works. Now I will check about Delta: if I get this correctly I might know the answer already.... ;-)

[lcrippa]

Confirm the code doesn't segfault now. I wonder which frequency array it used when it didn't segfault before (other drivers than test_delta). The weird delta behavior with the generic function is still present (as it should, rank-2 was not solving the issue). Thanks for looking into this

EDIT: it seems to persist independently on rank. Maybe something in the MPI splitting or in SC_COMMON?

[aamaricci]

Gotcha: the error is..... the integral with respect to the DOS to get Gloc.

As I suspected the is an incorrect cancellation in: \Delta = iw + \mu - H0 - \calG^{-1} = (as for self-consistency) = iw + \mu - H0 - \Sigma - Gloc^{-1}

Now, for \Delta to fall off as 1/iw it is necessary that Gloc^{-1} cancels exactly with iw+\mu-H0, leaving only \Sigma to rule the tail behavior. As you may know this kind of tail cancellations are quite tricky. In general they work in exact Dyson equations (as for instance that for the impurity OR the one for G_k) but not as much for the self-consistency (which looks like but it is not a Dyson equation). Because Gloc is not obtained to the required level of precision, this tail cancellation at some point does not hold anymore. [people doing QMC know this very well indeed]

Problem: \Delta_bethe=t*2Gloc (note this is already an incorrect way of writing the self-consistency as Gimp should be used) is different from \Delta above obtained using a generic self-consistency relation not restricted to the Bethe case.

I evaluate the error Y = abs|\Im\Delta_bethe - \Im\Delta|. This number is a linear function of iw Y = aw_n. For the parameter Lorenzo shared a=31e-6 \Re\Delta should be zero by p-h symmetry.

Check Hp1: the problem does not depend on the implementation of DMFT_TOOLS. I implemented in the code test_delta.f90 an equivalent but simplified version of the DMFT_TOOLS one. The problem persists.

Check Hp2: the problem does not depend on the inversion of Gloc appearing above. I implemented a version of the procedure returning Delta using 1/Gloc(a,a,i) rather than using the Lapack inversion procedure. The problem persists identical.

Check Hp3: the problem is not the number of Matsubara frequencies used in the tail. I tested using Lmats=84096. Nothing change and the error persists and also a=31e-6 identical

The is the integral for Gloc. Check Hp4: the integral for Gloc using Bethe DOS is misrepresented. I tested decreasing of 1/2 and increasing x10 the number of energy points sampling the Bethe DOS, finding the same linear behavior of Y but with: a = 11e-6 (Le=10000), a = 31e-6 (Le=1000), a= 9.5*1e-5 (Le=500)

Final test, solve the error: It is in principle impossible to get arbitrary precision in the integral for Gloc but we can use the exact Hilbert transform instead (in SciFor this is implemented as gfbethe and gfbether for, resp., matsubara and realaxis).

Result: the error is vanishingly small, reduced to noise of order 1d-13 around zero.

Cheers, A

test_delta_solved.tgz test_delta.pdf

[lcrippa]

Alright, that's a very thorough explanation, thanks a lot! (also sorry I got the mail only now). At least the thread wasn't completely useless since we got the rank-4 issue

Closing :)