在 Abinit 模守恒框架下实现电子局域函数（ELF）：理论推导与测试

动能密度不仅是 meta-GGA 泛函所需要的量，也是计算 ELF 的必要输入。前四节是理论部分，依次讨论平面波表示与对称性、电子密度与动能密度的对称化、ELF 的几种表述，以及动能密度的张量推广；最后两节是测试部分。为了与理论部分区分，测试部分的公式编号记作 (H.n)。

平面波表示下的波函数与对称性

平面波表示

在周期性体系中，波函数 $\psi$ 可以写成 Bloch 波函数的形式：

$$ \psi_{n\mathbf{k}}(\mathbf{r}) = u_{n\mathbf{k}}(\mathbf{r}) e^{i2\pi\mathbf{k}\cdot\mathbf{r}} \tag{1.1} $$

每个 Bloch 波函数由能带编号 $n$¹ 和波矢 $\mathbf{k}$ 标记。$u_{n\mathbf{k}}$ 项赋予 Bloch 波函数与体系相同的周期性：

$$ u_{n\mathbf{k}}(\mathbf{r}+\mathbf{R}) = u_{n\mathbf{k}}(\mathbf{r}) \tag{1.2} $$

其中 $\mathbf{R}$ 是实空间晶格矢量 $\mathbf{R}_{latt}$ 的任意线性组合。于是波函数满足

$$ \psi_{n\mathbf{k}}(\mathbf{r}+\mathbf{R}) = \psi_{n\mathbf{k}}(\mathbf{r}) e^{i2\pi\mathbf{k}\cdot\mathbf{R}} \tag{1.3} $$

$u_{n\mathbf{k}}$（从而波函数本身）可以用平面波展开：

$$ \begin{aligned} u_{n\mathbf{k}}(\mathbf{r}) &= \sum_{\mathbf{G}_\mathbf{k}} c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) e^{i2\pi\mathbf{G}_\mathbf{k}\cdot\mathbf{r}} && (1.4)\\ \psi_{n\mathbf{k}}(\mathbf{r}) &= \sum_{\mathbf{G}_\mathbf{k}} c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) e^{i2\pi(\mathbf{k}+\mathbf{G}_\mathbf{k})\cdot\mathbf{r}} && (1.5) \end{aligned} $$

其中 $e^{i2\pi\mathbf{G}_\mathbf{k}\cdot\mathbf{r}}$ 是平面波，$c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k})$ 是对应的权重，它们其实就是 $u_{n\mathbf{k}}(\mathbf{r})$ 的 Fourier 分量。这里的 $\mathbf{G}_\mathbf{k}$ 是倒格子基矢 $\mathbf{G}_{latt}$ 的任意线性组合²。关于 Abinit 代码的其他理论考虑可参见 \cite{gonze1wf}。

上式中 $\mathbf{G}_\mathbf{k}$ 带有下标 $\mathbf{k}$，原因在于：实际计算中为了避免对 $\mathbf{G}$ 矢量做无穷求和，对每个 $\mathbf{k}$ 只取由截断能 $E_{kin-cut}$ 限定的有限个 $\mathbf{G}$ 矢量。因此 $\mathbf{G}_\mathbf{k}$ 满足

$$ \frac{(2\pi)^2 |\mathbf{G}_\mathbf{k}+\mathbf{k}|^2}{2} < E_{kin-cut} \tag{1.6} $$

按照这个定义，$\mathbf{G}$ 矢量的个数（即平面波个数）对不同的 $\mathbf{k}$ 可以不同（见下图），这正是给 $\mathbf{G}$ 加下标的原因。

因此，严格来说波函数应显式写成

$$ \psi_{n\mathbf{k}}(\mathbf{r}) = \sum_{\mathbf{G}_\mathbf{k}}^{\mathrm{npw}_{\mathbf{k}}} c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) e^{i2\pi(\mathbf{k}+\mathbf{G}_\mathbf{k})\cdot\mathbf{r}} \tag{1.7} $$

最后，由于我们只在第一布里渊区内工作（见下图），波函数还满足一个关系：

$$ \begin{aligned} \psi_{n(\mathbf{k}+\mathbf{G})}(\mathbf{r}) &= \psi_{n\mathbf{k}}(\mathbf{r})\\ \text{从而}\quad u_{n(\mathbf{k}+\mathbf{G})}(\mathbf{r}) &= u_{n\mathbf{k}}(\mathbf{r}) e^{i2\pi\mathbf{G}\cdot\mathbf{r}} \end{aligned} $$

不过需要说明：当我们写两个波函数相等时，总应当考虑到它们的所有线性组合也同样成立。因此在这种情形下引入记号 $\overset{L.C.}{=}$³：

$$ \begin{aligned} \psi_{n(\mathbf{k}+\mathbf{G})}(\mathbf{r}) &\overset{L.C.}{=} \psi_{n\mathbf{k}}(\mathbf{r}) && (1.8)\\ u_{n(\mathbf{k}+\mathbf{G})}(\mathbf{r}) &\overset{L.C.}{=} u_{n\mathbf{k}}(\mathbf{r}) e^{i2\pi\mathbf{G}\cdot\mathbf{r}} && (1.9) \end{aligned} $$

对称性

对于第一性原理（ab initio）代码，内存和计算量往往是瓶颈，因此要尽量把计算压缩到最小代价。这就需要充分利用问题本身的对称性⁴以及体系的对称性（即晶胞可能具有的对称性）。

在实空间中，对称性由一组对称操作 $S_t$⁵ 定义 \cite{burns1990space}。这些算符作用在实空间矢量上：

$$ S_t(\mathbf{r}) = \mathbf{r}' \tag{1.10} $$

例如，$\mathbf{r}$ 指向某个原子，$\mathbf{r}'$ 则指向与之对称的原子。可能的对称操作包括：恒等、反演、（镜面）反射、旋转、旋转–反演⁶、平移、螺旋旋转⁷和滑移反射⁸。其中前几种——恒等、反演、反射、旋转和旋转–反演——称为点式（symmorphic）操作 $(S)$，可以用 $3\times 3$ 矩阵 $S_{\alpha\beta}$ 表示：

$$ \mathbf{r}' = S(\mathbf{r}) \qquad\Rightarrow\qquad r'_\alpha = \sum_\beta S_{\alpha\beta} r_\beta \tag{1.11} $$

对于平移操作以及与平移组合的操作（螺旋旋转和滑移反射），还需要一个平移矢量 $\mathbf{t}$。对于单纯的平移：

$$ \mathbf{r}' = \mathbf{r} + \mathbf{t} \qquad\Rightarrow\qquad r'_\alpha = r_\alpha + t_\alpha \tag{1.12} $$

对于螺旋旋转和滑移反射：

$$ \mathbf{r}' = S(\mathbf{r}) + \mathbf{t} \qquad\Rightarrow\qquad r'_\alpha = \sum_\beta S_{\alpha\beta} r_\beta + t_\alpha \tag{1.13} $$

出于实用考虑，最好把所有类型的对称操作统一到一种表示中。因此我们采用最一般的形式（即螺旋旋转和滑移反射的形式），定义广义算符 $S_\mathbf{t}$⁹：

$$ \mathbf{r}' = S_\mathbf{t}(\mathbf{r}) \qquad\Rightarrow\qquad r'_\alpha = \sum_\beta S_{\alpha\beta} r_\beta + t_\alpha \tag{1.14} $$

对于点式操作，$\mathbf{t}$ 为 $\vec{0}$；对于纯平移操作，矩阵 $S_{\alpha\beta}$ 为单位矩阵。把这些对称算符作用到波函数上，得到

$$ \psi'_{n\mathbf{k}} \overset{L.C.}{=} S_\mathbf{t}\left(\psi_{n\mathbf{k}}\right) \tag{1.15} $$

此时 $S_\mathbf{t}$ 作用的不再是实空间，而是波函数空间。把波函数看作波函数空间中的一个矢量，$S_\mathbf{t}$ 把它变换为另一个波函数矢量。我们的目标是弄清楚这个对称化后的波函数与原波函数之间是什么关系。由于波函数是实空间的函数，有

