自洽场计算、迭代对角化与 Hellmann-Feynman 力

本文是「第一性原理的微观计算模拟」系列的第 7 篇（共 11 篇），内容整理自同名书稿第 3 章，公式编号与原书一致。文中引用的文献依据原书参考文献清单整理并列于文末，按原书清单顺序从 1 开始编号。\nocite{*}

自洽场计算的实现

3.4.3 节已经给出了平面波-赝势框架下体系的总能表达式。从方程（3.458）不难看出，哈密顿矩阵元中的 Hartree 势 $V_{\mathrm{H}}(\boldsymbol{G})$ 必须通过 $\rho(\boldsymbol{G})$ 构建，但是 $\rho(\boldsymbol{G})$ 正是我们需要求解的物理量。这意味着必须已知 $\rho(\boldsymbol{G})$ 才能求解 $\rho(\boldsymbol{G})$。解决这个矛盾的一般方法是首先给定一个初始猜测，然后通过自洽场（self-consistent field，SCF）计算逐步逼近精确解。本节具体讨论自洽场计算的方法与过程。

自洽过程

图 3.11 为自洽场计算的流程图。

自洽场计算的基本步骤如下。

（1）初始化：首先对电荷密度 $\rho$ 进行一个合理的初始猜测。这通常可以通过将靠近原子核的位置上的电荷密度相叠加得到。初始电荷密度的选择对于自洽场计算的收敛速度和稳定性具有重要意义。

（2）构建势能：根据初始电荷密度 $\rho$ 计算出交换相关势 $V_{\mathrm{xc}}$ 和 Hartree 势 $V_{\mathrm{H}}$。这些势能与赝势一起构成系统的总势能。

（3）求解 Kohn-Sham 方程：使用总势能求解 Kohn-Sham 方程，得到一组 Kohn-Sham 本征值和本征函数，这些本征函数可以用于计算新的电荷密度。

（4）更新电荷密度：根据求得的本征函数计算新的电荷密度 $\rho'$。为了提高收敛稳定性，可以采用混合策略将新的电荷密度与旧的电荷密度相结合，如 $\rho_{\mathrm{new}}=\alpha\rho'+(1-\alpha)\rho$，其中 $\alpha$ 是介于 0 和 1 之间的混合系数。

（5）收敛检查：检查电荷密度、总能量或其他相关物理量的变化是否满足预设的收敛标准。如果满足收敛标准，则自洽场计算完成，可以得到体系的总能量和其他物理量。如果不满足收敛标准，则使用更新后的电荷密度 $\rho_{\mathrm{new}}$，返回步骤（2）继续迭代。

在实际计算的过程中，我们可能需要采用一些策略和算法来提升自洽场计算的收敛速度和稳定性，例如使用 Pulay 混合方法、Kerker 预条件等。

此外，根据具体问题的特性，我们可以选择适合的交换相关泛函、赝势和基组，以提高计算的精度和效率。

在计算过程中，最耗时的部分是矩阵对角化。最直接的算法是利用 LAPACK 中的标准库函数直接对角化哈密顿矩阵，从而得到体系的本征值和本征波函数。但是，这种直接对角化方法需要在计算过程中存储整个 $N$ 阶哈密顿矩阵，因此对内存的需求量很大。此外，直接对角化方法的计算量正比于哈密顿矩阵维数的三次方（即 $N^{3}$）。因此，直接对角化方法并不适合用于处理大型体系（例如原胞内原子数大于 20 的体系）。目前，采用平面波为基组并使用赝势的软件包大都使用所谓的迭代对角化方法，这种方法可以有效地克服直接对角化方法的这两个缺点。

通过自洽场计算，可以获取系统的电子结构信息，这为我们进一步研究材料的各种性质和现象提供了基础。

电荷密度更新

由 3.4.4.1 节可知，自洽场计算的每一步都需要迭代更新电荷密度。最直接的方案是用当前步得到的输出电荷密度 $\rho_i(\boldsymbol{r})$ 作为下一步的输入电荷密度。但是采取这种更新方式会导致自洽场计算不收敛。实际应用中常用的更新方法是将下一步的输入电荷表示为当前步的输入电荷以及输出电荷的线性叠加：

$$ \rho_{i+1}^{\mathrm{in}}=\beta\rho_i^{\mathrm{out}}+(1-\beta)\rho_i^{\mathrm{in}} \tag{3.463} $$

式中：$\beta$ 是一个经验参数。对于一般的体系，取 $\beta=0.3$ 可以保证自洽场计算收敛。但是对于自旋极化体系，$\beta$ 有可能需要取得很小，如 0.05 左右。因为新一轮计算所得的电荷强度仅有一小部分用于更新电荷，所以体系只能非常缓慢地向精确解逼近。在这种情况下，更好的选择是采用此前若干步的输入、输出电荷密度构建最佳的近似解，作为最新一步的输入电荷密度。这种方法称为 Pulay 电荷更新。其详细介绍请参看附录 A.8 节。

虽然平面波基组下电荷密度在动量空间中计算看起来更为直接，但效率更高的做法是通过快速傅里叶逆变换将本征矢 $\{c_{n,\boldsymbol{k}_i}(\boldsymbol{G})\}$ 变换为实空间网格点上的本征函数值 $\{\psi_{n,\boldsymbol{k}_i}(\boldsymbol{r})\}$，然后利用公式

$$ \rho(\boldsymbol r)=\sum_{\boldsymbol k}w_{\boldsymbol k} \sum_n f_{n\boldsymbol k}|\psi_{n\boldsymbol k}(\boldsymbol r)|^2, \qquad\sum_{\boldsymbol k}w_{\boldsymbol k}=1 \tag{3.464} $$

计算实空间各网格点上的电荷密度值。如在接下来的计算过程中需要，则再通过快速傅里叶变换计算 $\rho(\boldsymbol{G})$。

利用共轭梯度法求解广义本征值

在 DFT 理论框架下，对于弱相互作用体系，求解体系基态等价为优化下述泛函：

$$ E=\langle\phi\,|\,\hat{H}\,|\,\phi\rangle-\varepsilon_n\langle\phi\,|\,\hat{\mathrm{S}}\,|\,\phi\rangle \tag{3.465} $$

式中右端第二项来自于本征函数正交性的约束。在特定的函数基下看待此问题，本征函数 $|\,\phi\rangle$ 相当于矢量 $\boldsymbol{x}$，哈密顿算符 $\hat{H}$ 表征为一个矩阵 $\boldsymbol{A}$，$E$ 等于目标函数 $F$，则上述问题等价为优化一个二次函数。即使考虑 Kohn-Sham 方程 $\hat{H}\,|\,\phi\rangle=\varepsilon\hat{\mathrm{S}}\,|\,\phi\rangle$，其广义本征值问题仍然可以等同于函数优化，因此可以利用共轭梯度法求出最接近实际本征值的近似本征值及其相应的本征矢。这种迭代求解 Kohn-Sham 方程的做法有别于直接对角化矩阵以及 Car-Parrinello 动力学方法（统称为直接法）。

由 1.3.2 节的讨论可知，利用共轭梯度法优化目标函数时需要确定最速下降方向及共轭方向。因为 Kohn-Sham 方程本身的特点，还需要考虑正交化处理以及利用预处理技术提升收敛速度。对此分别加以介绍。

1. 最速下降方向

选定某条能带 $m$，根据方程（3.465），优化泛函 $E=\langle\phi_m\,|\,\hat{H}-\varepsilon_m\hat{\mathrm{S}}\,|\,\phi_m\rangle$。由此定义知，第 $i$ 次迭代时相应的残余矢量为

$$ |\,R(\phi_m^{i})\rangle=-(\hat{H}-\varepsilon_m\hat{\mathrm{S}})\,|\,\phi_m^{i}\rangle \tag{3.466} $$

即此处的最速下降方向与 $|\,\phi_m^{i}\rangle$ 正交。因此，此时的拉格朗日乘子可计算如下：

$$ \varepsilon_m^{i}=\frac{\langle\phi_m^{i}\,|\,\hat{H}\,|\,\phi_m^{i}\rangle}{\langle\phi_m^{i}\,|\,\hat{\mathrm{S}}\,|\,\phi_m^{i}\rangle} \tag{3.467} $$

