简介针对电子物理与凝聚态物理中周期性结构能带计算问题这份压缩包提供一套基于平面波展开PWE方法的能带曲线计算实现内含2个MATLAB脚本文件.m包体仅1KB尤其适合需要计算AB结构等周期势场能带的研究生、科研人员以及正在学习固体物理或计算材料课程的本科生。脚本覆盖完整计算流程周期势场定义、平面波基组构建、哈密顿量矩阵组装与特征值求解并最终提取能量绘制能带曲线代码精简、注释清晰便于对照理论公式自行修改参数验证不同周期排列下的能带变化也可作为课程作业、毕业设计或科研前期预研的起点。目前已有576人学习下载反映出该主题在能带理论数值模拟和材料设计中的实用价值。通过研读两个m文件读者可快速掌握PWE方法的核心步骤并迁移到光子晶体、声子晶体等其他周期体系对理解材料导电性、光学性质及超晶格设计具有直接帮助。1. PWE for band curve 到底是什么一次拉出整条能带曲线的平面波展开法拿到一个二维周期结构比如空气孔在硅板里排成三角格子你最想先知道的是它能不能开带隙带隙落在哪个频段想回答这两个问题多数人第一选择不是差分也不是有限元而是 PWE平面波展开Plane Wave Expansion。PWE 算能带曲线的思路很直接把周期结构里的场展开成一系列平面波的叠加再借 Bloch 定理把偏微分方程变成一个矩阵特征值问题对布里渊区路径上的每个 k 点扫一遍整条 band curve 就出来了。整个过程不用画网格、不用反复试边界条件十分钟内从几何参数换到一张能判断带隙的曲线图。这篇文章适合做光子晶体、声子晶体或超材料快速评估的人也适合刚接触能带计算但想把公式落到代码里的新手。接下来从方程怎么变成矩阵开始再给你能直接运行的 Python 片段最后把收敛性、参数和常见坑一起讲清楚。2. 从 Bloch 定理到广义本征值问题把波动方程翻译成矩阵2.1 波函数用 Bloch 波展开先看物理模型再写方程PWE 最舒服的入口是二维光子晶体晶格常数 a圆柱半径 r柱内介电常数 ε_c背景介电常数 ε_b柱子在背景里按周期排列。这四组参数决定后面所有结果。对于无源、无损耗、线性介质磁场的波动方程可以写成[ \nabla \times \left( \frac{1}{\varepsilon(\mathbf{r})} \nabla \times \mathbf{H} \right) \frac{\omega^2}{c^2} \mathbf{H} ]其中 ε(r) 是周期函数满足 ε(r)ε(rR)R 是任意晶格矢量。由于结构在 z 方向无限长二维问题可以按极化分开处理。TM 极化对应电场只有 z 分量方程退化为一个标量方程TE 极化对应磁场只有 z 分量方程形式略有不同。这个极化区分非常重要后面组装矩阵时要分别对待。Bloch 定理告诉我们周期结构里的本征场可以写成一个平面波因子乘以一个周期函数[ \mathbf{H}\mathbf{k}(\mathbf{r}) e^{i\mathbf{k}\cdot\mathbf{r}} \mathbf{u}\mathbf{k}(\mathbf{r}) ]而周期函数 u_k(r) 本身又可以展开成倒格矢 G 的平面波叠加。把这两步代进波动方程微分算符作用到指数上就会变成代数运算最终得到一个形如 (\sum_{G} M_{G,G} u_{G} (\omega/c)^2 u_G) 的矩阵方程。这就是 PWE 的全部核心把「解偏微分方程」换成「解矩阵特征值」。你不需要手推每个细节但必须记住后面代码里的每一条矩阵元都来自这一步。2.2 介电常数傅里叶系数直接展开还是逆展开矩阵方程里真正麻烦的是 ε(r) 的出现位置。展开场时方程里会出现 1/ε(r) 与场的乘积这个乘积在倒空间里会变成傅里叶系数的卷积。常见做法是对 1/ε 做傅里叶展开而不是对 ε 直接展开。原因很实际对 ε 直接展开时高介电对比结构的傅里叶级数收敛非常慢带隙附近容易出现伪带对 1/ε 展开则收敛快得多这是 PWE 从论文落到代码时最重要的一个选择。对二维正方晶格里的圆柱结构逆介电常数的傅里叶系数有解析式。令填充比 f πr²/a²则G0 项(\eta_0 f/\varepsilon_c (1-f)/\varepsilon_b)G≠0 项(\eta_G (1/\varepsilon_c - 1/\varepsilon_b) \cdot f \cdot \frac{2J_1(GR)}{GR})其中 J1 是第一类一阶贝塞尔函数。这个公式在很多教材和开源代码里都能对上建议作为自检点如果你用别的结构比如正方形柱或六边形孔解析式会变更稳妥的做法是在实空间单胞里离散采样然后做 FFT 得到傅里叶系数。FFT 方案对任意形状都成立但要注意介电常数突变处的 Gibbs 振荡通常用足够密的网格就能缓解。2.3 组装矩阵与求解器选择实对称矩阵交给 eigh二维正方晶格的倒格矢是 b1(2π/a,0) 和 b2(0,2π/a)任意倒格矢 G m·b1 n·b2。截断数 N 决定了平面波数量总数为 (2N1)²。所有 G 排列成一个列表后就可以对每个 k 点组装矩阵。TM 极化下矩阵元取范数乘积[ M_{G,G} |\mathbf{k}\mathbf{G}| \cdot |\mathbf{k}\mathbf{G}| \cdot \eta(\mathbf{G}-\mathbf{G}) ]TE 极化下矩阵元取点积[ M_{G,G} (\mathbf{k}\mathbf{G}) \cdot (\mathbf{k}\mathbf{G}) \cdot \eta(\mathbf{G}-\mathbf{G}) ]两个矩阵都是实对称矩阵。这个性质非常重要求解器应该选 scipy.linalg.eigh而不是通用的 eig。eigh 会利用对称性走 LAPACK 的专用路径速度快且数值稳定。只有在 N 很大、平面波数量上千、你只关心最低几条带的时候才需要考虑稀疏迭代求解器比如 scipy.sparse.linalg.eigsh。关于 eigsh 的坑第 5 章会专门讲这里先记住一个原则矩阵组装对了能带曲线就成功了一半。3. 用 PWE 计算二维光子晶体能带最小可复现代码与逐段拆解3.1 建立倒格矢与 G 列表先写参数区和倒空间基矢。这里统一取晶格常数 a1所有长度都以 a 为单位最后再换算成归一化频率。import numpy as np from scipy.special import j1 from scipy.linalg import eigh a 1.0 # 晶格常数归一化 r 0.2 * a # 圆柱半径 eps_c 13.0 # 柱体介电常数比如硅 eps_b 1.0 # 背景介电常数空气 N_G 10 # 每个方向的平面波截断数 # 正方晶格倒格矢 b1 2 * np.pi / a * np.array([1.0, 0.0]) b2 2 * np.pi / a * np.array([0.0, 1.0]) # 生成 (2N_G1)^2 个倒格矢 G_list [] for m in range(-N_G, N_G 1): for n in range(-N_G, N_G 1): G_list.append(m * b1 n * b2) G np.array(G_list)这段代码里最关键的是倒格矢的 2π 因子。很多人第一次写 PWE 把 2π 漏掉结果能带整体平移横轴纵轴全部对不上。G_list 的顺序不需要特殊处理只要后续组装矩阵时索引一致即可。N_G10 时平面波总数是 441矩阵维度 441×441对 scipy 的 eigh 来说非常轻松N_G15 时接近 1000 维仍然可接受。3.2 组装 TE/TM 矩阵从傅里叶系数到矩阵元接下来定义逆介电常数的傅里叶系数函数。注意 G0 的项必须单独处理因为贝塞尔函数公式在 G0 处没有定义。def inv_eps_fourier(Gdiff): f np.pi * r**2 / a**2 # 填充比 if np.linalg.norm(Gdiff) 0: return f / eps_c (1 - f) / eps_b x np.linalg.norm(Gdiff) * r return (1.0 / eps_c - 1.0 / eps_b) * f * 2.0 * j1(x) / x这个函数每次被调用都重新算一次填充比工程上可以提前算好。更高效的做法是预计算一个 eta_mat 矩阵把 G[i]-G[j] 对应的傅里叶系数一次填满nG len(G) eta_mat np.zeros((nG, nG)) for i in range(nG): for j in range(nG): eta_mat[i, j] inv_eps_fourier(G[i] - G[j])eta_mat 只依赖结构本身不依赖 k 点所以整个 k 路径上只用算一次。如果你要扫 50 个 k 点这一步能省掉大量重复计算。接下来对每个 k 点组装矩阵。以 TM 极化为例kfrac np.array([[0.0, 0.0], # Gamma [0.5, 0.0], # X [0.5, 0.5], # M [0.0, 0.0]]) # Gamma kvecs kfrac np.array([b1, b2]) def assemble_matrix(kvec, modeTM): M np.zeros((nG, nG)) for i in range(nG): for j in range(nG): kg_i kvec G[i] kg_j kvec G[j] if mode TM: M[i, j] eta_mat[i, j] * np.linalg.norm(kg_i) * np.linalg.norm(kg_j) else: M[i, j] eta_mat[i, j] * np.dot(kg_i, kg_j) return M这里 TM 和 TE 的唯一差别就是范数乘积换成点积。很多现成源码会把两个模式写到同一个函数里用参数切换。矩阵的对称性来自两个地方eta_mat 的对称性以及范数/点积对 i、j 交换不变。所以组装完成后可以直接放心交给 eigh。3.3 求解与能带曲线绘制把本征值换成归一化频率对 k 路径上的每个 k 点调用 eigh特征值默认升序排列省去自己排序的麻烦。本征值是 (ω/c)²画图前要开方并归一化。wlist [] for kvec in kvecs: M assemble_matrix(kvec, modeTM) w2 eigh(M, eigvals_onlyTrue) w np.sqrt(np.maximum(w2, 0.0)) # 防止数值误差产生的负值 wlist.append(w * a / (2 * np.pi)) # 转成 a/λ 单位 wbands np.array(wlist)k 路径的横轴不能直接用分数坐标要换算成倒空间的实际累计长度xcoord [0.0] for i in range(len(kvecs) - 1): xcoord.append(xcoord[-1] np.linalg.norm(kvecs[i 1] - kvecs[i])) xcoord np.array(xcoord)画图时取最低几条带即可import matplotlib.pyplot as plt band_num 8 for bi in range(band_num): plt.plot(xcoord, wbands[:, bi], b-) plt.xlim(0, xcoord[-1]) plt.ylim(0, 0.8) plt.xticks(xcoord, [Γ, X, M, Γ]) plt.ylabel(a/λ)这段代码跑出来的图就是你要的 band curve。如果发现曲线在某个 k 点突然跳变多半是能带交叉的地方 index 对不上这属于正常现象第 5 章会讲怎么处理。整个脚本不到 60 行结构上就是「参数区 → 倒格矢 → k 路径 → 矩阵组装 → 求解 → 画图」后续换三角晶格、换介质柱形状只需要改倒格矢和傅里叶系数两个函数。4. 截断数、填充比、归一化频率PWE 参数设定与收敛性检查4.1 截断数 N_G从趋势到收敛平面波截断数 N_G 是 PWE 里最核心的收敛参数。N_G5 时只有 121 个平面波计算很快但高频带和带隙位置误差很大N_G10 时有 441 个精度明显改善N_G15 时有 961 个已经能应付大多数光子晶体结构。矩阵维度增加一倍求解时间大约增加一个量级因为 eigh 的复杂度是维度的三次方。我一般会按两个阶段跑先用 N_G5 快速扫一遍 k 路径看能带的大致形状和带隙位置确认参数没写错然后再用 N_G10 或 15 跑最终结果。收敛性检查的简单做法是固定某个高对称点比如 Γ 点对比 N_G5、10、15 时第一条带和第二条带的本征频率变化小于 1% 就说明截断够了。带边频率往往比带中部对截断更敏感所以如果要用带隙边界做设计建议在带边处单独检查一次收敛。4.2 填充比与介电常数对比度怎么选填充比 fπr²/a² 直接决定带隙的位置和宽度。二维光子晶体最常见的经验是圆柱半径太小填充比过低带隙很难打开半径太大相邻柱体之间只剩很窄的背景通道PWE 的傅里叶级数收敛变差。对于正方晶格空气孔结构r 在 0.35a 到 0.45a 之间通常能拿到比较宽的带隙对于介质柱结构r 在 0.15a 到 0.2a 附近的组合更常用。这些数值不是死的但可以作为起点。介电常数对比度同样重要。硅ε≈12对空气ε1的对比度是 12:1是 PWE 能稳定处理的典型区间如果换成锗ε≈16或者更高对比材料PWE 会开始吃力伪带的概率明显上升。对比度低于 3:1 时往往没有完整带隙只有带隙很窄的局域效应。所以在跑参数扫描之前先确认你的结构在 PWE 的适用边界内。下表是常见参数的参考范围参数常用范围主要影响截断数 N_G520平面波数量收敛速度填充比 f0.050.5带隙位置与宽度介电常数对比度3:115:1带隙深度傅里叶收敛晶格类型正方 / 三角布里渊区路径与简并4.3 归一化坐标不统一量纲的曲线没法对比能带曲线的纵轴普遍用 a/λ 或 ωa/(2πc)两者等价。横轴不是简单的 0 到 1而是布里渊区路径上的累计倒空间距离。这两个细节是新手最常见的翻车点。如果纵轴直接画 ω/2πc 而不乘 a那么不同晶格常数的结果无法对比如果横轴用分数坐标直接画能带的斜率会被压缩或拉伸看起来形状和文献对不上。换算逻辑很简单求解器给出的是 (ω/c)²开方后得到角频率最后乘以 a 再除以 2π。写成代码就是上一章的w * a / (2 * np.pi)。横向的累计距离则用每段 k 点之间的欧氏距离累加。做完这两步你的能带图才能和论文里的曲线放在同一坐标系里比较。5. PWE 常见问题与避坑清单能带断裂、伪能带与求解器选择5.1 能带曲线断裂排序与能带交叉现象画出的曲线在某个 k 点附近突然断开或上下跳变带隙看起来忽开忽关。原因PWE 对每个 k 点独立求解本征值从小到大排序。当两条能带走近甚至交叉时固定按第 n 条连线就会在交叉点附近交换 index产生一条不连续的曲线。这不是你的代码错了而是能带本身的简并结构导致的。解决对于最低几条带如果只是要判断带隙通常排序已经够用。如果要做精细的能带追踪需要用特征向量判断相邻 k 点的模式是否连续再重新绑定 band index。更简单的工程折中是加密 k 路径采样让每条带在交叉区域的跳变足够窄视觉上仍然能看出趋势。对大多数需求来说eigh 的升序排列再加交叉点附近肉眼识别就够了。5.2 伪能带与高对比度从傅里叶系数找原因现象能带图里出现一堆几乎不随 k 变化的平带卡在带隙中间导致带隙显示不出来。原因最常见的是对 ε 直接展开而不是对 1/ε 展开导致傅里叶级数在高介电突变处出现严重振荡其次是截断数太小平面波数量不足以描述结构细节。伪带在介电常数对比度超过 15:1 时尤其明显。解决先把第 2 章公式里的 1/ε 展开用上这是最有效的修正。然后增大 N_G观察伪带是否逐渐消失。如果结构是尖锐的矩形柱或六边形孔解析傅里叶系数公式可能不好用改成 FFT 采样方案前先检查采样网格是否足够密。网格太稀会在介电常数突变处引入额外振荡伪带依旧存在。5.3 eigsh 求解低频本征值的坑现象当 N_G 很大、改用 scipy.sparse.linalg.eigsh 之后算出来的本征值全是高频或者干脆报 ArpackNoConvergence。原因eigsh 默认迭代目标是模最大的一组特征值。你的矩阵本征值都是正数模最大的自然是最高的那几条带而能带计算要的是最低频带。另一个常见原因是矩阵是 numpy 稠密数组传给 eigsh 前没有转成稀疏格式性能没有任何优势。解决求低频带时使用 shift-invert 模式让 ARPACK 把目标转到靠近 0 的特征值上。一种标准写法是from scipy.sparse.linalg import eigsh w2, _ eigsh(M, k8, sigma0.0, whichLM)sigma0 表示找零附近的特征值whichLM 表示在位移后取模最大实际得到的就是原问题的最小特征值。注意这里 M 要转换成 scipy.sparse 矩阵否则的还是稠密分解速度没有提升。如果你只是算 (2N1)² 在 900 以下的矩阵直接用 eigh 更省心不要为了省计算引入迭代求解器的复杂度。5.4 单位与 k 路径写错的英文式结果现象能带曲线形状对但带隙位置和文献差好几倍或者能带看起来完全不对称。原因倒格矢里的 2π/a 写成了 1/a横轴归一化忘了乘 a或者 k 路径用了错误的分数坐标。正方晶格的路径是 Γ-X-M-Γ三角晶格是 Γ-M-K-Γ如果把正方晶格的路径套到三角晶格上高对称点对不上能带必然畸形。解决先做一个均匀介质自检。把 ε_c 设成 ε_bPWE 应该退化成平面波色散 ωc|k|最低带在 Γ-X 中间某点的值应该精确等于该点到 Γ 的距离。这一条能同时验证单位、倒格矢和 k 路径三个环节值得每次跑例程前先用一次。6. 验证能带曲线的一条捷径用 TE/TM 两个极化方向交叉验算能带曲线的验证不一定要砸钱做实验把 TE 和 TM 两个极化分别算一遍本身就是成本最低的交叉检验。光子晶体里 TE 和 TM 的带隙经常错开同一个结构TM 在某个频段开带隙TE 可能在另一个频段开或者干脆没有。如果你的 TE 和 TM 算出来的带隙完全重合大概率是组装矩阵时把两个模式写成了一个比如忘了把范数乘积改成点积。反过来如果 TE/TM 带隙差得离谱比如其中一个几乎全是平带就要回去检查傅里叶系数有没有写错。第二个常用技巧是均匀介质极限验证。把柱体介电常数改成和背景一致相当于结构消失能带应该退化成一条直线归一化频率 a/λ 等于 k 矢量的归一化长度。在 Γ-X 路径上X 点对应 a/λ0.5。这个检查只要改一个参数就能跑能一次性暴露倒格矢 2π 因子、k 路径累计距离和纵轴归一化的问题。我做三角晶格光子晶体时第一次把倒格矢夹角算错能带整体往低频歪靠这个极限一对比立刻定位到倒格矢定义。第三个技巧是频率直方图法。把所有 k 点的本征频率收集起来画直方图带隙会表现为一段没有计数的空白区间。这个方法不用看能带曲线就能快速判断带隙边界尤其适合参数扫描时自动检测带隙。代码就两行all_w np.concatenate(wbands) plt.hist(all_w, bins200)直方图的谷底对应带隙的两条带边和 band curve 交叉验证后你可以放心把带隙位置写进论文或设计报告。我自己现在每次跑新结构都按照这个顺序来先均匀介质自检再 TE/TM 对比最后用直方图确认带隙。这套流程帮我挡掉了至少三次单位换算和极化写错的翻车希望这些习惯也能帮到你。本文还有配套的精品资源点击获取