Siesta关于E-k关系计算结果保存于systemlabel.bands文件内，故Matlab程序只需要读取此文档即可。经笔者测试，在siesta-2.0.2与siesta-3.2版本内均可完美运行。

function sband( filename, E_range)

%sband('systemlabel.bands',[RangeOfBandEnergy])

%Read the file 'systemlabel.bands' and plot the band structure from siesta.

%The first argument is a must, while the E_range is of your own choice.

%Edited by JackyTu, July 29, 2014.

%Please contact tuxingchen@pku.edu.cn

fid=fopen(filename);

frewind(fid);

% -------------------- Read --------------------

E_fermi=fscanf(fid,'%f',1);

kmin=fscanf(fid,'%f',1);

kmax=fscanf(fid,'%f',1);

Emin=fscanf(fid,'%f',1);

Emax=fscanf(fid,'%f',1);

nband=fscanf(fid,'%d',1);

spin=fscanf(fid,'%d',1);

k_num=fscanf(fid,'%d',1);

if(nargin<2)

E_range=[Emin,Emax];

end

kp=zeros(k_num,1);

e_data=zeros(nband*spin,k_num);

for i=1:k_num

kp(i)=fscanf(fid,'%f',1);

e_data(:,i)=fscanf(fid,'%f',nband*spin);

end

E_k=zeros(nband,k_num,2);

if spin==2

E_k(:,:,1)=e_data(1:nband,:);

E_k(:,:,2)=e_data((nband+1):end,:);

else

E_k(:,:,1)=e_data;

end

n=fscanf(fid,'%d',1);

fgetl(fid);

for j=1:n

str=fgetl(fid);

str_out=strsplit(strtrim(str));

xtick(j)=str2double(str_out{1});

xtickLabel{j}=str_out{2}(2:end-1);

end

fclose(fid);

% -------------------- Plot --------------------

figure;

if spin==2

subplot(1,2,1);

hold on;

for ii=1:nband

plot(kp,E_k(ii,:,1)-E_fermi,'b');

end

for jj=2:n-1

plot([xtick(jj),xtick(jj)],ylim,'k')

end

set(gca, 'xtick', xtick);

set(gca, 'xticklabel', xtickLabel);

ylabel('E-E_F(eV)');

axis([xtick(1),xtick(end),E_range]);

plot(xlim,[0,0],'--k');

hold off;

subplot(1,2,2);

hold on;

for ii=1:nband

plot(kp,E_k(ii,:,2)-E_fermi,'r');

end

for jj=2:n-1

plot([xtick(jj),xtick(jj)],ylim,'k')

end

set(gca, 'xtick', xtick);

set(gca, 'xticklabel', xtickLabel);

ylabel('E-E_F(eV)');

axis([xtick(1),xtick(end),E_range]);

plot(xlim,[0,0],'--k');

hold off;

else

hold on;

for ii=1:nband

plot(kp,E_k(ii,:,1)-E_fermi,'b');

end

for jj=2:n-1

plot([xtick(jj),xtick(jj)],ylim,'k')

end

set(gca, 'xtick', xtick);

set(gca, 'xticklabel', xtickLabel);

ylabel('E-E_F(eV)');

axis([xtick(1),xtick(end),E_range]);

plot(xlim,[0,0],'--k');

hold off;

end

end

下面是附有运行sband函数后的效果图(此处计算体系为BCC-Fe)：