主函数
posA=[0,0]; % 格点坐标
posB=[0,1]; % 格点坐标
a1=[sqrt(3),0]; % 晶格矢量
a2=[sqrt(3)/2,1.5]; % 晶格矢量
pos=[posA;posB]; % 格点所有坐标
a=1; % 晶格常数
n = 1000;
kx = linspace(-pi,pi,n);
ky = zeros(1,length(kx));
k=[kx;ky]';
Es=zeros(length(pos),length(k));
t=-1;
for i=1:length(kx)
H=Hamiltonian2Dk(pos,a1,a2,k(i,:),t,a);
[~,E]=eig(H);
Es(:,i)=diag(E);
end
plot(kx,real(Es),'.')
生成哈密顿量函数
function H=Hamiltonian2Dk(pos,a1,a2,k,t,a)
%pos 表示所有格点的坐标[x,y]
%a1,a2分别是两个晶格矢量
%t是最近邻hopping
%a是晶格常数
H =zeros(length(pos),length(pos));
eps=0.0001;
for i=1:length(pos)
for j=1:length(pos)
for n=-1:1:1
for m=-1:1:1
if abs(sqrt((pos(i,:)+n*a1+m*a2-pos(j,:))*(pos(i,:)+n*a1+m*a2-pos(j,:))')-a)<eps
H(i,j)=H(i,j)+t*exp(1i*k*(n*a1+m*a2)');
end
end
end
end
end