MATLAB手写紧束缚模型计算石墨烯能带结构-北京尧图网络

简介本资源是一套面向凝聚态物理与计算材料学初学者及研究者的石墨烯能带结构仿真MATLAB代码集聚焦于理解石墨烯电子性质的核心——线性色散关系与Dirac点特征。代码基于紧束缚模型完整覆盖晶格建模、布里渊区定义、薛定谔方程数值求解及能带图可视化全流程特别强化K点附近锥形能带分析助力掌握量子输运基础与二维材料模拟方法。压缩包共6个.m文件如graphene_NN.m、armchair_NNN.m等分别对应不同近邻阶数与边界条件的建模实现总大小仅2KB轻量易读、便于调试与拓展。已有1325人学习下载适合高校物理/材料专业本科生课程实践、研究生课题入门以及科研人员快速复现经典能带结果并开展参数敏感性分析。1. 用 MATLAB 快速计算并可视化石墨烯单层能带结构不是调用现成工具箱而是从紧束缚哈密顿量出发推导出可复现、可修改、可验证的完整代码链如果你在凝聚态物理、新型二维材料或计算材料学方向做研究大概率遇到过这个场景导师或审稿人问“石墨烯的狄拉克锥是怎么算出来的能不能把哈密顿量构造过程和布里渊区采样逻辑展示出来”——此时拿一个黑盒软件截图或调用某封装函数的脚本往往不够有说服力。本文聚焦标题中的核心动作“matlab 石墨烯代码_石墨烯能带”即在 MATLAB 环境中不依赖第三方能带计算软件如 VASP、Quantum ESPRESSO也不调用预编译的能带绘图工具箱而是从第一性原理紧束缚模型出发手写哈密顿量矩阵、定义六角布里渊区路径、对角化求解本征值、并绘制高对称点连线的能带色散曲线。整套流程完全基于 MATLAB 基础语法R2018a 及以上兼容所有变量命名直指物理含义如t -2.8; % eV, nearest-neighbor hopping矩阵维度与晶格自由度严格对应输出结果可与文献中经典图示如 Slonczewski-Weiss-McClure 模型下的双狄拉克锥逐点比对。适合刚接触能带理论的研究生快速建立物理图像也适合需要嵌入自定义参数如应变、外场、多层耦合的进阶用户二次开发。2. 构建石墨烯紧束缚哈密顿量从晶格矢量、原子位置到 k 空间哈密顿矩阵的完整推导石墨烯能带的核心在于其 AB 二原子基元与六角晶格对称性。MATLAB 中实现的关键不是“画图”而是准确表达哈密顿量随波矢 k 的依赖关系。这要求我们显式写出每个 k 对应的 2×2 哈密顿矩阵 H(k)而非一次性生成大矩阵再切片。以下步骤严格遵循固体物理惯例每行代码均可追溯至《Solid State Physics》Ashcroft Mermin第 10 章的数学表述。2.1 定义晶格与基元用向量运算生成真实空间坐标石墨烯晶格常数 a 2.46 Å但 MATLAB 计算中单位可归一化。重点在于正确构建两个子晶格 A 和 B 的相对位移% 晶格参数归一化单位a1 a 1; % 晶格常数 a1 a * [sqrt(3)/2, 1/2]; % 基矢1沿x轴偏转30度 a2 a * [sqrt(3)/2, -1/2];% 基矢2沿x轴偏转-30度 % 原子位置相对于原胞原点 delta1 [0, 0]; % A原子在原胞内位置 delta2 a/3 * [1, 0]; % B原子相对A的位移注意此处采用标准AB基元定义 % 验证最近邻距离应为 a/sqrt(3) ≈ 0.577a d_AB norm(delta2); % 输出应为 0.333...需修正——实际标准位移为 [1/3, 1/3] 在 a1,a2 坐标系下注意上述delta2初始设定有误这是新手最常踩的坑。石墨烯中 A、B 原子在原胞内的分数坐标应为 (0,0) 和 (1/3,1/3)以 a1,a2 为基。正确写法是delta2_frac [1/3, 1/3]; delta2 delta2_frac(1)*a1 delta2_frac(2)*a2; % 转换为笛卡尔坐标2.2 紧束缚哈密顿量 H(k) 的物理构成与 MATLAB 实现石墨烯最简紧束缚模型只保留最近邻跃迁 t ≈ −2.8 eV。H(k) 是 2×2 矩阵形式为H(k) [ 0 f(k) ] [ f*(k) 0 ]其中f(k) t * (exp(i*k·δ₁) exp(i*k·δ₂) exp(i*k·δ₃))δ₁,δ₂,δ₃ 是 A→B 的三个最近邻向量。在 MATLAB 中必须显式计算这三个向量% 三个最近邻向量从A指向B单位晶格常数a d1 a/3 * [1, 1]; % 注意此为笛卡尔坐标下的标准值需与a1,a2一致 d2 a/3 * [-1, 0]; d3 a/3 * [0, -1]; % 更严谨做法由基矢导出 d1_cart -a1/3 a2/3; % 标准推导结果 d2_cart a1/3 a2/3; d3_cart -a2/3; % 跃迁积分 teV t -2.8; % 对给定k [kx, ky]计算f(k) function f_val f_k(k, d1, d2, d3, t) phase1 exp(1i * dot(k, d1)); phase2 exp(1i * dot(k, d2)); phase3 exp(1i * dot(k, d3)); f_val t * (phase1 phase2 phase3); end2.3 布里渊区边界与高对称点坐标的 MATLAB 解析生成能带图必须沿特定路径Γ→M→K→Γ绘制。这些点的坐标需由倒格子基矢 b1,b2 精确计算不能手动输入近似值% 倒格子基矢由正格子基矢a1,a2通过公式 b_i 2π * ε_ij * a_j / (a1×a2) 得到 area cross([a1,0], [a2,0]); % 2D叉积模长 |a1x*a2y - a1y*a2x| b1 2*pi/area * [a2(2), -a2(1)]; % MATLAB中cross不直接支持2D改用分量式 b2 2*pi/area * [-a1(2), a1(1)]; % 高对称点分数坐标 → 笛卡尔坐标 Gamma [0, 0]; M [0.5, 0]; % M点b1/2 K [1/3, 1/3]; % K点b1/3 b2/3 % 转换为笛卡尔坐标 Gamma_cart Gamma(1)*b1 Gamma(2)*b2; M_cart M(1)*b1 M(2)*b2; K_cart K(1)*b1 K(2)*b2; % 验证K点是否满足 |K| 4π/(3a) —— 这是狄拉克点存在的几何前提 norm_K norm(K_cart); expected_K 4*pi/(3*a); fprintf(K点模长: %.6f, 理论值: %.6f\n, norm_K, expected_K); % 应接近2.3.1 生成 k 路径的线性插值序列能带计算本质是沿路径对每个 k 点求解 H(k) 的本征值。路径需足够密以分辨狄拉克锥尖锐性% 定义路径段Γ→M→K→Γ每段采样点数 N_GM 50; N_MK 60; N_KG 50; k_path []; % Γ→M 段 k_GM linspace(Gamma_cart, M_cart, N_GM); k_path [k_path; k_GM(1:end-1,:)]; % 去掉重复端点 % M→K 段 k_MK linspace(M_cart, K_cart, N_MK); k_path [k_path; k_MK(1:end-1,:)]; % K→Γ 段 k_KG linspace(K_cart, Gamma_cart, N_KG); k_path [k_path; k_KG]; % 总点数 N_total size(k_path, 1); % 例如 158 点2.4 对每个 k 点构造并求解哈密顿量避免循环低效用向量化预分配对 158 个 k 点逐个调用eig()是可行的但 MATLAB 中更高效的做法是预分配本征值数组并用循环确保逻辑清晰E zeros(N_total, 2); % 存储每个k点的2个本征值eV for i 1:N_total k k_path(i, :); f_val f_k(k, d1_cart, d2_cart, d3_cart, t); H [0, f_val; conj(f_val), 0]; % 2x2 Hermitian matrix eigvals eig(H); % 返回列向量 E(i, :) sort(real(eigvals)); % 排序确保价带/导带对应 end % 验证在K点本征值应趋近于0狄拉克点 fprintf(K点附近本征值取路径中第%d点: [%.6f, %.6f] eV\n, ... round(0.6*N_total), E(round(0.6*N_total),1), E(round(0.6*N_total),2));提示sort(real(eigvals))是关键。eig()返回顺序不固定而能带图要求价带负值在下、导带正值在上。real()是因为数值误差可能导致微小虚部1e-15必须清除。3. 绘制高精度石墨烯能带图坐标轴标注、能隙标记与物理量校准能带图的价值不仅在于形状更在于可被同行直接比对的定量信息。MATLAB 绘图必须包含能量零点设为狄拉克点、单位eV、高对称点标签及路径长度归一化。3.1 计算路径长度并归一化横坐标横轴不是 k_x 或 k_y而是沿路径的弧长 s单位为 Å⁻¹ 或归一化单位。这直接影响能带斜率即群速度的物理意义% 计算路径上每点到起点的累积距离欧氏距离 s zeros(N_total, 1); for i 2:N_total ds norm(k_path(i,:) - k_path(i-1,:)); s(i) s(i-1) ds; end % 归一化令 Γ→M 段总长为1.0便于文献对比 s_norm s / s(end) * 3; % 三段总长设为3每段约1单位 % 标记高对称点位置 s_Gamma 0; s_M s(N_GM); s_K s(N_GM N_MK); s_end s(end);3.2 绘制能带并添加物理标注使用plot()而非surf或mesh因为能带是线性路径上的函数figure(Position, [100, 100, 800, 500]); hold on; % 绘制两条能带分支 plot(s_norm, E(:,1), b-, LineWidth, 1.5); % 价带 plot(s_norm, E(:,2), r-, LineWidth, 1.5); % 导带 % 添加高对称点竖线与标签 line([s_Gamma, s_Gamma], ylim, Color, k, LineStyle, --); line([s_M, s_M], ylim, Color, k, LineStyle, --); line([s_K, s_K], ylim, Color, k, LineStyle, --); text(s_Gamma, ylim(2)*0.95, \Gamma, HorizontalAlignment,center,FontSize,12); text(s_M, ylim(2)*0.95, M, HorizontalAlignment,center,FontSize,12); text(s_K, ylim(2)*0.95, K, HorizontalAlignment,center,FontSize,12); % 坐标轴设置 xlabel(Wavevector Path, FontSize, 11); ylabel(Energy (eV), FontSize, 11); title(Graphene Band Structure (Tight-Binding Model), FontSize, 12); grid on; set(gca, FontSize, 10); % 关键物理标注狄拉克点能量 y_d 0; % 理论狄拉克点 line(xlim, [y_d, y_d], Color, g, LineStyle, :, LineWidth, 1); text(xlim(1)0.05*diff(xlim), y_d0.05, Dirac Point, Color, g); % 能隙标注石墨烯为零带隙但可显示最小间隔 min_gap min(E(:,2) - E(:,1)); text(xlim(1)0.05*diff(xlim), ylim(1)0.1*diff(ylim), ... sprintf(Min Gap %.3f eV, min_gap), FontSize, 10);3.2.1 验证狄拉克锥线性色散提取 K 点附近斜率真正的物理验证不是看图而是量化。在 K 点附近拟合 E(k) ∝ |k−K|斜率即费米速度 v_F% 找到K点索引最接近K_cart的点 [~, idx_K] min(sum(abs(k_path - repmat(K_cart, N_total, 1)), 2)); % 取K点附近20个点 half_win 10; idx_near max(1, idx_K-half_win) : min(N_total, idx_Khalf_win); k_dist sqrt(sum((k_path(idx_near,:) - repmat(K_cart, length(idx_near), 1)).^2, 2)); % 只取导带正值部分 E_near E(idx_near, 2); % 线性拟合E v_F * |k-K|v_F 单位为 eV·Å需转换为 m/s p polyfit(k_dist, E_near, 1); v_F_eVA p(1); % eV·Å v_F_ms v_F_eVA * 1.602e-19 / (1.0545718e-34 * 1e10); % 转换1 eV·Å 1.602e-19 J·1e-10 m → / ħ fprintf(Fermi velocity near K point: %.2e m/s (theoretical: 1.0e6 m/s)\n, v_F_ms);注意理论费米速度约为 1×10⁶ m/s。若计算结果偏差超过 20%需检查 d1_cart/d2_cart/d3_cart 是否严格满足 |d_i|a/√3 且夹角120°这是狄拉克锥出现的几何基础。3.3 与文献数据比对导出 CSV 并加载实测能带进行叠加为增强可信度可将计算结果导出并与 ARPES 实验数据如发表在Phys. Rev. Lett.上的石墨烯能带叠加% 导出为CSV含s_norm和两支能带 data_export [s_norm, E]; writematrix(data_export, graphene_band_structure.csv, Delimiter, ,); % 加载外部实验数据假设有 experiment_data.csv两列s_exp, E_exp % experiment readmatrix(experiment_data.csv); % hold on; plot(experiment(:,1), experiment(:,2), ko, MarkerSize, 3, MarkerFaceColor, k);4. 引入次近邻跃迁与自旋轨道耦合扩展模型以解析能带细节与拓扑特性仅最近邻模型给出完美的线性狄拉克锥但真实石墨烯存在微弱次近邻跃迁 t′ ≈ 0.1t 和自旋轨道耦合 λ_SO ≈ 10⁻³ eV它们分别导致能带弯曲和打开微小拓扑能隙。这部分代码展示了如何在不重构主框架的前提下增量式扩展哈密顿量体现 MATLAB 代码的模块化优势。4.1 添加次近邻跃迁 t′修正哈密顿量为 2×2 但含新相位项次近邻跃迁发生在同原子子晶格内A→A 或 B→B其哈密顿量变为H(k) [ ε_A(k) f(k) ] [ f*(k) ε_B(k) ]其中 ε_A(k) 2t′ * (cos(k·d₁′) cos(k·d₂′) cos(k·d₃′))dᵢ′ 是次近邻向量。在 MATLAB 中只需修改f_k函数并新增epsilon_k% 次近邻向量A→A共6个但对称性下可简化为3个独立向量 d1p a1; d2p a2; d3p a1-a2; % 标准选择 t_prime 0.1 * t; % t ≈ 0.28 eV % 修改后的f_k仍只含最近邻 function [f_val, eps_A, eps_B] H_elements(k, d1, d2, d3, d1p, d2p, d3p, t, t_prime) f_val t * (exp(1i*dot(k,d1)) exp(1i*dot(k,d2)) exp(1i*dot(k,d3))); % 同子晶格项ε_A ε_B 2t * Σ cos(k·d_i) cos_sum cos(dot(k,d1p)) cos(dot(k,d2p)) cos(dot(k,d3p)); eps_A 2 * t_prime * cos_sum; eps_B eps_A; % 假设A、B子晶格无差别 end % 主循环中更新H构造 for i 1:N_total k k_path(i, :); [f_val, eps_A, eps_B] H_elements(k, d1_cart, d2_cart, d3_cart, ... d1p, d2p, d3p, t, t_prime); H [eps_A, f_val; conj(f_val), eps_B]; eigvals eig(H); E(i, :) sort(real(eigvals)); end4.1.1 观察能带弯曲计算 Γ 点曲率有效质量次近邻项使能带在 Γ 点附近呈现抛物线形可拟合 E(k) ≈ ℏ²k²/(2m*) 计算有效质量% 取Γ点附近10个点k模长小 idx_Gamma 1:10; k_mag sqrt(sum(k_path(idx_Gamma,:).^2, 2)); % 只取导带正值 E_Gamma E(idx_Gamma, 2); % 抛物线拟合E C * k^2 p_quad polyfit(k_mag.^2, E_Gamma, 1); C p_quad(1); % 单位 eV·Å² % 有效质量 m* ħ²/(2C) 单位电子静止质量 m_e hbar_eVA 6.582e-16; % eV·s m_star_me (hbar_eVA^2) / (2 * C * (1e-10)^2) / 9.109e-31; fprintf(Effective mass at Γ point: %.3f m_e\n, m_star_me);4.2 引入自旋轨道耦合SOC升级为 4×4 哈密顿量以打开能隙SOC 项使石墨烯成为量子自旋霍尔绝缘体候选者。此时哈密顿量需在自旋空间张量积变为 4×4 矩阵。关键在于 SOC 项λ_SO * (σ_z ⊗ s_z)仅作用于次近邻且符号取决于子晶格% 自旋轨道耦合强度eV lambda_SO 1e-3; % 构造4x4哈密顿量基矢顺序 [A↑, A↓, B↑, B↓] % H [ ε_A 0 f 0 ] [ λ_SO 0 0 0 ] [ 0 0 0 0 ] % [ 0 ε_A 0 f ] [ 0 -λ_SO 0 0 ] [ 0 0 0 0 ] % [ f* 0 ε_B 0 ] [ 0 0 -λ_SO 0 ] [ 0 0 0 λ_SO ] % [ 0 f* 0 ε_B ] [ 0 0 0 λ_SO ] [ ... 同子晶格SOC项 ...] % 实际中SOC项为H_SOC i*λ_SO * Σ_{ij} (s_i × d_ij)_z * c_i† c_j导致K/K点符号相反 % 简化处理在K点附近H_SOC ≈ diag([λ_SO, -λ_SO, -λ_SO, λ_SO]) % 因此在K点能带隙 2*λ_SO % 计算K点本征值验证 [~, idx_K] min(sum(abs(k_path - repmat(K_cart, N_total, 1)), 2)); k_K k_path(idx_K, :); [f_val, eps_A, eps_B] H_elements(k_K, d1_cart, d2_cart, d3_cart, d1p, d2p, d3p, t, t_prime); % 4x4 H with SOC H4 zeros(4); H4(1,1) eps_A lambda_SO; H4(2,2) eps_A - lambda_SO; H4(3,3) eps_B - lambda_SO; H4(4,4) eps_B lambda_SO; H4(1,3) f_val; H4(3,1) conj(f_val); H4(2,4) f_val; H4(4,2) conj(f_val); eigvals4 eig(H4); gap_SOC min(diff(sort(real(eigvals4)))); % 最小能隙 fprintf(SOC-induced gap at K point: %.3e eV (target: %.3e eV)\n, gap_SOC, 2*lambda_SO);5. 参数敏感性分析与常见错误排查定位能带畸变、非零带隙或狄拉克点偏移的根本原因当计算出的能带不符合预期如狄拉克锥消失、带隙过大、K点不在零能问题往往不出在绘图而在哈密顿量构建的物理假设或数值实现细节。本章提供一套系统化的自查清单每项均对应可执行的 MATLAB 检查命令。5.1 晶格与倒格子一致性验证确保布里渊区几何正确错误的 a1/a2 会导致 b1/b2 错误进而使 K 点坐标失准% 检查正格子基矢是否正交石墨烯a1·a2应 a²/2 dot_a1a2 dot(a1, a2); a_sq a^2; fprintf(a1·a2 %.6f, expected a²/2 %.6f\n, dot_a1a2, a_sq/2); % 检查倒格子基矢是否满足 b_i · a_j 2π δ_ij check_b1a1 dot(b1, a1); check_b1a2 dot(b1, a2); check_b2a1 dot(b2, a1); check_b2a2 dot(b2, a2); fprintf(b1·a1%.3f (2π≈6.283), b1·a2%.3f (0)\n, check_b1a1, check_b1a2); fprintf(b2·a1%.3f (0), b2·a2%.3f (2π≈6.283)\n, check_b2a1, check_b2a2);5.2 最近邻向量 d1/d2/d3 的模长与夹角验证狄拉克锥存在的充要条件是三个 d_i 模长相等且互成120°d_vecs [d1_cart; d2_cart; d3_cart]; d_norms arrayfun((i) norm(d_vecs(i,:)), 1:3); angles zeros(3,1); for i 1:2 for j i1:3 cos_theta dot(d_vecs(i,:), d_vecs(j,:)) / (d_norms(i)*d_norms(j)); angles((i-1)*2j-i) acos(cos_theta) * 180/pi; % 度 end end fprintf(d_i norms: [%.6f, %.6f, %.6f] (should be equal)\n, d_norms); fprintf(d_i angles: [%.1f°, %.1f°, %.1f°] (should be 120°)\n, angles);5.3 能带计算结果的四项黄金检验运行完主程序后立即执行以下四条命令结果异常即表明模型有根本缺陷% 1. K点本征值是否接近0 [~, idx_K] min(sum(abs(k_path - repmat(K_cart, N_total, 1)), 2)); Ek_K E(idx_K, :); fprintf(K-point eigenvalues: [%.6f, %.6f] eV (should be ~[-0.001, 0.001])\n, Ek_K); % 2. Γ点是否为极值点价带顶、导带底 E_Gamma E(1, :); % Γ在路径起点 dE_dk_Gamma (E(2,2)-E(1,2)) / (s_norm(2)-s_norm(1)); % 导带在Γ点斜率 fprintf(Γ-point conduction band slope: %.6f eV/Å (should be ~0)\n, dE_dk_Gamma); % 3. M点是否为鞍点价带在此处最高 E_M E(N_GM, 2); % M点导带值 E_M_left E(N_GM-1, 2); E_M_right E(N_GM1, 2); fprintf(M-point is saddle? (E_M%.6f, left%.6f, right%.6f) → should be max\n, ... E_M, E_M_left, E_M_right); % 4. 路径总长是否匹配理论布里渊区尺寸 fprintf(Total path length s_max %.6f Å⁻¹ (theoretical Γ-M 2π/(3a)≈%.6f)\n, ... s_norm(end)/3, 2*pi/(3*a));5.3.1 当发现 K 点本征值偏离零时的三步定位法这是最常见问题按优先级排查检查 d1_cart/d2_cart/d3_cart 是否由 a1,a2 严格导出而非手动输入确认 f_k 计算中exp(1i*dot(k,d))的 k 是笛卡尔坐标不是分数坐标验证 K_cart 是否由 b1,b2 的 [1/3,1/3] 分数坐标生成而非 [1/3,0] 等错误组合。执行以下诊断代码% 手动计算K点处的f(k)应为0狄拉克点条件 k_test K_cart; f_test t * (exp(1i*dot(k_test,d1_cart)) exp(1i*dot(k_test,d2_cart)) exp(1i*dot(k_test,d3_cart))); fprintf(|f(K)| %.2e (should be 1e-12)\n, abs(f_test)); % 若不为零说明d_i或K_cart有误提示abs(f_test) 1e-10 即表示晶格几何定义存在系统性错误必须回溯d1_cart的生成逻辑而非调整t或绘图参数。本文还有配套的精品资源点击获取

