Fix Luttinger integral

- There was a double counting of the dw factor (silly me 🙈) - A very high wres value is needed to obtain sharp transitions > Updated the README and the showcase directory > Now plot.Uline() takes in also the Luttinger integral values, to plot them along the quasiparticle renormalization weights.
beddalumia · Jan 18, 2022 · 7856799 · 7856799
1 parent 8aa4104
commit 7856799
Show file tree

Hide file tree

Showing 6 changed files with 378 additions and 33 deletions.
diff --git a/.gitignore b/.gitignore
@@ -1,4 +1,4 @@
+output/
 *.fig
-*.gif
 *.asv
 *.sfit
diff --git a/README.md b/README.md
@@ -1,23 +1,33 @@
-# What MOTTlab is
+## What `MOTTlab` is
 A didactic/explorative implementation of (real axis) IPT-based dynamical mean-field theory for the half-filled Mott-Hubbard transition on the Bethe lattice, within MATLAB.
 
-# License and stuff
-This code has been implemented taking inspiration from these two didactic sources:
-1. http://www.physics.rutgers.edu/~haule/681/Perturbation.pdf ([local copy](docs/haule_IPTtheory_rutgers.pdf))
-2. https://www.cond-mat.de/events/correl19/manuscripts/rozenberg.pdf ([local copy](docs/rozenberg_review_julich.pdf))
-and the hands-on material given therein, in particular a tutorial-intended jupyter notebook provided by Óscar Nájera (available [on his cloud](http://mycore.core-cloud.net/index.php/s/oAz0lIWuBM90Gqt), or [locally](legacy/PYTHON/real_ipt-text_v3.ipynb)) under the BSD 3-Clause License. Here we provide a Matlab rewrite, plus some extensions (phase diagram loops, including convergence and self-mixing control, various post-processing and data analysis routines), again under BSD 3-Clause License. You can read more about allowed use of this code in the LICENSE file.
+## What you can get out of it
+Here we present a few examples of what you can obtain spending some time on the `MOTTlab`.
 
-A few fancy examples of what you can obtain running some simulation within the MOTTlab: visualizing the interaction-driven Mott transition at different relevant temperatures!
+> Characterize the interaction-driven quantum Mott transition through different physical markers, such as the _timeless classic_ quasiparticle renormalization weight _Z<sub>F</sub>_ , or the _quantized_ Luttinger integral _I<sub>L</sub>_, which gives some exotic hints of a topological interpretation.
+
+![Mott-Transition-Markers-zeroT](./showcase/ZvsL_beta1e6_wres2e15.svg) 
+
+> Truly visualize what happens at different relevant temperatures by means of the fanciest animations!
 
  Temperature  | DOS | SELF-ENERGY
 :-------------------------:|:-------------------------:|:-------------------------:
 `T -> 0: 2nd order MIT` |![Mott-Transition-AnimatedDOS-zeroT](./showcase/uDOS_zeroT.gif) | ![Mott-Transition-AnimatedSIGMA-zeroT](./showcase/uSigma_zeroT.gif) 
 `Intermediate T: 1st order MIT`  |![Mott-Transition-AnimatedDOS-intermediateT](./showcase/uDOS_beta50.gif) | ![Mott-Transition-AnimatedSIGMA-intermediateT](./showcase/uSigma_beta50.gif)
 `High T: supercritical MIT` |![Mott-Transition-AnimatedDOS-highT](./showcase/uDOS_beta1.gif) | ![Mott-Transition-AnimatedSIGMA-highT](./showcase/uSigma_beta1.gif)  
 
-# TODO
-- [ ] Insert a "restarting" protocol for full phase diagram spans. The gloc0=0 condition appers to be too unstable to obtain accurate UC1 lines. 
+## Licensing and legacy code
+This code has been implemented taking inspiration from these two didactic sources:
+1. http://www.physics.rutgers.edu/~haule/681/Perturbation.pdf ([local copy](docs/haule_IPTtheory_rutgers.pdf))
+2. https://www.cond-mat.de/events/correl19/manuscripts/rozenberg.pdf ([local copy](docs/rozenberg_review_julich.pdf))
+and the hands-on material given therein, in particular a tutorial-intended jupyter notebook provided by Óscar Nájera (available [on his cloud](http://mycore.core-cloud.net/index.php/s/oAz0lIWuBM90Gqt), or [locally](legacy/PYTHON/real_ipt-text_v3.ipynb)) under the BSD 3-Clause License. Here we provide an extensive MATLAB rewrite, augmented by basic phase-diagram workflows, with convergence and self-mixing control, various post-processing tools, finally some exploration of novel material, with explicit references to the relevant research papers. Everything again under BSD 3-Clause License. You can read more about allowed use of this code in the [LICENSE](./LICENSE) file.
+
+-------------
+
+## TODO
 - [x] Compute a MIvsFL marker based on the divergence of Im[Sigma(0)] (much sensible, very inexpensive), so to obtain sharper phase diagrams with respect to the Z-derived ones.  
  At the moment we implement the "strenght of correlations", S = norm[Sigma(0)-Sigma_HF], as defined in `PRL 114 185701`, with actual neat results: the marker is almost zero accross the whole FL phase and starts increasing fast in the Mott insulator.
-- [ ] Compute the Luttinger Integral, as defined in `PRB 102 081110(R)`. Since it appears to be quantized it could become the definitive _flag_ for phase diagrams; much better than Z or S for it is an _integer_.  
-UPDATE: Luttinger Theorem currently works for the FL phase *only*; the MI integer result for now is missing, for it arises from the nonanalytic Mott pole in the self-energy, which requires the most careful numerical treatment.  
+- [x] Compute the Luttinger Integral, as defined in `PRB 102 081110(R)`. Since it appears to be quantized it could become the definitive _flag_ for phase diagrams; much better than Z or S for it is an _integer_.  
+UPDATE: Luttinger Theorem currently works for very low temperatures only. It may well be an inherent limitation (in that case it would be a good marker for low temperature lines, but not phase diagrams). Also note that the quality of the first order step at the transition highly depends on the frequency resolution, so to have sharp boundaries you need to have at least `wres=2^14`, which makes the IPT solver far slower.
+- [ ] Insert a "restarting" protocol for lines and full phase diagram spans. The gloc0=0 condition appers to be too unstable to obtain accurate UC1 values. Furthermore this would most probably speed up a lot the convergence (lowering the required number of iterations), so to enable a relevant slowdown of the single iteration (e.g. bigger frequency resolution).
+
diff --git a/code/+phys/LuttingerIntegral.m b/code/+phys/LuttingerIntegral.m
@@ -1,22 +1,31 @@
-%% BSD 3-Clause License
-% 
-% Copyright (c) 2020, Gabriele Bellomia
-% All rights reserved.
-
-function I = LuttingerIntegral(w,sloc,gloc)
+function [I,Lfig] = LuttingerIntegral(w,sloc,gloc)
 %% Computes the Luttinger sum-rule, as defined in PRB 102 081110
+%
 %  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
-dw  = w(2)-w(1);
-ds  = diff(sloc);
+%
+%% BSD 3-Clause License
+% 
+%  Copyright (c) 2020, Gabriele Bellomia
+%  All rights reserved.
+                                                             global DoDEBUG
+ds  = diff(sloc); % dSigma -> already eliminates dw!
 g_  = gloc(1:end-1);
-integrand  = smoothdata(imag(g_.*ds));
-I = 1/pi*(dw*sum(integrand));
-%plot(w_,integrand);
+integrand  = imag(g_.*ds);
+I = 1/pi*sum(integrand); 
+I = abs(I); % J. Phys.: Condens. Matter 28 (2016) 025601
+                                                                 if DoDEBUG
+Lfig = figure("Name",'Luttinger integrand','Visible','off');
+plot(w(1:end-1),integrand);
+xlabel('\omega');
+ylabel('Im[G(\omega)\partial\Sigma/\partial\omega]');
+ylim([-0.1,0.2]);
+                                                                 end
 end
 
+
diff --git a/code/+plot/Uline.m b/code/+plot/Uline.m
@@ -1,10 +1,13 @@
-function line = Uline(Z,beta,Umin,Ustep,Umax)
+function line = Uline(Z,I,beta,Umin,Ustep,Umax)
     line = figure("Name",'U-driven MIT');
     U = Umin:Ustep:Umax;
     scatter(2.*U,Z,'k','filled'); hold on
+    scatter(2.*U,I,'r','filled'); % (Units: D=2t)
     title(sprintf('IPT  |  U-driven MIT at T = %f',1/beta))
     xlabel('$U/t$','Interpreter','latex')
-    ylabel('$Z = M/M^*$','Interpreter','latex')
+    legend({'$Z_F = 1 / ( 1 - Re[\partial_\omega\Sigma(0)] )$',...
+        '$I_L = 1/\pi\int d\omega Im[G(\omega)\partial_\omega\Sigma(\omega)]$'},...
+        'Interpreter','latex')
 end
 
 
diff --git a/code/main.m b/code/main.m
@@ -9,22 +9,23 @@
 %% INPUT: Physical Parameters 
 D    = 1;               % Bandwidth
 U    = 5;               % On-site Repulsion
-beta = 50;              % Inverse Temperature
+beta = 1e3;             % Inverse Temperature
 
 %% 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
-DoSPECTRAL   = 1;       % Controls plotting of spectral functions
+DoSPECTRAL   = 0;       % Controls plotting of spectral functions
 DoPLOT       = 1;       % Controls plotting of *all static* figures
 DoGIF        = 0;       % Controls plotting of *animated* figures
+DoDEBUG      = 0;       % Activates debug prints / plots / operations
 
 %% INPUT: Control Parameters
 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^12;           % Energy resolution in real-frequency axis
+wres  = 2^15;           % 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
@@ -73,7 +74,7 @@
         U = U + Ustep;
     end
     if(DoPLOT)
-        u_span = plot.Uline(Z,beta,Umin,Ustep,Umax)
+        u_span = plot.Uline(Z,I,beta,Umin,Ustep,Umax)
     end
     if(DoGIF && DoSPECTRAL)
         plot.spectral_gif(w,gloc,sloc,Umin:Ustep:Umax,1/beta,dt)
@@ -132,6 +133,3 @@
 
 
 
-
-
-
-Original file line number
+Diff line change
@@ -1,4 +1,4 @@
+    output/
     *.fig
-    *.gif
     *.asv
     *.sfit