单层磷烯10参数能带计算程序

坐标程序

function [x,y,z]=zigzag_pho(nx,ny)
u=0.08056;
v=0.10168;
a=3.314;
c=4.376;
b=10.48;
x1=zeros(1,4);
y1=zeros(1,4);
z1=zeros(1,4);
x1(1)=a/2;
x1(2)=0;
x1(3)=0;
x1(4)=a/2;
z1(1)=-v*b;
z1(2)=-v*b;
z1(3)=v*b;
z1(4)=v*b;
y1(1)=-(0.5-u)*c;
y1(2)=-u*c;
y1(3)=u*c;
y1(4)=(0.5-u)*c;
x2=x1;
y2=y1;
z2=z1;
for i=1:ny-1
    x2=[x2,x1];
    y2=[y2,y1+i*ones(1,4)*c];
    z2=[z2,z1];
end
x=x2;
y=y2;
z=z2;
n=length(x2);
for i=1:nx-1
    x=[x,x2+i*ones(1,n)*a];
    y=[y,y2];
    z=[z,z2];
end

哈密顿量程序

function H=Hamiltonian_NN_sp_phos(x,y,z)
N=length(x);
H=zeros(N,N);
for m=1:N
    for n=1:N
        d=sqrt((x(m)-x(n))^2+(y(m)-y(n))^2+(z(m)-z(n))^2);
        if abs(d-2.224)<0.01
            H(m,n)=-1.486;
        elseif abs(d-2.244)<0.01
            H(m,n)=3.729;
        elseif abs(d-3.314)<0.01
            H(m,n)=-0.252;
        elseif abs(d-3.334)<0.01
            H(m,n)=-0.071;
        elseif abs(d-3.475)<0.01
            H(m,n)=-0.019;
        elseif abs(d-4.245)<0.01
            H(m,n)=0.186;
        elseif abs(d-4.376)<0.01
            H(m,n)=-0.063;
        elseif abs(d-5.18)<0.01
            H(m,n)=0.101;
        elseif abs(d-5.37)<0.01
            H(m,n)=-0.042;
        elseif abs(d-5.49)<0.01
            H(m,n)=0.073;      
        end  
    end
end

计算程序

nx=3;  % 
ny=10; % 体系宽度（y方向的长度）
[x,y,z]=zigzag_pho(nx,ny);
plot3(x,y,z,'.','MarkerSize',20)
H=Hamiltonian_NN_sp_phos(x,y,z);
N=length(H);
HDL=H(N/3+1:N*2/3,1:N/3);
HD=H(N/3+1:N*2/3,N/3+1:N*2/3);
HDR=H(N/3+1:N*2/3,N*2/3+1:N);
n = length(HD);
dk = 0.001*pi*3/nx;
kx=-pi*3/nx:dk:pi*3/nx; % k空间路径
dN = length(kx);
Ek = zeros(n,dN);
for i = 1:dN
    Hk=HDL*exp(-1i*kx(i)*nx/3)+HD+HDR*exp(1i*kx(i)*nx/3);
    [~,E]=eig(Hk);
    Ek(:,i) = diag(E);
end
plot(kx,sort(real(Ek)),'.')
set(gca,'YLim',[-2.5 2.5]);
set(gca,'XLim',[-pi*3/nx pi*3/nx]);