简介这套资料围绕光子晶体陈数计算与COMSOL仿真面向具备一定光学或电磁场基础的研究生、工程师解决能带结构计算及陈数后处理中的具体问题。压缩包共含5个文件包括两篇PDF文档、一个mph模型文件、一个MATLAB脚本及一份txt说明整体约20MB可覆盖从理论原理到软件操作的完整链条。截至目前已有1042人学习浏览适合正在做光子晶体拓扑性质研究的读者参考。内容兼顾陈数概念讲解、COMSOL建模思路与MATLAB数据处理还涉及布里渊区分析和能带图的绘制其中mph模型可直观查看仿真设置m脚本方便修改参数并复用配合说明文档能帮助读者快速完成从电磁场模拟到陈数计算的实践显著缩短自行摸索周期。1. 在 Comsol 里算陈数到底要算哪个数很多做拓扑光子学或声子晶体的人把能带图跑出来后卡在同一句话上陈数的数怎么写。有人到全局定义里建了一个chern填0或者1发现什么都不变有人把某一支本征态的实部导出来做面积分得出一个连自己都不信的 0.03。问题不在 Comsol 用得不熟而是陈数根本就不是由单个波函数算出来的。它是布里渊区上 Berry 曲率的积分Comsol 在这件事里只负责提供每个 k 点的本征态真正的“数”要在 k 空间网格上做离散求和。本文就从公式写法、模型搭建、数值积分到结果验证把这条链路走通顺便指出那些让我自己调试到半夜的坑。2. 陈数公式怎么写从 Berry 曲率到布里渊区积分2.1 先看三行公式再谈 Comsol 里怎么写陈数严格说叫第一陈数Chern number在二维周期系统中的标准写法是A_n(k) i u_nk|∇_k|u_nk F_n(k) ∇_k × A_n(k) 二维下只有 F_xy 分量 C_n (1/2π) ∫_{BZ} F_n(k) d²k这里u_nk是第 n 条能带的周期布洛赫函数A是 Berry 连接也叫 Berry 联络F是 Berry 曲率积分区域是整个第一布里渊区BZ。说人话陈数不是一个位置处的性质它要求你在整个 k 空间上把相位变化“卷一圈”再积分。这个积分结果理论上一定是整数所以它才是拓扑不变量。那“陈数的数怎么写”这个问题落到 Comsol 里实际上有两层。第一层是数学表达式怎么写即上面的C (1/2π)∫F dk_x dk_y。第二层是在 Comsol 的表达式栏里怎么写。Comsol 默认不直接给你A或F这些物理量。如果你的模型是“参数化扫描”扫出来的 k 点解集最顺手的方式是在全局定义里用差分近似 Berry 连接再把曲率写成连接的空间差商。以下是我常用的一个对应关系表数学符号含义Comsol 表达式写法示意A_xBerry 连接的 x 分量-imag(conj(u1)*up(u1,kx))/dkx的离散形式A_yBerry 连接的 y 分量-imag(conj(u1)*up(u1,ky))/dky的离散形式F_xyBerry 曲率二维(A_y(kxdkx)-A_y(kx-dkx))/(2*dkx)-(A_x(kydky)-A_x(ky-dky))/(2*dky)C1第一陈数对F_xy在 BZ 上求和再乘1/(2π)实际放到后处理里注意up和down是 Comsol 结果节点里常用的“访问相邻网格点解”的算子。如果你用的是特征值研究加参数化扫描它们会按扫描参数的顺序取邻点。这里有个隐藏前提k 点的扫描顺序必须是规则的否则up/down取到的“邻点”在 k 空间里根本不邻接曲率就是错的。后面第五章专门说这个坑。2.2 为什么不能直接用波函数实部积分很多初学者误以为把 Comsol 计算得到的电场或位移场实部做面积分就能得到某种拓扑荷。这是错的。陈数里的被积函数是 Berry 曲率它不是波函数本身而是波函数在 k 空间里的“相位弯曲程度”。你可以这样直观理解给某个 k 点的波函数整体乘一个额外的全局相位e^{iθ(k)}波函数本身没有物理变化但 Berry 连接会改变Berry 曲率却不会变。这个性质叫作规范不变性。如果你直接对波函数积分得到的结果依赖你在建模时引入的随机相位完全没有意义。这就是为什么我在遇到“陈数算出来是 -0.35”的求助时第一反应永远是先问对方你确认你拿的是两分量本征态在各 k 点间的交叠而不是波函数本征值吗只有当公式里出现u(k₁)|u(k₂)这种内积项相位才被正确“缠绕”追踪。接下来我们就用 Comsol 造一个能提供这种内积的最小模型。3. 用 Comsol 搭最小拓扑模型k 空间扫描与特征值提取3.1 模型选择用“系数型 PDE”写 k·p 哈密顿量想算陈数先要有一个能对每个 k 点输出本征态的系统。两个常用路线一是用波动光学模块算光子晶体本征模二是用“系数型偏微分方程Coefficient Form PDE”接口直接把 k·p 哈密顿量写进去。我一般选后者因为它的本征态是显式的复向量提取内积不需要再处理电场分量归一化。最小的二维陈数模型是带质量项的 Dirac 哈密顿量H(kx,ky) vF * (kx*σx ky*σy) m * σz其中 σ 是泡利矩阵m 是质量项。当 m 为正时该系统在 k0 处打开一个带隙价带下支的陈数是 0m 为负时陈数变为 ±1。把哈密顿量写成 2×2 矩阵对应的就是两个因变量u1、u2的耦合本征方程H * [u1; u2] E * [u1; u2]在 Comsol 的“系数型 PDE”接口里你需要把标准形式改成特征值问题。设置要点是扩散系数c置为 0阻尼/质量项da设为恒等矩阵a矩阵填-H(kx,ky)。这样特征值求解器解出来的 λ 就是能量 E特征向量就是布洛赫函数的周期部分。具体模型树路径参考下表模型开发器位置设置内容全局定义 → 参数定义kx、ky作为扫描参数vF1m0.2等组件 → 系数型偏微分方程 (c)因变量设为两个u1、u2几何选 2D 单位正方形域网格任意稀疏都行因为方程里没有空间导数项研究 1 → 特征值指定研究维度为“参数化扫描”扫描kx、ky有人会问为什么几何网格可以随便划因为系数型 PDE 里c0没有空间耦合实际上每个节点上求解的是同一个 2×2 本征问题有限元只是提供了一个形式载体。你甚至可以把网格画成两个单元不影响结果。3.2 在 Comsol 里设置参数化扫描和波函数提取我习惯在“研究 1”里新建“参数化扫描”把kx和ky作为两个独立参数。假设布里渊区是[-π, π]每个方向取n_k21个点那么扫描参数写法可以像这样kx range(-pi, pi, 2*pi/n_k) ky range(-pi, pi, 2*pi/n_k)这里range的第三个参数是步长。Comsol 的“参数化扫描”支持用range()函数生成参数列表它会自动生成所有 kx 和 ky 的组合。这一步的意义是最终结果数据集里每一组 k 点都对应一组u1、u2的本征解。后面提取陈数时需要按(kx, ky)的顺序把本征向量逐个取出。数据导出时我建议用“派生值 → 全局计算”并把输出写成表格然后通过“文件 → 导出 → 数据”保存为 CSV。在导出设置里务必勾选“求解器应计算的表达式列表”把u1、u2的实部和虚部分别导出也就是四个量re(u1)、im(u1)、re(u2)、im(u2)。如果你只在结果里看场图CSV 里是拿不到这些原始复数分量的这一步很容易被忽略。3.3 不必在 Comsol 里积分导出到外部做数值积分更稳Comsol 的后处理本身可以做全域积分但陈数计算要对“离散 k 点之间的相位交叠”求和这不是有限元意义上的单元积分而是对 k 网格的拓扑和。在 Comsol 里用表达式硬算 Berry 曲率不是不行只是当 k 点数到 60×60 以上时后处理表达式很长而且up/down算子对参数化解集的处理速度会明显变慢。所以我一般只把 Comsol 当数据源导出每个 k 点的归一化本征向量然后回 Python 或 MATLAB 里算陈数。这样也方便后面换不同的 k 网格密度做收敛性验证。4. 在 Python 里把陈数的数算出来三角形相位累积法4.1 为什么要用三角形而不是四边形最常见的陈数数值算法是所谓的“三角形累积法”也叫 Fukui–Hatsugai–Suzuki 方法。核心思想是把二维布里渊区离散成若干个三角形对每个三角形的三个顶点(k1, k2, k3)计算相邻本征态的内积U12 u(k1)|u(k2) / |u(k1)|u(k2)| U23 u(k2)|u(k3) / |u(k2)|u(k3)| U31 u(k3)|u(k1) / |u(k3)|u(k1)|然后把U12 * U23 * U31的虚部对数加起来。每个三角形贡献的相位角正比于 Berry 曲率在这个小区域上的通量所有三角形累加后除以2π就是陈数。为什么不直接用矩形网格差分因为 Berry 连接在 k 空间的单值化U(1) 相位展开在三角形上做最自然。矩形的四个顶点内积会有一个多余的规范自由度处理不好容易漏掉拓扑贡献。三角形算法只需要相邻顶点内积边的方向固定相位累积天然闭合。以下是完整的可运行核心代码我从 CSV 读取数据算陈数import numpy as np import pandas as pd def normalize(v): # v 是复数向量长度等于波函数分量数 norm np.sqrt(np.vdot(v, v).real) return v / norm def read_comsol_csv(path): # Comsol 导出列名kx, ky, eig, re_u1, im_u1, re_u2, im_u2 df pd.read_csv(path) nk int(np.sqrt(len(df))) kx df[kx].to_numpy().reshape(nk, nk) ky df[ky].to_numpy().reshape(nk, nk) psi np.empty((nk, nk, 2), dtypecomplex) psi[..., 0] df[re_u1].to_numpy().reshape(nk, nk) 1j * df[im_u1].to_numpy().reshape(nk, nk) psi[..., 1] df[re_u2].to_numpy().reshape(nk, nk) 1j * df[im_u2].to_numpy().reshape(nk, nk) return kx, ky, psi def chern_triangle(kx, ky, psi): nk kx.shape[0] total 0.0 for ix in range(nk - 1): for iy in range(nk - 1): # 当前格点及其右、上、右上三个邻近格点 p1 (ix, iy) p2 (ix 1, iy) p3 (ix, iy 1) p4 (ix 1, iy 1) # 每个小矩形分成两个三角形 for tri in [(p1, p2, p3), (p2, p4, p3)]: v1 normalize(psi[tri[0]]) v2 normalize(psi[tri[1]]) v3 normalize(psi[tri[2]]) u12 np.vdot(v1, v2) u23 np.vdot(v2, v3) u31 np.vdot(v3, v1) phase np.angle(u12 * u23 * u31) total phase chern total / (2.0 * np.pi) return chern kx, ky, psi read_comsol_csv(comsol_export.csv) c1 chern_triangle(kx, ky, psi) print(fChern number {c1:.4f})这段代码的关键点有三个。第一np.vdot(v1, v2)做的是复数内积它自动对v1取共轭这正是量子力学里u(k1)|u(k2)的定义。第二每个内积后需要归一化但如果你在 Comsol 里已经归一化过本征向量这里可以跳过写成v1 psi[tri[0]]即可保留归一化步骤更保险。第三np.angle的结果范围是[-π, π]单个三角形贡献的相位必须远小于 π这要求 k 网格足够密如果网格太粗相位折叠会让最后结果变成随机小数。4.2 参数怎么给网格密度、质量和 BZ 范围在实际运行里我推荐先用小网格跑通再加大。比如先用n_k11验证程序不出错然后n_k31看陈数是否稳定收敛。下面是几个好用的参数起点参数建议值说明vF1速度项影响能带的斜率和带隙位置m0.2质量项必须非零才有陈数nk31每个方向采样数陈数对网格密度不敏感但 11 以下可能不够BZ 范围[-π, π]如果哈密顿量里 a1直接按标准 BZ你会在运行中发现一个有趣现象当m0.2时结果稳定在1.000或-1.000当m从正数连续调到负数时陈数直接跳变而不是连续过渡。这就是拓扑相变。如果算出来是非整数比如0.42或0.85不要急着调网格先按第四章的办法检查 k 点排序和三角剖分方向。5. 陈数结果不自洽时先查这三个位置5.1 带隙闭合点与 k 点排序第一个要查的是带隙是否真的在计算区域内完全打开。陈数只在绝热带隙不为零的前提下才是良定义的。你在扫参时如果m设为 0Dirac 点在 k0 处闭合陈数没有任何意义。所以我一般先画一条kx方向的能带图确认色散在布里渊区边界和最 Γ 点附近都有 gap再动手算陈数。第二个高发问题是参数化扫描的排序。Comsol 默认按第一个扫描参数优先变化生成解序列但导出的 CSV 行顺序不一定和你的kx、ky网格形状一致。如果顺序混乱我的reshape(nk, nk)会把 (kx, ky) 映射到错误的格点上三角形顶点方向反转得到的陈数符号会反过来。建议导出时顺手把kx、ky也导出来然后在 Python 里先做一次排序检查# 检查网格是否规则kx 沿行递增ky 沿列递增 dx np.diff(kx, axis1) dy np.diff(ky, axis0) print(dx.min(), dx.max(), dy.min(), dy.max())如果dx出现负值说明 kx 的方向在某一列反了如果dy全为 0 而dx有正有负同样是重排问题。最稳妥的办法是在 Comsol 里把扫描顺序固定为“外层 ky内层 kx”导出后按np.lexsort((kx_flat, ky_flat))重新索引数据再做 reshape。5.2 Wilson loop一个不依赖网格密度的验证技巧除了直接用三角形算法验证整数性质我还会在同一个数据集上算 Wilson loop。它的思路是固定一个ky沿着 kx 方向一圈把相邻点之间的贝里相位逐个累加当 ky 从 0 扫到 2π 时这个累计相位的变化绕了几圈就对应陈数。这个方法的直观意义是一条“相位缠绕组”经得起局部相位跳变的干扰。代码如下def wilson_loop(psi, ky_index): nk psi.shape[0] phase 0.0 for ix in range(nk): v1 normalize(psi[ky_index, ix]) v2 normalize(psi[ky_index, (ix 1) % nk]) phase np.angle(np.vdot(v1, v2)) return phase / (2.0 * np.pi) # 对每个 ky 行都算一圈 windings [wilson_loop(psi, iy) for iy in range(psi.shape[1])] print(np.median(windings), np.std(windings))当np.std(windings)小于0.05时说明陈数结果可信如果标准差很大通常意味着系统在某个 ky 处带隙接近闭合要么把m调大一点要么把 BZ 的离散范围缩到带隙稳定的区域。Wilson loop 的好处在于它不依赖每行之间的一致性只依赖沿 kx 闭合路径的相位累积因此它对 k 网格排序的要求比三角形算法更低非常适合作为交叉验证。最后一个小技巧在 Comsol 里保持波函数相位固定时尽量用特征值研究里的归一化选项并导出未经过场旋转的复数分量这样导出的数据在不同 k 点间才存在统一的参考相位你的陈数自然落在正确的整数上。本文还有配套的精品资源点击获取 COMSOL中计算陈数：从Berry曲率到布里渊区积分
简介这套资料围绕光子晶体陈数计算与COMSOL仿真面向具备一定光学或电磁场基础的研究生、工程师解决能带结构计算及陈数后处理中的具体问题。压缩包共含5个文件包括两篇PDF文档、一个mph模型文件、一个MATLAB脚本及一份txt说明整体约20MB可覆盖从理论原理到软件操作的完整链条。截至目前已有1042人学习浏览适合正在做光子晶体拓扑性质研究的读者参考。内容兼顾陈数概念讲解、COMSOL建模思路与MATLAB数据处理还涉及布里渊区分析和能带图的绘制其中mph模型可直观查看仿真设置m脚本方便修改参数并复用配合说明文档能帮助读者快速完成从电磁场模拟到陈数计算的实践显著缩短自行摸索周期。1. 在 Comsol 里算陈数到底要算哪个数很多做拓扑光子学或声子晶体的人把能带图跑出来后卡在同一句话上陈数的数怎么写。有人到全局定义里建了一个chern填0或者1发现什么都不变有人把某一支本征态的实部导出来做面积分得出一个连自己都不信的 0.03。问题不在 Comsol 用得不熟而是陈数根本就不是由单个波函数算出来的。它是布里渊区上 Berry 曲率的积分Comsol 在这件事里只负责提供每个 k 点的本征态真正的“数”要在 k 空间网格上做离散求和。本文就从公式写法、模型搭建、数值积分到结果验证把这条链路走通顺便指出那些让我自己调试到半夜的坑。2. 陈数公式怎么写从 Berry 曲率到布里渊区积分2.1 先看三行公式再谈 Comsol 里怎么写陈数严格说叫第一陈数Chern number在二维周期系统中的标准写法是A_n(k) i u_nk|∇_k|u_nk F_n(k) ∇_k × A_n(k) 二维下只有 F_xy 分量 C_n (1/2π) ∫_{BZ} F_n(k) d²k这里u_nk是第 n 条能带的周期布洛赫函数A是 Berry 连接也叫 Berry 联络F是 Berry 曲率积分区域是整个第一布里渊区BZ。说人话陈数不是一个位置处的性质它要求你在整个 k 空间上把相位变化“卷一圈”再积分。这个积分结果理论上一定是整数所以它才是拓扑不变量。那“陈数的数怎么写”这个问题落到 Comsol 里实际上有两层。第一层是数学表达式怎么写即上面的C (1/2π)∫F dk_x dk_y。第二层是在 Comsol 的表达式栏里怎么写。Comsol 默认不直接给你A或F这些物理量。如果你的模型是“参数化扫描”扫出来的 k 点解集最顺手的方式是在全局定义里用差分近似 Berry 连接再把曲率写成连接的空间差商。以下是我常用的一个对应关系表数学符号含义Comsol 表达式写法示意A_xBerry 连接的 x 分量-imag(conj(u1)*up(u1,kx))/dkx的离散形式A_yBerry 连接的 y 分量-imag(conj(u1)*up(u1,ky))/dky的离散形式F_xyBerry 曲率二维(A_y(kxdkx)-A_y(kx-dkx))/(2*dkx)-(A_x(kydky)-A_x(ky-dky))/(2*dky)C1第一陈数对F_xy在 BZ 上求和再乘1/(2π)实际放到后处理里注意up和down是 Comsol 结果节点里常用的“访问相邻网格点解”的算子。如果你用的是特征值研究加参数化扫描它们会按扫描参数的顺序取邻点。这里有个隐藏前提k 点的扫描顺序必须是规则的否则up/down取到的“邻点”在 k 空间里根本不邻接曲率就是错的。后面第五章专门说这个坑。2.2 为什么不能直接用波函数实部积分很多初学者误以为把 Comsol 计算得到的电场或位移场实部做面积分就能得到某种拓扑荷。这是错的。陈数里的被积函数是 Berry 曲率它不是波函数本身而是波函数在 k 空间里的“相位弯曲程度”。你可以这样直观理解给某个 k 点的波函数整体乘一个额外的全局相位e^{iθ(k)}波函数本身没有物理变化但 Berry 连接会改变Berry 曲率却不会变。这个性质叫作规范不变性。如果你直接对波函数积分得到的结果依赖你在建模时引入的随机相位完全没有意义。这就是为什么我在遇到“陈数算出来是 -0.35”的求助时第一反应永远是先问对方你确认你拿的是两分量本征态在各 k 点间的交叠而不是波函数本征值吗只有当公式里出现u(k₁)|u(k₂)这种内积项相位才被正确“缠绕”追踪。接下来我们就用 Comsol 造一个能提供这种内积的最小模型。3. 用 Comsol 搭最小拓扑模型k 空间扫描与特征值提取3.1 模型选择用“系数型 PDE”写 k·p 哈密顿量想算陈数先要有一个能对每个 k 点输出本征态的系统。两个常用路线一是用波动光学模块算光子晶体本征模二是用“系数型偏微分方程Coefficient Form PDE”接口直接把 k·p 哈密顿量写进去。我一般选后者因为它的本征态是显式的复向量提取内积不需要再处理电场分量归一化。最小的二维陈数模型是带质量项的 Dirac 哈密顿量H(kx,ky) vF * (kx*σx ky*σy) m * σz其中 σ 是泡利矩阵m 是质量项。当 m 为正时该系统在 k0 处打开一个带隙价带下支的陈数是 0m 为负时陈数变为 ±1。把哈密顿量写成 2×2 矩阵对应的就是两个因变量u1、u2的耦合本征方程H * [u1; u2] E * [u1; u2]在 Comsol 的“系数型 PDE”接口里你需要把标准形式改成特征值问题。设置要点是扩散系数c置为 0阻尼/质量项da设为恒等矩阵a矩阵填-H(kx,ky)。这样特征值求解器解出来的 λ 就是能量 E特征向量就是布洛赫函数的周期部分。具体模型树路径参考下表模型开发器位置设置内容全局定义 → 参数定义kx、ky作为扫描参数vF1m0.2等组件 → 系数型偏微分方程 (c)因变量设为两个u1、u2几何选 2D 单位正方形域网格任意稀疏都行因为方程里没有空间导数项研究 1 → 特征值指定研究维度为“参数化扫描”扫描kx、ky有人会问为什么几何网格可以随便划因为系数型 PDE 里c0没有空间耦合实际上每个节点上求解的是同一个 2×2 本征问题有限元只是提供了一个形式载体。你甚至可以把网格画成两个单元不影响结果。3.2 在 Comsol 里设置参数化扫描和波函数提取我习惯在“研究 1”里新建“参数化扫描”把kx和ky作为两个独立参数。假设布里渊区是[-π, π]每个方向取n_k21个点那么扫描参数写法可以像这样kx range(-pi, pi, 2*pi/n_k) ky range(-pi, pi, 2*pi/n_k)这里range的第三个参数是步长。Comsol 的“参数化扫描”支持用range()函数生成参数列表它会自动生成所有 kx 和 ky 的组合。这一步的意义是最终结果数据集里每一组 k 点都对应一组u1、u2的本征解。后面提取陈数时需要按(kx, ky)的顺序把本征向量逐个取出。数据导出时我建议用“派生值 → 全局计算”并把输出写成表格然后通过“文件 → 导出 → 数据”保存为 CSV。在导出设置里务必勾选“求解器应计算的表达式列表”把u1、u2的实部和虚部分别导出也就是四个量re(u1)、im(u1)、re(u2)、im(u2)。如果你只在结果里看场图CSV 里是拿不到这些原始复数分量的这一步很容易被忽略。3.3 不必在 Comsol 里积分导出到外部做数值积分更稳Comsol 的后处理本身可以做全域积分但陈数计算要对“离散 k 点之间的相位交叠”求和这不是有限元意义上的单元积分而是对 k 网格的拓扑和。在 Comsol 里用表达式硬算 Berry 曲率不是不行只是当 k 点数到 60×60 以上时后处理表达式很长而且up/down算子对参数化解集的处理速度会明显变慢。所以我一般只把 Comsol 当数据源导出每个 k 点的归一化本征向量然后回 Python 或 MATLAB 里算陈数。这样也方便后面换不同的 k 网格密度做收敛性验证。4. 在 Python 里把陈数的数算出来三角形相位累积法4.1 为什么要用三角形而不是四边形最常见的陈数数值算法是所谓的“三角形累积法”也叫 Fukui–Hatsugai–Suzuki 方法。核心思想是把二维布里渊区离散成若干个三角形对每个三角形的三个顶点(k1, k2, k3)计算相邻本征态的内积U12 u(k1)|u(k2) / |u(k1)|u(k2)| U23 u(k2)|u(k3) / |u(k2)|u(k3)| U31 u(k3)|u(k1) / |u(k3)|u(k1)|然后把U12 * U23 * U31的虚部对数加起来。每个三角形贡献的相位角正比于 Berry 曲率在这个小区域上的通量所有三角形累加后除以2π就是陈数。为什么不直接用矩形网格差分因为 Berry 连接在 k 空间的单值化U(1) 相位展开在三角形上做最自然。矩形的四个顶点内积会有一个多余的规范自由度处理不好容易漏掉拓扑贡献。三角形算法只需要相邻顶点内积边的方向固定相位累积天然闭合。以下是完整的可运行核心代码我从 CSV 读取数据算陈数import numpy as np import pandas as pd def normalize(v): # v 是复数向量长度等于波函数分量数 norm np.sqrt(np.vdot(v, v).real) return v / norm def read_comsol_csv(path): # Comsol 导出列名kx, ky, eig, re_u1, im_u1, re_u2, im_u2 df pd.read_csv(path) nk int(np.sqrt(len(df))) kx df[kx].to_numpy().reshape(nk, nk) ky df[ky].to_numpy().reshape(nk, nk) psi np.empty((nk, nk, 2), dtypecomplex) psi[..., 0] df[re_u1].to_numpy().reshape(nk, nk) 1j * df[im_u1].to_numpy().reshape(nk, nk) psi[..., 1] df[re_u2].to_numpy().reshape(nk, nk) 1j * df[im_u2].to_numpy().reshape(nk, nk) return kx, ky, psi def chern_triangle(kx, ky, psi): nk kx.shape[0] total 0.0 for ix in range(nk - 1): for iy in range(nk - 1): # 当前格点及其右、上、右上三个邻近格点 p1 (ix, iy) p2 (ix 1, iy) p3 (ix, iy 1) p4 (ix 1, iy 1) # 每个小矩形分成两个三角形 for tri in [(p1, p2, p3), (p2, p4, p3)]: v1 normalize(psi[tri[0]]) v2 normalize(psi[tri[1]]) v3 normalize(psi[tri[2]]) u12 np.vdot(v1, v2) u23 np.vdot(v2, v3) u31 np.vdot(v3, v1) phase np.angle(u12 * u23 * u31) total phase chern total / (2.0 * np.pi) return chern kx, ky, psi read_comsol_csv(comsol_export.csv) c1 chern_triangle(kx, ky, psi) print(fChern number {c1:.4f})这段代码的关键点有三个。第一np.vdot(v1, v2)做的是复数内积它自动对v1取共轭这正是量子力学里u(k1)|u(k2)的定义。第二每个内积后需要归一化但如果你在 Comsol 里已经归一化过本征向量这里可以跳过写成v1 psi[tri[0]]即可保留归一化步骤更保险。第三np.angle的结果范围是[-π, π]单个三角形贡献的相位必须远小于 π这要求 k 网格足够密如果网格太粗相位折叠会让最后结果变成随机小数。4.2 参数怎么给网格密度、质量和 BZ 范围在实际运行里我推荐先用小网格跑通再加大。比如先用n_k11验证程序不出错然后n_k31看陈数是否稳定收敛。下面是几个好用的参数起点参数建议值说明vF1速度项影响能带的斜率和带隙位置m0.2质量项必须非零才有陈数nk31每个方向采样数陈数对网格密度不敏感但 11 以下可能不够BZ 范围[-π, π]如果哈密顿量里 a1直接按标准 BZ你会在运行中发现一个有趣现象当m0.2时结果稳定在1.000或-1.000当m从正数连续调到负数时陈数直接跳变而不是连续过渡。这就是拓扑相变。如果算出来是非整数比如0.42或0.85不要急着调网格先按第四章的办法检查 k 点排序和三角剖分方向。5. 陈数结果不自洽时先查这三个位置5.1 带隙闭合点与 k 点排序第一个要查的是带隙是否真的在计算区域内完全打开。陈数只在绝热带隙不为零的前提下才是良定义的。你在扫参时如果m设为 0Dirac 点在 k0 处闭合陈数没有任何意义。所以我一般先画一条kx方向的能带图确认色散在布里渊区边界和最 Γ 点附近都有 gap再动手算陈数。第二个高发问题是参数化扫描的排序。Comsol 默认按第一个扫描参数优先变化生成解序列但导出的 CSV 行顺序不一定和你的kx、ky网格形状一致。如果顺序混乱我的reshape(nk, nk)会把 (kx, ky) 映射到错误的格点上三角形顶点方向反转得到的陈数符号会反过来。建议导出时顺手把kx、ky也导出来然后在 Python 里先做一次排序检查# 检查网格是否规则kx 沿行递增ky 沿列递增 dx np.diff(kx, axis1) dy np.diff(ky, axis0) print(dx.min(), dx.max(), dy.min(), dy.max())如果dx出现负值说明 kx 的方向在某一列反了如果dy全为 0 而dx有正有负同样是重排问题。最稳妥的办法是在 Comsol 里把扫描顺序固定为“外层 ky内层 kx”导出后按np.lexsort((kx_flat, ky_flat))重新索引数据再做 reshape。5.2 Wilson loop一个不依赖网格密度的验证技巧除了直接用三角形算法验证整数性质我还会在同一个数据集上算 Wilson loop。它的思路是固定一个ky沿着 kx 方向一圈把相邻点之间的贝里相位逐个累加当 ky 从 0 扫到 2π 时这个累计相位的变化绕了几圈就对应陈数。这个方法的直观意义是一条“相位缠绕组”经得起局部相位跳变的干扰。代码如下def wilson_loop(psi, ky_index): nk psi.shape[0] phase 0.0 for ix in range(nk): v1 normalize(psi[ky_index, ix]) v2 normalize(psi[ky_index, (ix 1) % nk]) phase np.angle(np.vdot(v1, v2)) return phase / (2.0 * np.pi) # 对每个 ky 行都算一圈 windings [wilson_loop(psi, iy) for iy in range(psi.shape[1])] print(np.median(windings), np.std(windings))当np.std(windings)小于0.05时说明陈数结果可信如果标准差很大通常意味着系统在某个 ky 处带隙接近闭合要么把m调大一点要么把 BZ 的离散范围缩到带隙稳定的区域。Wilson loop 的好处在于它不依赖每行之间的一致性只依赖沿 kx 闭合路径的相位累积因此它对 k 网格排序的要求比三角形算法更低非常适合作为交叉验证。最后一个小技巧在 Comsol 里保持波函数相位固定时尽量用特征值研究里的归一化选项并导出未经过场旋转的复数分量这样导出的数据在不同 k 点间才存在统一的参考相位你的陈数自然落在正确的整数上。本文还有配套的精品资源点击获取