FIX: bandwidth scaling (cfr. 81bd520)

HISTORY

Working on the lattices feature branch we found a severe problem in the code, namely that changing both U and D maintaining U/D fixed, DOES change the results (which of course makes no sense).

- Initially I assumed to have introduced the bug myself, while working on the generic-dos extension. But the problems was still present if the Bethe lattice was selected.

- So I turned to the +plot module (cfr. a32d24a), haunting for some silly error. No, the results *really* changed varying U and D together.

- Then I feared to have forgotten something important in the SOPT solver (the only place where U really enters and suspiciously D is totally absent). But no, it was fine.†

- At this point it was clear that the problem has been living unnoticed since the far past: in fact, it can be reproduced in the original notebook by Najera... BOOM!

So here I fix the notebook itself. 

> Maybe I should avoid doing this... isn't it meant to be a *static* resource? But what about future confusion if we diverge in such a core functionality. Also I may want the notebook to be useful to cloners, so yes, let's change it and make it clear somewhere in the markdown files: I could even "complete" it, by coding the exercises. :)

---------------------------

THE PROBLEM AND THE FIX

Najera had D=1 hardcoded in the <semi_circle_hiltrans()> function, so everything worked for itself. But if you make D actually a variable everything breaks, as in the Matlab code, for it is hardcoded also in another, somewhat hidden, place: a goddamn magic number!

    g0 = 1 / (w + eta - .25 * gloc)

which, for D=1, would stand for

    g0 = 1 / (w + eta - t**2 * gloc)

since for the Bethe lattice we have t = D/2

So we can fix the D ≠ 1 case by just changing it to

    g0 = 1 / (w + eta - .25 * D**2 * gloc)

- The same thing has been done in the dmft_loop() Matlab function. 
- The solution has been extensively tested in both Matlab's and Najera's code and it solves indeed the problem.
- A more general, lattice agnostic, solution has been implemented in the 'lattices' feature branch, where it is relevant. It will be merged here after the 'Bethe-only' version of the code is released, together with all the 'generic-lattices' material.
  > Cfr. 81bd520

---------------------------

FURTHER NOTES

†I want to comment extensively on this, for I feel today I learned something new... or at least improved my comprehension of the matter. I thought that the problem was that the second order diagram was proportional to U^2 instead of (U/D)^2. Afterall, perturbation theory is expected to come as an expansion on a small parameter and the proper small parameter here is indeed U/D or U/t, for sure not just U. Nevertheless I checked again the ref. by Haule and it was written as U^2.
Hence I thought a lot about the issue and finally realized that the 1/D^2 factor is already included in the diagram itself, for it comes from two nested convolutions of a spectral function, which values are indeed normalized depending on its support, i.e. D. If you increase D you decrease by D each value and so the result of the convolutions.

code/+phys/luttinger.m

@@ -1,60 +1,55 @@
-function [I,Lfig] = luttinger(w,sloc,gloc)
+function I = luttinger(w,sloc,gloc)
 %% 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
+%       w           : real valued array,    \omega domain
+%       sloc        : complex valued array, \Sigma(\omega)  function
 %       gloc        : complex valued array, G_{loc}(\omega) function
+%
 %  Output:
 %       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
 
-% 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);
-                                                                   if DEBUG
-Lfig = figure("Name",'Luttinger integrand','Visible','off');
-plot(w,integrand);
-xlabel('\omega');
-ylabel('Im[G(\omega)\partial\Sigma/\partial\omega]');
-ylim([-0.1,0.2]);
-                                                                   end
+    % 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);
+
 end
 
 
+

code/+phys/sopt.m

@@ -1,4 +1,4 @@
-function im_sloc = sopt(A,F,U)
+function iS = sopt(A,F,U)
 %% SOPT stands for Second Order Perturbation Theory:
 %  it computes the imag of the 2nd order diagram for the self-energy
 %
@@ -58,7 +58,7 @@
      D = D + flip(D); % flip{v(1:end)}=v(end:1)
 
   %% Imaginary part of the Self-Energy according to SOPT
-     im_sloc = -pi * U^2 * D;
+     iS = -pi * U^2 * D;
 
 end
 

code/dmft_loop.m

@@ -35,7 +35,7 @@
 %
 %% BSD 3-Clause License
 %
-%  Copyright (c) 2020, Gabriele Bellomia
+%  Copyright (c) 2020-2022, Gabriele Bellomia
 %  All rights reserved.
 
 if(~exist('pmode','var'))
@@ -55,7 +55,7 @@
     while LOOP
 
         % Weiss field from local Green's function
-            g0 = 1 ./ (w + eta - 0.25 .* gloc);
+            g0 = 1 ./ (w + eta - 0.25*D^2 .* gloc);
 
         % Spectral-function of Weiss field
             A0 = -imag(g0) ./ pi;

code/main.m

@@ -11,16 +11,16 @@
 end
 
 %% INPUT: Physical Parameters 
-D    = 1;               % Bandwidth
-U    = 0.1;             % On-site Repulsion
+D    = 4.0;             % Bandwidth
+U    = 4.0;             % On-site Repulsion
 beta = inf;             % Inverse Temperature
 
 %% INPUT: Boolean Flags
 MottBIAS     = 0;       % Changes initial guess of gloc (strongly favours Mott phase)
 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
+SPECTRAL     = 1;       % Controls plotting of spectral functions
 PLOT         = 1;       % Controls plotting of *all static* figures
 GIF          = 0;       % Controls plotting of *animated* figures
 UARRAY       = 0;       % Activates SLURM scaling of interaction values
@@ -32,11 +32,11 @@
 mloop = 1000;           % Max number of DMFT iterations 
 err   = 1e-5;           % Convergence threshold for self-consistency
 mix   = 0.10;           % Mixing parameter for DMFT iterations (=1 means full update)
-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.10;           % Hubbard U incremental step for phase diagrams
-Umax  = 6.00;           % Hubbard U maximum value for phase diagrams
+wres  = 2^13+1;         % Energy resolution in real-frequency axis
+wcut  = 6.00 * D;       % Energy cutoff in real-frequency axis
+Umin  = 0.00 * D;       % Hubbard U minimum value for phase diagrams
+Ustep = 0.10 * D;       % Hubbard U incremental step for phase diagrams
+Umax  = 6.00 * D;       % 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
 Tmax  = 5e-2;           % Temperature U maximum value for phase diagrams
@@ -94,8 +94,7 @@
         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),infig] = phys.luttinger(w,sloc{i},gloc{i});
-        plot.push_frame('luttok.gif',i,mloop,dt,infig);
+        I(i) = phys.luttinger(w,sloc{i},gloc{i});
         S(i) = phys.strcorrel(w,sloc{i});
     end
     if(PLOT)