%program to plot fcc bandstructures and DoS calculated by quantum espresso
%@Florian Mayer

bands_datname='bands_Ge.dat';
scf_datname='Ge.scf.out';
dos_datname='dos_Ge.dat';
%choose an order of symmetry points(the order must be the same as in QE!!!)
order='KWXGLWU';

%fcc - brilluin zone parameters
ao=10;
G=[0,0,0];
X=(2*pi)/ao.*[1,0,0];
L=(2*pi)/ao.*[1/2,1/2,1/2];
W=(2*pi)/ao.*[1,1/2,0];
K=(2*pi)/ao.*[3/4,3/4,0];
U=(2*pi)/ao.*[1,1/4,1/4];
vectors=[G',X',L',W',K',U'];
%calculate lenghts of vectors
G_l=0;
X_l=norm(X);
L_l=norm(L);
W_l=norm(W);
K_l=norm(K);
U_l=norm(L);
length=[X_l,L_l,W_l,K_l,U_l];
%calculate entire length
total_length=G_l+X_l+L_l+W_l+K_l+U_l;

%nbnb=12;
fid=fopen(bands_datname);
line1=fgets(fid);
nbnb=str2num(line1(15:16));
s=1;
if nbnb<=10
    while ~feof(fid)
    fgets(fid);
    lines=fgets(fid);
    lines_t=str2num(lines);
    bands(s,:)=lines_t;
    s=s+1;
    end
else
    while ~feof(fid)
    fgets(fid);
    lines1=fgets(fid);
    lines_t1=str2num(lines1);
    lines2=fgets(fid);
    lines_t2=str2num(lines2);
    bands(s,:)=[lines_t1,lines_t2];
    s=s+1;
    end
end
    
fclose(fid);
%search for fermi energy:
fid2=fopen(scf_datname);
while ~feof(fid2)
    line1=fgets(fid2);
    if numel(line1)>=22
    if strcmp(line1(1:21),'     highest occupied')==1;
        
        disp(line1);
        line1_ev=line1;
    end
    end
end
fermi_energy=str2num(line1_ev(34:end));
fclose(fid2);

%sort bands
bands=sort(bands,2);

[nn,mm]=size(bands);
% order of symmetry points
order_label=cell(1,numel(order));
for t=1:numel(order)
    if strcmp(order(t),'G')==1
        order_label{t}='G';
    else
    order_label{t}=order(t);
    end
end

for t=1:numel(order)
    switch order(t)
        case 'G'
            vectors_n(:,t)=vectors(:,1);
        case 'X'
            vectors_n(:,t)=vectors(:,2);
        case 'L'
            vectors_n(:,t)=vectors(:,3);
        case 'W'
            vectors_n(:,t)=vectors(:,4);
        case 'K'
            vectors_n(:,t)=vectors(:,5);
        case 'U'
            vectors_n(:,t)=vectors(:,6);
    end
end

for t=1:numel(order)-1
    length_order(t)=norm(vectors_n(:,t)-vectors_n(:,t+1));
end

total_length_order=sum(length_order);
re_scale=length_order./total_length_order;
npoints=(nn-1)/(numel(order)-1);
start=0;
for t=1:numel(length_order)
    %x_dat(t,:)=re_scale(t).*linspace(start,stop,npoints);
    stop=start+re_scale(t)*(npoints);
    x_dat(t,:)=start:re_scale(t):stop;
    start=stop;
    %start=stop;
    %stop=stop+1;
end
x_dat_n=[];
for t=1:numel(length_order)
    x_dat_n=[x_dat_n,x_dat(t,:)];
end
x_dat_n=x_dat_n(1:numel(length_order)*npoints+1);
ticks_n1=x_dat(:,end);
ticks_n=[0,ticks_n1'];
for t=1:numel(length)
   resc(t)=length(t)/total_length; 
end

ylmax=max(bands(1:end))+1.5;
ylmin=min(bands(1:end))-1.5;
figure(2)
for t=1:nbnb
plot(x_dat_n,bands(:,t),'k','LineWidth',2);
hold on;
end
for t=1:numel(order)-2
plot([ticks_n(t+1),ticks_n(t+1)],[ylmin,ylmax],'--','Color',[0.514,0.545,0.545]);
hold on;
end
hold on;
plot([min(x_dat_n),max(x_dat_n)],[fermi_energy,fermi_energy],'-.','Color',[0.514,0.545,0.545]);
hold on;
xlim([min(x_dat_n),max(x_dat_n)+1]);
set(0,'defaulttextInterpreter','tex')
hold on;
set(gca,'XTick',ticks_n,'XTickLabel',order_label,'FontSize',13);
hold on;
ylabel('E / eV');
ylim([ylmin,ylmax]);
%plot(1:nn,6.1*ones(1,nn),'r');
%xlim([0,nn]);
hold off;

% plotting density of states
dos_name=dos_datname;
dos_data=importdata(dos_name);

fermi_x=dos_data.data(dos_data.data(:,1)<=fermi_energy);
[n_fermi,~]=size(fermi_x);

figure(11)
plot(dos_data.data(:,1),dos_data.data(:,2),'k','LineWidth',1.5);
hold on;
ar=area(dos_data.data(1:n_fermi,1),dos_data.data(1:n_fermi,2));
set(ar,'FaceColor',[0.514,0.545,0.545]);
hold on;
xlabel('Energy [eV]');
ylabel('D(E) [unitcell^{-1} eV^{-1}]');
xlim([min(dos_data.data(:,1)),max(dos_data.data(:,1))]);
