DFT计算中SCF不收敛？一篇搞懂收敛机制与解决方案

SCF 收敛计算是如何实现的？

DFT 的核心是将复杂的多电子问题转化为单电子在有效势场中的运动问题（Kohn-Sham 方程）。但这里有一个"鸡生蛋，蛋生鸡"的死循环：你需要知道电子密度来计算有效势场，但你又需要知道有效势场才能求解方程得到电子密度。

在 KS-DFT 中，哈密顿算符 $H[\rho]$ 显式依赖电子密度 $\rho(\mathbf r)$，实际求解是通过 SCF 循环反复迭代 $\rho^{(n)}\to H[\rho^{(n)}]\to {\psi^{(n)}_i}\to \rho^{(n+1)}$，直到前后两步的能量或密度残差低于给定阈值。整个流程如下：

1. 初始猜测：软件首先根据原子的独立状态或预设模型，粗略"猜"一个初始的电子密度分布 $n_0(r)$。

2. 构建有效势：使用当前的电子密度，计算系统的有效势场 $V_{eff}$（包括外部势、库仑相互作用的 Hartree 势、以及交换关联势 $V_{xc}$）。

3. 求解 KS 方程：在计算出的有效势场下，求解薛定谔方程，得到单电子的波函数 $\psi_i$。

4. 计算新密度：根据求出的波函数，计算出系统"新的"电子密度 $n_{new}(r)$。

5. 收敛检查：将新的电子密度 $n_{new}(r)$（或系统总能量）与上一步的旧密度 $n_{old}(r)$ 进行对比。如果它们之间的差值小于设定的容忍度（阈值，例如 $10^{-5}$ eV），则宣布收敛，计算结束。

6. 混合与更新：如果没有达到收敛，软件不能直接把 $n_{new}(r)$ 作为下一步的输入（这会导致剧烈震荡）。大多数程序并不是直接用 $\rho^{(n+1)}{\text{out}}$ 替换 $\rho^{(n)}{\text{in}}$，而是做混合，例如 $\rho^{(n+1)}{\text{in}}=\rho^{(n)}{\text{in}}+\alpha(\rho^{(n)}{\text{out}}-\rho^{(n)}{\text{in}})$，或者用 Anderson、Pulay/DIIS、Broyden 这类拟牛顿方法利用多步历史残差构造"最佳"更新方向，以同时提高稳定性和收敛速度，然后将新旧密度按照一定比例混合，生成一个更合理的密度 $n_{in}^{(k+1)}$，然后回到第 2 步重新开始。

为什么有时候会"不收敛"？

SCF 不收敛的最核心物理现象叫做电荷震荡（Charge Sloshing）。

在迭代过程中，如果系统存在大量能量相近的空轨道和占据轨道（比如在费米能级附近），电子很容易在不同的空间区域（比如大分子的两端，或者固体表面的真空层）来回"跳跃"：

• 第 N 步：软件发现 A 端带正电，把电子移过去。
• 第 N+1 步：A 端电子太多了，排斥力骤增，软件又把电子赶回 B 端。

这种左右横跳会导致总能量和电子密度无法稳定下来，SCF 循环就会陷入死循环，直到达到最大迭代次数报错。

对金属或强非均匀体系，简单混合容易出现"charge sloshing"，即长波长电荷振荡导致 Hartree 能量在相邻步之间大幅波动，混合-纠正不断过冲，从而振荡甚至发散，因此需要配合合适的预条件器（如 Kerker 型）削弱长波分量的更新。

影响收敛的因素及解决策略

影响 SCF 收敛的因素主要分为物理体系特性和计算参数设置两大类。当遇到不收敛时，可以从以下几个方面进行排查和调整。

物理体系特性（"先天因素"）

金属体系与窄带隙半导体：费米面附近的电子态非常密集，极易发生电荷震荡。

• 对策：必须使用展宽技术（Smearing）。例如在 VASP 中调整 ISMEAR 和 SIGMA（增加 SIGMA 可以帮助收敛，但不能太大以免影响能量准确性）；Gaussian 中可以使用 Fermi 关键字。类似 DMol3 等程序的经验表明，适当的 Fermi 展宽或 smearing 有助于改善占据数对能量的敏感性，从而提升 SCF 的稳定性。

磁性体系（尤其是反铁磁/阻挫系统）：自旋向上和向下的电子密度很难找到最稳定的分布，容易卡在局部的亚稳态。