$$ \left(S_\mathbf{t}\left(\psi_{n\mathbf{k}}\right)\right)(\mathbf{r}) \overset{L.C.}{=} \psi'_{n\mathbf{k}}(\mathbf{r}) \overset{L.C.}{=} \psi_{n\mathbf{k}}\left(S_\mathbf{t}^{-1}(\mathbf{r})\right) \tag{1.16} $$

注意这与下面的操作并不相同：

$$ S_\mathbf{t}\left(\psi_{n\mathbf{k}}(\mathbf{r})\right) \overset{L.C.}{=} \psi'_{n\mathbf{k}}\left(S_\mathbf{t}(\mathbf{r})\right) \overset{L.C.}{=} \psi_{n\mathbf{k}}(\mathbf{r}) \tag{1.17} $$

在式 (1.16) 中，我们把对称算符的逆 $S_\mathbf{t}^{-1}$ 作用到实空间矢量 $\mathbf{r}$ 上，其定义为

$$ \mathbf{r}' = S_\mathbf{t}^{-1}(\mathbf{r}) \qquad\Rightarrow\qquad r'_\alpha = \sum_\beta S^{-1}_{\alpha\beta}\left(r_\beta - t_\beta\right) \tag{1.18} $$

现在来看对称算符作用后的波函数在 $\mathbf{r}+\mathbf{R}$ 处的取值：

$$ \left(S_\mathbf{t}\left(\psi_{n\mathbf{k}}\right)\right)(\mathbf{r}+\mathbf{R}) \overset{L.C.}{=} \psi'_{n\mathbf{k}}(\mathbf{r}+\mathbf{R}) \overset{L.C.}{=} \psi_{n\mathbf{k}}\left(S_\mathbf{t}^{-1}(\mathbf{r}+\mathbf{R})\right) \tag{1.19} $$

其中 $S_\mathbf{t}^{-1}(\mathbf{r}+\mathbf{R})$ 为

$$ \begin{aligned} \mathbf{r}'' = S_\mathbf{t}^{-1}(\mathbf{r}+\mathbf{R}) \qquad\Rightarrow\qquad r''_\alpha &= \sum_\beta S^{-1}_{\alpha\beta}\left(r_\beta + R_\beta - t_\beta\right)\\ r''_\alpha &= \sum_\beta S^{-1}_{\alpha\beta}\left(r_\beta - t_\beta\right) + \sum_\beta S^{-1}_{\alpha\beta} R_\beta\\ r''_\alpha &= r'_\alpha + \tilde{R}_\alpha \end{aligned} $$

$$ \begin{aligned} \mathbf{r}'' &= \mathbf{r}' + \tilde{\mathbf{R}}\\ \mathbf{r}'' &= S_\mathbf{t}^{-1}(\mathbf{r}) + S^{-1}(\mathbf{R}) \end{aligned} \tag{1.20} $$

于是，结合式 (1.3) 与上面的结果，有

$$ \begin{aligned} \left(S_\mathbf{t}\left(\psi_{n\mathbf{k}}\right)\right)(\mathbf{r}+\mathbf{R}) \overset{L.C.}{=} \psi_{n\mathbf{k}}(\mathbf{r}'') &\overset{L.C.}{=} \psi_{n\mathbf{k}}\left(S_\mathbf{t}^{-1}(\mathbf{r}) + S^{-1}(\mathbf{R})\right)\\ &\overset{L.C.}{=} \psi_{n\mathbf{k}}\left(S_\mathbf{t}^{-1}(\mathbf{r})\right) e^{i2\pi\mathbf{k}\cdot\left(S^{-1}(\mathbf{R})\right)} \end{aligned} \tag{1.21} $$

最后一步是改写相位因子 $e^{i2\pi\mathbf{k}\cdot\left(S^{-1}(\mathbf{R})\right)}$，把作用在 $\mathbf{R}$（实空间）上的 $S^{-1}$ 转移到 $\mathbf{k}$（倒空间）上。为此要用到 $S^{-1}$ 的转置：

$$ \begin{aligned} e^{i2\pi\mathbf{k}\cdot\left(S^{-1}(\mathbf{R})\right)} &= e^{i2\pi\sum_\alpha k_\alpha\left(\sum_\beta S^{-1}_{\alpha\beta}R_\beta\right)}\\ e^{i2\pi\mathbf{k}\cdot\left(S^{-1}(\mathbf{R})\right)} &= e^{i2\pi\sum_{\alpha\beta} k_\alpha S^{-1}_{\alpha\beta}R_\beta}\\ e^{i2\pi\mathbf{k}\cdot\left(S^{-1}(\mathbf{R})\right)} &= e^{i2\pi\sum_\beta\left(\sum_\alpha k_\alpha S^{-1,\mathfrak{t}}_{\beta\alpha}\right)R_\beta}\\ e^{i2\pi\mathbf{k}\cdot\left(S^{-1}(\mathbf{R})\right)} &= e^{i2\pi\left(S^{-1,\mathfrak{t}}(\mathbf{k})\right)\cdot\mathbf{R}} \end{aligned} \tag{1.22} $$

记新的矢量为 $\mathbf{k}'$：

$$ \mathbf{k}' = S^{-1,\mathfrak{t}}(\mathbf{k}) \qquad\Rightarrow\qquad k'_\beta = \sum_\alpha S^{-1,\mathfrak{t}}_{\beta\alpha} k_\alpha \tag{1.23} $$

注意，由 $\mathbf{k}$ 得到 $\mathbf{k}'$ 时用的是不含平移的对称算符，即只用到了对称操作的点式部分。最终得到

$$ \begin{aligned} \left(S_\mathbf{t}\left(\psi_{n\mathbf{k}}\right)\right)(\mathbf{r}+\mathbf{R}) \overset{L.C.}{=} \psi'_{n\mathbf{k}}(\mathbf{r}+\mathbf{R}) &\overset{L.C.}{=} \psi_{n\mathbf{k}}\left(S_\mathbf{t}^{-1}(\mathbf{r})\right) e^{i2\pi\mathbf{k}'\cdot\mathbf{R}}\\ \overset{L.C.}{=} \psi'_{n\mathbf{k}}(\mathbf{r}+\mathbf{R}) &\overset{L.C.}{=} \psi'_{n\mathbf{k}}(\mathbf{r}) e^{i2\pi\mathbf{k}'\cdot\mathbf{R}} \end{aligned} \tag{1.24} $$

这实际上说明，对称化后的波函数 $\psi'_{n\mathbf{k}}$ 就是 $\mathbf{k}'$ 处的 Bloch 波函数（参见式 (1.3)），因此

$$ \left(S_\mathbf{t}\left(\psi_{n\mathbf{k}}\right)\right)(\mathbf{r}) \overset{L.C.}{=} \psi'_{n\mathbf{k}}(\mathbf{r}) \overset{L.C.}{=} \psi_{n\mathbf{k}'}(\mathbf{r}) \tag{1.25} $$

Abinit 代码中只存储每个平面波的权重 $c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k})$，所以我们更关心 $c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k})$ 之间的关系：

$$ \begin{aligned} \left(S_\mathbf{t}\left(\psi_{n\mathbf{k}}\right)\right)(\mathbf{r}) &\overset{L.C.}{=} \psi_{n\mathbf{k}'}(\mathbf{r})\\ \psi_{n\mathbf{k}}\left(S_\mathbf{t}^{-1}(\mathbf{r})\right) &\overset{L.C.}{=} \psi_{n\mathbf{k}'}(\mathbf{r})\\ \sum_{\mathbf{G}_\mathbf{k}}^{\mathrm{npw}_{\mathbf{k}}} c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) e^{i2\pi(\mathbf{k}+\mathbf{G}_\mathbf{k})\cdot\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)} &\overset{L.C.}{=} \sum_{\mathbf{G}_{\mathbf{k}'}}^{\mathrm{npw}_{\mathbf{k}'}} c_{n\mathbf{k}'}(\mathbf{G}_{\mathbf{k}'}) e^{i2\pi(\mathbf{k}'+\mathbf{G}_{\mathbf{k}'})\cdot\mathbf{r}} \end{aligned} \tag{1.26} $$

