简介面向光子晶体与拓扑光学研究者的资源包聚焦陈数计算与仿真。内容涵盖光子晶体能带结构计算原理、陈数概念解析以及利用软件结合脚本进行仿真的完整方案。压缩包共五个文件包含两篇讲解文档分别介绍光子晶体建模原理与软件操作要点一个仿真模型可直接打开查看周期结构与边界条件设置一个脚本用于计算数据后处理并绘制布里渊区能带图另有文本文件作为补充说明整体大小约二十兆。目前已有1043人学习下载。资源内不仅提供理论讲解还给出可运行的仿真模型与脚本读者可结合交互流程从周期性结构建模、麦克斯韦方程频域求解到能带图绘制与陈数分析形成完整实践闭环。适合需要从零上手光子晶体数值计算的研究生与工程师。1. 陈数计算不神秘为什么用Comsol做拓扑不变量不少人第一次看到“陈数”这个词都会愣住然后问一句陈数怎么写是陈省身那个“陈”吗没错陈数Chern number以微分几何大师陈省身命名是个整数拓扑不变量在量子霍尔效应、拓扑绝缘体、拓扑光子学里反复出现。我这两年一直在用Comsol做电磁场仿真后来发现其实能把陈数计算也塞进Comsol工作流里把光子晶体的本征场导出来再用一小段Python算Berry曲率就能得到一堆能带各自的陈数。这个方法的好处是只用Comsol的波动光学模块和参数扫描不需要另建紧束缚模型适合验证连续介质里的拓扑性质。这篇文章就按我实际跑通的路径来写从物理公式到Comsol设置再到差分包络代码和那些折磨过我的坑。2. 先把物理模型写对从陈省身示性类到Berry曲率2.1 陈数的积分表达你真正要算的是哪个公式陈数定义在二维布里渊区上。对第n条能带公式写成C_n (1/2π) ∬ B_n(k) dk_x dk_y其中B_n(k)是Berry曲率由布洛赫函数u_{n,k}(r)的相位变化决定。严格写是B_n i (∂u*/∂k_x · ∂u/∂k_y - ∂u*/∂k_y · ∂u/∂k_x)但这需要波函数对k的导数数值上非常不稳。实际工程里没人直接求导而是用Fukui-Hatsugai-Suzuki离散化方法在k空间网格上用相邻格点波函数的内积来构造一个“链接变量”再取对数得到Berry相位最终求和得到整数陈数。这个办法的好处是规范无关不用管波函数整体相位怎么变。所以你要做的第一步不是写表达式而是先把“布洛赫波函数”从仿真里拿出来。光子晶体里布洛赫波是周期边界条件下的本征电场模式E_k(r)它满足E_k(r) u_k(r)e^{ik·r}。Comsol的Floquet周期条件可以直接把波向量k加进去求出来的本征场就是u_k(r)。这正好是上面公式里需要的对象。2.2 在Comsol中搭一个可扫参的哈密顿量以连续狄拉克模型为例在Comsol里算陈数需要先在二维布里渊区上得到一个完整的波函数网格。我用的是方形光子晶体晶胞边长a1 μm中央一个半径r0.3 μm的介质圆柱折射率3背景折射率1。用波动光学模块里的“电磁波频域”接口做本征频率求解。关键是周期边界条件。在“周期”边界条件里选择“Floquet周期”点对点绑定并定义“Floquet波向量”的k_x和k_y为两个全局参数。具体设置里我会把k_x kx2pi/ak_y ky2pi/a其中kx和ky是扫描参数单位是1。这样布里渊区范围就归一化成[-0.5, 0.5]取倒空间边界扫描间隔取0.05就得到21×21个点。注意这里算的是二维TE模式电场沿z方向所以本征方程是标量亥姆霍兹方程。用两个因变量不需要TE模式下只有一个电场分量解出来直接就是复数场。这比算矢量模式简单得多而且陈数计算只依赖一个分量的复相位。提示如果你用波段稍宽的软件包可能担心折射率色散。这里仿真是单频本征求解没有色散折射率直接写常数即可。3. 让Comsol输出波函数特征值求解与扫参设置3.1 参数化扫描k点最小可复现的研究设置在“研究”里添加“特征值”求解然后在旁边嵌套一个“参数扫描”。扫描参数选kx和ky从-0.5到0.5步长0.05。这样Comsol会循环求解本征势阱内的电磁模式每个(kx,ky)下返回一组本征频率和对应的场分布。我通常只关心第一条导带或者第一条在带隙内的能带。为了不让模式顺序乱跳需要在扫描里设置“特征值搜索基准”。对于二维晶格我一般搜eigenfrequency附近的值比如5 THz到10 THz之间。更稳的做法是先用一个固定点比如kxky0扫一次本征频谱拿到目标带的频率然后把扫描的初始区间缩窄到该频率附近。不然有些偏远的k点会解到别的带上去导出的场就串频了。参数扫描完成后把“每步求解”打开这样每个k点的数据都会保留在解集里。随后在“结果”里找到“数据集”它会生成一个带“参数维度”的求解数据里面索引就是(kx,ky)。这一步成功的话数据集节点下面能看到21×21条子解。3.2 导出波函数到文本数据格式与完整路径Comsol不能直接给你“每个k点一个文件”但可以在“结果”里用“求值组”把场值一次性导出。我用的方法右键“结果”选择“导出”添加“数据”选择上面的数据集然后设置“表达式”为电场z分量在频率域求解里通常叫ewfd.Ez复数值。输出格式选“文本”文件名写Ez_scan.txt。最关键的是在“文件类型”里选“x坐标、y坐标、Ez实部、Ez虚部”并且勾选“所有网格数据点”。这样导出的文本格式大致是x y z real(Ez) imag(Ez)... 每行对应网格上的一个点按点编号顺序排列但不同k点的数据会堆在一起。为了避免文件过大我一般把网格“细化”一次就行然后用“导入”到Python。这里有个坑导出的数据顺序是每个参数点一个块还是所有点混合不同Comsol版本有差别。老版本会重复输出“x y z”坐标只有场值属于当前参数步新6.x版本则会为每个参数点输出一遍完整坐标。我用了一个通用办法在导出文件里每个参数点的数据块行数等于网格自由节点数。所以只要知道网格节点数就能分块读取。或者在导出时把“参数值”也加进去这样每行都有(kx,ky)列。推荐后一种因为可以避免顺序假设。注意如果你只导出“被选择”的边界或域节点数会变分块时容易错位。所以我总是在“数据集”里选“完整数据集”导出全部体元素。4. 用Python把陈数算出来网格差分与数值积分4.1 Berry曲率的实空间差分公式拿到所有k点上的复电场E_k(x,y)后我要在每个k点计算四个相邻点之间的内积。定义链接变量U_x(k) ⟨E(k) | E(k kx̂)⟩ / |⟨E(k) | E(k kx̂)⟩|内积定义为对空间网格所有节点求和⟨E(k)|E(q)⟩ Σ_j conj(E_k(x_j,y_j)) * E_q(x_j,y_j)。注意这里只求和场点乘不需要再乘面积因子因为每个网格点代表同一个有限元积分权重。实际上严格有限元内积要乘质量矩阵但均匀网格近似下直接求和已经够用要精确的话可以在Comsol里用积分算子得到内积但导出后手动求和也能得到正确整数只是收敛稍慢。然后Berry曲率的离散化公式是F(k) log[ U_x(k) * U_y(k kx̂) / ( U_y(k) * U_x(k kŷ) ) ]注意这里的除法是复数除法。log取主值但最后对整条带累加时需要把每个F(k)的虚部累加再除以2π。累加过程中可能因为相位缠绕跳出主值区间所以要用“相位展开”技巧把当前累加值和下一个F(k)的虚部之差调整到(-π, π)内再累加。这个技巧当晚辈写代码时很容易漏。4.2 陈数积分与开路校验代码与参数下面这段Python代码可以直接用假设你已经从Comsol导出了数据并按参数点维度重构成了一个4维数组E[nky, nkx, Npts, 2]最后一维是实部和虚部。重构代码不贴了因为不同导出格式略有不同核心是计算环节import numpy as np def chern_number(E, nkx, nky): E: shape (nky, nkx, Npts) complex field values assumes k grids are uniform in [-0.5, 0.5] F np.zeros((nky, nkx), dtypecomplex) for iy in range(nky): for ix in range(nkx): # 四个角点索引周期性地包裹 k0 E[iy, ix] kx_plus E[iy, (ix1)%nkx] ky_plus E[(iy1)%nky, ix] kx_plus_ky_plus E[(iy1)%nky, (ix1)%nkx] # 计算链接变量 def U(a, b): val np.vdot(a, b) # np.vdot 取共轭求和 return val / np.abs(val) if np.abs(val) 1e-12 else 1.0 Ux1 U(k0, kx_plus) Uy2 U(kx_plus, kx_plus_ky_plus) Uy1 U(k0, ky_plus) Ux2 U(ky_plus, kx_plus_ky_plus) F[iy, ix] np.log( (Ux1 * Uy2) / (Uy1 * Ux2) ) # 相位展开累加 sum_F 0.0 for iy in range(nky): for ix in range(nkx): # 调整到 (-pi, pi) d np.angle(F[iy, ix]) sum_F d chern sum_F / (2.0 * np.pi) return chern这段代码里np.vdot(a,b)会自动共轭第一个参数再点乘正好符合内积定义。np.angle返回相位在[-π, π)如果直接累加遇到“相位跳变”会把陈数算成0。简单办法是检查相邻累加值的差如果差值超过π则加上或减去2π调整。我没有在代码里写因为实际调试时才发现留到避坑段落说。运行这个函数前先校验一个平凡情形把相邻k点的波函数设为同一个常数场此时所有内积模|⟨a|b⟩|1F0陈数自然为0程序不会报错。再跑一个已知的量子霍尔等价模型陈数应为1或-1。我用一个解析狄拉克模型验证过得到1。提示如果网格是21×21陈数收敛一般到±0.999左右不会恰好精确到1但足够确认拓扑数。继续加密到41×41会进一步逼近整数。5. 避坑算陈数最容易翻车的5个点5.1 现象陈数算出0.7或者1.3不是整数原因波函数整体相位在相邻k点间随机翻转。Comsol每次求解的相位原点由求解器决定不同k点之间没有关联但你计算的链接变量是规范不变的理论上不受影响。实际翻车是因为数值误差导致某个内积的模在1附近而分母接近0引入噪声。规范不变只是数学上不变离散化后大相位差会引起精度损失。解决先检查每个k点的场模与相邻点场模的重叠度。如果某个U的模偏离1超过1e-6说明网格分辨率不足或模式翻转。要保证在扫描布里渊区时相邻k点之间的场分布没有发生模式突变。把扫描步长从0.05减到0.02再看陈数是否稳定。5.2 现象简并点附近的Berry曲率剧烈跳动原因在布里渊区某些点两条能带交叉例如狄拉克点本征场在这些点附近非解析差分失效。特别是扫描网格刚好处在交叉位置F计算随机。解决在Comsol参数扫描里设置“跳过特定k点”不行Comsol没法自动跳过。我用过的方法是避开简并点把k网格稍微偏移半格例如取kx-0.475, -0.425 ... 而并非从-0.5开始。这样网格点不会正好落在简并点上。另外确定你算的是同一根能带可以观察频率在整个扫描里连续变化而无跳变。5.3 现象相位累加时陈数莫名其妙多一整个整数原因相位缠绕。离散F(k)的虚部可能因为取对数的主值在某个k点从π跳到-π累积时少加了2π最后陈数多1或减1。解决在累加前做一个相位展开例如用np.unwrap处理一维数组但二维展开麻烦。我写了一个简单循环每次拿当前累积和与F(k)的虚部差如果差大于π就减2π小于-π就加2π。这相当于强行把相邻(point)的相位差约束在(-π,π)。注意我们是对所有k点求和不是逐点独立所以需要按扫描顺序展开。5.4 现象导出文件太大内存被击穿原因21×21个k点每个网格上几千个节点导出全文是几百MB。直接全读进来会把16GB内存榨干。解决不要一次导出所有k点到一个文件。改成在Comsol里用“导出”里的“参数扫描”选项每个参数值一个文件文件名自动加索引例如Ez_scan1.txt,Ez_scan2.txt。Python端循环读文件逐个构造E数组。我用的块大小为每步一个文件内存峰值降了一个量级。5.5 现象布里渊区边界处链接变量计算错误原因当kx0.5和kx-0.5周期连接时Floquet波向量相差一个倒格矢波函数相位应该一致。但Comsol的Floquet条件在边界处可能引入了一个额外的相位因子e^{ikx * L}导致导出的场在周期方向上不完全周期。解决使用“自动周期”配对时Comsol会要求设置“点对点周期条件”但导出的u_k(r)是已经扣除了平面波的周期函数。大多数情况下没错但检查一下对同一个k在x最小和最大两个边界的场值应满足周期关系。如果不对需要用“比邻”周期边界条件并勾选“Floquet周期”并确认“波向量”表达式中k_x和k_y的单位与几何单位一致——因为几何尺度是微米波向量要用2π/a归一化否则边界条件求解的频率会发生折叠。6. 进阶扫描晶格参数验证拓扑相变6.1 用参数扫描观察陈数从0到1当你把波函数导出来以后其实很容易做“参数学习”改变介质柱半径r或折射率重新扫描整个布里渊区就能观察到陈数变化。比如在方形光子晶体里介质柱从细到粗系统会经历一次带隙闭合再打开陈数从0跳到1。做法是把Comsol“参数扫描”从原来的(kx,ky)扩展为(r,kx,ky)的三重扫描。计算量会大增我一般先固定r单独算一次陈数然后对感兴趣的r值逐点测试。如果想在一条命令里批量算可以给Comsol加Livelink for MATLAB或者用Python控制Comsol的Java API。不过直接改“参数扫描”列表是最简单的方式。唯一要注意的是每次改变r都必须重新划分网格因为圆柱边界移动了。在“网格”节点里勾选“每次参数扫描时重构网格”代价是扫描速度变慢但避免用旧网格去解新几何。6.2 Wilson loop验证另一种更稳的算法Fukui-Hatsugai-Suzuki方法已经够稳但它严格要求k网格均匀。如果网格不均匀或者你想检查结果Wilson loop是更通用的验证手段。Wilson loop本质是沿一条闭合路径累积Berry相位对固定的k_y沿k_x方向走过去计算每一步的链接变量乘积然后求相位。这个相位随k_y的缠绕数winding number就是陈数的投影。这段代码可以用上文的U_x和U_y重新组合def wilson_loop(E, nkx, nky): phases [] for iy in range(nky): prod 1.0 0j for ix in range(nkx): # 沿x方向渐进 Ux U(E[iy, ix], E[iy, (ix1)%nkx]) prod * Ux phase np.angle(prod) phases.append(phase) # 计算缠绕数 winding 0 for i in range(1, len(phases)): d phases[i] - phases[i-1] while d np.pi: d - 2*np.pi while d -np.pi: d 2*np.pi winding d return winding / (2*np.pi)如果Wilson loop给出的缠绕数和陈数结果一致就可以放心。我在实际项目里遇到网格较粗、Fukui方法有疑惑时用Wilson loop复核后最终都能得到同一个整数。最后留一个习惯每算完一个参数点我都会随手把k0处的模式频率和小k附近的三条能带频率打印出来确认带隙没有被意外打开或关闭。这个动作帮我抓住了两次因为几何参数输入错误导致“陈数处处为0”的诡异问题。希望这篇笔记能帮你在Comsol里把陈数算得干净利落。本文还有配套的精品资源点击获取