• 对策：提供高质量的初始磁矩猜测（如 VASP 中的 MAGMOM），或者先做非自旋极化计算得到波函数，再作为自旋极化计算的初始猜测。

包含强关联电子的体系（过渡金属、d/f 电子）：电子高度局域化，传统的 DFT 难以准确描述。

• 对策：使用 DFT+U 方法，并在初期设置较慢的混合速率。

大长条形分子或具有大真空层的体系（如二维材料）：偶极矩变化剧烈，长程静电相互作用容易导致跨越真空层的虚假电荷震荡。

• 对策：使用特定的混合预处理器（如 Kerker mixing），它能抑制长波长的电荷震荡。

体系物理本身"难"：金属、小带隙体系或强磁性、近简并态，都会加剧长波长密度振荡，使简单混合极不稳定。

展宽/占据处理不合适：对金属或小带隙体系，如果占据展宽（smearing）过小或者没有适当展宽，能量对占据数极其敏感，也会让 SCF 非常难收敛，多数量化软件都建议在这类体系开启合适的展宽帮助 SCF。

计算算法与参数设置（"后天调试"）

如果体系本身没问题，那就是软件参数没调好：

混合算法与权重：软件默认的混合步子可能迈得太大了。

• 对策：减小混合参数（降低新密度的权重）。在 VASP 中是减小 AMIX（例如从 0.4 降到 0.2 或更低）；在 Gaussian 中可以使用 SCF=QC（二次收敛算法，非常稳但很慢）或 SCF=XQC。混合参数过激会导致数值不稳定和振荡，这一行为在平面波 DFT 实现中普遍存在。

初始猜测太差：如果原子的初始摆放非常不合理（比如距离太近），或者直接从头猜的波函数太离谱。

• 对策：优化初始几何结构。如果是在做分子动力学或结构优化，尽量读取上一步计算好的波函数（如 VASP 的 ISTART=1, ICHARG=1；Gaussian 的 Guess=Read）。如果初始几何严重不合理或初始电荷密度离实际基态太远，SCF 很可能在局部极值附近徘徊或发散，实际经验中这在表面-吸附、缺陷、大超胞计算里尤为明显。

能级移动：有时候占据轨道和空轨道靠得太近导致震荡。

• 对策：强行在迭代过程中拉开 HOMO（最高占据轨道）和 LUMO（最低未占轨道）的能量差。Gaussian 中使用 SCF=VShift。

网格/K点精度不足：实空间网格太粗或 K 点采样不合理，导致积分误差大，引发数值不稳定。

• 对策：适当增加截断能（Cutoff energy）或加密 K 点网格。

混合参数过激：例如混合因子过大（ABACUS 中的 mixing_beta，越大表示每步密度变化幅度越大），在"好脾气"体系会加快收敛，但在困难体系反而容易导致数值不稳定和振荡。

预条件不当：预条件器如果没有有效抑制长波分量，或者参数选取得过强/过弱（如 Kerker 预条件的参数），很容易引起 charge sloshing 或收敛过慢，因此代码里一般提供若干预条件选项与参数供调节。

数值参数设置不合理：过大的能量截断、极端的 k 点设置或太"苛刻"的收敛阈值，加上有限的最大迭代步数，会使得 SCF 在给定 nstep 内达不到 toldfe/tolvrs 等阈值便被强制终止。

ABINIT 中 SCF 的实现与关键参数

ABINIT 把 SCF 迭代拆成两块：一块是密度/势的自洽（由 iscf、iprcel 等控制），另一块是每一步内部的波函数优化（由 wfoptalg、nnsclo 等控制）。在 ABINIT 里可以通过调整 iscf、iprcel、diemix/diemac/dielng、nstep、收敛阈值等一系列变量来改善收敛。

iscf：选择密度/势自洽算法

• iscf=2/3/4/5/7：在势上做简单混合、Anderson、共轭梯度或 Pulay 混合等
• iscf=12/13/15/17：在密度上做对应的混合算法
• 文档指出，范数守恒赝势默认 iscf=7（势上的 Pulay 混合），在 PAW 情形默认 iscf=17（密度上的 Pulay），这是速度与稳定性之间的折中；iscf=2 这类简单混合则"更安全但更慢"

iprcel 及 diemix/diemac/dielng：SCF 预条件器

