22. 化简 12：BFGS 结构弛豫¶

1. 化简位置¶

12 步化简链的最后一步。至此我们完成了所有电子结构部分：化简 1-10 构造了 KS 方程的有限维数值形式，化简 11 用 SCF 解出给定核位置下的电子基态和总能量 \(E(\{\vec{R}_I\}, \{\vec{a}_j\})\)。

但这还不是物理基态——核位置和晶格矢量也要优化到能量最低。这就是本步做的事：几何优化 (geometry optimization)，在 QE 里叫结构弛豫 (structural relaxation)。

2. 上一步的困难：SCF 只给一个点的能量¶

SCF 解出的是 给定核位置 下的电子基态。但是：

Born-Oppenheimer 势能面 \(E(\{\vec{R}_I\}, \{\vec{a}_j\})\) 是核坐标 + 晶格参数的函数
真实晶体的"基态结构"对应这个函数的最小值
初始结构（比如从实验值 / Materials Project 猜来的）一般不是最小值
我们要外层再做一个优化：找使能量最低的核位置和晶格

这是一个双层优化问题 (bilevel optimization)：

外层（化简 12）：min_{R, a} E(R, a)
    内层（化简 11）：E(R, a) = SCF({R, a}) 的收敛总能量

ML 类比直接到位：像 MAML（内层训练模型，外层调超参）、像 implicit differentiation 的结构。

3. 引入的假设 / 解法¶

3.1 Hellmann-Feynman 定理给梯度¶

外层优化需要梯度。幸运的是 Hellmann-Feynman 定理告诉我们：对已经在内层 SCF 收敛的波函数 \(\phi_i\)，能量对核坐标的梯度只需要对 H 里的显式参数依赖求导：

\[\vec{F}_I = -\frac{\partial E}{\partial \vec{R}_I} = -\langle \Psi | \frac{\partial \hat{H}}{\partial \vec{R}_I} | \Psi \rangle\]

注意关键点：不需要对 SCF 迭代本身反传。这是因为在驻点处，\(\partial E / \partial \phi_i \cdot \partial \phi_i / \partial R = 0\)——波函数对能量变分稳定。

ML 类比：这是 envelope theorem / implicit function theorem 的物理版本。类似你在 bilevel 里用 argmin 的 KKT 条件换掉对 argmin 的直接求导。

类似地有应力张量（能量对晶格应变的梯度）：

\[\sigma_{ij} = -\frac{1}{V} \frac{\partial E}{\partial \epsilon_{ij}}\]

QE 每次 SCF 收敛后都会输出 Forces 和 total stress。

Hellmann-Feynman 成立的前提：内层 SCF 必须严格收敛。如果 conv_thr 太松，\(\partial E/\partial \phi \cdot \partial \phi/\partial R \neq 0\) 的残差项会污染力 → BFGS 会走错方向。所以 vc-relax 时 conv_thr 一般比纯 SCF 更严 1-2 个量级。

3.2 BFGS 拟牛顿法¶

外层用 BFGS (Broyden-Fletcher-Goldfarb-Shanno)：

维护一个近似的逆 Hessian \(B_k^{-1}\)
每步根据 \((x_{k+1} - x_k, g_{k+1} - g_k)\) 做 rank-2 更新
更新方向 \(d_k = -B_k^{-1} g_k\)
收敛速度超线性（比最速下降快一个数量级）

对光滑的能量面（常温下远离相变的晶体都满足），BFGS 几步就够。

ML 读者应该很熟：PyTorch 的 torch.optim.LBFGS 就是内存受限版的 BFGS。QE 的 BFGS 不是 limited-memory（自由度通常只有几十到几百，不需要省内存）。

3.3 三种 calculation¶

值	变量	典型用途
`relax`	只动原子位置，晶格固定	已知晶格常数（实验值），优化内部坐标
`vc-relax`	variable-cell relax，同时动原子位置和晶格	从头预测晶格常数和原子位置
`md`	分子动力学	有限温度采样，不是静态优化

DFTBench 绝大部分 benchmark 用 vc-relax。

4. 引入的新概念¶

4.1 收敛判据¶

BFGS 在满足所有下述条件时停止：

