下面将为你提供一个用 MATLAB 计算并绘制 zigzag 石墨烯 ribbons 能带和态密度的代码示例。

下面将为你提供一个用 MATLAB 计算并绘制 zigzag 石墨烯 ribbons 能带和态密度的代码示例。

思路概述

1. 构建哈密顿量：根据 zigzag 石墨烯 ribbons 的结构特点，考虑最近邻跳跃相互作用构建哈密顿量矩阵。
2. 计算能带：在一定的动量空间内对哈密顿量进行对角化，得到本征值，即能带。
3. 计算态密度：通过对能带进行展宽处理来计算态密度。

4. 代码解释

   1. 参数设置：定义了条带长度 L、最近邻跳跃能 t、动量点数 Nk、能量范围 Emin 和 Emax、能量点数 NE 以及展宽因子 eta。
   2. 动量空间和能带计算：在动量空间 kx 内循环，对每个动量点构建哈密顿量矩阵 H，并进行对角化得到本征值，即能带 bands。
   3. 态密度计算：在能量网格 E_grid 上，通过对能带进行展宽处理（洛伦兹展宽）计算态密度 DOS。
   4. 绘图：使用 subplot 函数将能带图和态密度图绘制在同一窗口中。
   
   
   你可以根据需要调整参数，如条带长度 L、动量点数 Nk 等，以获得不同的结果。

% 参数设置

L = 20; % 条带长度

t = 1; % 最近邻跳跃能

Nk = 200; % 动量点数

Emin = -3*t; % 能量最小值

Emax = 3*t; % 能量最大值

NE = 200; % 能量点数

eta = 0.01; % 展宽因子

% 动量空间

kx = linspace(-pi, pi, Nk);

bands = zeros(2*L, Nk);

% 构建哈密顿量并计算能带

for ik = 1:Nk

    H = zeros(2*L);

    for n = 1:L

        % 本元胞内 A - B 子晶格跳跃

        if n < L

            H(2*n - 1, 2*n + 1) = -t;

            H(2*n + 1, 2*n - 1) = -t;

            H(2*n, 2*n + 2) = -t;

            H(2*n + 2, 2*n) = -t;

        end

        % 相邻元胞间 A - B 子晶格跳跃

        H(2*n - 1, 2*n) = -t * exp(1i * kx(ik));

        H(2*n, 2*n - 1) = -t * exp(-1i * kx(ik));

    end

    % 对角化哈密顿量

    [~, E] = eig(H);

    E = diag(E);

    [E_sorted, ~] = sort(E);

    bands(:, ik) = E_sorted;

end

% 计算态密度

E_grid = linspace(Emin, Emax, NE);

DOS = zeros(NE, 1);

for iE = 1:NE

    for ik = 1:Nk

        for n = 1:2*L

            DOS(iE) = DOS(iE) + eta / ((E_grid(iE) - bands(n, ik))^2 + eta^2);

        end

    end

end

DOS = DOS / (Nk * pi);

% 绘制能带图

figure;

subplot(1, 2, 1);

for n = 1:2*L

    plot(kx, bands(n, :));

    hold on;

end

xlabel('$k_x$', 'Interpreter', 'latex');

ylabel('Energy', 'Interpreter', 'latex');

title('Band Structure of Zigzag Graphene Ribbons', 'Interpreter', 'latex');

grid on;

% 绘制态密度图

subplot(1, 2, 2);

plot(DOS, E_grid);

xlabel('Density of States', 'Interpreter', 'latex');

ylabel('Energy', 'Interpreter', 'latex');

title('Density of States of Zigzag Graphene Ribbons', 'Interpreter', 'latex');

grid on;