全布里渊区第一性原理电子–声子耦合
完整三维硅 k/q 表示:以直接第一性原理参考检验 406 项回归查询、19 项独立电子 k 查询与 16 项独立联合 k/q 查询。
Felix Fei
全布里渊区第一性原理电子–声子耦合
完整三维硅 k/q 表示:以直接第一性原理参考检验 406 项回归查询、19 项独立电子 k 查询与 16 项独立联合 k/q 查询。
2026-10-20 · Mol2Mat · Felix Fei · revision 2
硅晶体研究里程碑
这一计算现已接入 Mol2Mat,可直接在浏览器工作台中使用。
电子–声子耦合计算连接三个对象:位于 的电子态、波矢为 的晶格振动,以及位于 的另一个电子态。准确计算每个对象只是其中一部分。当这些对象组合在一起,并从构建模型所用的网格点移向其他位置时,它们的相位、归一化、对称性和数值表示都必须保持一致。
我们已在 Mol2Mat 的金刚石结构硅四价带研究模型中完成这条计算链路。实现从自洽密度泛函理论出发,在完整的 8 × 8 × 8 电子网格与声子网格上构建声子动力学矩阵和复数电子–声子顶角,再通过 Fourier 插值计算连续 坐标处的响应。
最终模型在原有 5 % EPC 判据下通过了全部 406 项既有回归查询、19 项独立电子 查询和 16 项独立联合 查询。三组验证的最大误差分别为 4.2691 %、3.4066 % 和 4.5619 %。独立验证的预测值均在直接参考结果生成前冻结,整个验证过程中模型保持固定。
这一结果为指定的硅模型建立了可复现的三维动量空间第一性原理 EPC 计算。实现保留完整超软赝势顶角,确定响应所需的插值分辨率,并对弱耦合通道逐一验证。
从电子结构到复数散射振幅
电子结构计算采用 PBE 和 GBRV 硅超软赝势[1, 2],波函数与密度截断能分别为 20 和 100 Ha。金刚石结构原胞包含两个原子,对应六个笛卡尔位移坐标。电子子空间包含四条占据价带,由四个 Wannier 函数表示;计算采用非自旋极化设置,硅原子质量取 28.085 u。
对于原子 沿 方向的笛卡尔扰动,我们保留一个复数 4 × 4 导数矩阵,记为 ,单位为 Ha bohr⁻¹。对质量加权动力学矩阵进行对角化后,声子极化矢量与零点位移共同将该导数转化为逐模耦合:
其中,、 标记矩阵采用的电子基, 标记声子模式。存储的顶角采用 Wannier 基;对初态与末态电子端点分别施加幺正变换,即可得到相应的 Bloch 能带矩阵。将未加权的笛卡尔顶角与声子收缩分开保存,有助于明确插值过程和误差归因。
最终张量的形状为 (512, 512, 6, 4, 4),依次对应初态 、声子 、位移坐标、末态 Wannier 指标和初态 Wannier 指标。它包含 262,144 组 组合、25,165,824 个复数导数元素,以及 512 个复数 6 × 6 动力学矩阵。这两个网格决定插值分辨率;自洽电子积分则采用另行记录的求积方案。
这种拆分遵循 EPC 计算的物理结构。EPW [3] 采用局域表示,将电子、振动和耦合量插值到细密动量网格上;Perturbo [4] 则将这些第一性原理矩阵元用于输运与载流子动力学。在 Mol2Mat 中,此次里程碑是硅耦合基础模型的构建及独立数值验证,具体通过有限位移 JDFTx 响应和显式算符投影实现。
超软赝势顶角的一致构建
采用超软赝势时,电子本征问题包含重叠算符:
原子位移同时改变 和 。响应包含局域自洽势、非局域投影子、增广项和电子基的输运修正。两个电子端点也具有不同的 Bloch 相位。在对正负位移作差之前,每一项都必须采用相同的参考约定。
实现保存了投影矩阵 、、、、 和 。按照所存储的算符约定,中心化并经输运修正的矩阵为
平均势确定能量零点约定, 则提供双侧输运修正。组装前,代码检查恒等式 ,并检查参考 Hamiltonian 与重叠矩阵恒等式。原生投影使用已保存的自洽势,因此计算新的电子端点时可以复用已收敛响应。
中心差分通常沿选定单位方向采用 ±0.02 bohr 位移。原胞电子态嵌入 N 原胞超胞后,每个端点都带有 因子。因此,对投影矩阵块作差时,组装导数需恢复 因子。组合不同形状超胞的数据时,这一归一化尤其重要。
空间群输运同时作用于响应的各个部分:笛卡尔方向、原子置换、晶格平移偏移和两个电子端点的旋转。倒空间边界的重新标记同样属于算符定义的一部分。坐标越过布里渊区边界时,必须同步重标记平面波 G 指标;仅保留折返后的坐标数值还不够。
这些处理共同建立了一致的复数振幅。单个标量耦合强度会丢失检验这种一致性所需的相位和通道信息。
响应所需的插值分辨率
最终模型的形成过程揭示了一个重要区别:数值重构准确,并不意味着物理表示具有足够的表达能力。
早期 4³ 电子 Wannier–Fourier 表示能够以数值精度重现全部 64 个构建点,但在其他已保存 点上仍有较大误差。于是,我们检验了同一固定 Fourier 空间中的任意一组系数,能否足够准确地重现已有响应。
对每组电子矩阵,令 表示在 456 个已保存观测点上计算固定的 64 个基函数, 表示对应的直接响应, 表示到 列空间的投影。残差可以作正交分解:
第一项位于所选表示空间之外,调整空间内的系数无法消除它。通过薄 QR 分解进行投影,得到完整笛卡尔顶角的最小整体相对误差为 5.1642 %,最低声学模组为 8.6529 %,均超过 5 %。如果每个点的误差都不超过 5 %,整体 Frobenius 范数比值也应不超过 5 %。因此,这个固定表示无法在全部 456 个观测点上满足该条件。
我们还通过独立的有限位移对照,检查直接参考值是否足以解释历史上的大误差。在被审查点,将位移从 0.02 bohr 减半至 0.01 bohr,只使问题模式的参考矩阵变化 0.02742 %;冻结预测的误差则从 21.3081 % 变为 21.2993 %。这一配对对照将后续研究指向表示空间和采样分辨率。
电子网格提高到 8³ 后,五个固定 回归查询的最大误差从 18.1314 % 降至 2.0051 %。另外 392 个固定 回归点也全部通过,最大误差为 4.2691 %。
声子波矢方向仍需要独立提高分辨率。在 k8 × q4 下, = (2,3,5)/12 处八个非轴向查询的最大误差仍为 14.7329 %;完成 q8 后,同一组查询降至 3.6755 %。轴向 = (1/3,0,0) 查询也从 24.1940 % 改善至 3.2333 %。这里的坐标均为原胞倒格基下的分数坐标。
这一演进给出了计算困难的具体来源:电子 分辨率和声子 分辨率约束同一振幅的不同部分。提高一个坐标空间的分辨率后,另一个空间仍可能成为主导误差来源。
以 33 项新增自洽响应补齐 q8
q8 扩展复用了已有 q4 响应和轴向数据。十二组八原胞响应提供新增非轴向信息,一组保留的 24 原胞轴向响应再提供 16 个 节点。结合对称性,这些来源补齐了全部 512 个 节点:64 个保留的 q4 节点、16 个轴向节点和 432 个新增节点。
最初的正负位移计划包含 36 个 SCF 状态。晶体对称性将三个负位移状态映射到对应的正位移状态,把实际新增计算减少至 33 项。这些操作经过实际原子映射、带权电子积分网格和 FFT 网格的检查,再由已保存的正位移响应通过核验后的变换生成对应状态。
33 项自洽计算全部完成。每项输入均包含 16 个原子、32 条电子能带和 432 个显式带权电子积分点。收敛规则要求连续两个周期达到原有 10⁻⁹ SCF 残差阈值。后续用于独立联合验证的另外八项 SCF 属于单独的参考计算批次。
组装完成的笛卡尔顶角沿六个动量网格轴进行 Fourier 变换。连续坐标查询通过 和 的相位因子计算相应实空间系数:
实现中的相位因子对等价的最近周期镜像取平均。电子镜像的选择取决于 Wannier 中心之差;响应镜像的选择取决于原子相对于末态 Wannier 中心的位置;晶格度量决定哪些镜像最近。这些固定规则在保留粗网格相位恒等式的同时,确定节点之间的插值。动力学矩阵则采用相应的原子中心 空间构建。
声子组装中还出现了一个细微的对称性问题:简并模式极小的数值劈裂,会让原有模组比较器无法找到维数相同的匹配。解决方法是直接对每个新增动力学矩阵施加晶体的小群对称性。矩阵 Hermitian 化后,通过 Reynolds 平均施加所有保持 不变的操作:
原子位移表示 包含晶格偏移相位。在检查的 84 个 点中,最大相对矩阵修正为 0.02246 %,最大协变性残差为 3.22 × 10⁻¹⁶。所有公共点的模组维数均匹配;原始 q4 数值精确保留,原有分组容差保持固定。
对称性减少了所需响应计算的数量,也恢复了逐模比较所需的简并结构。
冻结预测与独立直接验证
最终评估分别检验三类问题。
既有回归。 406 项查询由 392 个已保存的固定 电子点、五个额外固定 对照、八个非轴向联合查询和一个轴向离网格查询组成。其直接结果在开发期间已经获得。最终 q8 模型通过全部 406 项。
独立电子 。 十九项新增电子查询经过与已暴露对称轨道的检查,覆盖晶体操作、时间反演和逆散射。预测在新增直接投影组装前冻结。由于 位于已有声子节点,这组验证检验新电子端点处的插值;十九项全部通过。
独立联合 。 十六项查询覆盖两组非轴向声子波矢 = (1,2,3)/12 和 = (1,4,7)/12,每组包含八项电子查询。八项新增正负位移 SCF 响应提供独立参考。独立参考采用 208 个固定势投影块构建。十六项查询均通过与参考生成前冻结预测的比较。
模型张量、Fourier 规则和 5 % 阈值始终固定。公共来源对照另外使用 0.5 % 判据,检查新参考组装与既有来源的一致性。
| 验证集合 | 通过查询数 | 仅电子部分最大误差 | 仅声子部分最大误差 | 联合最大误差 |
|---|---|---|---|---|
| 既有回归 | 406 / 406 | 4.2691 % | 0.6853 % | 4.2691 % |
| 独立电子 | 19 / 19 | 3.4066 % | 0.00001261 % | 3.4066 % |
| 独立联合 | 16 / 16 | 4.5619 % | 1.9479 % | 4.5619 % |
仅电子比较将预测顶角与参考动力学矩阵配对;仅声子比较将参考顶角与预测动力学矩阵配对;联合比较同时使用两项预测。这些误差来自同一耦合计算中的受控替换,不能直接相加。
这一结果结合了广泛的回归覆盖和 35 项前瞻性独立查询。Full-BZ 指连续三维 表示,数值精度由这些明确定义的测试集度量。
弱耦合通道的逐模验证
单个 EPC 范数不足以完整检验顶角。两个复数矩阵可以具有几乎相同的范数,却在电子通道空间中具有不同方向。较小的声学响应也可能来自更大原子贡献之间的抵消,从而放大最终模式中的相对误差。
因此,比较器保留完整原始顶角,并逐一检查声子模组,比较同规范下的矩阵差、耦合范数和奇异值谱。在简并模组内,将模式矩阵堆叠后得到的通道谱对模式基的幺正混合不变。预测与参考声子子空间通过最大化重叠的等维双射进行匹配;此外,也使用各自计算得到的模式评估耦合范数与通道谱。
例如,同规范矩阵误差定义为
单点评分取完整顶角误差,以及适用模组的矩阵、范数和通道谱误差中的最大值。无法获得有效模式匹配时,该点不能通过。零频或不稳定模式保留带符号本征值和诊断状态,不人为赋予正频率。
最大独立联合误差出现在 = (1,2,3)/12、 = (1,4,7)/12。其最弱模式的直接耦合范数为 74.1102 meV,预测为 73.6286 meV,范数误差仅 0.6498 %;但复数矩阵误差为 4.5619 %,通道谱误差为 2.1769 %。该模式只占原始笛卡尔顶角范数平方的 0.5431 %。其原子抵消比,即净模式振幅与两个原子振幅之和的比值,为 0.1333。
总强度图会给这一苛刻模式很低的权重,逐模比较则把它纳入报告精度的核心。
该查询的声子位移图提供了互补的实空间视角。最低模频率在直接计算中为 4.6445 THz,在冻结预测中为 4.6260 THz,两者本征矢重叠的模平方为 0.999994。振动图样的高度一致性与分项评分共同表明,较大的剩余误差位于电子顶角。
配套的交互式三维视图支持旋转、全部六个模式的选择、相位调整,并可导出 PNG 和 GLB。
这些检验针对后续散射计算所需的矩阵元。物理可观测量的收敛还取决于电子能量窗口、相空间积分、占据和能量守恒处理。保留完整复数顶角,为后续计算提供了比标量拟合摘要更扎实的基础。
八核心并行性能
本次硅 EPC 计算采用八核心预算。自洽响应和固定势投影可按独立任务并行处理;已有自洽势可以用于新的电子端点,减少重复计算。
我们在同一台机器上比较了两种核心分配。测试使用相同的八个已保存 SCF 状态,每个状态继续计算两个周期,耗时包含启动、检查点保存和结果输出:
| 配置 | 耗时 | 聚合驻留内存峰值 |
|---|---|---|
| 四进程 × 每进程两线程 | 704.67 s | 5.03 GiB |
| 八进程 × 每进程一线程 | 579.69 s | 9.59 GiB |
八个单线程任务的总耗时减少 17.74 %,内存峰值从 5.03 GiB 增至 9.59 GiB。两种配置的最终能量极差为 1.28 × 10⁻¹³ Ha。这组结果给出了明确的选择依据:内存充足时,可以增加同时运行的独立任务,以缩短等待时间。
线程分配也随剩余任务数调整:最多同时运行八个单线程任务,剩余四项或三项时每项使用两线程,剩余两项时每项四线程,最后一项可使用八线程。这样,批次接近完成时仍能利用可用核心。
连续动量查询与后续计算
这套模型支持连续三维 坐标处的顶角、声子模式和逐模耦合查询。研究者可以检查具体散射通道,将预测与直接参考逐项比较,并使用保存的矩阵复现误差。已有 441 项插值重放记录以 9.86 × 10⁻¹⁶ 的最大相对差异重现了冻结预测。
同一冻结模型现已接入浏览器的 FP64 WASM 计算工作流,用户可以执行连续 查询、查看声子模式与耦合、比较不同波矢,并导出复数矩阵。后续将通过外部矩阵对照、系统网格收敛和更大电子子空间检验其适用范围,在此基础上计算由这些顶角决定的散射率和输运性质。
在浏览器中使用
在 DFT 单点计算中选择电子–声子耦合,输入以原胞倒格矢为基底的电子波矢 和声子波矢 ,即可比较六个振动模式及不同波矢下的耦合强度。选中的模式与晶格振动动画联动;完整复数耦合矩阵可在详情中查看和导出。
浏览器使用上述 8³ 电子网格与声子网格构建的硅模型。首次使用需下载约 384.5 MiB 的模型数据,并校验文件大小和 SHA-256。下载完成后,查询由本机的 WebAssembly 内核执行。该内核由 C 编译而来,以双精度浮点数(FP64)计算复数 Fourier 插值,并复用相同 对应的中间结果。计算在独立 Worker 中进行,使页面在查询过程中保持响应。
插值得到的复数顶角和动力学矩阵用于计算声子频率、原子位移和逐模耦合。对于简并模式,界面同时给出模组耦合矩阵的 Frobenius 范数,便于比较不同波矢下的耦合强度。导出文件包含完整复数矩阵及对应的模型与内核版本。
我们在 Chrome 149 中运行这一实现,并将输出与 NumPy 参考程序逐项比较。全部 441 项查询中,原始顶角和动力学矩阵的最大相对差异分别为 2.11 × 10⁻¹⁵ 和 1.43 × 10⁻¹⁵;声子频率和模组范数也通过了相同参考的核验。再与直接第一性原理结果比较,406 项回归查询、19 项独立电子查询和 16 项独立联合查询全部通过原有的 5 % 判据。详细结果见浏览器计算验证摘要。
配套数据
配套表格包含全部 441 项查询评分、96 个独立联合模组比较、配对网格比较和固定表示空间误差下界。浏览器 WASM 验证摘要提供 441 项查询的移植核验结果。图表提供 SVG、PDF 与 PNG 格式。
Data and sources
- Browser WASM verification summary
- All 441 query scores
- Independent joint mode-group comparisons
- Matched grid comparisons
- Fixed-representation error bounds
- R. Sundararaman et al. JDFTx: Software for joint density-functional theory. SoftwareX 6, 278–284 (2017). doi:10.1016/j.softx.2017.10.006
- K. F. Garrity et al. Pseudopotentials for high-throughput DFT calculations. Computational Materials Science 81, 446–452 (2014). doi:10.1016/j.commatsci.2013.08.053
- H. Lee et al. Electron–phonon physics from first principles using the EPW code. npj Computational Materials 9, 156 (2023). doi:10.1038/s41524-023-01107-3
- J.-J. Zhou et al. Perturbo: A software package for ab initio electron–phonon interactions, charge transport and ultrafast dynamics. Computer Physics Communications 264, 107970 (2021). doi:10.1016/j.cpc.2021.107970