力收敛：每个原子 \(|\vec{F}_I| <\) forc_conv_thr（Ry/bohr）
能量收敛：最近两步 \(|E_{new} - E_{old}| <\) etot_conv_thr（Ry）
应力收敛（仅 vc-relax）：\(|\text{trace}(\sigma)| / 3 <\) press_conv_thr（kbar）

注意单位：QE 力的默认单位是 Ry/bohr（非 eV/Å）。

4.2 `cell_dofree`：对称性约束¶

vc-relax 里晶格有 6 个自由度（\(a, b, c, \alpha, \beta, \gamma\)），但对称性常常约束部分自由度。cell_dofree 选项：

值	含义
`'all'`	全 6 个自由度（默认）
`'ibrav'`	保持 Bravais 类型（cubic 只变 a）
`'xyz'`	只变 3 个长度，角度保持正交
`'xy'` / `'xz'` / `'yz'`	二维约束
`'2Dxy'`	二维体系（保持 z 方向真空层）
`'volume'`	只变体积（等比缩放）
`'shape'`	保持体积，只变形状

用对称性约束能减小搜索空间、加速收敛，并避免 BFGS 破坏晶体对称性（数值噪声会导致立方被"错误地"优化成正交）。

4.3 `cell_factor`¶

vc-relax 里晶格在变，平面波基组（依赖倒格矢）也在变——但 QE 在起始时分配一个固定大小的 G 向量空间。cell_factor（默认 1.2）表示预留 20% 的体积膨胀空间。如果优化中晶格变化超过 20%，会报错。

对大形变体系（如你之前 Ge 的情形）可能需要 cell_factor=1.5 甚至 2.0。

4.4 `press`：目标压强¶

默认 press=0.0（kbar），即大气压 0 压力。设置正值可做高压弛豫（模拟地核、金刚石压砧实验）。ML 视角：这是一个拉格朗日约束—— 对 \(E + PV\) 求极值。

5. 对应 QE 字段¶

`&CONTROL`¶

字段	含义
`calculation='vc-relax'`	变胞弛豫
`calculation='relax'`	固定晶格的弛豫
`etot_conv_thr`	总能量收敛阈值（Ry），默认 1e-4
`forc_conv_thr`	力收敛阈值（Ry/bohr），默认 1e-3

`&IONS`（原子位置优化）¶

字段	含义	典型值
`ion_dynamics`	`'bfgs'`（默认）/ `'damp'`（阻尼 MD）/ `'verlet'`	`'bfgs'`
`upscale`	SCF ethr 随 BFGS 步数的加严因子	100（默认）

`&CELL`（晶格优化，仅 vc-relax）¶

字段	含义	典型值
`cell_dynamics`	`'bfgs'`（默认）/ `'none'`	`'bfgs'`
`press`	目标压强（kbar）	0.0
`press_conv_thr`	压强收敛（kbar）	0.5
`cell_dofree`	晶格自由度约束	`'all'`
`cell_factor`	G 向量空间预留比	1.2

输出中的 BFGS 监控¶

BFGS Geometry Optimization

number of scf cycles    =   3
number of bfgs steps    =   2

energy              =  -15.84528512 Ry
new trust radius    =   0.050 bohr
new conv_thr        =   1.0E-10 Ry

Forces acting on atoms (cartesian axes, Ry/au):
     atom    1 type  1   force =  0.00012 -0.00008  0.00003
     ...
     Total force =     0.000157     Total SCF correction =     0.000002

     total   stress  (Ry/bohr**3)                   (kbar)      P=       0.32

每一步 BFGS 触发一次完整 SCF（内层 10-50 步）。

6. 对应 benchmark 条目¶

DFTBench 里 vc-relax 任务是主力评测，涉及字段：

输入模板：questions/vc_relax*.json
输出 CELL_PARAMETERS → 评分字段 a, b, c, α, β, γ（relative error）
输出 ATOMIC_POSITIONS → 评分字段 relaxed_structure
对称性分析（用 pymatgen / spglib 后处理）→ space_group, space_group_number, point_group, crystal_system（exact match）
最终 SCF 总能量 → total_energy_ev_per_fu

精度档（见 evaluate/compare.py 和 benchmark 文档）：

