Abinit SCF 波函数收敛过程：原理 + 代码实现详解
基于 Abinit 10.4.7 源码分析 核心模块：
src/94_scfcv/m_scfcv_core.F90、src/79_seqpar_mpi/m_vtorho.F90、src/79_seqpar_mpi/m_vtowfk.F90、各对角化器m_cgwf.F90/m_lobpcg2.F90/m_chebfi.F90/m_rmm_diis.F90、src/56_mixing/m_ab7_mixing.F90
引言：SCF 波函数收敛是什么？
DFT-Kohn-Sham 方程
$$\Bigl[ -\tfrac{1}{2}\nabla^2 + V_{\text{loc}}(\mathbf{r}) + V_{\text{nl}} + V_{\text{H}}[n] + V_{\text{xc}}[n] \Bigr],\psi_{n\mathbf{k}} = \epsilon_{n\mathbf{k}},\psi_{n\mathbf{k}}$$
是非线性特征值方程：哈密顿量本身依赖于通过 $\psi$ 构造出的电子密度 $n(\mathbf{r}) = \sum_{n\mathbf{k}} f_{n\mathbf{k}}|\psi_{n\mathbf{k}}|^2$。Abinit 用自洽场（SCF）迭代来求解：
固定 V → 解 H ψ = ε ψ → 由 ψ 算新密度 n → 算新势 V → ... 直到收敛每次迭代里要做两件相互嵌套的事：
- 外循环（SCF 大循环）：在当前 $V$ 下重新求一组 ${\psi_{n\mathbf{k}}}$，更新密度，再用混合算法把新密度 / 新势"安全地"代入下一轮。
- 内循环（波函数优化）：在每个 k 点上对固定的 $H_\mathbf{k}$，迭代极小化 Rayleigh 商 $\langle\psi|H|\psi\rangle/\langle\psi|\psi\rangle$，得到最低的若干本征对。
外循环的稳定性由"密度/势混合"决定；内循环的速度由"对角化器"决定。两者共同决定 SCF 收敛性。
输入参数总览
下表列出影响波函数收敛过程的所有关键变量（位置 abimkdocs/variables_abinit.py）。
| 输入变量 | 默认值 | 控制对象 | 含义 |
|---|---|---|---|
nstep | 30 | 外循环 | SCF 最大迭代次数 |
iscf | 7（NC）/ 17（PAW） | 外循环 | 混合方案（2 simple, 3/4 Anderson, 7 Pulay；+10 表示密度混合） |
tolvrs | 0 | 外循环 | 势残差平方收敛阈值 |
toldfe | 0 | 外循环 | 总能差收敛阈值 |
toldff | 0 | 外循环 | 最大原子力变化阈值 |
tolrff | 0 | 外循环 | 力的相对变化阈值 |
tolwfr | 0 | 外循环 | 波函数残差阈值（非自洽强制选） |
diemix | 1.0（NC）/ 0.7（PAW） | 外循环 | 混合因子 α |
diemac | 1e6 | 外循环 | 模型介电函数（介质常数） |
dielng | 1.08 | 外循环 | 介电模型的长波截断 |
npulayit | 7 | 外循环 | Pulay 混合保留的历史步数 |
iprcel | 0 | 外循环 | 介电预条件器类型 |
wfoptalg | 0（NC）/ 10（PAW）/ 114（kgb） | 内循环 | 对角化算法选择 |
nline | 4 | 内循环 | 单条带的最大 line minimization 次数 |
nnsclo | 0（auto） | 内循环 | 每个 SCF 步内的非自洽循环次数 |
nbdblock | 1 | 内循环 | LOBPCG / CG 的能带块大小 |
mdeg_filter | 6 | 内循环 | Chebyshev 滤波多项式最大阶 |
rmm_diis | 0 | 内循环 | 启用 RMM-DIIS（>0 在第 N 步后切换） |
tolrde | 0.005 | 内循环 | 单步 CG 改善阈值（早退条件） |
ortalg | -2 | 内循环 | 正交化算法选择 |
nbdbuf | 0 | 内循环 | 波函数收敛检查的"缓冲带"（顶部不算） |
mffmem | 1 | 内循环 | 是否把波函数缓存到磁盘 |
densfor_pred | 6 | 几何步 | 几何步间的密度外推方案 |
wfoptalg 取值（mod(wfoptalg, 10) 决定主算法）：
| 值 | 算法 | 适用 |
|---|---|---|
| 0 | 标准 CG（Payne 1992） | NC 默认 |
| 10 | CG + Kresse 改进（修正动能 + 改进预处理） | PAW 默认 |
| 1 / 111 | Chebyshev 滤波（ChebFi，Levitt 2015） | 大规模并行 |
| 4 / 14 / 114 | LOBPCG（Knyazev 2001） | 中等并行 / paral_kgb=1 默认 |
| 2 / 3 | Residual minimization w.r.t. shift | 非自洽 |
iscf 取值（mod(iscf,10) 决定混合算法，iscf>=10 表示混合密度而非势）：
| 值 | 算法 |
|---|---|
| 2 / 12 | 简单线性混合 |
| 3 / 13 | Anderson（用 1 步历史） |
| 4 / 14 | Anderson 2 阶（用 2 步历史） |
| 5 / 6 / 15 / 16 | CG 能量极小化 |
| 7 / 17 | Pulay（用 npulayit 步历史） |
整体调用链
m_gstate (顶层 SCF 驱动)
└─ scfcv (m_scfcv.F90)
└─ scfcv_core (m_scfcv_core.F90:263) ← SCF 外循环主体
├─ ab7_mixing_new (m_ab7_mixing.F90:161) ← 初始化混合器
│
└─ do istep = 1, nstep ← SCF 大循环
├─ vtorho (m_vtorho.F90) ← 由 V 求新 ρ
│ └─ do isppol, do ikpt ← 自旋 × k 点循环
│ └─ vtowfk (m_vtowfk.F90:177) ← 单 k 点波函数优化
│ ├─ cgwf (CG)
│ ├─ lobpcgwf2 (LOBPCG)
│ ├─ chebfi (Chebyshev 滤波)
│ ├─ rmm_diis (RMM-DIIS)
│ └─ subdiago (子空间对角化)
│
├─ newocc (m_occ.F90) ← 更新占据数 + Fermi 能级
│ (见 smearing 报告)
├─ rhotoxc → vhartre → setvtr ← 由 ρ 算新 V
├─ etotfor (m_scfcv_core.F90:2561) ← 算总能、力、应力
├─ scprqt (m_common.F90:183) ← 检查收敛
├─ if (quit) exit
└─ ab7_mixing_eval (m_ab7_mixing.F90:623)← 用历史信息生成下一轮 V/ρ外循环：SCF 主体（scfcv_core）
代码位置：src/94_scfcv/m_scfcv_core.F90:1034-2300。
SCF 主循环骨架
do istep = 1, max(1, nstep)
! ----- 1. 几何相关的一次性准备（仅 istep=1 或原子移动后）-----
if (moved_atm_inside == 1 .or. istep == 1) then
call symmetrize_xred(...)
call getcut(...) ! G 球截断
call getph(...) ! 结构因子
! PAW 时还要做 nhatgrid 等
end if
! ----- 2. 由当前 V 求新 ρ：vtorho 内含 k 点循环 -----
call vtorho(... eigen, occ, rhog, rhor, residm, ...)
! 返回最大波函数残差 residm
if (dtset%iscf < 0) exit ! 非自洽：vtorho 内已经迭代完，单步退出
! ----- 3. 在密度混合情况下，先算总能 -----
if (iscf >= 10) then
call etotfor(... etotal, deltae, ...)
end if
! ----- 4. 检查收敛 -----
call scprqt(... tollist, ... quit, conv_retcode)
if (quit == 1) exit
! ----- 5. 用历史信息混合（生成下一步的 V 或 ρ）-----
call ab7_mixing_eval(mix, ... istep, ...)
! ----- 6. 在势混合情况下，用混合后的 ρ 重算 V，再算总能 -----
if (iscf < 10) then
call rhotoxc(... rhor → vxc)
call hartre (rhog → vhartr)
call etotfor(... etotal, deltae, ...)
call scprqt(... quit ...) ! 再检查一次
end if
end do收敛判据：scprqt
位置：src/67_common/m_common.F90:183。从 tollist(1:8) 数组中读 5 个判据（详见之前的 SCF convergence 报告）：
tolwfr = tollist(2) ! 最大波函数残差
toldff = tollist(3) ! 最大力变化
toldfe = tollist(4) ! 总能差
tolvrs = tollist(6) ! 势/密度残差²
tolrff = tollist(7) ! 相对力变化只允许选一个（5 个里只有 1 个非零）。监控变量：
residm：所有 k 点和带的最大波函数残差 $\max_{n\mathbf{k}}|(H-\epsilon)\psi|^2$deltae = E_total(istep) - E_total(istep-1)res2：势/密度残差范数 $|V_{\text{new}}-V_{\text{old}}|^2$diffor：原子力最大变化
当对应判据连续满足 2 次（避免假阳性）时 quit = 1，跳出循环。
内循环：波函数优化（vtowfk）
代码位置：src/79_seqpar_mpi/m_vtowfk.F90:177-700。
分发器骨架
subroutine vtowfk(...)
wfoptalg = mod(dtset%wfoptalg, 100)
wfopta10 = mod(wfoptalg, 10)
...
! 把波函数过滤：动能 > huge 的 G 分量置零（防止"高频污染"）
do iband, do ipw
if (kinpw(ipw) > huge*1d-11) cwavef(:, ipw) = 0
end do
! 非自洽循环：每次 SCF 大步内可能多次调对角化器
do inonsc = 1, nnsclo_now
! 根据 wfopta10 派发到不同对角化器
if (wfopta10 == 4) then
call lobpcgwf2(...)
else if (wfopta10 == 1) then
call chebfi/chebfiwf2(...)
else ! wfopta10 = 0 或 10
if (.not. use_rmm_diis) then
call cgwf(...) ! 标准 CG
else
call rmm_diis(...) ! 切换到 RMM-DIIS
end if
end if
! 子空间对角化（Rayleigh-Ritz）：在已找到的 nband 维子空间中做精确对角化
if (do_subdiago) call subdiago(...)
! 正交化（部分算法跳过）
if (do_ortho) call zorthonormalize(...)
end do
end subroutinennsclo_now（非自洽循环次数）由 nnsclo 决定，自洽计算通常 = 1（每 SCF 步只调一次对角化器）；非自洽计算 = nstep（在固定 V 下反复优化波函数直到收敛）。
标准共轭梯度（CG，cgwf）
代码位置：src/67_common/m_cgwf.F90:488-1300。这是 Abinit 历史最悠久的算法，Teter–Allan–Payne 1992 的实现。
算法概要
对每个能带 $n$，逐一地求解：
$$\min_{|\psi_n\rangle;\perp;|\psi_1\rangle,\ldots,|\psi_{n-1}\rangle,;|\psi_n|=1} \langle\psi_n|H|\psi_n\rangle$$
用预条件 CG 在单位球面 + 正交补空间上做线搜索。
关键步骤（cgwf 第 488-1000 行）
对每个 iband = 1..nband：
- 归一化 $|\psi_n\rangle$
- 计算梯度 $|H\psi_n\rangle$：
call getghc(...)，这一步是整个 DFT 的 Hot Spot（每条 band 一次 FFT × 2 + 非局域投影） - 进入 line minimization 循环（最多
nline次）：
do iline = 1, nline
! a) Rayleigh quotient
chc = <ψ|H|ψ> ! 当前能量
lam0 = chc
! b) 残差向量与残差范数
|vresid> = |Hψ> - chc * |ψ>
resid(iband) = ‖vresid‖²
! c) 检查收敛 (tolwfr 是单 band 残差阈值)
if (resid(iband) < tolwfr) exit
! d) 投影到正交子空间：投掉前面 n-1 个 band
|direc> = |vresid> - Σ_{m<n} |ψ_m><ψ_m|vresid>
! e) Teter 多项式预处理（关键！）
call cg_precon(ψ, 0, ..., kinpw, ..., pcon, direc, ...)
! 预条件矩阵 K：在 G 空间对角，K(G) = (27+18x+12x²+8x³)/(27+18x+12x²+8x³+16x⁴)
! 其中 x = (k+G)²/⟨ψ|T|ψ⟩，⟨ψ|T|ψ⟩ 是当前能带的动能平均
|direc> ← K |direc>
! f) Polak-Ribière / Fletcher-Reeves: 加入历史方向
dotgg = <direc|H|ψ>
if (iline == 1) then
γ = 0
|conjgr> = |direc>
else
γ = dotgg / dotgp ! Fletcher-Reeves 系数
|conjgr> = |direc> + γ |conjgr>
end if
dotgp = dotgg
! g) 把 conjgr 投影到 ψ 的正交补 → direc
<ψ|conjgr> 投掉
! h) 在 span{ψ, direc} 平面内做精确二维 Rayleigh quotient 极小化
compute <ψ|H|ψ>, <direc|H|direc>, <ψ|H|direc>
tan(2θ) = 2 <ψ|H|direc> / (<ψ|H|ψ> - <direc|H|direc>)
ψ_new = cos(θ) * ψ + sin(θ) * direc
end do关键点
- Teter 多项式预处理：Abinit 招牌技术。预条件器 $K(G)$ 在小 $|G|$ 处近似 1（不改动低频信号），在大 $|G|$ 处近似 $1/T(G)$（压制高频）。代码
cg_precon在src/44_abitools/m_cgtools.F90:3651-3740。这一步把 CG 收敛速度从 $O(\kappa)$ 提到 $O(\sqrt{\kappa})$ 级别，是平面波 DFT 的标配。 - 解析的二维 line minimization：因为 Rayleigh quotient 在 2 维平面内是有理函数，可解析求其最优 θ（avoid bracket-and-search），效率最高。
- 单条带串行：CG 不能在 band 维度并行（每条带依赖前面的正交化结果）。这是其唯一劣势。
LOBPCG（Locally Optimal Block PCG）
代码位置：src/48_diago/m_lobpcg2.F90:340-700 和 src/79_seqpar_mpi/m_lobpcgwf.F90。
算法概要
把 $nband$ 条带分成大小为 blockdim 的若干块，每块一起优化（block iteration）：
对一块 $X \in \mathbb{C}^{N_{\text{pw}}\times b}$：
- 计算 $W = K\cdot R$，其中 $R = AX - BX\Lambda$（残差，$\Lambda$ 是 Rayleigh-Ritz 给出的本征值矩阵）
- 构造增广子空间 $V = [X, W, P]$（$P$ 是上一步的"惯性方向"）
- 在 $V$ 上做 Rayleigh-Ritz：求 $V^\dagger A V$ 与 $V^\dagger B V$ 的广义本征对，取最低 $b$ 个
- 更新 $X_{\text{new}}$，$P_{\text{new}} = $ 新 $X$ 在 $W,P$ 子空间的分量
代码骨架（m_lobpcg2.F90:456-700）
do iblock = 1, nblock
call lobpcg_getX0(...) ! 取当前块初始波函数
call lobpcg_orthoXwrtBlocks(...) ! 与之前所有 iblock 块 B-正交化
call getAX_BX(X, AX, BX) ! AX = H·X, BX = S·X
call xg_Borthonormalize(X, BX, AX) ! B-正交化 X
call xg_RayleighRitz(X, AX, BX, λ) ! 初始 RR，给出本征值估计 λ
do iline = 1, nline
! 1. 残差
W = AX - BX·diag(λ)
maxResidu = max ‖W_i‖
if (maxResidu < tolerance) exit
! 2. 预条件
W = K·W
! 3. 与之前块正交化
W = W - X_prev (X_prev^† B W)
! 4. 算 AW, BW
call getAX_BX(W, AW, BW)
! 5. Rayleigh-Ritz on [X, W, P] (or [X, W] in first iter)
if (iline == 1) then
call xg_Borthonormalize(XW, BXW, AXW)
call xg_RayleighRitz(X, AX, BX, λ, XW=XW, AW=AW, BW=BW, P, AP, BP, ...)
else
call xg_Borthonormalize(XWP, BXWP, AXWP)
call xg_RayleighRitz(X, AX, BX, λ, XW, AW, BW, P, AP, BP, XWP, ...)
end if
end do
end do关键点
- 块运算：所有矩阵-向量乘法变成矩阵-矩阵乘法（BLAS-3），缓存友好，能用 OpenMP/SIMD 加速。
- MPI 友好：不同块可以分到不同进程（
paral_kgb=1下），bandfft 转置（xgTransposer）使本征值问题在小子空间内做。 - 数值技巧：当残差降到 ~ 1e-27 时，$P$ 矩阵会出现 catastrophic cancellation；代码用
if minResidu < 1e-27重新做 $[X,W]$ 子空间的 RR 而不是 $[X,W,P]$。 - 推荐用法：
wfoptalg 114，配paral_kgb 1使用 band+FFT 并行。
Chebyshev 滤波（ChebFi）
代码位置：src/66_wfs/m_chebfi.F90:300-450。Levitt-Torrent 2015。
算法概要
不直接迭代单个本征对，而是用 Chebyshev 多项式 $T_m(H)$ 作为滤波器：对感兴趣谱区放大，其余压制。
设感兴趣区间 $[\lambda_{\min}, \mu]$（μ ≈ 当前最高占据本征值），让滤波器在区间外 $[\mu, \lambda_{\max}]$ 上小、区间内大：
$$\psi_{\text{new}} = T_m!\Bigl(\frac{2H - (\mu+\lambda_{\max})}{\lambda_{\max} - \mu}\Bigr);\psi_{\text{old}}$$
代码核心（m_chebfi.F90:340-440）
filter_center = (ecut + filter_low) / 2 ! μ 中心
filter_radius = (ecut - filter_low) / 2 ! μ 半径
do ideg = 1, mdeg_filter
! 应用 H
call getghc(cg, ghc, ...) ! H·ψ
! PAW 时还要应用 S^{-1}
if (paw) call apply_invovl(ghc, gsm1hc, ...)
else gsm1hc = ghc
! Chebyshev 递推：
! T_0(x) = 1, T_1(x) = x, T_{n+1}(x) = 2x T_n(x) - T_{n-1}(x)
if (ideg == 1) then
cg_next = 1/filter_radius * (gsm1hc - filter_center * cg)
else
cg_next = 2/filter_radius * (gsm1hc - filter_center * cg) - cg_prev
end if
cg_prev = cg ; cg = cg_next
end do
! 滤波完后做 Rayleigh-Ritz：在 nband 维子空间精确对角化
call subdiago(...)关键点
- 不做正交化也不做线搜索！整套 nband 条带"被同一个多项式作用"，对角化只在最后一步做。
- 多项式阶
mdeg_filter（默认 6）控制滤波质量。阶越高，滤波越锐利，但每步只多做几次 FFT，开销可控。 - 不可使用预条件：因为多项式作用是显式的。这导致它对
ecutsm不敏感，Pulay 应力无法修正（必须靠加大ecut）。 - 适合超大规模：每步 nband × mdeg_filter 次 FFT，可完美并行；当 LOBPCG 内部 RR 成为瓶颈时（>1000 bands），ChebFi 更优。
- 顶部几条 band 难收敛（多项式在尾部弱），代码强烈建议配
nbdbuf多算几条 band。
RMM-DIIS
代码位置：src/66_wfs/m_rmm_diis.F90，由 m_vtowfk.F90:465 调用。
Residual Minimization with Direct Inversion in Iterative Subspace（来自 VASP 工作）。在 CG 稳定后期切换可大幅加速。当 rmm_diis > 0 且 istep > 3 + rmm_diis 时启用：
use_rmm_diis = istep > 3 + dtset%rmm_diisDIIS 思想：保留若干个 (波函数, 残差) 历史对，求一组系数 $c_i$ 使线性组合 $\sum c_i R_i$ 范数最小，然后用同样系数组合波函数得到新解。
关键工具：getghc 与预条件器
getghc：$|H\psi\rangle$ 的核心算子
所有对角化器最贵的一步都是计算 $H|\psi\rangle$。Abinit 的 getghc 完成：
H |ψ> = T|ψ> <- G 空间对角，O(N_pw)
+ V_loc(r) ψ(r) <- 需要 FFT： r ←→ G
+ V_nl |ψ> <- 投影到原子基（nonlop），O(N_pw·N_proj)
+ V_Fock |ψ> <- 仅杂化泛函；最贵 O(N_pw² ln)时间分配（典型）：FFT ≈ 40%, nonlop ≈ 40%, 其他 20%。
Teter 预条件器
代码：src/44_abitools/m_cgtools.F90:3651-3740。在 G 空间对角，每个 G 用一个标量：
! 计算当前条带的平均动能
ek0 = Σ_G kinpw(G) |ψ(G)|²
! 对每个 G：
x = kinpw(G) / ek0
poly = 27 + x(18 + x(12 + x·8)) ! Teter 三阶多项式
fac = poly / (poly + 16·x⁴) ! 关键比率
pcon(G) = fac
! 应用到方向向量
vect(G) = (vect(G) - eval·ψ(G)) · fac数学上 $K(x)$ 在 $x\to 0$ 时 $\to 1$（保留低频），在 $x\to\infty$ 时 $\to (27+18x+...) / 16x^4 \to 0$（压制高频）。这把 Hessian 的条件数从 $E_{\text{cut}}/E_{\text{gap}}$（典型 1000）压到 $O(1)$。
密度/势混合算法
代码位置：src/56_mixing/m_ab7_mixing.F90。
混合的物理意义
设 SCF 输入势 $V_{\text{in}}^{(i)}$，由 KS 方程算得输出势 $V_{\text{out}}^{(i)} = F[V_{\text{in}}^{(i)}]$。SCF 不动点条件是 $V_{\text{out}} = V_{\text{in}}$。残差是 $R^{(i)} = V_{\text{out}}^{(i)} - V_{\text{in}}^{(i)}$。
朴素地令 $V_{\text{in}}^{(i+1)} = V_{\text{out}}^{(i)}$（不混合）会振荡发散（因为电子云对势变化反应过激）。需要"减速"：
$$V_{\text{in}}^{(i+1)} = V_{\text{in}}^{(i)} + \alpha,K,R^{(i)}$$
其中 $\alpha = $ diemix、$K = $ 模型介电预条件器。
七种混合算法（按 iscf 取值）
简单线性混合（iscf=2/12）
V_in[i+1] = V_in[i] + diemix * K_prec * R[i]代码（scfopt 第 1956-1960 行）：
f_fftgr(:, :, i_vstore) = vtrial(:, :) ! 存当前 V_in
vtrial(:, :) = vtrial(:, :) + f_fftgr(:, :, i_vrespc(1)) ! V_in + 预条件残差最稳定但收敛最慢。diemix 通常 0.3-1.0。
Anderson 1 阶（iscf=3/13）
代码（scfopt 第 1964-2018 行）。用上一步信息构造最优线性组合：
$$V_{\text{in}}^{(i+1)} = (1-\lambda)(V_{\text{in}}^{(i)} + R^{(i)}) + \lambda(V_{\text{in}}^{(i-1)} + R^{(i-1)})$$
其中
$$\lambda = \frac{\langle R^{(i)}, R^{(i)} - R^{(i-1)}\rangle}{|R^{(i)} - R^{(i-1)}|^2}$$
是直线段上极小化残差范数的最优系数。代码：
lambda = (resid_new - prod_resid) / (resid_new + resid_old - 2*prod_resid)Anderson 2 阶（iscf=4/14）
用 2 步历史，2×2 线性方程组。代码 m_ab7_mixing.F90:2021-2089：
det = aa1*aa2 - bb*bb
lambda = (aa2*cc1 - bb*cc2) / det
lambda2 = (aa1*cc2 - bb*cc1) / detPulay 混合（iscf=7/17，Abinit 默认）
代码 m_ab7_mixing.F90:2092-2200。用 npulayit（默认 7）步历史，求 $n$ 维线性方程组：
$$\min_{{\alpha_i}} \Bigl|\sum_{i=1}^n \alpha_i R^{(i)}\Bigr|^2,\qquad \sum \alpha_i = 1$$
构造 Pulay 矩阵 $A_{ij} = \langle R^{(i)}, R^{(j)}\rangle$，解约束线性方程组：
! 构造 A 矩阵
do ii = 1, niter
call dotprodm_v(... amat(ii, niter) ...) ! A_{i,n} = <R_i, R_n>
end do
! 求逆
call dgetrf(niter, niter, amatinv, ...)
call dgetri(niter, amatinv, ...)
! 取 alpha_i = A^{-1}[:, 1] 行和归一
alpha(:) = sum(amatinv, dim=2) / sum(amatinv)
! 新的混合势：
V_in[i+1] = Σ alpha_i (V_in[i] + R[i])Pulay 与 Anderson(npulayit) 等价：npulayit = 1 即 Anderson 1 阶，= 2 即 Anderson 2 阶。
CG 能量极小化（iscf=5/6/15/16）
把 $V$（或 $\rho$）当成无约束变量，用 CG 极小化 $E_{\text{tot}}[V]$。需要计算 $\partial E/\partial V$，开销大但理论上单调下降。仅 DEVELOP 用。
介电模型预条件 $K_{\text{prec}}$
在做混合前，先把残差 $R$ 用模型介电函数 $K = \varepsilon_{\text{model}}^{-1}$ 处理一遍（iprcel 控制类型）：
$$K(\mathbf{G}) = \frac{1}{\varepsilon_{\text{model}}(|\mathbf{G}|)}$$
最简单（iprcel=0）的 Thomas-Fermi 模型：
$$\varepsilon(G) = 1 + (\varepsilon_\infty - 1),\frac{G_{\text{TF}}^2}{G^2 + G_{\text{TF}}^2}$$
由 diemac（$\varepsilon_\infty$）和 dielng（Thomas-Fermi 长度的倒数）控制。对金属：diemac = 1e6；对半导体：diemac ≈ 10；对分子（真空）：diemac = 1。
SCF 一次完整迭代的逐步追踪
设输入：FCC Si, 4×4×4 k 点 → 10 IBZ 点, 8 bands, ecut=15 Ha, occopt=1, iscf=7, tolvrs=1e-12, wfoptalg=0 (CG), nstep=30。
Step 0: 初始化（istep=0 前）
- 读 WFK 或随机初始化 $\psi^{(0)}_{n\mathbf{k}}$
- 由 $\psi^{(0)}$ 算初始密度 $\rho^{(0)}$
- 算初始势 $V^{(0)}{\text{loc}} = V{\text{ion}} + V_H[\rho^{(0)}] + V_{xc}[\rho^{(0)}]$
ab7_mixing_new(mix, iscf=7, ...)分配 Pulay 历史空间
Step 1（istep=1）
进入 vtorho
- 对自旋通道
isppol=1（非极化） - 对每个 k 点
ikpt = 1..10（本进程负责的）：
在 k=1 处调 vtowfk
wfopta10 = 0→ 走 CG 分支cgwf
cgwf 内对每条 band（1..8）：
-
band 1:
- 归一化 $\psi$
getghc→ $|H\psi\rangle$（一次 FFT 正反 + nonlop）- 计算 $\langle H\rangle = \langle\psi|H|\psi\rangle$
- 残差 $|R\rangle = |H\psi\rangle - \langle H\rangle|\psi\rangle$，$|R|^2 = ?$
- 若 $|R|^2 < $
tolwfr_diago（默认 1e-22）→ skip - 否则 line minimization 循环（最多
nline=4次）：- 投影：$|R\rangle \to |R\rangle - \sum_{m<1} |\psi_m\rangle\langle\psi_m|R\rangle$（band 1 时无前导）
- Teter 预处理：$|D\rangle = K\cdot|R\rangle$
- CG 合并方向：$|C\rangle = |D\rangle + \gamma|C_{\text{prev}}\rangle$
- 投影 $|C\rangle \perp |\psi\rangle$
getghc算 $|HC\rangle$- 解析求 $\theta$：$\tan(2\theta) = 2\langle\psi|H|C\rangle/(\langle\psi|H|\psi\rangle - \langle C|H|C\rangle)$
- $|\psi_{\text{new}}\rangle = \cos\theta|\psi\rangle + \sin\theta|C\rangle$
-
band 2: 相同流程，但每步都先把 band 1 投掉
-
...
-
band 8
子空间对角化（Rayleigh-Ritz）
CG 收敛到正交基附近，但不严格在最优旋转上。subdiago 在 8 维子空间精确对角化：
$$H_{nm} = \langle\psi_n|H|\psi_m\rangle,\quad \text{对角化};H_{nm} \to \epsilon_n,;\psi_n^{\text{new}} = \sum_m c_{nm}\psi_m$$
累加 ρ 与能量
do iband = 1, 8
ρ(r) += wtk(1) · occ(iband) · |ψ_{iband,k=1}(r)|²
E_band += wtk(1) · occ(iband) · ε(iband)
end do循环回 k=2, 3, ..., 10 重复 1.2-1.5
vtorho 结束，得到 residm
residm = max over all (k, n) of resid_k(n)算总能、新势
call newocc(... → fermie, occ) ! 占据数 + Fermi 能（见 smearing 报告）
call rhotoxc(rho → vxc)
call hartre(rho → vhartr)
V_out = vpsp + vhartr + vxc
R = V_out - V_in
res2 = ‖R‖²
deltae = E_total[istep=1] - E_total[istep=0]检查收敛（scprqt）
if (tolvrs > 0 .and. res2 < tolvrs) then
tolvrs_ok = tolvrs_ok + 1
if (tolvrs_ok >= 2) quit = 1 ! 连续 2 次满足才退出
end ifistep=1 通常远未收敛，res2 ≈ 1e-2，继续。
Pulay 混合（ab7_mixing_eval）
istep=1 时还没有历史，退化为简单混合：
V_in[2] = V_in[1] + diemix * K_prec * R[1]存：f_fftgr(:, :, i_vrespc(1)) = R[1], f_fftgr(:, :, i_vtrial(1)) = V_in[1]。
Step 2（istep=2）
同样调 vtorho，但已经有 V_in[2]
- 用上一步收敛的 $\psi^{(1)}_{n\mathbf{k}}$ 作初值（很好的猜测）
- 每条 band 只需 1-2 次 line minimization 即可收敛到当前 V
Pulay 混合：1 步历史 = Anderson 1 阶
λ = <R[2], R[2]-R[1]> / ‖R[2]-R[1]‖²
V_in[3] = (1-λ)(V_in[2]+R[2]) + λ(V_in[1]+R[1])Step 3..7（istep=3-7）
每步 Pulay 历史增加，第 7 步用 7 步历史的 Pulay。res2 通常每步缩小 1-2 个数量级。
istep res2 deltae residm
1 1e-2 — 1e-3
2 3e-4 -0.012 2e-5
3 8e-7 -1.2e-4 3e-7
...
7 5e-13 -1e-10 2e-12Step 8（收敛退出）
res2 = 5e-13 < tolvrs = 1e-12 第二次满足 → quit = 1 → 退出 SCF 循环。
Step 后处理
- 写 WFK：${\psi_{n\mathbf{k}}}$ 完整谱
- 写 EIG：${\epsilon_{n\mathbf{k}}}$
- 算力、应力（如需）
- 报告
conv_retcode = 0
典型坑与诊断
SCF 不收敛振荡
- 金属体系 + diemix 过大 → 调小
diemix到 0.3-0.5 diemac错误：金属用 1e6，绝缘体用 6-10。给绝缘体设 1e6 会过预处理。- k 点不够密：金属需要 ≥ 16³ 配 smearing。
SCF 单调下降但极慢
- 切到 Pulay (
iscf 7/17) + 增大npulayit到 10-15 - 启用
iprcel 45用更好的预条件器 - 几何步间用
densfor_pred 6做密度外推
本征值求解器选择
| 场景 | 推荐 wfoptalg | 理由 |
|---|---|---|
| 小系统 (<100 原子)，单进程 | 0（CG） | 最稳定，开销低 |
| 中等系统，MPI ≤ 32 | 114（LOBPCG） | 块运算高效 |
| 大系统 (>500 原子)，MPI > 64 | 111（ChebFi） | RR 不再瓶颈 |
| PAW 计算 | 10（CG）/ 114 | 改进预条件器 |
| 难收敛/磁性 | 0 + nline 6-8 | CG 鲁棒性最好 |
tolwfr_diago vs tolwfr
tolwfr_diago：单 k 点内对角化器的目标残差。默认与 nstep 有关。tolwfr：作为 SCF 收敛判据（仅非自洽必填）。- 两者别混淆：自洽计算中 SCF 用
tolvrs，内部用tolwfr_diago控制对角化精度。
关键代码地图
| 模块 / 子程序 | 文件 | 行号 | 作用 |
|---|---|---|---|
scfcv_core | src/94_scfcv/m_scfcv_core.F90 | 263 | SCF 外循环主体 |
etotfor | 同上 | 2561 | 总能 + 力 |
scprqt | src/67_common/m_common.F90 | 183 | 收敛判据 |
vtorho | src/79_seqpar_mpi/m_vtorho.F90 | 全 | 由 V 求 ρ + k 点循环 |
vtowfk | src/79_seqpar_mpi/m_vtowfk.F90 | 177 | 单 k 点对角化器分发 |
cgwf | src/67_common/m_cgwf.F90 | 197 | 标准共轭梯度 |
lobpcgwf2 | src/79_seqpar_mpi/m_lobpcgwf.F90 | 82 | LOBPCG 高层 |
lobpcg_run | src/48_diago/m_lobpcg2.F90 | 340 | LOBPCG 核心算法 |
chebfi | src/66_wfs/m_chebfi.F90 | 102 | Chebyshev 滤波 |
rmm_diis | src/66_wfs/m_rmm_diis.F90 | — | RMM-DIIS |
getghc | src/66_wfs/m_getghc.F90 | — | $H\psi$ 算子（FFT + nonlop） |
cg_precon | src/44_abitools/m_cgtools.F90 | 3651 | Teter 多项式预条件 |
ab7_mixing_new | src/56_mixing/m_ab7_mixing.F90 | 161 | 混合器初始化 |
ab7_mixing_eval | 同上 | 623 | 混合主入口 |
scfopt | 同上 | 1868 | 势混合算法（包括 Pulay/Anderson） |
xg_RayleighRitz | src/45_xgTools/m_xg_ortho_RR.F90 | — | Rayleigh-Ritz 子空间对角化 |
SCF 收敛全景图（一图流总结）
┌──────────────────────────────────────────────────────────────────────────┐
│ SCF 外循环（scfcv_core） │
│ │
│ ┌──────────┐ ┌──────────┐ ┌──────────┐ ┌──────────┐ │
│ │ V_in │ → │ vtorho │ → │ etotfor │ → │ scprqt │ │
│ │ (势) │ │ 求 ψ, ρ │ │ E, F, σ │ │ 收敛？ │ │
│ └──────────┘ └──────────┘ └──────────┘ └──────────┘ │
│ ↑ │ │ │
│ │ │ │ │
│ │ ┌────────────────────────────────┐ │ │
│ │ │ 内循环（vtowfk） │ ↓ │
│ │ │ │ 收敛？ │
│ │ │ 对每个 (isppol, ikpt): │ / \ │
│ │ │ CG / LOBPCG / ChebFi │ 是 否 │
│ │ │ │ │ │ │
│ │ │ 每条/每块 band: │ ↓ │ │
│ │ │ getghc → 预条件 → RR │ 退出 │ │
│ │ └────────────────────────────────┘ │ │
│ │ ↓ │
│ │ ┌────────────────────────────────┐ ┌────────────────┐ │
│ │ │ ρ_out → vxc, vhartre → V_out │ │ ab7_mixing_eval│ │
│ │ └────────────────────────────────┘ │ Pulay/Anderson │ │
│ │ │ └────────────────┘ │
│ └──────────────────────────┴──────────────────────┘ │
│ │
└──────────────────────────────────────────────────────────────────────────┘
最贵的步骤：getghc (FFT × 2 + nonlop)
收敛速度由：预条件器 + 混合算法 决定
鲁棒性由：tolerance + 历史长度 决定报告生成于 2026-05-19，基于 Abinit 10.4.7。
References
Anderson, Donald G. "Iterative Procedures for Nonlinear Integral Equations." Journal of the ACM, vol. 12, no. 4, 1965, pp. 547–560. link
Gonze, Xavier, et al. "The Abinit Project: Impact, Environment and Recent Developments." Computer Physics Communications, vol. 248, 2020, p. 107042. link
Knyazev, Andrew V. "Toward the Optimal Preconditioned Eigensolver: Locally Optimal Block Preconditioned Conjugate Gradient Method." SIAM Journal on Scientific Computing, vol. 23, no. 2, 2001, pp. 517–541. link
Levitt, Antoine, and Marc Torrent. "Parallel Eigensolvers in Plane-Wave Density Functional Theory." Computer Physics Communications, vol. 187, 2015, pp. 98–105. link
Payne, M. C., et al. "Iterative Minimization Techniques for Ab Initio Total-Energy Calculations: Molecular Dynamics and Conjugate Gradients." Reviews of Modern Physics, vol. 64, no. 4, 1992, pp. 1045–1097. link
Pulay, Peter. "Convergence Acceleration of Iterative Sequences. The Case of SCF Iteration." Chemical Physics Letters, vol. 73, no. 2, 1980, pp. 393–398. link
Teter, Michael P., Michael C. Payne, and Douglas C. Allan. "Solution of Schrödinger's Equation for Large Systems." Physical Review B, vol. 40, no. 18, 1989, pp. 12255–12263. link
The Abinit Code. Version 10.4.7, 2024. link