DFT 计算中的优化算法与 Mixing 算法详解

核心计算步骤分类

DFT 计算步骤主要可分为两类：

1. 处理显式表达式：$B = f(A)$，已知 $A$ 即可直接求出 $B$。例如矩阵乘法、积分、插值、傅里叶变换等，均具有明确的表达式。
2. 处理隐式约束：$f(B) = A$，逆函数难以求解，甚至根本不存在显式表达。例如方程组、不动点问题、极小值问题等。

说明：不同 DFT 程序中涉及的显式表达式及其具体算法差异较大，难以一一列举。然而，处理隐式约束的方式在几乎所有 DFT 程序中都具有共性，因此本文将重点介绍隐式约束的处理方法。

在 DFT 计算流中，常见的隐式约束及相应的优化算法包括：

• 混合算法（Mixing）：物理量自洽，$\rho \to H(\rho) \to \psi \to \rho$
• 对角化算法（Diagonalization）：定态薛定谔方程，$H\psi - ES\psi = 0$
• 结构弛豫算法（Relaxation）：总能量极小，受力为零（即 $F=0$）

说明：DFT 计算的基本流程是一个自洽循环：哈密顿量 $H$ 依赖于电荷密度 $\rho$，波函数 $\psi$ 由 $H$ 解出，而 $\rho$ 又可由 $\psi$ 给出，形成循环。DFT 最终要达到的是经过一次循环后，电荷密度不再发生变化（即自洽）。

寻找 Kohn-Sham 方程的自洽解时，从试探电荷密度 $\rho^{\text{in}}$ 出发，执行：

$$\rho^{\text{in}} \to H \to \text{DM} \to \rho^{\text{out}}$$

隐式约束关系为：$\hat{\rho} \xleftarrow{H} \psi$

混合算法（Mixing）

一般情况下 $\rho^{\text{out}} \neq \rho^{\text{in}}$，直接将 $\rho^{\text{out}}$ 输入下一轮循环的效果并不理想。

Mixing 的核心思想：将 $\rho^{\text{out}}$ 与 $\rho^{\text{in}}$ 进行混合，给出下一步的 $\rho^{\text{in}}$：

$$\rho_n^{\text{in}} \xrightarrow{n^{\text{th}} \text{loop}} \rho_n^{\text{out}} \xrightarrow{\text{mix}} \rho_{n+1}^{\text{in}} \xrightarrow{n+1^{\text{th}} \text{loop}} \dots$$

注：在原子基程序中，除了电荷密度 $\rho$ 外，还支持密度矩阵（DM）和哈密顿量（H）的 mixing。虽然 mixing 操作可以放在哈密顿量或密度矩阵上，但公式形式是类似的，因此通常以电荷密度 $\rho$ 为例进行讲解。实际流程中，$\rho^{\text{in}}$ 和 $\rho^{\text{out}}$ 都要保存，最后经过 mix 给出一个新的 $\rho^{\text{in}}$，再输入到 $\rho \to H$ 的计算中。

Linear Mixing（线性混合）

$$\rho_{n+1}^{\text{in}} = \rho_n^{\text{in}} + G R_n$$

• 残差（Residual）：$R_n = \rho_n^{\text{out}} - \rho_n^{\text{in}}$
• 权重 $G$：人为指定，一般 $G \in (0,1)$，且与体系的介电常数 $\epsilon$ 有关（$G \sim 1/\epsilon$）。
• 特点：仅利用当前的 $\rho^{\text{in}}$ 和 $\rho^{\text{out}}$ 两个信息，收敛较慢。

说明：Linear mixing 仅比直接将 $\rho^{\text{out}}$ 作为下一步输入略好一点。权重 $G$ 的最佳值因体系而异（不同体系的介电函数不同），而在求解之前你不知道介电函数的具体形式，因此总体收敛很慢。而且如果权重给得不好，很多情况下完全无法收敛。因此这种方法在实际中基本不使用。

Pulay Mixing

说明：Pulay Mixing 又称 DIIS（Direct Inversion in the Iterative Subspace）。

