来源论文: https://arxiv.org/abs/2606.12782v1 生成时间: Jun 12, 2026 01:09

冲向热力学极限：周期性CCSD凝聚能与带隙的高密度布里渊区抽样深度解析

0. 执行摘要

在现代凝聚态物理与材料化学的计算模拟中，如何实现接近“化学精度”（~1 kcal/mol 或 0.04 eV）的电子结构预测始终是一个核心难题。尽管密度泛函理论（DFT）因其高计算效率而成为固态计算的主力军，但其在处理强相关效应、范德华力以及精确带隙预测时，常常受限于交换相关泛函的固有缺陷。相比之下，波函数理论（WFT）中的耦合簇理论（Coupled-Cluster with Singles and Doubles, CCSD）作为高精度量子化学的“金标准”，具有非微扰、尺寸一致性（Size-Extensivity）等优良性质，具备解决上述难题的理论潜力。

然而，将周期性 CCSD 及其激发态扩展（EOM-CCSD）应用于固态体系面临着极其严苛的计算挑战。其主要瓶颈源自布里渊区抽样（$k$ 点网格）的缓慢收敛性。由于周期性 CCSD 的计算量随 $k$ 点数量 $N_k$ 呈四次方甚至更高次方比例增长，传统的计算程序往往只能处理 $N_k \le 4^3 = 64$ 的小网格，导致有限尺寸效应（Finite-Size Errors, FSE）无法彻底消除，难以可靠地外推至热力学极限（Thermodynamic Limit, TDL）。

近期，来自 Emory 大学、Flatiron 研究所和哥伦比亚大学的研究团队（Shuhang Li, Huanchen Zhai, Francesco A. Evangelista 和 Timothy C. Berkelbach）在这一领域取得了里程碑式的进展。他们开发了一种全新的、基于分布式内存的并行化周期性 CCSD 和 EOM-CCSD 计算程序。该程序高度集成了空间群对称性、高效率的张量收缩优化、按需（On-the-fly）积分生成以及分阶段内存管理技术，能够在多达 12 个节点（共 1152 个 CPU 核心）上实现高度线性的并行扩展。利用该技术，研究团队首次在双 $\zeta$（DZ）基组下实现了高达 $6^3 = 216$ 个 $k$ 点的布里渊区抽样，并在三 $\zeta$（TZ）和四 $\zeta$（QZ）基组下分别实现了 $5^3$ 和 $4^3$ 的抽样。

基于这一高性能并行计算平台，他们系统研究了八种具有代表性的半导体和绝缘体（MgO, LiCl, LiF, LiH, BN, BP, Si, C）的凝聚能与带隙，将有限尺寸误差控制在 0.1 eV 以内，首次给出了 CCSD 理论级别下无偏的、收敛的基准数据（Benchmark）。结果表明：

1. 对于凝聚能，CCSD 相比实验值平均低估约 0.1–0.2 eV，而 MP2 则平均高估 0.24–0.30 eV；
2. 对于基本带隙，EOM-CCSD 相比实验值平均高估约 0.4 eV，但在岩盐结构（Rock-salt）体系（如 MgO, LiH, LiF, LiCl）中存在由于单激发表征度（$n_1^{IP}$）不足导致的系统性偏差；
3. 评估了以往基于小 $k$ 网格的外推方案，揭示了 GW 辅助外推在共价半导体中的局限性，并推荐了更为稳健的“二参数（AB）外推模型”与“复合间接带隙外推法”。

本博客将对该研究的核心科学问题、理论基础、技术难点、方法细节、计算基准数据以及分布式内存并行代码的设计架构进行全面且深度的技术剖析。

1. 核心科学问题、理论基础、技术难点与方法细节

1.1 核心科学问题：布里渊区抽样与热力学极限（TDL）

在周期性边界条件（PBC）下，固态体系的计算是在倒易空间（Reciprocal Space）中进行的，需要对第一布里渊区（First Brillouin Zone）进行离散化抽样，即引入 $k$ 点网格。理论上，只有当 $k$ 点网格趋于无穷密（即 $N_k \to \infty$）时，计算结果才能真正对应实际的无限大晶体，这一极限状态被称为热力学极限（TDL）。

对于 mean-field 方法（如 Hartree-Fock 或 DFT），由于其单粒子近似特性，通常可以通过使用极其密集的 $k$ 网格（如 $12^3$ 或 $16^3$）来轻松消除有限尺寸误差。然而，对于电子相关方法（Correlated Methods），由于双电子相互作用算符在倒易空间中的非局部特性，其有限尺寸效应的收敛极其缓慢。具体而言：

