FIX: Luttinger integral for odd wres (w=0 divergence)

The original definition, as given by Logan, Tucker and Galpin, prescribes
to integrate over the negative semiaxis:

    I = \frac{2}{\pi} \Im\int_{-\infty}^{0} dw G_loc(w) d\Sigma(w) / dw

Nevertheless we find most useful to exploit particle-hole symmetry and
thus integrate over the whole real-axis. Care has to be taken whenever
the frequency domain contains the origin, for therein resides a nasty
pole, arising from the derivative of the self-energy. Thus we divide
the integrand in (w < 0) and (w > 0) parts, and integrate the union.
Failing to exclude the (w = 0) point, if exist, leads to a diverging
uncontrolled result.

-------

Now the odd wres case should be totally safe.

ROADMAP.md

@@ -17,7 +17,7 @@
 - [x] Compute a mottness-marker based on the divergence of the scattering rate (Im[Sigma(0)]... very sensible, basically inexpensive), so to obtain sharper phase diagrams with respect to the Z-derived ones. It should also give insight about the supercritical phases: bad metal and pseudogap. [**this yet to check**]   
   > At the moment we implement the "strenght of correlations", S = norm[Sigma(0)-Sigma(∞)], as defined in `PRL 114 185701 (2015)`, with actual neat results: the marker is almost zero accross the whole FL phase and starts increasing very fast in the Mott insulator.
 
-- [x] Compute the Luttinger integral, as defined in `J. Phys.: Condens. Matter 28 (2016) 025601` and `PRB 102 081110(R) (2020)`. Since it appears to be quantized _at very low temperatures_, it could become the definitive **flag** for _quantum_ phase diagrams; much better than Z or S for it is an **integer**. [→ easier phase-boundary recognition!]  
+- [x] Compute the Luttinger integral, as defined in `PRB 90 075150 (2014)`. Since it appears to be quantized _at very low temperatures_, it could become the definitive **flag** for _quantum_ phase diagrams; much better than Z or S for it is an **integer**. [→ easier phase-boundary recognition!]  
     > Luttinger Theorem currently works for very low temperatures only. It might well be an inherent limitation, restricting its domain to the quantum Mott transition. Also note that to have a sharp first order step-up at the transition a very highly frequency resolution is needed, so to make the IPT solver performance-critical! [solved brilliantly with fast convolutions, see the `solver-optimization` section below]
 
 - [x] `SOLVER-OPTIMIZATION`: make the SOPT runs faster, by optimizing the needed convolutions. [implemented a pow2-optimized FFTW-based custom algorithm that actually greatly improves the cpu-time for the `wres=2^15` calculation: almost a x10 overall speedup!]

code/+phys/luttinger.m

@@ -1,31 +1,56 @@
 function [I,Lfig] = luttinger(w,sloc,gloc)
-%% Computes the Luttinger sum-rule, as defined in PRB 102 081110
+%% Computes the Luttinger sum-rule, as defined in PRB 90 075150
 %
 %  Input:
 %       w           : real valued array, \omega domain
 %       sloc        : complex valued array, \Sigma(\omega) function
 %       gloc        : complex valued array, G_{loc}(\omega) function
 %  Output:
-%       I = \frac{2}{\pi}\Im\int_{-\infty}^{0} dw G_loc(w) d\Sigma(w) / dw
-%         = \frac{1}{\pi}\Im\int_{-\infty}^{\infty}dwG_loc(w)d\Sigma(w)/dw
+%       I = \frac{1}{\pi}\Im\int_{-\infty}^{\infty}dwG_loc(w)d\Sigma(w)/dw
 %
 %% BSD 3-Clause License
 % 
 %  Copyright (c) 2022, Gabriele Bellomia
 %  All rights reserved.
                                                                global DEBUG
-% w = w(w<=0);          % Faster (half points to sum)
-% sloc = sloc(w<=0);    % but apparently less accurate
-% gloc = gloc(w<=0);    % --> better to exploit ph-sym
 