Pulay Mixing 是最小化残差法的一种。第 $n$ 步 mixing 时，利用过去 $N$ 步的信息 ${ \rho_m^{\text{in}}, R_m \mid m = n-N+1 \sim n }$，在 ${ \rho_m^{\text{in}} }$ 张成的空间中搜索，找到残差最小所对应的 $\rho_n^{\text{min}}$，再进行 linear mixing：

$$\rho_{n+1}^{\text{in}} = \rho_n^{\text{min}} + G R_n^{\text{min}}$$

• 子空间搜索：$\rho = \sum_i c_i \rho_i^{\text{in}}$（线性叠加），且满足归一化条件 $\sum_i c_i = 1$（因为总电荷是定值，不能改变总量）。

• 线性近似：$R(\rho) = \sum_i c_i R_i$（在接近收敛时，残差近似为线性关系，这是一个合理的假设）。

• 最小化残差：求 $\min_\rho \langle R(\rho) | R(\rho) \rangle$，解得系数：

  $$c_i = \frac{\sum_j A_{ij}^{-1}}{\sum_{ij} A_{ij}^{-1}}$$

  其中 $A_{ij} = \langle R_i | R_j \rangle$，$\langle \cdot | \cdot \rangle$ 是内积计算（最简单的内积是把残差的各个分量相乘再相加，也可以设计更高效的内积方式来提升收敛速率）。

• 最终更新公式：

  $$\rho_{n+1}^{\text{in}} = \sum_i c_i (\rho_i^{\text{in}} + G R_i)$$

• 优缺点：

  • 优点：收敛速度快，对 $G$ 的依赖小，是 DFT 程序中最广泛采用的算法。
  • 潜在缺陷：求 $A_{ij}^{-1}$ 时涉及矩阵求逆，当 ${ R_m }$ 线性相关时，矩阵奇异，数值上不稳定。

直观理解：假设收敛后的 $\rho$ 是某个红点（但我们不知道它在哪）。我们尝试了两次，知道了 $\rho_1$、$\rho_2$ 以及两者的残差。做法是：首先将两者的残差做线性组合，找到残差最小的组合（蓝线所示）；得到叠加系数后，将这个系数运用到两个电荷密度上，在它们连线上找到最佳点；最终从这个最佳点出发，沿它所对应的残差向前走一步，得到 $\rho_3$。可以发现 $\rho_3$ 更接近收敛值（红点）。

Broyden 2nd Mixing

说明：Broyden Mixing 是对 Pulay Mixing 的改进，通过在矩阵求逆时人为加入对角元来解决线性相关问题。

Broyden 2nd Mixing 为改进 Pulay mixing 的数值稳定性而生。对内积矩阵 $A$ 求逆时，加入对角元 $w_0$ 以改善 ${ R_i }$ 线性相关时的数值稳定性，但同时会降低收敛速度。

每个历史步 $\rho_i^{\text{in}}, R_i$ 都赋予可调权重 $w_i$：

说明：当残差线性相关时，矩阵求逆会出问题。人为加入对角元后，原本线性相关的矩阵变得不相关了，求逆操作可以正常进行。但这破坏了之前"残差最小"的条件，有可能会降低收敛速度。

总结与对比

• 电子步 Mixing：$\rho$ 的自由度很大（实空间或 K 空间格点数约 $10^3 \sim 10^5$），自由度越大越不易产生线性相关，故 Pulay mixing 最合适。
• 结构弛豫 Mixing：自由度较少（$3N$，$N$ 为原子数），更易出现线性相关，此时 Broyden 2nd mixing 的稳定性优势得以体现。
• 注：在 Linear mixing 与 Pulay mixing 的对比图中，横轴是迭代步数，纵轴是误差（取 log）。不同结构用 linear mixing 收敛行为差异很大；但用了 Pulay mixing 之后，基本重合，约 9 步左右都收敛了，非常稳定。

对角化算法（Diagonalization）

求解 $H \to \psi$（进而得到 DM）的算法，本质是将本征值问题转化为变分优化问题。

说明：求波函数这一步不一定非要用对角化方法，也有非对角化方法，但这里只讲用得最多的对角化方法。