• 基态相关能量（如 MP2 和 CCSD 相关能）：其有限尺寸误差通常表现为与体积成反比，即以 $N_k^{-1}$ 的渐近速度收敛。
• 激发态能（如带隙 $E_{gap}$）：由于外加电子或空穴与自身晶胞图像之间的长程库仑相互作用被介电常数筛选，带隙的有限尺寸误差收敛更为缓慢，其主导项与体积的三分之一次方成反比，即以 $N_k^{-1/3}$ 的渐近速度收敛。

由于高精度相关方法的计算复杂度极高（例如周期性 CCSD 的时间复杂度为 $\mathcal{O}(N_{atom}^6 N_k^4)$），以往的研究只能局限于 $N_k \le 3^3$ 或 $4^3$ 的网格。在如此稀疏的网格下，体系尚未进入渐近收敛区间，直接应用渐近外推公式往往会引入显著的数值不确定性。因此，如何在保障大基组计算的同时，突破计算瓶颈以获取 $N_k = 5^3$ 和 $6^3$ 的稠密网格数据，是确定固态 CCSD 本征物理精度、验证外推方案可靠性的关键科学问题。

1.2 理论基础：周期性 CCSD 与 EOM-CCSD

1.2.1 周期性 CCSD 理论

Coupled-Cluster 理论通过指数拟设将多电子波函数参数化：

其中 $|\Phi_0\rangle$ 为参考态（通常为 Hartree-Fock 决定式），$\hat{T} = \hat{T}_1 + \hat{T}_2$ 为簇算符（Cluster Operators）。在周期性晶体动量表象中，空间轨道带有晶体动量标记 $k$。簇算符定义为：

其中，带有撇号的求和号 $\sum^{\prime}$ 施加了严格的晶体动量守恒约束：

这里 $G$ 是倒易点阵矢量。这一守恒律极大地限制了非零 $t_2$ 振幅的规模：原本包含四个独立 $k$ 点的张量被降维为仅有三个独立的 $k$ 点，这在代码编写中可以通过动量守恒直接进行索引寻址。

1.2.2 周期性 EOM-CCSD 理论

为了计算体系的激发态（如电离能 IP 和电子亲和能 EA），引入运动方程（Equation-of-Motion）框架。激发态波函数表示为：

通过对相似变换后的哈密顿量 $\bar{H} = e^{-\hat{T}} H e^{\hat{T}}$ 进行对角化，可以求解激发能。对于带隙计算，我们特别关注：

1. 电离势（IP）：对应在 $N$ 电子系统基态上移除一个电子。IP-EOM-CCSD 在 1h（单空穴）和 2h1p（双空穴-单粒子）激发空间中求解：

   $$\hat{R}^-(k) = \sum_{i} r_{i_{k}} \hat{a}_{i_{k}} + \frac{1}{2} \sum_{k_i k_j k_a}^{\prime} \sum_{ija} r_{i_{k_i} j_{k_j}}^{a_{k_a}} \hat{a}_{a_{k_a}}^{\dagger} \hat{a}_{j_{k_j}} \hat{a}_{i_{k_i}}$$

2. 电子亲和能（EA）：对应在外加一个电子。EA-EOM-CCSD 在 1p（单粒子）和 2p1h（双粒子-单空穴）激发空间中求解：

   $$\hat{R}^+(k) = \sum_{a} r^{a_{k}} \hat{a}_{a_{k}}^{\dagger} + \frac{1}{2} \sum_{k_i k_a k_b}^{\prime} \sum_{iab} r_{i_{k_i}}^{a_{k_a} b_{k_b}} \hat{a}_{a_{k_a}}^{\dagger} \hat{a}_{b_{k_b}}^{\dagger} \hat{a}_{i_{k_i}}$$

基本带隙（Fundamental Band Gap）则定义为：

1.3 技术难点与应对策略

1. 极端的内存开销：周期性 CCSD 计算中，最庞大的数据结构是四脚双电子积分和 $t_2$ 振幅张量。对于包含 $N_{k}$ 个 $k$ 点、每个 $k$ 点有 $N_{occ}$ 个占据轨道和 $N_{vir}$ 个虚拟轨道的体系，振幅张量 $t_{ij}^{ab}(k_i, k_j, k_a)$ 占用的空间随 $N_k^3 N_{occ}^2 N_{vir}^2$ 剧烈增长。在常规双 $\zeta$（DZ）甚至三 $\zeta$（TZ）基组下，直接在内存中存储这些张量会轻易撑爆单机内存。

   • 应对策略：采用**三中心密度拟合（Density Fitting, DF）**技术，仅在内存中存储三中心 DF 积分 $L_{pq}^Q(k_p, k_q)$。在 CCSD 的张量收缩迭代中，涉及大量虚拟轨道索引的四中心积分（如最昂贵的 $VOVV$ 和 $VVVV$ 块）通过 DF 积分在线（On-the-fly）实时重构并直接参与收缩，避免了显式构建与存储。