到这里有两点需要说明。

第一，容易看出 $\mathrm{npw}_{\mathbf{k}} = \mathrm{npw}_{\mathbf{k}'}$，因为 $|\mathbf{k}| = |\mathbf{k}'|$（例如参见式 (1.6) 和前面的倒空间示意图）。事实上，$\mathbf{k}'$ 是把点式操作 $S^{-1,\mathfrak{t}}$ 作用于 $\mathbf{k}$ 得到的（即除带平移以外的所有对称操作，它们保持矢量的模不变）。这一点对于最终得到一一对应关系很重要。

第二，$\mathbf{G}_{\mathbf{k}'}$ 与 $\mathbf{G}_\mathbf{k}$ 是什么关系？其实就是（见下图）

$$ \mathbf{G}_{\mathbf{k}'} = \mathbf{G}'_\mathbf{k} = S^{-1,\mathfrak{t}}\left(\mathbf{G}_\mathbf{k}\right) \qquad\Rightarrow\qquad G_{\mathbf{k}',\beta} = G'_{\mathbf{k},\beta} = \sum_\alpha S^{-1,\mathfrak{t}}_{\beta\alpha} G_{\mathbf{k},\alpha} \tag{1.27} $$

于是有

$$ \sum_{\mathbf{G}_\mathbf{k}}^{\mathrm{npw}_{\mathbf{k}}} c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) e^{i2\pi\sum_{\alpha\beta}\left(k_\alpha+G_{\mathbf{k},\alpha}\right)S^{-1}_{\alpha\beta}\left(r_\beta-t_\beta\right)} = \sum_{\mathbf{G}'_\mathbf{k}}^{\mathrm{npw}_{\mathbf{k}'}} c_{n\mathbf{k}'}(\mathbf{G}'_\mathbf{k}) e^{i2\pi\sum_{\alpha\beta}\left(k_\alpha+G_{\mathbf{k},\alpha}\right)S^{-1,\mathfrak{t}}_{\beta\alpha}r_\beta} $$

由此得到

$$ \begin{aligned} c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) e^{-i2\pi\sum_{\alpha\beta}\left(k_\alpha+G_{\mathbf{k},\alpha}\right)S^{-1}_{\alpha\beta}t_\beta} &= c_{n\mathbf{k}'}(\mathbf{G}'_\mathbf{k})\\ c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) e^{-i2\pi\left(\mathbf{k}'+\mathbf{G}'_\mathbf{k}\right)\cdot\mathbf{t}} &= c_{n\mathbf{k}'}(\mathbf{G}'_\mathbf{k})\\ c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) &= c_{n\mathbf{k}'}(\mathbf{G}'_\mathbf{k}) e^{+i2\pi\left(\mathbf{k}'+\mathbf{G}'_\mathbf{k}\right)\cdot\mathbf{t}} \end{aligned} \tag{1.28} $$

或者等价地写成¹⁰

$$ c_{n\mathbf{k}'}(\mathbf{G}'_\mathbf{k}) = c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) e^{+i2\pi(\mathbf{k}+\mathbf{G}_\mathbf{k})\cdot\mathbf{t}} \tag{1.29} $$

回忆广义对称算符 $S_\mathbf{t}$ 的定义（式 (1.14)）。上面的结果意味着，对所有点式对称操作（即 $\mathbf{t} = \vec{0}$ 时），有简单关系

$$ c_{n\mathbf{k}'}(\mathbf{G}'_\mathbf{k}) = c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) \tag{1.30} $$

特别地，对于反演对称操作，有

$$ c_{n(\mathbf{k}'=-\mathbf{k})}\left(\left(\mathbf{G}'_\mathbf{k} = -\mathbf{G}_\mathbf{k}\right)\right) = c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) \tag{1.31} $$

对波函数而言（利用式 (1.25)）：

$$ \psi_{n\mathbf{k}}(-\mathbf{r}) \overset{L.C.}{=} \psi_{n(-\mathbf{k})}(\mathbf{r}) \tag{1.32} $$

作为对称性这一节的收尾，还可以加上一种此前没有考虑的对称性。它不是空间群对称性，而是问题本身的对称性——时间反演对称性（非磁性情形）：

$$ \psi_{n\mathbf{k}}(\mathbf{r}) \overset{L.C.}{=} \psi^{\ast}_{n(-\mathbf{k})}(\mathbf{r}) \tag{1.33} $$

它对 $c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k})$ 给出关系

$$ c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) = c^{\ast}_{n(-\mathbf{k})}(-\mathbf{G}_\mathbf{k}) \tag{1.34} $$

电子密度与动能密度

下面讨论两个物理量：电子密度和动能密度。之所以放在一起讲，是因为二者形式非常相似¹¹。

电子密度

电子密度 $n(\mathbf{r})$ 为

$$ n(\mathbf{r}) = \sum_\mathbf{k}\sum_n f(E_{n\mathbf{k}}) |\psi_{n\mathbf{k}}(\mathbf{r})|^2 \tag{2.1} $$

其中 $f(E)$ 是 Fermi-Dirac 分布。我们在零温下工作，把求和限制在 Fermi 能级 $E_f$ 及其以下的态上¹²：

$$ n(\mathbf{r}) = \sum_\mathbf{k}\sum_n |\psi_{n\mathbf{k}}(\mathbf{r})|^2 \tag{2.2} $$

再考虑到体系可能具有某些对称性，根据上一章的结论，可以进一步把求和约化为对一组不可约 $\mathbf{k}$ 矢量的求和，这对应于所谓的不可约布里渊区（$IBZ$），示例见下图。

因此，电子密度可以切分为若干对称部分，写成对称算符作用在 $IBZ$ 中定义的电子密度部分上之和：

$$ n(\mathbf{r}) = \sum_{S_\mathbf{t}} \left(S_\mathbf{t} n^{IBZ}\right)(\mathbf{r}) = \sum_{S_\mathbf{t}} n^{IBZ}\left(S_\mathbf{t}^{-1}(\mathbf{r})\right) \tag{2.3} $$

其中 $n^{IBZ}$ 是定义在 $IBZ$ 中的那部分电子密度：

$$ n^{IBZ}(\mathbf{r}) = \sum_{\mathbf{k}\in IBZ}\sum_n |\psi_{n\mathbf{k}}(\mathbf{r})|^2 = \sum_{\mathbf{k}\in IBZ}\sum_n n^{IBZ}_{n\mathbf{k}}(\mathbf{r}) \tag{2.4} $$

对于总电子密度，有

$$ \begin{aligned} n(\mathbf{r}) &= \sum_{S_\mathbf{t}}\left(S_\mathbf{t}\sum_{\mathbf{k}\in IBZ}\sum_n n^{IBZ}_{n\mathbf{k}}\right)(\mathbf{r}) && (2.5)\\ n(\mathbf{r}) &= \sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n n^{IBZ}_{n\mathbf{k}}\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)\\ n(\mathbf{r}) &= \sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n\left(\sum_{\mathbf{G}_\mathbf{k}} c^{\ast}_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) e^{-i2\pi(\mathbf{k}+\mathbf{G}_\mathbf{k})\cdot\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)}\sum_{\tilde{\mathbf{G}}_\mathbf{k}} c_{n\mathbf{k}}(\tilde{\mathbf{G}}_\mathbf{k}) e^{i2\pi\left(\mathbf{k}+\tilde{\mathbf{G}}_\mathbf{k}\right)\cdot\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)}\right)\\ n(\mathbf{r}) &= \sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n\left(\sum_{\mathbf{G}_\mathbf{k}}\sum_{\tilde{\mathbf{G}}_\mathbf{k}} c^{\ast}_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) c_{n\mathbf{k}}(\tilde{\mathbf{G}}_\mathbf{k}) e^{i2\pi\left(\tilde{\mathbf{G}}_\mathbf{k}-\mathbf{G}_\mathbf{k}\right)\cdot\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)}\right) && (2.6) \end{aligned} $$

