全布里渊区第一性原理电子–声子耦合

完整三维硅 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,可直接在浏览器工作台中使用。

电子–声子耦合计算连接三个对象:位于 k\mathbf k 的电子态、波矢为 q\mathbf q 的晶格振动,以及位于 k+q\mathbf k+\mathbf q 的另一个电子态。准确计算每个对象只是其中一部分。当这些对象组合在一起,并从构建模型所用的网格点移向其他位置时,它们的相位、归一化、对称性和数值表示都必须保持一致。

我们已在 Mol2Mat 的金刚石结构硅四价带研究模型中完成这条计算链路。实现从自洽密度泛函理论出发,在完整的 8 × 8 × 8 电子网格与声子网格上构建声子动力学矩阵和复数电子–声子顶角,再通过 Fourier 插值计算连续 k/q\mathbf k/\mathbf q 坐标处的响应。

最终模型在原有 5 % EPC 判据下通过了全部 406 项既有回归查询、19 项独立电子 k\mathbf k 查询和 16 项独立联合 k/q\mathbf k/\mathbf q 查询。三组验证的最大误差分别为 4.2691 %、3.4066 % 和 4.5619 %。独立验证的预测值均在直接参考结果生成前冻结,整个验证过程中模型保持固定。

这一结果为指定的硅模型建立了可复现的三维动量空间第一性原理 EPC 计算。实现保留完整超软赝势顶角,确定响应所需的插值分辨率,并对弱耦合通道逐一验证。

从电子结构到复数散射振幅

电子结构计算采用 PBE 和 GBRV 硅超软赝势[1, 2],波函数与密度截断能分别为 20 和 100 Ha。金刚石结构原胞包含两个原子,对应六个笛卡尔位移坐标。电子子空间包含四条占据价带,由四个 Wannier 函数表示;计算采用非自旋极化设置,硅原子质量取 28.085 u。

对于原子 κ\kappa 沿 α\alpha 方向的笛卡尔扰动,我们保留一个复数 4 × 4 导数矩阵,记为 Vκα(k,q)V^{\kappa\alpha}(\mathbf k,\mathbf q),单位为 Ha bohr⁻¹。对质量加权动力学矩阵进行对角化后,声子极化矢量与零点位移共同将该导数转化为逐模耦合:

gmnν(k,q)=∑καℏ2Mκωqνeκα,ν(q)Vmnκα(k,q).g_{mn\nu}(\mathbf{k},\mathbf{q})= \sum_{\kappa\alpha} \sqrt{\frac{\hbar}{2M_\kappa\omega_{\mathbf q\nu}}} e_{\kappa\alpha,\nu}(\mathbf q) V^{\kappa\alpha}_{mn}(\mathbf k,\mathbf q).
gmnν(k,q)=∑καℏ2Mκωqνeκα,ν(q)Vmnκα(k,q).\begin{aligned}g_{mn\nu}(\mathbf{k},\mathbf{q})&=\\ &\sum_{\kappa\alpha} \sqrt{\frac{\hbar}{2M_\kappa\omega_{\mathbf q\nu}}} e_{\kappa\alpha,\nu}(\mathbf q) V^{\kappa\alpha}_{mn}(\mathbf k,\mathbf q).\end{aligned}

其中,mm、nn 标记矩阵采用的电子基,ν\nu 标记声子模式。存储的顶角采用 Wannier 基;对初态与末态电子端点分别施加幺正变换,即可得到相应的 Bloch 能带矩阵。将未加权的笛卡尔顶角与声子收缩分开保存,有助于明确插值过程和误差归因。

最终张量的形状为 (512, 512, 6, 4, 4),依次对应初态 k\mathbf k、声子 q\mathbf q、位移坐标、末态 Wannier 指标和初态 Wannier 指标。它包含 262,144 组 k/q\mathbf k/\mathbf q 组合、25,165,824 个复数导数元素,以及 512 个复数 6 × 6 动力学矩阵。这两个网格决定插值分辨率;自洽电子积分则采用另行记录的求积方案。

