close all
clear all
% Parameters and constants for wire
N=10; % Number of atoms long the wire
M=10;% Number of atoms long the right-lead
L=10;% Number of atoms long the left-lead
tw=.75;
mu1=.45;
mu2=.25;
kT=0.025;
hbar=1.06e-34;
q=1.6e-19;
% Nearest-neighbor hopping energy
t_LL=1; % hopping energy for left-lead
t_RL=1; % hopping energy for right-lead
t_LW=.25;% hopping energy for wire-leftlead
t_WR=.25; % hopping energy for wire-rightlead
Enn=0; % on-site energy
eta=1i*.001;% Small imaginary number
% Hamiltoni
Hw=zeros (N,N);
Hw=(tw*diag (ones (1,N-1),1))+(tw*diag (ones (1,N-1),-1))+(Enn*diag (ones (1,N),0));
H_RL=zeros (M,M);
H_RL=(t_RL*diag (ones (1,M-1),1))+(t_RL*diag (ones (1,M-1),-1));
H_LL=zeros (L,L);
H_LL=(t_LL*diag (ones (1,L-1),1))+(t_LL*diag (ones (1,L-1),-1));
H_LW=zeros (L,N);
H_LW (1,1)=t_LW;
H_WR=zeros (N,M);
H_WR (N,N)=t_WR;
% Find DOS from the trace of imaginary Green function
%TT=[];
%E=[];
%II=[];
%for V=0:.1:.6
syms E
z=(E+eta)*eye (N);
g_R=inv ((E+eta)*eye (M)-H_RL);
g_L=inv ((E+eta)*eye (L)-H_LL);
sigma_R=H_WR*g_R*H_WR';
sigma_L=H_LW'*g_L*H_LW;
gr=inv (z-Hw-sigma_R-sigma_L);
f1=1./(1+exp((E-mu1)./kT))
f2=1./(1+exp((E-mu2)./kT))
DO1=i*(sigma_L-sigma_L');
DO2=i*(sigma_R-sigma_R');
A2=(gr*DO2*gr');
tic
T=trace(DO1*A2)
I=int((q./hbar)*T*(f1-f2),E,-100,100)
toc
%EE=[EE,E];
%%II=[II,I];
%end
%%plot(EE,II);
%xlabel ('Energy (ev)')
%%%ylabel ('Cuurent')