-ds  = diff(sloc); % dSigma -> already eliminates dw!
-g_  = gloc(1:end-1);
-integrand  = imag(g_.*ds);
+% The original definition, as given by Logan, Tucker and Galpin, prescribes
+% to integrate over the negative semiaxis:
+%
+%   I = \frac{2}{\pi} \Im\int_{-\infty}^{0} dw G_loc(w) d\Sigma(w) / dw
+%      
+% Nevertheless we find most useful to exploit particle-hole symmetry and
+% thus integrate over the whole real-axis. Care has to be taken whenever
+% the frequency domain contains the origin, for therein resides a nasty
+% pole, arising from the derivative of the self-energy. Thus we divide 
+% the integrand in (w < 0) and (w > 0) parts, and integrate the union.
+% Failing to exclude the (w = 0) point, if exist, leads to a diverging
+% uncontrolled result.
+
+wl = w(w<0);        % Build w < 0 patch
+gl = gloc(w<0);
+sl = sloc(w<0); 
+wl = wl(1:end-1);                   
+gl = gl(1:end-1);
+dl = diff(sl);
+
+wr = w(w>0);        % Build w > 0 patch
+gr = gloc(w>0); 
+sr = sloc(w>0);
+wr = wr(1:end-1); 
+gr = gr(1:end-1);       
+dr = diff(sr);    
+
+w  = [wl,wr];       % Join the patches
+ds = [dl,dr];
+g  = [gl,gr];
+
+integrand  = imag(g.*ds);
 I = 1/pi*sum(integrand);
-I = abs(I); % J. Phys.: Condens. Matter 28 (2016) 025601
+I = abs(I);
                                                                    if DEBUG
 Lfig = figure("Name",'Luttinger integrand','Visible','off');
-plot(w(1:end-1),integrand);
+plot(w,integrand);
 xlabel('\omega');
 ylabel('Im[G(\omega)\partial\Sigma/\partial\omega]');
 ylim([-0.1,0.2]);

code/main.m

@@ -1,6 +1,6 @@
 %% BSD 3-Clause License
 % 
-% Copyright (c) 2020, Gabriele Bellomia
+% Copyright (c) 2020-2022, Gabriele Bellomia
 % All rights reserved.
 
 clearvars; clc;
@@ -17,7 +17,7 @@
 
 %% INPUT: Boolean Flags
 MottBIAS     = 0;       % Changes initial guess of gloc (strongly favours Mott phase)
-ULINE        = 0;       % Takes and fixes the given beta value and performs a U-driven line
+ULINE        = 1;       % Takes and fixes the given beta value and performs a U-driven line
 TLINE        = 0;       % Takes and fixes the given U value and performs a T-driven line
 UTSCAN       = 0;       % Ignores both given U and beta values and builds a full phase diagram
 SPECTRAL     = 0;       % Controls plotting of spectral functions
@@ -35,7 +35,7 @@
 wres  = 2^15+1;           % Energy resolution in real-frequency axis
 wcut  = 6.00;           % Energy cutoff in real-frequency axis
 Umin  = 0.00;           % Hubbard U minimum value for phase diagrams
-Ustep = 0.09;           % Hubbard U incremental step for phase diagrams
+Ustep = 0.10;           % Hubbard U incremental step for phase diagrams
 Umax  = 6.00;           % Hubbard U maximum value for phase diagrams
 Tmin  = 1e-3;           % Temperature U minimum value for phase diagrams
 Tstep = 1e-3;           % Temperature incremental step for phase diagrams
@@ -94,7 +94,8 @@
         fprintf('< U = %f\n',U);
         [gloc{i},sloc{i}] = dmft_loop(gloc_0,w,D,U,beta,mloop,mix,err,'quiet');
         Z(i) = phys.zetaweight(w,sloc{i});
-        I(i) = phys.luttinger(w,sloc{i},gloc{i});
+        [I(i),infig] = phys.luttinger(w,sloc{i},gloc{i});
+        plot.push_frame('luttok.gif',i,mloop,dt,infig);
         S(i) = phys.strcorrel(w,sloc{i});
     end
     if(PLOT)