求解类型与相关程序包

• 矩阵对角化：稀疏矩阵（原子基，有限元程序）、稠密矩阵（平面波基程序）。
• $O(N)$ 算法：计算复杂度正比于原子数（准确说是轨道数），常用于原子基程序，如 Divide-conquer, Krylov subspace。但普遍缺点是精度不是非常高。
• 其它方法：Pole-expansion and selected inversion, Chebyshev polynomial。
• 相关程序包：
  • LAPACK + BLAS：线性代数包，可求解线性方程组、（广义）本征值问题、奇异值问题等，支持稀疏或稠密矩阵，包含多种求解器。第一性原理计算程序绝大多数都会调用这个包。
  • PEXSI：利用原子基程序中 $H$ 和 DM 的稀疏性，不对角化，而是选择性求逆，直接给出密度矩阵中的有用部分。
  • CHESS：利用 Chebyshev polynomial 的 filter 特性，过滤 $H$ 的高能部分，只保留占据态子空间，再在该子空间作对角化。

Davidson Scheme（DAV）

说明：Davidson Scheme 的共同思路是把本征值问题转化为变分优化问题——找到一个标量，使得当波函数满足本征值问题时，这个标量恰好也是极小的。选择不同的标量，就会导出不同的方法。

• 优化对象：能量期望 $\langle \psi_i | H | \psi_i \rangle$

• 优化目标：

  $$L_i = \langle \psi_i | H | \psi_i \rangle - \sum_j \lambda_{ij} (\langle \psi_j | S | \psi_i \rangle - \delta_{ij})$$

  （$S$ 是重叠矩阵，$|\psi_i\rangle$ 是试探波函数，有正交性限制 $\langle \psi_j | S | \psi_i \rangle = \delta_{ij}, j=1 \dots N$）

• 求解思路：最小化 $L_i$，可得极小值点 $|\psi_i^{\text{min}}\rangle$ 的表达式，但直接优化这个量表达式复杂、计算也复杂，不太可行。

• Preconditioner 解决方案：

  试探解 ${ |\psi_i\rangle, \epsilon_i }$ 的残差为：$|R\rangle = (H - \epsilon_i S) |\psi_i\rangle$。

  由 $|R\rangle$ 给出更新方向：$|\delta \psi_i\rangle = K |R\rangle$，使得 $|\delta \psi_i\rangle \approx |\psi_i^{\text{min}}\rangle - |\psi_i\rangle$。

  其中 $K$ 称为 preconditioner（预条件子），不同问题设定不同，一般依经验选取。

说明：残差和真正的误差有联系，但不一定恰好沿着误差方向。假设有一个算符 $K$ 作用到残差上，会给出更新方向（希望它沿着真正的误差方向更新）。$K$ 的严格形式我们不知道，不同问题合适的 $K$ 也不一样，但人们发现某种取法最好，就都用这种。

• 正交化（手动操作）：$K$ 给出的更新方向 $|\delta \psi_i\rangle$ 不一定满足正交性限制，需要手动正交化：

  $$|\delta \psi_i'\rangle = \left(I - \sum_{j<i} S |\psi_j\rangle \langle \psi_j|\right) |\delta \psi_i\rangle$$

  这相当于向其他向量之外的空间投影，保证新给出的向量和原来的子空间正交。

• VASP 中的 Preconditioner 示例：

  $$\langle G | K | G' \rangle = \delta_{GG'} \frac{27+18x+12x^2+8x^3}{27+18x+12x^2+8x^3+16x^4} \quad \text{且} \quad x = \frac{\hbar^2}{2m} \frac{G}{1.5 E_{\text{kin}}(R)}$$

  简化版：$\langle G | K | G' \rangle = \delta_{GG'} \max { G_{\text{cut}} / G^2, 1 }$

• 利用能量期望作更新：在二维子空间 ${ |\psi_i\rangle, |\delta \psi_i\rangle }$ 中对角化哈密顿量（2×2 算符，对角化非常容易），得到两个本征值，只取能量较低的态（希望能量越低越好）。用 $|\psi_i'\rangle$ 求出新的修正向量，构成三维子空间继续更新，直至 $N$ 步后开始更新下一个能级 $|\psi_{i+1}\rangle$。