可以直接用最后这个含平面波矢量集 $\mathbf{G}_\mathbf{k}$ 双重求和的式子计算电子密度，但从数值代价来看，更好的做法是从中凑出一个 Fourier 变换，然后借助高效的快速 Fourier 变换（FFT）工具。为此定义一组新的平面波矢量 $\tilde{\tilde{\mathbf{G}}}_\mathbf{k} = \tilde{\mathbf{G}}_\mathbf{k} - \mathbf{G}_\mathbf{k}$，它是一个更大的平面波矢量集¹³（见下图）。

电子密度变为

$$ \begin{aligned} n(\mathbf{r}) &= \sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n\left(\sum_{\mathbf{G}_\mathbf{k}}\sum_{\tilde{\tilde{\mathbf{G}}}_\mathbf{k}} c^{\ast}_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) c_{n\mathbf{k}}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k}+\mathbf{G}_\mathbf{k}) e^{i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)}\right)\\ n(\mathbf{r}) &= \sum_{S_\mathbf{t}}\sum_{\tilde{\tilde{\mathbf{G}}}_\mathbf{k}}\left(\sum_{\mathbf{k}\in IBZ}\sum_n\sum_{\mathbf{G}_\mathbf{k}} c^{\ast}_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) c_{n\mathbf{k}}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k}+\mathbf{G}_\mathbf{k})\right) e^{i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)}\\ n(\mathbf{r}) &= \sum_{S_\mathbf{t}}\sum_{\tilde{\tilde{\mathbf{G}}}_\mathbf{k}} \tilde{n}^{IBZ}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\ e^{i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)} && (2.7) \end{aligned} $$

把这个结果与式 (2.3) 比较，可以清楚地看到 $\tilde{n}^{IBZ}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})$ 就是 $n^{IBZ}\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)$ 的 Fourier 分量。此外，$\tilde{n}^{IBZ}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})$ 与对称算符 $S$ 无关，因此可以在任意一个对称操作下用 FFT 计算。实践中当然选恒等操作，也就是直接使用在 $IBZ$ 中算出的那部分电子密度 $n^{IBZ}(\mathbf{r})$。

接下来的思路是把作用于实空间的对称算符求和（$S_\mathbf{t}^{-1}(\mathbf{r})$）转移到倒空间中去¹⁴：

$$ \begin{aligned} n(\mathbf{r}) &= \sum_{\tilde{\tilde{\mathbf{G}}}_\mathbf{k}}\sum_{S_\mathbf{t}} \tilde{n}^{IBZ}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\ e^{i2\pi\sum_\alpha\tilde{\tilde{G}}_{\mathbf{k},\alpha}\left(\sum_\beta S^{-1}_{\alpha\beta}(r_\beta-t_\beta)\right)}\\ n(\mathbf{r}) &= \sum_{\tilde{\tilde{\mathbf{G}}}_\mathbf{k}}\sum_{S_\mathbf{t}} \tilde{n}^{IBZ}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\ e^{i2\pi\sum_{\alpha\beta}\tilde{\tilde{G}}_{\mathbf{k},\alpha}S^{-1,\mathfrak{t}}_{\beta\alpha}(r_\beta-t_\beta)}\\ n(\mathbf{r}) &= \sum_{\tilde{\tilde{\mathbf{G}}}_\mathbf{k}}\sum_{S_\mathbf{t}} \tilde{n}^{IBZ}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\ e^{i2\pi\left(S^{-1,\mathfrak{t}}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\right)\cdot(\mathbf{r}-\mathbf{t})} && (2.8) \end{aligned} $$

现在花点时间看看 $S^{-1,\mathfrak{t}}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})$ 这一项，并对照上面的 FFT box 示意图。事实上，把任意对称算符 $S^{-1,\mathfrak{t}}$ 作用到任意 $\tilde{\tilde{\mathbf{G}}}_\mathbf{k}$ 上，得到的另一个矢量 $\tilde{\tilde{\mathbf{G}}}'_\mathbf{k}$ 已经包含在 $\tilde{\tilde{\mathbf{G}}}_\mathbf{k}$ 矢量集中，因为平面波矢量集 $\tilde{\tilde{\mathbf{G}}}_\mathbf{k}$ 本身已经包含了全部对称性。因此，先对 $\tilde{\tilde{\mathbf{G}}}_\mathbf{k}$ 施加对称算符再对整个 $\tilde{\tilde{\mathbf{G}}}_\mathbf{k}$ 集求和，与不施加对称算符所得的和相同：

$$ \sum_{\tilde{\tilde{\mathbf{G}}}_\mathbf{k}} e^{i2\pi\left(S^{-1,\mathfrak{t}}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\right)\cdot(\mathbf{r}-\mathbf{t})} = \sum_{\tilde{\tilde{\mathbf{G}}}_\mathbf{k}} e^{i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot(\mathbf{r}-\mathbf{t})} \tag{2.9} $$

回到电子密度：

$$ \begin{aligned} n(\mathbf{r}) &= \sum_{\tilde{\tilde{\mathbf{G}}}_\mathbf{k}}\left(\sum_{S_\mathbf{t}} \tilde{n}^{IBZ}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\ e^{-i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\mathbf{t}}\right) e^{i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\mathbf{r}}\\ n(\mathbf{r}) &= \sum_{\tilde{\tilde{\mathbf{G}}}_\mathbf{k}} \tilde{n}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k}) e^{i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\mathbf{r}} \end{aligned} \tag{2.10} $$

最终结果表明：实空间中的总（即对称化后的）电子密度 $n(\mathbf{r})$¹⁵ 是倒空间中总（对称化）电子密度 $\tilde{n}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})$¹⁶ 的 Fourier 变换，后者定义为

$$ \tilde{n}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k}) = \sum_{S_\mathbf{t}} \tilde{n}^{IBZ}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\ e^{-i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\mathbf{t}} \tag{2.11} $$

相位因子 $e^{-i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\mathbf{t}}$ 称为非点式平移相位（nonsymmorphic translation phase）¹⁷。

动能密度

动能密度 $\tau(\mathbf{r})$ 与电子密度 $n(\mathbf{r})$ 非常相似，区别只在于要对波函数取梯度 $\nabla$：

$$ \tau(\mathbf{r}) = \frac{1}{2}\sum_\mathbf{k}\sum_n f(E_{n\mathbf{k}}) |\nabla\psi_{n\mathbf{k}}(\mathbf{r})|^2 \tag{2.12} $$

同样在零温下工作，把求和限制在 Fermi 能级 $E_f$ 及其以下的态上，得到

$$ \tau(\mathbf{r}) = \frac{1}{2}\sum_\mathbf{k}\sum_n |\nabla\psi_{n\mathbf{k}}(\mathbf{r})|^2 \tag{2.13} $$

我们沿用电子密度的处理流程，只是由于梯度的存在，推导会更复杂一些。但最终会看到，总动能密度同样可以由 $IBZ$ 中定义的动能密度、对称算符和 FFT 构造出来。和式 (2.5) 一样，先从下式开始：

