From 257fe7d11f088a0883b8ffb3dcb5f26e6953c15b Mon Sep 17 00:00:00 2001 From: Gabriele Bellomia Date: Wed, 3 Aug 2022 08:51:09 +0200 Subject: [PATCH] Refactor MFenergies Now the < >< > energy terms are computed directly in the fortran source, so we just read them (the explicit matlab computation is moved in a sub function, as a check that everything works). --- run/MFenergies.m | 221 +++++++++++++++++++++++++--------------------- run/PostMF.m | 4 +- src/mf_km_afm.f90 | 86 +++++++++++++----- 3 files changed, 191 insertions(+), 120 deletions(-) diff --git a/run/MFenergies.m b/run/MFenergies.m index 0e5159d..a2c534d 100644 --- a/run/MFenergies.m +++ b/run/MFenergies.m @@ -2,10 +2,9 @@ %% INPUT -doSmooth = 0 % true | false -> Smoothing on energy differences -doBands = false; % true | false --> plot mean-field bandstructure -doVector = false; % true | false --> save bands to vectorized file -doRaster = false; % true | false --> save bands to rasterized file +doSmooth = 0; % true | false --> Smoothing on energy differences +doBands = 0; % true | false --> plot mean-field bandstructure +doCheck = 1; % true | false --> check formula for < >< > terms mode = 'map'; % 'line' | 'map' --> many lines or phase-diagram @@ -33,100 +32,12 @@ fprintf(lineID{:}); fprintf('\n••••••••••••\n'); cd(lineID{:}); - clear('ids','ordpms','U_list'); - load('order_parameter_line.mat','ids','ordpms','U_list'); - Emfb(iSOI,:) = load('mean_field_gs.txt'); % mean-field bandstructure - try - Ekin(iSOI,:) = load('kinetic_energy.txt'); % dmft-like kinetic energy - catch - Ekin(iSOI,:) = NaN; - end - Npar = length(ordpms); - Npoints = length(U_list); - for iHubb = 1:Npoints - U = U_list(iHubb); - UDIR= sprintf('U=%f',U); - cd(UDIR); - % - Emag(iSOI,iHubb) = 0; - % - Sx = 0; - Sy = 0; - Sz = 0; - % - Rx = 0; - Ry = 0; - Rz = 0; - % - Pz = 0; - Ne = 2; % --> In the *unit-cell* - % - if Npar == 2 - Sz = ordpms{1}(iHubb); - Rz = ordpms{2}(iHubb); - elseif Npar == 4 - Sx = ordpms{1}(iHubb); - Sy = ordpms{2}(iHubb); - Rx = ordpms{3}(iHubb); - Ry = ordpms{4}(iHubb); - elseif Npar == 6 - Sx = ordpms{1}(iHubb); - Sy = ordpms{2}(iHubb); - Sz = ordpms{3}(iHubb); - Rx = ordpms{4}(iHubb); - Ry = ordpms{5}(iHubb); - Rz = ordpms{6}(iHubb); - end - % - fprintf('Computing "< >< >" term for U=%.2f ..',U); - Emag(iSOI,iHubb) = Emag(iSOI,iHubb) + (Sz+Rz)^2 - (Pz+Ne)^2 + (Sz-Rz)^2 - (Pz-Ne)^2; - Emag(iSOI,iHubb) = Emag(iSOI,iHubb) + (Sx+Rx)^2 + (Sy+Ry)^2 + (Sx-Rx)^2 + (Sy-Ry)^2; - Emag(iSOI,iHubb) = Emag(iSOI,iHubb) * U/32; - Emag(iSOI,iHubb) = Emag(iSOI,iHubb) + U/04; % Hartree-Fock correction (<-> Ne=0) - % - fprintf('.DONE!\n'); - %% Plotting Bands - %------------------------------------------------------------------ - if doBands - Eigenbands = load('EigenbandsKM.mf'); - id = sprintf('Mean-Field Bands [SOI=%f U=%f]',SOI,U); - figure("Name",id); - Ncell = length(Eigenbands); - Ev = Eigenbands(1:floor(Ncell/2),:); - Ec = Eigenbands(floor(Ncell/2)+1:end,:); - scatter(Ev(:,1),Ev(:,2),'r'); hold on - scatter(Ec(:,1),Ec(:,2),'b'); - % Title, legend, all of that - title(id); - xticks([0 2.418399 4.836798 7.245524]) - xlim([0,7.245524]); - xticklabels({'\Gamma','K','K`','\Gamma'}) - ylabel('\epsilon / t'); - ax = gca; ax.Box = 'on'; - % [These two lines ensure filling of the fig] - InSet = get(ax, 'TightInset'); - set(ax, 'Position', [InSet(1:2), 1-InSet(1)-InSet(3),... - 1-InSet(2)-InSet(4)]); - % PRINTING (if requested) - if doRaster - filename = sprintf('MF_bands_SOI=%f_U=%f.png',SOI,U); - fprintf('Printing %s..',filename); - print(gcf,filename,'-dpng','-r600'); - elseif doVector - filename = sprintf('MF_bands_SOI=%f_U=%f.pdf',SOI,U); - fprintf('Printing %s..',filename); - print(gcf,filename,'-dpdf','-fillpage'); - else - fprintf('SKIP printing!\n'); - end - fprintf('.DONE!\n'); - end - cd('..'); - end - cd('..'); - - Etot(iSOI,:) = Emfb(iSOI,:) + Emag(iSOI,:); % mean-field bands + order parameters - Epot(iSOI,:) = Etot(iSOI,:) - Ekin(iSOI,:); % should correspond to + Emfb(iSOI,:) = load('mf_bands_energy.txt'); % mean-field bandstructure + Ekin(iSOI,:) = load('kinetic_energy.txt'); % dmft-like kinetic energy + Epot(iSOI,:) = load('potential_energy.txt'); % matsubara potental energy + Emag(iSOI,:) = load('magnetic_energy.txt'); % magnetic ordering energy + Etot(iSOI,:) = Emfb(iSOI,:) + Emag(iSOI,:); % mean-field bands + order parameters + Upot(iSOI,:) = Etot(iSOI,:) - Ekin(iSOI,:); % should correspond to if strcmp(mode,'line') figure('Name',whichMF{iMF}) @@ -136,6 +47,18 @@ plot(U_list,Epot(iSOI,:)) end + if doBands + plot_bands(SOI) + end + + if doCheck + check_magnetic_energy(Emag(iSOI,:)); + end + + U_list = postDMFT.get_list('U'); + + cd('..'); + end cd('..') @@ -167,7 +90,7 @@ end -if size(AFMxy) == size(AFMz) +if strcmp(mode,'map') & size(AFMxy) == size(AFMz) figure('Name','GSenergy difference @ Hartree-Fock') title('Groundstate Energy Difference [AFMxy - AFMz]','Interpreter','latex') xlabel('$U/t$','Interpreter','latex') @@ -189,3 +112,103 @@ +%% contains + +function plot_bands(SOI) + doVector = false; % true | false --> save bands to vectorized file + doRaster = false; % true | false --> save bands to rasterized file + U_list = postDMFT.get_list('U'); + Npoints = length(U_list); + for iHubb = 1:Npoints + U = U_list(iHubb); + UDIR= sprintf('U=%f',U); + cd(UDIR); + %% Plotting Bands + Eigenbands = load('EigenbandsKMH.mf'); + id = sprintf('Mean-Field Bands [SOI=%f U=%f]',SOI,U); + figure("Name",id); + Ncell = length(Eigenbands); + Ev = Eigenbands(1:floor(Ncell/2),:); + Ec = Eigenbands(floor(Ncell/2)+1:end,:); + scatter(Ev(:,1),Ev(:,2),'r'); hold on + scatter(Ec(:,1),Ec(:,2),'b'); + % Title, legend, all of that + title(id); + xticks([0 2.418399 4.836798 7.245524]) + xlim([0,7.245524]); + xticklabels({'\Gamma','K','K`','\Gamma'}) + ylabel('\epsilon / t'); + ax = gca; ax.Box = 'on'; + % [These two lines ensure filling of the fig] + InSet = get(ax, 'TightInset'); + set(ax, 'Position', [InSet(1:2), 1-InSet(1)-InSet(3),... + 1-InSet(2)-InSet(4)]); + % PRINTING (if requested) + if doRaster + filename = sprintf('MF_bands_SOI=%f_U=%f.png',SOI,U); + fprintf('Printing %s..',filename); + print(gcf,filename,'-dpng','-r600'); + elseif doVector + filename = sprintf('MF_bands_SOI=%f_U=%f.pdf',SOI,U); + fprintf('Printing %s..',filename); + print(gcf,filename,'-dpdf','-fillpage'); + else + fprintf('SKIP printing!\n'); + end + fprintf('.DONE!\n'); + cd('..'); + end +end + +function Emag = check_magnetic_energy(Emag_ref) + load('order_parameter_line.mat','ids','ordpms','U_list'); + Npar = length(ordpms); + Npoints = length(U_list); + for iHubb = 1:Npoints + U = U_list(iHubb); + UDIR= sprintf('U=%f',U); + cd(UDIR); + % + Emag(iHubb) = 0; + % + Sx = 0; + Sy = 0; + Sz = 0; + % + Rx = 0; + Ry = 0; + Rz = 0; + % + Pz = 0; + Ne = 2; % --> In the *unit-cell* + % + if Npar == 2 + Sz = ordpms{1}(iHubb); + Rz = ordpms{2}(iHubb); + elseif Npar == 4 + Sx = ordpms{1}(iHubb); + Sy = ordpms{2}(iHubb); + Rx = ordpms{3}(iHubb); + Ry = ordpms{4}(iHubb); + elseif Npar == 6 + Sx = ordpms{1}(iHubb); + Sy = ordpms{2}(iHubb); + Sz = ordpms{3}(iHubb); + Rx = ordpms{4}(iHubb); + Ry = ordpms{5}(iHubb); + Rz = ordpms{6}(iHubb); + end + % + fprintf('Computing "< >< >" term for U=%.2f ..',U); + Emag(iHubb) = Emag(iHubb) + (Sz+Rz)^2 - (Pz+Ne)^2 + (Sz-Rz)^2 - (Pz-Ne)^2; + Emag(iHubb) = Emag(iHubb) + (Sx+Rx)^2 + (Sy+Ry)^2 + (Sx-Rx)^2 + (Sy-Ry)^2; + Emag(iHubb) = Emag(iHubb) * U/32; + Emag(iHubb) = Emag(iHubb) + U/04; % Hartree-Fock correction (<-> Ne=0) + % + fprintf('.DONE!\n'); + % + fprintf(' REF = %f\n DIFF = %f\n',Emag_ref(iHubb),Emag(iHubb)-Emag_ref(iHubb)); + % + cd('..'); + end +end \ No newline at end of file diff --git a/run/PostMF.m b/run/PostMF.m index cf3f8db..f32ddae 100644 --- a/run/PostMF.m +++ b/run/PostMF.m @@ -31,7 +31,9 @@ end [ids,ordpms] = postDMFT.order_parameter_line(U_list); fprintf('...DONE\n'); save('order_parameter_line.mat','ids','ordpms','U_list'); - mfens = postDMFT.custom_line('mean_field_gs.dat',U_list); + postDMFT.custom_line('mf_bands_energy.dat',U_list); + postDMFT.custom_line('potential_energy.dat',U_list); + postDMFT.custom_line('magnetic_energy.dat',U_list); try postDMFT.kinetic_line(); % Builds the 'kinetic_energy.txt' file catch diff --git a/src/mf_km_afm.f90 b/src/mf_km_afm.f90 index 931d535..31c50e1 100644 --- a/src/mf_km_afm.f90 +++ b/src/mf_km_afm.f90 @@ -7,7 +7,7 @@ program mf_km_2d integer :: Nparams !#{order_parameters} integer,parameter :: Norb=1,Nspin=2,Nlat=2,Nlso=Nlat*Nspin*Norb integer :: Nk,Nktot,Nkpath,Nkx,Nky,L,Ltail - integer :: unit_p,unit_e + integer :: unit_pms,unit_mfe,unit_mag integer :: i,j,k,ik,iorb,jorb,ispin,io,jo integer :: ilat,jlat integer :: ix,iy,iz @@ -163,18 +163,21 @@ program mf_km_2d call save_array("params.init",params) !Save used initial parameters, for reproducibility !SETUP I/O FILES - unit_p = free_unit() + unit_pms = free_unit() select case(Nparams) case(2) - open(unit_p,file="order_parameters_Sz_Rz.dat") + open(unit_pms,file="order_parameters_Sz_Rz.dat") case(4) - open(unit_p,file="order_parameters_Sx_Sy_Rx_Ry.dat") + open(unit_pms,file="order_parameters_Sx_Sy_Rx_Ry.dat") case(6) - open(unit_p,file="order_parameters_Sx_Sy_Sz_Rx_Ry_Rz.dat") + open(unit_pms,file="order_parameters_Sx_Sy_Sz_Rx_Ry_Rz.dat") end select ! - unit_e = free_unit() - open(unit_e,file="mean_field_gs.dat") + unit_mfe = free_unit() + open(unit_mfe,file="mf_bands_energy.dat") + ! + unit_mag = free_unit() + open(unit_mag,file="magnetic_energy.dat") !SELF-CONSISTENT LOOP converged=.false. ; iter=0 @@ -191,12 +194,17 @@ program mf_km_2d ! call end_loop end do - call save_array("params.restart",params) - close(unit_p) - close(unit_e) ! call save_array("Hmf_correction",Hmf_glob) call TB_write_hloc(Hmf_glob) !logging... + ! + call save_array("params.restart",params) + close(unit_pms) !Same content of restart + close(unit_mfe) !file but obs convention + ! + !Compute < >< > energy terms from ordpms + call mf_magnetic_energy(params) + close(unit_mag) !GREEN'S FUNCTIONS AND RELATED QUANTITIES if(withgf)then @@ -216,8 +224,8 @@ program mf_km_2d call dmft_print_gf_realaxis(Greal,"Greal",iprint=4) !Build mean-field self-energies do i=1,L - Smats(:,:,:,:,:,i) = lso2nnn(Hmf_glob) !We assume ∑(z) to be just a number, constant - Sreal(:,:,:,:,:,i) = lso2nnn(Hmf_glob) !for all z values: ∑(0) = ∑(∞) <=> Hmf = Htop + Smats(:,:,:,:,:,i) = dreal(lso2nnn(Hmf_glob)) !We assume ∑(z) to be just a number, constant + Sreal(:,:,:,:,:,i) = dreal(lso2nnn(Hmf_glob)) !for all z values: ∑(0) = ∑(∞) <=> Hmf = Htop enddo Sigma_lso = nnn2lso(Smats(:,:,:,:,:,1)) write(*,*) " " @@ -329,8 +337,8 @@ subroutine solve_MF_loc(order_parameters,mf_correction) enddo !Print MF-GS energy write(*,*)"∫_{BZ}Eᵥ(k)dk = "//str(Emf) - rewind(unit_e) - write(unit_e,"(F21.12)")Emf + rewind(unit_mfe) + write(unit_mfe,"(F21.12)")Emf ! !Update order parameters and print to stdout + file Sx = Sx + sum( GammaSx*rhoH ) @@ -345,20 +353,20 @@ subroutine solve_MF_loc(order_parameters,mf_correction) order_parameters = [Sz,Rz] write(*,*)"Sz Rz:" write(*,"(2F21.12)")Sz,Rz - rewind(unit_p) - write(unit_p,"(2F21.12)")Sz,Rz + rewind(unit_pms) + write(unit_pms,"(2F21.12)")Sz,Rz case(4) order_parameters = [Sx,Sy,Rx,Ry] write(*,*)"Sx Sy Rx Ry:" write(*,"(4F21.12)")Sx,Sy,Rx,Ry - rewind(unit_p) - write(unit_p,"(4F21.12)")Sx,Sy,Rx,Ry + rewind(unit_pms) + write(unit_pms,"(4F21.12)")Sx,Sy,Rx,Ry case(6) order_parameters = [Sx,Sy,Sz,Rx,Ry,Rz] write(*,*)"Sx Sy Sz Rx Ry Rz:" write(*,"(12F21.12)")Sx,Sy,Sz,Rx,Ry,Rz - rewind(unit_p) - write(unit_p,"(6F21.12)")Sx,Sy,Sz,Rx,Ry,Rz + rewind(unit_pms) + write(unit_pms,"(6F21.12)")Sx,Sy,Sz,Rx,Ry,Rz end select ! !Return mean-field correction to H(k) to parent scope @@ -390,6 +398,43 @@ function mf_Hk_correction(a) result(HkMF) HkMF = -HkMf * Uloc/4d0 ! end function mf_Hk_correction + ! + !Compute < >< > energy contributions from given order parameters + subroutine mf_magnetic_energy(a) + real(8),dimension(:) :: a !input array of order parameters + real(8) :: Emag !output < >< > "magnetic" energy + real(8) :: Sx,Sy,Sz,Rx,Ry,Rz !order parameters + ! + Sx = 0d0; Sy = 0d0; Sz = 0d0 + Rx = 0d0; Ry = 0d0; Rz = 0d0 + ! + Emag = 0d0 + ! + select case(Nparams) + case(2) + Sz = a(1) + Rz = a(2) + case(4) + Sx = a(1); Sy = a(2) + Rx = a(3); Ry = a(4) + case(6) + Sx = a(1); Sy = a(2); Sz = a(3) + Rx = a(4); Ry = a(5); Rz = a(6) + end select + ! + Emag = Emag + (Sx+Rx)**2 + (Sx-Rx)**2 + Emag = Emag + (Sy+Ry)**2 + (Sy-Ry)**2 + Emag = Emag + (Sz+Rz)**2 + (Sz-Rz)**2 + ! + Emag = Emag * Uloc/32d0 + ! + !Print magnetic energy + write(*,*)"< >< > = "//str(Emag) + rewind(unit_mag) + write(unit_mag,"(F21.12)")Emag + ! + end subroutine mf_magnetic_energy + !RESHAPE FUNCTIONS ! > nnn2lso: [Nlat,Nspin,Nspin,Norb,Norb] array -> [Nlso,Nlso] matrix @@ -436,6 +481,7 @@ function lso2nnn(Fin) result(Fout) enddo end function lso2nnn + !POTENTIAL ENERGY FROM MIGDAL-GALITSKIJ FORMULA IN MATSUBARA DOMAIN !The implementation is based on Matsubara formalism, moving from !Fetter-Walecka eq. 23.14 and trasforming to imaginary frequency