• Subspace Rotation：更新全部能级之后，在基组 ${ |\psi_j\rangle, j=1 \dots N_b }$ 下统一对角化哈密顿量，以提升精度。

• 优缺点：

  • 优点：VASP 中的默认方法，非常稳定。不论初始解什么样，因为总是向低能处搜索，不会遗漏本征解，求出的本征解一定从低到高依次排列。
  • 缺点：正交化步骤必不可少，但要计算很多内积，计算量较大。

Residual Minimization Method（RMM）

说明：RMM 全称 RMM-DIIS，也叫 Pulay Mixing for Diagonalization。为了解决 DAV 中正交化导致计算慢的问题，希望不要正交化。

• 优化对象：残差模长 $\langle R | R \rangle$（残差优化到 0 自动满足正交条件）。

  $$L_i = \langle R(\psi_i) | R(\psi_i) \rangle - \lambda_i (\langle \psi_i | S | \psi_i \rangle - 1)$$

  （无需考虑正交条件 $\langle \psi_j | S | \psi_i \rangle = \delta_{ij}$）

说明：当 $\psi$ 是严格的本征解时，残差一定为零；残差模长是正定的，所以零就是极小值点。当残差优化到零时，自动满足正交条件，不用显式地强行指定正交化，大大简化了计算。

• 求解与更新：已知过去 $N$ 步的试探波函数及残差 ${ |\psi_i^m\rangle, |R_i^m\rangle, m=1 \dots N }$。最小化残差，得本征方程：

  $$A \vec{c} = L_i \bar{S} \vec{c}$$

  其中 $A_{mn} = \langle R_i^m | R_i^n \rangle, \bar{S} = \langle \psi_i^m | S | \psi_i^n \rangle$。

  解得 $\vec{c}$，即可得 $\psi_i^{\text{min}}$ 和 $R_i^{\text{min}}$。

  更新公式：$|\psi_i^{n+1}\rangle = |\psi_i^{\text{min}}\rangle + G |R_i^{\text{min}}\rangle$（$G$ 是步长）。更新所有能带后，作 subspace rotation。

说明：这个优化会得到另一个本征方程，但维度非常低（与过去步数 $N$ 有关，$N$ 一般不大），非常好对角化。

• 优缺点：
  • 优点：无需正交化，计算量较小。无需其它能级信息，适合并行（每条能带可以分配给一个核，各自收敛到对应的能级上）。
  • 缺点：由于不需要其他能级的信息，不知道别人的状态，有可能和另一个线程结果一致（到同一个态上），或者跑到更高的态上去。初始化时需小心，需要基本覆盖正确的占据态空间才能启动 RMM。

经验补充：

• Pulay Mixing 又叫 DIIS。
• VASP 之所以能打败市面上 99% 的同类计算软件，一般认为有两点原因：
  1. 基于 PAW（Projector Augmented Wave）平面波的赝势做得非常好，那个赝势是有版权的，不允许随便改。
  2. 它的算法依赖于 RMM-DIIS。其中 RMM 就是刚才提到的基于残差优化的密度算法。
• VASP 有一个模式叫 fast：因为 RMM 如果初始选得差，能级会选错（漏掉能级或重复算能级），所以一般计算方法是：先用 Davidson 算法算 5 步（一般 3-5 步），算出初步结果后，立刻改用 RMM 算法加速计算。这样效果非常好，是 VASP 成功的秘诀之一。

结构弛豫算法（Relaxation）

说明：Relaxation 算法要解决的问题是：对于每一个原子结构，都能算出一个总能以及受力。注意，这里的受力不是用有限差分等方法算出来的，而是可以解析给出的，总能和受力可以同时得到。受力就是能量的梯度。我们要利用这两个信息去寻找什么位置下总能最低。

受力就是能量梯度。利用总能 $E_{\text{tot}}$ 和受力 $\vec{F}$，寻找 $\min_{\vec{x}} E_{\text{tot}}(\vec{x})$。