2. 空间群对称性导致的负载不均衡（Load Imbalance）：利用晶体空间群对称性（空间群阶数为 $G$），可以将整个布里渊区中的 $k$ 点四元组 $(k_i, k_j, k_a, k_b)$ 约化为不可约布里渊区（Irreducible Brillouin Zone, IBZ）内的四元组，使得计算量和存储量理论上可以下降为原来的 $1/G$。但在实际分布式计算中，由于一部分 $k$ 四元组本身就在 IBZ 内（直接读取，开销极小），而另一部分需要通过旋转矩阵从对称性相关的 IBZ 块中旋转得到（需要进行 $O(N_{occ}^2 N_{vir}^3)$ 的轨道旋转，开销极大），若采用平分 $k$ 空间任务的简单划分策略，会导致严重的计算负载失衡。

   • 应对策略：实现了一种基于干运行（Dry-run）和贪心最长处理时间优先（LPT）算法的动态负载均衡调度器。在正式计算前，先模拟评估各进程的几何旋转操作频数，为其分配“权重”，再利用最小堆（Min-heap）确保各进程的工作负载（FLOPs 估算值）与任务数量达到双重均衡。

3. 带隙收敛的 $N_k^{-1/3}$ 极慢收敛行为：如前所述，带隙的有限尺寸误差以 $N_k^{-1/3}$ 级数收敛，而直接带隙与间接带隙的收敛物理机制还存在差异。

   • 应对策略：发展了更为精细的多元渐近拟合外推公式，并针对间接带隙提出了“复合外推法”（详见 1.4 节）。

1.4 方法细节：高精度外推方案

1.4.1 基态能量外推

Hartree-Fock 基态能量采用三参数外推模型，引入二阶修正项以加速收敛：

而对于 MP2 和 CCSD 的相关能（Correlation Energy），由于相关空穴的局部性较好，采用更稳健的双参数 leading-order 外推：

1.4.2 带隙外推的三种模型

带隙随 $k$ 点网格密度的变化可通过以下通用公式描述：

根据保留项的不同，论文系统评估了三种拟合模型：

• Model ABC（三参数模型）：保留所有三项（$N_k^{-1/3}, N_k^{-2/3}, N_k^{-1}$）。这是最完备的模型，但需要极其稠密的 $k$ 点数据以防止过拟合。
• Model AB（二参数模型）：仅保留前两项（$N_k^{-1/3}, N_k^{-2/3}$）。
• Model A（单参数模型）：仅保留 leading-order 项（$N_k^{-1/3}$）。这在过去的低抽样密度研究中被广泛采用，但本工作证实，当抽样网格不够密时（如仅到 $4^3$），由于忽略了高阶项，Model A 会系统性地低估带隙 0.3–0.5 eV。

1.4.3 间接带隙的“复合外推法”（Composite Approach）

对于间接带隙体系（如 BN, BP, Si, C），价带顶（VBM）和导带底（CBM）位于倒易空间的不同 $k$ 点。计算这类带隙通常需要对 $k$ 网格进行整体平移（Shift），以使网格中刚好包含非 $\Gamma$ 点的 VBM 或 CBM。然而，$k$ 网格的平移会极大地破坏倒易空间的点群对称性，导致 IBZ 中的独立 $k$ 点大幅增加，计算量呈几何级数暴涨，使得大 $k$ 网格计算（如 $5^3$ 或 $6^3$）无法实施。

为了解决该痛点，论文提出了复合外推法（Composite Approach），其物理拟设如下：

• 第一项 $E_{direct, TDL}^{EOM}$：在价带顶 $k_{VBM}$ 处的直接带隙（对于本文研究的体系，其 VBM 均位于 $\Gamma$ 点）。由于计算直接带隙无需平移网格，能完美保留空间群对称性，因此可以使用最稠密的 $6^3$ 网格和完备的 Model ABC 进行外推，获取极高精度的直接带隙 TDL 值。