$$ \begin{aligned} \tau(\mathbf{r}) &= \frac{1}{2}\sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n \tau^{IBZ}_{n\mathbf{k}}\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)\\ \tau(\mathbf{r}) &= \frac{1}{2}\sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n \left|\nabla\psi_{n\mathbf{k}}\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)\right|^2\\ \tau(\mathbf{r}) &= \frac{1}{2}\sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n\sum_\gamma \left|\frac{\partial}{\partial\gamma}\psi_{n\mathbf{k}}\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)\right|^2\\ \tau(\mathbf{r}) &= \frac{1}{2}\sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n\sum_\gamma \left|\sum_{\mathbf{G}_\mathbf{k}} c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k})\frac{\partial}{\partial\gamma} e^{i2\pi(\mathbf{k}+\mathbf{G}_\mathbf{k})\cdot\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)}\right|^2\\ \tau(\mathbf{r}) &= \frac{1}{2}\sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n\sum_\gamma \left|\sum_{\mathbf{G}_\mathbf{k}} c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k})\frac{\partial}{\partial\gamma} e^{i2\pi\sum_{\alpha\beta}\left(k_\alpha+G_{\mathbf{k},\alpha}\right)S^{-1}_{\alpha\beta}(r_\beta-t_\beta)}\right|^2\\ \tau(\mathbf{r}) &= \frac{1}{2}\sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n\sum_\gamma \left|\sum_{\mathbf{G}_\mathbf{k}} c_{n\mathbf{k}}(\mathbf{G}_\mathbf{k})\left(i2\pi\sum_\alpha\left(k_\alpha+G_{\mathbf{k},\alpha}\right)S^{-1}_{\alpha\gamma}\right) e^{i2\pi(\mathbf{k}+\mathbf{G}_\mathbf{k})\cdot\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)}\right|^2\\ \tau(\mathbf{r}) &= \frac{1}{2}\sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n\sum_\gamma\\ &\quad\times\sum_{\mathbf{G}_\mathbf{k}}\sum_{\tilde{\mathbf{G}}_\mathbf{k}} c^{\ast}_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) c_{n\mathbf{k}}(\tilde{\mathbf{G}}_\mathbf{k}) (2\pi)^2\left(\sum_\alpha\left(k_\alpha+G_{\mathbf{k},\alpha}\right)S^{-1}_{\alpha\gamma}\right)\left(\sum_\alpha\left(k_\alpha+\tilde{G}_{\mathbf{k},\alpha}\right)S^{-1}_{\alpha\gamma}\right)\\ &\quad\times e^{i2\pi\left(\tilde{\mathbf{G}}_\mathbf{k}-\mathbf{G}_\mathbf{k}\right)\cdot\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)}\\ \tau(\mathbf{r}) &= \frac{1}{2}\sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n\sum_\gamma\\ &\quad\times\sum_{\mathbf{G}_\mathbf{k}}\sum_{\tilde{\tilde{\mathbf{G}}}_\mathbf{k}} c^{\ast}_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) c_{n\mathbf{k}}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k}+\mathbf{G}_\mathbf{k}) (2\pi)^2\left(\sum_\alpha\left(k_\alpha+G_{\mathbf{k},\alpha}\right)S^{-1}_{\alpha\gamma}\right)\\ &\quad\times\left(\sum_\alpha\left(k_\alpha+\tilde{\tilde{G}}_{\mathbf{k},\alpha}+G_{\mathbf{k},\alpha}\right)S^{-1}_{\alpha\gamma}\right) e^{i2\pi\left(\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\right)\cdot\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)} \end{aligned} \tag{2.14} $$

然后用与式 (2.9) 相同的技巧改写相位因子：

$$ \begin{aligned} \tau(\mathbf{r}) &= \frac{1}{2}\sum_{\tilde{\tilde{\mathbf{G}}}_\mathbf{k}}\sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n\sum_{\mathbf{G}_\mathbf{k}} c^{\ast}_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) c_{n\mathbf{k}}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k}+\mathbf{G}_\mathbf{k})\\ &\quad\times(2\pi)^2\sum_\gamma\left(\sum_\alpha\left(k_\alpha+G_{\mathbf{k},\alpha}\right)S^{-1}_{\alpha\gamma}\right)\left(\sum_\alpha\left(k_\alpha+\tilde{\tilde{G}}_{\mathbf{k},\alpha}+G_{\mathbf{k},\alpha}\right)S^{-1}_{\alpha\gamma}\right) e^{-i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\mathbf{t}}\\ &\quad\times e^{i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\mathbf{r}} \end{aligned} \tag{2.15} $$

形式几乎与式 (2.10) 相同，即

$$ \tau(\mathbf{r}) = \frac{1}{2}\sum_{\tilde{\tilde{\mathbf{G}}}_\mathbf{k}} \tilde{\tau}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\ e^{i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\mathbf{r}} \tag{2.16} $$

其中

$$ \tilde{\tau}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k}) = \sum_{S_\mathbf{t}} \tilde{\tau}^{IBZ}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\ e^{-i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\mathbf{t}} \tag{2.17} $$

但这一次，$\tau^{IBZ}\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)$ 的 Fourier 变换 $\tilde{\tau}^{IBZ}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})$ 看上去并不与 $S$ 无关：

$$ \begin{aligned} \tilde{\tau}^{IBZ}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k}) &= \sum_{\mathbf{k}\in IBZ}\sum_n\sum_{\mathbf{G}_\mathbf{k}} c^{\ast}_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) c_{n\mathbf{k}}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k}+\mathbf{G}_\mathbf{k})\\ &\quad\times(2\pi)^2\sum_\gamma\left(\sum_\alpha\left(k_\alpha+G_{\mathbf{k},\alpha}\right)S^{-1}_{\alpha\gamma}\right)\left(\sum_\alpha\left(k_\alpha+\tilde{\tilde{G}}_{\mathbf{k},\alpha}+G_{\mathbf{k},\alpha}\right)S^{-1}_{\alpha\gamma}\right) \end{aligned} \tag{2.18} $$

如果真是这样，我们就不能直接用 $IBZ$ 中算出的那部分动能密度 $\tau^{IBZ}(\mathbf{r})$ 来计算 $\tilde{\tau}^{IBZ}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})$，而必须对每个对称操作分别计算。幸运的是，可以利用 $S$ 矩阵的正交性证明式 (2.18) 与 $S$ 无关，从而等价于

$$ \begin{aligned} \tilde{\tau}^{IBZ}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\Big|_{S=1} &= \sum_{\mathbf{k}\in IBZ}\sum_n\sum_{\mathbf{G}_\mathbf{k}} c^{\ast}_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) c_{n\mathbf{k}}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k}+\mathbf{G}_\mathbf{k})\\ &\quad\times(2\pi)^2\sum_\gamma\left(k_\gamma+G_{\mathbf{k},\gamma}\right)\left(k_\gamma+\tilde{\tilde{G}}_{\mathbf{k},\gamma}+G_{\mathbf{k},\gamma}\right) \end{aligned} \tag{2.19} $$

把式 (2.18) 中含 $S$ 的三重求和单独拿出来改写：

