跳转至

14. 化简 4：Kohn-Sham 方程——从 N 体耦合到 N 个单体方程

本篇对应化简链的第 4 步。HK 定理（化简 3）告诉我们"找 \(\rho_0\) 就够了"，但没告诉我们 \(E[\rho]\) 长什么样。Kohn-Sham 1965 年的构造让 DFT 从一个存在性定理变成一个可在计算机上跑的算法。

1. 化简位置

化简 3：HK 定理         ρ_0 承载一切；E[ρ] 存在但形式未知
化简 4：Kohn-Sham 映射  ★ 本篇 ★  把 E[ρ] 落地成 N 个可解的单体方程
化简 5：XC 泛函近似     把残留的未知 E_xc[ρ] 用 LDA/GGA/… 模型近似

Kohn-Sham 的核心思想可以一句话概括：构造一个虚构的非相互作用电子系统，让它的密度等于真实系统的密度。 对虚构系统求解容易（单粒子方程），然后把密度拿出来用。

2. 上一步的困难

化简 3（HK）留下来的问题是：

$\(\rho_0 = \arg\min_\rho E[\rho], \qquad E[\rho] = T[\rho] + V_{ext}[\rho] + J[\rho] + E_{xc}[\rho]\)$ - \(V_{ext}[\rho]\) 和 \(J[\rho]\) 是简单积分，好算 - \(T[\rho]\) 呢？真实 N 电子动能不能用 ρ 显式写出来。早期 Thomas-Fermi 用 \(T^{TF}[\rho] \propto \int \rho^{5/3} d^3r\) 硬写，但误差巨大——甚至不能束缚原子成键，精度远远不够 - \(E_{xc}[\rho]\) 精确形式未知 核心症结：动能这个量，对波函数是个干净的微分算符（\(-\frac{1}{2}\nabla^2\)），但对密度没有干净表达。

ML 类比：如果你的 loss 有一项是输入梯度的某种统计量（比如 \(\int \|\nabla x\|^2\)），而你只能访问 \(x\) 的边缘分布而非 \(x\) 本身，就会很难算。Kohn-Sham 的解法相当于：偷偷保留一组辅助变量（KS 轨道）用来算梯度，对外仍然声称"我只用 ρ"。

3. 引入的假设/近似

核心构造（Kohn-Sham 映射）

假设：存在一个虚构的非相互作用电子系统，在某个有效势 \(V_{eff}(\vec{r})\) 下，其基态密度恰好等于真实相互作用系统的基态密度。

对虚构系统，电子之间没有相互作用，N 电子基态就是 N 个最低能单粒子态填满：

$\(\left[ -\tfrac{1}{2} \nabla^2 + V_{eff}(\vec{r}) \right] \phi_i^{KS}(\vec{r}) = \varepsilon_i \, \phi_i^{KS}(\vec{r})\)$ 这就是 Kohn-Sham 方程。密度直接由占据轨道构造：

\[\rho(\vec{r}) = \sum_{i=1}^{N} |\phi_i^{KS}(\vec{r})|^2\]

（自旋极化情况分开求和 \(\rho_\uparrow, \rho_\downarrow\)。）

有效势的构造

\(V_{eff}\) 怎么取？让虚构系统的密度等于真实系统密度、且满足 \(E[\rho]\) 的变分条件。答案是：

\[V_{eff}(\vec{r}) = V_{ext}(\vec{r}) + V_H(\vec{r}) + V_{xc}(\vec{r})\]

三项的物理含义：