• 第二项 $\Delta E_{CB, TDL}^{EOM}$：导带底的能量差值偏移，即：

  $$\Delta E_{CB} = E_{CB}(k_{CBM}) - E_{CB}(k_{VBM})$$

  物理研究表明，导带偏移量 $\Delta E_{CB}$ 随 $k$ 点网格的收敛速度极快（在 $N_k \ge 3^3$ 之后，其波动已小于 0.05 eV）。因此，我们可以使用较小的 $k$ 网格（如 $3^3$ 或 $4^3$）计算平移网格下的 $\Delta E_{CB}$，将其作为精确的修正项直接叠加上去。这一巧妙的分解方法既避开了对称性破缺带来的计算灾难，又保障了间接带隙的 TDL 预测精度。

2. 关键基准体系、计算数据与性能分析

2.1 研究的基准体系参数

研究团队精选了八种具有代表性的半导体和绝缘体，涵盖了不同的晶体结构：

• 金刚石结构：C, Si
• 闪锌矿结构（Zinc-Blende）：BN, BP
• 岩盐结构（Rock-Salt）：MgO, LiCl, LiF, LiH

下表汇总了各体系的实验晶格常数 $a$、价带顶（VBM）与导带底（CBM）在倒易空间中的相对坐标位置：

2.2 凝聚能基准计算结果分析

凝聚能（Cohesive Energy）定义为隔离原子总能与晶体单胞平均原子能量之差。下表汇总了利用 GTH 赝势（PP）及全电子冰冻核心近似（AE-FC）在不同理论级别（HF, MP2, CCSD）下得到的、外推至热力学极限与完备基组极限（TDL/CBS）的凝聚能，并与经零点能修正后的实验值进行了对比：

物理结果剖析

1. Hartree-Fock 的严重低估：HF 级别计算表现出极大的系统性偏差，平均绝对误差（MAE）高达 1.50 eV，这主要是由于 HF 忽略了电子动态相关效应，导致其无法正确描述晶体中的分散力与极化相互作用。
2. MP2 与 CCSD 的精细角逐：电子相关能的引入彻底改善了结果。MP2 表现出系统性高估（MAE 为 0.24–0.30 eV），而CCSD 表现出轻微的低估（MAE 降至 0.14 eV）。这表明 CCSD 作为高阶非微扰方法，有效克服了 MP2 理论中由于微扰过大引起的“过相关”（Overcorrelation）问题。
3. 赝势效应（PP vs AE-FC）：使用 GTH 赝势带来的系统误差非常微弱，绝大多数体系的 PP 凝聚能与 AE-FC 的差异在 0.05 eV 以内。唯独 BN 和 C 的偏差达到 0.1–0.2 eV，这与其极强的价电子局部化以及复杂的核电荷极化密切相关。

2.3 基本带隙基准计算结果分析

通过对 EOM-CCSD 计算结果施加 TDL 外推及基组不完备修正（CBS），研究团队获得了迄今为止最精准的固态 EOM-CCSD 带隙数据。下表对比了赝势（PP）和全电子冰冻核心（AE-FC）下的 EOM-CCSD 带隙与零点振动重整化（ZPR）修正后的实验带隙：

带隙物理本质分析

1. 系统性高估：在排除了有限尺寸误差和基组误差后，EOM-CCSD 带隙表现出系统性的轻微高估，MAE 约为 0.4 eV。这说明在不含三激发算符的情况下，EOM-CCSD 略微低估了体系在高激发状态下的极化自能筛选效应。
2. 岩盐结构体系的局限性：在 MgO, LiCl, LiF 和 LiH 等典型强离子晶体（岩盐结构）中，EOM-CCSD 的高估最为显著（最高达 0.97 eV）。
   • 深层物理根源：研究团队通过计算单激发表征度 $n_1^{IP}$（即电离算符中 1-hole 部分的振幅平方和）揭示了这一现象。在共价体系（C, Si, BN, BP）中，$n_1^{IP} \ge 95\%$，此时 EOM-CCSD 能完美适用。而在岩盐体系中，$n_1^{IP}$ 显著降低至 $93\% - 94\%$，这意味着双激发（如 shake-up 状态）和更高阶的电子相关效应在这些离子晶体中占有不可忽视的地位，必须引入三激发算符（如 EOM-CCSD(T)）才能彻底消除偏差。

2.4 并行计算性能展示

为了检验分布式代码的扩展性，研究团队以包含过渡金属原子的复杂体系——**金红石型二氧化钛（Rutile $TiO_2$）**为测试对象。计算采用 $2 \times 2 \times 3$ 的 $k$ 网格以及高精度的 GTH-DZVP-MOLOPT-SR 基组，并在配置有双路 96 核 AMD Genoa 处理器、1.5 TB 内存的超算节点上运行。下图展示了代码在 1 到 12 个节点（96 至 1152 核）范围内的强扩展性（Strong Scaling）表现：

  Speed-up (relative to 1 node)
   12 +-----------------------------------------------------------+ 
      |                                                   O Eff=0.94
   10 |                                          /--------/       | 
      |                                 O Eff=0.95                |
    8 |                        /--------/                         | 
      |               O Eff=0.98                                  |
    6 |      /--------/                                           | 
      |                                                           |
    4 |  /---/                                                    | 
      |                                                           |
    2 |O Eff=1.00                                                 | 
      +-----------------------------------------------------------+
      1               4               8              12
                            Number of Nodes

