二维石墨烯能带全数值计算代码

主函数

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