24小时热门版块排行榜

返回列表

wei.yc

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 3
帖子: 2
在线: 2.3小时
虫号: 3699170
注册: 2015-02-27
专业: 凝聚态物性 II ：电子结构

[交流] 关于有转角双层石墨烯14带能级连续性模型的程序求助~

本人想重复PNAS里的Moiré bands in twisted double-layer graphene（http://www.pnas.org/content/108/30/12233. full.pdf）文章里的结果（Fig.3A），但我一直是做实验的，编程能力不强，自己的程序（如下）与文章结果相差甚远，故想向各位请求指点，谢谢~

程序：
Program:

clear
w=0.11;
a=1.42;
t=3.16;
v=t*sqrt(3)*a/2;
kd=4*pi/(3*sqrt(3)*a);
theta=1.05*pi/180;
phi=2*pi/3;
ktheta=kd*sin(theta/2);
o=zeros(2);

%---------------------------
T1=[1,1;...
1,1];

T2=[exp(-1i*phi),1;...
exp(1i*phi),exp(-1i*phi)];

T3=[exp(1i*phi),1;...
exp(-1i*phi),exp(1i*phi)];

T4=[1,1;...
1,1];

T5=T2';

T6=T3';

%-------------0---------------
Kx1=linspace(0,0,800);
Ky1=linspace(-2,2,800);
Kx2=linspace(0,sqrt(3),692);
Ky2=-sqrt(3)*Kx2+2;
Kx3=linspace(sqrt(3),0,400);
Ky3=Kx3/sqrt(3)-2;
KX0=[Kx1,Kx2,Kx3];
KY0=[Ky1,Ky2,Ky3];

%-------------1---------------
Kx1=linspace(0,0,800);
Ky1=linspace(-4,0,800);
Kx2=linspace(0,sqrt(3),692);
Ky2=-sqrt(3)*Kx2;
Kx3=linspace(sqrt(3),0,400);
Ky3=Kx3/sqrt(3)-4;
KX1=[Kx1,Kx2,Kx3];
KY1=[Ky1,Ky2,Ky3];

%-------------2---------------
Kx1=linspace(sqrt(3),sqrt(3),800);
Ky1=linspace(-1,3,800);
Kx2=linspace(sqrt(3),2*sqrt(3),692);
Ky2=-sqrt(3)*Kx2+6;
Kx3=linspace(2*sqrt(3),sqrt(3),400);
Ky3=Kx3/sqrt(3)-2;
KX2=[Kx1,Kx2,Kx3];
KY2=[Ky1,Ky2,Ky3];

%-------------3----------------
Kx1=linspace(-sqrt(3),-sqrt(3),800);
Ky1=linspace(-1,3,800);
Kx2=linspace(-sqrt(3),0,692);
Ky2=-sqrt(3)*Kx2;
Kx3=linspace(0,-sqrt(3),400);
Ky3=Kx3/sqrt(3);
KX3=[Kx1,Kx2,Kx3];
KY3=[Ky1,Ky2,Ky3];

%-------------4------------------
Kx1=linspace(0,0,800);
Ky1=linspace(-2,2,800);
Kx2=linspace(0,sqrt(3),692);
Ky2=-sqrt(3)*Kx2+2;
Kx3=linspace(sqrt(3),0,400);
Ky3=Kx3/sqrt(3)-2;
KX4=[Kx1,Kx2,Kx3];
KY4=[Ky1,Ky2,Ky3];

%------------5---------------------
Kx1=linspace(sqrt(3),sqrt(3),800);
Ky1=linspace(-5,-1,800);
Kx2=linspace(sqrt(3),2*sqrt(3),692);
Ky2=-sqrt(3)*Kx2+2;
Kx3=linspace(2*sqrt(3),sqrt(3),400);
Ky3=Kx3/sqrt(3)-6;
KX5=[Kx1,Kx2,Kx3];
KY5=[Ky1,Ky2,Ky3];

%------------6---------------------
Kx1=linspace(-sqrt(3),-sqrt(3),800);
Ky1=linspace(-5,-1,800);
Kx2=linspace(-sqrt(3),0,692);
Ky2=-sqrt(3)*Kx2-4;
Kx3=linspace(0,-sqrt(3),400);
Ky3=Kx3/sqrt(3)-4;
KX6=[Kx1,Kx2,Kx3];
KY6=[Ky1,Ky2,Ky3];

%-----------------------------------
for m=1:length(KX1)
kx1=KX1(m);ky1=KY1(m);
kx2=KX2(m);ky2=KY2(m);
kx3=KX3(m);ky3=KY3(m);
kx4=KX4(m);ky4=KY4(m);
kx5=KX5(m);ky5=KY5(m);
kx6=KX6(m);ky6=KY6(m);
kx0=KX0(m);ky0=KY0(m);