• 理想的线性加速比：在最大规模的 12 节点计算中，代码斩获了 11.31 倍的加速比，并行效率（Parallel Efficiency）高达 94%。
• 耗时对比：传统的 PySCF 串行无对称性版本单次 CCSD 迭代耗时高达 41,442 秒；利用空间群对称性加速后降至 27,673 秒；而本工作通过极限优化张量收缩路径，单节点耗时缩短至 17,511 秒。当使用 12 节点并开启多任务并行调度（12 节点 $\times$ 4 任务/节点 $\times$ 24 线程/任务）时，单次 CCSD 迭代时间骤降至 437 秒，实现了近百倍的综合效率提升！

3. 代码实现细节、算法架构与复现指南

3.1 核心算法架构设计

周期性 CCSD 计算的最核心步骤在于反复迭代求解幅度方程，涉及大量形如 $W_{ab}^{cd} = \sum_{ij} t_{ij}^{cd} \langle ab||ij\rangle$ 的多维张量收缩（Tensor Contractions）。整个分布式内存并行程序的架构设计如下图所示：

                                +--------------------------+
                                |    Input Geometry & k    |
                                +-------------+------------+
                                              |
                                              v
                                +--------------------------+
                                |  Symmetry Analysis & G   |
                                |   Construct IBZ Mesh     |
                                +-------------+------------+
                                              |
                                              v
                                +--------------------------+
                                |     Dry-Run Estimator    | <--- Weights for symmetry 
                                |  Generate Task Schedule  |      rotations estimated
                                +-------------+------------+
                                              |
                                              v
                                +--------------------------+
                                |  Constrained LPT Greedy  |
                                |  Load Balancing Engine   |
                                +-------------+------------+
                                              |
                                 Distributed across MPI
                                              |
                      +-----------------------+-----------------------+
                      |                                               |
                      v                                               v
            +-------------------+                           +-------------------+
            |    MPI Rank 0     |                           |    MPI Rank N     |
            |                   |                           |                   |
            | Store IBZ-t2-amp  |                           | Store IBZ-t2-amp  |
            | Store 3-index DF  |                           | Store 3-index DF  |
            +---------+---------+                           +---------+---------+
                      |                                               |
                      +-----------------------+-----------------------+
                                              |
                                     K-Phase Buffer Split
                                              |
                                              v
                                +--------------------------+
                                |   MPI_Alltoallv Comm     |
                                |   On-the-fly DF Contract |
                                +--------------------------+

3.1.1 分布式空间群轨道的对称性旋转（Symmetry Rotations）

如前所述，为了节省存储空间，程序仅在内存中存储不可约布里渊区（IBZ）内的 $t_2$ 振幅块。当收缩过程中需要用到非不可约布里渊区（out-of-IBZ）的块时，程序会执行轨道旋转：

其中 $U^{\hat{R}}(k)$ 为单粒子晶体轨道的变换矩阵。这一步骤可以通过四步连续的矩阵乘法高效完成，每一步的时间复杂度被严格控制在 $\mathcal{O}(N_{occ}^2 N_{vir}^3)$。这一机制使得存储空间缩小了 $G$ 倍（对于高对称性晶体，如 MgO $G=48$），代价是增加了对称性重建的计算时间。为此，干运行估算器大显身手。

3.1.2 负载均衡干运行（Dry-Run）引擎

为了避免计算旋转操作造成的严重负载失衡，干运行阶段会遍历所有收缩任务。每个任务所需的对称旋转变换次数被定量记作“权重”（Weight）。随后，利用贪心的 LPT（Longest Processing Time First）策略进行分配：

1. 将所有任务按“权重”从大到小排序；
2. 依次将任务分派给当前累计权重之和最小的进程；
3. 限制每个进程分得的任务总数偏差不超过 1，从而兼顾了任务数均衡与 FLOPs 计算开销均衡。

3.1.3 基于 $K$ 阶段缓存分裂（$K$-Phase Buffer Splitting）的通信控制

为了防止在 MPI 通信中因瞬间缓冲巨大张量而导致物理内存溢出，本工作提出了一种创新的分阶段通信策略。每次迭代中，将该进程的任务拆分为 $K$ 个独立阶段。在每个阶段：