（$\vec{x}$ 为所有原子坐标构成的向量，$\vec{F} = \frac{\partial E_{\text{tot}}}{\partial \vec{x}}$）

说明：总能是电子态贡献的能量与晶格离子的库仑能之和，依赖于离子位置。总能极小等价于受力为零的点。

梯度下降法（Steepest Descent）

说明：最 naive 的想法——既然受力是能量的梯度，为什么不沿着受力方向去寻找能量极小？沿受力方向走，能量肯定降低。

Naïve 的想法，沿受力方向（$F_n$）搜索至极小值：

$$F_{n+1} \cdot \delta X_n = 0, \quad \delta X_n = X_{n+1} - X_n = \alpha_n F_n$$

（$\alpha_n$ 为步长），然后再沿 $F_{n+1}$ 继续搜索。

• 缺陷：势能面的曲率对应声子频率，受力 $F_n$ 在各个声子模式上均有投影。高频模式（陡峭方向）限制了搜索步长（受力大，步长不能太大，否则跨过极小点到另一边），低频模式（平缓方向）减缓了收敛速率（受力小，每次走不了多少）。局部的梯度不是最优的搜索方向，缺乏全局信息导致收敛很慢，而且有震荡。

直观理解：想象二维势能面，一个方向非常陡，另一个方向非常缓。受力既有陡方向的分量，也有缓方向的分量。陡方向给出很大的受力，导致步长不能太大；而缓方向受力小，收敛速度受限于此。结果就是在陡峭方向来回折返，真正有用的缓方向却每次走不了多少。

Damped Molecular Dynamics（阻尼分子动力学）

说明：把受力和能量想象成一个动力学系统——一堆原子在受力作用下运动。但只让它们运动起来不行（没有耗散），所以要引入阻尼，让能量慢慢收敛到最低点。

利用速度 $V$ 记录历史信息，优化方向为 $V$ 而非 $F$；引入阻尼 $\mu$ 减少震荡（有速度引入能量耗散）。

• 运动方程：

  $$a_n = \alpha F_n - \mu V_n$$ $$V_{n+1/2} = V_{n-1/2} + a_n$$ $$X_{n+1} = X_n + V_{n+1/2}$$

• 解得：

  $$V_{n+1/2} = \frac{(1-\mu/2)V_{n-1/2} + \alpha F_n}{1+\mu/2}, \quad X_{n+1} = X_n + V_{n+1/2}$$

说明：速度由加速度更新，而加速度包含之前的受力，所以速度把之前几步的受力、之前的更新都记录到了当前的 $V$ 中。这就是它比梯度下降好的原因。在合适的参数下，收敛路径非常平滑，慢慢跑到势能面最低点。

• 注：神经网络中常用的优化器 Momentum, Adam 等都是该方法的变种。神经网络中参数量非常巨大，在这种情况下反而这种方法效果比较好。
• 优缺点：可应对复杂势能面；但收敛性依赖于参数 $\alpha$ 和 $\mu$，需要调参（$\mu$ 太小可能还有震荡，$\mu$ 太大可能走不动或也震荡；$\alpha$ 太大也会出现类似梯度下降的问题）。

Conjugate Gradient（CG, 共轭梯度法）

说明：想要提升效率，一定要把之前各步的信息记录下来。CG 的想法是：第 $n+1$ 步的受力，不仅要和第 $n$ 步的更新方向垂直，还想和之前的都垂直。推导后发现，不需要把之前所有步都记录下来做正交化，只需要上一步的力、上一步的更新方向和这一步的力，三者就可以得到新的搜索方向，非常巧妙。

为提升搜索效率，$F_{n+1}$ 不仅应与 $\delta X_n$ 垂直，还应与 $\delta X_{n-1}, \delta X_{n-2} \dots$ 都垂直。

只需要 $F_n, F_{n-1}$ 和 $\Delta X_{n-1}$ 即可得到新的搜索方向：