(a) 电子端点 $k$ 与 $k+q$(a) 电子端点 $k$ 与 $k+q$
(b) 声子波矢 $q$(b) 声子波矢 $q$
图 1. 硅的倒空间采样,坐标单位为 2π/a2\pi/a。线框为面心立方晶格的 Wigner–Seitz 第一布里渊区,灰点为 8³ 构建网格的 512 个节点。(a) 35 项前瞻性查询使用八个初态 k\mathbf k 位置(蓝色)和 35 个末态 k+q\mathbf k+\mathbf q 位置(红色)。(b) 两组独立联合验证 q\mathbf q 波矢。倒格矢等价坐标折叠到同一区内。
下载论文版图表 · PDF

这种拆分遵循 EPC 计算的物理结构。EPW [3] 采用局域表示,将电子、振动和耦合量插值到细密动量网格上;Perturbo [4] 则将这些第一性原理矩阵元用于输运与载流子动力学。在 Mol2Mat 中,此次里程碑是硅耦合基础模型的构建及独立数值验证,具体通过有限位移 JDFTx 响应和显式算符投影实现。

超软赝势顶角的一致构建

采用超软赝势时,电子本征问题包含重叠算符:

Hψ=εSψ.H\psi=\varepsilon S\psi.

原子位移同时改变 HH 和 SS。响应包含局域自洽势、非局域投影子、增广项和电子基的输运修正。两个电子端点也具有不同的 Bloch 相位。在对正负位移作差之前,每一项都必须采用相同的参考约定。

实现保存了投影矩阵 HrawH_{\mathrm{raw}}、H0H_0、SS、S0S_0、TT 和 HUSPPH_{\mathrm{USPP}}。按照所存储的算符约定,中心化并经输运修正的矩阵为

A(u)=Hraw(u)−v‾(u)S(u)−T(u)E0−E0T(u)†,E0=H0−v‾0I.\mathcal A(u)=H_{\mathrm{raw}}(u)-\overline v(u)S(u) -T(u)E_0-E_0T(u)^\dagger, \qquad E_0=H_0-\overline v_0 I.
A(u)=Hraw(u)−v‾(u)S(u)−T(u)E0−E0T(u)†,E0=H0−v‾0I.\begin{aligned} \mathcal A(u)&=H_{\mathrm{raw}}(u)-\overline v(u)S(u)\\ &\quad-T(u)E_0-E_0T(u)^\dagger,\\ E_0&=H_0-\overline v_0 I. \end{aligned}

平均势确定能量零点约定,TT 则提供双侧输运修正。组装前,代码检查恒等式 Hraw−TH0−H0T†=HUSPPH_{\mathrm{raw}}-TH_0-H_0T^\dagger=H_{\mathrm{USPP}},并检查参考 Hamiltonian 与重叠矩阵恒等式。原生投影使用已保存的自洽势,因此计算新的电子端点时可以复用已收敛响应。

中心差分通常沿选定单位方向采用 ±0.02 bohr 位移。原胞电子态嵌入 N 原胞超胞后,每个端点都带有 1/N1/\sqrt N 因子。因此,对投影矩阵块作差时,组装导数需恢复 N/(2δ)N/(2\delta) 因子。组合不同形状超胞的数据时,这一归一化尤其重要。

空间群输运同时作用于响应的各个部分:笛卡尔方向、原子置换、晶格平移偏移和两个电子端点的旋转。倒空间边界的重新标记同样属于算符定义的一部分。坐标越过布里渊区边界时,必须同步重标记平面波 G 指标;仅保留折返后的坐标数值还不够。

这些处理共同建立了一致的复数振幅。单个标量耦合强度会丢失检验这种一致性所需的相位和通道信息。

响应所需的插值分辨率

最终模型的形成过程揭示了一个重要区别:数值重构准确,并不意味着物理表示具有足够的表达能力。

