| 查看: 347 | 回复: 0 | |||
[交流]
关于有转角双层石墨烯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 |
回复此楼» 猜你喜欢
材料调剂
已经有7人回复
-
化学工程085602 305分求调剂
已经有10人回复
-
289求调剂
已经有15人回复
-
291 求调剂
已经有7人回复
-
274求调剂
已经有14人回复
-
北京林业大学硕导招生广告
已经有3人回复
-
309求调剂
已经有5人回复
-
292求调剂
已经有8人回复
-
求调剂
已经有4人回复
-
一志愿 西北大学 总分282 英语一62 求调剂
已经有3人回复
» 本主题相关商家推荐: (我也要在这里推广)
点击这里搜索更多相关资源
5