$$\Delta X_n = \alpha_n (F_n + \gamma_n \Delta X_{n-1}), \quad \gamma_n = \frac{(F_n - F_{n-1}) \cdot F_n}{F_{n-1} \cdot F_{n-1}}$$

当总能 $E_{\text{tot}}$ 为严格的二次型时（$E_{\text{tot}} = \frac{1}{2} x^T B x$），以下关系严格成立：

$$F_{n+1} \cdot \Delta X_i = 0, \quad F_{n+1} \cdot F_i = 0, \quad \Delta X_{n+1} \cdot B \cdot \Delta X_i = 0, \quad (i=1 \dots n)$$

$B$ 为极小值点附近的 Hessian 矩阵：$B = \frac{\partial^2 E_{\text{tot}}}{\partial x^2}$。

说明：CG 本来就是为了解决正定二次型问题推导出来的。在 $N$ 维空间中，它想用 $N$ 次优化达到最低点。之所以叫"共轭"梯度法，是因为相邻两次搜索方向是中间夹着 Hessian 矩阵构成的（共轭正交）。

说明：如果总能是严格的二次型，上述正交关系严格成立。这意味着之前所做的搜索没有白做——上一步沿某个方向搜索过，下一步不会再朝这个方向走，受力总是与之前垂直。对于非严格二次型，虽然不可能恰好几步就到，但 CG 也是鲁棒的算法，对复杂地形也有效。在真实物理世界中，越接近极小值点就越接近二次型；即使在比较远的地方不是二次型，CG 也能帮助你进入到二次型定义的空间中。

$\alpha_n$ 由一维搜索给出（试探 + 外推）：朝搜索方向试探走一步，看新的一步梯度是多少；如果梯度有这个方向的分量，进行外推，再按外推结果去看，直至梯度不含搜索方向上的分量，line search 结束。

• 优缺点：方便稳定，几乎没有参数（$\alpha$ 由搜索给出，不用人为给定），无需二阶导数信息；但需要精确的一维搜索（每一步里头可能不止一步），收敛耗时不是最少。

Direct Inversion of Iterative Subspace（DIIS）

说明：DIIS 这个名字既可以指 mixing 中的方法，也可以用到 relaxation 中，因为思想是一样的。Mixing 中极小化的是残差，而 relaxation 中没有残差，有的是梯度（受力）。但在目标点处，残差为零，梯度也为零，所以可以把之前公式中的残差替换为受力（梯度），推导基本和 Pulay Mixing 一致。

DIIS 是拟牛顿法（Quasi-Newton）的一种。

• 优化对象：受力的模长 $|\vec{F}|^2$，即将 RMM 算法中的 $\langle R | R \rangle$ 替换为 $|\vec{F}|^2$，其余推导类似。

• 线性组合：$x = \sum_i c_i x_i, \quad \vec{F}(x) = \sum_i c_i \vec{F}_i$

• 归一化条件：没有明确的归一化条件（不像 mixing 中电荷总量一定），常用的为 $\sum_i c_i = 1$（为了和之前一致）。

• 更新方法：最小化 $|\vec{F}|^2$，得 $c_i, x^{\text{min}}, \vec{F}^{\text{min}}$。做一步梯度下降：

  $$X_{n+1} = X^{\text{min}} + G F^{\text{min}}$$

  （$G$ 是步长）。

• Broyden 2nd mixing（拟牛顿法的角度理解）：在 relaxation 优化时，线性相关问题有可能出现，所以矩阵求逆前加入对角项，以防止线性相关：

  $$C_{ij} = (w_0^2 I + \bar{A}){ij}^{-1}, \quad \bar{A}{ij} = \langle \Delta F_i | \Delta F_j \rangle, \quad \Delta F_i = F_{i+1} - F_i$$ $$(B_n)^{-1} = G - \sum_{ij} C_{ij} (G \Delta F_j + \delta X_j) \dots$$

  $B_n$ 是 Hessian 的近似值（Hessian 作用到梯度上）：$X_{n+1} = X_n + (B_n)^{-1} F_n$。