早期 4³ 电子 Wannier–Fourier 表示能够以数值精度重现全部 64 个构建点,但在其他已保存 k\mathbf k 点上仍有较大误差。于是,我们检验了同一固定 Fourier 空间中的任意一组系数,能否足够准确地重现已有响应。

对每组电子矩阵,令 AA 表示在 456 个已保存观测点上计算固定的 64 个基函数,YY 表示对应的直接响应,Π\Pi 表示到 AA 列空间的投影。残差可以作正交分解:

Y−AC4=(Y−ΠY)+(ΠY−AC4).Y-AC_4=(Y-\Pi Y)+(\Pi Y-AC_4).

第一项位于所选表示空间之外,调整空间内的系数无法消除它。通过薄 QR 分解进行投影,得到完整笛卡尔顶角的最小整体相对误差为 5.1642 %,最低声学模组为 8.6529 %,均超过 5 %。如果每个点的误差都不超过 5 %,整体 Frobenius 范数比值也应不超过 5 %。因此,这个固定表示无法在全部 456 个观测点上满足该条件。

原始整体误差与固定空间误差下界原始整体误差与固定空间误差下界
图 2. q\mathbf q = (1/4,0,0) 处的有限表示空间检验。灰条采用原始 4³ 系数,蓝条表示同一固定空间内可达到的最小误差。这里比较的是相同 456 个观测点上的整体矩阵误差,模式从一开始编号,虚线表示 5 %。该计算识别出指定 Wannier 规范和基下的表示能力限制。
下载论文版图表 · PDF

我们还通过独立的有限位移对照,检查直接参考值是否足以解释历史上的大误差。在被审查点,将位移从 0.02 bohr 减半至 0.01 bohr,只使问题模式的参考矩阵变化 0.02742 %;冻结预测的误差则从 21.3081 % 变为 21.2993 %。这一配对对照将后续研究指向表示空间和采样分辨率。

电子网格提高到 8³ 后,五个固定 q\mathbf q 回归查询的最大误差从 18.1314 % 降至 2.0051 %。另外 392 个固定 q\mathbf q 回归点也全部通过,最大误差为 4.2691 %。

声子波矢方向仍需要独立提高分辨率。在 k8 × q4 下,q\mathbf q = (2,3,5)/12 处八个非轴向查询的最大误差仍为 14.7329 %;完成 q8 后,同一组查询降至 3.6755 %。轴向 q\mathbf q = (1/3,0,0) 查询也从 24.1940 % 改善至 3.2333 %。这里的坐标均为原胞倒格基下的分数坐标。

(a) 固定 $q$,五个 $k$ 点($n=5$)(a) 固定 $q$,五个 $k$ 点($n=5$)
(b) 固定 $q$,已有 $k$ 点($n=392$)(b) 固定 $q$,已有 $k$ 点($n=392$)
(c) $q=(2,3,5)/12$ ($n=8$)(c) $q=(2,3,5)/12$ ($n=8$)
(d) $q=(1/3,0,0)$ ($n=1$)(d) $q=(1/3,0,0)$ ($n=1$)
图 3. 配对回归比较区分了电子网格与响应网格分辨率的作用。每个子图均保留相同的查询集合、直接参考和评分定义,虚线为原有 5 % 判据。这些已检查点用于诊断改进来源,独立验证见后文。
下载论文版图表 · PDF

这一演进给出了计算困难的具体来源:电子 k\mathbf k 分辨率和声子 q\mathbf q 分辨率约束同一振幅的不同部分。提高一个坐标空间的分辨率后,另一个空间仍可能成为主导误差来源。

以 33 项新增自洽响应补齐 q8

q8 扩展复用了已有 q4 响应和轴向数据。十二组八原胞响应提供新增非轴向信息,一组保留的 24 原胞轴向响应再提供 16 个 q\mathbf q 节点。结合对称性,这些来源补齐了全部 512 个 q\mathbf q 节点:64 个保留的 q4 节点、16 个轴向节点和 432 个新增节点。

