简介SBESemiconductor Bloch Equations是用于模拟二维材料中光与物质相互作用的Python仿真程序包面向光电子、凝聚态物理与量子计算领域的研究人员可计算解析能带结构、偶极矩、光激发电流及发射强度帮助理解超快光脉冲下的电子动力学行为。包内含128个文件以95个Python脚本为主要计算核心与接口模块另有26个npy参考数据文件、2个Jupyter Notebook交互示例、1个Markdown说明文档及若干配置文件压缩包整体约203MB。目前已有206人浏览学习。通过查看源码与运行示例读者可掌握从k空间能带解析到RK4/BDF求解器调参的完整模拟流程还能借助预计算结果复现论文图表非常适合作为科研入门和二次开发的基础工具。1. 二维光与物质相互作用仿真sbe 能帮你算到什么程度做超表面或者光栅设计的人大概都有过这种体验结构画好了材料参数查好了一打开商业仿真软件面对一堆模块和网格设置先懵半小时。我最早接触 sbe 这个 2D 光与物质相互作用仿真程序包就是被一个周期性纳米孔阵列的透射异常峰逼的——实验上测到某个波长有个尖锐的透射谷手算加近似模型都对不上需要把麦克斯韦方程组在二维空间老老实实解一遍。sbe 本质上就是这么一件事它把二维空间里的光与物质相互作用问题封装成几何、材料、光源、监视器几个核心对象用 Python 脚本驱动时域求解。它的价值不是功能花哨而是结构清晰、可改动、能快速验证想法。适合做超表面、衍射光栅、波导耦合这类器件的人也适合教学场景里让学生直观看见光在介质里到底经历了什么。2. 从麦克斯韦到数值网格sbe 的物理内核与选型逻辑2.1 为什么二维仿真就够用TE 与 TM 的解耦先说说 sbe 为什么只做二维。很多实际器件在第三个方向上是均匀的比如一段横截面不变的波导、一块无限大的多层膜、一个沿 z 方向完全相同的柱状纳米结构。这种情况下如果结构沿 z 方向没有变化电场和磁场在 z 方向上就没有偏导数麦克斯韦方程组会自动解耦成两组独立方程一组是 TE 模式电场只有 z 方向分量磁场在 xy 平面内旋转另一组是 TM 模式磁场只有 z 方向分量电场在 xy 平面内旋转。这两组方程各自独立演化互不干扰。我一开始用 sbe 的时候也犯过嘀咕会不会把三维效应丢掉了后来想明白你拿一个二维截面去算本质上是假设入射光波矢在 xy 平面内且结构沿 z 向无限延伸。很多实际场景恰好满足这个前提比如平面波斜入射到一维光栅、柱状纳米线的散射、多层膜反射率的计算。sbe 里你需要先选偏振态再选对应的波源类型这个顺序不要反。如果你要仿真的结构在 z 方向有明显的模式耦合或端面反射二维结果只能作为定性参考得老老实实上三维求解那不是 sbe 的适用范围。提示sbe 的二维求解默认结构沿 z 方向无限延伸。仿真周期结构时要用 Bloch 边界条件配合布洛赫波矢否则边界上会出现虚假反射。2.2 用 Python 做光与物质相互作用仿真生态优势与性能边界市面上不缺电磁仿真工具但 sbe 选择 Python 作为宿主语言是有道理的。你想想看一个光与物质相互作用的仿真流程通常是设计结构 → 扫描参数 → 算场分布 → 提取透射反射 → 画图看趋势。这个流程里几何构建、参数扫描、数据可视化、和机器学习优化算法对接全部由 Python 生态覆盖。numpy 负责数组运算matplotlib 负责看场图scipy 负责拟合和积分sbe 只需要把麦克斯韦方程的求解内核暴露成几个类和方法就行。性能方面要泼一盆冷水。二维 FDTD 的网格通常是百万量级时间步几千步Python 本身的循环效率撑不住所以 sbe 的核心时间步进循环是用编译型语言写的Python 层面只做封装和配置。这也是为什么你要先理解它的数据结构——如果你在 Python 层写个循环去逐点更新电场速度会慢到让你怀疑人生。sbe 的边界就在这里二维高效、三维不碰。用熟了之后你会发现它的定位更像一个「可编程的实验台」而不是一个万能的商业套件。2.3 Yee 网格与时间步进sbe 的离散化原理与 CFL 条件sbe 的时间域求解基于 FDTD 方法。FDTD 的祖师爷思想是把空间切成 Yee 网格——电场和磁场分量在空间上交错半个网格步长放置这样每个电场分量周围四个磁场分量每个磁场分量周围四个电场分量空间离散后麦克斯韦旋度方程变成一堆显式差分方程。时间上电场和磁场也交错半步推进所以电场更新需要上一时刻的磁场磁场更新需要上一时刻的电场就这样一步步把光场在结构里推着走。空间步长和时间步长不是随便取的这有物理约束。FDTD 有一个稳定性条件叫 CFL 条件光在一个时间步内走过的距离不能超过最小网格步长的某个比例。二维情况下CFL 条件可以写成 c·Δt ≤ 1 / √(1/Δx² 1/Δy²)实际使用时一般取一个安全系数。sbe 的初始化必然涉及这一步import numpy as np # 仿真区域单位用微米 x_min, x_max -2.0, 2.0 y_min, y_max -2.0, 2.0 res 120 # 每微米网格数这里约 8.3 nm 一个格子 nx int((x_max - x_min) * res) ny int((y_max - y_min) * res) dx (x_max - x_min) / nx dy (y_max - y_min) / ny # CFL 稳定性条件2D 情况下取安全系数 0.9 cfl 0.9 dt cfl / np.sqrt(1.0 / dx**2 1.0 / dy**2) print(f空间步长: dx{dx:.4e} um, dy{dy:.4e} um) print(f时间步长: dt{dt:.4e} um/c)这段代码里res决定了网格分辨率的直接表现——res120意味着每微米 120 个网格点结构特征尺寸至少要大于 45 个网格否则几何会被锯齿化。dt的推导来自 CFL 条件我这里取cfl0.9留 10% 的余量防止数值振荡。注意单位是微米光速在该单位制下取 299.792458这个细节在设置波源脉冲宽度时会用到后面讲。3. 把几何写入 sbe几何构建、材料色散与边界条件3.1 从设计图到网格几何对象与布尔运算sbe 的几何构建思路和很多仿真包类似先用简单几何体矩形、圆、多边形拼出结构再映射到网格上得到每个网格点的介质编号。我第一次用的时候以为要自己逐点赋值后来发现它有Box、Cylinder、Layer这些对象还支持布尔操作。比如做一根矩形波导加一个圆形腔可以把矩形和圆相加也可以做差集挖空。import sbe import numpy as np # 创建仿真对象 sim sbe.Simulation(domain[-2.5, 2.5, -2.5, 2.5], resolution120) # 定义基底一层厚度为 0.5 um 的介质 substrate sbe.geometry.Layer( y_range(-2.5, 0.0), materialSiN ) sim.add_geometry(substrate) # 定义波导矩形截面宽 0.6 um放在基片上方 wg sbe.geometry.Box( x_range(-0.3, 0.3), y_range(0.0, 0.22), materialSi ) sim.add_geometry(wg) # 挖一个空气孔半径为 0.12 um 的圆从波导中心穿过 hole sbe.geometry.Cylinder( center(0.0, 0.11), radius0.12, materialAir ) sim.add_geometry(hole, operationsubtract)Layer适合定义大面积的平面结构比如基底、包层。Box的x_range和y_range直接给坐标区间单位是微米。Cylinder的center参数是圆心坐标radius不用多说。operationsubtract表示从已有结构里挖掉这个圆柱体——这里模拟的是波导上挖空气孔。实际布局中用得最多的是add和subtract偶尔会用到intersect做重叠区域的材料替换。网格映射是自动完成的但分辨率一定要匹配结构尺寸。比如上面的空气孔半径只有 0.12 μm如果resolution60一个网格是 16.7 nm孔直径约 14 个网格还能看如果resolution30孔直径只剩 7 个网格边缘就会严重锯齿化等效半径产生偏移共振峰位置会飘。这是最常见的「结构对了但结果不对」的根源之一。3.2 材料模型介电常数、洛伦兹色散与损耗光与物质相互作用的仿真里材料模型决定了你看得到什么物理。sbe 内置了几种材料模型最基础的是常数介电常数适合工作波段窄、色散不明显的场景。稍微讲究一点用洛伦兹色散模型描述介质在某个频率附近的响应公式是 ε(ω) ε∞ ωp² / (ω0² - ω² - iγω)其中 ω0 是共振频率γ 是阻尼系数ωp 是等离子体频率。如果你要仿真金属纳米结构的等离激元共振这类色散模型是必须的因为金的介电常数在可见光波段变化剧烈用常数近似会直接导致共振峰位置错误。# 定义带色散的材料参数按常见介质近似给出 silver_model sbe.material.LorentzModel( eps_inf3.7, omega_p9.2, # 等离子体频率单位 rad/fs omega_00.0, # 无束缚共振Drude 型 gamma0.021 # 阻尼系数单位 rad/fs ) sim.set_material(Ag, silver_model) # 也可以读取已有材料库 sio2 sbe.material.get(SiO2) # 内置材料库eps_inf是高频介电常数对应的是远高于共振频率时的背景介电响应。omega_p和omega_0的单位是 rad/fs配合微米单位制使用时换算关系是 1 rad/fs 约等于 0.159 THz。gamma对应损耗取值越小共振越尖锐但也越容易在时域仿真中产生长时间的拖尾——这会直接拖慢收敛后面避坑章节会展开。金属材料用 Drude 模型时omega_00因为金属的自由电子没有恢复力介质材料则需要保留omega_0来描述束缚电子的共振吸收。3.3 边界条件PML 层与 Bloch 周期的选择边界条件是仿真成败的分水岭。sbe 支持三种常用边界PEC理想电导体、PML完美匹配层、Bloch 周期边界。PEC 只在特殊场景用比如全金属波导的截止特性分析PML 是开放区域的标配吸收向外传播的电磁波Bloch 周期边界用于周期性结构比如光栅或超表面阵列。PML 不是简单的吸收层它要求材料参数在层内渐变让波进入时不会产生反射。sbe 里 PML 的厚度和吸收剖面可以调常见的做法是sim.set_boundary(x_min, pml, thickness0.5, profilepolynomial, order3) sim.set_boundary(x_max, pml, thickness0.5, profilepolynomial, order3) sim.set_boundary(y_min, pml, thickness0.8, profilepolynomial, order3) sim.set_boundary(y_max, pml, thickness0.8, profilepolynomial, order3)thickness的单位是微米profilepolynomial表示吸收系数按多项式渐变order3对应三次多项式。层越厚吸收越好但计算量也越大。实际经验是PML 厚度至少要是最大波长的四分之一且order取 2 到 4 之间比较稳。order1吸收太慢边界反射明显order太高则层内参数梯度太陡数值反射反而增大——这是 FDTD 里一个非常典型的「过犹不及」现象。对周期性结构x 方向用 Bloch 周期边界可以只算一个周期单元来等效无限阵列sim.set_periodic(x, period0.8, bloch_vector[0.0, 0.0])bloch_vector是布洛赫波矢的 xy 分量平面波垂直入射时设为[0, 0]斜入射时需要按入射角和波长算。这里有个隐蔽的坑Bloch 边界要求网格在一个周期内严格对齐也就是说周期长度除以网格步长必须接近整数否则边界处会不连续。4. 让光跑起来波源配置、运行流程与数据采集4.1 波源类型平面波、高斯光束与点偶极子sbe 的波源接口支持三种常用激励平面波、高斯光束、点偶极子。三者对应不同的物理问题。平面波适合计算透射反射谱、能带结构高斯光束适合模拟聚焦光与结构的相互作用或空间分辨的响应点偶极子适合研究单粒子散射或自发辐射增强。# 平面波从 y 正方向入射沿 y 负方向传播 src sbe.source.PlaneWave( center_frequency350.0, # 单位 THz bandwidth150.0, # 脉冲频谱宽度单位 THz polarizationTM, # TM磁场沿 z 方向 amplitude1.0, delay50.0, # 脉冲延迟单位 fs direction[0.0, -1.0, 0.0] ) sim.add_source(src)center_frequency350THz 对应的真空波长约 857 nmbandwidth150THz 意味着脉冲在频域上覆盖约 300450 THz 范围相当于 6671000 nm 波段。脉冲越短频谱越宽但时域分辨率要求越高。delay50fs 是脉冲峰值到达源平面的时间给一个足够大的值让电场从零缓缓升起避免瞬间加源产生的数值冲击。polarizationTM和前面说的模式选择对应TM 模式下磁场只有 z 分量电场在 xy 平面内适合算薄膜吸收和光栅衍射。4.2 场监视器与数据采集光场算出来之后得知道某个位置的场随时间怎么演化。sbe 里的监视器分两类时间监视器和频域监视器。时间监视器记录某个点或某条线上的时域场值频域监视器在时域推进的同时做离散傅里叶变换直接得到频域复振幅。# 在结构前方 0.3 um 处放一条线监视器记录透过场的频域响应 monitor sbe.monitor.FrequencyDomainField( center(0.0, 0.8), size(1.2, 0.0), # 一条水平线 freq_range(280, 420), # THz field_componentEz ) sim.add_monitor(monitor) # 再放一个点监视器记录时域电场用于观察脉冲是否跑完 time_mon sbe.monitor.TimeDomainField( position(0.0, 1.2), field_componentEz ) sim.add_monitor(time_mon)FrequencyDomainField返回复数场数据包含了相位信息——这个在做近场到远场变换时需要用到。freq_range要和波源的频谱范围匹配超出波源带宽的频率点幅值会很小数值噪声占比高。TimeDomainField的作用是判断仿真是否该停了如果时域波形已经衰减到接近零说明仿真收敛可以结束如果振幅还在振荡那就得继续跑。4.3 运行流程与日志解读一切配置完毕开始推进sim.run(total_time1200.0) # 单位 fstotal_time的设置大有讲究。一个简单法则是入射光穿越整个仿真区域所需时间加上结构内部共振衰减时间再乘一个安全系数。比如区域高度 4 μm光在其中来回一次约 26.7 fs看起来 1200 fs 很充足但如果结构里有高 Q 值的共振模式能量在腔内反复振荡衰减时间可能达到几百甚至上千飞秒。这时只看时域波形是否衰减还不够要看频域监视器的输出是否随运行时间增加而稳定。跑完后的数据取用通常是这样transmission sim.get_transmission(monitor_id0, directiony) reflectance sim.get_reflectance(monitor_id1, directiony-)这两个方法返回的是频率数组透射率数组的元组可以直接拿来画谱线。sbe 内部会做归一化用同样的仿真区域但没有任何结构时的透过场作为参考透射率就是测得的透过场除以参考场再取模平方。这也是为什么你必须在加结构前先跑一次空场基准或者 sbe 内部自动保存基准。我开始的时候没有跑基准直接看绝对值谱线觉得奇怪后来意识到归一化对象搞错了。5. 避坑手册sbe 使用中的五个高频问题与排查5.1 装了包却起不来numpy 版本与编译产物冲突现象import sbe直接报错提示某个动态库加载失败或者出现undefined symbol之类的信息。原因sbe 的核心时间步进部分是用 C 扩展编译的编译时链接的 numpy C API 版本和你当前环境里的 numpy 不一致运行时会因为符号表对不上而崩溃。另一个常见诱因是跨 Python 版本安装——比如在 Python 3.11 环境下安装了一个为 3.9 编译的 wheel。解决先用pip卸载然后升级到某个稳定的组合。我自己的习惯是先建一个独立的虚拟环境指定 Python 3.10再安装 sbe 依赖固定版本python -m venv sbe_env source sbe_env/bin/activate pip install numpy1.21,2.0 scipy matplotlib pip install sbe如果下载的是源码包而不是预编译版那就需要本机有可用的 C 编译器。Linux 下缺 gcc 会直接编译失败报错信息会停在building sbe.core extension这一行后面跟一串找不到头文件的提示。装上python3-dev和gcc后重新编译即可。5.2 结构看着没问题结果却是全反射几何没被网格采样到现象明明在结构里加了一个细小的缝隙或薄层仿真结果却完全没有体现出它的存在透射谱看起来和没有缝隙时一模一样。原因特征尺寸小于网格步长时Yee 网格采样不到这个特征。比如resolution100也就是每个网格 10 nm一个宽度 5 nm 的缝隙在映射时被直接抹平了。更隐蔽的是几何映射算法在判断网格点归属时存在歧义薄层落在两种材料编号的交界缝里。解决先量一下最小特征尺寸和网格步长的比值至少要有 35 个网格支撑一个特征。如果特征确实太小且不可放宽就得提高分辨率并接受计算量增加。验证方法很简单把结构画出来看一眼sbe 里有sim.visualize_geometry()这类方法输出网格后的材料分布图肉眼看一遍就知道结构有没有被正确采样这一步别省。5.3 加 PML 后仍有反射层厚不够或吸收剖面参数不对现象场监视器记录的近场分布中边界附近出现明显的驻波条纹或者透射谱上叠加了周期性振荡的虚假纹波。原因PML 的厚度和吸收剖面参数不匹配。当入射角较大、PML 厚度不足或吸收系数增长过快时吸收层本身就像一层阻抗不匹配的介质产生反射。另外PML 参数是按垂直入射调校的掠入射时性能会下降。解决PML 厚度增加到最大波长的三分之一到一半order从 3 往下调到 2同时检查入射角范围。斜入射场景下需要更厚的 PML这是经验法则没有捷径。调试时可以临时关掉结构只放均匀介质加 PML跑一遍看边界反射到底来自 PML 还是来自结构散射这一步能快速定位问题。5.4 透射谱振荡不收敛脉冲没跑完或监视器频段超范围现象同一组参数把total_time从 500 fs 增加到 2000 fs透射谱却在某个波长附近明显变化或者谱线首尾出现极不合理的尖峰。原因时域脉冲还没完全穿过仿真区域频域监视器的傅里叶变换是被截断的等效于在时域上加了矩形窗频谱上就出现旁瓣振荡。另一种可能是监视器的freq_range覆盖了波源带宽之外的频率那些频率点的入射功率几乎为零透射率计算成了 0/0 的不稳定量。解决观察时域监视器的波形是否归零再结束仿真。如果波形有长尾说明结构里有高 Q 谐振需要用更长的total_time。对于频段问题把freq_range缩小到波源中心频率加减约 0.8 倍带宽的范围宁可少看一些频点也别让数值噪声污染谱线。每次调完参数后把监视器数据导出并重复检查几个特征频率确认幅值不再波动。5.5 内存暴涨或跑得很慢二维问题不该吃掉几十 GB现象网格规模看起来不大但运行几秒后内存占用飙升到几十 GB程序开始疯狂换页计算速度直线下降。原因常见的有三种。一是监视器数量过多且每个都存储整个频率范围的三维场快照数据量爆炸二是时间步推进的中间数组缓存没有复用每次循环都在重新分配内存三是在 Python 层对场数组做了循环操作触发了大量的临时对象生成。解决检查监视器去掉不需要的频域场监视器改用点监视器和线监视器。如果必须要保存全场数据降低保存频率或只在特定时间窗口保存。Python 侧不要自己写场更新的循环sbe 的全部数组运算都应该走向量化接口——你写for i in range(nx)逐点遍历的那一刻性能就已经毁了。必要时用sbe.enable_memory_optimization(True)开启内存复用。6. 从验证到进阶用解析解校核 sbe 结果与近场-远场变换6.1 用 Fabry-Perot 公式做基准验证拿到一个新工具第一件事不是直接算复杂结构而是找一个有解析解的场景验证它。最简单可靠的验证是用单层介质平板——厚度 d、折射率 n垂直入射时反射率和透射率可以用 Fabry-Perot 公式闭式表达。设空气-介质-空气三层结构反射系数 r 和透射系数 t 分别是r r₁₂ t₁₂ t₂₁ r₂₁ e^(i2δ) / (1 - r₂₁² e^(i2δ))其中 δ 2π n d / λ 是单程相位r₁₂、t₁₂ 是界面上的菲涅尔系数。这个公式我每次验证新环境都会跑一遍在 sbe 里建一个厚度 0.5 μm、折射率 1.5 的介质层用平面波垂直入射波长从 600 nm 扫到 1200 nm把透射率谱线和 Fabry-Perot 解析解叠加对比偏差通常在 1% 以内。如果偏差明显先怀疑网格分辨率再怀疑材料色散参数和归一化设置。这个验证习惯最大的价值是它把 sbe 的黑匣子拆开了一条缝。当后续算复杂超表面遇到反直觉的结果时你能确定问题出在结构设计上而不是数值内核上。我的经验是任何仿真工具都不值得盲目信任但一个通过了解析验证的工具至少能让你把精力集中在物理问题上。6.2 近场-远场变换与远场散射截面提取sbe 的频域监视器给出的是近场复振幅但很多实际需求是远场散射方向图比如纳米粒子散射到各个方向的光强分布。近场到远场变换的经典做法是基于等效原理在一个包围结构的外边界面上记录切向电场和磁场分量然后用 Stratton-Chu 积分把边界上的场源等效到远场。sbe 里这个变换不是自动完成的——至少需要你自己准备一个包围结构的闭合边界并导出边界上每个点的切向分量。实际实现大致是这个形态# 假设 monitor 已记录包围盒边界上的切向场 boundary_data sim.get_boundary_fields(monitor_idbox_surface) far_field sbe.postprocess.near_to_far( boundary_data, observation_anglesnp.linspace(0, 2*np.pi, 360), wavelength800.0, n_background1.0 ) # far_field 是 (角度, 复振幅) 的数组代码里的observation_angles是你要计算远场的角度数组单位弧度wavelength800.0对应垂直入射时的中心波长单位纳米。变换结果通常有两种输出远场电场幅值和远场强度。散射截面可以由远场强度对全空间积分再除以入射波强度得到。用这个接口时要注意包围盒必须完整包住所有非线性材料和结构,并且边界面上的场监视器分辨率不能太低——边界采样稀疏会导致积分不准确远场方向图出现锯齿状伪影。如果是规则形状的散射体比如圆形介质柱还可以用米氏散射的解析解交叉验证近远场变换的中间结果。对于金属纳米粒子先算吸收截面、散射截面再比较两者之和与消光截面是否守恒误差在几个百分点以内就说明流程可靠。这个守恒关系是一块很好的试金石,一旦发现消光不等于吸收加散射优先回头检查场监视器是否遗漏了损耗区域的内部场分布。从那以后我每次换机器、换环境、换版本拿到 sbe 的第一件事就是跑一遍单层平板的 Fabry-Perot 验证再跑一个圆形柱子的米氏散射对照。这两步跑通了我才敢把时间花在复杂结构上。这套习惯帮我挡掉了至少三次「仿真结果和实验对不上」的危机——其中一次是 PML 参数在更新版本后的默认值变动另一次是材料库的色散模型被替换成了更精确的版本导致共振峰移动。希望帮到你。本文还有配套的精品资源点击获取