• 优缺点：在极小值点附近收敛较快；但未利用能量信息，只最小化受力，受力非单调时可能向错误的方向优化（比如在解离曲线左侧，往力减小的方向优化反而是错误的方向）。因此可能需要先用别的方法预优化到势阱里面，再切换到这个方法。

BFGS

说明：拟牛顿法的意思是，我们不知道 Hessian 矩阵 $B$ 是什么样的，但想通过有限几步搜索给出的信息，大致给出 Hessian 矩阵 $B$ 的一个近似值 $B_n$。所有拟牛顿法都可以写成：通过 $B_n$ 作用到梯度上给出更新。如果是二次型且 $B$ 恰好是精确的 $B$，那直接下一步就得到精确的解。但实际上不是严格的二次型，$B$ 也是近似出来的，之后还要迭代。不同的拟牛顿法差异在于如何更新 Hessian 矩阵。

BFGS 是拟牛顿法的一种：不断更新 Hessian 近似值 $B_n$，并由 $B_n$ 给出更新量：

$$\Delta X_n = - G B_n F_n$$

• Rank-2 Update：

  $$B_{n+1} = B_n + a u u^T + b v v^T$$

  其中：

  $$u = \Delta X_n, \quad v = B_n \Delta X_n$$ $$a = \frac{1}{\Delta F_n \cdot \Delta X_n}, \quad b = \frac{1}{\Delta X_n \cdot B_n \cdot \Delta X_n}$$

说明：BFGS 中的 B 就是 Broyden。"Rank-2 update" 指的是在上一步的 $B$ 中添入两个秩为 1 的矩阵（$uu^T$ 和 $vv^T$），总共添加了一个秩为 2 的修正项。

• 优缺点：收敛速度快，$B_n$ 总是正定的；但远离平衡位置时可能失败。

拟牛顿法综合比较

说明：从 Quantum ESPRESSO 软件手册中可以看到四种拟牛顿法（DIIS, BFGS, EF, RF）的对比。EF 和 RF 是比较新的算法，看起来更快。

常见的几种拟牛顿法包括：EF, BFGS, RF, DIIS。

经验补充（VASP 中的 IBRION 参数）：

• VASP 中调控结构弛豫的参数叫 IBRION，有很多模式可以选择：
  • IBRION = 1：RMM-DIIS（类似电子步中的 RMM-DIIS 思路）
  • IBRION = 2：Conjugate Gradient (CG)
  • IBRION = 3：Damped Molecular Dynamics
  • 还有其他模式（4, 5, 6 等）
• 经验法则：
  • 如果结构本身就在平衡点附近，用 IBRION = 1（RMM-DIIS）是最快的
  • 如果结构偏离平衡点比较远，用 IBRION = 2（CG）比较稳定

总结：算法之间的联系与分类

说明：以下是对所有方法之间区别和联系的梳理。

按 Hessian 矩阵是否已知分类

1. Hessian 已知（如对角化中能量关于波函数的二阶导矩阵）：

   • Davidson Scheme：做正交化，最后在更小的子空间（如 2×2 或 3×3）中求本征值。更广泛的名字叫 Rayleigh-Ritz 方法。
   • RMM-DIIS：转到拟牛顿法策略上，优化残差模长。

2. Hessian 未知（如结构弛豫问题、RMM 中残差模长关于波函数的二阶导）：

   • 拟牛顿法：尝试构造一个 Hessian 出来，当做真正的 Hessian 去优化。配合 line search（不一定要找到真正的极小值点，只需要找到比当前更好的点就行，一般会出现步长系数）。
   • 共轭梯度法（CG）：依赖于每一步都是极小值点，作为递推中的条件，去证明下一步的更新。要求每一步的 line search 要找到这个方向的极小值。

按步长是否人为指定分类

• 步长计算出来：如共轭梯度法中的 $\alpha$ 由 line search 给出
• 步长人为固定：如 linear mixing 中的 $G$

补充：RMM-DIIS 的思路同时对电子步和离子步都适用。电子步指的是 Pulay Mixing（电荷密度 mixing），离子步指的是用 DIIS 做结构弛豫（ionic relaxation），思路类似，所以都叫 RMM-DIIS。