项 表达式 物理
外势 $V_{ext}(\vec{r}) = -\sum_I \frac{Z_I}{ \vec{r}-\vec{R}_I
Hartree 势 $V_H(\vec{r}) = \int \frac{\rho(\vec{r}')}{ \vec{r}-\vec{r}'
XC 势 \(V_{xc}(\vec{r}) = \frac{\delta E_{xc}[\rho]}{\delta \rho(\vec{r})}\) 交换-关联势，化简 5 要近似的对象

能量分解

对应地，总能量泛函可以写成：

\[E[\rho] = T_s[\rho] + V_{ext}[\rho] + J[\rho] + E_{xc}[\rho]\]

注意 \(T[\rho]\)（真实动能）被替换成 \(T_s[\rho]\)（non-interacting kinetic energy）：

\[T_s[\rho] = \sum_{i=1}^{N} \langle \phi_i^{KS} | -\tfrac{1}{2}\nabla^2 | \phi_i^{KS} \rangle\]

两者不相等！它们的差值被 归入 \(E_{xc}\)：

\[E_{xc}[\rho] \equiv \underbrace{(T[\rho] - T_s[\rho])}_{\text{动能关联修正}} + \underbrace{(V_{ee}[\rho] - J[\rho])}_{\text{交换 + 库仑关联}}\]

这就是为什么 \(E_{xc}\) 叫做交换-关联能——它打包了两种"非平凡"的量子多体效应。

代价

  1. 引入虚构系统：KS 轨道 \(\{\phi_i^{KS}\}\) 和真实电子波函数没有直接物理对应，它们只是"用来造出正确密度的辅助数学对象"
  2. \(E_{xc}\) 仍然未知：化简 4 没解决 \(E_{xc}\) 问题，只是把所有"复杂动能关联 + 电-电关联"都塞到这一项里。需要化简 5 进一步近似
  3. 自洽循环：\(V_{eff}\) 依赖 \(\rho\)，\(\rho\) 又由 KS 方程的解给出。鸡生蛋问题——需要化简 11 的 SCF 迭代解决
  4. \(\varepsilon_i\) 的物理解释有限：KS 本征值不是真实电子能级，严格只有最高占据轨道能量 \(\varepsilon_{HOMO}\) 对应真实负电离能（Janak 定理/Koopmans-like），其他都是近似解读

4. 引入的新概念

概念 符号 含义
KS 轨道 / Kohn-Sham orbital \(\phi_i^{KS}(\vec{r})\) 虚构非相互作用系统的单粒子波函数，构成正交归一集合
KS 本征值 \(\varepsilon_i\) KS 方程的本征值，常被（近似地）当作电子能级画能带
有效势 / effective potential \(V_{eff}(\vec{r})\) 外势 + Hartree + XC 三者之和
Hartree 势 \(V_H(\vec{r})\) 经典电-电平均场
XC 势 / exchange-correlation potential \(V_{xc}(\vec{r})\) \(E_{xc}\) 对 ρ 的泛函导数
非相互作用动能 \(T_s[\rho]\) KS 轨道给出的动能，非真实动能
占据数 / occupation \(f_i\) 每个 KS 轨道被占据的程度（绝缘体 0 或 2，金属用 smearing 取分数）
HOMO / LUMO — Highest Occupied / Lowest Unoccupied Molecular Orbital；固体里对应价带顶 VBM 和导带底 CBM

5. 对应 QE 字段

化简 4 的每一件事在 QE 里都有具体对应：

化简 4 的概念 QE 中的位置
解 KS 方程 每个 SCF 迭代的核心操作，矩阵对角化（diagonalization='david' 默认 Davidson，或 'cg'）
KS 轨道 \(\phi_i^{KS}\) 存到 outdir/<prefix>.save/wfc*.dat（或 .hdf5），每个 k 点一个文件
KS 本征值 \(\varepsilon_i\) SCF 输出里的 k = ... 段落：每个 k 点下列出 nbnd 个本征值
占据数 \(f_i\) 绝缘体固定；金属由 occupations='smearing' + smearing=... + degauss 决定
有效势 \(V_{eff}\) 隐式构造；可以用 pp.x + plot_num=1 导出可视化
Hartree 势 \(V_H\) 隐式；能量项 hartree contribution 输出
XC 势/能 \(V_{xc}, E_{xc}\) 由赝势文件里的 <PP_HEADER functional="PBE"> 字段隐式选择；化简 5 展开
非相互作用动能 \(T_s\) 输出 kinetic energy = ... Ry（注意：这是 \(T_s\)，不是真实动能）
轨道数 nbnd &SYSTEM nbnd = ...，默认 \(N/2\)（绝缘体）或更多（金属需空带）
Fermi 能级 金属/smearing 下输出 the Fermi energy is ...

一次 SCF 迭代 = 解一次 KS 方程 + 更新密度 + 混合。化简 4 的数学结构直接决定了 pw.x 运行时 CPU 做的大部分工作。

6. 对应 benchmark 条目

benchmark 字段 与化简 4 的关系
total_energy_ev_per_fu \(E[\rho_0] = T_s + V_{ext} + J + E_{xc}\) 在收敛密度处的值
band_gap（若评分涉及） 用 KS 本征值的 CBM - VBM 估算（严格来说低估真实带隙，见化简 5）
relaxed_structure 弛豫过程中每个构型都要解一次 KS 方程得到 \(E(\{\vec{R}_I\})\) 和力
conv_thr &ELECTRONS conv_thr 控制 KS 方程 SCF 停机阈值；ground_truth 生成要求足够紧

7. ML 类比

类比 1：Mean-field / variational inference

Kohn-Sham 是相互作用系统 → 非相互作用辅助系统的映射。这和变分推断里"用一个易处理分布 \(q\) 去逼近难处理后验 \(p\)"高度神似：

Kohn-Sham 变分推断
真实 N 电子系统（耦合，难） 真后验 $p(z
虚构非相互作用系统（单体，易） 变分分布 \(q_\phi(z)\)（易）
约束：密度 ρ 相等 约束：最小化 KL(q || p)
自洽通过 \(V_{eff}[\rho]\) 实现 自洽通过参数 \(\phi\) 更新
残差放进 \(E_{xc}\) ELBO 的 gap 放进损失

类比 2：把 N 体问题分解成 N 个单体 + 共享环境

一个更直白的类比：想象 N 个 agent 要优化一个含两两交互项的目标函数。Mean-field 办法是：每个 agent 不看其他 agent 的具体位置，只看所有 agent 的密度 ρ 形成的平均场，在这个平均场里做自己的优化。

真实问题:    N 个 agent × N-1 个两两交互项  (O(N²))
KS 分解:     N 个 agent × 1 个共享环境 V_eff[ρ]  (O(N) 环境 eval + O(N³) 对角化)
                                         +
                        ρ = Σ |φ_i|²  的自洽一致性要求

共享环境 \(V_{eff}[\rho]\) 是所有 agent"共同造出来"又"共同响应"的——这正是 mean-field 游戏的经典结构。

类比 3：KS 轨道类似"特征表示（representation）"

KS 轨道 \(\{\phi_i^{KS}\}\) 本身不是物理可观测量，但由它们构造出来的密度 ρ 是。这类似于：

  • 神经网络的隐藏层激活（hidden activations）本身不直接有物理意义
  • 但由它们构造出的输出（predictions）是可测可比的

KS 轨道是"密度的高效参数化"。可以把 Kohn-Sham 方程看成：寻找一组正交归一的 3 维函数 \(\{\phi_i\}\)，让它们构成的密度 \(\sum_i |\phi_i|^2\) 正好最小化能量泛函。

类比 4：降维的代价——丢弃了相关性信息

从 N 电子波函数 \(\Psi\) 到 N 个独立 KS 轨道 \(\{\phi_i^{KS}\}\)，相当于丢弃了电子之间的量子关联。但真实电子显然是关联的（两个电子不可能同时出现在同一点，除了自旋不同）。

这个关联被补偿到 \(E_{xc}\) 里。所以 \(E_{xc}\) 不只是"交换能"（源自反对称性），还包含"关联能"（源自关联效应）。

8. 典型取值 / 常见坑

坑 1：KS 本征值 ≠ 真实电子能级

常见误用：把 PBE 算出的 KS 本征值之差当作真实带隙，然后和实验比。这是系统性低估： - Si 实验带隙 ≈ 1.17 eV - PBE KS gap ≈ 0.5-0.6 eV（低估约 50%）

原因：KS 系统是虚构的非相互作用系统，其"LUMO"不对应真实电子加一的激发能。严格处理需要 GW 近似、DFT+U、杂化泛函（见化简 5 的 HSE06）或 TDDFT。

坑 2：自洽迭代不一定收敛

\(V_{eff}[\rho] \to \phi_i \to \rho \to V_{eff} \to \ldots\) 是一个固定点迭代。原始直接代入几乎肯定发散，需要混合（mixing）技术（Broyden、Anderson、Pulay）——这是化简 11 的主题。

典型难收敛体系： - 金属：Fermi 面附近态密度大，\(\rho\) 对 \(V_{eff}\) 微扰敏感 - 过渡金属氧化物：强关联 + 多个能量相近的磁构型 - 大超胞 + 低对称性：相空间大，混合步长要小

坑 3：nbnd 要开够

nbnd 决定求多少个 KS 本征值。默认只开占据轨道： - 绝缘体：\(N_{val}/2\) 个轨道（非自旋极化） - 金属/smearing：需要 extra bands（通常 20% 额外），否则 Fermi level 附近填充不稳

如果要算非占据态（能带图的导带、DOS 的高能部分），必须在 nscf 里手动加大 nbnd。

坑 4：占据方式选错导致金属被算成绝缘体（或反之）

  • 绝缘体用 occupations='fixed'（默认），要求 nbnd 恰好使能隙分明
  • 金属必须用 occupations='smearing'，否则 SCF 一般不收敛（Fermi 面上态跳跃）
  • 常见 smearing：smearing='mv'（Marzari-Vanderbilt，固体推荐）、smearing='gaussian'、smearing='fd'（Fermi-Dirac）
  • degauss 典型 0.01-0.02 Ry；太大扭曲能量，太小难收敛

坑 5：spin 要按体系设好

  • 非磁体系：nspin=1（KS 方程一套，每个轨道占据 0 或 2）
  • 共线磁性：nspin=2（KS 方程两套，\(\phi^\uparrow\) 和 \(\phi^\downarrow\)），并设初始磁矩 starting_magnetization
  • 非共线 / SOC：nspin=4 + lspinorb=.true.，KS 轨道变成两分量旋量

典型数值感受（Si 体系）

一次 Si vc-relax 的 SCF 迭代： - ecutwfc = 40 Ry，原胞 2 原子 8 价电子 - nbnd 默认 4（价带）；要看导带要手动加大 - 10-20 次 Davidson 对角化内可达 conv_thr = 1e-8 - KS 本征值前 4 个在 -5 ~ +5 eV 范围，第 5 个（LUMO）~0.6 eV 高于 HOMO → PBE 带隙约 0.6 eV

下一步阅读