最初的正负位移计划包含 36 个 SCF 状态。晶体对称性将三个负位移状态映射到对应的正位移状态,把实际新增计算减少至 33 项。这些操作经过实际原子映射、带权电子积分网格和 FFT 网格的检查,再由已保存的正位移响应通过核验后的变换生成对应状态。

33 项自洽计算全部完成。每项输入均包含 16 个原子、32 条电子能带和 432 个显式带权电子积分点。收敛规则要求连续两个周期达到原有 10⁻⁹ SCF 残差阈值。后续用于独立联合验证的另外八项 SCF 属于单独的参考计算批次。

组装完成的笛卡尔顶角沿六个动量网格轴进行 Fourier 变换。连续坐标查询通过 k\mathbf k 和 q\mathbf q 的相位因子计算相应实空间系数:

V(k,q)=∑Re,Rpe2πi(k⋅Re+q⋅Rp)V~(Re,Rp).V(\mathbf k,\mathbf q)= \sum_{\mathbf R_e,\mathbf R_p} e^{2\pi i(\mathbf k\cdot\mathbf R_e+\mathbf q\cdot\mathbf R_p)} \widetilde V(\mathbf R_e,\mathbf R_p).

实现中的相位因子对等价的最近周期镜像取平均。电子镜像的选择取决于 Wannier 中心之差;响应镜像的选择取决于原子相对于末态 Wannier 中心的位置;晶格度量决定哪些镜像最近。这些固定规则在保留粗网格相位恒等式的同时,确定节点之间的插值。动力学矩阵则采用相应的原子中心 q\mathbf q 空间构建。

声子组装中还出现了一个细微的对称性问题:简并模式极小的数值劈裂,会让原有模组比较器无法找到维数相同的匹配。解决方法是直接对每个新增动力学矩阵施加晶体的小群对称性。矩阵 Hermitian 化后,通过 Reynolds 平均施加所有保持 q\mathbf q 不变的操作:

Dsym(q)=1∣Gq∣∑s∈GqWs(q)DH(q)Ws(q)†.D_{\mathrm{sym}}(\mathbf q)= \frac{1}{|G_{\mathbf q}|}\sum_{s\in G_{\mathbf q}} W_s(\mathbf q)D_H(\mathbf q)W_s(\mathbf q)^\dagger.

原子位移表示 WW 包含晶格偏移相位。在检查的 84 个 q\mathbf q 点中,最大相对矩阵修正为 0.02246 %,最大协变性残差为 3.22 × 10⁻¹⁶。所有公共点的模组维数均匹配;原始 q4 数值精确保留,原有分组容差保持固定。

对称性减少了所需响应计算的数量,也恢复了逐模比较所需的简并结构。

冻结预测与独立直接验证

最终评估分别检验三类问题。

既有回归。 406 项查询由 392 个已保存的固定 q\mathbf q 电子点、五个额外固定 q\mathbf q 对照、八个非轴向联合查询和一个轴向离网格查询组成。其直接结果在开发期间已经获得。最终 q8 模型通过全部 406 项。

独立电子 k\mathbf k。 十九项新增电子查询经过与已暴露对称轨道的检查,覆盖晶体操作、时间反演和逆散射。预测在新增直接投影组装前冻结。由于 q\mathbf q 位于已有声子节点,这组验证检验新电子端点处的插值;十九项全部通过。

独立联合 k/q\mathbf k/\mathbf q。 十六项查询覆盖两组非轴向声子波矢 q\mathbf q = (1,2,3)/12 和 q\mathbf q = (1,4,7)/12,每组包含八项电子查询。八项新增正负位移 SCF 响应提供独立参考。独立参考采用 208 个固定势投影块构建。十六项查询均通过与参考生成前冻结预测的比较。

模型张量、Fourier 规则和 5 % 阈值始终固定。公共来源对照另外使用 0.5 % 判据,检查新参考组装与既有来源的一致性。