$$ \begin{aligned} &\sum_\gamma\left(\sum_\alpha\left(k_\alpha+G_{\mathbf{k},\alpha}\right)S^{-1}_{\alpha\gamma}\right)\left(\sum_\alpha\left(k_\alpha+\tilde{\tilde{G}}_{\mathbf{k},\alpha}+G_{\mathbf{k},\alpha}\right)S^{-1}_{\alpha\gamma}\right) && (2.20)\\ ={}& \sum_{\gamma\alpha\alpha'} S^{-1}_{\alpha\gamma}S^{-1}_{\alpha'\gamma}\left(k_\alpha+G_{\mathbf{k},\alpha}\right)\left(k_{\alpha'}+\tilde{\tilde{G}}_{\mathbf{k},\alpha'}+G_{\mathbf{k},\alpha'}\right) && (2.21)\\ ={}& \sum_{\alpha\alpha'}\left(\sum_\gamma S^{-1}_{\alpha\gamma}S^{-1}_{\alpha'\gamma}\right)\left(k_\alpha+G_{\mathbf{k},\alpha}\right)\left(k_{\alpha'}+\tilde{\tilde{G}}_{\mathbf{k},\alpha'}+G_{\mathbf{k},\alpha'}\right) && (2.22) \end{aligned} $$

这里可以利用：所有点式对称矩阵 $S$ 都是正交矩阵，其逆矩阵等于转置矩阵，因此

$$ \begin{aligned} S^{-1} &= S^{\mathfrak{t}} && (2.23)\\ SS^{\mathfrak{t}} &= 1 && (2.24)\\ \text{或}\quad\sum_\gamma S_{\gamma\alpha}S_{\gamma\alpha'} &= \delta_{\alpha\alpha'} && (2.25)\\ \text{亦即}\quad\sum_\gamma S^{-1}_{\alpha\gamma}S^{-1}_{\alpha'\gamma} &= \delta_{\alpha\alpha'} && (2.26) \end{aligned} $$

于是

$$ \begin{aligned} &\sum_\gamma\left(\sum_\alpha\left(k_\alpha+G_{\mathbf{k},\alpha}\right)S^{-1}_{\alpha\gamma}\right)\left(\sum_\alpha\left(k_\alpha+\tilde{\tilde{G}}_{\mathbf{k},\alpha}+G_{\mathbf{k},\alpha}\right)S^{-1}_{\alpha\gamma}\right) && (2.27)\\ ={}& \sum_\alpha\left(k_\alpha+G_{\mathbf{k},\alpha}\right)\left(k_\alpha+\tilde{\tilde{G}}_{\mathbf{k},\alpha}+G_{\mathbf{k},\alpha}\right) && (2.28) \end{aligned} $$

它与 $S$ 无关，因此动能密度的对称化可以采用与电子密度相同的技术¹⁸。

电子局域函数（ELF）

在前面讨论了电子密度和动能密度这两个重要的密度之后，就可以构造一个在化学键分析中常用的函数：电子局域函数（Electron Localization Function，$ELF$）。在量子化学拓扑学中，研究化学键通常有两种方法：一是 Bader 提出的分子中的原子（Atoms in Molecules，$AIM$）理论，二是 Silvi 和 Savin 基于 Becke 与 Edgecombe 给出的 $ELF$ 表达式所提出的 $ELF$ 拓扑分析。

Savin 表述

$ELF$ 的表述要借助一个类似动能密度的新物理量，称之为“Pauli 动能密度”，记为 $D(\mathbf{r})$，定义为

$$ D(\mathbf{r}) = \tau(\mathbf{r}) - \frac{1}{8}\frac{|\nabla n(\mathbf{r})|^2}{n(\mathbf{r})} \tag{3.1} $$

其中 $|\nabla n(\mathbf{r})|^2$ 是电子密度梯度的模方，$\frac{1}{8}\frac{|\nabla n(\mathbf{r})|^2}{n(\mathbf{r})}$ 实际上就是 Weizsäcker 动能密度。当 Pauli 动能密度 $D(\mathbf{r})$ 趋于零时，找到局域电子的概率增大。对于 $D(\mathbf{r})$ 不为零的情形，需要一个参考值来判断电子是否算得上局域。我们采用 Thomas-Fermi 动能密度 $D^0(\mathbf{r})$ 作为参考：

$$ D^0(\mathbf{r}) = \frac{3}{10}\left(3\pi^2\right)^{2/3} n^{5/3}(\mathbf{r}) = C_F n^{5/3}(\mathbf{r}) \tag{3.2} $$

其中 $C_F \sim 2.871$ 为 Fermi 常数。电子局域函数于是定义为

$$ ELF(\mathbf{r}) = \frac{1}{1+\left(\dfrac{D(\mathbf{r})}{D^0(\mathbf{r})}\right)^2} \tag{3.3} $$

按此定义，$ELF(\mathbf{r})$ 是一个介于 0 和 1 之间的无量纲量。当 $ELF(\mathbf{r}) \sim 1$ 时，该区域的电子是局域的；当 $ELF(\mathbf{r}) \sim \frac{1}{2}$ 时，意味着 $D(\mathbf{r}) = D^0(\mathbf{r})$，即电子的局域程度并不比均匀电子气更高。

自旋相关的 ELF

如果考虑自旋相关的密度，就可以构造自旋相关的 $ELF$，对每个自旋分量都由同样的一般公式给出：

$$ ELF_\sigma(\mathbf{r}) = \frac{1}{1+\left(\dfrac{D_\sigma(\mathbf{r})}{D^0_\sigma(\mathbf{r})}\right)^2} \tag{3.4} $$

事实上，Becke 和 Edgecombe 给出的 ELF 原始表述正是建立在同自旋电子的电子对密度之上的。由于 Pauli 原理，可以认为同自旋电子对之间的关联强于异自旋电子对。此时需要如下定义的自旋相关密度：

$$ \begin{aligned} n_\sigma(\mathbf{r}) &= \sum_\mathbf{k}\sum_{n_\sigma} f(E_{n_\sigma\mathbf{k}}) |\psi_{n_\sigma\mathbf{k}}(\mathbf{r})|^2 && (3.5)\\ \tau_\sigma(\mathbf{r}) &= \frac{1}{2}\sum_\mathbf{k}\sum_{n_\sigma} f(E_{n_\sigma\mathbf{k}}) |\nabla\psi_{n_\sigma\mathbf{k}}(\mathbf{r})|^2 && (3.6) \end{aligned} $$

它们满足¹⁹

$$ \begin{aligned} n(\mathbf{r}) &= n_{\sigma=\uparrow}(\mathbf{r}) + n_{\sigma=\downarrow}(\mathbf{r}) && (3.7)\\ \tau(\mathbf{r}) &= \tau_{\sigma=\uparrow}(\mathbf{r}) + \tau_{\sigma=\downarrow}(\mathbf{r}) && (3.8) \end{aligned} $$

Becke–Edgecombe 表述与 Savin 表述的关系

Becke 和 Edgecombe（B-E）对每个自旋的 ELF 表述，基于下面的 Pauli 动能密度和 Thomas-Fermi 动能密度公式：

$$ \begin{aligned} D^{\mathrm{B-E}}_\sigma(\mathbf{r}) &= 2\tau_\sigma(\mathbf{r}) - \frac{1}{4}\frac{|\nabla n_\sigma(\mathbf{r})|^2}{n_\sigma(\mathbf{r})} && (3.9)\\ D^{0\ \mathrm{B-E}}_\sigma(\mathbf{r}) &= C_F\left(2n_\sigma(\mathbf{r})\right)^{5/3} && (3.10) \end{aligned} $$

式 (3.1) 和式 (3.2) 给出的 Savin（S）表述与自旋无关，对于闭壳层体系它与 Becke–Edgecombe 表述一致。事实上，对这类体系有

$$ \begin{aligned} \frac{1}{2}n(\mathbf{r}) &= n_{\sigma=\uparrow}(\mathbf{r}) = n_{\sigma=\downarrow}(\mathbf{r}) && (3.11)\\ \frac{1}{2}\tau(\mathbf{r}) &= \tau_{\sigma=\uparrow}(\mathbf{r}) = \tau_{\sigma=\downarrow}(\mathbf{r}) && (3.12) \end{aligned} $$

$$ \begin{aligned} D^{\mathrm{S}}(\mathbf{r}) &= \tau(\mathbf{r}) - \frac{1}{8}\frac{|\nabla n(\mathbf{r})|^2}{n(\mathbf{r})}\\ &= 2\tau_\sigma(\mathbf{r}) - \frac{1}{8}\frac{|\nabla 2n_\sigma(\mathbf{r})|^2}{2n_\sigma(\mathbf{r})}\\ &= 2\tau_\sigma(\mathbf{r}) - \frac{1}{4}\frac{|\nabla n_\sigma(\mathbf{r})|^2}{n_\sigma(\mathbf{r})} = D^{\mathrm{B-E}}_\sigma(\mathbf{r}) && (3.13)\\ D^{0\ \mathrm{S}}(\mathbf{r}) &= C_F n^{5/3}(\mathbf{r})\\ &= C_F\left(2n_\sigma(\mathbf{r})\right)^{5/3} = D^{0\ \mathrm{B-E}}_\sigma(\mathbf{r}) && (3.14) \end{aligned} $$

按照 Becke–Edgecombe 基于同自旋电子对密度的 ELF 表述精神，Savin 表述因此只对闭壳层体系严格成立。不过在没有自旋信息（自旋相关密度）的情况下，Savin 的解释对开壳层体系仍然有用。

Kohout–Savin 表述

由于上述两种表述（B-E 与 S）之间存在容易引起误解的因子 2，Kohout 和 Savin（K-S）按照自旋密度泛函理论，重新定义了自旋相关情形下的 Pauli 动能密度和 Thomas-Fermi 动能密度 \cite{kohout1996atomic}：

$$ \begin{aligned} D^{\mathrm{K-S}}_\sigma(\mathbf{r}) &= \frac{D^{\mathrm{B-E}}_\sigma(\mathbf{r})}{2} = \tau_\sigma(\mathbf{r}) - \frac{1}{8}\frac{|\nabla n_\sigma(\mathbf{r})|^2}{n_\sigma(\mathbf{r})} && (3.15)\\ D^{0\ \mathrm{K-S}}_\sigma(\mathbf{r}) &= \frac{D^{0\ \mathrm{B-E}}_\sigma(\mathbf{r})}{2} = 2^{2/3}C_F n_\sigma^{5/3}(\mathbf{r}) && (3.16) \end{aligned} $$

他们还给出了总 ELF 的新表述。由于这次它建立在自旋相关量之上，原则上对开壳层体系同样适用：

$$ \begin{aligned} D^{\mathrm{K-S}}(\mathbf{r}) &= D^{\mathrm{K-S}}_{\sigma=\uparrow}(\mathbf{r}) + D^{\mathrm{K-S}}_{\sigma=\downarrow}(\mathbf{r})\\ &= \tau_{\sigma=\uparrow}(\mathbf{r}) + \tau_{\sigma=\downarrow}(\mathbf{r}) - \frac{1}{8}\frac{|\nabla n_{\sigma=\uparrow}(\mathbf{r})|^2}{n_{\sigma=\uparrow}(\mathbf{r})} - \frac{1}{8}\frac{|\nabla n_{\sigma=\downarrow}(\mathbf{r})|^2}{n_{\sigma=\downarrow}(\mathbf{r})}\\ &= \tau(\mathbf{r}) - \frac{1}{8}\left(\frac{|\nabla n_{\sigma=\uparrow}(\mathbf{r})|^2}{n_{\sigma=\uparrow}(\mathbf{r})} + \frac{|\nabla n_{\sigma=\downarrow}(\mathbf{r})|^2}{n_{\sigma=\downarrow}(\mathbf{r})}\right) && (3.17)\\ D^{0\ \mathrm{K-S}}(\mathbf{r}) &= D^{0\ \mathrm{K-S}}_{\sigma=\uparrow}(\mathbf{r}) + D^{0\ \mathrm{K-S}}_{\sigma=\downarrow}(\mathbf{r})\\ &= 2^{2/3}C_F\left(n^{5/3}_{\sigma=\uparrow}(\mathbf{r}) + n^{5/3}_{\sigma=\downarrow}(\mathbf{r})\right) && (3.18) \end{aligned} $$

动能密度的张量推广

前面（式 (2.13)）已经给出了动能密度的定义，现在把它推广为张量形式。事实上，动能密度可以看作动能密度张量 $\overline{\tau}$（一个 $3\times 3$ 矩阵）的迹，张量的每个元素定义为

$$ \tau_{\alpha\beta}(\mathbf{r}) = \sum_\mathbf{k}\sum_n\left(\frac{\partial}{\partial\alpha}\psi_{n\mathbf{k}}(\mathbf{r})\right)^{\ast}\left(\frac{\partial}{\partial\beta}\psi_{n\mathbf{k}}(\mathbf{r})\right) \tag{4.1} $$

按此定义，动能密度可表示为

$$ \tau(\mathbf{r}) = \frac{1}{2}\sum_\alpha \tau_{\alpha\alpha}(\mathbf{r}) = \frac{1}{2}\mathrm{Tr}\left[\overline{\tau}(\mathbf{r})\right] \tag{4.2} $$

动能密度张量的单个元素在许多发展中都可能有用，尤其是单个对角元²⁰。和前面处理其他密度一样，下面来看在利用对称性、只在 $IBZ$ 中计算时，如何对动能密度张量的元素做对称化：

$$ \begin{aligned} \tau_{\alpha\beta}(\mathbf{r}) &= \sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n \tau^{IBZ}_{\alpha\beta\ n\mathbf{k}}\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)\\ \tau_{\alpha\beta}(\mathbf{r}) &= \sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n\left(\frac{\partial}{\partial\alpha}\psi_{n\mathbf{k}}\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)\right)^{\ast}\left(\frac{\partial}{\partial\beta}\psi_{n\mathbf{k}}\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)\right)\\ \tau_{\alpha\beta}(\mathbf{r}) &= \sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n\\ &\quad\times\sum_{\mathbf{G}_\mathbf{k}}\sum_{\tilde{\mathbf{G}}_\mathbf{k}} c^{\ast}_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) c_{n\mathbf{k}}(\tilde{\mathbf{G}}_\mathbf{k}) (2\pi)^2\left(\sum_{\alpha'}\left(k_{\alpha'}+G_{\mathbf{k},\alpha'}\right)S^{-1}_{\alpha'\alpha}\right)\left(\sum_{\beta'}\left(k_{\beta'}+\tilde{G}_{\mathbf{k},\beta'}\right)S^{-1}_{\beta'\beta}\right)\\ &\quad\times e^{i2\pi\left(\tilde{\mathbf{G}}_\mathbf{k}-\mathbf{G}_\mathbf{k}\right)\cdot\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)}\\ \tau_{\alpha\beta}(\mathbf{r}) &= \sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n\\ &\quad\times\sum_{\mathbf{G}_\mathbf{k}}\sum_{\tilde{\tilde{\mathbf{G}}}_\mathbf{k}} c^{\ast}_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) c_{n\mathbf{k}}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k}+\mathbf{G}_\mathbf{k}) (2\pi)^2\left(\sum_{\alpha'}\left(k_{\alpha'}+G_{\mathbf{k},\alpha'}\right)S^{-1}_{\alpha'\alpha}\right)\\ &\quad\times\left(\sum_{\beta'}\left(k_{\beta'}+\tilde{\tilde{G}}_{\mathbf{k},\beta'}+G_{\mathbf{k},\beta'}\right)S^{-1}_{\beta'\beta}\right) e^{i2\pi\left(\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\right)\cdot\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)}\\ \tau_{\alpha\beta}(\mathbf{r}) &= \sum_{\tilde{\tilde{\mathbf{G}}}_\mathbf{k}}\sum_{S_\mathbf{t}}\sum_{\mathbf{k}\in IBZ}\sum_n\sum_{\mathbf{G}_\mathbf{k}} c^{\ast}_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) c_{n\mathbf{k}}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k}+\mathbf{G}_\mathbf{k})\\ &\quad\times(2\pi)^2\left(\sum_{\alpha'}\left(k_{\alpha'}+G_{\mathbf{k},\alpha'}\right)S^{-1}_{\alpha'\alpha}\right)\left(\sum_{\beta'}\left(k_{\beta'}+\tilde{\tilde{G}}_{\mathbf{k},\beta'}+G_{\mathbf{k},\beta'}\right)S^{-1}_{\beta'\beta}\right) e^{-i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\mathbf{t}}\\ &\quad\times e^{i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\mathbf{r}} \end{aligned} \tag{4.3} $$