%--------------0-----------------
if kx0==0
      if ky0>=0
         thetak0=pi/2;
      else
         thetak0=-pi/2;
      end
elseif ky0>=0
      thetak0=acos(kx0/sqrt(kx0^2+ky0^2));
else
      thetak0=-acos(kx0/sqrt(kx0^2+ky0^2));
end

%--------------1-------------------
   if kx1==0
      thetak1=-pi/2;
else
      thetak1=-acos(kx1/sqrt(kx1^2+ky1^2));
end

%--------------2-------------------
if ky2==0
      thetak2=0;
elseif ky1>0;
      thetak2=acos(kx2/sqrt(kx2^2+ky2^2));
else
      thetak2=-acos(kx2/sqrt(kx2^2+ky2^2));
end

%---------------3--------------------
if ky3==0
      thetak3=pi;
elseif ky3>0
      thetak3=acos(kx3/sqrt(kx3^2+ky3^2));
else
      thetak3=-acos(kx3/sqrt(kx3^2+ky3^2));
end

%---------------4--------------------
if kx4==0
      if ky4>=0
         thetak4=pi/2;
      else
         thetak4=-pi/2;
      end
elseif ky4>=0
      thetak4=acos(kx4/sqrt(kx4^2+ky4^2));
else
      thetak4=-acos(kx4/sqrt(kx4^2+ky4^2));
end

%--------------5---------------------
thetak5=-acos(kx5/sqrt(kx5^2+ky5^2));

%--------------6---------------------
thetak6=-acos(kx6/sqrt(kx6^2+ky6^2));

%--------------k---------------------
k0=sqrt(kx0^2+ky0^2)*ktheta;
k1=sqrt(kx1^2+ky1^2)*ktheta;
k2=sqrt(kx2^2+ky2^2)*ktheta;
k3=sqrt(kx3^2+ky3^2)*ktheta;
k4=sqrt(kx4^2+ky4^2)*ktheta;
k5=sqrt(kx5^2+ky5^2)*ktheta;
k6=sqrt(kx6^2+ky6^2)*ktheta;

%--------------H-----------------
hk0=-v*k0*[0,exp(1i*(thetak0-theta/2));...
         exp(-1i*(thetak0-theta/2)),0];
hk1=-v*k1*[0,exp(1i*(thetak1+theta/2));...
         exp(-1i*(thetak1+theta/2)),0];
hk2=-v*k2*[0,exp(1i*(thetak2+theta/2));...
         exp(-1i*(thetak2+theta/2)),0];
hk3=-v*k3*[0,exp(1i*(thetak3+theta/2));...
         exp(-1i*(thetak3+theta/2)),0];
hk4=-v*k4*[0,exp(1i*(thetak4+theta/2));...
         exp(-1i*(thetak4+theta/2)),0];
hk5=-v*k5*[0,exp(1i*(thetak5+theta/2));...
         exp(-1i*(thetak5+theta/2)),0];
hk6=-v*k6*[0,exp(1i*(thetak6+theta/2));...
         exp(-1i*(thetak6+theta/2)),0];

%---------------------------------
   H= [hk0,w*T1,w*T2,w*T3,w*T4,w*T5,w*T6;...
         w*T1',hk1,o,o,o,o,o;...
         w*T2',o,hk2,o,o,o,o;...
         w*T3',o,o,hk3,o,o,o;...
         w*T4',o,o,o,hk4,o,o;...
         w*T5',o,o,o,o,hk5,o;...
         w*T6',o,o,o,o,o,hk6];
%--------------------------------
E0=eig(H);
E1(m)=E0(1);E2(m)=E0(2);E3(m)=E0(3);E4(m)=E0(4);
E5(m)=E0(5);E6(m)=E0(6);E7(m)=E0(7);E8(m)=E0(8);
E9(m)=E0(9);E10(m)=E0(10);E11(m)=E0(11);E12(m)=E0(12);
E13(m)=E0(13);E14(m)=E0(14);
end
figure
K=[1:1:1892];
plot(K,E1,'-b');hold on
plot(K,E2,'-g');hold on
plot(K,E3,'-r');hold on
plot(K,E4,'-c');hold on
plot(K,E5,'-m');hold on
plot(K,E6,'-y');hold on
plot(K,E7,'-k');hold on
plot(K,E8,'-k');hold on
plot(K,E9,'-y');hold on
plot(K,E10,'-m');hold on
plot(K,E11,'-c');hold on
plot(K,E12,'-r');hold on
plot(K,E13,'-g');hold on
plot(K,E14,'-b');hold on

回复此楼