表 1. 独立与回归查询的最大相对误差
验证集合 通过查询数 仅电子部分最大误差 仅声子部分最大误差 联合最大误差
既有回归 406 / 406 4.2691 % 0.6853 % 4.2691 %
独立电子 k\mathbf k 19 / 19 3.4066 % 0.00001261 % 3.4066 %
独立联合 k/q\mathbf k/\mathbf q 16 / 16 4.5619 % 1.9479 % 4.5619 %

仅电子比较将预测顶角与参考动力学矩阵配对;仅声子比较将参考顶角与预测动力学矩阵配对;联合比较同时使用两项预测。这些误差来自同一耦合计算中的受控替换,不能直接相加。

(a) 回归验证($n=406$)(a) 回归验证($n=406$)
(b) 独立电子 $k$($n=19$)(b) 独立电子 $k$($n=19$)
(c) 独立联合 $k/q$($n=16$)(c) 独立联合 $k/q$($n=16$)
图 4. 每组通过验证的查询均包含在图中。横坐标为比较器返回的单点最大误差,纵坐标为误差不超过该值的查询比例。仅电子与联合曲线高度重叠,各组最大值相同,但部分联合查询的分项评分不同。独立电子 k\mathbf k 组的声子误差接近零,源于其 q\mathbf q 位于已有声子节点。竖直虚线表示 5 %。
下载论文版图表 · PDF

这一结果结合了广泛的回归覆盖和 35 项前瞻性独立查询。Full-BZ 指连续三维 k/q\mathbf k/\mathbf q 表示,数值精度由这些明确定义的测试集度量。

弱耦合通道的逐模验证

单个 EPC 范数不足以完整检验顶角。两个复数矩阵可以具有几乎相同的范数,却在电子通道空间中具有不同方向。较小的声学响应也可能来自更大原子贡献之间的抵消,从而放大最终模式中的相对误差。

因此,比较器保留完整原始顶角,并逐一检查声子模组,比较同规范下的矩阵差、耦合范数和奇异值谱。在简并模组内,将模式矩阵堆叠后得到的通道谱对模式基的幺正混合不变。预测与参考声子子空间通过最大化重叠的等维双射进行匹配;此外,也使用各自计算得到的模式评估耦合范数与通道谱。

例如,同规范矩阵误差定义为

ϵV=∥Vpred−Vref∥F∥Vref∥F.\epsilon_V=\frac{\|V_{\mathrm{pred}}-V_{\mathrm{ref}}\|_F} {\|V_{\mathrm{ref}}\|_F}.

单点评分取完整顶角误差,以及适用模组的矩阵、范数和通道谱误差中的最大值。无法获得有效模式匹配时,该点不能通过。零频或不稳定模式保留带符号本征值和诊断状态,不人为赋予正频率。

最大独立联合误差出现在 k\mathbf k = (1,2,3)/12、q\mathbf q = (1,4,7)/12。其最弱模式的直接耦合范数为 74.1102 meV,预测为 73.6286 meV,范数误差仅 0.6498 %;但复数矩阵误差为 4.5619 %,通道谱误差为 2.1769 %。该模式只占原始笛卡尔顶角范数平方的 0.5431 %。其原子抵消比,即净模式振幅与两个原子振幅之和的比值,为 0.1333。

总强度图会给这一苛刻模式很低的权重,逐模比较则把它纳入报告精度的核心。

(a) 耦合范数对照(a) 耦合范数对照
(b) 矩阵与范数误差(b) 矩阵与范数误差
图 5. 十六项联合查询产生 96 个稳定、非简并模组。(a) 以 meV 表示的 Frobenius 耦合范数,两种顶角均使用固定的直接参考声子模式和频率。(b) 同一组数据的相对范数误差与相对模式矩阵误差。颜色与标记形状区分两组 q\mathbf q 波矢。(a) 的对角线表示完全一致;(b) 显示仅比较范数会低估部分矩阵误差。
下载论文版图表 · PDF