同样可以凑出一个 Fourier 变换：

$$ \tau_{\alpha\beta}(\mathbf{r}) = \sum_{\tilde{\tilde{\mathbf{G}}}_\mathbf{k}} \tilde{\tau}_{\alpha\beta}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\ e^{i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\mathbf{r}} \tag{4.4} $$

其中

$$ \tilde{\tau}_{\alpha\beta}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k}) = \sum_{S_\mathbf{t}} \tilde{\tau}^{IBZ}_{\alpha\beta}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\ e^{-i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\mathbf{t}} \tag{4.5} $$

但这一次，$\tau^{IBZ}_{\alpha\beta}\left(S_\mathbf{t}^{-1}(\mathbf{r})\right)$ 的 Fourier 变换 $\tilde{\tau}^{IBZ}_{\alpha\beta}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})$ 确实与 $S$ 有关：

$$ \begin{aligned} \tilde{\tau}^{IBZ}_{\alpha\beta}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\Big|_{\forall S} &= \sum_{\mathbf{k}\in IBZ}\sum_n\sum_{\mathbf{G}_\mathbf{k}} c^{\ast}_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) c_{n\mathbf{k}}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k}+\mathbf{G}_\mathbf{k})\\ &\quad\times(2\pi)^2\left(\sum_{\alpha'}\left(k_{\alpha'}+G_{\mathbf{k},\alpha'}\right)S^{-1}_{\alpha'\alpha}\right)\left(\sum_{\beta'}\left(k_{\beta'}+\tilde{\tilde{G}}_{\mathbf{k},\beta'}+G_{\mathbf{k},\beta'}\right)S^{-1}_{\beta'\beta}\right) \end{aligned} \tag{4.6} $$

这里已经没有正交关系可以帮忙了。不过，对给定的一对 $\{\alpha,\beta\}$，任意 $S$ 下的 $\tilde{\tau}^{IBZ}_{\alpha\beta}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})$ 都可以改写为 $S=1$ 时各量的求和。事实上，对于恒等操作 $S=1$ 有

