一维原子链体系态密度的计算

使用大体系实空间哈密顿量结合KPM方法稀疏矩阵，计算态密度

首先需要生成指定的稀疏矩阵

N=10^6;
e = ones (N,1); % e似乎必须写成n*1的形式
H = spdiags ([ e,zeros(N,1), e ],[-1,0,1],N,N);
E = (-1.0 : 0.01 : 1.0) * (2.2);
E_max = max(E) * 1.1; % to be safe
M = 500; % number of moments in the Chebyshev expansion
phi=create_state(N);
tic;
DOS=find_dos(M,E_max,E,H,phi);
toc;
plot(E,DOS,'linewidth', 3)

给定稀疏矩阵计算态密度

function DOS=find_dos(M,E_max,E,H,phi)
% M 是切比雪夫多项式展开项数
% E_max 是计算的最高上限能量
% E 是列向量，为计算的能量范围
% H 是稀疏矩阵形势下的哈密顿矩阵
% phi是随机列向量
% DOS是计算返回的态密度列向量
H=H/E_max;   %标准化处理
E=E/E_max;
phi_1=phi;
phi_2=H*phi_1;
DOS=zeros(length(E),1);  %态密度列向量
CnDOS=zeros(M,1);        % CuDOS(m,1) = "Tr[T_m(H)]" 哈密顿量部分的切比雪夫多项式
CnDOS(1)=1;
CnDOS(2)=phi'*phi_2;
% 哈密顿量部分的各阶切比雪夫多项式
for m=3:M
    temp=phi_2;
    phi_2=2.0*H*phi_2-phi_1;
    phi_1=temp;
    CnDOS(m)=phi'*phi_2; 
end
% M项雅可比核函数
for m=0:M-1 % Apply a damping factor (Jackson kernel)
    a=1.0/(M+1);
    g=(1-m*a)*cos(pi*m*a)+a*sin(pi*m*a)/tan(pi*a);
    CnDOS(m+1)=CnDOS(m+1)*g;
end
% 对采样的能量部分和哈密顿量的切比雪夫多项式求和
for ne=1:length(E) % Perform the Chebyshev Summation
    T1=1.0;
    T2=E(ne);
    DOS(ne)=CnDOS(2)*T2;
    for m=3:M
        T_temp=T2;
        T2=2.0*E(ne)*T2-T1;
        T1=T_temp;
        DOS(ne)=DOS(ne)+CnDOS(m)*T2;
    end
    DOS(ne)=(DOS(ne)*2.0+CnDOS(1))*2.0/pi/sqrt(1.0-E(ne)*E(ne))/E_max; 
end  
DOS=real(DOS);

尝生随机归一化的一维向量

function initial_state = create_state(N)
% 生成所需的N维复数归一化随机列向量
random_phase = rand(N, 1) * 2 * pi;
initial_state = cos(random_phase) + sin(random_phase) * 1i;   
initial_state = initial_state / norm(initial_state); % length = 1