这个值也是第 $m$ 个本征值的当前最佳估计值。

2. 正交化

正交化的要求源自以下矛盾：共轭梯度法是一种无约束的优化方法，但是如果要进行一系列能量本征值的求解，那么要求分属不同本征值的本征矢彼此正交。可以通过对最速下降方向进行正交化处理而将约束优化问题转化为无约束优化问题。因为每条能带的本征矢都与其他的本征矢正交，所以假设已经将指标小于 $m$ 的所有能带优化完毕，那么第 $m$ 条能带的本征值应该是满足与 $m-1$ 个本征矢正交的最小的本征值。因此，第 $i$ 次迭代的最速下降方向应该与 $m-1$ 个本征矢正交。利用格拉姆-施密特正交化方案实现，即

$$ |\zeta_m^i\rangle=|R(\phi_m^i)\rangle-\sum_{n\lt m}|\phi_n\rangle\langle\phi_n|\hat S|R(\phi_m^i)\rangle \tag{3.468} $$

3. 预处理

从理论上讲，在 $N$ 步之内对能带 $m$ 的优化可以结束，但是可以进行一系列操作来提高优化效率。从数学上讲，预处理等于对矩阵 $\boldsymbol{A}$ 进行相似变换，改善其条件数，使得尽量多的本征值简并。这种操作之所以可以提高效率，是因为如果最速下降方向是当前误差（即当前解与精确解之差）乘以一个常数的话，那么沿着最速下降方向移动适当的距离就可以非常精确地到达精确解。

设当前步骤下最优解与精确解之间的差别为 $|\,\delta\phi_m^{i}\rangle$，这个量可以用体系的本征值展开为

$$ |\,\delta\phi_m^{i}\rangle=\sum_{n}c_{n,m}\,|\,\phi_n\rangle \tag{3.469} $$

因此第 $m$ 条能带的精确解可以写为 $|\,\phi_m\rangle=|\,\phi_m^{i}\rangle+|\,\delta\phi_m^{i}\rangle$，将其代入方程（3.466），可得

$$ |\,R(\phi_m^{i})\rangle=-(\hat{H}-\varepsilon_m\hat{\mathrm{S}})\,|\,\phi_m\rangle+(\hat{H}-\varepsilon_m\hat{\mathrm{S}})\,|\,\delta\phi_m\rangle \tag{3.470} $$

在比较接近精确解的时候，式（3.470）右端的第一项可以忽略，仅考虑第二项即可。将式（3.469）代入式（3.470），可得

$$ |\,R(\phi_m^{i})\rangle=(\hat{H}-\varepsilon_m\hat{\mathrm{S}})\,|\,\delta\phi_m\rangle=\sum n(\varepsilon_n-\varepsilon_m^{i})c_{n,m}\,|\,\phi_n\rangle $$

可以看出，如果 $n$ 个能级彼此简并，则最速下降方向是 $|\,\delta\phi_m^{i}\rangle$ 的常数倍。因此，如前所述，沿着最速下降方向移动适当的距离就可以非常精确地到达精确解。这可以大大提升共轭梯度法的收敛速度。而预处理可以通过乘以一个预处理矩阵 $\boldsymbol{K}$ 得以实现，$\boldsymbol{K}$ 取决于在计算中所采用的基函数。

以平面波基为例，对于 $G$ 较高的平面波，动能项为主要项，因此，如果要构造一个简并度比较高的变换，那么令计算最简便的矩阵 $\boldsymbol{K}$ 是一个对角矩阵，对角元是动能的倒数。但是对于 $G$ 较低的平面波，动能项不占优势，因此 $\boldsymbol{K}$ 应该趋近于 1。一般取下面的表达式：

$$ K_{m,n}=\frac{27+18x+12x^{2}+8x^{3}}{27+18x+12x^{2}+8x^{3}+16x^{4}}\delta_{mn} \tag{3.471} $$

式中：$x$ 为平面波动能与 $|\,\phi_m^{i}\rangle$ 动能的比值。

同时考虑最速下降方向的正交性与预处理，可以构造最速下降方向为

$$ |\,\eta_m^{i}\rangle=\boldsymbol{K}\,|\,\zeta_m^{i}\rangle $$

但是乘以矩阵 $\boldsymbol{K}$ 会破坏正交性，因此需要特别对 $|\,\eta_m^{i}\rangle$ 再进行一次格拉姆-施密特正交化：

$$ |\eta_m^{\prime i}\rangle=|\eta_m^i\rangle -\sum_{n\le m}|\phi_n\rangle\langle\phi_n|\hat S|\eta_m^i\rangle \tag{3.472} $$

将 $|\,\eta_m'^{\,i}\rangle$ 作为最速下降方向。注意式（3.472）右端第三项中的 $|\,\phi_n\rangle$ 没有上标，表明第 $n$ 个能带以下的能带均已优化到精确解。

4. 共轭方向

确定最速下降方向之后，可以依照经典的共轭梯度方法构造共轭方向：

