%Graphene
clear all
tic
%Constants (all MKS, except energy which is in eV)
hbar=1.06e-34;q=1.6e-19;qh=q/hbar;a=1e-9;
 
%inputs
t0=-2.5;D=1e-12;ctr=0;zplus=i*1e-6;
NW=25;NL=1;L=zeros(NL);R=L;L(1,1)=1;R(NL,NL)=1;
config=1;%1 for armchair, 2 for zigzag edge
 
%Hamiltonian
al=t0*[0 1 0 0;1 0 1 0;0 1 0 1;0 0 1 0];
if config==1
bL=t0*[0 0 0 0;0 0 0 0;0 0 0 0;1 0 0 0];
bW=t0*[0 0 0 0;1 0 0 0;0 0 0 1;0 0 0 0];end
if config==2
bW=t0*[0 0 0 0;0 0 0 0;0 0 0 0;1 0 0 0];
bL=t0*[0 0 0 0;1 0 0 0;0 0 0 1;0 0 0 0];end
 
n=4;% al=4;bW=-1;bL=-1;n=1;
alpha=kron(eye(NW),al)+kron(diag(ones(1,NW-1),+1),bW)+kron(diag(ones(1,NW-1),-1),bW');
%  alpha=alpha+kron(diag(ones(1,1),1-NW),bW)+kron(diag(ones(1,1),NW-1),bW');% for CNT's
 
sigB=zeros(NW*NL*n);siginB=zeros(NW*NL*n);
ii=0;for EE=-0.6:+0.02:0.6
 for BB=[0:.1:0]
ii=ii+1;
ig0=(EE+zplus)*eye(NW*n)-alpha;
if ii==1
gs1=inv(ig0);gs2=inv(ig0);end
 
beta=kron(diag(exp(i*qh*BB*a*a*[1:1:NW])),bL);
%beta=kron(eye(NW),bL);
H=kron(eye(NL),alpha);if NL>1
H=H+kron(diag(ones(1,NL-1),+1),beta)+kron(diag(ones(1,NL-1),-1),beta');end
 
    change=1;
    while change>4e-5
        Gs=inv(ig0-beta'*gs1*beta);
        change=sum(sum(abs(Gs-gs1)))/(sum(sum(abs(gs1)+abs(Gs))));
        gs1=0.5*Gs+0.5*gs1;
    end
    sig1=beta'*gs1*beta;sig1=kron(L,sig1);gam1=i*(sig1-sig1');
    change=1;
    while change>4e-5
        Gs=inv(ig0-beta*gs2*beta');
        change=sum(sum(abs(Gs-gs2)))/(sum(sum(abs(gs2)+abs(Gs))));
        gs2=0.5*Gs+0.5*gs2;
    end
    sig2=beta*gs2*beta';sig2=kron(R,sig2);gam2=i*(sig2-sig2');
 
    change=1;
    while change>1e-5
        G=inv((EE*eye(NW*NL*n))-H-sig1-sig2-sigB);sigBnew=D*G;
        change=ctr*(sum(sum(abs(sigBnew-sigB))))/(sum(sum(abs(sigBnew+sigB))));
        sigB=sigB+0.5*ctr*(sigBnew-sigB);
    end
    change=1;
    while change>1e-5
        Gn=G*(gam1+siginB)*G';siginBnew=D*Gn;
        change=ctr*(sum(sum(abs(siginBnew-siginB))))/(sum(sum(abs(siginBnew+siginB))));
        siginB=siginB+0.5*ctr*(siginBnew-siginB);
    end
    GN(:,ii)=real(diag(Gn));A(:,ii)=real(diag(i*(G-G')));
    T(ii)=real(trace(gam2*Gn));B(ii)=BB;E(ii)=EE;EE
end
end
toc
hold on
%%
%imagesc(A)
%colormap(gray)
figure(1)
h=plot(T,E,'b');
set(h,'linewidth',[3.0])
h=plot(T,E,'b+');
set(gca,'Fontsize',[24]);
grid on