matlab 绘制石墨烯能带图(最近邻)(笔记)

刚学习lobster,石墨烯狄拉克点附近的能带是C-pz轨道构成的，比较简单可以入门哈。为了作对比我用紧束缚画了一个石墨烯的能带。我只考虑了C-C的最近邻（without SOC），当然次近邻和三近邻稍微改动也可以加进去。

                                                                       石墨烯俯视图

                                                                 石墨烯能带图

代码如下

%%Nie-Wei Wang,2024-1-8
clear all
%%%%%%%%晶格矢量 正格矢量 
a1=[2.4682260730186987   0.0000000000000000   0.0000000000000000];
a2=[-1.2341130365093496  2.1375473857986873  0.0000000000000000];
a3=[ 0.0000000000000006  0.0000000000000011   10.0000000000000018];
% C1=[0.0000018561517373 -0.0000018561604471  0.5000000000000000];
% C2=[0.6666681678482611  0.3333318611604476  0.5000000000000000];
%%%%%%%%没错就是POSCAR粘贴过来的
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
V=dot(a1,cross(a2,a3));%  %%%晶胞体积

b1=2*pi*cross(a2,a3)/V; %%%%倒格矢量
b2=2*pi*cross(a3,a1)/V;
b3=2*pi*cross(a1,a2)/V;

a1=a1(1:2);                     %晶格矢量a1 二维
a2=a2(1:2);                     %晶格矢量a2
 
b1=b1(1:2);                     %%%倒空间晶格矢量b1，二维
b2=b2(1:2);                     %%%倒格矢b2
 
gamma=[0 0];                    %%%高对称点gamma 
M=1/2*b1;                       %%%高对称点M    
K=1/3*b1+1/3*b2;                %%%高对称点K   


kn=20;     %两个高对称点之间分成有kn个点
% Gamma-M段  每个高对称点都有两个 150和151都是M点
kx(1:kn)=linspace(gamma(1),M(1),kn);
ky(1:kn)=linspace(gamma(2),M(2),kn);
% M-K段
kx(kn+1:2*kn)=linspace(M(1),K(1),kn);
ky(kn+1:2*kn)=linspace(M(2),K(2),kn);
% K-Gamma段
kx(2*kn+1:3*kn)=linspace(K(1),0,kn);
ky(2*kn+1:3*kn)=linspace(K(2),0,kn);
k=zeros(3*kn,2);
k(:,1)=kx;                            %%存储kx
k(:,2)=ky; 
Nw=2; %%%%能带条数需要修改
k_new=zeros(1,3*kn);    %%%为了沿着高对称点画出能带需要将k做调整,按照路径的模长走
energy=zeros(Nw,kn);    %%存放本征值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:3*kn   %%%%高对称路径
    if i<=kn
        kv=[kx(i)-gamma(1),ky(i)-gamma(2)];
        k_new(i)=norm(kv);                        %gamma->M      
    elseif i>kn && i<=2*kn
        kv=[kx(i)-M(1),ky(i)-M(2)];
        k_new(i)=norm(M-gamma)+norm(kv);           %M->K 
    else 
        kv=[kx(i)-K(1),ky(i)-K(2)];
        k_new(i)=norm(kv)+norm(M-gamma)+norm(K-M); %%%K->gamma 
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%路径不要动
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
efermi=0;
e1=0;
e2=0;
t=3;  %跃迁积分

Nw=2; %两条能带

h0=[e1,0;0,e2];
h1=[0,t;0 0];
%%%%%%%%%%%%%%%%最近邻
R1=-a1;
R2=-(a1+a2);
R3=-R1;
R4=-R2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:3*kn 
%%%%%%%%%%%%%%%normal   
% hsum=h0+(exp(1i*dot(k(i,:),R1))+exp(1i*dot(k(i,:),R2)))*h01...
%        +(exp(1i*dot(k(i,:),-R1))+exp(1i*dot(k(i,:),-R2)))*h01'
%%%%%%%%%%%%%%simple 其实是对normal的简化，
h=(1+exp(1i*dot(k(i,:),R1))+exp(1i*dot(k(i,:),R2)))*t;%%最近邻
hsum=h0+[0,h;h',0];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
energy(:,i)=eig(hsum);
end
 
%%%%%%%plot 画图
for m=1:Nw                                 %%%  原胞中有几个轨道就有几条能带
    plot(k_new,energy(m,:)-efermi,'color','b','LineWidth',1)  %%
    hold on
end
 
 
%%%%设置xticks布里渊区高对称点
xticks([k_new(1) k_new(kn+1)  k_new(2*kn+1)  k_new(3*kn)]);
xticklabels({'\Gamma','M','K','\Gamma'})
 
%%%辅助线
xline(k_new(kn+1),'--');
xline(k_new(2*kn+1),'--');
yline(0,'--');
%grid on
 
%%%
xlim([k_new(1) k_new(3*kn)])
ylim([-10,10])
ylabel("Energy")
box on

上次写好了找不到了......，今天重新写了一个。