$$ \begin{aligned} \tilde{\tau}^{IBZ}_{\alpha\beta}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\Big|_{S=1} &= \sum_{\mathbf{k}\in IBZ}\sum_n\sum_{\mathbf{G}_\mathbf{k}} c^{\ast}_{n\mathbf{k}}(\mathbf{G}_\mathbf{k}) c_{n\mathbf{k}}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k}+\mathbf{G}_\mathbf{k})\\ &\quad\times(2\pi)^2\left(k_\alpha+G_{\mathbf{k},\alpha}\right)\left(k_\beta+\tilde{\tilde{G}}_{\mathbf{k},\beta}+G_{\mathbf{k},\beta}\right) \end{aligned} \tag{4.7} $$

因此

$$ \tilde{\tau}^{IBZ}_{\alpha\beta}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\Big|_{\forall S} = \sum_{\alpha',\beta'} \tilde{\tau}^{IBZ}_{\alpha'\beta'}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\Big|_{S=1} S^{-1}_{\alpha'\alpha}S^{-1}_{\beta'\beta} \tag{4.8} $$

这意味着，要对张量的某一个元素做对称化，需要知道 $IBZ$ 中全部九个元素。最终得到

$$ \tau_{\alpha\beta}(\mathbf{r}) = \sum_{\tilde{\tilde{\mathbf{G}}}_\mathbf{k}}\left(\sum_{S_\mathbf{t}}\left(\sum_{\alpha',\beta'} \tilde{\tau}^{IBZ}_{\alpha'\beta'}(\tilde{\tilde{\mathbf{G}}}_\mathbf{k})\Big|_{S=1} S^{-1}_{\alpha'\alpha}S^{-1}_{\beta'\beta}\right) e^{-i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\mathbf{t}}\right) e^{i2\pi\tilde{\tilde{\mathbf{G}}}_\mathbf{k}\cdot\mathbf{r}} \tag{4.9} $$

测试：孤立 H 原子

这一部分来自配套的测试报告，用解析结果检验上述 ELF 实现。

计算采用 Fermi-Amaldi 交换关联泛函（$ixc = 20$），不考虑自旋极化（该泛函不支持自旋极化）。

单个 H 原子的波函数是 $1s$ 原子轨道。做解析推导²¹时采用球谐函数形式：

$$ \psi = \varphi_{1s}(r,\theta,\phi) = \sqrt{\frac{Z^3}{\pi a_0^3}} e^{-Z\frac{|\mathbf{r}|}{a_0}} \tag{H.1} $$

其中 $Z$ 为原子序数，$a_0$ 为 Bohr 半径。对 H 原子（$Z = 1$）可以得到：

• 电子密度

  $$ n(\mathbf{r}) = |\psi|^2 = |\varphi_{1s}(r,\theta,\phi)|^2 = \frac{1}{\pi a_0^3} e^{-\frac{2|\mathbf{r}|}{a_0}} \tag{H.2} $$

• 动能密度

  $$ \tau(\mathbf{r}) = \frac{1}{2}|\nabla\psi|^2 = \frac{1}{2}|\nabla\varphi_{1s}(r,\theta,\phi)|^2 = \frac{1}{2\pi a_0^5} e^{-\frac{2|\mathbf{r}|}{a_0}} \tag{H.3} $$

• 电子密度梯度的模方

  $$ |\nabla n(\mathbf{r})|^2 = \left|\nabla|\psi|^2\right|^2 = \left|\nabla|\varphi_{1s}(r,\theta,\phi)|^2\right|^2 = \left|\frac{-2}{\pi a_0^4} e^{-\frac{2|\mathbf{r}|}{a_0}}\right|^2 = \frac{4}{\pi^2 a_0^8} e^{-\frac{4|\mathbf{r}|}{a_0}} \tag{H.4} $$

• Weizsäcker 动能密度

  $$ \frac{1}{8}\frac{|\nabla n(\mathbf{r})|^2}{n(\mathbf{r})} = \frac{1}{8}\frac{\frac{4}{\pi^2 a_0^8} e^{-\frac{4|\mathbf{r}|}{a_0}}}{\frac{1}{\pi a_0^3} e^{-\frac{2|\mathbf{r}|}{a_0}}} = \frac{1}{2\pi a_0^5} e^{-\frac{2|\mathbf{r}|}{a_0}} \tag{H.5} $$

• Thomas-Fermi 动能密度

  $$ \frac{3}{10}\left(3\pi^2\right)^{2/3} n^{5/3}(\mathbf{r}) = 2.871 \times \left(\frac{1}{\pi a_0^3} e^{-\frac{2|\mathbf{r}|}{a_0}}\right)^{5/3} \tag{H.6} $$

• ELF²²

  $$ ELF(\mathbf{r}) = \frac{1}{1+\left(\dfrac{\tau(\mathbf{r})-\frac{1}{8}\frac{|\nabla n(\mathbf{r})|^2}{n(\mathbf{r})}}{2.871\times n^{5/3}(\mathbf{r})}\right)} = \frac{1}{1+\left(\dfrac{0}{2.871\times n^{5/3}(\mathbf{r})}\right)} = 1 \tag{H.7} $$

可以看到，对单个氢原子，$ELF$ 应当处处等于 1，因为此时动能密度与 Weizsäcker 动能密度相等²³（见式 (H.3) 和式 (H.5)）。

标准测试

使用的标准输入文件如下：

acell 3*30
ecut 100
diemac 1.0d0
diemix 0.5d0
iscf 3
ixc 20
kpt 3*0.25
natom 1
nband 1
nkpt 1
nline 3
nsppol 1
nstep 6
nsym 8
ntypat 1
occ 1
rprim 100 010 001
symrel
1 0 0   0 1 0   0 0 1
-1 0 0  0 1 0   0 0 1
1 0 0   0-1 0   0 0 1
-1 0 0  0-1 0   0 0 1
1 0 0   0 1 0   0 0-1
-1 0 0  0 1 0   0 0-1
1 0 0   0-1 0   0 0-1
-1 0 0  0-1 0   0 0-1
tnons 24*0
tolwfr 1.0d-14
typat 1
wtk 1
znucl 1
xred 3*0
prtelf 1 #output a _ELF file.

下面的图给出 ABINIT 计算结果与前述解析公式的对比。首先是对 acell 参数的收敛性：

然后是对 ecut 参数的收敛性。可以发现，至少对氢原子而言，$ELF$ 的收敛对 acell 参数比对 ecut 参数更敏感。例如，只用 10 Ha 的 ecut、但采用 10 Bohr 的盒子，就已经得到处处为 1 的结果。

测试：孤立 Li 原子

氢原子用于测试 $ELF$ 有些特殊，因此我们又对另一个孤立原子做了测试。这里选用锂（Li），因为 $ELF$ 可用来显示孤立原子的壳层结构，而对 Li 来说这一结构非常简单——实际上只有单个 $s$ 壳层。为此我们做了全电子计算²⁴。

首先采用裸（bare）赝势：

然后采用 fhi 赝势：