$$ |\,\varphi_m^{i}\rangle=|\,\eta_m'^{\,i}\rangle+\gamma_m^{i}\,|\,\varphi_m^{i-1}\rangle,\quad\gamma_m^{i}=\frac{\langle\eta_m'^{\,i}\,|\,\zeta_m^{i}\rangle}{\langle\eta_m'^{\,i-1}\,|\,\zeta_m^{i-1}\rangle} $$

应当注意，共轭方向中的 $\gamma_m^{i}$ 不仅仅有一种表达式。比如 Dyutiman Das 采用了 $\gamma_m^{i}$ 的另外一种形式：

$$ \gamma_m^{i}=\frac{(\langle\eta_m'^{\,i}\,|-\langle\eta_m'^{\,i-1}\,|)\,|\,\eta_m'^{\,i}\rangle}{(\langle\eta_m'^{\,i}\,|-\langle\eta_m'^{\,i-1}\,|)\,|\,\varphi_m^{i}\rangle} \tag{3.473} $$

这些表达式在没有外约束的条件下是彼此等价的，但是因为本征矢正交性的限制，由不同的 $\gamma_m^{i}$ 构造出来的 $|\,\varphi_m^{i}\rangle$ 并不相同，很难说哪一种效率更高，需要在具体的问题中通过测试确定。

当前的共轭方向 $|\,\varphi_m^{i}\rangle$ 还需要与当前第 $m$ 个能带的本征矢正交，并归一化。因此，最终的共轭方向的表达式为

$$ \begin{cases} |\,\varphi_m''^{\,i}\rangle=|\,\varphi_m^{i}\rangle-\langle\phi_m^{i}\,|\,\hat S\,|\,\varphi_m^{i}\rangle\,|\,\phi_m^{i}\rangle\\ |\,\varphi_m'^{\,i}\rangle=\dfrac{|\,\varphi_m''^{\,i}\rangle}{\sqrt{\langle\varphi_m''^{\,i}\,|\,\hat S\,|\,\varphi_m''^{\,i}\rangle}} \end{cases} \tag{3.474} $$

而以 $|\,\varphi_m'^{\,i}\rangle$ 最终的共轭梯度方向作为优化方向。

5. 一维搜索

确定了优化方向后，需要沿优化方向求出目标函数的最优解，而相应的本征矢 $|\,\phi_m^{i+1}\rangle$ 相当于当前本征矢 $|\,\phi_m^{i}\rangle$ 和优化方向 $|\,\varphi_m'^{\,i}\rangle$ 的一个线性叠加。考虑到 $|\,\phi_m^{i+1}\rangle$ 和 $|\,\phi_m^{i}\rangle$ 的模方应该相等，因此可写为 ${|\,\phi_m^{i+1}\rangle=\cos\theta\,|\,\phi_m^{i}\rangle+\sin\theta\,|\,\varphi_m'^{\,i}\rangle}$，优化参数 $\theta$ 即可。前面说过，对于二次正定的函数，优化步长有解析的形式。假设采用经验赝势方法，则可写出 $\theta_{\min}$ 的解析式为

$$ \tan(2\theta)=\frac{2\langle\varphi_m'^{\,i}\,|\,\hat{H}\,|\,\phi_m^{i}\rangle}{\langle\phi_m^{i}\,|\,\hat{H}\,|\,\phi_m^{i}\rangle-\langle\varphi_m'^{\,i}\,|\,\hat{H}\,|\,\varphi_m'^{\,i}\rangle} \tag{3.475} $$

如果采用严格的第一性原理计算，那么 $\theta_{\min}$ 虽然仍有解析形式，但是需要考虑实空间内交换关联能及 Hartree 项的积分。

另外一种方法则需要算出 $\theta=0$ 时的函数值及一阶导数值，以及取另一个 $\theta$ 值（通常取 $\pi/300$）时的函数值，具体的步骤如下。

首先将能量 $E$ 写为关于 $\theta$ 的三角级数：

$$ E(\theta)=E_0+\sum_{n}[A_n\cos(2n\theta)+B_n\sin(2n\theta)] \tag{3.476} $$

Payne 和 Joannopoulods 指出，这个级数中 $n\gt 1$ 的项均可以省略。因此式（3.476）简化为 $E(\theta)=E_0+A_1\cos2\theta+B_1\sin2\theta$。因此，如果要确定 $\theta_{\min}$，我们需要先求出 $E_0$、$A_1$ 及 $B_1$ 三个参数的值。可以通过下述三个方程求解：

$$ \begin{gathered} E_0=\frac{E\left(\dfrac{\pi}{300}\right)-\dfrac{1}{2}\left.\dfrac{\partial E}{\partial\theta}\right|_{\theta=0}-E(0)\cos\dfrac{2\pi}{300}}{1-\cos\dfrac{2\pi}{300}}\\ A_1=\frac{E(0)-E\left(\dfrac{\pi}{300}\right)+\dfrac{1}{2}\left.\dfrac{\partial E}{\partial\theta}\right|_{\theta=0}}{1-\cos\dfrac{2\pi}{300}}\\ B_1=\frac{1}{2}\left.\frac{\partial E}{\partial\theta}\right|_{\theta=0} \end{gathered} $$

式中

$$ \left.\frac{\partial E}{\partial\theta}\right|_{\theta=0}=\langle\varphi_m'^{\,i}\,|\,\hat{H}\,|\,\phi_m^{i}\rangle+\langle\phi_m^{i}\,|\,\hat{H}\,|\,\varphi_m'^{\,i}\rangle \tag{3.477} $$

则极值点 $\theta_s=\dfrac{1}{2}\arctan\left(\dfrac{B_1}{A_1}\right)$，在区间 $\left[0,\dfrac{\pi}{2}\right]$ 中的 $\theta_s$ 即为所求的 $\theta_0$。

以上所介绍的这两种方法的计算量相差无几。

综上所述，利用共轭梯度法求解 Kohn-Sham 方程的本征值具体步骤如下：

（1）预设收敛判据 $\tau$ 及 $\lambda$，初始化 $N$ 个本征矢 $|\,\phi\rangle$（如利用随机数作为系数），设 $j=0$，$|\,\phi^{j}\rangle=|\,\phi\rangle$，以原子电荷分布的叠加作为初始电荷密度 $\rho_j^{\mathrm{in}}$；

（2）选择最低的能带 $m$，设 $i=0$，根据式（3.467）求出在上述 $\rho_j^{\mathrm{in}}$ 下的期待值 $\varepsilon_m^{i}$，并由式（3.466）计算残余矢量 $|\,R(\phi_m^{i})\rangle$；

（3）依次按照式（3.468）至式（3.472）对 $|\,R(\phi_m^{i})\rangle$ 进行操作，并通过式（3.473）和式（3.474）构造归一化的共轭梯度方向 $|\,\varphi_m'^{\,i}\rangle$；

（4）进行一维搜索，利用方程（3.475）或者上述两种方法计算优化的本征矢 $|\,\phi_m^{i+1}\rangle$，计算出 $\varepsilon_m^{i+1}$ 和 $\|\,|\,R(\phi_m^{i+1})\rangle\,\|$，若 $\|\,|\,R(\phi_m^{i+1})\rangle\,\|\leqslant\tau$，优化结束，设 $m=m+1$，移到下一条能带，否则设 $i=i+1$，回到步骤（2）；

（5）重复上述过程，直至算法收敛或者 $m\geqslant N$，计算总能 $E$ 和能量变化值 $\Delta E^{j}$，若 $\Delta E^{j}\leqslant\lambda$，全部计算结束，转到步骤（6），否则设 $j=j+1$，利用 $\phi^{j}$ 更新哈密顿矩阵以及电荷密度 $\rho_j^{\mathrm{in}}$，回到步骤（1）；

（6）计算并保存电荷密度、本征矢、总能等各种信息。

对最速下降方向以及共轭方向进行约束的方法不只上面一种。如果正交化操作不仅仅是对这些能带进行，而是将 $n$ 的取值遍历所有能带指标，则利用类似的算法可以得到同样的本征能级，但是不能保证所得的本征矢是正确的。实际上，这些本征矢是 Kohn-Sham 本征矢的线性叠加。因此，对金属而言，需要在共轭梯度法求解过程结束之后再进行子空间转动（subspace rotation）这一步骤，即以占据数不为零的所有能带对应的本征矢张开子空间，计算哈密顿矩阵 $\boldsymbol{H}$ 以及交叠矩阵 $\boldsymbol{S}$，再进行矩阵的直接对角化。用所得的本征矢 $\{|\,B\rangle\}$ 乘以此前得到的本征矢 $\{|\,\phi^{j}\rangle\}$，即可得到最后的结果。

值得注意的是，在优化全部 $N$ 个本征矢的过程中，电荷密度 $\rho_j^{\mathrm{in}}$ 是保持不变的，只有在所有 $N$ 个本征矢优化结束之后才更新 $\rho_j^{\mathrm{in}}$。因此，上述方法是迭代算法。另一种可能的方法是直接使用共轭梯度法对体系进行求解，这个过程与上述算法大体一致。该方法与上述算法主要的区别在于，对第 $m$ 条能带优化结束之后，需要先更新 $\rho^{\mathrm{in}}$，然后移至下一条能带。然而，在实际计算中，特别是在使用平面波基方法的情况下，最初的几轮优化更新 $\rho^{\mathrm{in}}$ 可能会导致系统严重偏离基态（这是因为本征矢的初始化使用了随机数），因此这种方法在实践中存在可行性问题。一种解决办法是，在所有能带的优化过程中保持初始的 $\rho^{\mathrm{in}}$ 不变，从第二步开始随时更新 $\rho^{\mathrm{in}}$。为了提高效率，所设置的共轭梯度法的收敛判据 $\tau$ 和 $\lambda$ 的值会随着优化的进行而逐渐减小，而不是固定的。

迭代对角化方法

哈密顿矩阵 $\boldsymbol{H}$（或者重叠矩阵 $\boldsymbol{S}$）的维数通常很大。直接对角化之后，总共有 $N$ 个本征值和本征矢，而被电子占据的能带只占其中一小部分。因此，为了避免对高维矩阵的直接对角化，我们可以专注于最低的 $n$ 个能带的精度。一种合理的方法是针对给定的能带，首先给出近似的本征值和本征矢，将其代入本征值方程（或广义本征值方程），以得到改进的结果。使上述过程迭代进行，直至结果的精度达到要求。这就是迭代对角化方法的基本思想。常见的迭代对角化方法包括 Lanczos 方法、Davidson 方法和残差矢量最小化方法（RMM-DIIS）等。在本节中，我们将首先介绍迭代对角化方法的基本理论，然后详细讨论 RMM-DIIS 方法，最后探讨实际应用中可以提高计算效率的因素。迭代对角化方法的核心思想是逐步改进对目标能带的本征值和本征矢的近似，以在有限的迭代步骤内获得所需精度。这些方法通常利用初始近似解和目标矩阵的一些性质来构建一个较小的子空间，在这个子空间中进行对角化，然后通过更新近似解和子空间来进行迭代。

总之，迭代对角化方法为求解大规模矩阵本征问题提供了一种高效的解决方案。通过选择合适的方法和技巧，可以在有限的迭代步骤内获得满足精度要求的结果，从而提高计算效率和可靠性。

基本理论

绝大多数的迭代对角化方法都会定义或者构建三组 $N$ 维矢量。第一组 $\{|\,\varphi_i\rangle\}$ 有 $N$ 个元素，是希尔伯特空间的基函数（或称坐标轴），如平面波基等，哈密顿矩阵 $\boldsymbol{H}$ 和重叠矩阵 $\boldsymbol{S}$ 都可以由此得出；第二组 $\{|\,x_i\rangle\}$ 也有 $N$ 个元素，是一组完备基，可以张开整个希尔伯特空间中的任意矢量；第三组 $\{|\,b_i\rangle\}$ 只有 $N_0$（$N_0\ll N$）个元素，是 $N_0$ 维子空间基函数，因为只要求 $\{|\,b_i\rangle\}$ 张开可以包含 $n$ 条最低能带的 $N_0$ 维子空间，所以每个 $|\,b_i\rangle$ 中只需要前 $N_0$ 个元素准确。启动迭代对角化方法时，首先选定一个 $N_0$ 阶哈密顿矩阵 $\boldsymbol{H}_0$（$\boldsymbol{H}$ 的一部分），然后利用直接对角化方法如 Cholesky 分解法等求解 $\boldsymbol{H}_0$ 的本征值和本征矢，若精度不够，则利用所得结果更新 $\{|\,b_i\rangle\}$，并在 $\{|\,b_i\rangle\}$ 张开的子空间中重新求解下列方程：

$$ \varXi\,|\,c\rangle=\varepsilon\varOmega\,|\,c\rangle \tag{3.478} $$

有

$$ \begin{cases} \varXi_{ij}=\langle b_i|\hat H|b_j\rangle,\\ \varOmega_{ij}=\langle b_i|\hat S|b_j\rangle. \end{cases} \tag{3.479} $$

$\varepsilon_k$ 为当前哈密顿矩阵 $\boldsymbol{H}$ 的第 $k$ 个本征值的近似值，相应的本征矢 $|\,a_k\rangle$ 为

$$ |a_k\rangle=\sum_i c_i^{(k)}|b_i\rangle,\qquad \varXi\boldsymbol c^{(k)}=\varepsilon_k\varOmega\boldsymbol c^{(k)} \tag{3.480} $$

由计算结果构建 $|\,x_i\rangle$ 以及 $|\,b_i\rangle$ 的方法不同，导致了迭代对角化方法的不同。但是这些方法均遵循一个原则，即更新后的本征矢应当尽量靠近体系的精确解。为了完成这个任务，首先定义残余矢量为

$$ |R(A^{\mathrm c},E^{\mathrm c})\rangle =(\hat H-E^{\mathrm c}\hat S)|A^{\mathrm c}\rangle \tag{3.481} $$

式中：$A^{\mathrm{c}}$ 和 $E^{\mathrm{c}}$ 分别为当前本征矢和本征值的近似估算值。$|\,R\rangle$ 的模 $(\langle R\,|\,R\rangle/\langle A^{\mathrm{c}}\,|\,\hat{\mathrm{S}}\,|\,A^{\mathrm{c}}\rangle)^{1/2}$ 反映了当前结果至精确值的“距离”。而 $E^{\mathrm{c}}$ 的计算非常直接：

$$ E^{\mathrm{c}}=\frac{\langle A^{\mathrm{c}}\,|\,\hat{H}\,|\,A^{\mathrm{c}}\rangle}{\langle A^{\mathrm{c}}\,|\,\hat{\mathrm{S}}\,|\,A^{\mathrm{c}}\rangle} \tag{3.482} $$

下一轮迭代，相当于在当前的本征矢近似值 $|\,A^{\mathrm{c}}\rangle$ 上叠加一个 $|\,\delta A\rangle$。最理想的结果是更新后的矢量 $|\,A^{\mathrm{c}}\rangle+|\,\delta A\rangle$ 就是精确的本征矢，此时残余矢量为零，且

$$ |R(A^{\mathrm c}+\delta A,E^{\mathrm c})\rangle =|R(A^{\mathrm c},E^{\mathrm c})\rangle +(\hat H-E^{\mathrm c}\hat S)|\delta A\rangle=0 \tag{3.483} $$

由此可得最理想的 $|\,\delta A\rangle$ 应为

$$ |\,\delta A\rangle=-(\hat{H}-E^{\mathrm{c}}\hat{\mathrm{S}})^{-1}\,|\,R(|\,A^{\mathrm{c}}\rangle,E^{\mathrm{c}})\rangle \tag{3.484} $$

但是因为需要对一个 $N$ 阶矩阵求逆，所以通过式（3.484）计算 $|\,\delta A\rangle$ 并不现实。若利用完备基组 $|\,x_j\rangle$ 将 $|\,\delta A\rangle$ 展开，代入式（3.483），并与 $\langle x_i\,|$ 做内积，则有

$$ \langle x_i\,|\,R\rangle+\sum_{j}\langle x_i\,|\,(\hat{H}-E^{\mathrm{c}}\hat{\mathrm{S}})\,|\,x_j\rangle\langle x_j\,|\,\delta A\rangle=0 \tag{3.485} $$

式（3.485）并无法减小直接计算 $|\,\delta A\rangle$ 所需的计算量。为了解决这个问题，一般采用所谓对角近似，即只保留式（3.485）中 $i=j$ 的项，因此有

$$ |\,\delta A\rangle=-\sum_{i}{}'\frac{\langle x_i\,|\,R\rangle\,|\,x_i\rangle}{\langle x_i\,|\,\hat{H}-E^{\mathrm{c}}\hat{\mathrm{S}}\,|\,x_i\rangle} \tag{3.486} $$

其中求和符号上的“$'$”表示剔除任何分母小于某个阈值 $\delta$ 的项。这个定义确保了当第 $k$ 个本征矢的残余矢量 $|\,R_k\rangle$ 为 $\boldsymbol{0}$ 时，相应的 $|\,\delta A_k\rangle$ 也为 $\boldsymbol{0}$。

利用式（3.486）更新本征矢的近似值 $|\,A\rangle$ 后，再将其代入式（3.482）求得更新后的本征值近似值 $E^{\mathrm{new}}$。若相应的 $R^{\mathrm{new}}$ 仍然比较大，则重复上述过程。这就构成了大多数迭代对角化方法的基本算法。此外，应当注意，到目前为止，算法对本征值的优化是串行的，即每次只优化一个本征值，结束之后再优化下一个。

RMM-DIIS 方法

RMM-DIIS 方法是一种在量子化学和凝聚态物理计算中广泛应用的迭代对角化方法。其核心思想是在每一轮迭代过程中，找出一个最佳线性组合以便最小化残差矢量的范数。RMM-DIIS 方法的优势在于，其能够加快收敛速度，尤其对于那些难以收敛的问题，同时它保持了较低的计算成本。在实际应用中，可以通过采取适当的预处理技术、控制子空间的尺寸，以及采用选择性的收敛标准等方式，提升 RMM-DIIS 方法的计算效率。

这个方法是由 Wood 和 Zunger 在 1984 年正式提出的\cite{wood1985method}。在此之前，Pulay 为了提高自洽场计算中电荷更新的精确度，提出过一个与 RMM-DIIS 方法非常类似的算法，这类算法称为迭代子空间直接求逆（direct inversion in the iterative subspace，DIIS）法，或者 Pulay 电荷更新法\cite{pulay1980convergence}，在 9.8 节中将简要介绍该方法。在本节中，我们主要介绍 RMM-DIIS 方法。

对于第 $j$ 个本征矢和本征值，选取 $\{|\,x_i\rangle\}$ 和 $\{|\,b_i\rangle\}$ 分别为

$$ \{|\,x_i\rangle\}=\{|\,a_j^{0}\rangle,j=1,2,\cdots,N_0\}+\{|\,e_j\rangle,j=N_0+1,\cdots,N\} \tag{3.487} $$

$$ \{|\,b_i\rangle\}[p=0/1/2/\cdots]=[\,|\,\delta A_j^{(0)}\rangle/|\,\delta A_j^{(1)}\rangle/|\,\delta A_j^{(2)}\rangle/\cdots] \tag{3.488} $$

式中：$|\,a_j^{0}\rangle$ 表示哈密顿矩阵 $\boldsymbol{H}_0$ 的一套本征矢；$|\,\delta A_j^{(0)}\rangle$ 是一个 $N$ 维矢量，前 $N_0$ 个元素由 $|\,a_j^{0}\rangle$ 给出，之后的 $N-N_0$ 个元素为 0；$p$ 代表迭代的次数，“/”表示每一轮会增加的子空间基函数 $|\,b_i\rangle$。例如，首次迭代，$|\,b_i\rangle=|\,\delta A_j^{(0)}\rangle$，而下一次迭代，$\{|\,b_i\rangle\}$ 有两个矢量——$|\,\delta A_j^{(0)}\rangle$ 和 $|\,\delta A_j^{(1)}\rangle$，依次类推。可见，RMM-DIIS 方法包含了全部迭代信息。其中用来优化本征值和本征矢的子空间由迭代产生的 $|\,\delta A_j^{(p)}\rangle$ 张开，又称为迭代子空间。将式（3.487）代入式（3.486），可得 $\delta A$ 的各个分量：

$$ \langle e_k\,|\,\delta A\rangle=-\sum_{i\in(1,N_0)}{}'\frac{{\langle a_i^{0}\,|\,R(E^{\mathrm{old}})\rangle\langle e_k\,|\,a_i^{0}\rangle}}{(\lambda_i^{0}-E^{\mathrm{old}})\langle a_i^{0}\,|\,\hat{\mathrm{S}}\,|\,a_i^{0}\rangle}-\sum_{i\in(N_0+1,N)}{}'\frac{\langle e_k\,|\,R(E^{\mathrm{old}})\rangle\delta_{ik}}{(H_{ii}-E^{\mathrm{old}}S_{ii})} \tag{3.489} $$

式中：$\lambda_i^{0}$ 是矩阵 $\boldsymbol{H}_0$ 的第 $i$ 个本征值。设现在开始进行第 $m$ 次迭代，已经有了 $m$ 个迭代子空间的基函数 $\{\delta A^{(0)},\delta A^{(1)},\cdots,\delta A^{(m-1)}\}$，同时有 $E^{\mathrm{old}}=E^{(m-1)}$，则根据式（3.489）得到 $|\,\delta A^{(m)}\rangle$。根据这些信息构建第 $m$ 轮中本征矢 $|\,A^{\mathrm{new}}\rangle$ 的最佳估计值，将其在迭代子空间中展开，有

$$ |\,A^{\mathrm{new}}\rangle=\sum_{i=0}^{m}\alpha_i\,|\,\delta A^{(i)}\rangle \tag{3.490} $$

$|\,A^{\mathrm{new}}\rangle$ 最佳意味着相应的残余矢量 $|\,R(A^{\mathrm{new}},E^{\mathrm{old}})\rangle$ 的模方 $\rho^{2}$ 最小。因此 $\rho^{2}$ 对于任意系数 $\alpha_i$ 的偏导都应为 0，即

$$ \frac{\partial\rho^{2}}{\partial\alpha_k^{*}}=\frac{\partial}{\partial\alpha_k^{*}}\frac{\displaystyle\sum_{r,s=0}^{m}\alpha_r^{*}\alpha_s\langle\delta A^{(r)}\,|\,\hat{H}-E^{\mathrm{old}}\hat{\mathrm{S}}\,|\,(\hat{H}-E^{\mathrm{old}}\hat{\mathrm{S}})\delta A^{(s)}\rangle}{\displaystyle\sum_{r,s=0}^{m}\alpha_r^{*}\alpha_s\langle\delta A^{(r)}\,|\,\hat{\mathrm{S}}\,|\,\delta A^{(s)}\rangle}=0 \tag{3.491} $$

$k$ 从 0 遍历到 $m$，则形如式（3.491）的 $m+1$ 个方程联立，求解 $\rho^{2}$ 的极小值等价于求解下列广义本征值方程：

$$ \boldsymbol{P}\,|\,\alpha\rangle=\rho^{2}\boldsymbol{Q}\,|\,\alpha\rangle \tag{3.492} $$

矩阵 $\boldsymbol{P}$ 和 $\boldsymbol{Q}$ 的矩阵元如下：

$$ \begin{cases} P_{rs}=\langle(\hat H-E^{\mathrm{old}}\hat S)\delta A^{(r)}\,|\, (\hat H-E^{\mathrm{old}}\hat S)\delta A^{(s)}\rangle,\\ Q_{rs}=\langle\delta A^{(r)}|\hat S|\delta A^{(s)}\rangle. \end{cases} \tag{3.493} $$

因为 $\boldsymbol{P}$ 与 $\boldsymbol{Q}$ 均为 $m+1$ 阶方阵，所以可以利用直接对角化求得最小的本征值，将相应的本征矢 $|\,\alpha\rangle$ 代入式（3.490），就得到最优的 $|\,A^{\mathrm{new}}\rangle$。

得到 $|\,A^{\mathrm{new}}\rangle$ 之后可以更新对本征值的近似

$$ E^{\mathrm{new}}=\frac{\langle A^{\mathrm{new}}\,|\,\hat{H}\,|\,A^{\mathrm{new}}\rangle}{\langle A^{\mathrm{new}}\,|\,\hat{\mathrm{S}}\,|\,A^{\mathrm{new}}\rangle} \tag{3.494} $$

以及残余矢量

$$ |R^{\mathrm{new}}\rangle =\frac{(\hat H-E^{\mathrm{new}}\hat S)|A^{\mathrm{new}}\rangle} {\sqrt{\langle A^{\mathrm{new}}|\hat S|A^{\mathrm{new}}\rangle}} \tag{3.495} $$

若 $\|\,|\,R^{\mathrm{new}}\rangle\,\|^{2}$ 大于收敛判据，则再重复上述过程，直至收敛为止。

Wood-Zunger 提出的上述算法在处理大型矩阵时会遇到收敛困难的问题，而且计算效率也比较低，所以当前流行的 RMM-DIIS 方法更接近于 Pulay 的方法\cite{kresse1996efficiency,kresse1996efficient}。二者的主要区别在于在 RMM-DIIS 方法中，首次迭代时，子空间基函数 $|\,b_j\rangle=|\,\delta A_j^{(0)}\rangle$ 不再要求预先求解 $\boldsymbol{H}_0$ 以得到一个初始的本征矢，而是利用共轭梯度或者最速下降法非自洽地求解整个哈密顿矩阵 $\boldsymbol{H}$，以得到的第 $j$ 个本征矢 $|\,A_j^{0}\rangle$ 作为 $|\,b_j\rangle$，相应的本征值记为 $E^{\mathrm{old}}$。并且在首次迭代中，不再用方程（3.489）构建当前步的叠加矢量 $|\,\delta A\rangle$，而是利用 $|\,A_j^{0}\rangle$ 对应的残余矢量 $|\,R(|\,A_j^{0}\rangle,E^{\mathrm{old}})\rangle$ 构建新的基函数，记为 $|\,A_j^{1}\rangle$，即

$$ |\,A_j^{1}\rangle=|\,A_j^{0}\rangle+\lambda\boldsymbol{K}\,|\,R(|\,A_j^{0}\rangle,E^{\mathrm{old}})\rangle \tag{3.496} $$

式中：$\boldsymbol{K}$ 是预处理矩阵，矩阵元由方程（3.471）给出；$\lambda$ 是步长，通常取值范围为 $[0.3,1]$。由式（3.481）得到 $|\,A_j^{1}\rangle$ 相应的残余矢量 $|\,R_j^{1}\rangle=|\,R_j^{1}(|\,A_j^{1},E^{\mathrm{old}})\rangle$。再根据式（3.490）写出本次迭代中本征矢的最佳估计值：

$$ |\,A_j^{M,\mathrm{new}}\rangle=\sum_{i=0}^{M}\alpha_i\,|\,A_j^{i}\rangle,\quad M=1 \tag{3.497} $$

设残余矢量是一个线性算符，则有

$$ |\,R_j^{M,\mathrm{new}}\rangle=\sum_{i=0}^{M}\alpha_i\,|\,R_j^{i}\rangle \tag{3.498} $$

对式（3.498）求 $\|\,|\,R_j^{\mathrm{new}}\rangle\,\|^{2}$ 的极小值，同样可以得到关于 $\{\alpha_i\}$ 的广义本征值方程（3.492），其中

$$ \begin{cases} P_{rs}=\langle R_j^{r}\,|\,R_j^{s}\rangle\\ Q_{rs}=\langle A_j^{r}\,|\,\hat{\mathrm{S}}\,|\,A_j^{s}\rangle \end{cases} \tag{3.499} $$

将求得的 $\{\alpha_i\}$ 代入式（3.497），得到本次迭代下本征矢最优近似解 $|\,A_j^{M,\mathrm{new}}\rangle$，再由式（3.494）和式（3.495）得到 $E_j^{\mathrm{new}}$ 和 $|\,R_j^{M,\mathrm{new}}\rangle$，如果算法未收敛，则设 $E_j^{\mathrm{old}}=E_j^{\mathrm{new}}$，$M=M+1$，并且增加一个新的迭代子空间基函数

$$ |\,A_j^{M+1}\rangle=|\,A_j^{M,\mathrm{new}}\rangle+\lambda\boldsymbol{K}\,|\,R_j^{M,\mathrm{new}}\rangle \tag{3.500} $$

再重复上述过程，直至算法收敛为止。

计算效率优化

从前面的讨论可以看到，迭代对角化方法需要计算 $\boldsymbol{H}\psi$（$N$ 维列向量）。因此，采用迭代对角化方法可以大大地降低计算时对内存空间的要求。而且，快速傅里叶变换允许在实空间和倒空间之间相互切换，对于 $\boldsymbol{H}\psi$ 每一项都选取最高效的方法进行计算。哈密顿算符 $\hat{H}$ 已经在 3.4.3.1 节中由式（3.380）给出。为讨论方便，这里将赝势分解为局域赝势与 KB 非局域赝势两部分，然后重新写出 $\hat{H}$：

$$ \hat{H}=-\frac{\boldsymbol{\nabla}^{2}}{2}+\hat{V}_{\mathrm{H}}+\hat{V}_{\mathrm{ps}}^{\mathrm{loc}}+\hat{\mu}_{\mathrm{xc}}+\hat{V}_{\mathrm{ps}}^{\mathrm{nl}} \tag{3.501} $$

利用平面波将 $\hat{H}$ 和 $\psi$ 分别展开为矩阵 $\boldsymbol{H}$ 和矢量 $c_{n,\boldsymbol{k}_i,\boldsymbol{G}}$ 之后，可以看到哈密顿矩阵的各个组成部分可以分别和矢量相乘。动能项 $\hat{T}_{\mathrm{e}}$ 非常简单，因为在倒空间中只有对角项。而势能算符中，$\hat{V}_{\mathrm{H}}$、$\hat{V}_{\mathrm{ps}}^{\mathrm{loc}}$ 和 $\hat{V}_{\mathrm{ps},l}^{\mathrm{nl}}$ 分别由式（3.384）、式（3.393）和式（3.397）直接在倒空间里给出。而 $\hat{\mu}_{\mathrm{xc}}$ 可以在实空间的格点上计算，再经由快速傅里叶变换得到在倒空间中的表达式。具体写出各项 $\boldsymbol{G}$ 分量的矩阵表达式如下：

$$ \boldsymbol T_{\mathrm e}\,c_{n,\boldsymbol k_i}(\boldsymbol G)=\frac12|\boldsymbol k_i+\boldsymbol G|^2c_{n,\boldsymbol k_i}(\boldsymbol G) \tag{3.502} $$

$$ \boldsymbol{V}_{\mathrm{H}}\boldsymbol{c}_{n,\boldsymbol{k}_i}(\boldsymbol{G})=\sum_{\boldsymbol{G}'}V_{\mathrm{H}}(\boldsymbol{G}-\boldsymbol{G}')c_{n,\boldsymbol{k}_i}(\boldsymbol{G}') \tag{3.503} $$

$$ \boldsymbol{V}_{\mathrm{ps}}^{\mathrm{loc}}\boldsymbol{c}_{n,\boldsymbol{k}_i}(\boldsymbol{G})=\sum_{\boldsymbol{G}'}V_{\mathrm{ps}}^{\mathrm{loc}}(\boldsymbol{G}-\boldsymbol{G}')c_{n,\boldsymbol{k}_i}(\boldsymbol{G}') \tag{3.504} $$

$$ \hat{\mu}_{\mathrm{xc}}\boldsymbol{c}_{n,\boldsymbol{k}_i}(\boldsymbol{G})=\sum_{\boldsymbol{G}'}\mu_{\mathrm{xc}}(\boldsymbol{G}-\boldsymbol{G}')c_{n,\boldsymbol{k}_i}(\boldsymbol{G}') \tag{3.505} $$

$$ \boldsymbol{V}_{\mathrm{ps}}^{\mathrm{nl}}\boldsymbol{c}_{n,\boldsymbol{k}_i}(\boldsymbol{G})=\sum_{l=0}^{l_{\max}}\sum_{m=-l}^{l}\sum_{s=1}^{P_s}\sum_{I=1}^{N_s}\beta_{lm}^{s}\mathrm{e}^{-\mathrm{i}\boldsymbol{G}\cdot\boldsymbol{R}_{I,s}}f_{lm}^{s,*}(\boldsymbol{k}_i+\boldsymbol{G})F_{I,n}^{lm,s}(\boldsymbol{k}_i+\boldsymbol{G}) \tag{3.506} $$

其中式（3.506）用到了式（3.397）至式（3.399）及式（3.462）。至此，我们已经得到了 $\boldsymbol{H}\psi$ 各元素的计算公式。仔细考察式（3.503）至式（3.505），可知当 $\boldsymbol{G}=\boldsymbol{G}'$ 时需要考虑限制条件 $V_{\mathrm{H}}(0)=V_{\mathrm{ps}}^{\mathrm{loc}}(0)=0$。这在一定程度上增加了计算的复杂性，而且对 $\boldsymbol{G}'$ 的求和也显得比较繁杂。有一些软件包，如 CPMD¹ 等，就采用了另一种办法计算这三个方程，即首先在实空间内计算 $\hat{V}\psi_{n,\boldsymbol{k}_i}$，然后利用快速傅里叶变换直接得到相应的 $\boldsymbol{G}$ 分量\cite{kohanoff2006electronic}。为了避免 $V_{\mathrm{H}}(0)$ 与 $V_{\mathrm{ps}}^{\mathrm{loc}}(0)$ 的发散，该方法利用了 3.4.3.3 节中引入的附加电荷分布 $\rho_{\mathrm{aux}}(\boldsymbol{r})$，并定义 $V_{\mathrm{es}}^{\mathrm{loc}}$：

$$ V_{\mathrm{es}}^{\mathrm{loc}}(\boldsymbol{G})=\frac{4\pi}{|\,\boldsymbol{G}\,|^{2}}[\rho(\boldsymbol{G})+\rho_{\mathrm{aux}}(\boldsymbol{G})]+\sum_{s}S^{s}(\boldsymbol{G})\left[V_{\mathrm{ps}}^{\mathrm{loc},s}(\boldsymbol{G})+\frac{4\pi Z^{s}}{|\,\boldsymbol{G}\,|^{2}\varOmega_{\mathrm{cell}}}\mathrm{e}^{-|\boldsymbol{G}|^{2}\sigma^{2}/4}\right] \tag{3.507} $$

不难证明，$V_{\mathrm{es}}^{\mathrm{loc}}(0)=0$，所以不需考虑发散项。利用快速傅里叶逆变换将 $V_{\mathrm{es}}^{\mathrm{loc}}(\boldsymbol{G})$ 与实空间格点上的势函数值 $V_{\mathrm{es}}^{\mathrm{loc}}(\boldsymbol{r})$ 相关联，则有

$$ \hat{V}^{\mathrm{loc}}\boldsymbol{c}_{n,\boldsymbol{k}_i}(\boldsymbol{G})=\frac{1}{\varOmega_{\mathrm{cell}}}\int_{\varOmega_{\mathrm{cell}}}[V_{\mathrm{es}}^{\mathrm{loc}}(\boldsymbol{r})+\mu_{\mathrm{xc}}(\boldsymbol{r})]\psi_{n,\boldsymbol{k}_i}(\boldsymbol{r})\mathrm{e}^{-\mathrm{i}\boldsymbol{G}\cdot\boldsymbol{r}}\,\mathrm{d}\boldsymbol{r} \tag{3.508} $$

图 3.12 展示了迭代对角化方法中 $\boldsymbol{H}\psi$ 的构建过程。可以看到，计算过程同时用到了实空间和倒空间中的积分。两个空间之间通过快速傅里叶变换及其逆变换联系了起来。

实际上，通过对 $\boldsymbol{H}\psi$ 的不同解读可以给出另一种对求解 Kohn-Sham 方程方法的不同理解。例如，将其视为波函数在由基函数张开的希尔伯特空间中受到的广义力，则可以利用 Car-Parrinello 动力学方法（CPMD）求解本征态；将其视为线性空间内对矢量的形变操作，则可利用 RMM-DIIS 方法求解；而将其视为二次函数的梯度，则可利用最优化方法（如共轭梯度法等）进行求解。

Hellmann-Feynman 力

第 $\alpha$ 个原子受力 $\boldsymbol{F}_\alpha$ 的普遍表达式是

$$ \boldsymbol{F}_\alpha=-\frac{\partial E_{\mathrm{tot}}}{\partial\boldsymbol{R}_\alpha} \tag{3.509} $$

式中：$\boldsymbol{R}_\alpha$ 为原子 $i$ 的坐标。从原则上讲，可以用有限差分的方法近似求解式（3.509），但是这种方法的效率和精度都比较差，因此在实际工作中都是采用解析表达式求解 $\boldsymbol{F}_\alpha$，而其理论基础最早由 Hellmann 与 Feynman 提出\cite{hellmann1937einfuhrung,feynman1939forces}。他们指出，在基组完备以及本征波函数严格正确的条件下，$\boldsymbol{F}_\alpha^{\mathrm{HF}}$ 即为该原子在体系中受到的静电力。证明过程如下：

$$ \boldsymbol F_\alpha=-\frac{\partial E_{\mathrm{tot}}}{\partial\boldsymbol R_\alpha} =\boldsymbol F_\alpha^{\mathrm{el}}+\boldsymbol F_\alpha^{\mathrm{ion}},\qquad \boldsymbol F_\alpha^{\mathrm{ion}}=-\frac{\partial E_{\mathrm{II}}}{\partial\boldsymbol R_\alpha} \tag{3.510} $$

这里 $|\Psi\rangle$ 为归一化的电子多体本征态；原子受力分为电子贡献与离子间贡献。对电子能量求导有

$$ \begin{aligned} \boldsymbol F_\alpha^{\mathrm{el}} &=-\frac{\partial}{\partial\boldsymbol R_\alpha} \langle\Psi|\hat H_{\mathrm{el}}|\Psi\rangle\\ &=-\bigl[\langle\partial_\alpha\Psi|\hat H_{\mathrm{el}}|\Psi\rangle +\langle\Psi|\partial_\alpha\hat H_{\mathrm{el}}|\Psi\rangle +\langle\Psi|\hat H_{\mathrm{el}}|\partial_\alpha\Psi\rangle\bigr]\\ &=-\langle\Psi|\partial_\alpha\hat H_{\mathrm{el}}|\Psi\rangle, \quad\hat H_{\mathrm{el}}|\Psi\rangle=E_{\mathrm{el}}|\Psi\rangle, \quad\partial_\alpha\langle\Psi|\Psi\rangle=0. \end{aligned} \tag{3.511} $$

式（3.511）中波函数导数项相消的原因是本征方程和归一化条件；本征值和一般会随核位置改变。若使用自洽的单电子有效哈密顿量，还须在总能中保留 Hartree 与交换关联的双计数修正。其形式为

$$ \hat{H}=-\frac{1}{2}\boldsymbol{\nabla}^{2}+V_{\mathrm{ee}}+V_{\mathrm{ext}}+V_{\mathrm{xc}} $$

在固定密度的显式核坐标导数中，局域外势贡献可由密度积分给出；若存在非局域赝势，还须对其投影算符求导。因此局域部分为

$$ \boldsymbol F_\alpha^{\mathrm{el}} =-\int\rho(\boldsymbol r)\frac{\partial V_{\mathrm{ext}}(\boldsymbol r;\boldsymbol R)}{\partial\boldsymbol R_\alpha}\,\mathrm d\boldsymbol r =\int\rho(\boldsymbol r)\boldsymbol\nabla_{\boldsymbol r}v_\alpha(\boldsymbol r-\boldsymbol R_\alpha)\,\mathrm d\boldsymbol r \quad\text{（局域外势）} \tag{3.512} $$

而 $\boldsymbol{F}_\alpha^{\mathrm{ion}}$ 的计算比较简单，下面直接给出 $\boldsymbol{F}_\alpha$ 的公式：

$$ \boldsymbol{F}_\alpha^{\mathrm{HF}}=\int\rho(\boldsymbol{r})\left.\frac{\mathrm{d}v_\alpha}{\mathrm{d}\boldsymbol{r}'}\right|_{\boldsymbol{r}'=\boldsymbol{r}-\boldsymbol{R}_\alpha}\mathrm{d}\boldsymbol{r}+\sum_{\beta,\beta\neq\alpha}\frac{Z_\alpha Z_\beta(\boldsymbol{R}_\alpha-\boldsymbol{R}_\beta)}{|\,\boldsymbol{R}_\alpha-\boldsymbol{R}_\beta\,|^{3}} \tag{3.513} $$

因此，原子 $\alpha$ 所受的力即为静电力。方程（3.513）称为 Hellmann-Feynman 表达式，它是第一性原理动力学以及体系弛豫的理论基础。然而，在实际应用中，以上述 Hellmann-Feynman 表达式的形式直接获取精确值通常是不可行的，这主要出于两个原因。首先，我们无法采用严格意义上的完备基组来构建哈密顿矩阵，因为完备基组由无穷多个基函数组成，而在计算过程中我们需要通过截断方法选取有限数量的基函数。其次，我们无法达到“完全自洽”的条件，即对角化哈密顿矩阵后得到完全精确的本征波函数。这两方面原因造成的误差均会影响原子受力的结果，因此需要分别对这两方面的因素进行修正。

再次写出体系的总能表达式

$$ E_{\mathrm{tot}}=T_{\mathrm s}[\rho]+E_{\mathrm H}[\rho]+E_{\mathrm{xc}}[\rho]+E_{\mathrm{Ie}}[\rho]+E_{\mathrm{II}} \tag{3.514} $$

其中对动能泛函 $T_0$ 需要做特别考虑。根据 Kohn-Sham 方程，有 $\left(-\dfrac{1}{2}\boldsymbol{\nabla}^{2}+V_{\mathrm{eff}}(\boldsymbol{r})\right)\psi_i=\varepsilon_i\psi_i$，而 $T=\displaystyle\sum_{i}\langle\psi_i\,|-\frac{1}{2}\boldsymbol{\nabla}^{2}\,|\,\psi_i\rangle$。因此由这两个方程可得

$$ T[\rho]=\sum_{i}n_i\varepsilon_i-\sum_{i}\langle\psi_i\,|\,V_{\mathrm{eff}}(\boldsymbol{r})\,|\,\psi_i\rangle \tag{3.515} $$

不难看出，原子坐标 $\boldsymbol{R}_\alpha$ 对 $E_{\mathrm{tot}}$ 的影响并不仅限于 $E_{\mathrm{II}}$ 和 $E_{\mathrm{Ie}}$，它的变化也会影响本征函数 $\psi_i$，从而导致电子密度 $\rho(\boldsymbol{r})$ 乃至本征能级 $\varepsilon_i$ 的变化。因此，为求得原子受力的普遍表达式，对方程（3.514）求关于 $\boldsymbol{R}_\alpha$ 的全微分\cite{bendt1983simultaneous,srivastava1987theory}：

$$ \boldsymbol F_\alpha=-\frac{\mathrm dE_{\mathrm{tot}}}{\mathrm d\boldsymbol R_\alpha} =-\frac{\mathrm d}{\mathrm d\boldsymbol R_\alpha} \left[T_s[\rho]+E_{\mathrm H}[\rho]+E_{\mathrm{xc}}[\rho] +E_{\mathrm{Ie}}[\rho,\{\boldsymbol R_I\}]+E_{\mathrm{II}}(\{\boldsymbol R_I\})\right] \tag{3.516} $$

式（3.516）右端最后一项在整体负号作用下为 $\boldsymbol F_\alpha^{\mathrm{ion}}=-\mathrm dE_{\mathrm{II}}/\mathrm d\boldsymbol R_\alpha$。 以下将总能各项对核坐标求导；本征值及密度一般随核坐标变化，不能视为常数。

设 $\hat{H}^{0}=-\dfrac{1}{2}\boldsymbol{\nabla}^{2}+V_{\mathrm{eff}}(\boldsymbol{r})$，且 $\psi_i$ 是严格满足 $\hat{H}^{0}\psi_i=\varepsilon_i\psi_i$ 的本征函数，而 $\hat{H}^{0}$ 的表达式中仅 $V_{\mathrm{eff}}$ 与 $\boldsymbol{R}_\alpha$ 有关。由本征方程可得

$$ \frac{\mathrm d\varepsilon_i}{\mathrm d\boldsymbol R_\alpha} =\left\langle\psi_i\left|\frac{\mathrm d\hat H^0}{\mathrm d\boldsymbol R_\alpha}\right|\psi_i\right\rangle \quad\text{（归一化的精确本征态）} \tag{3.517} $$

第二项的计算比较直接：

$$ \begin{aligned} \frac{\mathrm d}{\mathrm d\boldsymbol R_\alpha}\sum_i n_i\langle\psi_i|V_{\mathrm{eff}}|\psi_i\rangle ={}&\sum_i n_i\langle\psi_i|\partial_\alpha V_{\mathrm{eff}}|\psi_i\rangle\\ &+\int V_{\mathrm{eff}}(\boldsymbol r)\,\partial_\alpha\rho(\boldsymbol r)\,\mathrm d\boldsymbol r \end{aligned} \tag{3.518} $$

Hartree 和交换关联泛函对核位置没有显式依赖，但它们通过电子密度产生隐式依赖；外势及非局域赝势还具有显式位置依赖。因此

$$ \begin{aligned} \frac{\mathrm d}{\mathrm d\boldsymbol R_\alpha} \bigl(E_{\mathrm H}[\rho]+E_{\mathrm{xc}}[\rho]+E_{\mathrm{Ie}}[\rho,\boldsymbol R]\bigr) ={}&\int\bigl(V_{\mathrm H}+V_{\mathrm{xc}}+V_{\mathrm{ext}}\bigr) \partial_\alpha\rho\,\mathrm d\boldsymbol r\\ &+\int\rho(\boldsymbol r)\,\partial_\alpha V_{\mathrm{ext}}(\boldsymbol r)\,\mathrm d\boldsymbol r +\left.\partial_\alpha E_{\mathrm{Ie}}^{\mathrm{nl}}\right|_\rho \end{aligned} \tag{3.519} $$

式中

$$ V_{\mathrm{KS}}=V_{\mathrm{ee}}+V_{\mathrm{xc}}+V_{\mathrm{ext}} $$

将式（3.517）至式（3.519）代入式（3.516），并设 $\psi_i$ 由基函数 $\{\chi_j\}$ 展开（$\psi_i=\displaystyle\sum_{j}a_{ij}\chi_j$），则有

$$ \begin{aligned} \boldsymbol{F}_\alpha&=-2\sum_{i}n_i\mathrm{Re}\sum_{j}a_{ij}\langle\frac{\mathrm{d}\chi_j}{\mathrm{d}\boldsymbol{R}_\alpha}\,|\,\hat{H}^{0}-\varepsilon_i\,|\,\psi_i\rangle-\int(V_{\mathrm{KS}}-V_{\mathrm{eff}})\frac{\mathrm{d}\rho(\boldsymbol{r})}{\mathrm{d}\boldsymbol{R}_\alpha}\,\mathrm{d}\boldsymbol{r}+\boldsymbol{F}_\alpha^{\mathrm{el}}+\boldsymbol{F}_\alpha^{\mathrm{ion}}\\ &=\boldsymbol{F}_\alpha^{\mathrm{IBS}}+\boldsymbol{F}_\alpha^{\mathrm{NSF}}+\boldsymbol{F}_\alpha^{\mathrm{HF}} \end{aligned} \tag{3.520} $$

式（3.520）右端第一项是对非完备基的修正，也被称为 Pulay 力\cite{pulay1969initio}；第二项是对非自洽场计算的修正。可以看到：若迭代精度无限高，即 $V_{\mathrm{eff}}$ 严格等于 $V_{\mathrm{KS}}$，则第二项 $\boldsymbol{F}_\alpha^{\mathrm{NSF}}$ 为零；若本征函数严格满足 Hellmann-Feynman 条件，则第一项 $\boldsymbol{F}_\alpha^{\mathrm{IBS}}$ 为零。此外，若基函数 $\chi_j$ 是平面波，$\boldsymbol{F}_\alpha^{\mathrm{IBS}}$ 也为零，因为 $\chi_j$ 不依赖于原子坐标，所以 $\mathrm{d}\chi_j/\mathrm{d}\boldsymbol{R}_\alpha=0$。方程（3.520）只是普遍表达式，对于具体的计算方法，需要依照具体情况提出最适合的表达式，例如文献\cite{savrasov1992full}给出了 LMTO 方法中原子受力的表达式，在这里给出平面波-赝势框架下的 Hellmann-Feynman 力表达式\cite{ihm1979momentum}：

$$ \begin{aligned} \boldsymbol{F}_k={}&\sum_{j,j\neq k}\frac{Z_jZ_k(\boldsymbol{R}_k-\boldsymbol{R}_j)}{|\,\boldsymbol{R}_k-\boldsymbol{R}_j\,|^{3}}-\mathrm{i}\varOmega_{\mathrm{cell}}\sum_{i,l,\boldsymbol{G},\boldsymbol{G}'}(\boldsymbol{G}'-\boldsymbol{G})\mathrm{e}^{\mathrm{i}(\boldsymbol{G}'-\boldsymbol{G})\cdot\boldsymbol{R}_k}\\ &\times\psi^{*}(\boldsymbol{k}_i+\boldsymbol{G})\psi(\boldsymbol{k}_i+\boldsymbol{G}')V_{\mathrm{ps},l,\boldsymbol{k}_i+\boldsymbol{G},\boldsymbol{k}_i+\boldsymbol{G}'} \end{aligned} \tag{3.521} $$

系列导航

1. 分子轨道理论与 Hartree-Fock 方法
2. 均匀电子气、基组选取与超越 Hartree-Fock 近似
3. 密度泛函理论：从托马斯-费米模型到 LDA+U
4. 赝势：正交化平面波、模守恒赝势与超软赝势
5. 平面波-赝势方法：布里渊区积分
6. 平面波-赝势框架下的体系总能与 Ewald 求和
7. 自洽场计算、迭代对角化与 Hellmann-Feynman 力（本文）
8. 缀加平面波方法及其线性化
9. 过渡态搜索：拖曳法、NEB 方法与 Dimer 方法
10. 电子激发谱与准粒子近似：GW 方法与 Bethe-Salpeter 方程
11. 第一性原理计算的应用实例：缺陷、表面与合金相图