1. 仅通过 MPI_Alltoallv 向目标进程请求该阶段所需的特定 $t_2$ 振幅块；
2. 进行局部收缩计算，累加贡献；
3. 立即释放接收缓存区，然后进入下一阶段。 这一设计允许用户通过调整阶段数 $K$，在通信带宽开销与峰值内存占用之间寻找最优平衡。当 $K$ 增大时，单次通信缓冲大小缩小为 $1/K$，峰值内存压力得到极大缓解，使超大型计算成为可能。

3.2 源码复现与环境配置指南

3.2.1 依赖环境安装

复现本论文的计算需要配置支持分布式并行及 GPU/CPU 混合调度的量子化学计算平台。核心依赖包括：

• Python $\ge 3.8$
• PySCF $\ge 2.3$ (量子化学引擎核心)
• mpi4py $\ge 3.1.0$ (Python 包装的 MPI 接口)
• 兼容 MPI 3.0+ 标准的并行环境 (如 OpenMPI 4.1+ 或 Intel MPI)
• MKL / OpenBLAS 高性能线性代数库

使用 Conda 配置环境的命令如下：

conda create -n peri_ccsd python=3.10
conda activate peri_ccsd
conda install -c conda-forge numpy scipy h5py mpi4py pytest
pip install pyscf

3.2.2 编译与克隆高性能计算库

克隆作者的开源仓库并将其路径加入 Python 的环境变量中：

git clone https://github.com/shuhangli98/mpi_periodic_eom_ccsd.git
cd mpi_periodic_eom_ccsd
export PYTHONPATH=$PWD:$PYTHONPATH

3.2.3 快速复现测试脚本 (以 MgO 为例)

以下脚本用于提交多节点并行任务（保存为 run_mgo.py）：

import numpy as np
from pyscf.pbc import mingo, scf, gto
# 导入作者并行包中的高层接口
from mpi_periodic_eom_ccsd import MPI_CCSD, MPI_EOM_IP, MPI_EOM_EA
from mpi4py import MPI

comm = MPI.COMM_WORLD
rank = comm.Get_rank()

# 1. 定义 MgO 晶胞参数（高精度计算需使用更完备的基组和赝势）
cell = gto.Cell()
cell.atom = '''
Mg 0.0000000000 0.0000000000 0.0000000000
O  2.1065000000 2.1065000000 2.1065000000
'''
cell.a = np.array([
    [0.0, 2.1065, 2.1065],
    [2.1065, 0.0, 2.1065],
    [2.1065, 2.1065, 0.0]
])
cell.basis = 'gth-cc-pvdz'
cell.pseudo = 'gth-hf-rev'
cell.exp_to_discard = 0.05  # 去除弥散低指数高斯基，消除线性相关
cell.build()

# 2. 布里渊区离散化 3x3x3 网格采样
kmesh = [3, 3, 3]
kpts = cell.make_kpts(kmesh, wrap_around=True)

# 3. 执行周期性 Hartree-Fock 计算
kmf = scf.KRHF(cell, kpts=kpts).density_fit()
kmf.exxdiv = 'ewald'  # 启用 Ewald 修正以降低有限尺寸效应
if rank == 0:
    print("Starting Mean-Field HF Calculation...")
kmf.kernel()

# 4. 执行分布式基态 CCSD 计算
# 控制 task splitting 数，优化峰值内存
phase_k = 4
cc_solver = MPI_CCSD(kmf, phase_k=phase_k, comm=comm)
if rank == 0:
    print("Starting Distributed CCSD Calculation...")
cc_solver.kernel()

# 5. 执行激发态 EOM-CCSD 计算带隙
if rank == 0:
    print("Starting EOM-CCSD Calculations...")
# 求解 Gamma 点处的电离电位 (IP) 和电子亲和能 (EA)
ip_solver = MPI_EOM_IP(cc_solver, nroots=1, kpt_idx=0)
ip_energy = ip_solver.kernel()

ea_solver = MPI_EOM_EA(cc_solver, nroots=1, kpt_idx=0)
ea_energy = ea_solver.kernel()

if rank == 0:
    fundamental_gap = (ip_energy + ea_energy) * 27.2114 # 转换为 eV
    print(f"=== MgO Gamma-Point Direct Gap: {fundamental_gap:.4f} eV ===")

通过 MPI 调度多节点并行执行（例如在 SLURM 超算集群中申请 4 个节点，每个节点 24 进程）：

srun -N 4 -n 96 python run_mgo.py