该查询的声子位移图提供了互补的实空间视角。最低模频率在直接计算中为 4.6445 THz,在冻结预测中为 4.6260 THz,两者本征矢重叠的模平方为 0.999994。振动图样的高度一致性与分项评分共同表明,较大的剩余误差位于电子顶角。

(a) 直接参考: 4.6445 THz(a) 直接参考: 4.6445 THz
(b) 冻结预测: 4.6260 THz(b) 冻结预测: 4.6260 THz
图 6. 最大独立联合误差查询处的最低模式,在包含 16 个硅原子的 2 × 2 × 2 原胞显示窗口中呈现。位移遵循 Re⁡[eκexp⁡(2πiq⋅R−iϕ)]\operatorname{Re}[e_\kappa\exp(2\pi i\mathbf q\cdot\mathbf R-i\phi)],ϕ=0\phi=0。两个本征矢采用相同的相位约定与相机方向。显示位移幅度为 0.20 Å,箭头再统一放大四倍;淡色轮廓表示平衡位置。底层计算保留原始归一化本征矢和频率。
下载论文版图表 · PDF

配套的交互式三维视图支持旋转、全部六个模式的选择、相位调整,并可导出 PNG 和 GLB。

这些检验针对后续散射计算所需的矩阵元。物理可观测量的收敛还取决于电子能量窗口、相空间积分、占据和能量守恒处理。保留完整复数顶角,为后续计算提供了比标量拟合摘要更扎实的基础。

八核心并行性能

本次硅 EPC 计算采用八核心预算。自洽响应和固定势投影可按独立任务并行处理;已有自洽势可以用于新的电子端点,减少重复计算。

我们在同一台机器上比较了两种核心分配。测试使用相同的八个已保存 SCF 状态,每个状态继续计算两个周期,耗时包含启动、检查点保存和结果输出:

表 2. 相同八个 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。这组结果给出了明确的选择依据:内存充足时,可以增加同时运行的独立任务,以缩短等待时间。

线程分配也随剩余任务数调整:最多同时运行八个单线程任务,剩余四项或三项时每项使用两线程,剩余两项时每项四线程,最后一项可使用八线程。这样,批次接近完成时仍能利用可用核心。

连续动量查询与后续计算

这套模型支持连续三维 k/q\mathbf k/\mathbf q 坐标处的顶角、声子模式和逐模耦合查询。研究者可以检查具体散射通道,将预测与直接参考逐项比较,并使用保存的矩阵复现误差。已有 441 项插值重放记录以 9.86 × 10⁻¹⁶ 的最大相对差异重现了冻结预测。

同一冻结模型现已接入浏览器的 FP64 WASM 计算工作流,用户可以执行连续 k/q\mathbf k/\mathbf q 查询、查看声子模式与耦合、比较不同波矢,并导出复数矩阵。后续将通过外部矩阵对照、系统网格收敛和更大电子子空间检验其适用范围,在此基础上计算由这些顶角决定的散射率和输运性质。

在浏览器中使用

在 DFT 单点计算中选择电子–声子耦合,输入以原胞倒格矢为基底的电子波矢 k\mathbf k 和声子波矢 q\mathbf q,即可比较六个振动模式及不同波矢下的耦合强度。选中的模式与晶格振动动画联动;完整复数耦合矩阵可在详情中查看和导出。

浏览器使用上述 8³ 电子网格与声子网格构建的硅模型。首次使用需下载约 384.5 MiB 的模型数据,并校验文件大小和 SHA-256。下载完成后,查询由本机的 WebAssembly 内核执行。该内核由 C 编译而来,以双精度浮点数(FP64)计算复数 Fourier 插值,并复用相同 q\mathbf q 对应的中间结果。计算在独立 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

  1. R. Sundararaman et al. JDFTx: Software for joint density-functional theory. SoftwareX 6, 278–284 (2017). doi:10.1016/j.softx.2017.10.006
  2. 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
  3. 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
  4. 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