NEGF程序计算态密度，电子密度，透射率

计算电子态密度主程序，体系的坐标和哈密顿量在之前的博客里面有

需要用Pade展开费米函数

ny=48; 
nx=3;   
% [x,y]=zigzagP(N,number);
[x,y]=zigzagP(ny,nx);                             %生成zigzag石墨烯的坐标
% plot(x,y,'.','MarkerSize',20); 
t1=-2.7;
H0=Hamiltonian_NN_graphene(x,y,t1,1);             %最近邻hooping
e=0.01;
m=length(H0);
Hab=Hab(m,e);
H0=H0+Hab;
hL=H0(1:ny,1:ny);                                 %取出最左边一列的哈密顿量，用于半无穷表面格林函数的计算
hR=H0(ny*(nx-1)+1:ny*nx,ny*(nx-1)+1:ny*nx);       %取出最右边的一列的哈密顿量
hLD=H0(1:ny,ny+1:2*ny);                           %取出最左与右边一列的耦合矩阵
hDL=hLD';
hRD=H0((nx-1)*ny+1:nx*ny,ny+1:2*ny);              %取出最右边与左边的耦合矩阵
hDR=hRD';
HD=H0(ny+1:ny*(nx-1),ny+1:(nx-1)*ny);             %取出device区域的哈密顿量
HLD=H0(1:ny,ny+1:ny*(nx-1));                      %左边与device区域的耦合矩阵
HRD=H0(ny*(nx-1)+1:ny*nx,ny+1:(nx-1)*ny);         %右边与device区域的耦合矩阵
HDL=HLD';
HDR=HRD';
density=zeros(length(HD));                        %初始化密度矩阵
%导入pade系数
file1=fopen('Rp.txt','r');                        %Pade展开费米函数的系数
file2=fopen('zp.txt','r');                        %Pade展开费米函数的留数
Rp=fscanf(file1,'%f');
zp=fscanf(file2,'%f');
fclose(file1);
fclose(file2);
kT=0.038,u=0; %温度与化学势
%准备工作完成，开始计算                            
number =length(zp);                               %pade系数项数
for j=1:20 
gL=NEGFSe(zp(j)*i*kT+u,hL,hRD,100);               %计算左表面格林函数
gR=NEGFSe(zp(j)*i*kT+u,hR,hLD,100);               %计算右表面格林函数
sigmaL=HDL*gL*HLD;                                %计算左自能
sigmaR=HDR*gR*HRD;                                %计算右自能
Gd=((zp(j)*i*kT+u)*eye(length(HD))-HD-sigmaL-sigmaR)^(-1);%device区域格林函数
density=density+2*pi*i*Rp(j)*Gd;%密度矩阵叠加
end
density=(1/2)*eye(length(HD))+imag(density)/(-pi);%得到密度矩阵
Density=diag(density);%取出对角部分电子密度

其中计算半无穷表面格林函数的程序

function g=NEGFSe(E,Ha,Hb,N); 
%给定一个E值，计算对应的半无穷格林函数g
delta=1.0*10^(-3);
e=E+i*delta;
A(:,:,1)=Hb(:,:); B(:,:,1)=Ha(:,:);
C(:,:,1)=Hb(:,:)'; D(:,:,1)=Ha(:,:);
for n=1:N
    A(:,:,n+1)=A(:,:,n)*inv(e*eye(length(Ha))-B(:,:,n))*A(:,:,n);
    B(:,:,n+1)=B(:,:,n)+C(:,:,n)*inv(e*eye(length(Ha))-B(:,:,n))*A(:,:,n) + A(:,:,n)*inv((e*eye(length(Ha))-B(:,:,n)))*C(:,:,n);
    C(:,:,n+1)=C(:,:,n)*inv(e*eye(length(Ha))-B(:,:,n))*C(:,:,n);
    D(:,:,n+1)=D(:,:,n)+A(:,:,n)*inv(e*eye(length(Ha))-B(:,:,n))*C(:,:,n);
end 
g=inv(e*eye(length(Ha))-D(:,:,n+1));

加入AB格子后有电荷密度波CDW

下面的代码是计算实空间的局域态密度LDOS的，局域态密度为格林函数的负虚部/pi