• iprcel 设置势/密度残差的预条件类型，iprcel=0 为默认，一般就能工作
• 对强非均匀体系，文档推荐尝试 iprcel=45 这类更高级预条件，并通过 diemix/diemac/dielng 定义广义 Kerker 型预条件参数，可以显著加快收敛

收敛阈值：toldfe/tolvrs/tolwfr

• toldfe：能量差收敛标准
• tolvrs：势或密度残差的范数阈值
• tolwfr：波函数残差阈值
• 这些阈值与 nstep 一起决定 SCF 何时停止，论坛中不少"SCF 没到阈值但提前停止"的报错就是 nstep 不够配合过严的阈值导致

nstep

单次 GS 计算中允许的最大 SCF 循环次数，若在 nstep 内达不到 toldfe/tolvrs 等收敛条件，就报"不收敛"并停止。

波函数优化相关：nnsclo、wfoptalg、nbdbuf

• nnsclo 控制每个 SCF 步内部为求解本征值/波函数所需的迭代收敛性，与 iscf 配合决定整体效率和稳定性
• wfoptalg 则选择具体波函数优化算法（不同共轭梯度/直接最小化等），与上面密度/势混合部分相互独立
• 在一些扩展分子动力学（EXTFPMD）场景中，有开发者经验表明适当增加"带缓冲" nbdbuf 可以明显改善潜在残差的 SCF 收敛，避免在给定步数内收敛停滞

ABINIT 中改善 SCF 收敛的常用"处方"

结合上面的机制，可以在 Abinit 里按以下思路系统地排查和调参：

先选一个稳健的 SCF 算法

• 对一般体系，保持默认 iscf=7（范数守恒）或 iscf=17（PAW）即可，这对应势/密度上的 Pulay 混合，一般兼顾稳健和效率
• 若遇到剧烈振荡，可以临时改用更"保守"的简单混合如 iscf=2 或密度简单混合 iscf=12 让 SCF 稍微"慢一点，但更稳"

调小混合强度、增强预条件

• 减小与混合相关的参数（如 diemix 等效于降低步长）能直接压制振荡
• Kerker 等预条件通过在 G 空间对残差加权来抑制长波长分量，可以显著改善金属体系的 SCF 收敛，这是 CP2K 与 ABACUS 等实现中重点强调的技巧
• 对强非均匀体系或表面/缺陷超胞，尝试把 iprcel 设为 45 这类更复杂的广义 Kerker 预条件，文档明确指出这类选择对非均匀体系收敛有显著帮助

合理设置收敛阈值与 nstep

• 如果你要求非常严的 toldfe=1e-11 或 tolvrs<1e-12，务必相应增大 nstep（比如 100–200），否则论坛中那种"nstep 不够、最大能量差刚好略高于阈值就停"的情形就会出现
• 作为调试，可以先用稍松的阈值让 SCF 找到一个"还不错"的参考态，再在此基础上收紧阈值进行高精度计算

对金属/小带隙体系使用展宽

• 类似 DMol3 等程序的经验表明，适当的 Fermi 展宽或 smearing 有助于改善占据数对能量的敏感性，从而提升 SCF 的稳定性
• 在 Abinit 中，你可以通过相应的占据控制和电子温度参数（如选择合适的 occopt 方案并设置电子温度）来起到类似效果

检查模型和波函数收敛

• 如果是缺陷或大超胞，先在较小体系或较松参数下进行预收敛，再迁移电荷密度和波函数到目标计算，可显著减少 SCF 步数与发散概率，这一点在 Abinit 论坛关于缺陷 SCF 不收敛的讨论中多次被提及
• 对于 DFPT/EXTFPMD 等波函数更敏感的场景，可以适当增大 nbdbuf（带缓冲）和相应的波函数迭代参数，社区测试显示这能明显改善势残差的收敛情况

总结建议

遇到 SCF 不收敛，不要盲目增加最大迭代次数（NELM 或 MaxCycle）。如果前 30 步能量还在剧烈上下跳动，算 1000 步通常也是徒劳的。

最常规的"救火"顺序是：

1. 检查结构是否合理
2. 增加 Smearing 展宽
3. 减小电荷混合参数（AMIX 或 diemix）
4. 尝试更高级的收敛算法（如 QC）
5. 调整预条件参数（如 iprcel=45）
6. 合理设置收敛阈值与步数

参考文献