Skip to content

Commit

Permalink
Refactor MFenergies
Browse files Browse the repository at this point in the history 
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).
  • Loading branch information
@beddalumia
beddalumia committed Aug 3, 2022
1 parent 189c4c2 commit 257fe7d
Show file tree
Hide file tree
Showing 3 changed files with 191 additions and 120 deletions.
221 changes: 122 additions & 99 deletions run/MFenergies.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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 <Uloc*nup*ndw>
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 <Uloc*nup*ndw>

if strcmp(mode,'line')
figure('Name',whichMF{iMF})
Expand All @@ -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('..')
Expand Down Expand Up @@ -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')
Expand All @@ -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
4 changes: 3 additions & 1 deletion run/PostMF.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading

0 comments on commit 257fe7d

Please sign in to comment.