本文是「第一性原理的微观计算模拟」系列的第 5 篇（共 11 篇），内容整理自同名书稿第 3 章，公式编号与原书一致。文中引用的文献依据原书参考文献清单整理并列于文末，按原书清单顺序从 1 开始编号。\nocite{*}
在 3.2.3 节中我们详细介绍了 Kohn-Sham 方程，但是距离 Kohn-Sham 方程的具体求解尚有一定距离。一般可以选择三类基函数来展开波函数。第一类是平面波函数，其在空间中没有固定的参考点。第二类是局域波函数，例如原子轨道或者高斯基函数等。第三类则是混合基组，即将平面波函数“缀加”于局域波函数上作为基函数。选取不同的基函数，则相应的 Kohn-Sham 方程形式、哈密顿矩阵元表达式及总能的表达式会有显著差异。在本节中我们以平面波-赝势框架下 Kohn-Sham 方程具体求解过程为例，对第一性原理计算程序的若干关键点进行详细讨论。
布里渊区积分——特殊 $\boldsymbol{k}$ 点
在各种周期性边界条件下的第一性原理计算方法中，往往会涉及在布里渊区积分的问题，例如总能、电荷密度分布，以及金属体系中费米面的确定等等。为了提高计算效率，需要寻找一种高效的积分方法，可以通过较少的 $\boldsymbol{k}$ 点运算取得较高的精度。这些 $\boldsymbol{k}$ 点称为平均值点或者特殊点，而这种方法就称为特殊 $\boldsymbol{k}$ 点法。
特殊 $\boldsymbol{k}$ 点法基本思想
Chadi 和 Cohen 最早提出了这种特殊 $\boldsymbol{k}$ 点法的数学基础\cite{chadi1973special}。考虑一个光滑周期性函数 $g(\boldsymbol{k})$，周期为 $\boldsymbol{G}$，可以将其展开为如下傅里叶级数：
$$ g(\boldsymbol{k})=g_0+\sum_{m=1}g_m\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{R}_m} \tag{3.296} $$
式中：$\boldsymbol{R}_m$ 是与倒格矢 $\boldsymbol{G}$ 相应的晶体格子，其对称性用对称点群 $G$ 来描述。假设另有一个拥有体系全部对称性的函数 $f(\boldsymbol{k})$ 满足条件
$$ f(\mathrm{T}\boldsymbol{k})=f(\boldsymbol{k}),\quad\forall\,\mathrm{T}\in G $$
则可以将 $f(\boldsymbol{k})$ 用 $g(\boldsymbol{k})$ 展开为
$$ f(\boldsymbol{k})=\frac{1}{n_G}\sum_{i}g(T_i\boldsymbol{k})=g_0+\sum_{m=1}^{\infty}\sum_{i}\frac{1}{n_G}g_m\mathrm{e}^{\mathrm{i}T_i\boldsymbol{k}\cdot\boldsymbol{R}_m} \tag{3.297} $$
式中：$n_G$ 是点群 $G$ 的阶数。设 $g_0=f_0$，将式（3.297）的求和顺序重新调整可以得到
$$ f(\boldsymbol{k})=f_0+\sum_{m=1}\frac{g_m}{n_G}\sum_{T_i\in G}\mathrm{e}^{\mathrm{i}\boldsymbol{k}T_i^{-1}\cdot\boldsymbol{R}_m}=f_0+\sum_{m=1}f_m\sum_{|\boldsymbol{R}|=C_m}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{R}}=f_0+\sum_{m=1}f_m\mathrm{A}_m(\boldsymbol{k}) \tag{3.298} $$
式中：$C_m$ 是距离原点第 $m$ 近邻的球半径，按升序排列，$C_m\leqslant C_{m+1}$。注意限制条件 $C_m\leqslant|\,\boldsymbol{R}\,|\leqslant C_{m+1}$ 具有球对称性，也即高于 $G$ 的对称性，所以满足限制条件的格点集合 $\{\boldsymbol{R}\}$ 并不一定可以通过 $G$ 中的操作联系起来。方程（3.298）中的函数 $A_m$ 满足下列条件：
$$ \begin{cases} \dfrac{\varOmega}{(2\pi)^{3}}\displaystyle\int_{\mathrm{BZ}}A_m(\boldsymbol{k})\,\mathrm{d}\boldsymbol{k}=0,\quad\forall\,(m\gt 0,m\in\mathbf{Z})\\ \dfrac{\varOmega}{(2\pi)^{3}}\displaystyle\int_{\mathrm{BZ}}A_m(\boldsymbol{k})A_n(\boldsymbol{k})\,\mathrm{d}\boldsymbol{k}=\mathrm{N}_n\delta_{nm}\\ A_m(\boldsymbol{k}+\boldsymbol{G})=A_m(\boldsymbol{k})\\ A_m(T_i\boldsymbol{k})=A_m(\boldsymbol{k})\\ A_m(\boldsymbol{k})A_n(\boldsymbol{k})=\displaystyle\sum_{j}a(j,m,n)A_j(\boldsymbol{k}) \end{cases} \tag{3.299} $$
式中：$\boldsymbol{G}$ 是倒格矢；$N_n$ 是满足条件 $|\,\boldsymbol{R}\,|=C_n$ 的格点数。
式（3.299）中后四个方程分别表明函数 $A_m(\boldsymbol{k})$ 在第一布里渊区内的正交性、周期性、体系对称性和完备性，第一个方程则给出了对 $A_m(\boldsymbol{k})$ 的要求。对特殊 $\boldsymbol{k}$ 点法而言，前两个方程更为重要。
注意到式（3.299）中的求和从 $m=1$ 开始，因此需要对 $m=0$ 的情况进行单独定义。定义 $A_0(\boldsymbol{k})=1$，则函数 $f(\boldsymbol{k})$ 的平均值为
$$ \bar{f}=\frac{\varOmega}{(2\pi)^{3}}\int_{\mathrm{BZ}}f(\boldsymbol{k})\,\mathrm{d}\boldsymbol{k}=f_0 \tag{3.300} $$
由方程（3.298）可知，如果存在 $\boldsymbol{k}_0$，满足
$$ A_m(\boldsymbol{k}_0)=0,\quad\forall\,(m\gt 0,m\in\mathbf{Z}) \tag{3.301} $$
那么立刻可以得到 $f=f_0=f(\boldsymbol{k}_0)$，这样的 $\boldsymbol{k}_0$ 点即为平均值点。但是满足上述条件的 $\boldsymbol{k}$ 点并不是普遍存在的，所以需要构建满足一定条件的集合 $\{\boldsymbol{k}\}$，利用这些点上函数值的加权平均计算 $f_0$。也即
$$ \begin{cases} \displaystyle\sum_{i=1}^{n}\alpha_iA_m(\boldsymbol{k}_i)=0,\quad m=1,2,\cdots,N\\ \displaystyle\sum_{i}\alpha_i=1 \end{cases} \tag{3.302} $$
式中 $N$ 可以取有限值。
利用方程（3.298），可以得到
$$ \sum_{i=1}^{n}\alpha_i f(\boldsymbol k_i) =f_0+\sum_{m=N+1}^{\infty}f_m\sum_{i=1}^{n}\alpha_i A_m(\boldsymbol k_i) \tag{3.303} $$
根据方程（3.303），有
$$ f_0=\sum_{i=1}^{n}\alpha_i f(\boldsymbol k_i) -\sum_{m=N+1}^{\infty}f_m\sum_{i=1}^{n}\alpha_i A_m(\boldsymbol k_i) \tag{3.304} $$
考虑到 $f_m$ 随 $m$ 的增大而迅速减小的性质，可以近似地得到 $f(\boldsymbol{k})$ 的平均值，即
$$ \bar f=f_0\approx\sum_{i=1}^{n}\alpha_i f(\boldsymbol k_i) \tag{3.305} $$
而将方程（3.304）右端的第二项作为可控误差。因此，如果可以找到一组 $\boldsymbol{k}$ 点，使得集合中的 $\boldsymbol{k}$ 点尽量少，而且这些 $\boldsymbol{k}$ 点在 $N$ 尽量大的情况下满足方程（3.303），则我们进行布里渊区积分的时候就可以尽可能快地得到精度较高的结果。这正是特殊 $\boldsymbol{k}$ 点法的要点所在。反过来讲，这也表明进行具体计算的时候我们需要对计算精度进行测试，也即保证所取 $\boldsymbol{k}$ 点使得式（3.304）右端第二项足够小。
Chadi-Cohen 方法
在 3.4.1.1 节我们讨论了 $\boldsymbol{k}$ 点的可行性。Chadi 和 Cohen 提出了一套可以得出这些特殊 $\boldsymbol{k}$ 点的方法\cite{chadi1973special}。首先找出两个特殊 $\boldsymbol{k}$ 点——$\boldsymbol{k}_1$、$\boldsymbol{k}_2$，二者分别在 $\{N_1\}$ 和 $\{N_2\}$ 的情况下满足
$$ A_m(\boldsymbol{k})=0 $$
然后通过这两个 $\boldsymbol{k}$ 点构造新的 $\boldsymbol{k}$ 点集合：
$$ \boldsymbol{k}_i=\boldsymbol{k}_1+T_i\boldsymbol{k}_2 $$
且权重 $\alpha_i=\dfrac{1}{n_G}$。下面证明 $\boldsymbol{k}_i$ 在 $\{N_1\}\cup\{N_2\}$ 的情况下仍然满足方程（3.303）。
根据 $\boldsymbol{k}_1$ 和 $\boldsymbol{k}_2$ 的定义可知，对于 $m\in\{N_1\}$ 和 $m\in\{N_2\}$，有
$$ A_m(\boldsymbol{k}_1)A_m(\boldsymbol{k}_2)=0 $$
即
$$ \left(\sum_{|\boldsymbol{R}|=C_m}\mathrm{e}^{\mathrm{i}\boldsymbol{k}_1\cdot\boldsymbol{R}}\right)\left(\sum_{|\boldsymbol{R}|=C_m}\mathrm{e}^{\mathrm{i}\boldsymbol{k}_2\cdot\boldsymbol{R}}\right)=0 \tag{3.306} $$
由式（3.306）可进行如下推导：
$$ \begin{aligned} &\left(\sum_{|\boldsymbol{R}|=C_m}\mathrm{e}^{\mathrm{i}\boldsymbol{k}_1\cdot\boldsymbol{R}}\right)\left(\sum_{i}\mathrm{e}^{\mathrm{i}\boldsymbol{k}_2\cdot T_i\boldsymbol{R}}\right)=\left(\sum_{|\boldsymbol{R}|=C_m}\mathrm{e}^{\mathrm{i}\boldsymbol{k}_1\cdot\boldsymbol{R}}\right)\left(\sum_{l}\mathrm{e}^{\mathrm{i}T_l\boldsymbol{k}_2\cdot\boldsymbol{R}}\right)=\sum_{l}\sum_{|\boldsymbol{R}|=C_m}\mathrm{e}^{\mathrm{i}(\boldsymbol{k}_1+T_l\boldsymbol{k}_2)\cdot\boldsymbol{R}}=0\\ &\Rightarrow\sum_{l}A_m(\boldsymbol{k}_1+T_l\boldsymbol{k}_2)=0\\ &\Rightarrow\sum_{l}A_m(\boldsymbol{k}_l)=0 \end{aligned} $$
因此可以用这种方法产生一系列 $\boldsymbol{k}$ 点，用于计算布里渊区内的积分。如果此时的精度不够，则利用同样的方法继续生成新的 $\boldsymbol{k}$ 点集合，从而改进精度，即
$$ \boldsymbol{k}_{ii}=\boldsymbol{k}_i+T_i\boldsymbol{k}_3 $$
式中：$\boldsymbol{k}_3$ 为在 $m\in\{N_3\}$ 情况下满足 $A_m(\boldsymbol{k})=0$ 的特殊 $\boldsymbol{k}$ 点。
利用点群对称性，可将整个网格约化至不可约布里渊区。设点群阶数为 $n_G$，代表点 $\boldsymbol{k}_i$ 的稳定子（波矢群）阶数为 $n_i$，则其星形轨道含 $\nu_i=n_G/n_i$ 个等价点。若原网格各点等权，约化后的归一化权重等于星形轨道的相对大小：
$$ \omega_{\boldsymbol{k}_i}=\frac{\nu_i}{\displaystyle\sum_j\nu_j}, \qquad \nu_i=\frac{n_G}{n_i} \tag{3.307} $$
若从其他构造得到未归一化权重 $\alpha_i$，则可写为
$$ \omega_{\boldsymbol{k}_i}=\frac{\alpha_i}{\displaystyle\sum_{j}\alpha_j} \tag{3.308} $$
Monkhorst-Pack 方法
上述 Chadi-Cohen 方法非常巧妙，但是在具体应用中，必须首先确定 2～3 个性能比较好的 $\boldsymbol{k}$ 点，由此构建出的 $\boldsymbol{k}$ 点集合才拥有比较高的效率和精度。因此，对于每一个具体问题，在计算之前都必须经过相当多的对称性分析。对程序编写而言，这是一个相当烦琐的任务。Monkhorst 和 Pack 提出了一种简单的产生 $\boldsymbol{k}$ 点网格的方法，同时又可使方程（3.303）得到满足，这就是通常所说的 Monkhorst-Pack 方法\cite{monkhorst1976special}。
晶体中的格点 $\boldsymbol{R}$ 总可以表示为 $\boldsymbol{R}=R_1\boldsymbol{a}_1+R_2\boldsymbol{a}_2+R_3\boldsymbol{a}_3$，其中 $\boldsymbol{a}_i$ 是实空间三个方向上的基矢。Monkhorst 和 Pack 建议按如下方法划分布里渊区：
$$ u_r=(2r-q-1)/(2q),\quad1\leqslant r\leqslant q \tag{3.309} $$
将 $\boldsymbol{k}$ 点写为分量形式，则可得到如下表达式：
$$ \boldsymbol{k}_{prs}=u_p\boldsymbol{b}_1+u_r\boldsymbol{b}_2+u_s\boldsymbol{b}_3 \tag{3.310} $$
式中：$\boldsymbol{b}_1$、$\boldsymbol{b}_2$、$\boldsymbol{b}_3$ 是倒空间的基矢。与 Chadi-Cohen 方法相似，Monkhorst-Pack 方法定义函数 $A_m$ 为
$$ \begin{cases} A_m(\boldsymbol{k})=\dfrac{1}{\sqrt{\mathrm{N}_m}}\displaystyle\sum_{|\boldsymbol{R}|=C_m}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{R}},&m\gt 1\\ A_1(\boldsymbol{k})=1 \end{cases} \tag{3.311} $$
则相应于式（3.299）中的 $\dfrac{\varOmega}{(2\pi)^{3}}\displaystyle\int_{\mathrm{BZ}}A_m(\boldsymbol{k})A_n(\boldsymbol{k})\,\mathrm{d}\boldsymbol{k}$，可以计算方程（3.309）所生成的离散化网格点上相同的量：
$$ S_{mn}(q)=\frac{1}{q^{3}}\sum_{p,r,s=1}^{q}A_m^*(\boldsymbol{k}_{prs})A_n(\boldsymbol{k}_{prs}) =\frac{1}{\sqrt{N_mN_n}}\sum_{a=1}^{N_m}\sum_{b=1}^{N_n}\prod_{j=1}^{3}W_j^{ab}(q),\quad N_1=1 \tag{3.312} $$
式中
$$ W_j^{ab}(q)=\frac{1}{q}\sum_{r=1}^{q} \exp\!\left[\frac{\mathrm{i}\pi(2r-q-1)}{q}(R_j^{b}-R_j^{a})\right] \tag{3.313} $$
注意到 $\boldsymbol{R}_j^{a}$ 和 $\boldsymbol{R}_j^{b}$（$j=1,2,3$）都是整数，因此可以算出
$$ W_j^{ab}(q)=\begin{cases} (-1)^{(q+1)d_j/q},&d_j\in q\mathbb Z,\\ 0,&d_j\notin q\mathbb Z, \end{cases}\qquad d_j=R_j^b-R_j^a \tag{3.314} $$
第三种情况来自等比数列求和为零，并非 $W_j^{ab}(q)$ 是奇函数。引入限制条件
$$ \begin{cases}|\,\boldsymbol{R}_j^{a}\,|\lt q/2\\|\,\boldsymbol{R}_j^{b}\,|\lt q/2\end{cases} \tag{3.315} $$
则可得
$$ S_{mn}(q)=\delta_{mn} $$
也即在满足方程上述限制条件的前提下，$A_m$ 在 $\boldsymbol{k}$ 点网格上是正交的。与 Chadi-Cohen 方法类似，将函数 $f(\boldsymbol{k})$ 用 $A_m$ 展开，有
$$ f(\boldsymbol{k})=\sum_{m=1}f_mA_m(\boldsymbol{k}) \tag{3.316} $$
同时左乘 $A_m^{*}(\boldsymbol{k})$ 并在布里渊区内积分，可得
$$ f_m=\frac{\varOmega}{(2\pi)^{3}}\int_{\mathrm{BZ}}A_m^{*}(\boldsymbol{k})f(\boldsymbol{k})\,\mathrm{d}\boldsymbol{k} \tag{3.317} $$
因为 $A_1(\boldsymbol{k})=1$，所以由方程（3.317）可得
$$ f=\int_{\mathrm{BZ}}f(\boldsymbol{k})\,\mathrm{d}\boldsymbol{k}=\frac{8\pi^{3}}{\varOmega}f_1 \tag{3.318} $$
忽略前面的常数因子，可以看到 Monkhorst-Pack 方法中 $f$ 的表达式与 Chadi-Cohen 方法中的完全一样。
方程（3.318）虽然表明函数 $f(\boldsymbol{k})$ 的积分值可以用 $f_1$ 准确地给出，但是我们无法得到 $f_1$ 的精确值。因此仍然需要用上述 $\boldsymbol{k}$ 点网格得到 $f_1$，以及更普遍的 $f_m$ 的近似值 $\tilde{f}_m$：
$$ \tilde{f}_m=\frac{1}{q^{3}}\sum_{j=1}^{q^{3}}f(\boldsymbol{k}_j)A_m^{*}(\boldsymbol{k}_j) \tag{3.319} $$
相应地，函数 $f(\boldsymbol{k})$ 的近似值 $\tilde{f}(\boldsymbol{k})$ 可表示为
$$ \tilde{f}(\boldsymbol{k})=\sum_{m=1}\tilde{f}_mA_m(\boldsymbol{k}) \tag{3.320} $$
将恒等式\cite{gu1989solid}
$$ \frac{1}{(2\pi)^{3}}\int_{\mathrm{BZ}}f(\boldsymbol{k})=\lim_{V\to\infty}\frac{1}{V}\sum_{j}f(\boldsymbol{k}_j) $$
式（3.319）对完整网格求和，每点的权重已经由 $1/q^3$ 给出；若改用不可约布里渊区，则须采用式（3.324）的星形轨道重数。式（3.311）的 $A_m$ 对归一化布里渊区积分正交归一，而对未归一化的体积积分满足
$$ \int_{\mathrm{BZ}}A_m^{*}A_n\,\mathrm{d}\boldsymbol{k}=\frac{8\pi^{3}}{\varOmega}\delta_{mn} \tag{3.321} $$
利用这种 $\boldsymbol{k}$ 点网格近似布里渊区积分所产生的误差可按以下公式计算：
$$ \begin{aligned} \varepsilon_{\mathrm{BZ}} &=\int_{\mathrm{BZ}}[f(\boldsymbol k)-\tilde f(\boldsymbol k)]\,\mathrm d^3\boldsymbol k =\frac{(2\pi)^3}{\varOmega}(f_1-\tilde f_1)\\ &=-\frac{(2\pi)^3}{\varOmega} \sum_{m\gt 1}f_m\left[\frac1{q^3}\sum_{p,r,s=1}^q A_m(\boldsymbol k_{prs})\right] =-\frac{(2\pi)^3}{\varOmega}\sum_{m\gt 1}f_m S_{1m}(q). \end{aligned} \tag{3.322} $$
式中
$$ S_{1m}(q)=\frac{1}{\sqrt{N_m}} \sum_{|\boldsymbol R|=C_m}\prod_{j=1}^{3} \left[\frac1q\sum_{r=1}^{q} e^{\mathrm i\pi(2r-q-1)R_j/q}\right] \tag{3.323} $$
这里的 $S_{1m}(q)$ 与式（3.312）定义相同，因为 $A_1=1$。每个方括号可按式（3.314）计算；只有满足离散网格混叠条件的格矢才对误差有贡献。同一星形壳层内不能预设每个格矢的贡献都相同。
式（3.322）说明网格积分的误差由发生混叠的 Fourier 分量决定。对于足够平滑的周期函数，增加 $q$ 通常能改善积分，但收敛速度仍应以实际计算检验；费米面不连续处尤其不能仅凭此展开假定快速收敛。
但是根据方程（3.319）可知，$\tilde{f}_1$ 的计算量与 $q^{3}$ 成正比。如果 $q$ 值取得比较大，那么所需计算的 $\boldsymbol{k}$ 点数目就会非常大，如何提高 Monkhorst-Pack 方法的效率呢？如果考虑体系的对称性，则 $\boldsymbol{k}$ 点的数目会大大减少。重新写出 $f_1$ 的表达式如下：
$$ \tilde f_1=\frac{1}{q^{3}}\sum_{j=1}^{P(q)}\nu_jf(\boldsymbol{k}_j), \qquad\sum_{j=1}^{P(q)}\nu_j=q^3 \tag{3.324} $$
式中 $\nu_j=n_G/n_j$ 为代表点的星形轨道重数，$P(q)$ 为不可约网格中不等价点的数目。边界等价点只计一次，权重总和等于完整网格点数 $q^3$。高对称点的稳定子较大，因此轨道重数较小。文献\cite{monkhorst1976special}中给出了偶数 $q$ 时 BCC 和 FCC 两种格子的 $P(q)$：
对于 BCC 格子，
$$ P(q)=\begin{cases}q(q+4)(q+8)/192,&q/2\text{ 为偶数}\\(q+2)(q+4)(q+6)/192,&q/2\text{ 为奇数}\end{cases},\quad q\in2\mathbb N \tag{3.325} $$
对于 FCC 格子，
$$ P(q)=\begin{cases}q(q+2)(q+4)/96,&q/2\text{ 为偶数}\\(q+2)(q^{2}+4q+12)/96,&q/2\text{ 为奇数}\end{cases},\quad q\in2\mathbb N \tag{3.326} $$
可以看出，即使对于较大的 $q$ 值，$P(q)$ 也是比较小的，因此 Monkhorst-Pack 方法效率是比较高的。
Monkhorst–Pack 网格在倒格矢分数坐标中生成，适用于六角、单斜等晶格；选点及权重应符合晶格对称性\cite{chadi1977special}。以六角格子为例，Pack 指出 $\boldsymbol{k}$ 点网格应按以下公式生成\cite{pack1977special}：
$$ u_p=(p-1)/q_a,\quad u_r=(r-1)/q_a,\quad p,r=1,\ldots,q_a \tag{3.327} $$
$$ u_s=(2s-q_c-1)/2q_c,\quad s\in[1,q_c] \tag{3.328} $$
即 $a$ 轴和 $c$ 轴分别设置。相应地，$P(q)$ 的大小可按下式计算：
$$ P_a(q_a)=(\alpha+1)(3\alpha+\beta)+\delta_{\beta0},\quad \alpha=\lfloor q_a/6\rfloor,\quad\beta=q_a-6\alpha $$
$$ P_c(q_c)=\begin{cases}q_c/2,&q_c\text{ 为偶数}\\(q_c+1)/2,&q_c\text{ 为奇数}\end{cases} $$
对上述网格，面内分量包含零；沿 $c$ 方向只有 $q_c$ 为奇数时 $u_s=0$，故仅此时网格包含 $\varGamma$ 点。是否另加整体偏移，应由所选网格和对称性决定。
网格约化使用使积分函数不变的点群操作；若体系具有反演或镜面对称，也须一并考虑。倒格矢平移用于识别等价波矢。即使同属一种晶系，具体材料的对称性也可能不同，因此应按实际对称群确定权重。
Chadi-Cohen 方法的应用实例
1. $\boldsymbol{k}$ 点集合的生成
Cunningham\cite{cunningham1974special} 对于二维情况依照 Chadi-Cohen 方法分别生成了 $\boldsymbol{k}$ 点集合。我们选择长方格子和正方格子这两种情况进行具体的分析。
1）长方格子
实空间和倒空间的基矢及格点坐标分别为
$$ \begin{gathered} \boldsymbol{a}_1=a(1,0),\quad\boldsymbol{a}_2=a(0,\beta)\ (\text{其中}\ \beta\lt 1),\quad\boldsymbol{R}=a(l,n\beta)\\ \boldsymbol{b}_1=(2\pi/a)(1,0),\quad\boldsymbol{b}_2=(2\pi/a)(0,1/\beta),\quad\boldsymbol{K}=(2\pi/a)(k,n/\beta) \end{gathered} $$
选择
$$ \boldsymbol{k}_1^{0}=(\pi/a)[1/2,1/(2\beta)],\quad\boldsymbol{k}_2^{0}=(\pi/a)[1/4,1/(4\beta)] $$
前者保证 $l$ 或 $n$ 为奇数时 $A_m(\boldsymbol{k})=0$，而后者保证 $l/2$ 或 $n/2$ 为奇数时 $A_m(\boldsymbol{k})=0$。该长方格子的对称操作为 $\{E,c_2,\sigma_{\mathrm{v}}^{1},\sigma_{\mathrm{v}}^{2}\}$。按照 Chadi-Cohen 方法，可以构建 $\boldsymbol{k}_i$ 点如下：
$$ \begin{cases} \boldsymbol{k}_1=\boldsymbol{k}_1^{0}+E\boldsymbol{k}_2^{0}=[1/2,1/(2\beta)]+[1/4,1/(4\beta)]=[3/4,3/(4\beta)]\\ \boldsymbol{k}_2=\boldsymbol{k}_1^{0}+c_2\boldsymbol{k}_2^{0}=[1/2,1/(2\beta)]+[-1/4,-1/(4\beta)]=[1/4,1/(4\beta)]\\ \boldsymbol{k}_3=\boldsymbol{k}_1^{0}+\sigma_{\mathrm{v}}^{1}\boldsymbol{k}_2^{0}=[1/2,1/(2\beta)]+[-1/4,1/(4\beta)]=[1/4,3/(4\beta)]\\ \boldsymbol{k}_4=\boldsymbol{k}_1^{0}+\sigma_{\mathrm{v}}^{2}\boldsymbol{k}_2^{0}=[1/2,1/(2\beta)]+[1/4,-1/(4\beta)]=[3/4,1/(4\beta)] \end{cases} \tag{3.329} $$
每个 $\boldsymbol{k}$ 点的权重 $\alpha_i=1/4$。
2）正方格子
在上述情况下，令 $\beta=1$，则长方格子转变为正方格子。两种情况最主要的不同是布里渊区不可约部分有了变化。从式（3.329）可以看出，在正方格子中 $\beta=1$，$\boldsymbol{k}_3$ 和 $\boldsymbol{k}_4$ 重合。因此只有三个不同的 $\boldsymbol{k}$ 点，每个 $\boldsymbol{k}$ 点的权重分别为 $\alpha_1=\alpha_2=1/4$，$\alpha_3=1/2$，而且 $\displaystyle\sum_{i=1}^{3}\alpha_i=1$。
2. 利用特殊 $\boldsymbol{k}$ 点计算电荷密度
将 Bloch 函数用 Wannier 函数展开，有\cite{chadi1973electronic}
$$ \varPsi_{\boldsymbol{k}}(\boldsymbol{r})=\frac{1}{\sqrt{N}}\sum_{m}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{R}_m}a(\boldsymbol{r}-\boldsymbol{R}_m) \tag{3.330} $$
则在给定 $\boldsymbol{k}$ 点的电荷密度为
$$ \rho_{\boldsymbol{k}}(\boldsymbol{r})=\varPsi_{\boldsymbol{k}}^{*}(\boldsymbol{r})\varPsi_{\boldsymbol{k}}(\boldsymbol{r})=\frac{1}{N}\sum_{mn}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot(\boldsymbol{R}_m-\boldsymbol{R}_n)}a(\boldsymbol{r}-\boldsymbol{R}_m)a^{*}(\boldsymbol{r}-\boldsymbol{R}_n) \tag{3.331} $$
而
$$ \rho(\boldsymbol{r})=\int_{\mathrm{BZ}}\rho_{\boldsymbol{k}}(\boldsymbol{r})\,\mathrm{d}\boldsymbol{k} \tag{3.332} $$
将 $\rho_{\boldsymbol{k}}(\boldsymbol{r})$ 的表达式（3.331）写为
$$ \rho_{\boldsymbol{k}}(\boldsymbol{r})=\frac{1}{N}\sum_{m}|\,a(\boldsymbol{r}-\boldsymbol{R}_m)^{2}\,|+\frac{1}{N}\sum_{j}{}'\sum_{m}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{R}_j}a(\boldsymbol{r}-\boldsymbol{R}_m)a^{*}(\boldsymbol{r}+\boldsymbol{R}_j-\boldsymbol{R}_m) \tag{3.333} $$
式中：求和符号上的撇号（$'$）表明 $\boldsymbol{R}_j\neq\boldsymbol{0}$ 而且 $\boldsymbol{R}_j=\boldsymbol{R}_m-\boldsymbol{R}_n$。因此，考虑到对称性，$\rho_{\boldsymbol{k}}(\boldsymbol{r})$ 又可表示为
$$ \begin{aligned} \rho_{\boldsymbol{k}}(\boldsymbol{r})&=\frac{1}{n_G}\sum_{T_i}\rho_{T_i\boldsymbol{k}}(\boldsymbol{r})\\ &=\frac{1}{Nn_G}\sum_{T_i}\sum_{m}|\,a(\boldsymbol{r}-\boldsymbol{R}_m)\,|^{2}+\frac{1}{Nn_G}\sum_{j}{}'\sum_{m}\sum_{T_i}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{R}_j}a(\boldsymbol{r}-\boldsymbol{R}_m)a^{*}(\boldsymbol{r}+\boldsymbol{R}_j-\boldsymbol{R}_m) \end{aligned} \tag{3.334} $$
式（3.334）右端第一项与 $T_i$ 和 $\boldsymbol{k}$ 无关，相当于 Chadi-Cohen 方法中的 $f_0$，而第二项因为是对所有的 $j$ 求和，因此可以写成如下形式：
$$ \begin{aligned} F(\boldsymbol{r})&=\frac{1}{Nn_G}\sum_{j}{}'\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{R}_j}\sum_{m}\sum_{T_i}a(\boldsymbol{r}-\boldsymbol{R}_m)a^{*}(\boldsymbol{r}-T_i\boldsymbol{R}_j-\boldsymbol{R}_m)\\ &=\frac{1}{Nn_G}\sum_{j}{}'\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{R}_j}\sum_{m}S_m(\boldsymbol{r}) \end{aligned} \tag{3.335} $$
式中：$S_m(\boldsymbol{r})$ 与 $\boldsymbol{R}_j$ 无关，且随 $|\,T_i\boldsymbol{R}_j+\boldsymbol{R}_m\,|$ 的增大而减小，相当于 $f_m$。因此 $\rho_{\boldsymbol{k}}(\boldsymbol{r})$ 可写为
$$ \rho_{\boldsymbol{k}}(\boldsymbol{r})=f_0+\sum_{m}\sum_{|\boldsymbol{R}_j|=C_m}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{R}_j}f_m=f_0+\sum_{m}A_m(\boldsymbol{k})f_m \tag{3.336} $$
如果存在 $\boldsymbol{k}_0$，满足 $A_m(\boldsymbol{k}_0)=0$，$m=1,2,\cdots$，则可得
$$ \rho(\boldsymbol{r})=f_0=\frac{1}{N}\sum_m \bigl|a(\boldsymbol{r}-\boldsymbol{R}_m)\bigr|^{2} =\rho_{\boldsymbol{k}_0}(\boldsymbol{r}) \tag{3.337} $$
但是普遍来讲，这样的 $\boldsymbol{k}_0$ 并不存在。例如，在 FCC 格子中考虑第一、二、三近邻，写出 $A_m(\boldsymbol{k})$：
$$ \begin{cases} \cos k_x\cos k_y+\cos k_x\cos k_z+\cos k_y\cos k_z=0\\ \cos2k_x+\cos2k_y+\cos2k_z=0\\ \cos2k_x\cos k_y\cos k_z+\cos k_x\cos2k_y\cos k_z+\cos k_x\cos k_y\cos2k_z=0 \end{cases} \tag{3.338} $$
不存在单独的 $\boldsymbol{k}_0$ 点同时满足上述三个方程。因此，需要寻找一系列特殊的 $\boldsymbol{k}$ 点，满足
$$ \sum_{i=1}^{n}\sum_{|\boldsymbol{R}_j|=C_m}\alpha_i\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{R}_j}=\sum_{i=1}^{n}\alpha_iA_m(\boldsymbol{k}_i)=0,\quad\sum_{i}\alpha_i=1 \tag{3.339} $$
则 $\rho(\boldsymbol{r})=\displaystyle\sum_{i}\alpha_i\rho_{\boldsymbol{k}_i}(\boldsymbol{r})$。
Chadi 和 Cohen\cite{chadi1973electronic} 采用 $\boldsymbol{k}_1=(0.5,0,0)$、$\boldsymbol{k}_2=(1.0,0.5,0)$ 和 $\boldsymbol{k}_3=(0.5,0.5,0)$ 三个 $\boldsymbol{k}$ 点计算 $\rho(\boldsymbol{r})$ 取得了较好的结果：
$$ \rho(\boldsymbol{r})=\frac{1}{4}\rho_{\boldsymbol{k}_1}(\boldsymbol{r})+\frac{1}{2}\rho_{\boldsymbol{k}_2}(\boldsymbol{r})+\frac{1}{4}\rho_{\boldsymbol{k}_3}(\boldsymbol{r}) $$
而如果采用 $\boldsymbol{k}_1=(0.75,0.25,0.25)$ 和 $\boldsymbol{k}_2=(0.25,0.25,0.25)$，则可以改进计算结果：
$$ \rho(\boldsymbol{r})=\frac{3}{4}\rho_{\boldsymbol{k}_1}(\boldsymbol{r})+\frac{1}{4}\rho_{\boldsymbol{k}_2}(\boldsymbol{r}) $$
布里渊区积分——四面体方法
四面体方法是一种在第一性原理计算程序中广泛使用的积分方法，尤其适合用于计算总能量、态密度和各能带电子占据数等。与特殊 $\boldsymbol{k}$ 点法相比，四面体方法具有一定的优势和特点。
首先，四面体方法的基本原理是将布里渊区内的积分区域划分为许多小的四面体，然后在每个四面体内部进行积分。这种方法能够更加准确地描述能带结构的细节，特别是在处理金属或半金属体系时，能够有效地降低对 $\boldsymbol{k}$ 点采样密度的要求。
其次，四面体方法对于具有复杂能带结构的材料具有较高的计算效率。这是因为在四面体方法中，能带的贡献是通过局部积分得到的，而不是通过全局积分。因此，四面体方法能够减少不必要的计算量，从而提高计算速度。
最后，四面体方法在处理弱束缚和强关联体系时表现出较好的性能。这是由于四面体方法能够更好地捕捉能带结构中的局部特征，如能带交叉点和能带窄化现象，从而为研究这些体系提供一个更精确的描述。
以下我们将对四面体方法进行详细的讨论和分析。
总能量
四面体方法涉及各 $\boldsymbol{k}$ 点上的能级值，通常用 $E$ 表示。因此，为了避免符号混淆，此处的总能量用 $F$ 表示。则其期待值 $\langle F\rangle$ 在倒空间中计算如下：
$$ \langle F\rangle=\frac{1}{V_G}\sum_{n}\int_{V_G}\mathrm{d}^{3}k\,F_n(\boldsymbol{k})f[E_n(\boldsymbol{k})] \tag{3.340} $$
式中：$V_G$ 为第一布里渊区的体积；$f(\varepsilon)$ 是费米-狄拉克分布函数；$F_n(\boldsymbol{k})$ 为哈密顿算符 $\hat{H}$ 在第 $n$ 条能带中指定 $\boldsymbol{k}$ 点上的值。为了计算这个积分，四面体方法将第一布里渊区分解成若干个小的四面体，然后在各个四面体中分别计算积分，再对所有四面体求和，即
$$ \frac{1}{V_G}\int_{\mathrm{BZ}}F(\boldsymbol k)\,\mathrm d^3\boldsymbol k =\frac{1}{V_G}\sum_{T}\int_T F(\boldsymbol k)\,\mathrm d^3\boldsymbol k \tag{3.341} $$
而在每个四面体内 $F_n(\boldsymbol{k})$ 采用线性函数 $f(\boldsymbol{k})=a_0+a_1k_x+a_2k_y+a_3k_z$ 近似。$f(\boldsymbol{k})$ 的线性系数由边界条件 $f(\boldsymbol{k}_i)=F_{n,i}$（$i=1,2,3,4$）确定，$i$ 是四面体的四个顶点。这样，$F_n(\boldsymbol{k})$ 在每个四面体中的积分为
$$ \frac{1}{V_G}\int_{V_T}F(\boldsymbol{k})\,\mathrm{d}^{3}k=\frac{V_T}{V_G}\sum_{j=1}^{4}\frac{F_{n,j}}{4} \tag{3.342} $$
式中：$F_{n,j}$ 为四面体顶点的函数值。可以证明，式（3.342）中 1/4 这个权重来源于积分\cite{molenaar1982extended}
$$ \omega_i=\frac1{V_{T_0}}\int_{T_0}\lambda_i(\boldsymbol k)\,\mathrm d^3\boldsymbol k =\frac14,\qquad \sum_{i=1}^4\lambda_i=1,\quad V_{T_0}=\frac16 \tag{3.343} $$
将顶点分别为 $(0,0,0)$，$(1,0,0)$，$(0,1,0)$，$(0,0,1)$ 的参考四面体称为 $T_0$，$\lambda_i$ 为对应顶点的重心坐标；它们的平均值均为 $1/4$。式（3.342）适用于整个四面体均被占据的情形；若费米面穿过四面体，应在占据子区域上积分重心坐标以得到各顶点权重。式（3.342）也是现有算法中最基本的表达式。在不产生误解的前提下，此后的讨论中如无必要，将省略能带指数 $n$。
更加精确的计算表明\cite{molenaar1982extended,zaharioudakis2004tetrahedron}，利用线性函数近似 $F(\boldsymbol{k})$ 有些情况下并不是最理想的选择，这时对其更好的近似应该取为两个线性函数的商：
$$ F(\boldsymbol{k})\simeq\frac{p(\boldsymbol{k})}{g(\boldsymbol{k})} =\frac{a_0+a_1k_x+a_2k_y+a_3k_z} {b_0+b_1k_x+b_2k_y+b_3k_z} $$
$g(\boldsymbol{k})$ 的系数由 $g(\boldsymbol{k}_i)=E_i$（$i=1,2,3,4$）确定，其中 $E_i$ 是第 $i$ 个顶点上选定的分母函数值；为了匹配被积函数的顶点值，须取 $p(\boldsymbol{k}_i)=F_iE_i$。当分母取自能量时，必须明确能量零点及分母的物理来源。经推导可得，这时四面体的积分可表示为
$$ \frac1{V_G}\int_T F(\boldsymbol k)\,\mathrm d^3\boldsymbol k =\frac{V_T}{V_G}\sum_{i=1}^{4}F_i\omega_i \tag{3.344} $$
式中
$$ \omega_i=\frac{1}{\prod_{k\neq i}\left(1-\dfrac{E_k}{E_i}\right)}+\sum_{k\neq i}\frac{1}{\prod_{l\neq k}\left(1-\dfrac{E_l}{E_k}\right)}\frac{\ln\dfrac{E_k}{E_i}}{\dfrac{E_k}{E_i}-1} \tag{3.345} $$
式（3.345）及下列实对数表达式要求各 $E_i\ne0$，且在四面体内线性分母 $g(\boldsymbol{k})$ 不穿过零点；若四个顶点同号，此条件自动满足。若分母过零，这些权重不能作为普通实值体积分直接使用，应按具体响应函数的主值或复频率规定另行处理。
考虑到简并情况，Zaharioudakis 给出了比较完整的权重表达式\cite{zaharioudakis2004tetrahedron}：
（1）当 $E_1\lt E_2\lt E_3\lt E_4$ 时，同方程（3.345）；
（2）当 $E_1=E_2\lt E_3\lt E_4$ 时，
$$ \begin{aligned} \omega_1=\omega_2={}&\dfrac{5}{2\left(1-\dfrac{\mathrm{E}_{3}}{\mathrm{E}_{1}}\right)\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{1}}\right)}-\dfrac{1}{\left(1-\dfrac{\mathrm{E}_{3}}{\mathrm{E}_{1}}\right)^{2}\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{1}}\right)}-\dfrac{1}{\left(1-\dfrac{\mathrm{E}_{3}}{\mathrm{E}_{1}}\right)\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{1}}\right)^{2}}\\&+\dfrac{1}{\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{3}}\right)^{3}\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{3}}\right)}\dfrac{\ln\dfrac{\mathrm{E}_{3}}{\mathrm{E}_{1}}}{\dfrac{\mathrm{E}_{3}}{\mathrm{E}_{1}}}+\dfrac{1}{\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{4}}\right)^{3}\left(1-\dfrac{\mathrm{E}_{3}}{\mathrm{E}_{4}}\right)}\dfrac{\ln\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{1}}}{\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{1}}} \end{aligned} \tag{3.346} $$
$$ \begin{aligned} \omega_3={}&\dfrac{1}{\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{3}}\right)^{2}\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{3}}\right)}+\dfrac{\dfrac{\mathrm{E}_{3}}{\mathrm{E}_{1}}}{\left(1-\dfrac{\mathrm{E}_{3}}{\mathrm{E}_{1}}\right)^{2}\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{1}}\right)}+\dfrac{1}{\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{4}}\right)^{2}\left(1-\dfrac{\mathrm{E}_{3}}{\mathrm{E}_{4}}\right)^{2}}\dfrac{\ln\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{3}}}{\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{3}}}\\&+\left[\dfrac{3}{\left(1-\dfrac{\mathrm{E}_{3}}{\mathrm{E}_{1}}\right)^{2}\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{1}}\right)}-\dfrac{1}{\left(1-\dfrac{\mathrm{E}_{3}}{\mathrm{E}_{1}}\right)^{2}\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{1}}\right)^{2}}-\dfrac{2}{\left(1-\dfrac{\mathrm{E}_{3}}{\mathrm{E}_{1}}\right)^{3}\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{1}}\right)}\right]\dfrac{\ln\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{3}}}{\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{3}}} \end{aligned} \tag{3.347} $$
$\omega_4$ 与 $\omega_3$ 形式相同，只需要将式（3.347）中的 $\mathrm{E}_3$ 和 $\mathrm{E}_4$ 位置互换即可。
（3）当 $E_1\lt E_2=E_3\lt E_4$ 时，
$$ \begin{aligned} \omega_1={}&\dfrac{1}{\left(1-\dfrac{\mathrm{E}_{2}}{\mathrm{E}_{1}}\right)^{2}\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{1}}\right)}+\dfrac{\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{2}}}{\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{2}}\right)^{2}\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{2}}\right)}+\dfrac{1}{\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{4}}\right)^{2}\left(1-\dfrac{\mathrm{E}_{2}}{\mathrm{E}_{4}}\right)^{2}}\dfrac{\ln\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{1}}}{\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{1}}}\\&+\left[\dfrac{3}{\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{2}}\right)^{2}\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{2}}\right)}-\dfrac{1}{\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{2}}\right)^{2}\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{2}}\right)^{2}}-\dfrac{2}{\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{2}}\right)^{3}\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{2}}\right)}\right]\dfrac{\ln\dfrac{\mathrm{E}_{2}}{\mathrm{E}_{1}}}{\dfrac{\mathrm{E}_{2}}{\mathrm{E}_{1}}} \end{aligned} \tag{3.348} $$
$$ \begin{aligned} \omega_2=\omega_3={}&\dfrac{5}{2\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{2}}\right)\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{2}}\right)}-\dfrac{1}{\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{2}}\right)^{2}\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{2}}\right)}-\dfrac{1}{\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{2}}\right)\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{2}}\right)^{2}}\\&+\dfrac{1}{\left(1-\dfrac{\mathrm{E}_{2}}{\mathrm{E}_{1}}\right)^{3}\left(1-\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{1}}\right)}\dfrac{\ln\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{2}}}{\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{2}}}+\dfrac{1}{\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{4}}\right)\left(1-\dfrac{\mathrm{E}_{2}}{\mathrm{E}_{4}}\right)^{3}}\dfrac{\ln\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{2}}}{\dfrac{\mathrm{E}_{4}}{\mathrm{E}_{2}}} \end{aligned} \tag{3.349} $$
$\omega_4$ 与 $\omega_1$ 形式相同，只需要将式（3.348）中的 $\mathrm{E}_1$ 和 $\mathrm{E}_4$ 位置互换即可。
（4）当 $E_1\lt E_2\lt E_3=E_4$ 时，
$$ \begin{aligned} \omega_1={}&\dfrac{1}{\left(1-\dfrac{\mathrm{E}_{3}}{\mathrm{E}_{1}}\right)^{2}\left(1-\dfrac{\mathrm{E}_{2}}{\mathrm{E}_{1}}\right)}+\dfrac{\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{3}}}{\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{3}}\right)^{2}\left(1-\dfrac{\mathrm{E}_{2}}{\mathrm{E}_{3}}\right)}+\dfrac{1}{\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{2}}\right)^{2}\left(1-\dfrac{\mathrm{E}_{3}}{\mathrm{E}_{2}}\right)^{2}}\dfrac{\ln\dfrac{\mathrm{E}_{2}}{\mathrm{E}_{1}}}{\dfrac{\mathrm{E}_{2}}{\mathrm{E}_{1}}}\\&+\left[\dfrac{3}{\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{3}}\right)^{2}\left(1-\dfrac{\mathrm{E}_{2}}{\mathrm{E}_{3}}\right)}-\dfrac{1}{\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{3}}\right)^{2}\left(1-\dfrac{\mathrm{E}_{2}}{\mathrm{E}_{3}}\right)^{2}}-\dfrac{2}{\left(1-\dfrac{\mathrm{E}_{1}}{\mathrm{E}_{3}}\right)^{3}\left(1-\dfrac{\mathrm{E}_{2}}{\mathrm{E}_{3}}\right)}\right]\dfrac{\ln\dfrac{\mathrm{E}_{3}}{\mathrm{E}_{1}}}{\dfrac{\mathrm{E}_{3}}{\mathrm{E}_{1}}} \end{aligned} \tag{3.350} $$
$$ \begin{aligned} \omega_3=\omega_4={}&\dfrac{5}{2\left(1-\dfrac{E_{1}}{E_{3}}\right)\left(1-\dfrac{E_{2}}{E_{3}}\right)}-\dfrac{1}{\left(1-\dfrac{E_{1}}{E_{3}}\right)^{2}\left(1-\dfrac{E_{2}}{E_{3}}\right)}-\dfrac{1}{\left(1-\dfrac{E_{1}}{E_{3}}\right)\left(1-\dfrac{E_{2}}{E_{3}}\right)^{2}}\\&+\dfrac{1}{\left(1-\dfrac{E_{2}}{E_{1}}\right)\left(1-\dfrac{E_{3}}{E_{1}}\right)^{3}}\dfrac{\ln\dfrac{E_{1}}{E_{3}}}{\dfrac{E_{1}}{E_{3}}}+\dfrac{1}{\left(1-\dfrac{E_{1}}{E_{2}}\right)\left(1-\dfrac{E_{3}}{E_{2}}\right)^{3}}\dfrac{\ln\dfrac{E_{2}}{E_{3}}}{\dfrac{E_{2}}{E_{3}}} \end{aligned} \tag{3.351} $$
$\omega_2$ 与 $\omega_1$ 形式相同，只需要将式（3.350）中的 $E_1$ 和 $E_2$ 位置互换即可。
（5）当 $E_1=E_2=E_3\lt E_4$ 时，
$$ \omega_1=\omega_2=\omega_3=\dfrac{11}{6\left(1-\dfrac{E_{4}}{E_{1}}\right)}-\dfrac{5}{2\left(1-\dfrac{E_{4}}{E_{1}}\right)^{2}}+\dfrac{1}{\left(1-\dfrac{E_{4}}{E_{1}}\right)^{3}}+\dfrac{1}{\left(1-\dfrac{E_{1}}{E_{4}}\right)^{4}}\dfrac{\ln\dfrac{E_{4}}{E_{1}}}{\dfrac{E_{4}}{E_{1}}} \tag{3.352} $$
$$ \begin{aligned} \omega_4={}&\left[\dfrac{3}{\left(1-\dfrac{E_{4}}{E_{1}}\right)^{2}}-\dfrac{6}{\left(1-\dfrac{E_{4}}{E_{1}}\right)^{3}}+\dfrac{3}{\left(1-\dfrac{E_{4}}{E_{1}}\right)^{4}}\right]\dfrac{\ln\dfrac{E_{1}}{E_{4}}}{\dfrac{E_{1}}{E_{4}}}+\dfrac{1}{\left(1-\dfrac{E_{1}}{E_{4}}\right)^{3}}\\&+\dfrac{5}{2}\dfrac{\dfrac{E_{4}}{E_{1}}}{\left(1-\dfrac{E_{4}}{E_{1}}\right)^{2}}-2\dfrac{\dfrac{E_{4}}{E_{1}}}{\left(1-\dfrac{E_{4}}{E_{1}}\right)^{3}} \end{aligned} \tag{3.353} $$
（6）当 $E_1\lt E_2=E_3=E_4$ 时，
$$ \begin{aligned} \omega_1={}&\left[\dfrac{3}{\left(1-\dfrac{E_{1}}{E_{2}}\right)^{2}}-\dfrac{6}{\left(1-\dfrac{E_{1}}{E_{2}}\right)^{3}}+\dfrac{3}{\left(1-\dfrac{E_{1}}{E_{2}}\right)^{4}}\right]\dfrac{\ln\dfrac{E_{2}}{E_{1}}}{\dfrac{E_{2}}{E_{1}}}+\dfrac{1}{\left(1-\dfrac{E_{2}}{E_{1}}\right)^{3}}\\&+\dfrac{5}{2}\dfrac{\dfrac{E_{1}}{E_{2}}}{\left(1-\dfrac{E_{1}}{E_{2}}\right)^{2}}-2\dfrac{\dfrac{E_{1}}{E_{2}}}{\left(1-\dfrac{E_{1}}{E_{2}}\right)^{3}} \end{aligned} \tag{3.354} $$
$$ \omega_2=\omega_3=\omega_4=\dfrac{11}{6\left(1-\dfrac{E_{1}}{E_{2}}\right)}-\dfrac{5}{2\left(1-\dfrac{E_{1}}{E_{2}}\right)^{2}}+\dfrac{1}{\left(1-\dfrac{E_{1}}{E_{2}}\right)^{3}}+\dfrac{1}{\left(1-\dfrac{E_{2}}{E_{1}}\right)^{4}}\dfrac{\ln\dfrac{E_{1}}{E_{2}}}{\dfrac{E_{1}}{E_{2}}} \tag{3.355} $$
（7）当 $E_1=E_2\lt E_3=E_4$ 时，
$$ \omega_1=\omega_2=\dfrac{5}{2\left(1-\dfrac{E_{3}}{E_{1}}\right)^{2}}-\dfrac{2}{\left(1-\dfrac{E_{3}}{E_{1}}\right)^{3}}+\dfrac{\dfrac{E_{1}}{E_{3}}}{\left(1-\dfrac{E_{1}}{E_{3}}\right)^{3}}+\left[\dfrac{3}{\left(1-\dfrac{E_{1}}{E_{3}}\right)^{3}}-\dfrac{3}{\left(1-\dfrac{E_{1}}{E_{3}}\right)^{4}}\right]\dfrac{\ln\dfrac{E_{3}}{E_{1}}}{\dfrac{E_{3}}{E_{1}}} \tag{3.356} $$
$$ \omega_3=\omega_4=\frac{5}{2\left(1-\dfrac{E_1}{E_3}\right)^{2}}-\frac{2}{\left(1-\dfrac{E_1}{E_3}\right)^{3}}+\frac{\dfrac{E_3}{E_1}}{\left(1-\dfrac{E_3}{E_1}\right)^{3}}+\left[\frac{3}{\left(1-\dfrac{E_3}{E_1}\right)^{3}}-\frac{3}{\left(1-\dfrac{E_3}{E_1}\right)^{4}}\right]\frac{\ln\dfrac{E_1}{E_3}}{\dfrac{E_1}{E_3}} \tag{3.357} $$
（8）当 $E_1=E_2=E_3=E_4$ 时，
$$ \omega_1=\omega_2=\omega_3=\omega_4=\frac{1}{4} \tag{3.358} $$
态密度
四面体方法的提出，最早就是为了求解形如
$$ I(E)=\int_{E(\boldsymbol{k})=E}F(\boldsymbol{k})\,|\,\boldsymbol{\nabla}E(\boldsymbol{k})\,|^{-1}\,\mathrm{d}S \tag{3.359} $$
的积分式。当 $F(\boldsymbol{k})\equiv1$ 时，积分式前乘以因子 $1/V_G$，式（3.359）就成为态密度的定义式。除了在物理上的重要性以外，讨论态密度有助于直观地理解四面体方法的几何意义。
Lehmann 和 Taut\cite{lehmann1972numerical} 指出，设四个顶点的能量本征值满足条件 $E_1\lt E_2\lt E_3\lt E_4$，在四面体内能量按照线性函数展开，与 3.4.2.1 节的 $g(\boldsymbol{k})$ 相同：
$$ E(\boldsymbol{k})=E_1+\boldsymbol{b}\cdot(\boldsymbol{k}-\boldsymbol{k}_1) \tag{3.360} $$
式中：$\boldsymbol{b}=\displaystyle\sum_{i=1}^{3}(E_{i+1}-E_1)\boldsymbol{r}_i$，$\boldsymbol{r}_i\cdot\boldsymbol{k}_j=\delta_{ij}$，其中 $\boldsymbol{k}_j=\boldsymbol{k}_{j+1}-\boldsymbol{k}_1$。因此可以按照倒格矢与正格矢的关系式写出 $\boldsymbol{r}_i$，有
$$ \boldsymbol{r}_1=\frac{\boldsymbol{k}_2\times\boldsymbol{k}_3}{\boldsymbol{k}_1\cdot(\boldsymbol{k}_2\times\boldsymbol{k}_3)},\quad\boldsymbol{r}_2=\frac{\boldsymbol{k}_3\times\boldsymbol{k}_1}{\boldsymbol{k}_1\cdot(\boldsymbol{k}_2\times\boldsymbol{k}_3)},\quad\boldsymbol{r}_3=\frac{\boldsymbol{k}_1\times\boldsymbol{k}_2}{\boldsymbol{k}_1\cdot(\boldsymbol{k}_2\times\boldsymbol{k}_3)} $$
式中分母是以 $\boldsymbol{k}_i$ 为顶点的四面体体积的六倍，即 $6V_T$。因此该四面体对态密度 $D_T(\mathrm{E})$ 的贡献为
$$ D_T(E)=\frac{1}{V_G}\frac{\mathrm{dS}(E)}{|\,\boldsymbol{b}\,|}=\begin{cases}0,&E\leqslant E_1\\\dfrac{1}{V_G}\dfrac{f_1}{|\,\boldsymbol{b}\,|},&E_1\leqslant E\leqslant E_2\\\dfrac{1}{V_G}\dfrac{f_1-f_2}{|\,\boldsymbol{b}\,|},&E_2\leqslant E\leqslant E_3\\\dfrac{1}{V_G}\dfrac{f_4}{|\,\boldsymbol{b}\,|},&E_3\leqslant E\leqslant E_4\\0,&E\geqslant E_4\end{cases} \tag{3.361} $$
图 3.9 等能面 S(E) 在四面体中的截面
其中函数 $f$ 是等能面 $S(E)$ 在四面体内的截面面积，如图 3.9 所示。因此容易得出 $f/|\boldsymbol{b}|$ 的表达式为
$$ \begin{cases} \dfrac{f_1}{|\boldsymbol{b}|}=3V_T\dfrac{(E-E_1)^{2}}{(E_2-E_1)(E_3-E_1)(E_4-E_1)}\\ \dfrac{f_2}{|\boldsymbol{b}|}=3V_T\dfrac{(E-E_2)^{2}}{(E_2-E_1)(E_3-E_2)(E_4-E_2)}\\ \dfrac{f_4}{|\boldsymbol{b}|}=3V_T\dfrac{(E-E_4)^{2}}{(E_4-E_1)(E_4-E_2)(E_4-E_3)} \end{cases} \tag{3.362} $$
对于态密度，Jepsen 和 Andersen\cite{jepsen1971electronic} 还提出了更直观简单的计算方法，之后由 Blöchl 对他们所提方法进行了改进\cite{blochl1994improved}。可以将每个四面体看作容器，给定等能面 $S(E)$ 之后，该四面体对电子态数目 $n(E)$ 的贡献等于 $\varepsilon\leqslant E$ 的包络体积对第一布里渊区体积的比值。而态密度 $D_T(E)$ 可以定义为 $\mathrm{d}n/\mathrm{d}E$。通过简单的三角锥体积计算，可得
$$ n_T(E)=\frac{V_T}{V_G}\begin{cases} 0,&E\le E_1,\\ \dfrac{(E-E_1)^3}{(E_2-E_1)(E_3-E_1)(E_4-E_1)},&E_1\lt E\le E_2,\\ \dfrac{(E-E_1)^3}{(E_2-E_1)(E_3-E_1)(E_4-E_1)} -\dfrac{(E-E_2)^3}{(E_2-E_1)(E_3-E_2)(E_4-E_2)},&E_2\lt E\le E_3,\\ 1-\dfrac{(E_4-E)^3}{(E_4-E_1)(E_4-E_2)(E_4-E_3)},&E_3\lt E\le E_4,\\ 1,&E\gt E_4. \end{cases} \tag{3.363} $$
将式（3.363）对 $E$ 求导，即得方程（3.361）。对于简并情况，式（3.361）至式（3.363）并不适用。因此对于式（3.363）中的第三种情况，电子态数目方程可等效地写为\cite{jepsen1971electronic}
$$ \begin{aligned} n(E)={}&\frac{V_T}{V_G}\frac{1}{(E_3-E_1)(E_4-E_1)}\left[(E_2-E_1)^{2}+3(E_2-E_1)(E-E_2)+3(E-E_2)^{2}\right.\\ &\left.-\frac{E_3-E_1+E_4-E_2}{(E_3-E_2)(E_4-E_2)}(E-E_2)^{3}\right],\quad E_2\leqslant E\leqslant E_3 \end{aligned} \tag{3.364} $$
而相应地，该情况下的态密度 $D_T(E)$ 可写为
$$ \begin{aligned} D_T(E)={}&\frac{V_T}{V_G}\frac{1}{(E_3-E_1)(E_4-E_1)}\left[3(E_2-E_1)+6(E-E_2)\right.\\ &\left.-3\frac{(E_3-E_1+E_4-E_2)(E-E_2)^{2}}{(E_3-E_2)(E_4-E_2)}\right],\quad E_2\leqslant E\leqslant E_3 \end{aligned} $$
在传统的四面体方法中，首先通过对称群的操作找出第一布里渊区的不可约部分，然后在等间距的 $\boldsymbol{k}$ 点网格上将其手动划分为若干四面体。这种划分有一定的随意性，因为给定一组 $\boldsymbol{k}$ 点网格，可以有很多种不同的方法将该网格划分为若干互不重叠的四面体，进而进行计算。这种随意性是否会造成计算上的误差甚至错误？Kleinman 对此做了讨论\cite{kleinman1983error}。通过一个简单的例子，他证明在 $\boldsymbol{k}$ 点比较稀疏的情况下，不同的划分会造成高达 14% 的误差，而划分正确时会得出与特殊 $\boldsymbol{k}$ 点法相同的结果。在 $\boldsymbol{k}$ 点足够密集时，这种划分上的随意性对计算结果的影响可以忽略不计。其原因在于，在四面体方法中，除个别 $\boldsymbol{k}$ 点外，绝大多数 $\boldsymbol{k}$ 点由多个四面体共享，因此实际上在积分中参与了超过一次的计算。边界及靠近边界的 $\boldsymbol{k}$ 点在不同的划分中所归属的四面体数量不尽相同，从而使得各自对积分值的贡献有差异。而在网格中间区域的 $\boldsymbol{k}$ 点则不存在这个问题。因此，在 $\boldsymbol{k}$ 点密集时，可以忽略边界处 $\boldsymbol{k}$ 点的贡献，但在 $\boldsymbol{k}$ 点较少时，这样做则可能影响最终结果。Kleinman 的研究促使 Jepsen 和 Andersen 重新审视四面体方法，并提出了一种通用的、适用于编程的四面体划分法\cite{jepsen1984error}。这为后来 Blöchl 的工作奠定了基础，进一步推动了四面体方法在第一性原理计算中的应用和发展。
Kleinman 的工作的重要性还在于第一次明确指出了利用四面体方法，同样可以将 $\langle F\rangle$ $=\dfrac{1}{V_G}\displaystyle\sum_{n}\int_{V_G}\mathrm{d}^{3}k\,F_n(\boldsymbol{k})f(E_n(\boldsymbol{k}))$ 中的积分计算转化为各个 $\boldsymbol{k}$ 点上的被积函数值的加权求和：
$$ \langle F\rangle\simeq\sum_n\sum_i F_n(\boldsymbol k_i) \sum_{T\ni\boldsymbol k_i}w_{T,i}^{(n)}(E_F),\qquad w_{T,i}^{(n)}(E_F)=\frac1{V_G} \int_{T\cap\{E_n(\boldsymbol k)\le E_F\}}\lambda_{T,i}(\boldsymbol k)\, \mathrm d^3\boldsymbol k \tag{3.365} $$
内层求和遍历以 $\boldsymbol{k}_i$ 为顶点的四面体；$\lambda_{T,i}$ 为该顶点的重心坐标函数。完整占据的四面体满足 $w_{T,i}=V_T/(4V_G)$，但部分占据时各顶点的权重一般不同，不能统一用 $V_T^{\mathrm{occ}}/(4V_G)$ 代替。占据体积可由式（3.363）求得。在前面工作的基础上，Blöchl 给出了有较大改进的四面体方法的普适算法\cite{blochl1994improved}，其主要特点体现在以下三个方面：四面体的自动划分；各 $\boldsymbol{k}$ 点的权重；对金属体系的 Blöchl 修正。
1. 四面体的自动划分
前面已经指出，为了减少需要计算的四面体数目，首先需要找出第一布里渊区的不可约部分。这样的策略有一个副作用，即对不可约部分的四面体划分几乎不可避免地要进行人工干预，不利于编程求解。因此 Blöchl 提出，首先利用 Monkhorst-Pack（MP）方法\cite{monkhorst1976special}在第一布里渊区内生成等距的 $\boldsymbol{k}$ 点网格，然后给每个 $\boldsymbol{k}$ 点编号：
$$ \mathrm{N}=1+\frac{l-l_0}{2}+(n_1+1)\left[\frac{m-m_0}{2}+(n_2+1)\frac{n-n_0}{2}\right] \tag{3.366} $$
式中：$(l,m,n)$ 是该 $\boldsymbol{k}$ 点沿倒格矢 $\boldsymbol{b}_1$、$\boldsymbol{b}_2$、$\boldsymbol{b}_3$ 的序数的 2 倍；$n_i$ 是在三个方向上的 $\boldsymbol{k}$ 点数；$(l_0,m_0,n_0)$ 是 MP 方法中 $\varGamma$ 点的偏移量，有偏移则为 1，否则为 0。编号之后建立标识数组，其位置与该位置储存的元素值相同，例如，在第一个位置存储 1，在第二个位置存储 2，依次类推。然后从第一个位置开始，利用对称群的操作矩阵对每个 $\boldsymbol{k}$ 点坐标进行操作，再与数组中其他 $\boldsymbol{k}$ 点的坐标进行比较，如果彼此相同且后者的编号大于前者，则将后者的元素值改为前者的。这样对全部数组操作完毕之后，可以立即挑出所有不可约 $\boldsymbol{k}$ 点——只有当 $\boldsymbol{k}_i$ 点为不可约 $\boldsymbol{k}$ 点时，其编号才与其存储位置相同。之后，为了计算方便，可以对所有这些不可约 $\boldsymbol{k}$ 点按存储位置的顺序重新编号，即从 1 到 $\boldsymbol{k}_{\mathrm{irr}}(\max)$。数组中的各个元素也相应地改为新的编号（或名称）。这样整个第一布里渊区中的 $\boldsymbol{k}$ 点都可用不可约 $\boldsymbol{k}$ 点标记。
下一步讨论四面体的自动划分过程。以下八组坐标代表的 $\boldsymbol{k}$ 点构成平行六面体：
$$ \begin{gathered} (l,m,n)-1,\quad(l+2,m,n)-2,\quad(l,m+2,n)-3,\quad(l,m,n+2)-4,\\ (l+2,m+2,n)-5,\ (l+2,m,n+2)-6,\ (l,m+2,n+2)-7,\ (l+2,m+2,n+2)-8 \end{gathered} $$
为了尽量减小插值引起的误差，可以取此平行六面体中最短的体对角线作为等体积的六个四面体的公共对角线。将平行六面体的顶点依次设为 1～8，则可以采用下面六组途径确定这六个四面体的各个顶点：
$$ \begin{gathered} 3\to5\to6\to7,\quad3\to6\to7\to8,\quad3\to4\to6\to8\\ 1\to3\to5\to6,\quad1\to2\to3\to6,\quad2\to3\to4\to6 \end{gathered} $$
整个分解过程如图 3.10 所示。为简单起见，图中以立方体为例。对每个平行六面体重复上述过程，可以将整个第一布里渊区划分为体积相等的若干个四面体，每个四面体的顶点可用标识数组中的不可约 $\boldsymbol{k}$ 点标记。将这四个顶点的标号按升序排列，则每个四面体的简并度可以轻易得出。因此，这个过程保证了可以只用不可约 $\boldsymbol{k}$ 点上的信息进行整个第一布里渊区的积分，而无须考虑如何划定其不可约部分。上述过程可以避免 Kleinman 所说的计算误差，而且整个过程可以通过程序自动实现而无须人工干预。
2. 各个 $\boldsymbol{k}$ 点的权重计算
四面体方法的基本过程是，先算出每个四面体对积分值的贡献，再对所有四面体求和。
图 3.10 k 点网格中的四面体划分
因此如前文所述，绝大多数 $\boldsymbol{k}$ 点参与了多次计算。同时采用这种算法时也需要知道每个不可约 $\boldsymbol{k}$ 点上的被积函数值，对大规模计算而言难免会有存储方面的困难。如同 Kleinman 所指出的，积分 $\langle F\rangle$ 可以表示为 $\boldsymbol{k}$ 点的加权求和，即
$$ \langle F\rangle=\sum_{i,n}F_n(\boldsymbol{k}_i)\omega_{ni} \tag{3.367} $$
与标准的四面体方法相同，在每个四面体内，$F_n(\boldsymbol{k})$ 用线性函数 $f_n(\boldsymbol{k})$ 近似。$f_n(\boldsymbol{k})$ 也可写为
$$ f_n(\boldsymbol{k})=\sum_{i\in T}F_n(\boldsymbol{k}_i)\lambda_{T,i}(\boldsymbol{k}), \qquad\boldsymbol{k}\in T \tag{3.368} $$
将其代入方程（3.367），可得每个不可约 $\boldsymbol{k}$ 点的积分权重 $\omega_{ni}$：
$$ \omega_{ni}=\frac{1}{V_G}\sum_{T\ni\boldsymbol k_i} \int_{T\cap\{E_n(\boldsymbol k)\le E_F\}} \lambda_{T,i}(\boldsymbol k)\,\mathrm d^3\boldsymbol k \tag{3.369} $$
式（3.368）的 $\lambda_{T,i}(\boldsymbol k)$ 是四面体内的线性重心坐标，顶点 $i$ 处为 1，其余三个顶点处为 0。将其代入式（3.369），即得到式（3.365）的顶点权重；若把部分占据四面体内的权重函数都取为常数 $1/4$，得到的只是额外的常数权重近似。为计算线性插值下的权重 $\omega_{ni}$，首先必须计算体系的费米能 $E_{\mathrm{F}}$：利用方程（3.363）计算能量小于给定 $E$ 的状态数 $n(E)$，直至填充电子数与体系总电子数相符，此时的 $E$ 即为 $E_{\mathrm{F}}$。据此可以得出各四面体的顶点权重（下文省略能带指数 $n$）：
当 $E_{\mathrm{F}}\lt E_1$ 时，
$$ \omega_1=\omega_2=\omega_3=\omega_4=0 \tag{3.370} $$
当 $E_1\lt E_{\mathrm{F}}\lt E_2$ 时
$$ \begin{cases} \omega_1=C_0\left[4-(E_{\mathrm{F}}-E_1)\left(\dfrac{1}{E_2-E_1}+\dfrac{1}{E_3-E_1}+\dfrac{1}{E_4-E_1}\right)\right]\\ \omega_2=C_0\dfrac{E_{\mathrm{F}}-E_1}{E_2-E_1}\\ \omega_3=C_0\dfrac{E_{\mathrm{F}}-E_1}{E_3-E_1}\\ \omega_4=C_0\dfrac{E_{\mathrm{F}}-E_1}{E_4-E_1} \end{cases} \tag{3.371} $$
式中
$$ C_0=\frac{V_T}{4V_G}\frac{(\mathrm{E}_{\mathrm{F}}-\mathrm{E}_1)^{3}}{(\mathrm{E}_2-\mathrm{E}_1)(\mathrm{E}_3-\mathrm{E}_1)(\mathrm{E}_4-\mathrm{E}_1)} $$
当 $E_2\lt E_{\mathrm{F}}\lt E_3$ 时
$$ \begin{cases} \omega_1=\mathrm{C}_1+(\mathrm{C}_1+\mathrm{C}_2)\dfrac{\mathrm{E}_3-\mathrm{E}_{\mathrm{F}}}{\mathrm{E}_3-\mathrm{E}_1}+(\mathrm{C}_1+\mathrm{C}_2+\mathrm{C}_3)\dfrac{\mathrm{E}_4-\mathrm{E}_{\mathrm{F}}}{\mathrm{E}_4-\mathrm{E}_1}\\ \omega_2=\mathrm{C}_1+\mathrm{C}_2+\mathrm{C}_3+(\mathrm{C}_2+\mathrm{C}_3)\dfrac{\mathrm{E}_3-\mathrm{E}_{\mathrm{F}}}{\mathrm{E}_3-\mathrm{E}_2}+\mathrm{C}_3\dfrac{\mathrm{E}_4-\mathrm{E}_{\mathrm{F}}}{\mathrm{E}_4-\mathrm{E}_2}\\ \omega_3=(\mathrm{C}_1+\mathrm{C}_2)\dfrac{\mathrm{E}_{\mathrm{F}}-\mathrm{E}_1}{\mathrm{E}_3-\mathrm{E}_1}+(\mathrm{C}_2+\mathrm{C}_3)\dfrac{\mathrm{E}_{\mathrm{F}}-\mathrm{E}_2}{\mathrm{E}_3-\mathrm{E}_2}\\ \omega_4=(\mathrm{C}_1+\mathrm{C}_2+\mathrm{C}_3)\dfrac{\mathrm{E}_{\mathrm{F}}-\mathrm{E}_1}{\mathrm{E}_4-\mathrm{E}_1}+\mathrm{C}_3\dfrac{\mathrm{E}_{\mathrm{F}}-\mathrm{E}_2}{\mathrm{E}_4-\mathrm{E}_2} \end{cases} \tag{3.372} $$
式中
$$ \begin{cases} \mathrm{C}_1=\dfrac{V_T}{4V_G}\dfrac{(\mathrm{E}_{\mathrm{F}}-\mathrm{E}_1)^{2}}{(\mathrm{E}_4-\mathrm{E}_1)(\mathrm{E}_3-\mathrm{E}_1)}\\ \mathrm{C}_2=\dfrac{V_T}{4V_G}\dfrac{(\mathrm{E}_{\mathrm{F}}-\mathrm{E}_1)(\mathrm{E}_{\mathrm{F}}-\mathrm{E}_2)(\mathrm{E}_3-\mathrm{E}_{\mathrm{F}})}{(\mathrm{E}_4-\mathrm{E}_1)(\mathrm{E}_3-\mathrm{E}_2)(\mathrm{E}_3-\mathrm{E}_1)}\\ \mathrm{C}_3=\dfrac{V_T}{4V_G}\dfrac{(\mathrm{E}_{\mathrm{F}}-\mathrm{E}_2)^{2}(\mathrm{E}_4-\mathrm{E}_{\mathrm{F}})}{(\mathrm{E}_4-\mathrm{E}_2)(\mathrm{E}_3-\mathrm{E}_2)(\mathrm{E}_4-\mathrm{E}_1)} \end{cases} $$
当 $E_3\lt E_{\mathrm{F}}\lt E_4$ 时，
$$ \begin{cases} \omega_1=\dfrac{V_T}{4V_G}-C_4\dfrac{E_4-E_{\mathrm{F}}}{E_4-E_1}\\ \omega_2=\dfrac{V_T}{4V_G}-C_4\dfrac{E_4-E_{\mathrm{F}}}{E_4-E_2}\\ \omega_3=\dfrac{V_T}{4V_G}-C_4\dfrac{E_4-E_{\mathrm{F}}}{E_4-E_3}\\ \omega_4=\dfrac{V_T}{4V_G}-C_4\left[4-\left(\dfrac{1}{E_4-E_1}+\dfrac{1}{E_4-E_2}+\dfrac{1}{E_4-E_3}\right)(E_4-E_{\mathrm{F}})\right] \end{cases} \tag{3.373} $$
式中
$$ C_4=\frac{V_T}{4V_G}\frac{(E_4-E_{\mathrm{F}})^{3}}{(E_4-E_1)(E_4-E_2)(E_4-E_3)} $$
当 $E_{\mathrm{F}}\gt E_4$ 时，
$$ \omega_1=\omega_2=\omega_3=\omega_4=\frac{V_T}{4V_G},\quad E_{\mathrm{F}}\gt E_4 \tag{3.374} $$
因此，对于给定的 $\boldsymbol{k}_i$，可以找出其自身和等价点所属的四面体，通过分析各四面体内的能级分布情况，由式（3.370）至式（3.374）计算出积分权重 $\omega_{ni}$。
3. Blöchl 修正
在四面体方法中，通过对被积函数进行线性插值，可以得到简单的计算公式，同时可在很大程度上保持计算精度。然而，对于费米面形状复杂且部分填充的过渡金属体系，我们仍然需要考虑线性近似可能带来的误差以及相关的修正问题。对于这些系统的精确计算和理解，需要在计算方法和数值技巧上做出适当的调整，以减小近似所产生的影响并提高计算精度。
首先讨论误差。在四面体中，真实的被积函数 $F(\boldsymbol{k})$ 总会有正曲率和负曲率的部分，而线性函数 $f(\boldsymbol{k})$ 则会高估正曲率区间的函数值，同时低估负曲率区间的函数值。对于绝缘体和半导体，由于积分区域覆盖整个四面体，这两部分误差在很大程度上可以相互抵消。但是，对于金属体系，由于导带部分填充，即积分区域仅覆盖四面体的一部分而非全部，高估和低估的误差部分之间可能存在显著的不平衡，这将导致误差大大增加，具体表现为总能量随 $\boldsymbol{k}$ 点数目收敛的速度较慢。
再考虑对上述误差的修正。设一个普遍的二次函数 $X(\boldsymbol{k})$ 在以四面体三条边为坐标轴的坐标系内可以表示为
$$ X(\boldsymbol{k})=\bar{a}+\sum_{i}\bar{b}_ik_i+\frac{1}{2}\sum_{i,j}k_i\bar{c}_{ij}k_j,\quad i,j=1,2,3 $$
相应的线性差值函数 $x(\boldsymbol{k})$ 则为
$$ x(\boldsymbol{k})=\bar{a}+\sum_{i}\left(\bar{b}_i+\frac{1}{2}\bar{c}_{ii}\right)k_i $$
则四面体中的误差为
$$ \delta\langle X\rangle_T=\int_{V_T}\mathrm{d}^{3}k\cdot\frac{1}{2}\left(\sum_{ij}k_i\bar{c}_{ij}k_j-\sum_{i}\bar{c}_{ii}k_i\right)=V_T\cdot\frac{1}{40}\left(\sum_{i\neq j}\bar{c}_{ij}-3\sum_{i}\bar{c}_{ii}\right) \tag{3.375} $$
其中最后一步利用了公式\cite{abramowitz1972handbook}
$$ \frac{1}{V_T}\int_{V_T}X(\boldsymbol{k})\,\mathrm{d}^{3}\boldsymbol{k}=\frac{1}{40}\sum F_{\mathrm{v}}+\frac{9}{40}\sum F_{\mathrm{f}}+O(4) $$
式中：$F_{\mathrm{v}}$ 和 $F_{\mathrm{f}}$ 分别指四面体顶点处及面心处的函数值。利用方程（3.375）可以计算整个第一布里渊区（必须是全部区域，仅包含不可约部分的求和无效）中的误差总值 $\delta\langle X\rangle$。具体的计算过程可以在文献\cite{blochl1994improved}中找到，这里不赘述，仅给出结果。利用高斯定理以及有限差分近似，最终可得
$$ \delta\langle X\rangle=\sum_{T}D_T(E_{\mathrm{F}})\cdot\frac{1}{40}\sum_{i=1}^{4}X_i\sum_{j=1}^{4}(E_j-E_i) \tag{3.376} $$
则相应的修正权重 $\mathrm{d}\omega_i$ 为
$$ \mathrm{d}\omega_i=\frac{\mathrm{d}\delta\langle X\rangle}{\mathrm{d}X_i}=\sum_{T}\frac{1}{40}D_T(E_{\mathrm{F}})\sum_{j=1}^{4}(E_j-E_i) \tag{3.377} $$
在利用方程（3.367）计算形如式（3.359）的积分时，对于金属体系，计入上述修正项会有效地改善被积函数的收敛性。