不光滑的态密度

背景

在固体物理或相关物理教材可知，对于如下的非常经典的抛物线型色散关系

\[H = \bm{k}^2 \]

其对应的态密度色散关系为:

一维时, 与 \(E^{-1/2}\) 成正比;
二维时，其是个常数；
三维时, 与 \(E^{1/2}\) 成正比.

态密度公式

(参考下面该篇博士论文附录，先码住有时间我再展开写)

Green函数方法计算

参考以前的笔记balabala

附件

(1) Sunhee_Lee 博士论文中的matlab代码

点击查看代码

function bulk_DOS_test(dim,gv,nk)
% bulk_DOS\_test(dim,gv,nk)
% Comparison numerical result vs. analytical solution of E=k^2
% dim: dimension (1/2/3) of the parabolic band E=k^2.
% gv: Broadening factor in (eV). Set to (0~10e-3)(eV) depending on simulation conditions.
% nk: Number of k-points between the interval k=[0 1].
% Note: setting this value too large might cause memory error,
% especially in 2D and 3D cases.
% Spin degeneracy is set to 2
close all;
if dim==3
    % 3D
    kx=0:1/nk:1;lkx=length(kx);
    ky=0:1/nk:1;lky=length(ky);
    kz=0:1/nk:1;lkz=length(kz);
    [kxx kyy kzz]=meshgrid(kx,ky,kz);
    E=(kxx.^2+kyy.^2+kzz.^2)+0.1;
    Egrid=0.0:1e-3:0.6;
    E=reshape(E,[lkx*lky*lkz 1]);
    kxx=reshape(kxx,[lkx*lky*lkz 1]);
    kyy=reshape(kyy,[lkx*lky*lkz 1]);
    kzz=reshape(kzz,[lkx*lky*lkz 1]);
    DOS=zeros(size(Egrid));
    
    for eidx=1:length(E)
        m_factor=8;
        if abs(kxx(eidx,1)+0)<1e-3 || abs(kxx(eidx,1)-1)<1e-3
            m_factor = m_factor/2;
        end
        if abs(kyy(eidx,1)+0)<1e-3 || abs(kyy(eidx,1)-1)<1e-3
            m_factor = m_factor/2;
        end
        if abs(kzz(eidx,1)+0)<1e-3 || abs(kzz(eidx,1)-1)<1e-3
            m_factor = m_factor/2;
        end
        DOS = DOS + 2/sqrt(2*pi*gv*gv)*exp(-(E(eidx,1)-Egrid).^2 ...
            /(2*gv*gv))*m_factor/(nk*nk*nk);
    end
    figure(1);
    set(gca,'Fontsize',24);
    plot(Egrid,DOS/(8*pi*pi*pi),'b');hold on;
    X=Egrid;
    Y=1/(2*pi*pi)*sqrt(X-0.1);
    plot(X,Y,'r')
    axis([0 0.6 0 max(DOS)/(8*pi*pi*pi)*1.2])
    xlabel('Energy(eV)');
    ylabel('DOS(#/eV/m^{3})');
end

if dim==2
    %% 2D
    kx=0:1/nk:1;
    ky=kx;
    [kxx kyy ]=meshgrid(kx,ky);
    lkx=length(kx);
    lky=length(ky);
    E=(kxx.^2+kyy.^2)+0.1;
    Egrid=0.0:1e-3:0.6;
    E=reshape(E,[lkx*lky 1]);
    kxx=reshape(kxx,[lkx*lky 1]);
    kyy=reshape(kyy,[lkx*lky 1]);
    DOS=zeros(size(Egrid));
    for eidx=1:length(E)
        m_factor=4;
        if abs(kxx(eidx,1)+0)<1e-3 || abs(kxx(eidx,1)-1)<1e-3
            m_factor = m_factor/2;
        end
        if abs(kyy(eidx,1)+0)<1e-3 || abs(kyy(eidx,1)-1)<1e-3
            m_factor = m_factor/2;
        end
        DOS = DOS + 2/sqrt(2*pi*gv*gv)*exp(-(E(eidx,1)-Egrid).^2/(2*gv*gv))*m_factor/(nk*nk);
    end
    figure(2);
    set(gca,'Fontsize',24);
    plot(Egrid,DOS/(4*pi*pi),'b');hold on;
    X=Egrid;
    for idx=1:length(X)
        if X(idx)<0.1
            Y(idx)=0;
        else
            Y(idx)=1/(2*pi);
        end
    end
    plot(X,Y,'r')
    axis([0 0.6 0 max(DOS)/(2*pi*2*pi)*1.2])
    xlabel('Energy(eV)');
    ylabel('DOS(#/eV/m^{2})');