档位	`forc_conv_thr`	`press_conv_thr`	能量精度要求
low_acc	1e-3	0.5 kbar	20 meV/atom
medium_acc	1e-4	0.1 kbar	10 meV/atom
high_acc	1e-5	0.05 kbar	1 meV/atom

注意：benchmark ground_truth 当前全空（见项目级 README），所以这些阈值目前只是形式上存在。

7. ML 类比速查¶

DFT 概念	ML 对应
SCF + BFGS 双层结构	bilevel optimization / MAML（inner 训模型，outer 调超参）
BFGS	拟牛顿法，= PyTorch 的 `torch.optim.LBFGS` 的一般版本
Hellmann-Feynman 定理	envelope theorem / implicit function theorem（不需要对 argmin 求导）
力 \(\vec{F}_I\)	损失对参数的梯度
`forc_conv_thr`	梯度范数的 early stopping 阈值
`etot_conv_thr`	loss 变化的 early stopping 阈值
`cell_dofree`	参数约束（类似 weight tying 或对称性正则）
陷入局部极小	非凸优化的经典陷阱
BFGS 走错方向（SCF 未收敛）	梯度噪声导致优化器发散
trust radius（QE 自适应步长）	trust-region method

8. 典型取值与常见坑¶

典型参数¶

精度/场景	conv_thr (SCF)	forc_conv_thr	press_conv_thr	etot_conv_thr
粗调	1e-6	1e-3	0.5 kbar	1e-4
标准	1e-8	1e-4	0.1 kbar	1e-5
高精度	1e-10	1e-5	0.05 kbar	1e-6

典型 BFGS 步数¶

初始结构接近最小（如从实验值出发）：3-5 步
中等偏离（Materials Project 猜初值）：5-15 步
大偏离（瞎猜）：20+ 步，可能卡住

你之前跑的 Si（ibrav=2, 实验值初猜）3 步就收敛，典型。Ge 出现过 6% 体积压缩问题——那不是 BFGS 的锅，而是初始结构 + 泛函选择导致 BFGS 走到了非物理局部极小。

常见坑¶

BFGS 陷入非物理局部极小
症状：弛豫后晶格收缩/膨胀严重（>5%），空间群错
原因：初始结构太差，或 XC 泛函（如 LDA）本身有系统偏差
解决：从 experimental 结构初猜；换泛函（PBE → PBEsol）；对比 Materials Project 的 structure_data
SCF 阈值太松 → 力不准 → BFGS 振荡
症状：BFGS 能量下降后反弹，或力一直在 1e-3 附近徘徊
原因：conv_thr 设得不够严，Hellmann-Feynman 残差项污染力
解决：conv_thr ≤ forc_conv_thr / 100（经验法则）
cell_factor 太小
症状：运行中报 "G-vectors exceed limit" 错误
解决：cell_factor=1.5 或更大
压强一直不收敛
原因：DFT 计算的压强绝对精度就 \(\sim\) 1 kbar（来自 XC 泛函和 ecutrho 的误差），press_conv_thr=0.01 纯属浪费
解决：不要把 press_conv_thr 设得比 ecut / 泛函本身的精度还严
对称性被破坏
症状：cubic 材料弛豫成 orthorhombic（三个轴长度略有差异）
原因：BFGS 数值噪声没有被对称性约束
解决：cell_dofree='ibrav'（保持 Bravais 类型）；或 post-process 时用 spglib 做对称化
ion_positions vs cell 自由度耦合
vc-relax 每一步 BFGS 同时更新原子位置和晶格，有时能量曲线不光滑
可以分两阶段：先 relax 优化原子位置，再 vc-relax
结构弛豫后要重跑 SCF
vc-relax 结束时的 SCF 是在旧基组下（G 向量数随晶格变化但 cell_factor 预留的空间固定）
严格报告能量/能带时，用最终结构重新跑 calculation='scf'（和后续 nscf/bands）

下一步阅读¶

至此 12 步化简链全部走完。下一篇转到"DFT 输出端"：

30-observables.md — DFT 能算的物性（能量、能带、DOS、磁矩、弹性、介电、声子……）

如果想复习 SCF 内层：

21-scf-iteration.md — 本篇依赖的内层循环