本文是「第一性原理的微观计算模拟」系列的第 3 篇（共 11 篇），内容整理自同名书稿第 3 章，公式编号与原书一致。文中引用的文献依据原书参考文献清单整理并列于文末，按原书清单顺序从 1 开始编号。\nocite{*}
密度泛函理论（density functional theory，DFT）是一种研究如何将复杂的多体问题严格转化为相对简单的单体问题的理论。与 Hartree-Fock 方法不同，DFT 方法并不关心电子的具体组态，它关注的是体系基态所对应的空间中的电荷密度 $\rho(\boldsymbol{r})$。DFT 方法在当今的量子物理计算、材料模拟和化学反应研究中占据着主导地位，具有较高的精度和较快的计算速度。这一理论分支的发展非常迅速，已经广泛应用于各个领域。
DFT 的核心思想是将多体问题转化为一个单体问题，从而避免了直接处理多电子波函数的困难。在 DFT 框架下，电子的基态能量可以通过电子密度 $\rho(\boldsymbol{r})$ 来表示，而不是通过复杂数量的多电子波函数。这使得 DFT 方法在处理大量电子体系时变得相对简单和高效。
DFT 的基础是 Hohenberg-Kohn 定理和 Kohn-Sham 方程。Hohenberg-Kohn 定理表明，一个多电子体系的基态密度可以唯一确定其基态性质，从而建立了电子密度与多电子哈密顿量之间的一一对应关系。Kohn-Sham 方程将电子间的相互作用转化为一系列单电子方程，这些方程包含了一个有效的外势，包括电子-核吸引势、库仑排斥势以及交换关联势。
解 Kohn-Sham 方程得到的是 Kohn-Sham 轨道，从这些轨道可以获得电子密度，并进一步计算得到体系的基态能量。
尽管 DFT 方法具有很多优点，但它的精确性仍然受限于交换关联势的近似。为了提高 DFT 方法的精度，研究人员发展了许多不同的交换关联泛函，例如局域密度近似（LDA）泛函、广义梯度近似（GGA）泛函和杂化泛函等。这些泛函在不同的应用领域和体系中表现出各自的优缺点，用户需要根据实际问题选择合适的泛函。
在本节中我们将会详细讨论 DFT 的基本理论及具体实现。
托马斯-费米-狄拉克近似
1927 年，托马斯和费米各自独立地提出了将体系能量写为仅显含电子密度 $\rho(\boldsymbol{r})$ 的表达式的方法。这是密度泛函理论的第一次尝试，称为托马斯-费米理论。其主要的思想是利用均匀电子气的解析结果，将体系总能各项，如动能项、Hartree 项等表示成如下形式的方程：
$$ E_i=\int\varepsilon_i[\rho(\boldsymbol{r})]\rho(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r} \tag{3.150} $$
式中：$\varepsilon_i[\rho(\boldsymbol{r})]$ 是均匀电子气模型下各项的能量“密度”（相对于 $\boldsymbol{r}$ 终点处的电子密度 $\rho(\boldsymbol{r})$ 而言）。此外，式（3.150）表示体系能量只与该点处的 $\rho(\boldsymbol{r})$ 有关。为使计算简便，在下面的讨论中设体系的体积为 1。
对于动能项，根据第 2 章中自由电子气模型，有
$$ T=\frac{3}{5}\varepsilon_{\mathrm{F}}\rho \tag{3.151} $$
其中，$\varepsilon_{\mathrm{F}}=\dfrac{\hbar^{2}k_{\mathrm{F}}^{2}}{2m}$ 是费米动能。而
$$ \rho=\frac{1}{3\pi^{2}}\left(\frac{2m}{\hbar^{2}}\right)^{3/2}\varepsilon_{\mathrm{F}}^{3/2} \tag{3.152} $$
由式（3.151）和式（3.152）可得，每个电子的平均动能为
$$ t[\rho]=\frac{T}{\rho}=\frac{3}{5}\frac{\hbar^{2}}{2m}(3\pi^{2})^{2/3}\rho^{2/3} \tag{3.153} $$
代入式（3.150）可得
$$ T[\rho]=C_1\int\rho(\boldsymbol{r})^{5/3}\,\mathrm{d}\boldsymbol{r} \tag{3.154} $$
显然，取原子单位制，并令 ${\hbar=m_{\mathrm e}=e=4\pi\varepsilon_0=1}$，可知
$$ C_1=\frac{3(3\pi^{2})^{2/3}}{10}\approx2.871 \tag{3.155} $$
Hartree 项及电子-核相互作用能的表达式分别为
$$ E_{\mathrm{H}}=\frac{1}{2}\iint\frac{\rho(\boldsymbol{r})\rho(\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}' \tag{3.156} $$
$$ E_{\mathrm{ext}}=\int V_{\mathrm{ext}}(\boldsymbol{r})\rho(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r} \tag{3.157} $$
托马斯和费米在工作中忽略了多体体系的交换关联作用。狄拉克对托马斯-费米理论予以发展，提出电子的交换能密度应满足 $\varepsilon_{\mathrm{x}}\propto\rho^{1/3}$。同样，根据均匀电子气模型，体系的交换能为
$$ E_{\mathrm{x}}=C_2\int\rho(\boldsymbol{r})^{4/3}\,\mathrm{d}\boldsymbol{r} \tag{3.158} $$
在原子单位制下，有
$$ C_2=-\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}\approx-0.739 \tag{3.159} $$
因此，体系的总能可以表示为
$$ E_{\mathrm{TFD}}=2.871\int\rho(\boldsymbol{r})^{5/3}\,\mathrm{d}\boldsymbol{r}+\frac{1}{2}\iint\frac{\rho(\boldsymbol{r})\rho(\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}'+\int V_{\mathrm{ext}}(\boldsymbol{r})\rho(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r}-0.739\int\rho(\boldsymbol{r})^{4/3}\,\mathrm{d}\boldsymbol{r} \tag{3.160} $$
式（3.160）称为托马斯-费米-狄拉克（TFD）近似。可以看到，$E_{\mathrm{TFD}}$ 仅与体系的电子密度分布 $\rho(\boldsymbol{r})$ 有关，可表示为 $\rho(\boldsymbol{r})$ 的泛函。
根据方程（3.160）及约束条件，有
$$ \int\rho(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r}=N \tag{3.161} $$
通过拉格朗日乘子法求解体系的基态能量及相应的电子密度分布：
$$ \frac{\delta\left[E_{\mathrm{TFD}}[\rho]-\mu\left(\displaystyle\int\rho(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r}-N\right)\right]}{\delta\rho(\boldsymbol{r})}=0 \tag{3.162} $$
即可得 TFD 方程
$$ \frac{5}{3}C_1\rho(\boldsymbol r)^{2/3}+V_{\mathrm{ext}}(\boldsymbol r) +\int\frac{\rho(\boldsymbol r')}{|\boldsymbol r-\boldsymbol r'|}\,\mathrm d\boldsymbol r' +\frac{4}{3}C_2\rho(\boldsymbol r)^{1/3}-\mu=0 \tag{3.163} $$
其中拉格朗日乘子 $\mu$ 是电子的化学势，即费米能。在已知 $\mu$ 和 $V_{\mathrm{ext}}(\boldsymbol{r})$ 的情况下，可以反解该方程，得到基态电子密度分布。可以看到，由于积分项及非线性项的存在，这个任务并不容易完成。
根据 TFD 近似，体系的基态能量原则上可以通过对单一函数 $\rho(\boldsymbol{r})$ 的变分求得。相对 Hartree-Fock 方法需要求解 $N$ 个联立的方程组而言，TFD 近似要简单得多，且在碱金属体系的计算上得到了理想的结果。但是因为方程（3.160）的导出是依据均匀电子气模型，且没有考虑电子的关联作用，所以 TFD 近似方法用于处理成键方向性较强的体系（如离子键或共价键体系等）时效果并不理想。严格的密度泛函理论及其实践算法是在 TFD 近似提出三十多年后由 Hohenberg、Kohn 及 Sham 给出的，我们将在下面几节中进行详细的讨论。
Hohenberg-Kohn 定理
Hohenberg–Kohn 定理确立了基态电子密度与外势的关系及能量变分原理，是现代密度泛函理论的基础；构造可求解的单粒子方程还需要 Kohn–Sham 辅助体系\cite{hohenberg1964inhomogeneous}。该定理主要由两部分组成。
定理 3.1 任意一个由相互作用粒子组成的体系所感受到的外势 $V_{\mathrm{ext}}(\boldsymbol{r})$，除了常数因子外，唯一地由该体系的基态电子密度分布 $\rho^{0}(\boldsymbol{r})$ 确定。
推论 3.1 既然 $V_{\mathrm{ext}}(\boldsymbol{r})$ 决定了体系的哈密顿量 $H(\boldsymbol{r})$，而 $V_{\mathrm{ext}}(\boldsymbol{r})$ 又由 $\rho^{0}(\boldsymbol{r})$ 确定，那么体系的多电子基态波函数 $\varPsi^{0}$ 完全由 $\rho^{0}(\boldsymbol{r})$ 确定，是 $\rho^{0}(\boldsymbol{r})$ 的泛函。
定理 3.2 对于可由某外势的基态表示的试探密度，可以定义与外势无关的普适泛函。若扩大定义域至 $N$ 可表象密度，可采用受约束搜索 $F[\rho]=\min_{\Psi\to\rho}\langle\Psi|\hat T+\hat V_{\mathrm{ee}}|\Psi\rangle$。若给定 $V_{\mathrm{ext}}(\boldsymbol{r})$，仅当电子密度分布 $\rho(\boldsymbol{r})$ 为该体系的基态电子密度分布 $\rho^{0}(\boldsymbol{r})$ 时，泛函 $E[\tilde{\rho}(\boldsymbol{r})]$ 最小，且给出体系的基态能量。
首先证明定理 3.1。用反证法，设有外势差别不仅仅为一个常数的两个外势场（外势分别为 $V_{\mathrm{ext}}^{0}(\boldsymbol{r})$ 和 $V_{\mathrm{ext}}^{1}(\boldsymbol{r})$），它们所对应的基态电子密度分布为 $\rho^{0}(\boldsymbol{r})$。这两个不同的外势场规定了体系的两个哈密顿量 $H^{0}$ 和 $H^{1}$，相应地有两个不同的基态波函数 $\psi^{0}$ 以及 $\psi^{1}$。当基态非简并时（简并情况可以通过引入一个微弱势场加以消除），有
$$ \begin{aligned} E^{0}&=\langle\psi^{0}|H^{0}|\psi^{0}\rangle \lt \langle\psi^{1}|H^{0}|\psi^{1}\rangle\\ &=E^{1}+\int[V_{\mathrm{ext}}^{0}(\boldsymbol r)-V_{\mathrm{ext}}^{1}(\boldsymbol r)]\rho^{1}(\boldsymbol r)\,\mathrm d\boldsymbol r \end{aligned} \tag{3.164} $$
同理，有
$$ \begin{aligned} E^{1}&=\langle\psi^{1}\,|\,H^{1}\,|\,\psi^{1}\rangle\lt \langle\psi^{0}\,|\,H^{1}\,|\,\psi^{0}\rangle=\langle\psi^{0}\,|\,H^{0}\,|\,\psi^{0}\rangle+\langle\psi^{0}\,|\,H^{1}-H^{0}\,|\,\psi^{0}\rangle\\ &=E^{0}+\int[V_{\mathrm{ext}}^{1}(\boldsymbol{r})-V_{\mathrm{ext}}^{0}(\boldsymbol{r})]\rho^{0}(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r} \end{aligned} \tag{3.165} $$
将式（3.164）和式（3.165）相加，则有不等式 $E^{0}+E^{1}\lt E^{0}+E^{1}$，这显然不成立。因此假设错误，故 $V_{\mathrm{ext}}^{0}(\boldsymbol{r})$ 与 $V_{\mathrm{ext}}^{1}(\boldsymbol{r})$ 只可能相差一个常数。定理 3.1 得证。
再证明定理 3.2。根据定理 3.1 及推论 3.1，多体波函数 $\psi^{0}$ 是电子密度分布 $\rho^{0}(\boldsymbol{r})$ 的泛函。而波函数的确定意味着体系的所有性质，如动能、电子间相互作用等均可确定。因此，可以认为体系动能和电子间相互作用也可以表示为 $\rho^{0}(\boldsymbol{r})$ 的泛函，分别记为 $T[\rho(\boldsymbol{r})]$ 和 $E_{\mathrm{ee}}[\rho(\boldsymbol{r})]$。 这里省去基态上标时，须将密度的定义域限定为可由某外势基态表示的密度；若采用受约束搜索，则可推广至 $N$ 可表象密度。由此，可以设泛函 $F[\rho(\boldsymbol{r})]$ 为
$$ F[\rho(\boldsymbol{r})]=T[\rho(\boldsymbol{r})]+E_{\mathrm{ee}}[\rho(\boldsymbol{r})]=\langle\psi\,|\,\hat{T}+\hat{V}_{\mathrm{ee}}\,|\,\psi\rangle \tag{3.166} $$
这是泛函的普遍形式。该泛函仅取决于 $\rho(\boldsymbol{r})$，而与外势 $V_{\mathrm{ext}}(\boldsymbol{r})$ 无关。
对于任意 $V_{\mathrm{ext}}(\boldsymbol{r})$，可以定义 Hohenberg-Kohn 能量泛函 $E^{\mathrm{HK}}[\rho(\boldsymbol{r}),V_{\mathrm{ext}}(\boldsymbol{r})]$ 为
$$ E^{\mathrm{HK}}[\rho(\boldsymbol{r}),V_{\mathrm{ext}}(\boldsymbol{r})]=T[\rho(\boldsymbol{r})]+E_{\mathrm{ee}}[\rho(\boldsymbol{r})]+\int V_{\mathrm{ext}}(\boldsymbol{r})\rho(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r}+E_{\mathrm{II}}(\{\boldsymbol{R}_{\mathrm{I}}\}) \tag{3.167} $$
若已给定 $V_{\mathrm{ext}}^{0}(\boldsymbol{r})$，其相应的基态电子密度为 $\rho^{0}(\boldsymbol{r})$，则 Hohenberg-Kohn 能量泛函等于哈密顿量 $H^{0}$ 对基态多体波函数 $\psi^{0}(\boldsymbol{r})$ 的期待值，即
$$ E^{\mathrm{HK}}[\rho^{0},V_{\mathrm{ext}}^{0}]=\langle\psi^{0}\,|\,H^{0}\,|\,\psi^{0}\rangle \tag{3.168} $$
也即体系的基态能量。设对于另一个电子密度 $\rho'(\boldsymbol{r})$（外势为 $V'(\boldsymbol{r})$ 的外势场的基态电子密度），相应地有基态多体波函数 $\psi'(\boldsymbol{r})$。类似于定理 3.1 的证明过程，有
$$ E'=\langle\psi'\,|\,H^{0}\,|\,\psi'\rangle\gt \langle\psi^{0}\,|\,H^{0}\,|\,\psi^{0}\rangle \tag{3.169} $$
因此，$E^{\mathrm{HK}}[\rho^{0},V_{\mathrm{ext}}^{0}]$ 的最小值仅在电子密度（相对于 $V^{0}(\boldsymbol{r})$）为基态电子密度 $\rho^{0}(\boldsymbol{r})$ 时才能取得。
Hohenberg–Kohn 变分原理表明，若普适泛函 $F[\rho]$ 已知，可通过对允许的电子密度求能量最小值而获得基态能量。这明显比 Hartree-Fock 方法来得简单。但是到目前为止，对于相互作用电子气，人们还不知道 $F[\rho(\boldsymbol{r})]$ 的具体形式。
Kohn-Sham 方程
Hohenberg-Kohn 定理从理论上保证了体系基态能量泛函的存在性与唯一性，且指出了该能量泛函对电荷密度 $\rho(\boldsymbol{r})$ 求变分达到的极值即为体系基态能量。但是如 3.2.2 节所述，与外势场无关的普适泛函 $F[\rho(\boldsymbol{r})]$ 的具体形式未知，因此 Hohenberg-Kohn 定理不能直接用于求解问题。实际运用密度泛函理论时，需要用到 Kohn-Sham（KS）方程\cite{kohn1965self}。将电子密度表示为
$$ \rho(\boldsymbol{r})=\sum_{j=1}^{N}\phi_j^{*}(\boldsymbol{r})\phi_j(\boldsymbol{r}) \tag{3.170} $$
$\{\phi_j\}$ 为互为正交的一组波函数（共 $N$ 个）。Kohn–Sham 构造引入一个无相互作用辅助体系，并要求它的基态密度 $\rho^{0}(\boldsymbol r)$ 等于相互作用体系的目标密度。下述轨道方程的建立需假定所考虑的密度可由合适的无相互作用局域势表示；Hohenberg–Kohn 定理本身并不保证任意试探密度都满足这个条件。辅助体系的占据轨道给出
$$ \rho^{0}(\boldsymbol{r})=\sum_{j}\phi_j^{*}(\boldsymbol{r})\phi_j(\boldsymbol{r}) $$
因为 $\rho(\boldsymbol{r})$ 与 $\rho^{0}(\boldsymbol{r})$ 相等，所以方程（3.170）总是成立的。上述讨论也称为 Kohn-Sham 拟设。
具体写出外势的表达式，即
$$ V_{\mathrm{ext}}(\boldsymbol{r})=\sum_{l}V(\boldsymbol{r}-\boldsymbol{R}_{\mathrm{I}}) \tag{3.171} $$
并且引入两个已知的泛函：无相互作用电子气的动能泛函 $T_0[\rho]$ 和电子间库仑相互作用（又称 Hartree 项）$E_{\mathrm{H}}[\rho]$。采用原子单位制 ${\hbar=m_{\mathrm e}=e=4\pi\varepsilon_0=1}$，则有
$$ \begin{gathered} T_0[\rho]=\sum_j\langle\phi_j|{-\tfrac12\nabla^2}|\phi_j\rangle,\\ E_{\mathrm H}[\rho]=\frac12\iint\frac{\rho(\boldsymbol r)\rho(\boldsymbol r')}{|\boldsymbol r-\boldsymbol r'|}\,\mathrm d\boldsymbol r\,\mathrm d\boldsymbol r' =\frac12\sum_{ij}\langle\phi_i\phi_j|r_{12}^{-1}|\phi_i\phi_j\rangle \end{gathered} \tag{3.172} $$
$T_0[\rho]+E_{\mathrm{H}}[\rho]$ 与方程（3.166）显然并不一致，两者之间的差别可以归结为描述多体相互作用的交换关联泛函 $E_{\mathrm{xc}}[\rho]$，且
$$ E_{\mathrm{xc}}[\rho]=T[\rho]-T_0[\rho]+E_{\mathrm{ee}}[\rho]-E_{\mathrm{H}}[\rho] \tag{3.173} $$
因此，借助式（3.172）及式（3.173），方程（3.167）可以重新写为
$$ \begin{aligned} E^{\mathrm{HK}}[\rho,V_{\mathrm{ext}}]={}&\sum_j\langle\phi_j|-\tfrac12\nabla^2+V_{\mathrm{ext}}|\phi_j\rangle +\frac12\sum_{i,j}\langle\phi_i\phi_j|r_{12}^{-1}|\phi_i\phi_j\rangle\\ &+E_{\mathrm{xc}}[\rho]+E_{\mathrm{II}}. \end{aligned} \tag{3.174} $$
引入所谓的交换关联势 $V_{\mathrm{xc}}$，即
$$ V_{\mathrm{xc}}=\frac{\delta E_{\mathrm{xc}}}{\delta\rho} \tag{3.175} $$
可以将 $E_{\mathrm{xc}}[\rho]$ 随 $\rho$ 的小量变化写为
$$ \delta E_{\mathrm{xc}}=\int V_{\mathrm{xc}}(\boldsymbol r)\,\delta\rho(\boldsymbol r)\,\mathrm d\boldsymbol r,\qquad \delta\rho=\sum_j\bigl(\delta\phi_j^*\phi_j+\phi_j^*\delta\phi_j\bigr) \tag{3.176} $$
将式（3.176）代入方程（3.174），考虑到占据轨道必须两两正交且归一，满足约束条件
$$ \langle\phi_i|\phi_j\rangle=\delta_{ij}\qquad(i,j=1,\ldots,N) \tag{3.177} $$
则利用 1.3.7 节中介绍的带限制条件的拉格朗日乘子法，对方程（3.174）关于 $\langle\phi_j\,|$ 求变分极值，记相应的约束泛函为 $L_{\mathrm{KS}}$，可以得到
$$ {L_{\mathrm{KS}}}=E^{\mathrm{HK}}[\{\phi_i\}]-\sum_{i,j}\Lambda_{ij} \bigl(\langle\phi_i|\phi_j\rangle-\delta_{ij}\bigr), \qquad \Lambda_{ij}=\Lambda_{ji}^{*}. $$
$$ \begin{aligned} \frac{\delta{L_{\mathrm{KS}}}}{\delta\phi_j^*(\boldsymbol r)} &=\left[-\frac{\nabla^2}{2}+V_{\mathrm{ext}}(\boldsymbol r) +V_{\mathrm H}(\boldsymbol r)+V_{\mathrm{xc}}(\boldsymbol r)\right]\phi_j(\boldsymbol r) -\sum_i\Lambda_{ji}\phi_i(\boldsymbol r)=0,\\ V_{\mathrm H}(\boldsymbol r)&=\int\frac{\rho(\boldsymbol r')}{|\boldsymbol r-\boldsymbol r'|}\,\mathrm d\boldsymbol r'. \end{aligned} \tag{3.178} $$
通过对厄米乘子矩阵 $\boldsymbol\Lambda$ 作幺正对角化，得到正则轨道；式（3.178）于是化为单轨道方程。
对式（3.178）重新进行整理，可得
$$ \left(-\frac{\boldsymbol{\nabla}^{2}}{2}+V_{\mathrm{KS}}(\boldsymbol{r})\right)\phi_j(\boldsymbol{r})=\varepsilon_j\phi_j(\boldsymbol{r}) \tag{3.179} $$
式中
$$ V_{\mathrm{KS}}(\boldsymbol{r})=V_{\mathrm{ext}}(\boldsymbol{r})+V_{\mathrm{H}}(\boldsymbol{r})+V_{\mathrm{xc}}(\boldsymbol{r}) \tag{3.180} $$
方程（3.179）即著名的 Kohn-Sham 方程。假设已经得到一组 Kohn-Sham 方程的本征值 $\{\varepsilon_j\}$，则可以将体系基态总能表示为
$$ E_0=\sum_{j}\varepsilon_j-\frac{1}{2}\iint\frac{\rho(\boldsymbol{r}')\rho(\boldsymbol{r})}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}'+E_{\mathrm{xc}}[\rho(\boldsymbol{r})]-\int V_{\mathrm{xc}}(\boldsymbol{r})\rho(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r}+E_{\mathrm{II}} \tag{3.181} $$
方程（3.181）右端第一项 $\displaystyle\sum_{j}\varepsilon_j$ 称为能带结构能（band structure energy），而后三项称为冗余项（double counting，d. c.）。
交换关联能概述
Kohn-Sham 方程最重要的一个特点就是将所有未知的、难以求得的多体项的贡献全部包含在交换关联能 $E_{\mathrm{xc}}$ 中了，所以这一项的精确度直接决定了 Kohn-Sham 方程的计算精度。由 3.2.3 节的讨论可知，在严格意义上，DFT 理论中的交换关联能 $E_{\mathrm{xc}}$ 与 Hartree-Fock 近似下对应的项 $E_{\mathrm{xc}}^{\mathrm{HF}}$ 并不相同，因为前者还包含了相互作用电子气的动能泛函的修正。 由泛函的定义可知，精确交换关联能包含相互作用动能修正、交换能及库仑关联能：
$$ E_{\mathrm{xc}}[\rho]=\bigl(T[\rho]-T_0[\rho]\bigr)+\bigl(E_{\mathrm{ee}}[\rho]-E_{\mathrm H}[\rho]\bigr) \tag{3.182} $$
为了将动能泛函的修正包括进去，一般采用耦合常数积分法。设一个电子密度为 $\rho(\boldsymbol{r})$ 的体系，电子间的相互作用 $E_{\mathrm{ee}}$ 正比于 $e^{2}$。现在假设该体系的电子所带电荷可以在 $[0,1]$ 之间变化，设为 $\sqrt{\lambda}e$，则 $E_{\mathrm{ee}}\propto\lambda e^{2}$，$\lambda$ 称为耦合常数。显然，若 $\lambda=0$，则体系无相互作用电子气，若 $\lambda=1$，则体系为真实的物理体系。参照式（3.57），可以将 $E_{\mathrm{xc}}$ 表示为
$$ E_{\mathrm{xc}}=\frac{1}{2}\int\rho(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r}\int\frac{\bar{\rho}_{\mathrm{xc}}(\boldsymbol{r},\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}' \tag{3.183} $$
式中
$$ \bar{\rho}_{\mathrm{xc}}(\boldsymbol{r},\boldsymbol{r}')=\int_{0}^{1}\rho_{\mathrm{xc}}(\boldsymbol{r},\boldsymbol{r}',\lambda)\,\mathrm{d}\lambda=\rho(\boldsymbol{r}')\left[\int_{0}^{1}g(\boldsymbol{r},\boldsymbol{r}',\lambda)\,\mathrm{d}\lambda-1\right]=\rho(\boldsymbol{r}')[\bar{g}(\boldsymbol{r},\boldsymbol{r}')-1] \tag{3.184} $$
在均匀电子气模型中，$\rho_{\mathrm{xc}}$ 仅是电子间距 $r=|\,\boldsymbol{r}-\boldsymbol{r}'\,|$ 的函数。到目前为止，虽然给出了 $E_{\mathrm{xc}}$ 的一个合理的近似，但是仍然需要知道 $\rho_{\mathrm{xc}}(r,\lambda)$ 随耦合常数 $\lambda$ 的变化关系才能进行定量计算。为了具体给出定量计算的方法，需要讨论均匀电子气模型下的方程（3.2）。首先利用玻恩-奥本海默近似忽略原子核的动能项，然后设正电荷也以同样的密度 $\rho_0$ 在空间均匀分布，则方程（3.2）中的哈密顿量为（原子单位制下）
$$ \hat{H}=-\sum_{i}\frac{\boldsymbol{\nabla}_i^{2}}{2}+\frac{1}{2}\sum_{i\neq j}\frac{1}{|\,\boldsymbol{r}_i-\boldsymbol{r}_j\,|}-\frac{1}{2}\iint\frac{\rho_0^{2}}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}' \tag{3.185} $$
现在将位置坐标 $\boldsymbol{r}$ 以电子密度参数 $r_{\mathrm{s}}a_0$（见式（3.92））为单位进行约化（$\tilde{\boldsymbol{r}}=\boldsymbol{r}/(r_{\mathrm{s}}a_0)$），则式（3.185）右端的三项分别为
$$ \begin{gathered} \sum_{i}\frac{\boldsymbol{\nabla}_i^{2}}{2}=\left(\frac{1}{r_{\mathrm{s}}a_0}\right)^{2}\sum_{i}\frac{\tilde{\boldsymbol{\nabla}}_i^{2}}{2}\\ \frac{1}{2}\sum_{i\neq j}\frac{1}{|\,\boldsymbol{r}_i-\boldsymbol{r}_j\,|}=\frac{1}{2r_{\mathrm{s}}a_0}\sum_{i\neq j}\frac{1}{\tilde{\boldsymbol{r}}_i-\tilde{\boldsymbol{r}}_j}\\ \frac{1}{2}\iint\frac{\rho_0^{2}}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}'=\frac{1}{2}\sum_{i}\frac{1}{\rho_0}\int\frac{\rho_0^{2}}{|\,\boldsymbol{r}-\boldsymbol{r}_i\,|}\,\mathrm{d}\boldsymbol{r}=\frac{1}{2r_{\mathrm{s}}a_0}\frac{3}{4\pi}\sum_{i}\int\frac{1}{|\,\tilde{\boldsymbol{r}}-\tilde{\boldsymbol{r}}_i\,|}\,\mathrm{d}\tilde{\boldsymbol{r}} \end{gathered} $$
由此可得
$$ \hat{H}=\left(\frac{1}{r_{\mathrm{s}}a_0}\right)^{2}\sum_{i}\left[-\frac{\tilde{\boldsymbol{\nabla}}_i^{2}}{2}+\frac{1}{2}r_{\mathrm{s}}a_0\left(\sum_{i\neq j}\frac{1}{|\,\tilde{\boldsymbol{r}}_i-\tilde{\boldsymbol{r}}_j\,|}-\frac{3}{4\pi}\int\frac{1}{|\,\tilde{\boldsymbol{r}}-\tilde{\boldsymbol{r}}_i\,|}\,\mathrm{d}\tilde{\boldsymbol{r}}\right)\right] \tag{3.186} $$
式（3.186）表明，耦合常数可以用 $r_{\mathrm{s}}a_0$，即电子密度来表示。因此，原则上可以通过模拟不同密度下的均匀电子气\cite{gorigiorgi2000analytic}，得到相应的 $\rho_{\mathrm{xc}}(\boldsymbol{r},\boldsymbol{r}',\lambda)$ 或者 $g_{\mathrm{xc}}(\boldsymbol{r},\boldsymbol{r}',\lambda)$（二者通过式（3.61）相互联系），从而求得 $E_{\mathrm{xc}}$。对于电子非均匀分布的体系，如原子、分子、固体等，情况显然更为复杂。为了使问题可解，通常需要预先做某种假设，以便尽可能地利用均匀电子气的解析或模拟结果。我们将在下几节中做具体介绍。
联系 3.1.5 节以及 3.1.6 节中的讨论可知，原则上体系的交换能是可以通过解析形式给出的，但是采用这种做法计算量过大，且精度不会显著提高（因为关联能并没有对应的精确解，所以即使其他各项均得到精确解，关联能部分的误差也无法被部分抵消）。因此，在实际应用中，研究人员通常寻求合适的近似方法来处理交换关联能。
局域密度近似
与理想化的均匀电子气模型不同，实际体系中的电荷分布往往呈现出非常明显的起伏及各向异性。为了使用均匀电子气的结果，最简单的办法是将 $E_{\mathrm{xc}}[\rho]$ 的求解视为各个离散的 $\boldsymbol{r}$ 终点处仅由局域电荷密度 $\rho(\boldsymbol{r})$ 决定的交换关联能密度 $\varepsilon_{\mathrm{xc}}[\rho(\boldsymbol{r})]$ 的加权求和，而权重就是 $\rho(\boldsymbol{r})$，也即
$$ E_{\mathrm{xc}}[\rho]=\int\varepsilon_{\mathrm{xc}}[\rho(\boldsymbol{r}),\boldsymbol{r}]\rho(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r} \tag{3.187} $$
这种处理方法称为局域密度近似（local density approximation，LDA）。LDA 中 $\varepsilon_{\mathrm{xc}}[\rho]$ 分为交换能密度 $\varepsilon_{\mathrm{x}}[\rho]$ 和关联能密度 $\varepsilon_{\mathrm{c}}[\rho]$ 两部分。$\varepsilon_{\mathrm{x}}[\rho]$ 一般采用均匀电子气结果（见式（3.109））给出。而 $\varepsilon_{\mathrm{c}}[\rho]$ 则通常没有严格的解析解，其表达式主要基于 Ceperley 和 Alder 对均匀电子气的量子蒙特卡罗（QMC）模拟\cite{ceperley1980ground}。在实际应用中，为了避免大量的计算，通常用拟合的函数形式来近似 $\varepsilon_{\mathrm{c}}$，其中最常见的几种形式（能量单位均为 Hartree）如下。
1. Perdew-Zunger（PZ）函数\cite{perdew1981self}
$$ \varepsilon_{\mathrm{c}}^{\mathrm{PZ}}(r_{\mathrm{s}})=\begin{cases}A\ln r_{\mathrm{s}}+B+Cr_{\mathrm{s}}\ln r_{\mathrm{s}}+Dr_{\mathrm{s}}&(r_{\mathrm{s}}\leqslant1)\\\gamma/(1+\beta_1\sqrt{r_{\mathrm{s}}}+\beta_2r_{\mathrm{s}})&(r_{\mathrm{s}}\gt 1)\end{cases} \tag{3.188} $$
式中：$A=0.0311$，$B=-0.048$，$C=0.002$，$D=-0.0116$。其中 $A$ 与 $B$ 与方程（3.112）的 PZ 高密度分支一致。$r_{\mathrm{s}}\leqslant1$ 是高密度极限，不考虑自旋极化的情况已经在 3.1.11 节中的关联能部分讨论过了；而电子气密度较低时，$r_{\mathrm{s}}\gt 1$，在 Perdew-Zunger 函数中，$\gamma=-0.1423$，$\beta_1=1.0529$，$\beta_2=0.3334$。
2. Vosko-Wilk-Nusair（VWN）函数\cite{vosko1980influence,vosko1980accurate}
$$ \begin{aligned} \varepsilon_{\mathrm c}^{\mathrm{VWN}}(r_{\mathrm s})=A\Bigg\{& \ln\!\frac{x^2}{X(x)}+\frac{2b}{Q}\arctan\!\frac{Q}{2x+b}\\ &-\frac{bx_0}{X(x_0)}\left[ \ln\!\frac{(x-x_0)^2}{X(x)} +\frac{2(b+2x_0)}{Q}\arctan\!\frac{Q}{2x+b} \right]\Bigg\},\\ x&=\sqrt{r_{\mathrm s}},\qquad X(x)=x^2+bx+c,\qquad Q=\sqrt{4c-b^2}. \end{aligned} \tag{3.189} $$
对于自旋非极化的情况，有 $A=0.0310907$、$x_0=-0.10498$、$b=3.72744$、$c=12.9352$。
还有另外一些 LDA 下的函数形式，这里不多做介绍，具体函数形式请参考文献\cite{perdew1992accurate}。
自旋极化情况
首先定义自旋极化分布：
$$ \zeta=\frac{\rho^{\uparrow}(\boldsymbol{r})-\rho^{\downarrow}(\boldsymbol{r})}{\rho^{\uparrow}(\boldsymbol{r})+\rho^{\downarrow}(\boldsymbol{r})} \tag{3.190} $$
对于自旋极化情况，一般将交换关联能密度表示为完全非极化情况（$\zeta=0$）与完全极化情况（$\zeta=1$）的插值。Barth 和 Hedin 指出，对于交换能密度，可以采取如下方式\cite{vonbarth1972local}：
$$ \varepsilon_{\mathrm{x}}(r_{\mathrm{s}},\zeta)=\varepsilon_{\mathrm{x}}^{\mathrm{U}}+(\varepsilon_{\mathrm{x}}^{\mathrm{P}}-\varepsilon_{\mathrm{x}}^{\mathrm{U}})f(\zeta) \tag{3.191} $$
其中上标 U 和 P 分别代表 $\zeta=0$ 和 $\zeta=1$ 的情况，而 $f(\zeta)$ 为
$$ f(\zeta)=\frac{(1+\zeta)^{4/3}+(1-\zeta)^{4/3}-2}{2^{4/3}-2} \tag{3.192} $$
Perdew 与 Zunger 给出了完全自旋极化均匀电子气关联能的分段参数化\cite{perdew1981self}：
$$ \varepsilon_{\mathrm{c}}^{\mathrm{P}}(r_{\mathrm{s}})=\begin{cases}0.01555\ln r_{\mathrm{s}}-0.0269+0.0007r_{\mathrm{s}}\ln r_{\mathrm{s}}-0.0048r_{\mathrm{s}},&r_{\mathrm{s}}\leqslant1\\-0.0843/(1+1.3981\sqrt{r_{\mathrm{s}}}+0.2611r_{\mathrm{s}}),&r_{\mathrm{s}}\gt 1\end{cases} \tag{3.193} $$
对于部分自旋极化的情况，Vosko、Wilk 和 Nusair 给出了包含自旋刚度的插值形式\cite{vosko1980influence,vosko1980accurate,kohanoff2006electronic}：
$$ \varepsilon_{\mathrm{c}}(r_{\mathrm{s}},\zeta)=\varepsilon_{\mathrm{c}}^{\mathrm{U}}+\left[\frac{f(\zeta)}{f''(0)}\right](1-\zeta^{4})\alpha_{\mathrm{c}}(r_{\mathrm{s}})+f(\zeta)\zeta^{4}[\varepsilon_{\mathrm{c}}^{\mathrm{P}}-\varepsilon_{\mathrm{c}}^{\mathrm{U}}] \tag{3.194} $$
其中 $\varepsilon_{\mathrm{c}}^{\mathrm{P}}$ 中的参数分别为
$$ A^{\mathrm{P}}=0.01553535,\quad x_0^{\mathrm{P}}=-0.325,\quad b^{\mathrm{P}}=7.06042,\quad c^{\mathrm{P}}=18.0578 $$
而 $\alpha_{\mathrm{c}}(r_{\mathrm{s}})$ 也取方程（3.189）所示的形式，其中四个参数分别为
$$ A^{\alpha}=-1/(6\pi^{2}),\quad x_0^{\alpha}=-0.0047584,\quad b^{\alpha}=1.13107,\quad c^{\alpha}=13.0045 $$
最后给出局域自旋密度近似（local spin density approximation，LSDA）下的交换关联能 $E_{\mathrm{xc}}^{\mathrm{LSDA}}$ 的表达式：
$$ E_{\mathrm{xc}}^{\mathrm{LSDA}}[\rho^{\uparrow},\rho^{\downarrow}]=\int\rho(\boldsymbol{r})[\varepsilon_{\mathrm{x}}(r_{\mathrm{s}},\zeta)+\varepsilon_{\mathrm{c}}(r_{\mathrm{s}},\zeta)]\,\mathrm{d}\boldsymbol{r} \tag{3.195} $$
LDA 下的交换关联势
LDA 下的交换关联势有非常简单的形式。由 $V_{\mathrm{xc}}$ 的定义式（3.175）及 LDA 下 $E_{\mathrm{xc}}$ 的表达式（式（3.187）），可得
$$ V_{\mathrm{xc}}(\boldsymbol{r})=\frac{\delta E_{\mathrm{xc}}[\rho(\boldsymbol{r})]}{\delta\rho(\boldsymbol{r})}=\varepsilon_{\mathrm{xc}}[\rho,\boldsymbol{r}]+\rho(\boldsymbol{r})\frac{\delta\varepsilon_{\mathrm{xc}}[\rho,\boldsymbol{r}]}{\delta\rho(\boldsymbol{r})} \tag{3.196} $$
与关于 $\varepsilon_{\mathrm{xc}}$ 的讨论类似，一般也将 $V_{\mathrm{xc}}(\boldsymbol{r})$ 分为 $V_{\mathrm{x}}(\boldsymbol{r})$ 和 $V_{\mathrm{c}}(\boldsymbol{r})$ 两部分。由前面几节的讨论可知，实际计算中需要将 $V_{\mathrm{xc}}$ 表示成 $r_{\mathrm{s}}$ 的函数。其中交换能部分 $V_{\mathrm{x}}(r_{\mathrm{s}})$ 比较简单，根据方程（3.109），有
$$ V_{\mathrm{x}}=\frac{4}{3}\varepsilon_{\mathrm{x}}(r_{\mathrm{s}})\propto\rho^{1/3} \tag{3.197} $$
与 $X_\alpha$ 方法的表达式（3.143）等价。
而由式（3.92）及式（3.109）可以直接得到 $V_{\mathrm{c}}(r_{\mathrm{s}})$ 的表达式：
$$ V_{\mathrm{c}}(r_{\mathrm{s}})=\varepsilon_{\mathrm{c}}-\frac{r_{\mathrm{s}}}{3}\frac{\mathrm{d}\varepsilon_{\mathrm{c}}}{\mathrm{d}r_{\mathrm{s}}} \tag{3.198} $$
由此不难得到 PZ 函数及 VWN 函数形式的交换关联势。
考虑自旋自变量 $\sigma$ 的情况，直接给出 $V_{\mathrm{xc}}$ 的表达式：
$$ V_{\mathrm{xc}}^{\sigma}(\boldsymbol r)=\frac{\delta E_{\mathrm{xc}}[\rho^\uparrow,\rho^\downarrow]}{\delta\rho^\sigma(\boldsymbol r)} =\varepsilon_{\mathrm{xc}}(\rho^\uparrow,\rho^\downarrow) +\rho(\boldsymbol r)\frac{\partial\varepsilon_{\mathrm{xc}}(\rho^\uparrow,\rho^\downarrow)}{\partial\rho^\sigma},\qquad\rho=\rho^\uparrow+\rho^\downarrow \tag{3.199} $$
LDA 的特性概述
应该指出的是，交换能本质上是一个非局域的函数。也就是说，其泛函值取决于全空间的电子密度，而非取决于空间某点的局域电子密度。这一点通过比较方程（3.183）和方程（3.187）就可看出：
$$ \varepsilon_{\mathrm{xc}}[\rho]=\frac{1}{2}\int\frac{\bar{\rho}_{\mathrm{xc}}(\boldsymbol{r},\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}' \tag{3.200} $$
均匀电子气只是一个可以解析求解的特例。对于电子气分布较均匀的体系，如简单金属等，LDA 显然是非常合理的。但是对于以共价键为主的晶体或分子体系，LDA 的效果往往要差一些。除了前面已经讨论过的假设之外，LDA 方法还定义交换关联空穴 $\bar{\rho}_{\mathrm{xc}}^{\mathrm{LDA}}$ 为
$$ \bar\rho_{\mathrm{xc}}^{\mathrm{LDA}}(\boldsymbol r,\boldsymbol r') =\rho(\boldsymbol r)\left[\bar g^{\mathrm{hom}}(|\boldsymbol r-\boldsymbol r'|;\rho(\boldsymbol r))-1\right] \tag{3.201} $$
与严格的定义方程（3.184）相比，LDA 采用 $\boldsymbol{r}$ 终点处而非 $\boldsymbol{r}'$ 终点处的电子密度，而且电子对关联函数 $g(\boldsymbol{r},\boldsymbol{r}')$ 采用了密度为 $\rho(\boldsymbol{r})$ 的均匀电子气的结果。这表明在 LDA 下，交换关联能是 $\boldsymbol{r}$ 终点处的电子密度与 $\boldsymbol{r}$ 终点处的交换关联空穴之间的局域相互作用，而且由式（3.201）可得
$$ \int\bar\rho_{\mathrm{xc}}^{\mathrm{LDA}}(\boldsymbol r,\boldsymbol r')\,\mathrm d\boldsymbol r'=-1 \tag{3.202} $$
因为在每一个 $\boldsymbol{r}$ 终点处，关联函数 $\bar{g}^{\mathrm{hom}}$ 都对应着一个密度为 $\rho(\boldsymbol{r})$ 的均匀电子气体系。这样，式（3.202）中的积分实际上就是对均匀电子气的交换关联空穴的积分，因此 LDA 满足交换关联空穴的求和要求。因为只考虑局域的电荷密度信息，所以 LDA 是一个比较粗略的近似，但是其在实际使用中却取得了非常好的效果，式（3.202）能够成立是其中一个很重要的原因。
但是，在实际使用 LDA 过程中会出现以下问题：
（1）LDA 在计算中往往会高估结合能、低估晶格常数，因此用 LDA 计算的弹性常数往往比实验值高大约 10%。
（2）因为无法正确处理电子跃迁时产生的交换关联势的突变，所以 LDA 会低估半导体或绝缘体的带隙及介电常数。
（3）因为没有考虑 $E_{\mathrm{xc}}[\rho]$ 的非局域效应，所以无法有效处理范德瓦尔斯力。
（4）对于强关联体系，如过渡金属氧化物等，LDA 无法获得令人满意的结果。
尽管 LDA 在很多情况下可以提供合理的近似，但在一些特定场景中，它的局限性仍然较为明显。为了解决这些问题，研究人员在不断探索更加精确和可靠的交换关联势近似方法。
广义梯度近似
在实际的固体体系中，电子云在晶体中的分布并不均匀，因此对 LDA 的一个自然改进便是在交换关联项中引入电子密度的梯度以及更高阶的导数项。然而，由于实际体系中 $|\,\boldsymbol{\nabla}\rho\,|$ 通常较大，且最初提出的泛函形式并不满足 $\varepsilon_{\mathrm{xc}}$ 准则，因此早期的尝试并未取得成功。经过不断地尝试和改进，研究人员提出了一系列方案，能够很好地处理梯度项并尽可能地满足上述准则。这些方案统称为广义梯度近似（generalized gradient approximation，GGA）。
考虑电子密度梯度的修正，可以将 $E_{\mathrm{xc}}$ 表示为
$$ E_{\mathrm{xc}}[\rho]=\int\mathrm{d}\boldsymbol{r}\,\rho(\boldsymbol{r})\varepsilon_{\mathrm{xc}}^{\mathrm{hom}}F_{\mathrm{xc}}[\rho^{\uparrow}(\boldsymbol{r}),\rho^{\downarrow}(\boldsymbol{r}),|\,\boldsymbol{\nabla}\rho^{\uparrow}(\boldsymbol{r})\,|,|\,\boldsymbol{\nabla}\rho^{\downarrow}(\boldsymbol{r})\,|,\cdots] \tag{3.203} $$
式中：$F_{\mathrm{xc}}$ 称为增效函数，包含非局域、非均匀项对均匀电子气结果的修正。
对于交换能，因为其只存在于相同自旋态中，所以可以分解成两项，即
$$ E_{\mathrm{x}}[\rho^{\uparrow},\rho^{\downarrow}]=\frac{1}{2}(E_{\mathrm{x}}[2\rho^{\uparrow}]+E_{\mathrm{x}}[2\rho^{\downarrow}]) \tag{3.204} $$
式（3.204）表示自旋标度关系，由 Oliver 和 Perdew 提出\cite{oliver1979spin}。该方程右端的两项均为非自旋极化的结果，因此可以将 $F_{\mathrm{x}}$ 表示为总电子密度和总电子密度各阶导数的函数。为方便起见，引入无量纲约化梯度 $s_1$ 及带符号的约化拉普拉斯量 $s_2$，有
$$ s_1=\frac{|\boldsymbol\nabla\rho|}{2k_{\mathrm F}\rho},\qquad s_2=\frac{\nabla^2\rho}{(2k_{\mathrm F})^2\rho},\qquad k_{\mathrm F}=(3\pi^2\rho)^{1/3} \tag{3.205} $$
GGA 的理论推导均涉及多体理论以及复杂的公式推导，在这里不拟详细讨论，具体可参考文献\cite{ma1968correlation,kleinman1988gradient,perdew1996comparison}。$F_{\mathrm{x}}(s)$ 做泰勒展开精确到 $O(\boldsymbol{\nabla}^{6}\rho)$ 的普遍表达式为\cite{kohanoff2006electronic,svendsen1996gradient}
$$ F_{\mathrm{x}}(s)=1+\frac{10}{81}s_1^{2}+\frac{146}{2025}s_2^{2}-\frac{73}{405}s_1^{2}s_2+Ds_1^{4}+O(\boldsymbol{\nabla}^{6}\rho) \tag{3.206} $$
式中四阶梯度系数 $D$ 依赖所选的交换能量密度规范，不能普遍设为零；这一展开也不等同于任意 GGA 泛函的定义。
关联能的计算更加困难，到目前为止还没有人提出比较普遍的表达式。因此一般的做法是尽量使 $\varepsilon_{\mathrm{c}}$ 满足交换关联能的准则，且在 $s\to0$ 的情况下回归到 LDA 的形式。
GGA 泛函
目前 DFT 计算中最为常见的 GGA 泛函包括 Becke-Lee-Yang-Parr（BLYP）、Perdew-Wang（PW91）和 Perdew-Burke-Ernzerhof（PBE）泛函（能量单位均为 Hartree）。
1. BLYP\cite{becke1988density,lee1988development} 泛函
BLYP 组合了 Becke 1988 交换与 LYP 关联；对于自旋非极化密度，前者的每电子交换能写为：
$$ \varepsilon_{\mathrm{x}}^{\mathrm{B88}}=\varepsilon_{\mathrm{x}}^{\mathrm{LDA}} \left(1+\frac{\beta}{2^{1/3}A_x}\frac{x^{2}}{1+6\beta x\sinh^{-1}x}\right) \tag{3.207} $$
式中：$\beta=0.0042$，$A_x=(3/4)(3/\pi)^{1/3}$，$x=2(6\pi^{2})^{1/3}s_1=2^{1/3}|\,\boldsymbol{\nabla}\rho(\boldsymbol{r})\,|/\rho(\boldsymbol{r})^{4/3}$。
$$ e_{\mathrm c}^{\mathrm{LYP}}(\boldsymbol r) =-\frac{a}{1+d\rho^{-1/3}}\left\{\rho+b\rho^{-2/3}\left[C_{\mathrm{F}}\rho^{5/3}-2t_{\mathrm{W}}+\frac{1}{9}\left(t_{\mathrm{W}}+\frac{1}{2}\boldsymbol{\nabla}^{2}\rho\right)\right]\mathrm{e}^{-c\rho^{-1/3}}\right\} \tag{3.208} $$
式（3.208）为闭壳层情形的 LYP 关联能体密度，$E_{\mathrm c}^{\mathrm{LYP}}=\int e_{\mathrm c}^{\mathrm{LYP}}(\boldsymbol r)\,\mathrm d^3\boldsymbol r$；它不是每电子关联能，也不能直接用于任意自旋极化密度。这里 $t_{\mathrm W}=\frac18\bigl(|\nabla\rho|^2/\rho-\nabla^2\rho\bigr)$，$C_{\mathrm{F}}=\frac3{10}(3\pi^{2})^{2/3}$，$a=0.04918$，$b=0.132$，$c=0.2533$，$d=0.349$。
2. PW91\cite{perdew1992atoms} 泛函
在 PW91 泛函中，将 $\varepsilon_{\mathrm{x}}$ 表示为
$$ \begin{aligned} \varepsilon_{\mathrm{x}}^{\mathrm{PW91}}(r_{\mathrm s}) &=\varepsilon_{\mathrm{x}}^{\mathrm{hom}}(r_{\mathrm s})F_{\mathrm{x}}^{\mathrm{PW91}}(s_1),\\ F_{\mathrm{x}}^{\mathrm{PW91}}(s_1) &=\frac{1+0.19645s_1\sinh^{-1}(7.7956s_1) +(0.2743-0.1508e^{-100s_1^2})s_1^2} {1+0.19645s_1\sinh^{-1}(7.7956s_1)+0.004s_1^4}. \end{aligned} \tag{3.209} $$
可以看出，在 $s_1$ 较小时，有
$$ F_{\mathrm{x}}\simeq1+0.1234s_1^{2}+O(|\,\boldsymbol{\nabla}\rho\,|^{4}) $$
即方程（3.206）。
而关联能表示为
$$ E_{\mathrm{c}}^{\mathrm{PW91}}[\rho^{\uparrow},\rho^{\downarrow}]=\int\rho(\boldsymbol{r})[\varepsilon_{\mathrm{c}}^{\mathrm{hom}}(r_{\mathrm{s}},\zeta)+H(r_{\mathrm{s}},t,\zeta)]\,\mathrm{d}\boldsymbol{r} \tag{3.210} $$
式中：$t=|\,\boldsymbol{\nabla}\rho\,|/(2\phi(\zeta)\boldsymbol{k}_{\mathrm{s}}\rho)$，其中 $\boldsymbol{k}_{\mathrm{s}}$ 是 Thomas-Fermi 屏蔽波矢，$|\,\boldsymbol{k}_{\mathrm{s}}\,|=\sqrt{4k_{\mathrm{F}}/\pi}$，而 $\phi(\zeta)=[(1+\zeta)^{2/3}+(1-\zeta)^{2/3}]/2$。
方程（3.210）中的第二项 $H=H_0+H_1$，其中
$$ H_0=\phi^3(\zeta)\frac{\beta^2}{2\alpha} \ln\!\left[1+\frac{2\alpha}{\beta} \frac{t^2+At^4}{1+At^2+A^2t^4}\right] \tag{3.211} $$
式中
$$ A=\frac{2\alpha/\beta} {\exp\!\left[-\frac{2\alpha\varepsilon_{\mathrm c}^{\mathrm{hom}}(r_{\mathrm s},\zeta)}{\phi^3(\zeta)\beta^2}\right]-1} \tag{3.212} $$
其中：$\alpha=0.09$，$\beta=\nu C_{\mathrm{c}}(0)=(16/\pi)(3\pi^{2})^{1/3}\times0.004235$。而
$$ H_1=\nu[C_{\mathrm{c}}(r_{\mathrm{s}})-C_{\mathrm{c}}(0)-3C_{\mathrm{x}}/7]\phi^{3}(\zeta)t^{2}\times\exp[-100\phi^{4}(\zeta)(k_{\mathrm{s}}^{2}/k_{\mathrm{F}}^{2})t^{2}] \tag{3.213} $$
式中：$C_{\mathrm{x}}=-0.001667$；$C_{\mathrm{c}}(r_{\mathrm{s}})$ 由文献\cite{rasolt1986exchange}给出，且有
$$ C_{\mathrm{c}}(r_{\mathrm{s}})=10^{-3}\frac{2.568+ar_{\mathrm{s}}+br_{\mathrm{s}}^{2}}{1+cr_{\mathrm{s}}+dr_{\mathrm{s}}^{2}+10br_{\mathrm{s}}^{3}} \tag{3.214} $$
其中：$a=23.266$，$b=7.389\times10^{-3}$，$c=8.723$，$d=0.472$。
3. PBE\cite{perdew1996generalized} 泛函
PBE 泛函的形式与 PW91 泛函类似，但是其形式更简单，而且在实际应用中取得了比较好的效果，因此是目前材料计算里被广泛采用的一种 GGA 泛函。PBE 泛函理论规定
$$ \varepsilon_{\mathrm{x}}^{\mathrm{PBE}}=\varepsilon_{\mathrm{x}}^{\mathrm{hom}}(r_{\mathrm{s}})F_{\mathrm{x}}(r_{\mathrm{s}},\zeta,s_1)=\varepsilon_{\mathrm{x}}^{\mathrm{hom}}\left(1+\kappa-\frac{\kappa}{1+\mu s_1^{2}/\kappa}\right) \tag{3.215} $$
其中：$\mu\simeq0.21915$，$\kappa=0.804$。$F_{\mathrm{x}}$ 的这个形式保证了方程（3.215）在 $s_1\to0$ 的情况下回到 LDA，而且满足式（3.204）。
与 PW91 泛函类似，PBE 泛函也将 $E_{\mathrm{c}}$ 表示为两项之和，即
$$ E_{\mathrm{c}}^{\mathrm{PBE}}[\rho^{\uparrow},\rho^{\downarrow}]=\int\rho(\boldsymbol{r})[\varepsilon_{\mathrm{c}}^{\mathrm{hom}}(r_{\mathrm{s}},\zeta)+H(r_{\mathrm{s}},t,\zeta)]\,\mathrm{d}\boldsymbol{r} \tag{3.216} $$
式中
$$ H(r_{\mathrm{s}},t,\zeta)=\frac{e^{2}}{a_0}\gamma\phi^{3}(\zeta)\times\ln\left[1+\frac{\beta}{\gamma}t^{2}\left(\frac{1+At^{2}}{1+At^{2}+A^{2}t^{4}}\right)\right] \tag{3.217} $$
而
$$ A=\frac{\beta}{\gamma}\frac{1}{\exp[-\varepsilon_{\mathrm{c}}^{\mathrm{hom}}(r_{\mathrm{s}},\zeta)/(\gamma\phi^{3}(\zeta)e^{2}/a_0)]-1} \tag{3.218} $$
方程（3.217）及方程（3.218）中，$\beta\simeq0.066725$，$\gamma=(1-\ln2)/\pi^{2}\simeq0.031091$。$t$、$\phi(\zeta)$、$k_{\mathrm{s}}$ 均与 PW91 泛函中的定义相同。注意 PBE 泛函中，$e$ 和 $a_0$ 均取原子单位。
通常情况下，由于考虑了对电子密度梯度的修正，GGA 泛函的计算结果要比 LDA 泛函的精确，这一点的主要体现是原子能量、晶体结合能、体系键长、键角等的值在 GGA 中可以更接近实验结果。但是也存在例外情况。例如，在描述金属和氧化物的表面能时，GGA 方法反而不如 LDA 方法合理。
PW91 泛函和 PBE 泛函被广泛运用在各种晶体性质计算中。在计算各种分子在贵金属表面的吸附能时，这两种形式的 GGA 泛函都有高估吸附能的趋势。为了解决这个问题，Hammer 提出了 Revised-Perdew-Burke-Ernzerhof（RPBE）泛函\cite{hammer1999improved}。但是最近的研究表明，为了精确地计算分子在贵金属表面的吸附能，需要引入自能修正或者范德瓦尔斯修正，这里不做详细讨论。
GGA 下的交换关联势
因为 GGA 包含了密度梯度，所以其交换关联势的计算也相对复杂。设当电子密度改变 $\delta\rho$、密度梯度改变 $\delta\boldsymbol{\nabla}\rho=\boldsymbol{\nabla}\delta\rho$ 时，交换关联能改变 $\delta E_{\mathrm{xc}}[\rho,\boldsymbol{\nabla}\rho]$，则
$$ \delta E_{\mathrm{xc}}=\int\left[ \left(\varepsilon_{\mathrm{xc}}+\rho\frac{\partial\varepsilon_{\mathrm{xc}}}{\partial\rho}\right)\delta\rho +\rho\frac{\partial\varepsilon_{\mathrm{xc}}}{\partial(\boldsymbol\nabla\rho)}\cdot\boldsymbol\nabla(\delta\rho) \right]\mathrm d\boldsymbol r \tag{3.219} $$
式（3.219）右端方括号中的前两项即 LDA 中的交换关联势；对最后一项做分部积分，可得\cite{martin2004electronic}
$$ V_{\mathrm{xc}}=\varepsilon_{\mathrm{xc}}+\rho\frac{\partial\varepsilon_{\mathrm{xc}}}{\partial\rho} -\boldsymbol\nabla\!\cdot\!\left(\rho\frac{\partial\varepsilon_{\mathrm{xc}}}{\partial(\boldsymbol\nabla\rho)}\right) \tag{3.220} $$
通常交换能部分与关联能部分遵循各自独立的方程，因此（3.220）也相应地分为 $V_{\mathrm{x}}$ 和 $V_{\mathrm{c}}$ 两部分。对于考虑自旋自变量 $\sigma$ 的情况，直接给出 $V_{\mathrm{xc}}^{\sigma}$ 的表达式：
$$ V_{\mathrm{xc}}^\sigma=\frac{\partial f_{\mathrm{xc}}}{\partial\rho^\sigma} -\boldsymbol\nabla\!\cdot\!\frac{\partial f_{\mathrm{xc}}}{\partial(\boldsymbol\nabla\rho^\sigma)},\qquad f_{\mathrm{xc}}=\rho\varepsilon_{\mathrm{xc}}(\rho^\uparrow,\rho^\downarrow,\boldsymbol\nabla\rho^\uparrow,\boldsymbol\nabla\rho^\downarrow) \tag{3.221} $$
混合泛函
在电子结构的计算中，自关联项和交换项没有办法抵消，因此引起了较大的误差。为了解决这个问题，在混合泛函方法中，交换关联势的表达式加入了部分 Hartree-Fock 的精确交换。这类泛函包括 B3LYP、HSE03\cite{heyd2003hybrid}、HSE06\cite{heyd2006erratum} 等。
此混合泛函在基于局域基组的量子化学计算中得到了广泛的应用，但是在基于平面波的程序中，由于非局域的交换项计算比较困难，因此应用受到一定的限制。最新发展的 HSE06 等平面波基的混合泛函，通过将交换项分解成长程部分和短程部分，并仅对短程部分加入精确交换势，从而使计算量减小，为混合泛函在平面波基的密度泛函程序中的应用铺平了道路。
混合泛函中，体系的交换关联能往往表示为
$$ E_{\mathrm{xc}}^{\mathrm{HF}}=\alpha E_{\mathrm{x}}^{\mathrm{HF}}+(1-\alpha)E_{\mathrm{x}}^{\mathrm{DFT}}+E_{\mathrm{c}}^{\mathrm{DFT}} \tag{3.222} $$
式中：$\alpha$ 是一个可调参数。在 HSE03 和 HSE06 中，一般取 $\alpha=0.25$。
强关联与 LDA＋$U$ 方法
含有部分占据且较局域的 $d$ 或 $f$ 电子的体系，如某些过渡金属氧化物和稀土化合物，可能受到较强的局域电子相互作用影响。半局域 LDA/GGA 对其中一些绝缘体会低估带隙，甚至错误预测为金属；也有真实的强关联金属，不能把此类体系一概描述为半导体或绝缘体。局域轨道与巡游态可同时存在，具体电子结构取决于材料和所处相。
对局域相互作用重要的材料，可在适当模型和参数条件下使用 DFT＋$U$ 等方法改善描述；GW 用于准粒子激发，但单次 GW 本身并不能普遍解决所有强关联问题。DFT＋$U$ 的结果依赖相关子空间、$U/J$ 参数及双计数形式，其计算代价通常低于更完整的多体方法。我们将在本节中对 LDA＋$U$ 方法进行简要的介绍。
LDA＋$U$ 的理论推导需要用到二次量子化的知识，这超出了本书的范围。因此，这里只给出重要的结果并做相关讨论。以 $d$ 轨道为例，在不考虑自旋极化的情况下，可以将体系中的电子分为两个亚系统，分别是局域性较强的 $d$ 电子与离域性较强的 $s$ 电子和 $p$ 电子。后者的相互作用可以用 LDA 描述，而 $d$ 电子与 $d$ 电子间的库仑相互作用（简称 $d$-$d$ 库仑作用）则写为
$$ E_{d\text{-}d}=\frac{U}{2}\sum_{i\neq j}n_in_j \tag{3.223} $$
式中：$n_i$ 和 $n_j$ 分别为第 $i$ 个和第 $j$ 个 $d$ 轨道上的电子占据数；$U$ 为库仑参数，取正值。此时体系的总能可以表示为
$$ E_{\mathrm{tot}}=E_{\mathrm{LDA}}+E_{d\text{-}d}-E_{\mathrm{d.\,c.}} \tag{3.224} $$
式中：右端第一项就是普通的 LDA 近似下体系基态总能，参见方程（3.181）；第二项由式（3.223）给出；第三项是冗余项，这是因为 $E_{\mathrm{LDA}}$ 中已经包含了 $d$-$d$ 库仑作用，因此需要将这部分重复计入的能量作为冗余项排除。Anisimov 等人假设 LDA 中，$d$-$d$ 库仑作用只与总的 d 轨道占据数 $N$ 相关，因此可以将 $E_{\mathrm{d.\,c.}}$ 写为\cite{anisimov1997first}
$$ E_{\mathrm{d.\,c.}}=UN(N-1)/2 \tag{3.225} $$
式中 $N=\displaystyle\sum_{i}n_i$。由此可以得到 LDA＋$U$ 方法下第 $i$ 个 $d$ 轨道的本征值 $\varepsilon_i$ 为
$$ \varepsilon_i=\frac{\partial E_{\mathrm{tot}}}{\partial n_i}=\varepsilon_{i,\mathrm{LDA}}+U\left(\frac{1}{2}-n_i\right) \tag{3.226} $$
式（3.226）表明，与 LDA 的结果相比，被占据的 $d$ 轨道能量下移 $U/2$，未被占据的 $d$ 轨道能量上移 $U/2$。因此引入 $U$ 有助于改进被低估的能隙。但是因为能量表达式改变，所以相应的 Kohn-Sham 方程中的 $V_{\mathrm{eff}}$ 和哈密顿矩阵都要做相应的修改。一般而言，将 $d$ 轨道或者 $f$ 轨道用一组正交的局域轨道基 $|\,i,nlm,\sigma\rangle$ 展开，其中 $i$ 表示格点，$nlm$ 为轨道基的量子数，而 $\sigma$ 代表自旋。为了简化推导，在这里认为只有特定的 $nl$ 轨道需要利用 $U$ 来准确描述，因此，只有磁量子数 $m$ 可以变化。由此可定义格点 $i$ 上的密度矩阵元 $n_{i,mm'}^{\sigma}$：
$$ n_{i,mm'}^{\sigma}=\sum_{n\boldsymbol{k}}f_{n\boldsymbol{k}\sigma} \langle i,nlm,\sigma\,|\,\psi_{n\boldsymbol{k}\sigma}\rangle \langle\psi_{n\boldsymbol{k}\sigma}\,|\,i,nlm',\sigma\rangle \tag{3.227} $$
这里 $|i,nlm,\sigma\rangle$ 为待修正的局域轨道，$|\psi_{n\boldsymbol{k}\sigma}\rangle$ 为单粒子态，$f_{n\boldsymbol{k}\sigma}$ 为其占据数；若布里渊区权重未并入 $f$，则求和中还须乘相应权重。密度矩阵的两个轨道指标可以不同。写出普遍的 LSDA＋$U$ 的总能\cite{anisimov1997first,liechtenstein1995density}：
$$ E^{\mathrm{LSDA+U}}[\rho^{\sigma}(\boldsymbol{r})\{n_i^{\sigma}\}]=E^{\mathrm{LSDA}}+E^{U}[\{n_i^{\sigma}\}]-E_{\mathrm{d.\,c.}}[\{n_i^{\sigma}\}] \tag{3.228} $$
式中：$\rho^{\sigma}(\boldsymbol{r})$ 是自旋态为 $\sigma$ 的电子密度。右端第一项 $E^{\mathrm{LSDA}}$ 由式（3.181）、式（3.195）及式（3.199）给出。第二项为
$$ E^{U}[\{n_i^{\sigma}\}]=\frac12\sum_i\sum_{abcd}\sum_{\sigma\sigma'} n_{i,ba}^{\sigma}n_{i,dc}^{\sigma'} \left[\langle ac|V_{\mathrm{ee}}|bd\rangle -\delta_{\sigma\sigma'}\langle ac|V_{\mathrm{ee}}|db\rangle\right] \tag{3.229} $$
式中 $a,b,c,d$ 分别遍历同一相关壳层的局域轨道；$\langle ac|V_{\mathrm{ee}}|bd\rangle$ 是两电子库仑积分，交换项仅在 $\sigma=\sigma'$ 时出现。 $V_{\mathrm{ee}}$ 为处于 $\boldsymbol{r}(r,\theta,\phi)$ 和 $\boldsymbol{r}'(r',\theta',\phi')$ 终点的两个点电荷之间的库仑相互作用，用球谐函数展开为
$$ V_{\mathrm{ee}}(\boldsymbol{r},\boldsymbol{r}')=\frac{1}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}=\sum_{l}\frac{4\pi}{2l+1}\frac{r_{\lt }^{l}}{r_{\gt }^{l+1}}\sum_{m=-l}^{l}\mathrm{Y}_l^{m}(\theta,\phi)\mathrm{Y}_l^{m*}(\theta',\phi') \tag{3.230} $$
式中 $r_\lt =\min(r,r')$、$r_\gt =\max(r,r')$。
式（3.229）中的第一项积分可以写为
$$ \begin{aligned} &\langle m,m''\,|\,V_{\mathrm{ee}}\,|\,m',m'''\rangle\\ & =\iint\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}' R_l^{*}(r)\mathrm{Y}_l^{m*}(\theta,\phi) R_l(r)\mathrm{Y}_l^{m'}(\theta,\phi)V_{\mathrm{ee}} R_l^{*}(r')\mathrm{Y}_l^{m''*}(\theta',\phi') R_l(r')\mathrm{Y}_l^{m'''}(\theta',\phi')\\ &=\sum_{k=0}^{2l}a_k(m,m',m'',m''')F^{k} \end{aligned} \tag{3.231} $$
式（3.231）将径向积分提出为同一组 $F^k$，以壳层内各 $m$ 共用径向函数 $R_l(r)$ 为前提。 $F^{k}$ 包括了径向函数积分，称为屏蔽 Slater 积分\cite{judd1963operator}；$a_k$ 称为 Gaunt 系数，有
$$ a_k(m,m',m'',m''')=\sum_{q=-k}^{k}\frac{4\pi}{2k+1}\langle lm\,|\,\mathrm{Y}_k^{q}\,|\,lm'\rangle\langle lm''\,|\,\mathrm{Y}_k^{q*}\,|\,lm'''\rangle \tag{3.232} $$
式（3.229）中其余两项积分也可以类似地写成上述形式。描述 $d$ 电子，需要 $F^{0}$、$F^{2}$ 及 $F^{4}$，描述 $f$ 电子，则还需要 $F^{6}$。稍后讨论哈密顿矩阵元时，我们还要再讨论 $F^{k}$。
LSDA＋$U$ 总能表达式中的冗余项 $E_{\mathrm{d.\,c.}}$ 为
$$ E_{\mathrm{d.\,c.}}[\{n_i^{\sigma}\}] =\sum_i\left\{\frac{U}{2}N_i(N_i-1) -\frac{J}{2}\sum_{\sigma}N_i^{\sigma}(N_i^{\sigma}-1)\right\} \tag{3.233} $$
式中
$$ N_i^{\sigma}=\operatorname{tr}n_i^{\sigma}=\sum_m n_{i,mm}^{\sigma}, \qquad N_i=N_i^{\uparrow}+N_i^{\downarrow} $$
其中 $\sigma$ 表示 ↑ 或 ↓。
式（3.233）中出现了两个参数 $U$ 和 $J$，它们分别表示局域轨道的平均库仑作用与平均 Hund 交换作用。可以看到，$J$ 的存在部分抵消了 $U$ 所描述的排斥作用，这一点我们在讨论 Hartree-Fock 方程的时候已经发现了。
将式（3.229）至式（3.233）代入 LSDA＋$U$ 的总能表达式，并对 $n_{i,mm'}^{\sigma}$ 求变分，可以得到作用在格点 $i$（或称第 $i$ 个 $d$ 轨道或 $f$ 轨道）的 Kohn-Sham 有效势 $V_{i,\mathrm{eff}}^{\sigma}$：
$$ V_{i,\mathrm{eff}}^{\sigma}=V_{\mathrm{KS}}^{\mathrm{LSDA}}+\sum_{mm'}|\,i,nlm,\sigma\rangle V_{i,mm'}^{\sigma}\langle i,nlm',\sigma\,| \tag{3.234} $$
$V_{\mathrm{KS}}^{\mathrm{LSDA}}$ 即为通常 LSDA 近似下的 Kohn-Sham 有效势（见方程（3.180）、方程（3.199）），而附加的一项则代表了 $d$ 电子或 $f$ 电子相互作用的影响，其中 $V_{i,mm'}^{\sigma}$ 可写为
$$ \begin{aligned} V_{i,mm'}^{\sigma} ={}&\frac{\partial(E^U-E_{\mathrm{d.\,c.}})} {\partial n_{i,m'm}^{\sigma}}\\ ={}&\sum_{cd,\sigma'}n_{i,dc}^{\sigma'} \left[\langle mc|V_{\mathrm{ee}}|m'd\rangle -\delta_{\sigma\sigma'}\langle mc|V_{\mathrm{ee}}|dm'\rangle\right]\\ &-\delta_{mm'}\left[U\left(N_i-\frac12\right) -J\left(N_i^{\sigma}-\frac12\right)\right]. \end{aligned} \tag{3.235} $$
至此，我们构建了包含自旋极化的 LDA＋$U$ 理论框架，但是屏蔽 Slater 积分并未给出，而且在实际工作中需要设定的是 $U$ 和 $J$，所以必须给出 $F^{k}$ 和 $U$、$J$ 之间的关系。对于 $d$ 电子，有\cite{degroot1990xray}
$$ U=F^{0},\quad J=\frac{F^{2}+F^{4}}{14},\quad\frac{F^{4}}{F^{2}}=0.625 $$
对于 $f$ 电子，计算程序 VASP 采用
$$ U=F^{0},\quad J=\frac{286F^{2}+195F^{4}+250F^{6}}{6435},\quad\frac{F^{4}}{F^{2}}=0.668,\quad\frac{F^{6}}{F^{2}}=0.494 $$