end
if dim==1
    %% 1D
    kx=0:1/nk:1;
    Egrid=-0.1:1e-3:0.6;
    kxx=kx';
    E=(kxx.^2)+0.1;
    DOS=zeros(size(Egrid));
    for eidx=1:length(E)
        m_factor=2;
        if abs(kxx(eidx,1)+0)<5e-6 || abs(kxx(eidx,1)-1)<5e-6
            m_factor = m_factor/2;
        end
        DOS = DOS + 2/sqrt(2*pi*gv*gv)*exp(-(E(eidx,1)-Egrid).^2 ...
            /(2*gv*gv))*m_factor/nk;
    end
    figure(3);
    set(gca,'Fontsize',24);
    plot(Egrid,DOS/(2*pi),'b');hold on;
    
    X=1e-8:1e-3:1.0;
    Y=1/(1*pi)./sqrt(X-0.1);
    plot(X,Y,'r')
    axis([0 0.6 0 max(DOS)/(2*pi)*1.2])
    xlabel('Energy(eV)');
    ylabel('DOS(#/eV/m)');
end

(2) Green函数方法

点击查看代码 dos_1D.m

clear all;
tic

eta = 5e-3;
E = -0.1:1e-3:0.6;
N_E = length(E);

N_kx = 200;
kx = linspace(0,1,N_kx);
dkx = kx(2) - kx(1);

DOS = zeros(N_E,1);
for i_E = 1:N_E
    for i_kx = 1:N_kx
        Ham =  kx(i_kx)^2 + 0.1;
        Gr = inv(E(i_E)+1j*eta - Ham);
        DOS(i_E) = DOS(i_E) - imag(Gr)/pi*dkx/(2*pi);
    end
end
figure;
plot(E,DOS);
xlim([0,0.6])
toc

点击查看代码 dos_2D.m

clear all;
tic

eta = 1e-4;

N_E = 500;
E = linspace(-1,1,N_E);

N_kx = 300;
kx = linspace(-pi,pi,N_kx);
dkx = kx(2) - kx(1);
N_ky = 300;
ky = linspace(-pi,pi,N_ky);
dky = ky(2) - ky(1);


DOS = zeros(N_E,1);
for i_E = 1:N_E
    for i_kx = 1:N_kx
        for i_ky = 1:N_ky
            Ham =  kx(i_kx)^2 + ky(i_ky)^2;
            Gr = inv(E(i_E)+1j*eta - Ham);
            DOS(i_E) = DOS(i_E) - imag(Gr)/pi*dkx*dky/(4*pi^2);
        end
    end
end
figure;
plot(E,DOS);
toc

点击查看代码 dos_3D.m

clear all;
tic

eta = 1e-3;

N_E = 500;
E = linspace(-1,1,N_E);

N_kx = 50;
kx = linspace(-pi,pi,N_kx);
dkx = kx(2) - kx(1);
N_ky = 50;
ky = linspace(-pi,pi,N_ky);
dky = ky(2) - ky(1);
N_kz = 50;
kz = linspace(-pi,pi,N_kz);
dkz = kz(2) - kz(1);

DOS = zeros(N_E,1);
for i_E = 1:N_E
    for i_kx = 1:N_kx
        for i_ky = 1:N_ky
            for i_kz = 1: N_kz
                Ham =  sin(kx(i_kx))^2 + sin(ky(i_ky))^2 + sin(kz(i_kz))^2;
                %Ham = v*(kx(i_kx)*sigma_x + ky(i_ky)*sigma_y) + delta/2*sigma_z + Ef*sigma_0 + m*(kx(i_kx)^2 + ky(i_ky)^2)*sigma_0;
                Gr = inv(E(i_E)+1j*eta - Ham);
                DOS(i_E) = DOS(i_E) + sum(diag(-imag(Gr))/pi)*dkx*dky*dkz/(8*pi^3);
            end
        end
    end
end
figure;
plot(E,DOS);

toc

参考文献

[1] Sunhee Lee PhD thesis. Purdue University. https://engineering.purdue.edu/gekcogrp/publications/theses/PhD_11_2011_Sunhee_Lee_PhD_Thesis_main.pdf