Comsol计算陈数：从布洛赫波函数到Fukui格点规范实现

简介这是一份面向光子晶体研究者与Comsol/MATLAB用户的陈数计算专题资料内容聚焦光子晶体能带结构模拟与拓扑陈数计算解决从Comsol建模到MATLAB数据后处理的核心问题。压缩包共5个文件内含两篇PDF说明文档、一个mph模型源文件、一段m计算脚本以及txt配套文本总体约20MBPDF侧重原理解读与文件说明mph与m代码则提供可运行的三维建模和能带提取示例。已有1043人学习下载。资料兼顾概念与实操既讲解光子晶体周期性结构、光子禁带、布里渊区划分等基础也演示利用频域仿真求解麦克斯韦方程获取电磁场分布并通过LiveLink将结果导入MATLAB绘制能带图最终用于陈数分析。对设计光子晶体光纤、光子晶体激光器等器件或从事拓扑光子学研究的读者具有直接参考价值。1. 陈数计算分享这个拓扑整数到底怎么从Comsol里“写”出来陈数Chern number在拓扑能带理论里算是最重要的一个整数二维周期结构每一条能带都带一个 C∈ZC 不为 0 意味着边界上必然出现单向输运的边界态。Comsol 用户把能带图扫得很漂亮卡就卡在“陈数的数怎么写”——翻遍后处理菜单也找不到一个叫 Berry 曲率的变量。这篇文章就把两件事讲透用 Comsol 把布洛赫波函数导出来再在后处理脚本里用格点规范场把陈数算出来同时解释为什么能带打开、带隙清晰陈数照样可能算成 0.5。适合已经在 Comsol 里跑过周期结构本征模、现在想算拓扑数和边界态的工程师、研究生。下面按我实际跑通的一条路径来写所有参数都给出可以直接改的值。2. 把陈数拆成能上机的公式先从 Berry 曲率写到格点相位求和2.1 陈数的标准写法与为什么不能直接在 Comsol 后处理里积分布洛赫态写成ψ_{n,k}(r) e^{ik·r} u_{n,k}(r)其中u是周期因子注意它跟 Comsol 直接解出的场ψ差了一个平面波因子。Berry 曲率的标准写法是Ω_n(k) i ( ⟨∂_kx u | ∂_ky u⟩ − ⟨∂_ky u | ∂_kx u⟩ )陈数就是 Berry 曲率在布里渊区上的积分C_n (1/2π) ∬_{BZ} Ω_n(k) dk_x dk_y这个式子看着简洁直接拿到 Comsol 里却写不出来原因有三个。第一Comsol 解出的布洛赫波函数相位是不固定的每次重新求解都会乘一个任意的e^{iχ(k)}直接用数值微分算∂_k u微商会被人工相位噪声淹没。第二你需要在相邻k点拿到同一物理能带上的波函数但特征值求解器默认按频率排序简并点附近很容易串带。第三有限元节点值的内积并不是简单的逐点相乘后续会专门讲。所以业内做陈数几乎都用“格点规范场”写法也就是 Fukui–Hatsugai–Suzuki 方法把布里渊区切成N_x × N_y个小格子每个格子边上定义U1(i,j) ⟨u(i,j)|u(i1,j)⟩ / |⟨u(i,j)|u(i1,j)⟩|U2(i,j) ⟨u(i,j)|u(i,j1)⟩ / |⟨u(i,j)|u(i,j1)⟩|然后陈数写成C (1/2πi) Σ_ij ln[ U1(i,j) · U2(i1,j) · conj(U1(i,j1)) · conj(U2(i,j)) ]这才是“陈数的数”在后处理代码里的标准写法四个 U 围成一个小闭合回路取整体辐角把所有小格子的辐角累加再除以2π。直接用这段公式避开了 Berry 曲率里所有的规范依赖问题。2.2 Comsol 四个必须对齐的配置Floquet 边界、复值解、扫描顺序和导出列要把上面的公式跑起来Comsol 模型里有四个配置不能错错了后面代码再漂亮也算不出整数。配置项具体做法出错后果布洛赫波矢全局定义kx、ky单位用rad/m扫描区间用[-π/a, π/a)布里渊区闭合失败Floquet 周期边界在周期边界节点下选 Floquet源/目标边界引用kx、ky只能跑 Γ 点特征值研究求解类型选“特征值”确认输出复本征值和复场波函数没有相位信息参数化扫描外层参数kx内层参数ky导出行顺序和代码重排对不上导出表达式同时导出u和imag(u)不要只导默认实部陈数恒等于 0扫描网格我一般取Nx Ny 24起步。kx -π/a (0:Nx-1) * 2π/(a*Nx)意思是最后一个点不取π/a让首尾天然闭合。如果你用linspace(-π/a,π/a,Nx)那样的写法布里渊区边界会被重复采样后面mod周期索引会出现错位。特征值求解时要多求两条能带。只求目标能带很容易搅进非物理的零模或数值模式多算两条带做模式连接时才有余量。对于光子晶体、声子晶体物理场接口本身自带复数求解如果用 PDE 模式自己写方程记得确认“允许复数本征值”没有被关掉。这些配置对齐后求解器会给你一组带kx/ky/mode/freq/u/im(u)的记录下一步就是导数据。3. 把 Comsol 波函数倒进 MATLAB从导出列到陈数积分的最小脚本3.1 导出文件的行列布局一列复数场分量怎么摆在 Comsol 的“导出”节点里把数据集指到特征值求解的解。导出的文本通常是这种结构# kx ky mode freq re_u im_u 0.0000 0.0000 1 2.1345e11 1.0243e-1 0.0000e0 0.0000 0.0000 2 3.2011e11 3.2210e-2 4.1103e-1 ... 0.2618 0.0000 1 2.1410e11 5.1110e-2 8.0030e-2同一组(kx, ky)会有mode1,2,3,...多行每一行是该波矢下的一个本征模。这里最容易犯的错是把mode列当成能带编号它只是求解器输出顺序不是物理能带顺序。物理上同一条能带要靠后面讲的模式连接来重新排。读取时按kx、ky重排成数组。以下代码假设数据列依次是 kx、ky、mode、freq、实部、虚部% 读取Comsol导出的文本跳过1行表头 raw readmatrix(chern_data.txt, NumHeaderLines, 1); kx raw(:,1); ky raw(:,2); mode raw(:,3); u raw(:,5) 1i*raw(:,6); % 拼成复数波函数 band 1; % 先取第1条物理能带 sel mode band; % 这一句只在模式已经对齐时可用 kxb kx(sel); kyb ky(sel); ub u(sel); kx_list unique(kxb); ky_list unique(kyb); Nx numel(kx_list); Ny numel(ky_list); Psi zeros(Nx, Ny, 1); % 先存单分量波函数 for q 1:numel(kxb) ix find(kx_list kxb(q), 1); iy find(ky_list kyb(q), 1); Psi(ix, iy, 1) ub(q); end如果你是弹性波、压电或多分量 PDE就把每个分量的实部/虚部都导出来在循环里拼成Psi(ix,iy,:)的一个复数向量再传入后面的函数。unique要求kx/ky网格是规则的如果扫描里混入了失败解数组重排前先检查numel(kx_list)*numel(ky_list)是否等于实际行数。3.2 Fukui 格子规范把陈数的数用四行相位凑出来重排好波函数后主菜就是一页纸的陈数函数。这个函数不计算 Berry 曲率只计算相邻k点之间的 U(1) 相位链接function C chern_from_comsol(Psi) % Psi: 复数波函数数组 size(Nx, Ny, Ndof) % 输出: 整数陈数 Nx size(Psi, 1); Ny size(Psi, 2); ix_next [(2:Nx), 1]; % kx 方向周期相邻 iy_next [(2:Ny), 1]; % ky 方向周期相邻 U1 zeros(Nx, Ny); U2 zeros(Nx, Ny); for i 1:Nx for j 1:Ny v0 Psi(i, j, :); v0 v0(:); v1 Psi(ix_next(i), j, :); v1 v1(:); v2 Psi(i, iy_next(j), :); v2 v2(:); O1 v0 * v1; % 复内积 O2 v0 * v2; U1(i, j) O1 / abs(O1); % 只保留相位丢掉模长 U2(i, j) O2 / abs(O2); end end F U1 .* U2(ix_next, :) .* conj(U1(:, iy_next)) .* conj(U2); C round( sum(angle(F), all) / (2*pi) ); end这段代码的逻辑就三句话。第一U1/U2是相邻波矢上波函数的内积再归一化它把求解器随机乘的e^{iχ}消掉了这就是为什么不需要手动固定规范。第二F把四个方向链接合成一个闭合回路angle(F)是每格 Berry 相位。第三全部小格相位累加后除以2πround只负责清掉1e-13量级的数值尾巴。调用时直接说C chern_from_comsol(Psi)。如果算出来是0.9998取整后1正常如果算出来是0.63或1.35别急着 round先去看是不是模式串了、网格不够细、或者内积近似出了问题。4. 陈数算成 0.5 的翻车现场五个让结果变成废数的坑4.1 能带在交叉处互换积分给出 0.63现象同一个模型跑了三遍陈数分别是 0.63、1.02、0.5怎么改参数都不稳定。原因Comsol 特征值求解器按频率大小返回本征模。能带交叉或反交叉时mode序号在某个k点突然跳变你取的第 1、2、3 条“带”已经不是同一个物理能带。解决在算法里做模式连接。具体做法是用上一个k点的波函数跟当前点所有候选模式做内积取模最大的候选作为同一条带的延续。模式连接做完陈数通常立刻回到整数。这也是我不建议按modeband直接筛数据的原因第 5 章会给可操作步骤。4.2 导出只带实部陈数恒等于 0现象所有U1/U2都等于 1angle(F)全为 0C 永远为 0。原因Comsol 默认后处理视图只显示实数部分导出节点如果没有显式加imag(u)你拿到的是一组纯实波函数。相位信息全丢陈数自然写不出来。解决导出表达式中单独加一列imag(u)。用 PDE 多分量时记得是imag(u1)、imag(u2)不是imag(u)一个量。检查导出文件里波形是否带明显的复数虚部。4.3 布里渊区两端点重合陈数对 N 不稳定现象Nx20时 C0.87Nx30时 C0.92算出来的数离整数越来越近但就是不落在整点上。原因扫描范围写成了[-π/a, π/a]闭区间。边界点-π/a和π/a在布里渊区里是同一个点等于是重复采样。重复点破坏周期相邻关系整个闭合回路被多算或少算了一个格点。解决扫描区间用半开区间。比如kx -pi/a (0:Nx-1) * 2*pi/(a*Nx)。代码里的ix_next[(2:Nx) 1]已经默认了周期闭合扫的时候别再把同一个点扫两次。4.4 高对称点简并相位因子乱跳现象C 在 0 和 1 之间反复横跳angle(F)在某个局部区域出现类似棋盘格的乱值。原因在 Γ、K、M 这类高对称点上两条能带简并。简并子空间里的本征矢可以任意旋转单条能带的 Berry 相位在简并点附近根本不唯一Fukui 公式在简并点上给不出稳定结果。解决把网格整体偏移半个格点避开高对称点。具体做法是kx -pi/a (i-0.5) * 2*pi/(a*Nx)让采样点落在高对称点之间的位置。如果简并不是孤立的而是整条能带简并单能带陈数本身没有定义需要做多能带的 Wilson loop那不是几行脚本能解决的问题。4.5 节点值直乘当内积网格一改符号就翻现象用一样的网格算某个模型 C1加密网格后符号变 −1甚至出现 0.5 附近的数。原因有限元的自由度是节点值严格的两个波函数内积应该是u1ᵀ M u2其中 M 是质量矩阵。直接用节点值逐点相乘等于把 M 当成单位阵。网格很细时逼近得还可以但粗网格或者波函数剧烈振荡时会出问题。解决有三个常见做法。第一用 LiveLink for MATLAB 从 Comsol 里拿质量矩阵做完整内积第二把解插值到规则采样网格再算逐点内积能减小质量矩阵带来的偏差第三先做收敛性检查在 Nx24、32、40 三种网格上各跑一次陈数如果对 N 不敏感就可以接受。陈数本身是拓扑整数抗噪声能力很强但如果出现基础积分权重错误收敛性检查是能揪出来的。5. 把结果钉死用模式连接和边界态验证陈数而不是只看整数5.1 模式连接让能带编号在扫描里不换脸拿到整数 C 之前我建议先跑一遍模式连接至少能确认你的“第 1 条带”从头到尾是同一条物理能带。做法很简单对每个(i,j)格点先拿上一个k点已经确认的波函数v_old跟当前格点所有候选模式算内积G(m) | ⟨v_old | v_m(k)⟩ |模最大的m就是这条能带的延续。把所有格点都这样扫一遍重排后的Psi再丢进chern_from_comsol。如果重排前后 C 不一致说明之前mode列表确实串带了。这个步骤的关键参数是内积阈值G_max如果小于 0.5说明当前格点已经不在同一物理带附近要么扫描步长太大要么存在简并。这时候加密kx/ky网格比怀疑公式更有效。高对称点附近建议用 0.5 格偏移的网格来处理。5.2 边界态穿越体带隙的次数是一台“验钞机”陈数是整数量子数单看一个整数也可能算错符号。我一般的验证习惯是把同一几何改成一维带状结构x 方向仍用 Floquet 条件扫kxy 方向取 20 个周期并截断成有限宽边界保留开放。算完特征频率后画f(kx)能带图。这时你会看到体带围出的带隙里出现一条或几条边界态色散。边界态从下体带穿越带隙、连到上体带的次数就是整条介质边界处的拓扑陈数。比如 C1边界态必须跨带隙一次C2要跨两次。方向对不上说明符号写反了通常回去检查U1和U2的定义顺序即可。我还习惯画出每个小格子的angle(F)热图陈数非 0 时热图上曲率会集中在几个动量点附近而不是均匀铺开看到热图就基本知道这个整数是从哪里挣来的。这套流程我从最初算光子晶体、声子晶体一路用到现在最大的教训就是陈数不能靠肉眼看出能带边缘的绕向也不能只看后处理里的某个“趋势”必须把布洛赫波函数的相位老老实实导出来用格点规范场凑成闭合回路。把模式连接和边界态验证两个习惯一起用上即便不依赖别人现成的代码包也能在 Comsol 里稳定复现别人文章里的陈数。希望帮到你。本文还有配套的精品资源点击获取