ny=48; 
nx=3;   
[x,y]=zigzagP(ny,nx);                             %生成zigzag石墨烯的坐标
% plot(x,y,'.','MarkerSize',20); 
t1=-2.7;
e=0.1;
H0=Hamiltonian_NN_graphene(x,y,t1, 1);            %最近邻hooping
m=length(H0);
Hab=Hab(m,e);
H0=H0+Hab;
hL=H0(1:ny,1:ny);                                 %取出最左边一列的哈密顿量，用于半无穷表面格林函数的计算
hR=H0(ny*(nx-1)+1:ny*nx,ny*(nx-1)+1:ny*nx);       %取出最右边的一列的哈密顿量
hLD=H0(1:ny,ny+1:2*ny);                           %取出最左与右边一列的耦合矩阵
hDL=hLD';
hRD=H0((nx-1)*ny+1:nx*ny,ny+1:2*ny);              %取出最右边与左边的耦合矩阵
hDR=hRD';
HD=H0(ny+1:ny*(nx-1),ny+1:(nx-1)*ny);             %取出device区域的哈密顿量
HLD=H0(1:ny,ny+1:ny*(nx-1));                      %左边与device区域的耦合矩阵
HRD=H0(ny*(nx-1)+1:ny*nx,ny+1:(nx-1)*ny);         %右边与device区域的耦合矩阵
HDL=HLD';
HDR=HRD';
delta = 0.0001;                                     %计算用的小虚部
E=-0.05:0.01:0.05;
m=length(HD);
n=length(E);
Density=zeros(m,n);
for i=1:n
     gL=NEGFSe(E(i)+1i*delta,hL,hRD,100);              %左表面格林函数
     gR=NEGFSe(E(i)+1i*delta,hR,hLD,100);              %右表面格林函数
     sigmaL=HDL*gL*HLD;                                %左自能
    sigmaR=HDR*gR*HRD;                                %右自能
    Gd=(E(i)*eye(length(HD))-HD-sigmaL-sigmaR)^(-1);  %device区域格林函数
    Density(:,i)=-imag(diag(Gd))/pi;
end
Zs=0;
for i=1:n
    Zs=Density(:,i)+Zs;
end
plot(Zs)

最后是计算体系透射率的

ny=48; 
nx=3;   
[x,y]=zigzagP(ny,nx);                             %生成zigzag石墨烯的坐标
% plot(x,y,'.','MarkerSize',20); 
t1=-2.7;
e=0.1;
H0=Hamiltonian_NN_graphene(x,y,t1, 1);            %最近邻hooping
m=length(H0);
Hab=Hab(m,e);
H0=H0+Hab;
hL=H0(1:ny,1:ny);                                 %取出最左边一列的哈密顿量，用于半无穷表面格林函数的计算
hR=H0(ny*(nx-1)+1:ny*nx,ny*(nx-1)+1:ny*nx);       %取出最右边的一列的哈密顿量
hLD=H0(1:ny,ny+1:2*ny);                           %取出最左与右边一列的耦合矩阵
hDL=hLD';
hRD=H0((nx-1)*ny+1:nx*ny,ny+1:2*ny);              %取出最右边与左边的耦合矩阵
hDR=hRD';
HD=H0(ny+1:ny*(nx-1),ny+1:(nx-1)*ny);             %取出device区域的哈密顿量
HLD=H0(1:ny,ny+1:ny*(nx-1));                      %左边与device区域的耦合矩阵
HRD=H0(ny*(nx-1)+1:ny*nx,ny+1:(nx-1)*ny);         %右边与device区域的耦合矩阵
HDL=HLD';
HDR=HRD';
delta = 0.0001;                                     %计算用的小虚部
E=-1:0.01:1;                                    %计算能量范围
n=length(E);
T=zeros(n,1);
for i=1:n
      gL=NEGFSe(E(i)+1i*delta,hL,hRD,100);              %左表面格林函数
      gR=NEGFSe(E(i)+1i*delta,hR,hLD,100);              %右表面格林函数
      sigmaL=HDL*gL*HLD;                                %左自能
      sigmaR=HDR*gR*HRD;                                %右自能
      Gd=(E(i)*eye(length(HD))-HD-sigmaL-sigmaR)^(-1);  %device区域格林函数
      GammaL=-2*imag(sigmaL);                           %左关联函数
      GammaR=-2*imag(sigmaR);                           %右关联函数
      T(i,1)=trace(GammaL*Gd*GammaR*conj(Gd));
end
plot(real(T),E)