• 开源代码仓库: https://github.com/shuhangli98/mpi_periodic_eom_ccsd
• 基准原始数据集: https://github.com/shuhangli98/supporting_data

4. 关键引用文献与局限性批判评述

4.1 关键参考文献

1. PySCF 核心架构: Q. Sun, et al. PySCF: the Python-based simulations of chemistry framework, WIREs Comput. Mol. Sci. 8, e1340 (2018).
   • 评述：奠定了高效率高分子晶体密度拟合与轨道操作的 Python-C 混合平台。
2. 周期性波函数理论奠基: J. McClain, Q. Sun, T. C. Berkelbach, et al. Gaussian-Based Coupled-Cluster Theory for the Ground-State and Band Structure of Solids, J. Chem. Theory Comput. 13, 1209 (2017).
   • 评述：首次在晶体动量动量守恒下系统实现了周期性高斯基组的 CCSD 方程。
3. 有限尺寸渐近外推行为研究: E. Moerman, et al. Finite-Size Effects in Periodic EOM-CCSD for Ionization Energies and Electron Affinities: Convergence Rate and Extrapolation to the Thermodynamic Limit, J. Chem. Theory Comput. 21, 1865 (2025).
   • 评述：对固态体系下激发能收敛特征作出了系统的解析数学推导，是本工作 Model ABC 理论的最直接支撑。
4. GW-EOM 融合尝试: E. Moerman, A. Gallo, T. Gruber, F. Hummel, et al. Exploring the accuracy of the equation-of-motion coupled-cluster band gap of solids, Phys. Rev. B 111, L121202 (2025).
   • 评述：提出了利用 GW 计算辅助外推 EOM-CCSD 激发能的重要拟设。

4.2 局限性深度评论

尽管本工作在并行规模以及基准精度上取得了惊艳的成就，但作为科学探讨，其仍存在以下不可忽视的方法学与技术性局限：

1. 缺乏三激发（Triple Excitations）算符：精度封顶

尽管 CCSD 是相关能计算的黄金标准，但在分子量子化学中，只有包含非微扰三激发修正的 $\text{CCSD(T)}$ 才能真正被称为“金标准”（实现 $\le 0.05$ eV 的精度）。本工作也清晰表明，由于缺乏三激发，EOM-CCSD 在预测强极化的岩盐结构体系时，会系统性高估带隙约 0.4–0.9 eV。如何在本并行的框架下引入极高成本的 $\mathcal{O}(N_{atom}^7 N_k^5)$ 规模三激发计算，是决定波函数理论能否在固态中超越 GW/BSE 近似的物理关键。

2. “复合间接带隙外推法”的对称性红利不可推广性

本工作提出的复合外推法成功避开了平移网格对称性破缺对计算资源的吞噬，这要求价带顶（VBM）刚好位于 $\Gamma$ 点。对于绝大多数典型半导体（如 Si, Ge, BN），其 VBM 确实在 $\Gamma$ 点，此时该法堪称完美。但如果面临一些更为奇特的拓扑材料或过渡金属卤化物，其 VBM 本身就在非高对称性的特殊 $k$ 点，复合外推法将失效，计算仍不得不回退至“平移网格对称性破缺”的窘境。

3. 弥散基组（Diffuse Functions）导致的线性相关性梦魇

在周期性边界条件下，使用包含长程弥散基函数的 correlation-consistent 基组（如 cc-pVTZ）极易引发严重的线性相关（Linear Dependency）问题。为了稳定计算，作者不得不采用剔除重叠矩阵中特征值低于 $10^{-6}$ 的本征向量、并手动删去指数小于 0.05 的高斯弥散原子的权宜之计。这一机制虽然保证了数值计算的收敛，但人为引入了断点。微小的参数截断变动往往会引起总能量数毫电子伏特的波动，这对实现极高精度的基态热力学性质预测带来了一定的不确定性。

5. 补充讨论：有限尺寸效应与 GW 辅助外推的物理机理

5.1 有限尺寸效应（FSE）的深层物理剖析

为了更直观地理解为何带隙的有限尺寸误差收敛慢于总能量，我们需要从微观的多体自能（Self-Energy）和库仑自相互作用（Coulomb Self-interaction）来剖析。

在有限胞计算中，当我们在导带中加入一个电子（对应计算 EA）或在价带中移走一个电子（对应计算 IP）时，这个外加的电子或空穴不仅在原胞内发生极化，还会与由于周期性边界条件产生的、位于其它超胞内的**自身镜像（Images）**发生长程静电排斥与自相互作用。这一排斥能的大小正比于马德隆常数（Madelung Constant），对能量的贡献与超胞的特征线度成反比，即 $L^{-1} \sim N_k^{-1/3}$。因此，这一物理项是激发态能量缓慢收敛的核心元凶。

