简介本资源面向凝聚态物理、光学器件设计及电磁仿真方向的研究生与科研人员聚焦狄拉克半金属这一前沿量子材料的光学建模基础工作——提供计算其频率依赖型相对介电常数实部与虚部的关键脚本与数据。包内共2个文件1个MATLAB源码文件.m 1个文本数据文件.txt总大小仅4KB轻量但高度实用MATLAB脚本用于基于能带结构解析狄拉克点并数值计算介电常数实部文本文件则直接给出经理论推导或拟合得到的虚部频域数据二者可无缝导入CST Studio Suite开展电磁响应仿真。已有853人学习下载适用于构建光吸收谱、反射率模型及新型光电探测器的前期参数准备阶段。读者可即取即用快速启动狄拉克半金属在太赫兹至红外波段的光学特性仿真研究显著降低从理论到仿真的技术门槛。1. 项目概述从理论到代码的狄拉克半金属光学性质探索最近在整理凝聚态物理和计算材料学相关的项目时狄拉克半金属这个主题反复被提及尤其是其独特的光学响应特性在光电探测、太赫兹技术等领域展现出了巨大的潜力。很多同学和同行在入门时往往卡在了如何将抽象的物理模型转化为可计算、可可视化的代码这一环。理论公式看得懂但一到用MATLAB实现就无从下手。这其实是一个典型的“最后一公里”问题我们掌握了物理图像但缺乏将其“落地”的工具和路径。这个项目就是聚焦于解决这个问题。我们将以“狄拉克半金属的光学性质”为核心手把手带你用MATLAB完成从理论模型构建、关键物理量计算到结果可视化分析的全过程。它不仅仅是一份代码更是一个完整的计算物理研究框架。无论你是物理、材料专业的研究生需要完成课程作业或开展初步科研还是对拓扑材料计算感兴趣的工程师这个项目都能为你提供一个清晰的、可复现的起点。我们会从最基础的紧束缚模型出发推导出光学电导率的Kubo公式并最终用MATLAB画出那标志性的、与频率成线性关系的光学电导谱让你直观地理解为什么狄拉克半金属对特定波段的光如此“透明”又如此“敏感”。2. 核心物理模型与计算思路拆解要计算光学性质首先必须明确我们描述的对象是什么。对于狄拉克半金属其低能有效模型通常可以用一个二维或三维的狄拉克哈密顿量来描述。这里我们以一个最简单的、各向同性的三维狄拉克半金属为例这也是很多理论研究的起点。2.1 狄拉克半金属的低能有效哈密顿量在动量空间其哈密顿量可以写为 [ \hat{H}(\mathbf{k}) \hbar v_F (\sigma_x k_x \sigma_y k_y \sigma_z k_z) ] 其中( \hbar ) 是约化普朗克常数( v_F ) 是费米速度一个关键材料参数例如在Cd3As2中约为 ( 1.5 \times 10^6 m/s )( \sigma_{x,y,z} ) 是泡利矩阵( \mathbf{k} (k_x, k_y, k_z) ) 是相对于狄拉克点的动量。这个哈密顿量的本征值能谱为 ( E_{\pm}(\mathbf{k}) \pm \hbar v_F |\mathbf{k}| )呈现典型的线性色散关系即狄拉克锥。我们的计算目标——光学电导率描述的是材料在交变电场下的电流响应。对于各向同性的系统我们通常计算其对角分量 ( \sigma_{xx}(\omega) )。这里Kubo公式是我们的核心武器。2.2 光学电导率的Kubo公式框架在松原-久保Kubo线性响应理论的框架下零温度下的光学电导率实部可以表示为 [ \text{Re} \sigma_{xx}(\omega) \frac{\pi e^2}{\hbar \omega} \int \frac{d^3k}{(2\pi)^3} \sum_{n,m} |\langle m|\hat{v}_x|n \rangle|^2 \delta(E_n - E_m - \hbar\omega) [f(E_n) - f(E_m)] ] 这个公式看起来复杂但物理图像非常清晰求和与积分对全布里渊区这里用动量积分近似的所有动量点求和并对所有能带n, m求和。对于两带模型导带和价带n和m就是“”和“-”。速度矩阵元( |\langle m|\hat{v}_x|n \rangle|^2 ) 是连接初态|n和末态|m的速度算符的矩阵元平方。速度算符由哈密顿量对动量的导数给出( \hat{v}_x (1/\hbar) \partial \hat{H} / \partial k_x )。能量守恒( \delta(E_n - E_m - \hbar\omega) ) 是狄拉克δ函数保证光子能量 ( \hbar\omega ) 等于初末态能级差这是光学跃迁必须满足的条件。费米分布差( [f(E_n) - f(E_m)] ) 在零温下简化为 ( \Theta(E_F - E_m) - \Theta(E_F - E_n) )其中 ( \Theta ) 是阶跃函数( E_F ) 是费米能级。这保证了跃迁发生在费米面附近初态被占据而末态空着。对于本征狄拉克半金属费米能级正好在狄拉克点上并且考虑带内跃迁intraband即nm贡献通常被忽略需要特别处理主要贡献来自带间跃迁interband即从价带“-”到导带“”。经过一系列推导涉及泡利矩阵的运算和动量空间积分可以得到一个非常简洁且重要的解析结果 [ \text{Re} \sigma_{xx}(\omega) \frac{e^2}{16\hbar} \frac{\omega}{v_F} \quad \text{(对于 } \hbar\omega 2|E_F| \text{)} ]这就是狄拉克半金属光学电导率的标志性线性关系。我们的MATLAB程序最终目标就是要通过数值计算复现这个线性依赖关系并理解其背后的物理。注意实际数值计算中我们不会直接使用这个最终解析式而是会回到Kubo公式的离散形式进行数值积分这样既能验证解析解又能为处理更复杂模型如各向异性、有质量项、考虑杂质散射等打下基础。这是从学习到科研的关键一步。3. MATLAB实现从公式到代码的完整流程接下来我们进入实操环节。我将把整个计算分解为几个清晰的模块并给出详细的MATLAB代码和注释。假设我们的费米能级 ( E_F 0 )本征情况计算一定频率范围内的光学电导率。3.1 参数定义与动量空间网格化首先我们需要设置物理参数和计算参数。动量空间的积分需要离散化选择一个合适的网格范围和分辨率至关重要。%% 参数设置 clear; close all; clc; % 物理常数 e 1.602e-19; % 元电荷单位 C hbar 1.0546e-34; % 约化普朗克常数单位 J*s eps0 8.854e-12; % 真空介电常数 % 材料参数 (以典型狄拉克半金属为例) vF 1.5e6; % 费米速度单位 m/s Ef 0.0; % 费米能级设定在狄拉克点单位 eV Ef_J Ef * e; % 转换为焦耳 % 计算参数 % 动量空间网格由于线性色散能量截断对应动量截断 E_cutoff 1.0; % 能量截断单位 eV。决定我们积分多大的动量空间区域。 k_cutoff E_cutoff * e / (hbar * vF); % 对应的最大波矢 Nk 150; % 每个动量方向的网格点数。分辨率越高结果越平滑但计算越慢。 % 生成三维动量网格。这里采用均匀网格对于各向同性模型是合适的。 [kx, ky, kz] ndgrid(linspace(-k_cutoff, k_cutoff, Nk), ... linspace(-k_cutoff, k_cutoff, Nk), ... linspace(-k_cutoff, k_cutoff, Nk)); dk (2*k_cutoff)/(Nk-1); % 动量网格间距 V_k (2*k_cutoff)^3; % 动量空间总体积参数选择心得E_cutoff不宜过小否则会丢失高能跃迁贡献也不宜过大否则会引入远离狄拉克点的区域那里低能有效模型可能已经失效。通常取0.5-2 eV是一个合理的范围具体看所研究材料的能带宽度。Nk是精度与计算成本的权衡。对于三维积分计算量随Nk三次方增长。Nk150在普通台式机上计算可能需要几分钟。初次调试可用Nk50快速查看趋势。这里使用了ndgrid生成三维网格它会返回三个Nk x Nk x Nk的矩阵便于后续向量化计算比三层循环快得多。3.2 能带计算与速度矩阵元对于每一个动量点(kx, ky, kz)我们需要计算其对应的能量本征值和本征态用于计算速度矩阵元。%% 计算能带和速度矩阵元 % 初始化数组 E_plus zeros(size(kx)); % 导带能量 E_minus zeros(size(kx)); % 价带能量 % 泡利矩阵 sigma_x [0,1;1,0]; sigma_y [0,-1i;1i,0]; sigma_z [1,0;0,-1]; % 由于哈密顿量是2x2矩阵我们可以解析求解本征值和本征矢。 % 对于 H hbar*vF*(sigma·k)本征值 E ± hbar*vF * |k| k_norm sqrt(kx.^2 ky.^2 kz.^2); E_plus hbar * vF * k_norm; % 导带能量单位 J E_minus -hbar * vF * k_norm; % 价带能量单位 J % 计算速度算符 v_x (1/hbar) * dH/dkx vF * sigma_x % 速度矩阵元 psi_minus| v_x |psi_plus % 对于任意k点哈密顿量 H d*(sigma·n)其中 d hbar*vF*|k|, n k/|k| % 其本征态旋量部分是 n 在泡利矩阵矢量方向上的本征旋量。 % 我们可以直接利用泡利矩阵的性质计算速度矩阵元的模平方。 % 经过推导对于带间跃迁 (从价带到导带)有 % |psi_| v_x |psi_-|^2 vF^2 * (1 - (kx/|k|)^2) / 2 % 这个表达式是各向同性的体现。 v_matrix_element_sq (vF^2) * (1 - (kx ./ k_norm).^2) / 2; % 注意当 |k| 0 时kx./k_norm 会出现 NaN我们需要处理这个奇点。 % 物理上|k|0的点狄拉克点对积分的贡献为零因为态密度为零我们可以将其设为0。 v_matrix_element_sq(k_norm 0) 0;关键点解析向量化计算所有操作 (k_norm,E_plus,v_matrix_element_sq) 都是对Nk x Nk x Nk的大数组进行的没有使用循环这充分利用了MATLAB的矩阵运算优势速度极快。解析推导的价值我们直接使用了速度矩阵元模平方的解析表达式避免了数值对角化哈密顿量再计算本征态和矩阵元的繁琐过程这大大提高了计算效率也减少了数值误差。这是处理简单模型时的常用技巧。奇点处理k_norm0处的奇点必须处理否则会导致NaN非数污染整个结果。根据物理图像将其设为零是合理的。3.3 光学电导率的数值积分实现这是最核心的一步。我们将Kubo公式中的连续积分离散为求和。δ函数需要用数值方法近似这里我们采用洛伦兹展宽或高斯展宽来模拟有限的能级寿命或计算中的平滑需求。%% 计算光学电导率 (带间跃迁贡献) % 频率范围 (从太赫兹到近红外) omega_THz linspace(0.1, 100, 300); % 单位THz omega omega_THz * 2*pi * 1e12; % 转换为角频率单位 rad/s % 展宽参数 (模拟准粒子寿命或实验分辨率) gamma 0.01 * e / hbar; % 展宽单位 rad/s。对应约 0.01 eV 的能量展宽。 % 使用洛伦兹型函数近似δ函数L(x) (1/pi) * (gamma/2) / (x^2 (gamma/2)^2) % 但更常用的是L(E1-E2-hbarω) (1/pi) * gamma / ((E1-E2-hbarω)^2 gamma^2) Re_sigma_xx zeros(size(omega)); % 初始化光学电导率数组 % 零温费米分布差 % 对于带间跃迁 (从价带m-到导带n)f(E-) - f(E) 1 - 0 1 (因为Ef0E-0, E0) % 所以费米分布差项恒为1。如果Ef不为0则需要用阶跃函数判断。 f_factor ones(size(kx)); % 本例中Ef0所有k点满足E-EfE % 数值积分对动量空间所有点求和 % 注意这是一个三重求和但我们已经将kx,ky,kz向量化所以可以用点乘和求和函数。 % 计算所有k点的带间能级差 DeltaE E_plus - E_minus 2 * hbar * vF * |k| DeltaE E_plus - E_minus; % 单位 J % 预计算一些常量避免在循环中重复计算 prefactor (e^2 * dk^3) / (hbar * V_k); % 积分测度和常系数 % Kubo公式中的因子是 pi/(hbar * omega) * |v|^2 * δ(ΔE - hbarω) * f_factor % 我们将其拆解。 fprintf(开始计算光学电导率共有 %d 个频率点...\n, length(omega)); tic; % 开始计时 for i 1:length(omega) hbar_omega hbar * omega(i); % 洛伦兹展宽的δ函数近似 % 注意这里我们计算的是 Re σ(ω)公式中已有π因子而洛伦兹函数积分面积为1。 % 一种更直接的写法是δ(ΔE - hbarω) ≈ (1/pi) * gamma / ((ΔE - hbarω)^2 gamma^2) lorentzian (1/pi) * gamma ./ ((DeltaE - hbar_omega).^2 gamma.^2); % 被积函数|v|^2 * f_factor * lorentzian integrand v_matrix_element_sq .* f_factor .* lorentzian; % 数值积分求和并乘以体积元 dk^3再乘以系数 % 注意我们的动量网格是均匀的每个格点的权重是 dk^3 % 公式中的 ∫d^3k/(2π)^3 ≈ (1/(2π)^3) * Σ * dk^3 integral_value sum(integrand(:)) * dk^3 / ((2*pi)^3); % 计算 Re σ_xx(ω) Re_sigma_xx(i) (pi * e^2 / (hbar * omega(i))) * integral_value; end t_cal toc; fprintf(计算完成耗时 %.2f 秒。\n, t_cal);实现细节与技巧展宽参数gamma的选择这是数值计算中的关键“艺术”。gamma太小δ函数太尖锐需要极高的动量网格分辨率Nk才能捕捉否则结果会充满噪声gamma太大会过度平滑化结果掩盖真实的物理特征如激子峰等。通常gamma应远小于所关心的特征能量尺度如hbar*omega的范围。这里设为0.01 eV是一个合理的起始值。循环与向量化频率循环是不可避免的因为每个ω对应不同的δ函数中心。但在每个频率点内对动量点的计算是完全向量化的integrand是整个数组的运算这保证了核心计算的高效。积分权重dk^3 / ((2*pi)^3)这个因子至关重要。dk^3是离散化的体积元(2π)^3来源于动量空间积分测度的定义d^3k/(2π)^3。忘记这个因子会导致结果在量纲和数值上完全错误。单位制一致性全程使用国际单位制SI。能量从eV转换到J乘以e频率从THz转换到rad/s。确保所有物理量在运算前单位统一是物理计算编程不出错的基础。3.4 结果可视化与理论对比计算完成后我们需要用图形直观展示结果并与理论预期进行对比。%% 结果可视化 figure(Position, [100, 100, 1200, 500]); % 子图1光学电导率谱 subplot(1,2,1); plot(omega_THz, real(Re_sigma_xx), b-, LineWidth, 2); xlabel(频率 (THz), FontSize, 12); ylabel(Re \sigma_{xx}(\omega) (S/m), FontSize, 12); % S/m 是电导率单位 Siemens per meter title(狄拉克半金属光学电导率 (带间跃迁), FontSize, 14); grid on; set(gca, FontSize, 11); % 在同一图中叠加理论线性关系 % 理论公式Re σ_xx (e^2/(16*hbar)) * (omega/vF) % 注意单位e (C), hbar (J*s), vF (m/s), omega (rad/s) - 结果单位是 S/m sigma_theory (e^2/(16*hbar)) * (omega/vF); hold on; plot(omega_THz, sigma_theory, r--, LineWidth, 1.5); legend(数值计算, 理论解析 (线性), Location, northwest); hold off; % 子图2在双对数坐标下观察幂律关系 subplot(1,2,2); loglog(omega_THz, real(Re_sigma_xx), bo, MarkerSize, 4); hold on; loglog(omega_THz, sigma_theory, r-, LineWidth, 1.5); xlabel(频率 (THz), FontSize, 12); ylabel(Re \sigma_{xx}(\omega) (S/m), FontSize, 12); title(双对数坐标下的光学电导率, FontSize, 14); grid on; legend(数值计算, 理论解析 (斜率1), Location, northwest); set(gca, FontSize, 11); % 添加注释标出线性区的斜率 % 在双对数坐标中直线代表幂律关系。斜率为1对应线性关系。 text(10, 1e-3, 斜率 ≈ 1, FontSize, 12, Color, r);可视化要点双图对比左图线性坐标展示整体形状右图双对数坐标用于精确判断幂律关系直线斜率即为指数。这是分析标度行为的标准方法。理论曲线叠加将解析公式计算出的曲线以虚线形式叠加可以直观评估数值计算的准确性。在合适的gamma和Nk下两条曲线应该基本重合。单位标注坐标轴标注清晰的单位THz, S/m是科研绘图的基本要求。电导率σ的国际单位是西门子每米S/m1 S 1 Ω^{-1}。4. 关键参数影响分析与计算优化在实际操作中你会发现计算结果严重依赖于几个关键参数。理解它们的影响是调出可靠结果的关键。4.1 展宽参数gamma的选取策略gamma不是一个真实的物理参数除非特意引入散射而是一个数值工具。它的选取原则是下限受限于动量网格间距dk。粗略估计gamma应大于由于离散化造成的能级最小差异。可以先用一个稍大的gamma如0.05 eV计算然后逐步减小观察谱线何时开始出现剧烈的锯齿状噪声。出现噪声时的gamma值就是当前网格分辨率下的下限。上限应小于你想要分辨的谱特征宽度。例如如果你关心0.1 eV附近的细微结构gamma最好小于0.01 eV。实操建议进行一组参数扫描。固定Nk100分别用gamma [0.005, 0.01, 0.02, 0.05] eV进行计算并绘图。你会看到随着gamma减小曲线在低频区hbarω接近0的峰值变得更尖锐这是带间跃迁阈值行为但高频部分可能噪声增加。选择一个能使曲线平滑且保留物理特征的折中值。4.2 动量网格分辨率Nk与截断E_cutoff的权衡这是一个计算精度与成本的平衡问题。Nk的影响Nk直接决定积分精度。增加Nk可以降低由离散化带来的误差允许使用更小的gamma来揭示更精细的结构但计算量以 (O(Nk^3)) 增长。一个实用的检查方法是将Nk提高一倍例如从80到160比较两次结果。如果主要频段内的电导率相对变化小于5%则可以认为当前的Nk已足够。E_cutoff的影响它定义了积分区域。对于光学电导率高频部分hbarω大的贡献来自大动量态。如果E_cutoff设置过小高频部分的计算结果会偏低因为高能态的跃迁被截断了。你需要确保E_cutoff大于你感兴趣的频率范围对应的最大能量hbar * max(omega)。通常设置E_cutoff 2 * max(hbar*omega)是一个安全的经验法则。重要心得在正式运行长时间计算前务必先用低分辨率参数如Nk50,E_cutoff0.5eV快速跑一遍。这能在几分钟内帮你确认代码逻辑是否正确、参数设置是否合理、结果趋势是否符合预期避免浪费数小时甚至数天时间在错误的计算上。4.3 费米能级Ef非零情况的处理当费米能级不在狄拉克点掺杂或门压调控时计算会发生变化。主要修改在费米分布差因子f_factor和跃迁条件。修改f_factorf_factor heaviside(Ef_J - E_minus, 0) - heaviside(Ef_J - E_plus, 0)。这里heaviside是阶跃函数第二个参数0定义了在0点的值。这个表达式确保了只有初态价带能量低于Ef且末态导带能量高于Ef的跃迁才被允许。修改DeltaE的筛选在带间跃迁公式中DeltaE仍然是E_plus - E_minus。但f_factor会自动将不满足占据条件的跃迁权重设为零。理论公式的变化此时解析公式变为Re σ_xx(ω) (e^2/16hbar) * (ω/vF) * Θ(ħω - 2|Ef|)。即存在一个吸收边ħω 2|Ef|只有当光子能量大于两倍费米能级时带间跃迁才会发生。你的数值结果应该清晰地显示出这个吸收阈值。在代码中实现这一点只需修改f_factor的计算部分并注意在可视化时理论曲线也要乘以阶跃函数heaviside(hbar*omega - 2*abs(Ef_J), 0)。5. 常见问题排查与性能优化技巧即使按照步骤操作你也可能会遇到一些问题。这里列出一些典型情况及解决方法。5.1 计算结果为NaN或Inf原因1动量空间原点奇点。在计算v_matrix_element_sq时k_norm作为分母出现。当k_norm0时会出现0/0或除以零的情况。解决在计算v_matrix_element_sq后立即添加一行v_matrix_element_sq(k_norm 0) 0;。物理上该点贡献为零。原因2频率omega包含零。在计算Re_sigma_xx时公式中有1/omega项。如果omega数组的第一个值是0会导致除以零。解决定义频率范围时从一个小正数开始例如linspace(1e-2, 100, 300)。从零开始本身物理上也不合理静态极限需单独处理。5.2 计算速度过慢三维动量空间积分是计算瓶颈。除了使用向量化还有以下优化手段利用对称性对于各向同性系统光学电导率只与|k|有关。可以将三维积分转化为对k模的一维积分计算量从 (O(Nk^3)) 降至 (O(Nk))。这需要将积分测度从d^3k变为4π k^2 dk并相应地修改矩阵元的角向平均表达式。这是大幅提速的最有效方法。并行计算频率循环for i 1:length(omega)是独立的非常适合用parfor并行。确保你的MATLAB安装了Parallel Computing Toolbox并在循环前用parpool启动工作进程。parpool(local, 4); % 启动4个工作进程 parfor i 1:length(omega) % ... 循环体内的计算 ... end注意循环内的变量需要满足parfor的使用规则如切片变量。降低精度需求在调试和寻找趋势阶段果断使用较小的Nk如50和较大的gamma。5.3 数值结果与理论曲线偏差大低频区不匹配在ω - 0时数值结果可能不为零而理论预测为零对于本征情况。这通常是gamma过大导致的。gamma相当于给δ函数一个宽度使得在ħω gamma的范围内也有非零的跃迁概率。尝试减小gamma同时适当增加Nk。高频区斜率偏离1在双对数坐标下高频区直线斜率明显大于或小于1。检查E_cutoff如果E_cutoff不够大高频部分积分不完整会导致斜率下降。增大E_cutoff。检查gamma的影响过大的gamma会平滑掉高频部分的线性行为也可能影响斜率。尝试在保证曲线平滑的前提下减小gamma。检查单位确保理论曲线sigma_theory和数值曲线Re_sigma_xx使用了完全一致的物理常数和单位。最好在代码开头统一定义所有常数。整体幅值偏差数值结果的整体幅值y轴刻度与理论值有系统性差异。检查积分权重因子反复核对dk^3 / ((2π)^3)这个因子是否正确。这是最常见的错误来源。检查vF值确认数值计算和理论曲线中使用的费米速度vF是同一个值。5.4 内存不足错误当Nk较大时如200kx, ky, kz三个矩阵每个都有200^3 8e6个元素以双精度8字节存储每个矩阵约64 MB三个就是192 MB。后续计算的中间变量如E_plus,integrand也会占用类似大小的内存很容易超过数GB。解决使用单精度浮点数。在参数定义后添加kx single(kx); ky single(ky); kz single(kz);。单精度变量占用内存减半且对于此类计算精度通常足够。此外及时用clear清除不再需要的大变量。这个基于MATLAB的狄拉克半金属光学性质计算框架为你打开了一扇门。你可以在此基础上进行各种拓展研究各向异性模型如Na3Bi只需修改哈密顿量引入有限温度修改费米分布函数f_factor考虑电子-声子散射在gamma中引入频率依赖关系甚至计算其他光学性质如折射率、反射率。最关键的是通过这个亲手实现的过程那些书本上的公式和图表对你而言不再是黑箱而是你可以掌控、可以探究其背后每一个细节的鲜活对象。计算物理的魅力正在于此。本文还有配套的精品资源点击获取 
长河君 · 2026/9/5 12:36:32 · 阅读 1.2k MATLAB实现狄拉克半金属光学电导率计算：从Kubo公式到数值模拟
简介本资源面向凝聚态物理、光学器件设计及电磁仿真方向的研究生与科研人员聚焦狄拉克半金属这一前沿量子材料的光学建模基础工作——提供计算其频率依赖型相对介电常数实部与虚部的关键脚本与数据。包内共2个文件1个MATLAB源码文件.m 1个文本数据文件.txt总大小仅4KB轻量但高度实用MATLAB脚本用于基于能带结构解析狄拉克点并数值计算介电常数实部文本文件则直接给出经理论推导或拟合得到的虚部频域数据二者可无缝导入CST Studio Suite开展电磁响应仿真。已有853人学习下载适用于构建光吸收谱、反射率模型及新型光电探测器的前期参数准备阶段。读者可即取即用快速启动狄拉克半金属在太赫兹至红外波段的光学特性仿真研究显著降低从理论到仿真的技术门槛。1. 项目概述从理论到代码的狄拉克半金属光学性质探索最近在整理凝聚态物理和计算材料学相关的项目时狄拉克半金属这个主题反复被提及尤其是其独特的光学响应特性在光电探测、太赫兹技术等领域展现出了巨大的潜力。很多同学和同行在入门时往往卡在了如何将抽象的物理模型转化为可计算、可可视化的代码这一环。理论公式看得懂但一到用MATLAB实现就无从下手。这其实是一个典型的“最后一公里”问题我们掌握了物理图像但缺乏将其“落地”的工具和路径。这个项目就是聚焦于解决这个问题。我们将以“狄拉克半金属的光学性质”为核心手把手带你用MATLAB完成从理论模型构建、关键物理量计算到结果可视化分析的全过程。它不仅仅是一份代码更是一个完整的计算物理研究框架。无论你是物理、材料专业的研究生需要完成课程作业或开展初步科研还是对拓扑材料计算感兴趣的工程师这个项目都能为你提供一个清晰的、可复现的起点。我们会从最基础的紧束缚模型出发推导出光学电导率的Kubo公式并最终用MATLAB画出那标志性的、与频率成线性关系的光学电导谱让你直观地理解为什么狄拉克半金属对特定波段的光如此“透明”又如此“敏感”。2. 核心物理模型与计算思路拆解要计算光学性质首先必须明确我们描述的对象是什么。对于狄拉克半金属其低能有效模型通常可以用一个二维或三维的狄拉克哈密顿量来描述。这里我们以一个最简单的、各向同性的三维狄拉克半金属为例这也是很多理论研究的起点。2.1 狄拉克半金属的低能有效哈密顿量在动量空间其哈密顿量可以写为 [ \hat{H}(\mathbf{k}) \hbar v_F (\sigma_x k_x \sigma_y k_y \sigma_z k_z) ] 其中( \hbar ) 是约化普朗克常数( v_F ) 是费米速度一个关键材料参数例如在Cd3As2中约为 ( 1.5 \times 10^6 m/s )( \sigma_{x,y,z} ) 是泡利矩阵( \mathbf{k} (k_x, k_y, k_z) ) 是相对于狄拉克点的动量。这个哈密顿量的本征值能谱为 ( E_{\pm}(\mathbf{k}) \pm \hbar v_F |\mathbf{k}| )呈现典型的线性色散关系即狄拉克锥。我们的计算目标——光学电导率描述的是材料在交变电场下的电流响应。对于各向同性的系统我们通常计算其对角分量 ( \sigma_{xx}(\omega) )。这里Kubo公式是我们的核心武器。2.2 光学电导率的Kubo公式框架在松原-久保Kubo线性响应理论的框架下零温度下的光学电导率实部可以表示为 [ \text{Re} \sigma_{xx}(\omega) \frac{\pi e^2}{\hbar \omega} \int \frac{d^3k}{(2\pi)^3} \sum_{n,m} |\langle m|\hat{v}_x|n \rangle|^2 \delta(E_n - E_m - \hbar\omega) [f(E_n) - f(E_m)] ] 这个公式看起来复杂但物理图像非常清晰求和与积分对全布里渊区这里用动量积分近似的所有动量点求和并对所有能带n, m求和。对于两带模型导带和价带n和m就是“”和“-”。速度矩阵元( |\langle m|\hat{v}_x|n \rangle|^2 ) 是连接初态|n和末态|m的速度算符的矩阵元平方。速度算符由哈密顿量对动量的导数给出( \hat{v}_x (1/\hbar) \partial \hat{H} / \partial k_x )。能量守恒( \delta(E_n - E_m - \hbar\omega) ) 是狄拉克δ函数保证光子能量 ( \hbar\omega ) 等于初末态能级差这是光学跃迁必须满足的条件。费米分布差( [f(E_n) - f(E_m)] ) 在零温下简化为 ( \Theta(E_F - E_m) - \Theta(E_F - E_n) )其中 ( \Theta ) 是阶跃函数( E_F ) 是费米能级。这保证了跃迁发生在费米面附近初态被占据而末态空着。对于本征狄拉克半金属费米能级正好在狄拉克点上并且考虑带内跃迁intraband即nm贡献通常被忽略需要特别处理主要贡献来自带间跃迁interband即从价带“-”到导带“”。经过一系列推导涉及泡利矩阵的运算和动量空间积分可以得到一个非常简洁且重要的解析结果 [ \text{Re} \sigma_{xx}(\omega) \frac{e^2}{16\hbar} \frac{\omega}{v_F} \quad \text{(对于 } \hbar\omega 2|E_F| \text{)} ]这就是狄拉克半金属光学电导率的标志性线性关系。我们的MATLAB程序最终目标就是要通过数值计算复现这个线性依赖关系并理解其背后的物理。注意实际数值计算中我们不会直接使用这个最终解析式而是会回到Kubo公式的离散形式进行数值积分这样既能验证解析解又能为处理更复杂模型如各向异性、有质量项、考虑杂质散射等打下基础。这是从学习到科研的关键一步。3. MATLAB实现从公式到代码的完整流程接下来我们进入实操环节。我将把整个计算分解为几个清晰的模块并给出详细的MATLAB代码和注释。假设我们的费米能级 ( E_F 0 )本征情况计算一定频率范围内的光学电导率。3.1 参数定义与动量空间网格化首先我们需要设置物理参数和计算参数。动量空间的积分需要离散化选择一个合适的网格范围和分辨率至关重要。%% 参数设置 clear; close all; clc; % 物理常数 e 1.602e-19; % 元电荷单位 C hbar 1.0546e-34; % 约化普朗克常数单位 J*s eps0 8.854e-12; % 真空介电常数 % 材料参数 (以典型狄拉克半金属为例) vF 1.5e6; % 费米速度单位 m/s Ef 0.0; % 费米能级设定在狄拉克点单位 eV Ef_J Ef * e; % 转换为焦耳 % 计算参数 % 动量空间网格由于线性色散能量截断对应动量截断 E_cutoff 1.0; % 能量截断单位 eV。决定我们积分多大的动量空间区域。 k_cutoff E_cutoff * e / (hbar * vF); % 对应的最大波矢 Nk 150; % 每个动量方向的网格点数。分辨率越高结果越平滑但计算越慢。 % 生成三维动量网格。这里采用均匀网格对于各向同性模型是合适的。 [kx, ky, kz] ndgrid(linspace(-k_cutoff, k_cutoff, Nk), ... linspace(-k_cutoff, k_cutoff, Nk), ... linspace(-k_cutoff, k_cutoff, Nk)); dk (2*k_cutoff)/(Nk-1); % 动量网格间距 V_k (2*k_cutoff)^3; % 动量空间总体积参数选择心得E_cutoff不宜过小否则会丢失高能跃迁贡献也不宜过大否则会引入远离狄拉克点的区域那里低能有效模型可能已经失效。通常取0.5-2 eV是一个合理的范围具体看所研究材料的能带宽度。Nk是精度与计算成本的权衡。对于三维积分计算量随Nk三次方增长。Nk150在普通台式机上计算可能需要几分钟。初次调试可用Nk50快速查看趋势。这里使用了ndgrid生成三维网格它会返回三个Nk x Nk x Nk的矩阵便于后续向量化计算比三层循环快得多。3.2 能带计算与速度矩阵元对于每一个动量点(kx, ky, kz)我们需要计算其对应的能量本征值和本征态用于计算速度矩阵元。%% 计算能带和速度矩阵元 % 初始化数组 E_plus zeros(size(kx)); % 导带能量 E_minus zeros(size(kx)); % 价带能量 % 泡利矩阵 sigma_x [0,1;1,0]; sigma_y [0,-1i;1i,0]; sigma_z [1,0;0,-1]; % 由于哈密顿量是2x2矩阵我们可以解析求解本征值和本征矢。 % 对于 H hbar*vF*(sigma·k)本征值 E ± hbar*vF * |k| k_norm sqrt(kx.^2 ky.^2 kz.^2); E_plus hbar * vF * k_norm; % 导带能量单位 J E_minus -hbar * vF * k_norm; % 价带能量单位 J % 计算速度算符 v_x (1/hbar) * dH/dkx vF * sigma_x % 速度矩阵元 psi_minus| v_x |psi_plus % 对于任意k点哈密顿量 H d*(sigma·n)其中 d hbar*vF*|k|, n k/|k| % 其本征态旋量部分是 n 在泡利矩阵矢量方向上的本征旋量。 % 我们可以直接利用泡利矩阵的性质计算速度矩阵元的模平方。 % 经过推导对于带间跃迁 (从价带到导带)有 % |psi_| v_x |psi_-|^2 vF^2 * (1 - (kx/|k|)^2) / 2 % 这个表达式是各向同性的体现。 v_matrix_element_sq (vF^2) * (1 - (kx ./ k_norm).^2) / 2; % 注意当 |k| 0 时kx./k_norm 会出现 NaN我们需要处理这个奇点。 % 物理上|k|0的点狄拉克点对积分的贡献为零因为态密度为零我们可以将其设为0。 v_matrix_element_sq(k_norm 0) 0;关键点解析向量化计算所有操作 (k_norm,E_plus,v_matrix_element_sq) 都是对Nk x Nk x Nk的大数组进行的没有使用循环这充分利用了MATLAB的矩阵运算优势速度极快。解析推导的价值我们直接使用了速度矩阵元模平方的解析表达式避免了数值对角化哈密顿量再计算本征态和矩阵元的繁琐过程这大大提高了计算效率也减少了数值误差。这是处理简单模型时的常用技巧。奇点处理k_norm0处的奇点必须处理否则会导致NaN非数污染整个结果。根据物理图像将其设为零是合理的。3.3 光学电导率的数值积分实现这是最核心的一步。我们将Kubo公式中的连续积分离散为求和。δ函数需要用数值方法近似这里我们采用洛伦兹展宽或高斯展宽来模拟有限的能级寿命或计算中的平滑需求。%% 计算光学电导率 (带间跃迁贡献) % 频率范围 (从太赫兹到近红外) omega_THz linspace(0.1, 100, 300); % 单位THz omega omega_THz * 2*pi * 1e12; % 转换为角频率单位 rad/s % 展宽参数 (模拟准粒子寿命或实验分辨率) gamma 0.01 * e / hbar; % 展宽单位 rad/s。对应约 0.01 eV 的能量展宽。 % 使用洛伦兹型函数近似δ函数L(x) (1/pi) * (gamma/2) / (x^2 (gamma/2)^2) % 但更常用的是L(E1-E2-hbarω) (1/pi) * gamma / ((E1-E2-hbarω)^2 gamma^2) Re_sigma_xx zeros(size(omega)); % 初始化光学电导率数组 % 零温费米分布差 % 对于带间跃迁 (从价带m-到导带n)f(E-) - f(E) 1 - 0 1 (因为Ef0E-0, E0) % 所以费米分布差项恒为1。如果Ef不为0则需要用阶跃函数判断。 f_factor ones(size(kx)); % 本例中Ef0所有k点满足E-EfE % 数值积分对动量空间所有点求和 % 注意这是一个三重求和但我们已经将kx,ky,kz向量化所以可以用点乘和求和函数。 % 计算所有k点的带间能级差 DeltaE E_plus - E_minus 2 * hbar * vF * |k| DeltaE E_plus - E_minus; % 单位 J % 预计算一些常量避免在循环中重复计算 prefactor (e^2 * dk^3) / (hbar * V_k); % 积分测度和常系数 % Kubo公式中的因子是 pi/(hbar * omega) * |v|^2 * δ(ΔE - hbarω) * f_factor % 我们将其拆解。 fprintf(开始计算光学电导率共有 %d 个频率点...\n, length(omega)); tic; % 开始计时 for i 1:length(omega) hbar_omega hbar * omega(i); % 洛伦兹展宽的δ函数近似 % 注意这里我们计算的是 Re σ(ω)公式中已有π因子而洛伦兹函数积分面积为1。 % 一种更直接的写法是δ(ΔE - hbarω) ≈ (1/pi) * gamma / ((ΔE - hbarω)^2 gamma^2) lorentzian (1/pi) * gamma ./ ((DeltaE - hbar_omega).^2 gamma.^2); % 被积函数|v|^2 * f_factor * lorentzian integrand v_matrix_element_sq .* f_factor .* lorentzian; % 数值积分求和并乘以体积元 dk^3再乘以系数 % 注意我们的动量网格是均匀的每个格点的权重是 dk^3 % 公式中的 ∫d^3k/(2π)^3 ≈ (1/(2π)^3) * Σ * dk^3 integral_value sum(integrand(:)) * dk^3 / ((2*pi)^3); % 计算 Re σ_xx(ω) Re_sigma_xx(i) (pi * e^2 / (hbar * omega(i))) * integral_value; end t_cal toc; fprintf(计算完成耗时 %.2f 秒。\n, t_cal);实现细节与技巧展宽参数gamma的选择这是数值计算中的关键“艺术”。gamma太小δ函数太尖锐需要极高的动量网格分辨率Nk才能捕捉否则结果会充满噪声gamma太大会过度平滑化结果掩盖真实的物理特征如激子峰等。通常gamma应远小于所关心的特征能量尺度如hbar*omega的范围。这里设为0.01 eV是一个合理的起始值。循环与向量化频率循环是不可避免的因为每个ω对应不同的δ函数中心。但在每个频率点内对动量点的计算是完全向量化的integrand是整个数组的运算这保证了核心计算的高效。积分权重dk^3 / ((2*pi)^3)这个因子至关重要。dk^3是离散化的体积元(2π)^3来源于动量空间积分测度的定义d^3k/(2π)^3。忘记这个因子会导致结果在量纲和数值上完全错误。单位制一致性全程使用国际单位制SI。能量从eV转换到J乘以e频率从THz转换到rad/s。确保所有物理量在运算前单位统一是物理计算编程不出错的基础。3.4 结果可视化与理论对比计算完成后我们需要用图形直观展示结果并与理论预期进行对比。%% 结果可视化 figure(Position, [100, 100, 1200, 500]); % 子图1光学电导率谱 subplot(1,2,1); plot(omega_THz, real(Re_sigma_xx), b-, LineWidth, 2); xlabel(频率 (THz), FontSize, 12); ylabel(Re \sigma_{xx}(\omega) (S/m), FontSize, 12); % S/m 是电导率单位 Siemens per meter title(狄拉克半金属光学电导率 (带间跃迁), FontSize, 14); grid on; set(gca, FontSize, 11); % 在同一图中叠加理论线性关系 % 理论公式Re σ_xx (e^2/(16*hbar)) * (omega/vF) % 注意单位e (C), hbar (J*s), vF (m/s), omega (rad/s) - 结果单位是 S/m sigma_theory (e^2/(16*hbar)) * (omega/vF); hold on; plot(omega_THz, sigma_theory, r--, LineWidth, 1.5); legend(数值计算, 理论解析 (线性), Location, northwest); hold off; % 子图2在双对数坐标下观察幂律关系 subplot(1,2,2); loglog(omega_THz, real(Re_sigma_xx), bo, MarkerSize, 4); hold on; loglog(omega_THz, sigma_theory, r-, LineWidth, 1.5); xlabel(频率 (THz), FontSize, 12); ylabel(Re \sigma_{xx}(\omega) (S/m), FontSize, 12); title(双对数坐标下的光学电导率, FontSize, 14); grid on; legend(数值计算, 理论解析 (斜率1), Location, northwest); set(gca, FontSize, 11); % 添加注释标出线性区的斜率 % 在双对数坐标中直线代表幂律关系。斜率为1对应线性关系。 text(10, 1e-3, 斜率 ≈ 1, FontSize, 12, Color, r);可视化要点双图对比左图线性坐标展示整体形状右图双对数坐标用于精确判断幂律关系直线斜率即为指数。这是分析标度行为的标准方法。理论曲线叠加将解析公式计算出的曲线以虚线形式叠加可以直观评估数值计算的准确性。在合适的gamma和Nk下两条曲线应该基本重合。单位标注坐标轴标注清晰的单位THz, S/m是科研绘图的基本要求。电导率σ的国际单位是西门子每米S/m1 S 1 Ω^{-1}。4. 关键参数影响分析与计算优化在实际操作中你会发现计算结果严重依赖于几个关键参数。理解它们的影响是调出可靠结果的关键。4.1 展宽参数gamma的选取策略gamma不是一个真实的物理参数除非特意引入散射而是一个数值工具。它的选取原则是下限受限于动量网格间距dk。粗略估计gamma应大于由于离散化造成的能级最小差异。可以先用一个稍大的gamma如0.05 eV计算然后逐步减小观察谱线何时开始出现剧烈的锯齿状噪声。出现噪声时的gamma值就是当前网格分辨率下的下限。上限应小于你想要分辨的谱特征宽度。例如如果你关心0.1 eV附近的细微结构gamma最好小于0.01 eV。实操建议进行一组参数扫描。固定Nk100分别用gamma [0.005, 0.01, 0.02, 0.05] eV进行计算并绘图。你会看到随着gamma减小曲线在低频区hbarω接近0的峰值变得更尖锐这是带间跃迁阈值行为但高频部分可能噪声增加。选择一个能使曲线平滑且保留物理特征的折中值。4.2 动量网格分辨率Nk与截断E_cutoff的权衡这是一个计算精度与成本的平衡问题。Nk的影响Nk直接决定积分精度。增加Nk可以降低由离散化带来的误差允许使用更小的gamma来揭示更精细的结构但计算量以 (O(Nk^3)) 增长。一个实用的检查方法是将Nk提高一倍例如从80到160比较两次结果。如果主要频段内的电导率相对变化小于5%则可以认为当前的Nk已足够。E_cutoff的影响它定义了积分区域。对于光学电导率高频部分hbarω大的贡献来自大动量态。如果E_cutoff设置过小高频部分的计算结果会偏低因为高能态的跃迁被截断了。你需要确保E_cutoff大于你感兴趣的频率范围对应的最大能量hbar * max(omega)。通常设置E_cutoff 2 * max(hbar*omega)是一个安全的经验法则。重要心得在正式运行长时间计算前务必先用低分辨率参数如Nk50,E_cutoff0.5eV快速跑一遍。这能在几分钟内帮你确认代码逻辑是否正确、参数设置是否合理、结果趋势是否符合预期避免浪费数小时甚至数天时间在错误的计算上。4.3 费米能级Ef非零情况的处理当费米能级不在狄拉克点掺杂或门压调控时计算会发生变化。主要修改在费米分布差因子f_factor和跃迁条件。修改f_factorf_factor heaviside(Ef_J - E_minus, 0) - heaviside(Ef_J - E_plus, 0)。这里heaviside是阶跃函数第二个参数0定义了在0点的值。这个表达式确保了只有初态价带能量低于Ef且末态导带能量高于Ef的跃迁才被允许。修改DeltaE的筛选在带间跃迁公式中DeltaE仍然是E_plus - E_minus。但f_factor会自动将不满足占据条件的跃迁权重设为零。理论公式的变化此时解析公式变为Re σ_xx(ω) (e^2/16hbar) * (ω/vF) * Θ(ħω - 2|Ef|)。即存在一个吸收边ħω 2|Ef|只有当光子能量大于两倍费米能级时带间跃迁才会发生。你的数值结果应该清晰地显示出这个吸收阈值。在代码中实现这一点只需修改f_factor的计算部分并注意在可视化时理论曲线也要乘以阶跃函数heaviside(hbar*omega - 2*abs(Ef_J), 0)。5. 常见问题排查与性能优化技巧即使按照步骤操作你也可能会遇到一些问题。这里列出一些典型情况及解决方法。5.1 计算结果为NaN或Inf原因1动量空间原点奇点。在计算v_matrix_element_sq时k_norm作为分母出现。当k_norm0时会出现0/0或除以零的情况。解决在计算v_matrix_element_sq后立即添加一行v_matrix_element_sq(k_norm 0) 0;。物理上该点贡献为零。原因2频率omega包含零。在计算Re_sigma_xx时公式中有1/omega项。如果omega数组的第一个值是0会导致除以零。解决定义频率范围时从一个小正数开始例如linspace(1e-2, 100, 300)。从零开始本身物理上也不合理静态极限需单独处理。5.2 计算速度过慢三维动量空间积分是计算瓶颈。除了使用向量化还有以下优化手段利用对称性对于各向同性系统光学电导率只与|k|有关。可以将三维积分转化为对k模的一维积分计算量从 (O(Nk^3)) 降至 (O(Nk))。这需要将积分测度从d^3k变为4π k^2 dk并相应地修改矩阵元的角向平均表达式。这是大幅提速的最有效方法。并行计算频率循环for i 1:length(omega)是独立的非常适合用parfor并行。确保你的MATLAB安装了Parallel Computing Toolbox并在循环前用parpool启动工作进程。parpool(local, 4); % 启动4个工作进程 parfor i 1:length(omega) % ... 循环体内的计算 ... end注意循环内的变量需要满足parfor的使用规则如切片变量。降低精度需求在调试和寻找趋势阶段果断使用较小的Nk如50和较大的gamma。5.3 数值结果与理论曲线偏差大低频区不匹配在ω - 0时数值结果可能不为零而理论预测为零对于本征情况。这通常是gamma过大导致的。gamma相当于给δ函数一个宽度使得在ħω gamma的范围内也有非零的跃迁概率。尝试减小gamma同时适当增加Nk。高频区斜率偏离1在双对数坐标下高频区直线斜率明显大于或小于1。检查E_cutoff如果E_cutoff不够大高频部分积分不完整会导致斜率下降。增大E_cutoff。检查gamma的影响过大的gamma会平滑掉高频部分的线性行为也可能影响斜率。尝试在保证曲线平滑的前提下减小gamma。检查单位确保理论曲线sigma_theory和数值曲线Re_sigma_xx使用了完全一致的物理常数和单位。最好在代码开头统一定义所有常数。整体幅值偏差数值结果的整体幅值y轴刻度与理论值有系统性差异。检查积分权重因子反复核对dk^3 / ((2π)^3)这个因子是否正确。这是最常见的错误来源。检查vF值确认数值计算和理论曲线中使用的费米速度vF是同一个值。5.4 内存不足错误当Nk较大时如200kx, ky, kz三个矩阵每个都有200^3 8e6个元素以双精度8字节存储每个矩阵约64 MB三个就是192 MB。后续计算的中间变量如E_plus,integrand也会占用类似大小的内存很容易超过数GB。解决使用单精度浮点数。在参数定义后添加kx single(kx); ky single(ky); kz single(kz);。单精度变量占用内存减半且对于此类计算精度通常足够。此外及时用clear清除不再需要的大变量。这个基于MATLAB的狄拉克半金属光学性质计算框架为你打开了一扇门。你可以在此基础上进行各种拓展研究各向异性模型如Na3Bi只需修改哈密顿量引入有限温度修改费米分布函数f_factor考虑电子-声子散射在gamma中引入频率依赖关系甚至计算其他光学性质如折射率、反射率。最关键的是通过这个亲手实现的过程那些书本上的公式和图表对你而言不再是黑箱而是你可以掌控、可以探究其背后每一个细节的鲜活对象。计算物理的魅力正在于此。本文还有配套的精品资源点击获取