计算简单体系的电导

求解简单体系电导MATLAB程序如下:

clear all; 
tic
W = 30;
L = 100; %% run 4h
t = 1;
eta = 1e-10;
Ham = zeros(W*L,W*L);
SigmaL = zeros(W*L,W*L);
SigmaR = zeros(W*L,W*L);
GammaL = zeros(W*L,W*L);
GammaR = zeros(W*L,W*L);
H00 = zeros(W,W);
H01 = zeros(W,W);
for i = 1:W
    H00(i,i) = 4*t;
    H01(i,i) = t;
    if i < W
        H00(i,i+1) = t;
        H00(i+1,i) = t;
    end
end

for j = 1:L
    Ham(j*W-W+1:j*W,j*W-W+1:j*W) = H00;
    if j < L
        Ham(j*W-W+1:j*W,j*W+1:j*W+W) = H01;
        Ham(j*W+1:j*W+W,j*W-W+1:j*W) = H01';
    end
end

N_omega = 300;
omega = linspace(0,1,N_omega);
Cond = zeros(N_omega,1);
for i = 1:N_omega
    Gs = surface_green_function(H00,H01,omega(i));
    SigmaL(1:W,1:W) = H01*Gs*H01';
    SigmaR(W*L-W+1:W*L, W*L-W+1:W*L) = H01*Gs*H01';
    GammaL = 1j*(SigmaL-SigmaL');
    GammaR = 1j*(SigmaR-SigmaR');
    Gr = inv((omega(i) + 1j*eta)*eye(W*L) - Ham - SigmaL - SigmaR); 
    Ga = Gr';
    Cond(i) = trace(GammaL*Gr*GammaR*Ga);   
end
figure;
plot(omega,real(Cond));
xlabel('Energy (t)');
ylabel('Conducatance (e^2/h)')
toc

　　

 

Dyson 迭代优化代码:

clear all;
tic
W = 30;
L = 100; %% 17s
t = 1;
eta = 1e-10;
SigmaL = zeros(W,W);SigmaR = zeros(W,W);
GammaL = zeros(W,W);GammaR = zeros(W,W);
H00 = zeros(W,W);H01 = zeros(W,W);
for i = 1:W
    H00(i,i) = 4*t;
    H01(i,i) = t;
    if i < W
        H00(i,i+1) = t;
        H00(i+1,i) = t;
    end
end

N_omega = 300;
omega = linspace(0,1,N_omega);
Cond = zeros(N_omega,1);
for i = 1:N_omega
    Gs = surface_green_function(H00,H01,omega(i));
    SigmaL = H01*Gs*H01';SigmaR = H01*Gs*H01';
    G11 = inv((omega(i) + 1j*eta)*eye(W)- H00 - H01*Gs*H01');
    Gjj_j = G11;G1j_j = G11;
    for j = 2:L-1 
       %G_{N+1,N+1}^{N+1} = inv(omega - H_{N+1,N+1} - H01'*G_{N,N}^{N}*H01)
       Gj1j1_j1 = inv((omega(i) + 1j*eta)*eye(W)- H00 - H01'*Gjj_j*H01);
       %G_{i,N+1}^{N+1} = G_{i,N}^{N}*V_{N}*G_{N+1,N+1}^{N+1}.
       G1j1_j1 = G1j_j*H01*Gj1j1_j1;
       Gjj_j = Gj1j1_j1;G1j_j = G1j1_j1;
    end
    GLL_L = inv((omega(i) + 1j*eta)*eye(W) - H00 - H01'*Gjj_j*H01 - H01*Gs*H01');
    G1L_L = G1j_j*H01*GLL_L;
    GammaL = 1j*(SigmaL-SigmaL');GammaR = 1j*(SigmaR-SigmaR');
    Cond(i) = trace(GammaL*G1L_L*GammaR*G1L_L');
end
plot(omega,real(Cond));
xlabel('Energy [t]');
ylabel('Conducatance [e^2/h]')
toc

电导结果为:

 

 

surface_green_function.m 程序如下:

function Gs = surface_green_function(H00,H01,omega)
    eta = 1e-10;
    [m, n] = size(H00); 
    % Initial
    alphai_1 = H01*((omega+1j*eta)*eye(m)-H00)^(-1)*H01;
    betai_1 = H01'*((omega+1j*eta)*eye(m)-H00)^(-1)*H01';
    ei_1s = H00 + H01*((omega+1j*eta)*eye(m)-H00)^(-1)*H01';
    ei_1 = H00 + H01*((omega+1j*eta)*eye(m)-H00)^(-1)*H01' + H01'*((omega+1j*eta)*eye(m)-H00)^(-1)*H01;

    for i = 1:100
        alphai = alphai_1*((omega+1j*eta)*eye(m)-ei_1)^(-1)*alphai_1;
        betai = betai_1*((omega+1j*eta)*eye(m)-ei_1)^(-1)*betai_1;
        ei = ei_1 + alphai_1*((omega+1j*eta)*eye(m)-ei_1)^(-1)*betai_1 + betai_1*((omega+1j*eta)*eye(m)-ei_1)^(-1)*alphai_1;
        eis = ei_1s + alphai_1 *((omega+1j*eta)*eye(m)-ei_1)^(-1)*betai_1;
        % exit the loop
        if sum(sum(abs(alphai))) <= 1e-8
            break
        end
        % iterative
        alphai_1 = alphai;
        betai_1 = betai;
        ei_1s = eis;
        ei_1 = ei;  
    end
    Gs = ((omega+1j*eta)*eye(m)-eis)^(-1);
end