From 785679981c77411b1881f9f30bec6cc3c3d6872a Mon Sep 17 00:00:00 2001 From: Gabriele Bellomia Date: Tue, 18 Jan 2022 08:22:47 +0100 Subject: [PATCH] Fix Luttinger integral MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 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. --- .gitignore | 2 +- README.md | 32 ++- code/+phys/LuttingerIntegral.m | 31 ++- code/+plot/Uline.m | 7 +- code/main.m | 14 +- showcase/ZvsL_beta1e6_wres2e15.svg | 325 +++++++++++++++++++++++++++++ 6 files changed, 378 insertions(+), 33 deletions(-) create mode 100644 showcase/ZvsL_beta1e6_wres2e15.svg diff --git a/.gitignore b/.gitignore index 4403b64..8fef34d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,4 @@ +output/ *.fig -*.gif *.asv *.sfit diff --git a/README.md b/README.md index cd98d32..85ae5e2 100644 --- a/README.md +++ b/README.md @@ -1,13 +1,14 @@ -# 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 _ZF_ , or the _quantized_ Luttinger integral _IL_, 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 :-------------------------:|:-------------------------:|:-------------------------: @@ -15,9 +16,18 @@ A few fancy examples of what you can obtain running some simulation within the M `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 index 4e6c6b1..cccb4e2 100644 --- a/code/+phys/LuttingerIntegral.m +++ b/code/+phys/LuttingerIntegral.m @@ -1,10 +1,6 @@ -%% 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 @@ -12,11 +8,24 @@ % 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 index e50b662..7b1d6d9 100644 --- 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 index 31910ba..23878b0 100644 --- 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 @@ - - - diff --git a/showcase/ZvsL_beta1e6_wres2e15.svg b/showcase/ZvsL_beta1e6_wres2e15.svg new file mode 100644 index 0000000..176eda9 --- /dev/null +++ b/showcase/ZvsL_beta1e6_wres2e15.svg @@ -0,0 +1,325 @@ + + +