简介平面波展开法PWE是一种基于布里渊区平面波基底展开的数值方法通过求解周期势场中的薛定谔方程获得能量本征值从而绘制电子能量与动量关系的能带曲线是计算晶体、超晶格等AB结构电子性质的核心手段。压缩包内共含2个.m脚本文件整体仅1KB文件分工明确一个用于定义周期势场并构建哈密顿量矩阵另一个负责求解特征值并绘制能带图适用于AB结构等周期性排列体系的能带计算。目前已有576人学习浏览代码精简、可读性好适合凝聚态物理、材料科学方向的研究生和科研人员作为平面波展开法的入门实践模板。借助这套程序可快速修改势场参数、平面波截断数等关键设置直观对比不同周期结构的能带差异理解带隙形成、导电类型等物理机理为后续材料设计或半导体性能预测提供便捷的数值验证基础。1. 平面波展开算能带先搞清楚它在算什么做周期结构仿真的人迟早会撞到「算能带」这件事。光子晶体、声子晶体、超表面、周期性 AB 堆叠结构凡是周期性排列的东西都要回答一个问题在这个周期势场里允许哪些能量/频率传播PWE平面波展开法就是解决这个问题最经典的一种数值手段——把波函数拆成一组平面波的线性组合用布洛赫定理把无穷大周期结构压扁到第一布里渊区最后求解一个哈密顿量矩阵的特征值问题得到的就是能带曲线。这份 MATLAB 资源核心就是两个脚本Untitled.m、Untitled4.m一个搭模型、一个扫波矢算特征值并绘图。适合刚接触能带计算、想从零跑通一条曲线、又不想一上来去啃商业软件的人。2. 平面波展开的骨架布洛赫波、倒格矢与哈密顿量矩阵2.1 为什么偏偏是平面波布洛赫定理给的理由周期结构的本征问题核心难点是势场 (V(r)) 具有空间周期性(V(r R) V(r))其中 (R) 是晶格平移矢量。直接在这个势场里解薛定谔方程数值上很麻烦因为我们要在一个无穷大的空间里描述波函数。布洛赫定理给了关键的一刀周期势场中的本征态一定可以写成[ \psi_k(r) e^{ik\cdot r} u_k(r) ]的形式其中 (u_k(r)) 和势场一样具有晶格周期性。把 (u_k(r)) 再按倒格矢展开成傅里叶级数波函数就变成了一组平面波 (e^{i(kG)\cdot r}) 的线性叠加(G) 是倒格矢。这一步的价值在于原本一个「无限空间里的微分方程」变成了一个「倒空间里、以 (G) 为指标的代数问题」。平面波基底天然满足周期性不需要处理边界条件这就是 PWE 在周期结构里比有限差分、有限元更顺手的根本原因——它把周期性放到基底里而不是放到网格上。这里要强调一个容易懵的点虽然最终算的是倒空间的矩阵但物理图像还在实空间。平面波展开并不会丢失任何实空间信息因为势场的傅里叶系数 (V(G)) 就是从实空间周期势场做积分得到的(V(G) \frac{1}{\Omega}\int_{\Omega} V(r)e^{-iG\cdot r} dr)。你给 PWE 输入什么势场形状它就忠实反映在哈密顿量矩阵的哪些位置上有非零元。换句话说PWE 的精度上限由你建模的势场决定而不是由算法本身决定——这个认识能帮你避免后面很多玄学问题。2.2 哈密顿量矩阵怎么搭从公式到代码映射以电子能带为例单电子薛定谔方程在平面波基底下的矩阵元是[ H_{GG} \frac{\hbar^2}{2m}|kG|^2\delta_{GG} V(G-G) ]等式右边第一项是动能项只在矩阵对角元上出现量级随 (|kG|^2) 增长第二项是势场项由势场傅里叶系数平移 ((G-G)) 决定。代码里要做的事非常机械生成一组倒格矢 (G_1, G_2, ..., G_N)通常按 (|G_x|, |G_y|) 截断比如取 (-N_{max} \le n_x, n_y \le N_{max})对每个波矢 (k)用上面的公式装配出 (N\times N) 的哈密顿量矩阵(N) 是平面波总数用eig()对这个矩阵做特征值分解得到的特征值从小到大排就是当前 (k) 点各能带能量沿布里渊区边界的高对称路径比如 (\Gamma \to X \to M \to \Gamma)逐点重复 2-3把所有特征值连成曲线就是能带图。第一项对角元是确定性的第二项是卷积结构所以在 MATLAB 里最实用的是用两层循环先把势场的傅里叶系数表做好再用下标映射方式填充矩阵。下面给一段我惯用的核心骨架直接对应Untitled.m里最可能长成什么样的部分% 平面波数截断设置一维示例按实际需求扩展到二维 Nmax 5; % 每维倒格矢个数截断 G (-Nmax:Nmax); % 一维倒格矢序列总平面波数 N2*Nmax1 N length(G); % 构造势场傅里叶系数 V(G-G)方阱势势垒高度 V0、宽度 a、晶格常数 d V0 1.0; a 0.3; d 1.0; Vmat zeros(N, N); for i 1:N for j 1:N dG G(i) - G(j); % 倒格矢差值决定矩阵元位置 if dG 0 Vmat(i,j) V0 * a / d; % G0 分量即平均势 else % 方阱势的傅里叶系数解析式 Vmat(i,j) V0 * a / d * sinc(dG * a / (2*pi)); end end end这段代码的逻辑分三层第一层是截断参数Nmax它直接控制矩阵维度Nmax 越大、矩阵越大能带越收敛但计算量按平方甚至立方涨第二层是用双层循环对每一对 ((G, G)) 计算势场矩阵元这里sinc函数来自方阱势的矩形轮廓积分换成其他势场形状就换成对应的解析式或fft数值积分第三层是矩阵装配方式注意势场项依赖差值 (G-G) 而不是绝对位置这是卷积结构在倒空间的典型表现。实际运行中你不需要每次都重写这段通常把势场参数抽到脚本头部方便扫参数。2.3 扫波矢路径归一化决定了能带图横轴长什么样矩阵装配完成后剩下的工作就是沿高对称路径扫描波矢 (k)。这一步的坑主要集中在「横轴长度」上能带图横轴不是均匀放 n 个点就完事而是要按照高对称点在倒空间里的实际距离加权。% 一维情况第一布里渊区 [-pi/d, pi/d]扫描 Gamma 到 X kpath linspace(0, pi/d, 50); % 默认归一Gamma 在 0X 在 pi/d Elist zeros(length(kpath), N); for ik 1:length(kpath) k kpath(ik); H diag( (k G).^2 ) Vmat; % 动能项用约化单位 hbar^2/2m1 Elist(ik,:) sort(eig(H)); % 排序能量本征值 end plot(kpath, Elist, LineWidth, 1.2);注意sort(eig(...))这一步不是可选项。eig()返回的特征值顺序在相邻 (k) 点间不保证连续如果直接画曲线会在带交叉处出现锯齿跳跃。先排序再画是确保曲线连续性的最低成本手段。另外动能项写成(k G).^2意味着你已经把手伸进了约化单位制(\hbar^2/2m1)后面把所有能量换算回物理单位时记得乘回系数——这是新手最容易忽略的比例因子问题我在第 5 章统一展开。3. 拆开两个 MATLAB 脚本从势场定义到能带曲线的完整链路3.1 Untitled.m 大概干了什么主程序与参数区拿到这份资源里两个 MATLAB 脚本站到地面上看最靠谱的判断方式是Untitled.m和Untitled4.m是 MATLAB 默认的未命名文件标题里面装的极可能一个是主计算脚本、一个是辅助函数或备选版本。这类命名下最常见的拆分是主脚本负责定义物理参数、调用核心计算、绘图另一个脚本放着势场定义或者某一特定结构的算例。我拆这类资源的一般顺序是先打开主脚本通常行数多、含figure和plot的那个看参数区再回头看辅助脚本里的函数签名。主脚本参数区通常包含这几组量晶格常数 (d)或二维下的基矢 (a_1, a_2)、势场参数V0 势垒高度、ab 结构的介质占比、截断数 Nmax、波矢路径上的高对称点坐标、以及单位制开关。建议你用同一种单位写完整套参数——要么全用无量纲要么全用 eV 与 Å 配搭不要在中间混进去一个厘米克秒单位后面换算会让你怀疑人生。3.2 势场 特征值计算Untitled4.m 里最可能躺着的辅助函数辅助脚本里最值得细看的是势场傅里叶系数函数。AB 结构两种材料交替的介电常数或势能轮廓是一个矩形波矩形波傅里叶系数是解析的不需要做数值积分代码里大概率是类似这样的实现function VG fourier_coeff_AB(Va, Vb, frac_a, d, G) % Va, Vb: 两种材料的势能/介电常数; frac_a: A层厚度占周期的比例 % d: 晶格常数; G: 需要计算的倒格矢(可以是数组) Vmean Va * frac_a Vb * (1 - frac_a); % 平均势G0 时的值 VG zeros(size(G)); for i 1:length(G) if G(i) 0 VG(i) Vmean; else % 矩形函数傅里叶系数的标准形式 VG(i) (Va - Vb) * sin(pi * frac_a * G(i)) / (pi * G(i)); % 注意这个形式隐含了 G*d 的约化具体比例因子看单位设定 end end end你拿这份资源对照时重点核对三个位置frac_a的默认值是否和 AB 结构的描述一致比如 0.5 一半一半d是否出现在sin的自变量里矩形波傅里叶系数要求无量纲参数 (\pi \cdot frac_a \cdot G \cdot d)漏乘 (d) 是高频错误以及Va - Vb这个差值项前面的符号——它决定带隙出现在 Brillouin 区边界的上还是下虽然不影响定性结论但会影响你对照文献时的匹配度。这个函数把结构参数压缩成了一组倒空间系数算是整个脚本里信息密度最高的一段争论最少、值得逐行读。3.3 绘制能带图画到能用的标准能带图画出来容易画到「能用」有两条隐藏标准第一条是横轴高对称点标注。如果脚本里只写了plot(E, -)而没有按路径坐标标注 (\Gamma)、X、M 的位置图基本只适合自己看过程不适合放进论文或交付。第二条是纵轴量纲能带图纵轴可能是归一化频率光子晶体常见或能量电子能带常见脚本里通常放在最后一节一个ylabel就能看出作者用的是哪套约定。% 绘制能带图隐藏要点高对称点位置标注 Gamma_idx 1; X_idx length(kpath); set(gca, XTick, [kpath(Gamma_idx), kpath(X_idx)]); set(gca, XTickLabel, {\Gamma, X}); ylabel(Energy (eV)); xlabel(Wave vector); grid on; title(Band structure by PWE (AB structure));set(gca,...)这两行是在给横轴做「语义标记」——把数值 tick 换成高对称点符号。做这个动作的注意点是如果路径中间还插了别的点比如 (\Gamma \to X \to M) 三段路径你要在linspace拼接处记录每个高对称点在总数组里的索引而不是像我这段简化代码里取首尾。凡是能带图画完之后带着一条时间线从头标到尾的脚本基本就是完成度比较高的资源。4. AB 周期结构怎么建模倒格子选取与材料参数设置4.1 AB 结构在 PWE 里本质是「一个势场函数」AB 结构听起来像晶体学概念但在 PWE 代码里它的全部物理内涵都浓缩在势场傅里叶系数里。A 原子和 B 原子交替排列宏观上等价于一个周期性的势能轮廓比如一个周期内 A 材料占据 (\mathrm{frac}_a) 的宽度、势能为 (V_a)B 材料占据剩余宽度、势能为 (V_b)。这样建模有一个连带好处当你把 A/B 的宽度比从 0.5 调到 0.3不需要改代码结构只改frac_a一个标量就能重新算一条能带。参数设置的实操建议是先写一个结构参数表再进代码调试。表格里的量要包括晶格常数 (d)决定布里渊区边界位置 (\pi/d)、势垒高度或介电常数比决定带隙宽度、占空比frac_a决定带隙中心位置、以及平面波截断数 Nmax决定数值收敛程度。下面是我惯用的参数表模板参数含义常见取值调试策略d晶格常数1.0无量纲/ 5 Å只影响坐标尺度不影响带隙定性frac_aA 层占空比0.3 ~ 0.7调节带隙宽度与位置Va, Vb两种材料的势能/介电常数相差越大带隙越宽先设差值小验证程序再放大Nmax每维倒格矢截断5 ~ 15从小到大扫描画收敛曲线单位制约化常数(\hbar^2/2m1) 或 eV·Å中途不要改改一次错一次4.2 二维扩展倒格子基矢怎么从这个文件里改出来如果这份脚本只处理一维而你手头要算的是二维光子晶体或二维超晶格那么改动集中在两处倒格矢生成方式和矩阵维度。一维里的G (-Nmax:Nmax)在二维变成% 二维倒格矢枚举正方格子倒格子基矢 b1, b2 NG Nmax; [Gx, Gy] meshgrid(-NG:NG, -NG:NG); Gvec [Gx(:), Gy(:)]; % 每行一个倒格矢 (Gx, Gy) N size(Gvec, 1); % 总平面波数 (2*NG1)^2这段代码把一维的倒格矢序列换成二维网格点集矩阵维度从 (2N_{max}1) 变成 ((2N_{max}1)^2)。一个直观的数字对比Nmax5 时一维只要 11 个平面波二维要 121 个矩阵就变成 121×121eig()的时间还扛得住Nmax 拉到 9 就到 361 维单点特征值求解已经明显变慢所以二维计算里 Nmax 的增长是平方代价建议从 3 开始跑。如果你的结构是三角格子或六角格子倒格矢基矢要换成对应形式但代码结构不变——都是把 (b_1, b_2) 组合成网格亮点。4.3 势场不是方阱时用 FFT 数值化替代解析公式矩形势场的傅里叶系数有解析式但实际结构常常更复杂——渐变折射率、高斯势、粗糙界面。解析公式失效后正确姿势是回到实空间把势场离散到网格上用fft数值求出傅里叶系数。这个做法在脚本改造上非常省事在第一个脚本的势场函数里把解析计算替换成% 实空间离散势场 - 倒空间傅里叶系数 Ngrid 256; % 实空间采样点数取 2 的幂 x linspace(0, d, Ngrid, false); Vreal V0 * (ab_waveform(x, param)); % 任意周期势场波形 Vfft fft(Vreal) / Ngrid; % 归一化后的傅里叶系数 G_idx ...; % 从 fft 输出映射到倒格矢序列 VG_for_G Vfft( G_idx 1 ); % 取需要的倒格矢分量这里最要紧的一步是fft输出的索引映射MATLAB 的fft把零频放在第一个元素正频率在前半段、负频率在后半段直接取下标会让人绕晕。我的习惯是先把实空间坐标 x 摆成从 0 到 d不含 d因为周期性已经说明问题了确保采样点数取 2 的幂上面。另外提醒一句fft数值系数和前面解析式里的系数差一个归一化因子 (\frac{1}{N_{grid}})代码里容易来回差几倍跟能带图对照时先看 G0 分量的值是否等于平均势场不等就说明归一化出了问题。5. PWE 最常见的五个坑截断数、波矢路径与单位制排查5.1 带边对不上文献先怀疑截断数不够现象算出来的能带低频段和文献重合高频段数值偏高或抖动。原因高频能带对应的波函数在倒空间需要更多平面波分量来展开Nmax 太小则高能态没被充分描述。解决做收敛性测试——固定结构参数Nmax 取 3、5、7、9分别算第一带隙上下边界画一条「带隙 vs Nmax」曲线看数值是否趋于水平。我见过有人把 Nmax 设成 3 就直接拿去做参数优化最后白跑了几百组参数全错在同一条低频带上。5.2 路径扫描出现锯齿线特征值排序没做现象能带画出来像一堆断线尤其带交叉点附近来回跳。原因eig()返回的特征值顺序在相邻 k 点之间没有连续性保证是数学上的排序问题而非物理问题。解决每个 k 点特征值算完先排序再存一句Elist(ik,:) sort(eig(H));就能解决。更复杂的方案是按各 k 点间波函数重叠做「能带追踪」但大部分场景排序就足够。5.3 带隙位置偏移且能量差整数倍单位制混了现象带隙中心和理论值一致但绝对值差了一个或几个 (2\pi) 系数。原因动能项 (\frac{\hbar^2}{2m}|kG|^2) 里的约化单位没统一有人把 (k) 用 (1/d) 表示、倒格矢却用 (2\pi/d) 表示混在一起就差系数。解决全脚本通读一遍单位把长度统一到 (\mathrm{Å})、能量统一到 (\mathrm{eV})动能前置系数写进代码常量禁止在循环内做临时换算。我自己的习惯是把系数定义放在参数区第一行后面所有地方只引用变量。5.4 本应有带隙的 AB 结构画出了交叉对称性保护没破现象AB 占空比调到非 0.5能带在布里渊区边界还是硬交叉。原因结构和网格保留过高对称性比如 A、B 层宽度相等时势场具有中心反转对称某些简并被对称性保护着不是计算错误。解决检查frac_a是否确实不等于 0.5如果结构本来就对称比如 0.5/0.5那交叉是物理不是 bug。要验证计算无误可以给势场加一个非对称扰动比如 Va 和 Vb 微调 5%看交叉是否分裂——分裂了说明程序对对称性敏感物理上是对的没分裂才需要查矩阵装配。5.5 计算结果与「常识」不符先查波矢路径再查势场现象算出的能带结构完全看不出半导体特征或者能带数量和理论预期差一倍。原因波矢扫描路径错了或者势场傅里叶系数符号写反。解决先打印几个关键 k 点(\Gamma) 点和布里渊区边界的本征值手动用一维 Kronig-Penney 模型解析结果对照能对上全路径就基本无错。凡是脚本能跑通但物理不对的九成问题出在路径或系数符号上不要先怀疑算法本身——PWE 是成熟方法跑不通只能是实现问题。提示排查顺序建议固定为「截断数 → 单位制 → 路径 → 对称性 → 势场系数」。每次改一个变量画一条新曲线对比不要同时改两个参数不然出问题分不清是哪个改动引起的。6. 把能带曲线验证到可信再用三条检查路径第一条检查路径是「低能极限对照」固定势场很弱V0 趋近 0PWE 算出来的能带应该退化成自由电子抛物线 (\frac{\hbar^2}{2m}k^2)。这个检查能一次性验证动能项系数、单位制、波矢扫描是否正确——如果自由电子极限都画不成抛物线后面算任何真实结构都是白费功夫。第二条检查路径是「带隙边界解析解对照」对于一维周期方阱势可以用 Kronig-Penney 模型的转移矩阵法在另一个独立脚本里算出前几条带边能量和 PWE 结果画在同一张图上。带边误差控制在 1% 以内才算收敛通过。实际经验是 Nmax 取 7 到 9 时一维问题基本收敛取 5 时误差大约百分之几可以接受但不适合做高精度对比。第三条检查路径是「结构对称性验证」把 A/B 材料参数交换Va↔VbPWE 算出的能带必须完全不变再把波矢 k 扫成 -k能带要关于 (\Gamma) 点对称。这两条是程序内部一致性的试金石能排除掉绝大多数矩阵装配和路径设置的问题。结构做对称性验证通过后才算真正把这份资源吃透改造成二维或者加损耗都有可靠的基线。说句实在话我第一次拿这类脚本的时候连eig()要排序都不知道画出来的能带锯齿满天飞还以为是程序写得不对其实只是少了一行排序。从那以后我每次拿到新的能带计算脚本都强制自己先跑一遍自由电子极限、再跑一遍 Kronig-Penney 对照、最后验对称性三步全过才敢往项目里用。整套流程下来你踩过的坑都会变成可复用的检查清单希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/31 9:43:09 · 阅读 1.2k 平面波展开法计算能带：从布洛赫定理到MATLAB实现
简介平面波展开法PWE是一种基于布里渊区平面波基底展开的数值方法通过求解周期势场中的薛定谔方程获得能量本征值从而绘制电子能量与动量关系的能带曲线是计算晶体、超晶格等AB结构电子性质的核心手段。压缩包内共含2个.m脚本文件整体仅1KB文件分工明确一个用于定义周期势场并构建哈密顿量矩阵另一个负责求解特征值并绘制能带图适用于AB结构等周期性排列体系的能带计算。目前已有576人学习浏览代码精简、可读性好适合凝聚态物理、材料科学方向的研究生和科研人员作为平面波展开法的入门实践模板。借助这套程序可快速修改势场参数、平面波截断数等关键设置直观对比不同周期结构的能带差异理解带隙形成、导电类型等物理机理为后续材料设计或半导体性能预测提供便捷的数值验证基础。1. 平面波展开算能带先搞清楚它在算什么做周期结构仿真的人迟早会撞到「算能带」这件事。光子晶体、声子晶体、超表面、周期性 AB 堆叠结构凡是周期性排列的东西都要回答一个问题在这个周期势场里允许哪些能量/频率传播PWE平面波展开法就是解决这个问题最经典的一种数值手段——把波函数拆成一组平面波的线性组合用布洛赫定理把无穷大周期结构压扁到第一布里渊区最后求解一个哈密顿量矩阵的特征值问题得到的就是能带曲线。这份 MATLAB 资源核心就是两个脚本Untitled.m、Untitled4.m一个搭模型、一个扫波矢算特征值并绘图。适合刚接触能带计算、想从零跑通一条曲线、又不想一上来去啃商业软件的人。2. 平面波展开的骨架布洛赫波、倒格矢与哈密顿量矩阵2.1 为什么偏偏是平面波布洛赫定理给的理由周期结构的本征问题核心难点是势场 (V(r)) 具有空间周期性(V(r R) V(r))其中 (R) 是晶格平移矢量。直接在这个势场里解薛定谔方程数值上很麻烦因为我们要在一个无穷大的空间里描述波函数。布洛赫定理给了关键的一刀周期势场中的本征态一定可以写成[ \psi_k(r) e^{ik\cdot r} u_k(r) ]的形式其中 (u_k(r)) 和势场一样具有晶格周期性。把 (u_k(r)) 再按倒格矢展开成傅里叶级数波函数就变成了一组平面波 (e^{i(kG)\cdot r}) 的线性叠加(G) 是倒格矢。这一步的价值在于原本一个「无限空间里的微分方程」变成了一个「倒空间里、以 (G) 为指标的代数问题」。平面波基底天然满足周期性不需要处理边界条件这就是 PWE 在周期结构里比有限差分、有限元更顺手的根本原因——它把周期性放到基底里而不是放到网格上。这里要强调一个容易懵的点虽然最终算的是倒空间的矩阵但物理图像还在实空间。平面波展开并不会丢失任何实空间信息因为势场的傅里叶系数 (V(G)) 就是从实空间周期势场做积分得到的(V(G) \frac{1}{\Omega}\int_{\Omega} V(r)e^{-iG\cdot r} dr)。你给 PWE 输入什么势场形状它就忠实反映在哈密顿量矩阵的哪些位置上有非零元。换句话说PWE 的精度上限由你建模的势场决定而不是由算法本身决定——这个认识能帮你避免后面很多玄学问题。2.2 哈密顿量矩阵怎么搭从公式到代码映射以电子能带为例单电子薛定谔方程在平面波基底下的矩阵元是[ H_{GG} \frac{\hbar^2}{2m}|kG|^2\delta_{GG} V(G-G) ]等式右边第一项是动能项只在矩阵对角元上出现量级随 (|kG|^2) 增长第二项是势场项由势场傅里叶系数平移 ((G-G)) 决定。代码里要做的事非常机械生成一组倒格矢 (G_1, G_2, ..., G_N)通常按 (|G_x|, |G_y|) 截断比如取 (-N_{max} \le n_x, n_y \le N_{max})对每个波矢 (k)用上面的公式装配出 (N\times N) 的哈密顿量矩阵(N) 是平面波总数用eig()对这个矩阵做特征值分解得到的特征值从小到大排就是当前 (k) 点各能带能量沿布里渊区边界的高对称路径比如 (\Gamma \to X \to M \to \Gamma)逐点重复 2-3把所有特征值连成曲线就是能带图。第一项对角元是确定性的第二项是卷积结构所以在 MATLAB 里最实用的是用两层循环先把势场的傅里叶系数表做好再用下标映射方式填充矩阵。下面给一段我惯用的核心骨架直接对应Untitled.m里最可能长成什么样的部分% 平面波数截断设置一维示例按实际需求扩展到二维 Nmax 5; % 每维倒格矢个数截断 G (-Nmax:Nmax); % 一维倒格矢序列总平面波数 N2*Nmax1 N length(G); % 构造势场傅里叶系数 V(G-G)方阱势势垒高度 V0、宽度 a、晶格常数 d V0 1.0; a 0.3; d 1.0; Vmat zeros(N, N); for i 1:N for j 1:N dG G(i) - G(j); % 倒格矢差值决定矩阵元位置 if dG 0 Vmat(i,j) V0 * a / d; % G0 分量即平均势 else % 方阱势的傅里叶系数解析式 Vmat(i,j) V0 * a / d * sinc(dG * a / (2*pi)); end end end这段代码的逻辑分三层第一层是截断参数Nmax它直接控制矩阵维度Nmax 越大、矩阵越大能带越收敛但计算量按平方甚至立方涨第二层是用双层循环对每一对 ((G, G)) 计算势场矩阵元这里sinc函数来自方阱势的矩形轮廓积分换成其他势场形状就换成对应的解析式或fft数值积分第三层是矩阵装配方式注意势场项依赖差值 (G-G) 而不是绝对位置这是卷积结构在倒空间的典型表现。实际运行中你不需要每次都重写这段通常把势场参数抽到脚本头部方便扫参数。2.3 扫波矢路径归一化决定了能带图横轴长什么样矩阵装配完成后剩下的工作就是沿高对称路径扫描波矢 (k)。这一步的坑主要集中在「横轴长度」上能带图横轴不是均匀放 n 个点就完事而是要按照高对称点在倒空间里的实际距离加权。% 一维情况第一布里渊区 [-pi/d, pi/d]扫描 Gamma 到 X kpath linspace(0, pi/d, 50); % 默认归一Gamma 在 0X 在 pi/d Elist zeros(length(kpath), N); for ik 1:length(kpath) k kpath(ik); H diag( (k G).^2 ) Vmat; % 动能项用约化单位 hbar^2/2m1 Elist(ik,:) sort(eig(H)); % 排序能量本征值 end plot(kpath, Elist, LineWidth, 1.2);注意sort(eig(...))这一步不是可选项。eig()返回的特征值顺序在相邻 (k) 点间不保证连续如果直接画曲线会在带交叉处出现锯齿跳跃。先排序再画是确保曲线连续性的最低成本手段。另外动能项写成(k G).^2意味着你已经把手伸进了约化单位制(\hbar^2/2m1)后面把所有能量换算回物理单位时记得乘回系数——这是新手最容易忽略的比例因子问题我在第 5 章统一展开。3. 拆开两个 MATLAB 脚本从势场定义到能带曲线的完整链路3.1 Untitled.m 大概干了什么主程序与参数区拿到这份资源里两个 MATLAB 脚本站到地面上看最靠谱的判断方式是Untitled.m和Untitled4.m是 MATLAB 默认的未命名文件标题里面装的极可能一个是主计算脚本、一个是辅助函数或备选版本。这类命名下最常见的拆分是主脚本负责定义物理参数、调用核心计算、绘图另一个脚本放着势场定义或者某一特定结构的算例。我拆这类资源的一般顺序是先打开主脚本通常行数多、含figure和plot的那个看参数区再回头看辅助脚本里的函数签名。主脚本参数区通常包含这几组量晶格常数 (d)或二维下的基矢 (a_1, a_2)、势场参数V0 势垒高度、ab 结构的介质占比、截断数 Nmax、波矢路径上的高对称点坐标、以及单位制开关。建议你用同一种单位写完整套参数——要么全用无量纲要么全用 eV 与 Å 配搭不要在中间混进去一个厘米克秒单位后面换算会让你怀疑人生。3.2 势场 特征值计算Untitled4.m 里最可能躺着的辅助函数辅助脚本里最值得细看的是势场傅里叶系数函数。AB 结构两种材料交替的介电常数或势能轮廓是一个矩形波矩形波傅里叶系数是解析的不需要做数值积分代码里大概率是类似这样的实现function VG fourier_coeff_AB(Va, Vb, frac_a, d, G) % Va, Vb: 两种材料的势能/介电常数; frac_a: A层厚度占周期的比例 % d: 晶格常数; G: 需要计算的倒格矢(可以是数组) Vmean Va * frac_a Vb * (1 - frac_a); % 平均势G0 时的值 VG zeros(size(G)); for i 1:length(G) if G(i) 0 VG(i) Vmean; else % 矩形函数傅里叶系数的标准形式 VG(i) (Va - Vb) * sin(pi * frac_a * G(i)) / (pi * G(i)); % 注意这个形式隐含了 G*d 的约化具体比例因子看单位设定 end end end你拿这份资源对照时重点核对三个位置frac_a的默认值是否和 AB 结构的描述一致比如 0.5 一半一半d是否出现在sin的自变量里矩形波傅里叶系数要求无量纲参数 (\pi \cdot frac_a \cdot G \cdot d)漏乘 (d) 是高频错误以及Va - Vb这个差值项前面的符号——它决定带隙出现在 Brillouin 区边界的上还是下虽然不影响定性结论但会影响你对照文献时的匹配度。这个函数把结构参数压缩成了一组倒空间系数算是整个脚本里信息密度最高的一段争论最少、值得逐行读。3.3 绘制能带图画到能用的标准能带图画出来容易画到「能用」有两条隐藏标准第一条是横轴高对称点标注。如果脚本里只写了plot(E, -)而没有按路径坐标标注 (\Gamma)、X、M 的位置图基本只适合自己看过程不适合放进论文或交付。第二条是纵轴量纲能带图纵轴可能是归一化频率光子晶体常见或能量电子能带常见脚本里通常放在最后一节一个ylabel就能看出作者用的是哪套约定。% 绘制能带图隐藏要点高对称点位置标注 Gamma_idx 1; X_idx length(kpath); set(gca, XTick, [kpath(Gamma_idx), kpath(X_idx)]); set(gca, XTickLabel, {\Gamma, X}); ylabel(Energy (eV)); xlabel(Wave vector); grid on; title(Band structure by PWE (AB structure));set(gca,...)这两行是在给横轴做「语义标记」——把数值 tick 换成高对称点符号。做这个动作的注意点是如果路径中间还插了别的点比如 (\Gamma \to X \to M) 三段路径你要在linspace拼接处记录每个高对称点在总数组里的索引而不是像我这段简化代码里取首尾。凡是能带图画完之后带着一条时间线从头标到尾的脚本基本就是完成度比较高的资源。4. AB 周期结构怎么建模倒格子选取与材料参数设置4.1 AB 结构在 PWE 里本质是「一个势场函数」AB 结构听起来像晶体学概念但在 PWE 代码里它的全部物理内涵都浓缩在势场傅里叶系数里。A 原子和 B 原子交替排列宏观上等价于一个周期性的势能轮廓比如一个周期内 A 材料占据 (\mathrm{frac}_a) 的宽度、势能为 (V_a)B 材料占据剩余宽度、势能为 (V_b)。这样建模有一个连带好处当你把 A/B 的宽度比从 0.5 调到 0.3不需要改代码结构只改frac_a一个标量就能重新算一条能带。参数设置的实操建议是先写一个结构参数表再进代码调试。表格里的量要包括晶格常数 (d)决定布里渊区边界位置 (\pi/d)、势垒高度或介电常数比决定带隙宽度、占空比frac_a决定带隙中心位置、以及平面波截断数 Nmax决定数值收敛程度。下面是我惯用的参数表模板参数含义常见取值调试策略d晶格常数1.0无量纲/ 5 Å只影响坐标尺度不影响带隙定性frac_aA 层占空比0.3 ~ 0.7调节带隙宽度与位置Va, Vb两种材料的势能/介电常数相差越大带隙越宽先设差值小验证程序再放大Nmax每维倒格矢截断5 ~ 15从小到大扫描画收敛曲线单位制约化常数(\hbar^2/2m1) 或 eV·Å中途不要改改一次错一次4.2 二维扩展倒格子基矢怎么从这个文件里改出来如果这份脚本只处理一维而你手头要算的是二维光子晶体或二维超晶格那么改动集中在两处倒格矢生成方式和矩阵维度。一维里的G (-Nmax:Nmax)在二维变成% 二维倒格矢枚举正方格子倒格子基矢 b1, b2 NG Nmax; [Gx, Gy] meshgrid(-NG:NG, -NG:NG); Gvec [Gx(:), Gy(:)]; % 每行一个倒格矢 (Gx, Gy) N size(Gvec, 1); % 总平面波数 (2*NG1)^2这段代码把一维的倒格矢序列换成二维网格点集矩阵维度从 (2N_{max}1) 变成 ((2N_{max}1)^2)。一个直观的数字对比Nmax5 时一维只要 11 个平面波二维要 121 个矩阵就变成 121×121eig()的时间还扛得住Nmax 拉到 9 就到 361 维单点特征值求解已经明显变慢所以二维计算里 Nmax 的增长是平方代价建议从 3 开始跑。如果你的结构是三角格子或六角格子倒格矢基矢要换成对应形式但代码结构不变——都是把 (b_1, b_2) 组合成网格亮点。4.3 势场不是方阱时用 FFT 数值化替代解析公式矩形势场的傅里叶系数有解析式但实际结构常常更复杂——渐变折射率、高斯势、粗糙界面。解析公式失效后正确姿势是回到实空间把势场离散到网格上用fft数值求出傅里叶系数。这个做法在脚本改造上非常省事在第一个脚本的势场函数里把解析计算替换成% 实空间离散势场 - 倒空间傅里叶系数 Ngrid 256; % 实空间采样点数取 2 的幂 x linspace(0, d, Ngrid, false); Vreal V0 * (ab_waveform(x, param)); % 任意周期势场波形 Vfft fft(Vreal) / Ngrid; % 归一化后的傅里叶系数 G_idx ...; % 从 fft 输出映射到倒格矢序列 VG_for_G Vfft( G_idx 1 ); % 取需要的倒格矢分量这里最要紧的一步是fft输出的索引映射MATLAB 的fft把零频放在第一个元素正频率在前半段、负频率在后半段直接取下标会让人绕晕。我的习惯是先把实空间坐标 x 摆成从 0 到 d不含 d因为周期性已经说明问题了确保采样点数取 2 的幂上面。另外提醒一句fft数值系数和前面解析式里的系数差一个归一化因子 (\frac{1}{N_{grid}})代码里容易来回差几倍跟能带图对照时先看 G0 分量的值是否等于平均势场不等就说明归一化出了问题。5. PWE 最常见的五个坑截断数、波矢路径与单位制排查5.1 带边对不上文献先怀疑截断数不够现象算出来的能带低频段和文献重合高频段数值偏高或抖动。原因高频能带对应的波函数在倒空间需要更多平面波分量来展开Nmax 太小则高能态没被充分描述。解决做收敛性测试——固定结构参数Nmax 取 3、5、7、9分别算第一带隙上下边界画一条「带隙 vs Nmax」曲线看数值是否趋于水平。我见过有人把 Nmax 设成 3 就直接拿去做参数优化最后白跑了几百组参数全错在同一条低频带上。5.2 路径扫描出现锯齿线特征值排序没做现象能带画出来像一堆断线尤其带交叉点附近来回跳。原因eig()返回的特征值顺序在相邻 k 点之间没有连续性保证是数学上的排序问题而非物理问题。解决每个 k 点特征值算完先排序再存一句Elist(ik,:) sort(eig(H));就能解决。更复杂的方案是按各 k 点间波函数重叠做「能带追踪」但大部分场景排序就足够。5.3 带隙位置偏移且能量差整数倍单位制混了现象带隙中心和理论值一致但绝对值差了一个或几个 (2\pi) 系数。原因动能项 (\frac{\hbar^2}{2m}|kG|^2) 里的约化单位没统一有人把 (k) 用 (1/d) 表示、倒格矢却用 (2\pi/d) 表示混在一起就差系数。解决全脚本通读一遍单位把长度统一到 (\mathrm{Å})、能量统一到 (\mathrm{eV})动能前置系数写进代码常量禁止在循环内做临时换算。我自己的习惯是把系数定义放在参数区第一行后面所有地方只引用变量。5.4 本应有带隙的 AB 结构画出了交叉对称性保护没破现象AB 占空比调到非 0.5能带在布里渊区边界还是硬交叉。原因结构和网格保留过高对称性比如 A、B 层宽度相等时势场具有中心反转对称某些简并被对称性保护着不是计算错误。解决检查frac_a是否确实不等于 0.5如果结构本来就对称比如 0.5/0.5那交叉是物理不是 bug。要验证计算无误可以给势场加一个非对称扰动比如 Va 和 Vb 微调 5%看交叉是否分裂——分裂了说明程序对对称性敏感物理上是对的没分裂才需要查矩阵装配。5.5 计算结果与「常识」不符先查波矢路径再查势场现象算出的能带结构完全看不出半导体特征或者能带数量和理论预期差一倍。原因波矢扫描路径错了或者势场傅里叶系数符号写反。解决先打印几个关键 k 点(\Gamma) 点和布里渊区边界的本征值手动用一维 Kronig-Penney 模型解析结果对照能对上全路径就基本无错。凡是脚本能跑通但物理不对的九成问题出在路径或系数符号上不要先怀疑算法本身——PWE 是成熟方法跑不通只能是实现问题。提示排查顺序建议固定为「截断数 → 单位制 → 路径 → 对称性 → 势场系数」。每次改一个变量画一条新曲线对比不要同时改两个参数不然出问题分不清是哪个改动引起的。6. 把能带曲线验证到可信再用三条检查路径第一条检查路径是「低能极限对照」固定势场很弱V0 趋近 0PWE 算出来的能带应该退化成自由电子抛物线 (\frac{\hbar^2}{2m}k^2)。这个检查能一次性验证动能项系数、单位制、波矢扫描是否正确——如果自由电子极限都画不成抛物线后面算任何真实结构都是白费功夫。第二条检查路径是「带隙边界解析解对照」对于一维周期方阱势可以用 Kronig-Penney 模型的转移矩阵法在另一个独立脚本里算出前几条带边能量和 PWE 结果画在同一张图上。带边误差控制在 1% 以内才算收敛通过。实际经验是 Nmax 取 7 到 9 时一维问题基本收敛取 5 时误差大约百分之几可以接受但不适合做高精度对比。第三条检查路径是「结构对称性验证」把 A/B 材料参数交换Va↔VbPWE 算出的能带必须完全不变再把波矢 k 扫成 -k能带要关于 (\Gamma) 点对称。这两条是程序内部一致性的试金石能排除掉绝大多数矩阵装配和路径设置的问题。结构做对称性验证通过后才算真正把这份资源吃透改造成二维或者加损耗都有可靠的基线。说句实在话我第一次拿这类脚本的时候连eig()要排序都不知道画出来的能带锯齿满天飞还以为是程序写得不对其实只是少了一行排序。从那以后我每次拿到新的能带计算脚本都强制自己先跑一遍自由电子极限、再跑一遍 Kronig-Penney 对照、最后验对称性三步全过才敢往项目里用。整套流程下来你踩过的坑都会变成可复用的检查清单希望帮到你。本文还有配套的精品资源点击获取