MATLAB手写紧束缚模型计算石墨烯能带结构

相关推荐

SpringBoot+Vue人事系统实战：解决HR数据断点与事务一致性

AWS无服务器应用开发指南：从Lambda到SAM的架构与实践

缓存更新模式与高并发落地：从Cache Aside到TEI实践

更多资讯

OpenClaw 安全过滤调用内容安全模型报 401？TaoToken 这样给通道地址

AI写论文实用指南：工具选择、技巧方法与合规注意事项全解析

配额超限 429？TaoToken 这样改 Codex 的限流重试

eVTOL毕设从开题到答辩，我常用的AI工具组合

拒绝“憋论文”：用书匠策AI把课程论文从“苦役”变成“可控工程”

Claude Code 并行开 5 个进程，改走 TaoToken 行不行？

今日更新

AWS无服务器应用开发指南：从Lambda到SAM的架构与实践

SpringBoot+Vue人事系统实战：解决HR数据断点与事务一致性

MATLAB手写紧束缚模型计算石墨烯能带结构

本周精选

AI驱动的自动化专利分析平台技术解析

LeetCode-Go 题解：1252. Cells with Odd Values in a Matrix（奇数值单元格计数）双解法剖析

在线答疑系统Java毕设实战：Spring Boot前后端分离与状态流转

本月热门

基于YOLO的焊缝缺陷检测：从数据标注到工业部署全流程实践

分布式服务的决策记录

创业产品的首版范围

读完文章，想聊聊您的网站？