相比之下，对于基态总能量（计算凝聚能的核心），体系是呈电中性的，且不存在额外注入的局域电荷，电子的相关效应表现为局部“相关空穴（Correlation Hole）”的屏蔽。相关空穴的相互作用强度随距离呈偶极-偶极（$\sim R^{-6}$）或更高阶衰减，因而其有限尺寸效应通过三维空间积分后表现为与体积成反比（$\sim N_k^{-1}$）。

这两种截然不同的物理本源，决定了带隙外推时必须包含 $N_k^{-1/3}$ 及其高阶多项式，也使得在小网格下直接应用带隙外推非常危险。

5.2 GW-EOM-234 辅助外推机制及其在共价体系中的局限

为了在低至 $3^3$ 或 $4^3$ 的稀疏 $k$ 网格下预测 EOM-CCSD 的 TDL 极限，前人（Moerman 等，Ref. 19）曾提出过一种非常巧妙的 GW 辅助外推方案（GW-EOM-234）。其物理逻辑基于：由于 GW 近似与 EOM-CCSD 的有限尺寸效应都源于动力学屏蔽（Screening）和长程库仑相互作用，两者的带隙有限尺寸收敛曲线在数学结构上应该高度相似，只是在振幅上相差一个比例系数：

其中比例系数 $b$ 可以通过在极小 $k$ 网格下拟合 $E^{CC}$ 与 $E^{GW}$ 的比值轻松锁定。由于 GW 计算开销极低，可以在极其密集的网格（如 $12^3$）下进行外推获取极为可靠的 $E_{gap}^{GW}(\infty)$。一旦 $b$ 确定，我们就能用低阶 $k$ 点的 $E_{gap}^{CC}(N_k)$ 外推出极高精度的 $E_{gap}^{CC}(\infty)$，被称为 GW-EOM-234 方案。

本论文作者利用他们首次获得的 $6^3$ 稠密 CCSD 网格数据，首次对这一假说的可靠性进行了无偏的、严苛的审视（参见 Table IV 数据）：

• 在强离子体系（MgO, LiCl, LiF）中，GW-EOM 拟合表现极其出色。其预测与理论 TDL 的偏差在 0.11–0.13 eV 以内。这是因为在高度极化、电荷局部化的离子晶体中，极化筛选主要由局域偶极决定，各阶外推系数（$A/A', B/B', C/C'$）表现出了高度的各向同性比例特征（如 MgO 中，各阶比值均稳定在 1.35–1.92 之间）。
• 但在共价半导体（BN, BP, Si, C）中，该方案出现了灾难性的失败。其外推偏差飙升至 0.39–0.55 eV（严重高估）。
  • 深层物理机制缺陷：研究表明，在高度去域化的共价体系中，极化筛选具有强烈的方向性（各向异性）。这导致 GW 的二阶项、三阶项系数比例与 CCSD 严重脱节。如 Table IV 所示，对于 LiCl 而言，一阶与二阶系数比例较为平稳，而三阶系数比例高达 35.27；对于 LiF 其三阶比例甚至变号为 -2.39。这一物理本源的不匹配，使得仅通过稀疏网格拟合出一个 isotropic 比例系数 $b$ 的 GW-EOM-234 方案，在面对高各向异性的共价晶体时，会引入极其严重的外推伪影。

因此，本论文给出的一条极其宝贵的避坑指南是：在面对非离子型的共价半导体体系时，千万不要迷信 GW 辅助外推，而应老老实实地通过高性能并行计算程序，基于最少 $3^3 - 5^3$ 的无偏 CCSD 散点数据，直接采用 Model AB 或复合外推法进行纯粹的几何外推。

5.3 展望：迈向智能与百亿亿次时代的固态量子化学

本项研究标志着高精度固态波函数理论正式迈入了“大 $k$ 点、消除有限尺寸效应”的成熟期。未来的发展将聚焦于两个重要维度：

1. 异构计算（GPU）加速：周期性 CCSD 算法高度契合 GPU 的大规模密集群积（GEMM）特性。将本工作中的分布式 MPI 调度与 GPU 的按需三中心 DF 张量收缩融合，有望将计算时间进一步缩短一个数量级。
2. 机器学习辅助外推：通过引入物理约束的神经网络（PINN）来拟合倒易空间电子密度分布，替代简单的多项式外推，或许可以使我们在仅有 $2^3$ 和 $3^3$ 网格下，就能精确预测出由于电荷极化产生的全部有限尺寸修正项，彻底颠覆传统的拟合游戏规则。