首页

四面体积分法

星期三, 2009-05-06 02:31 — satchel1979

四面体积分法的基本公式
 

 第一原理计算方法中，基于平面波展开的周期性方法因其高性能得到了越来越广泛的应用。这一类方法的实现中有一个关键性的步骤，即算符的期待值或平均值的计算需要在体系的第一布里渊区（或其不可约部分）进行积分。比较常用的数值积分方法有两大类：特殊点法(special points)以及解析四面体积分法 (analytical tetrahedron method)。其中后者（以下简称为四面体法）是一种更为普适的方法，除了可以应用于计算绝缘体和金属体系的算符期待值之外，还可以给出体系的谱函数，也即可以计算体系的动态响应性质。

 体积分

 在这一节里首先给出四面体法的基本概念和一些重要的公式。对于算符$ F $，其期待值$ \langle F\rangle $在倒空间中计算如下：

$ \langle F\rangle = \frac{1}{V_G}\sum_{n}\int_{V_G}d^3k F_n(\vec{k})f[E_n(\vec{k})]. (1) $

 上式中$ V_G $是第一布里渊区的体积，$ f(\epsilon) $是费米-狄拉克分布函数，而$ F_n(\vec{k}) $是算符$ F $在第$ n $条能带中指定$ \vec{k} $点上的值。 为计算这个积分，四面体方法将第一布里渊区分解成若干个小的四面体，然后在各个四面体中分别计算积分，再对所有四面体求和，即

$ \frac{1}{V_G}\int_{V_G}F(\vec{k})d^3 k = \frac{1}{V_G}\sum^{N}_{j=1}\int_{V_T}F(\vec{k})d^3 k. (2) $
 
 而在每个四面体内$ F_n(\vec{k}) $采用线性函数$ f(\vec{k})=a_0+a_1k_x+a_2k_y+a_3k_z $近似。$ f $的线性系数由边界条件$ f(\vec{k_i})=F_{n,i} (i=1,2,3,4) $决定，$ i $是四面体的四个顶点。这样，$ F_n(\vec{k}) $在每个四面体中的积分为
 
 $ \frac{1}{V_G}\int_{V_T}F(\vec{k})d^3 k = \frac{V_T}{V_G}\sum^4_{j=1}\frac{F_{n,j}}{4}. (3) $
 
 其中$ F^j_n $为四面体顶点的函数值。可以证明$ \frac{1}{4} $这个权重来源于积分[1]

$ \omega_j = \int_{T_0}k_xd^3k, j = 1,2,3 (4) $
 
 其中$ T_0 $是顶点分别为(0,0,0),(1,0,0),(0,1,0),(0,0,1)的参考四面体。而$ \omega_4 = 1-\sum_{j=1}^3\omega_j $。更严格的说，考虑$ f(E) $，上式中的前置系数应改为$ V_T^{occ} $，也就是这个四面体小于费米能级$ E_F $的部分的体积。上式也是现有算法中最基本的表达式。在不产生误解的前提下，此后的讨论中如无必要，能带指数$ n $将不再出现。

 更严格的工作指出[1,2]，利用线性函数近似$ F(\vec{k}) $有些情况下并不是最理想的选择，这时对其更好的近似应该取为两个线性函数的商
  $ F(\vec{k})=\frac{f(\vec{k})}{g(\vec{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(\vec{k}) $的系数由边界条件$ g(\vec{k_i})=E_i (i=1,2,3,4) $决定，其中$ E_i $是第$ i $个顶点上的能量本征值。相应的$ f(\vec{k}) $在顶点处等于$ F_i\cdot E_i $。经推导可得，这时四面体的积分可表示为

$ \int_{V_T}F(\vec{k})d^3 k = \frac{V_T}{V_G}\sum^4_{i=1}F_i\omega_i (5) $
 
 其中$ \omega_i $为
 
 $ \omega_i=\frac{1}{\prod_{k\neq i}(1-\frac{E_k}{E_i})}+\sum_{k\neq i}\frac{1}{\prod_{l\neq k}(1-\frac{E_l}{E_k})}\frac{\ln\frac{E_k}{E_i}}{\frac{E_k}{E_i}-1}. (6) $
 
 考虑到简并情况，Zaharioudakis给出了比较完整的权重表示[2]：
1. $ E_1 $ < $ E_2 $ < $ E_3 $ < $ E_4 $ 同方程(6)
2. $ E_1 $ = $ E_2 $ < $ E_3 $ < $ E_4 $
      
       $ \omega_1=\omega_2=\frac{5}{2(1-\frac{E_3}{E_1})(1-\frac{E_4}{E_1})}-\frac{1}{(1-\frac{E_3}{E_1})^2(1-\frac{E_4}{E_1})}-\frac{1}{(1-\frac{E_3}{E_1})(1-\frac{E_4}{E_1})^2}+\frac{1}{(1-\frac{E_1}{E_3})^3(1-\frac{E_4}{E_3})}\frac{\ln\frac{E_3}{E_1}}{\frac{E_3}{E_1}}+\frac{1}{(1-\frac{E_1}{E_4})^3(1-\frac{E_3}{E_4})}\frac{\ln\frac{E_4}{E_1}}{\frac{E_4}{E_1}}. (7) $
      
       $ \omega_3=\frac{1}{(1-\frac{E_1}{E_3})^2(1-\frac{E_4}{E_3})}+\frac{\frac{E_3}{E_1}}{(1-\frac{E_3}{E_1})^2(1-\frac{E_4}{E_1})}+\frac{1}{(1-\frac{E_1}{E_4})^2(1-\frac{E_3}{E_4})^2}\frac{\ln\frac{E_4}{E_3}}{\frac{E_4}{E_3}}+[\frac{3}{(1-\frac{E_3}{E_1})^2(1-\frac{E_4}{E_1})}-\frac{1}{(1-\frac{E_3}{E_1})^2(1-\frac{E_4}{E_1})^2}-\frac{2}{(1-\frac{E_3}{E_1})^3(1-\frac{E_4}{E_1})}]\frac{\ln\frac{E_1}{E_3}}{\frac{E_1}{E_3}}. (8) $
       
       $ \omega_4 $与$ \omega_3 $形式相同，只需要将上式中的$ E_3 $和$ E_4 $互换即可。
3. $ E_1 $ <  $ E_2 $ = $ E_3 $ < $ E_4 $
      
       $ \omega_2=\omega_3=\frac{5}{2(1-\frac{E_1}{E_2})(1-\frac{E_4}{E_2})}-\frac{1}{(1-\frac{E_1}{E_2})^2(1-\frac{E_4}{E_2})}-\frac{1}{(1-\frac{E_1}{E_2})(1-\frac{E_4}{E_2})^2}+\frac{1}{(1-\frac{E_2}{E_1})^3(1-\frac{E_4}{E_1})}\frac{\ln\frac{E_1}{E_2}}{\frac{E_1}{E_2}}+\frac{1}{(1-\frac{E_1}{E_4})(1-\frac{E_2}{E_4})^3}\frac{\ln\frac{E_4}{E_2}}{\frac{E_4}{E_2}}. (9) $

       $ \omega_1=\frac{1}{(1-\frac{E_2}{E_1})^2(1-\frac{E_4}{E_1})}+\frac{\frac{E_1}{E_2}}{(1-\frac{E_1}{E_2})^2(1-\frac{E_4}{E_2})}+\frac{1}{(1-\frac{E_1}{E_4})^2(1-\frac{E_2}{E_4})^2}\frac{\ln\frac{E_4}{E_1}}{\frac{E_4}{E_1}}+[\frac{3}{(1-\frac{E_1}{E_2})^2(1-\frac{E_4}{E_2})}-\frac{1}{(1-\frac{E_1}{E_2})^2(1-\frac{E_4}{E_2})^2}-\frac{2}{(1-\frac{E_1}{E_2})^3(1-\frac{E_4}{E_2})}]\frac{\ln\frac{E_2}{E_1}}{\frac{E_2}{E_1}}. (10) $
  
       $ \omega_4 $与$ \omega_1 $形式相同，只需要将上式中的$ E_1 $和$ E_4 $互换即可。
4. $ E_1 $ < $ E_2 $ < $ E_3 $ = $ E_4 $
   
       $ \omega_3=\omega_4=\frac{5}{2(1-\frac{E_1}{E_3})(1-\frac{E_2}{E_3})}-\frac{1}{(1-\frac{E_1}{E_3})^2(1-\frac{E_2}{E_3})}-\frac{1}{(1-\frac{E_1}{E_3})(1-\frac{E_2}{E_3})^2}+\frac{1}{(1-\frac{E_2}{E_1})(1-\frac{E_3}{E_1})^3}\frac{\ln\frac{E_1}{E_3}}{\frac{E_1}{E_3}}+\frac{1}{(1-\frac{E_1}{E_2})(1-\frac{E_3}{E_2})^3}\frac{\ln\frac{E_2}{E_3}}{\frac{E_2}{E_3}}. (11) $

       $ \omega_1=\frac{1}{(1-\frac{E_3}{E_1})^2(1-\frac{E_2}{E_1})}+\frac{\frac{E_1}{E_3}}{(1-\frac{E_1}{E_3})^2(1-\frac{E_2}{E_3})}+\frac{1}{(1-\frac{E_1}{E_2})^2(1-\frac{E_3}{E_2})^2}\frac{\ln\frac{E_2}{E_1}}{\frac{E_2}{E_1}}+[\frac{3}{(1-\frac{E_1}{E_3})^2(1-\frac{E_2}{E_3})}-\frac{1}{(1-\frac{E_1}{E_3})^2(1-\frac{E_2}{E_3})^2}-\frac{2}{(1-\frac{E_1}{E_3})^3(1-\frac{E_2}{E_3})}]\frac{\ln\frac{E_3}{E_1}}{\frac{E_3}{E_1}}. (12) $

       $ \omega_2 $与$ \omega_1 $形式相同，只需要将上式中的$ E_1 $和$ E_2 $互换即可。
5. $ E_1 $ = $ E_2 $ = $ E_3 $ < $ E_4 $
      
      $ \omega_1=\omega_2=\omega_3=\frac{11}{6(1-\frac{E_4}{E_1})}-\frac{5}{2(1-\frac{E_4}{E_1})^2}+\frac{1}{(1-\frac{E_4}{E_1})^3}+\frac{1}{(1-\frac{E_1}{E_4})^4}\frac{\ln\frac{E_4}{E_1}}{\frac{E_4}{E_1}}. (13) $

      $ \omega_4=[\frac{3}{(1-\frac{E_4}{E_1})^2}-\frac{6}{(1-\frac{E_4}{E_1})^3}+\frac{3}{(1-\frac{E_4}{E_1})^4}]\frac{\ln\frac{E_1}{E_4}}{\frac{E_1}{E_4}}+\frac{1}{(1-\frac{E_1}{E_4})^3}+\frac{5}{2}\frac{\frac{E_4}{E_1}}{(1-\frac{E_4}{E_1})^2}-2\frac{\frac{E_4}{E_1}}{(1-\frac{E_4}{E_1})^3}. (14) $
     
6. $ E_1 $ < $ E_2 $ = $ E_3 $ = $ E_4 $
      
      $ \omega_1=[\frac{3}{(1-\frac{E_1}{E_2})^2}-\frac{6}{(1-\frac{E_1}{E_2})^3}+\frac{3}{(1-\frac{E_1}{E_2})^4}]\frac{\ln\frac{E_2}{E_1}}{\frac{E_2}{E_1}}+\frac{1}{(1-\frac{E_2}{E_1})^3}+\frac{5}{2}\frac{\frac{E_1}{E_2}}{(1-\frac{E_1}{E_2})^2}-2\frac{\frac{E_1}{E_2}}{(1-\frac{E_1}{E_2})^3}. (15) $
      
      $ \omega_2=\omega_3=\omega_4=\frac{11}{6(1-\frac{E_1}{E_2})}-\frac{5}{2(1-\frac{E_1}{E_2})^2}+\frac{1}{(1-\frac{E_1}{E_2})^3}+\frac{1}{(1-\frac{E_2}{E_1})^4}\frac{\ln\frac{E_1}{E_2}}{\frac{E_1}{E_2}}. (16) $
     
7. $ E_1 $ = $ E_2 $ < $ E_3 $ = $ E_4 $
      
      $ \omega_1=\omega_2=\frac{5}{2(1-\frac{E_3}{E_1})^2}-\frac{2}{(1-\frac{E_3}{E_1})^3}+\frac{\frac{E_1}{E_3}}{(1-\frac{E_1}{E_3})^3}+[\frac{3}{(1-\frac{E_1}{E_3})^3}-\frac{3}{(1-\frac{E_1}{E_3})^4}]\frac{\ln\frac{E_3}{E_1}}{\frac{E_3}{E_1}}. (17) $
     
      $ \omega_3=\omega_4=\frac{5}{2(1-\frac{E_1}{E_3})^2}-\frac{2}{(1-\frac{E_1}{E_3})^3}+\frac{\frac{E_3}{E_1}}{(1-\frac{E_3}{E_1})^3}+[\frac{3}{(1-\frac{E_3}{E_1})^3}-\frac{3}{(1-\frac{E_3}{E_1})^4}]\frac{\ln\frac{E_1}{E_3}}{\frac{E_1}{E_3}}. (18) $
     
8. $ E_1 $ = $ E_2 $ = $ E_3 $ = $ E_4 $
      
      $ \omega_1=\omega_2=\omega_3=\omega_4=\frac{1}{4}. (19) $
 

 态密度

 四面体法的提出，最早是为了解决形如

 $ I(E) = \int_{E(\vec{k})=E}F(\vec{k})|\nabla E(\vec{k})|^{-1}dS. (20) $
 
 类的积分。当$ F(\vec{k})\equiv1 $，积分式前乘以因子$ 1/V_G $，上式就是态密度的定义。除了在物理上的重要性以外，讨论态密度有助于直观的理解四面体法的几何意义。

 Lehmann和Taut指出[3]，设四个顶点的能量本征值满足条件$ E_1 $ < $ E_2 $ < $ E_3 $ < $ E_4 $，在四面体内能量按照线性函数展开，与上节的$ g(\vec{k}) $相同：

 $ E(\vec{k})=E_1+\vec{b}\cdot(\vec{k}-\vec{k_1}). (21) $

 其中$ \vec{b}=\sum^3_{i=1}[E_{i+1}-E_1]\vec{r_i} $以及$ \vec{r_i}\cdot\vec{k_j}=\delta_{ij} $。上式中$ \vec{k_j}=\vec{k_{j+1}}-\vec{k_1} $，因此可以按照倒格矢与正格矢的关系式写出$ \vec{r_i} $如下：

$ \vec{r_1}=\frac{\vec{k_2}\times\vec{k_3}}{\vec{k_1}\cdot(\vec{k_2}\times\vec{k_3})};\vec{r_2}=\frac{\vec{k_3}\times\vec{k_1}}{\vec{k_1}\cdot(\vec{k_2}\times\vec{k_3})}; \vec{r_3}=\frac{\vec{k_1}\times\vec{k_2}}{\vec{k_1}\cdot(\vec{k_2}\times\vec{k_3})}. (22) $

 分母是以$ \vec{k_i} $为顶点的四面体体积的六倍$ 6V_T $。因此该四面体对于态密度$ D_T(E) $的贡献为

  $ D_T(E)=1/V_G\frac{dS(E)}{|\vec{b}|}=0, E \leq E_1 $
                   $ =1/V_G\frac{f_1}{|\vec{b}|}, E_1 \leq E \leq E_2 $
                   $ =1/V_G\frac{f_1-f_2}{|\vec{b}|}, E_2 \leq E \leq E_3 $
                   $ =1/V_G\frac{f_4}{|\vec{b}|}, E_3 \leq E \leq E_4 $
                   $ =0, E $ > $ E_4 (23) $

 其中函数$ f $的几何意义是等能面$ S(E) $在四面体内的截面面积。因此容易得出$ f/|\vec{b}| $的表达式为
 
  $ f_1/|\vec{b}|=3V_T\frac{(E-E_1)^2}{(E_2-E_1)(E_3-E_1)(E_4-E_1)}, $
  $ f_2/|\vec{b}|=3V_T\frac{(E-E_2)^2}{(E_2-E_1)(E_3-E_2)(E_4-E_2)}, $
  $ f_4/|\vec{b}|=3V_T\frac{(E-E_4)^2}{(E_4-E_1)(E_4-E_2)(E_4-E_3)}. (24) $

 实际上对于态密度，Jepsen和Andersen[4]还提出了更直观简单的计算方法，Blöchl[5]也采用了这种方法，并做了改进。可以将每个四面体看作是容器，给定等能面$ S(E) $之后，该四面体对于电子态数目$ n(E) $的贡献等于$ \epsilon $ $ \leq $ $ E $填充的体积相对于第一布里渊区体积的比值。而态密度$ D_T(E) $可以定义为$ dn/dE $。通过简单的三角锥体积计算，可以得到

  $ n(E)=0, E\leq E_1; $
         $ =\frac{V_T}{V_G}\frac{(E-E_1)^3}{(E_2-E_1)(E_3-E_1)(E_4-E_1)}, E_1 \leq E \leq E_2; $
         $ =\frac{V_T}{V_G}[\frac{(E-E_1)^3}{(E_2-E_1)(E_3-E_1)(E_4-E_1)}- \frac{(E-E_1)^3}{(E_2-E_1)(E_3-E_2)(E_4-E_2)}], E_2 \leq E \leq E_3; $
         $ =\frac{V_T}{V_G}[(1-\frac{(E_4-E)^3}{(E_4-E_1)(E_4-E_2)(E_4-E_3)}], E_3 \leq E \leq E_4; $
         $ =\frac{V_T}{V_G}, E $ > $ E_4 $                          (25)

 上式对于$ E $求导，即得方程(23)。对于简并情况，方程(23)-(25)并不适用。因此对于上式的第三种情况，可等效的写为[4]
  当$ E_2 \leq E \leq E_3 $时，
  $ n(E)=\frac{V_T}{V_G}\frac{1}{(E_3-E_1)(E_4-E_1)}[(E_2-E1)^2+3(E_2-E_1)(E-E_2)+3(E-E_2)^2-\frac{E_3-E_1+E_4-E_2}{(E_3-E_2)(E_4-E_2)}(E-E_2)^2] $
 
 而相应的，该情况下的态密度$ D_T(E) $也可写为

  $ D_T(E)=\frac{V_T}{V_G}\frac{1}{(E_3-E_1)(E_4-E_1)}[3(E_2-E_1)+6(E-E_2)-3\frac{(E_3-E_1+E_4-E_2)(E-E_2)^2}{(E_3-E_2)(E_4-E_2)}], E_2 \leq E \leq E_3 $

图1：等能面截面图图1：等能面截面图 

 

 其他函数的面积分

 考虑方程(20)，对于任意函数$ F(\vec{k}) $的加权积分可以利用四面体法通过类似于态密度的讨论进行计算[3]。将$ F(\vec{k}) $在四面体内也利用线性函数近似，则类似于方程(21)，$ F(\vec{k}) $可写为

 $ F(\vec{k}) = f_0+\vec{f}\cdot(\vec{k}-\vec{k_1}) (26) $

 其中$ \vec{f}=\sum^3_{i=1}[F_{i+1}-F_1]\vec{r_i} $。因此每个四面体对于方程(20)的贡献为

  $ i(E)=1/V_G\frac{dS(E)\cdot\vec{k}}{|\vec{b}|}=0, E \leq E_1 $
        $ =1/V_G\frac{f_0s_0}{|\vec{b}|}, E_1 \leq E \leq E_2 $
        $ =1/V_G\frac{f_0-f_1}{|\vec{b}|}\frac{s_1f_1-s_2f_2}{f_1-f_2}, E_2 \leq E \leq E_3 $
        $ =1/V_G\frac{f_4s_4}{|\vec{b}|},E_3 \leq E \leq E_4 $
        $ =0, E $ > $ E_4 (27) $

 其中$ s_i $是密度随$ k $呈线性变化的截面$ f_i $的重心：
 
 $ s_i(E)=k_i+\frac{E-E_i}{3}\sum_{j=1,j\neq i}^{4}\frac{k_j-k_i}{E_j-E_i}. (28) $
 
 从上面几个公式也可以看出四面体法的一个优点，就是在计算公式中不出现被积函数的导数值。这一点大大的简化了计算，使得编写程序更为容易。

从格林函数的角度统一思考
 

 上一部分中我们分别给出了体积分和面积分的四面体法计算公式。这两类积分所以重要，是因为固体物理中绝大多数的响应性质的计算公式都可以归结到其中（例如有效质量张量、电导率以及磁化率等等）。重新写出更为普遍的体积分式如下

 $ J(E)={\rm P}\int_{V_G}\frac{F(\vec{k})}{E-E(\vec{k})}d^3k. (29) $

其中$ {\rm P} $代表积分主值。这样，方程(20)和(29)分别对应于体系的格林函数$ G(z) $的虚部和实部。因此，利用格林函数方法可以统一得到这两者的计算公式。现简要介绍如下：

 在给定基函数的表象下，预解算符可表示为[6]

 $ (z-H)^{-1}=\sum_n\sum_{\vec{k}}|\psi_n(\vec{k})\rangle\langle\psi_n(\vec{k})|/[z-E_n(\vec{k})] (30) $

 将上式中对$ k $点求和改为对第一布里渊区积分，并且将谱投影算符$ |\psi_n(\vec{k})\rangle\langle\psi_n(\vec{k})| $的矩阵元记为$ F_n(\vec{k}) $，则可以得出预解算符的矩阵元为

 $ R(z)=1/V_G\sum_n\int_{V_G}F_n(\vec{k})/[z-E_n(\vec{k})]d^3k (31) $

 与第一部分相同，利用$ k $点网格将第一布里渊区划分为若干四面体，则对于给定的能带指数$ n $，每个四面体对于预解矩阵元$ R(z) $的贡献为

 $ \sum_{i=1}^{4}r_i^n(z)F_n(\vec{k_i})V_T/V_G (32) $
 

 其中复权重$ r_i^n(z) $遵循方程

 $ r_i^n(z)=\frac{(z-E_i^n)^2}{\Pi_{k(\neq i)}(E^n_k-E^n_i)}+\sum_{j(\neq i)}[\frac{(z-E^n_j)^3}{\Pi_{k(\neq i)}(E^n_k-E^n_i)}\frac{\ln[(z-E^n_j)/(z-E^n_i)]}{E^n_i-E^n_j}] (33) $

 实宗量的格林函数$ G(E) $可以看作是预解矩阵元$ R(z) $沿虚轴正半轴接近实轴的极限，即

 $ G(E)=R(E+i0^+) (34) $

 相应的复权重可以写为

 $ r_i(E+i0^+)=d_i(E)-i\pi c_i(E) (35) $

 利用恒等式
 $ \sum_{j(\neq i)}[\frac{(z-E_j)^3}{\Pi_{k(\neq j)(E_k-E_j)}}\frac{1}{E_i-E_j}]=-\frac{(z-E_i)^2}{\Pi_{k(\neq i)}(E_k-E_i)}\sum_{j(\neq i)}\frac{z-E_j}{E_i-E_j} (36) $

 代入方程(33)，就可以得出$ d_i(E) $和$ c_i(E) $分别为

  $ d_i(E)=\frac{(E-E_i)^2}{\Pi_{k(\neq i)}(E_k-E_i)}[1+\sum_{j(\neq i)}\frac{E-E_j}{E_i-E_k}\ln|E-E_i|]+\sum_{j(\neq i)}[\frac{(E-E_j)^3}{\Pi_{k(\neq j)}(E_k-E_j)}\frac{\ln|E-E_j|}{E_i-E_j}]; $

  $ c_i(E)=\frac{(E-E_i)^2}{\Pi_{k(\neq i)}(E_k-E_i)}\sum_{j(\neq i)}\frac{E-E_j}{E_i-E_j}\theta(E-E_i)+\sum_{j(\neq i)}[\frac{(E-E_j)^3}{\Pi_{k(\neq j)}(E_k-E_j)}\frac{\theta(E-E_j)}{E_i-E_j}] (37) $

 其中$ \theta(x) $是Heaviside单位阶跃函数。因此结合方程(32)、(34)、(35)可得给定能带中，每个四面体对于$ G(E) $的贡献为
 
 $ G(E)=\sum^4_{i=1}[d_i(E)-i\pi c_i(E)]F_iV_T/V_G (38) $
 
 格林函数理论表明[7]
$ I(E)=-\frac{1}{\pi} {\rm Im}G(E) $
$ J(E)={\rm Re}G(E) (39) $

方程将(37)中的$ d_i(E) $宗量$ E $取为0，则可得到方程(5)，方程(20)中的$ F_{\vec{k}}\equiv1 $，则可得方程(23)。一般的面积分可表为$ \frac{V_T}{V_G}\sum_{i=1}^4c_i(E)F_i $。

 通用的四面体划分法及改进的四面体法

 传统的四面体法中，首先通过对称群的操作找出第一布里渊区的不可约部分。然后在等间距的$ k $点网格上将其手动划分为若干四面体。这种划分有一定的随意性，因为给定一组$ k $点网格， 可以有很多种不同的方法将其围成若干互不重叠的四面体，进而进行计算。这种随意性是否会造成计算上的误差甚至错误？ Kleinman对此做了讨论[8]，利用一个简单的例子，他证明在$ k $点比较稀疏的情况下，不同的划分会造成高达14%的误差，而正确的划分会给出与特殊点法相同的结果。而在$ k $点足够稠密时，这种划分上的随意性对计算结果的影响可以忽略不计。其原因在于四面体法中，除个别$ k $点以外，绝大多数$ k $点由多个四面体共享，因此在积分中实际上参与超过一次计算，边界及靠近边界的$ k $点在不同的划分中所归属的四面体个数不尽相同，从而使得各自对于积分值的贡献有所差别，而网格中间区域的$ k $点不存在这个问题。因此$ k $点稠密时，边界处$ k $点贡献可以忽略，但是$ k $点较少时这种差别则会影响最终结果。Kleinman的工作使得Jepsen和Andersen对四面体法重新思考，从而提出了通用的、适用于编程的四面体划分法[9]，为后来Blöchl的工作奠定了基础。

 Kleinman的工作的另一个重要性在于第一次明确指出了利用四面体法，同样可以将方程(1)所示积分转化为各个$ k $点上的被积函数值的加权求和：

 $ \langle F\rangle = \frac{1}{4V_G}\sum_n\sum_iF_n(\vec{k_i})(\sum_{T\ni \vec{k_i}}V_T^{occ}) (40) $

 其中圆括号中的求和遍历所有以$ \vec{k_i} $为一个顶点的四面体，而$ V_T^{occ} $则是这个四面体中满足$ E_n(\vec{k}) $ $ \leq $ $ E_F $的体积，具体公式参见方程(25)。实际上，这一点也可以从方程(3)直接推出。

 在前面工作的基础上，Blöchl给出了包含较大改进的四面体法的普适算法[5]，有三个主要特点：四面体的自动划分、各$ k $点的权重以及对于金属体系的Blöchl修正。下面分别介绍如下。

 四面体的自动划分

 前面已经指出，为了减少需要计算的四面体数目，首先需要找出第一布里渊区的不可约部分。这样的策略有一个副作用，即对不可约部分的四面体划分几乎不可避免的要进行人工干预，不利于编程求解。因此Blöchl提出，利用Monkhorst-Pack方法(MP)[10]首先在第一布里渊区内生成等距的$ k $点网格。然后给每个$ k $点编号：

 $ N=1+\frac{l-l_0}{2}+(n_1+1)(\frac{m-m_0}{2}+(n_2+1)\frac{n-n_0}{2}) (41) $

 其中$ (l,m,n) $分别是该$ k $点沿倒格矢$ b_1,b_2,b_3 $的序数的二倍，$ n_i $分别是在三个方向上的$ k $点数。$ (l_0,m_0,n_0) $是MP方法中$ \Gamma $点的偏移量，有偏移则为1，否则为0。编号之后建立标识数组，其位置与该位置储存的元素值相同，例如，第一个位置存储“1”，第二个位置存储“2”，依此类推。然后从第一个位置开始，利用对称群的操作矩阵对每个$ k $点坐标进行操作，再与数组中其他$ k $点的坐标进行比较，如果彼此相同且后者的编号大于前者，即将后者的元素值改为前者。这样对全部数组操作完毕之后，立即可以挑出所有不可约$ k $点——只有当$ k_i $为不可约$ k $点时，其编号才与其存储位置相同。之后，为了计算方便，可以对所有这些不可约$ k $点按存储位置的顺序重新编号，即从“1”到“$ k_{irr}({\rm max}) $”。数组中的各个元素也相应的改为新的编号（或叫名称）。这样整个第一布里渊区中的$ k $点都可用不可约$ k $点标记。

 下一步讨论四面体的自动划分过程。以下八组坐标代表的$ k $点构成平行六面体：
 
  $ (l,m,n)-1,(l+2,m,n)-2,(l,m+2,n)-3,(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 $
 
 为了尽量减小插值引起的误差，可以取此平行六面体中最短的体对角线作为等体积的六个四面体的公共对角线，设为$ 1-8 $，则可以采用下面6组途径确定这六个四面体的各个顶点：

  $ 1\rightarrow2\rightarrow5\rightarrow8;\quad 1\rightarrow2\rightarrow6\rightarrow8;\quad 1\rightarrow3\rightarrow5\rightarrow8; $
  $ 1\rightarrow3\rightarrow7\rightarrow8;\quad 1\rightarrow4\rightarrow6\rightarrow8;\quad 1\rightarrow4\rightarrow7\rightarrow8. $

 对每个平行六面体重复上述过程，可以将整个第一布里渊区划分为体积相等的若干四面体，每个四面体的顶点可用标识数组中的不可约$ k $点标记。将这四个顶点的标号按升序排列，则每个四面体的简并度可以轻易得出。因此，这个过程保证了可以只用不可约$ k $点上的信息进行整个第一布里渊区的积分，而无须考虑如何划定其不可约部分。

 上述过程可以避免Kleinman所说的计算误差，而且整个过程可以通过程序自动实现而无须人工干预。

 

 各个$ k $点的权重计算

 四面体法的基本过程是算出每个四面体对积分值的贡献再对所有四面体求和。因此在本节开始时已经指出，绝大多数$ k $参与了多次计算。同时这种算法也需要同时知道每个不可约$ k $点上的被积函数值。因此对于大规模计算而言难免会有存储方面的困难。如同Kleinman指出，方程(1)可以被表示为$ k $点的加权求和

 $ \langle F\rangle = \sum_{i,n}F_n(\vec{k_i})\omega_{ni} (42) $

 与标准的四面体法相同，在每个四面体内，$ F_n(\vec{k}) $用线性函数$ f_n(\vec{k}) $近似。等价于方程(26)，$ f_n(\vec{k}) $也可写为

 $ f_n(\vec{k})=\sum_iF_n(\vec{k_i})\omega_i(\vec{k_i}) (43) $

 将其带入方程(41)，可得每个不可约$ k $点的积分权重$ \omega_{ni} $：

 $ \omega_{ni}=\frac{1}{V_G}\int_{V_G}\omega_i(\vec{k})f[E_n(\vec{k})] (44) $

 若取$ \omega_i(\vec{k}) $为常数1/4，则根据方程(44)很容易得到方程(40)。事实上方程(40)非常简单实用。但是考虑到问题的完整性以及与修正项的自洽问题，我们在这里还是重复一下 Blöch的计算。将$ \omega_i(\vec{k}) $作为线性函数，满足在$ k_i $以及等价点$ \{k^{*}_i\} $上为1，而在其他$ k $点均为0（当然在两$ k $点间$ \omega_i(\vec{k}) $成线性变化）。为计算这种情况下的积分权重$ \omega_{ni} $，首先必须计算出体系的费米能级$ E_F $：利用方程(25)计算能量小于给定$ E $的状态数$ n(E) $，直到$ \sum_ng_nn(E_n) $等于体系的总电子数（$ g_n $为该能级的简并度），此时的$ E $即为所求$ E_F $。据此可以得出$ \omega_{ni} $（省略能带指数$ n $）：

当$ E_F $ < $ E_1 $时，
$ \omega_1=\omega_2=\omega_3=\omega_4=0. (45) $

当E_1 < E_F < E_2时，
  $ \omega_1=C_0[4-(E_F-E_1)(\frac{1}{E_2-E_1}+\frac{1}{E_3-E_1}+\frac{1}{E_4-E_1})], $
  $ \omega_2=C_0\frac{E_F-E_1}{E_2-E_1}, $
  $ \omega_3=C_0\frac{E_F-E_1}{E_3-E_1}, $
  $ \omega_4=C_0\frac{E_F-E_1}{E_4-E_1}. (46) $

 其中$ C_0 $为
 
  $ C_0=\frac{V_T}{4V_G}\frac{(E_F-E_1)^3}{(E_2-E_1)(E_3-E_1)(E_4-E_1)}. $
 

当$ E_2 $ < $ E_F $ < $ E_3 $时，
  $ \omega_1=C_1+(C_1+C_2)\frac{E_3-E_F}{E_3-E_1}+(C_1+C_2+C_3)\frac{E_4-E_F}{E_4-E_1}, $
  $ \omega_2=C_1+C_2+C_3+(C_2+C_3)\frac{E_3-E_F}{E_3-E_2}+C_3\frac{E_4-E_F}{E_4-E_2}, $
  $ \omega_3=(C_1+C_2)\frac{E_F-E_1}{E_3-E_1}+(C_2+C_3)\frac{E_F-E_2}{E_3-E_2}, $
  $ \omega_4=(C_1+C_2+C_3)\frac{E_F-E_1}{E_4-E_1}+C_3\frac{E_F-E_2}{E_4-E_2}. (47) $
 
 其中
  $ C_1=\frac{V_T}{4V_G}\frac{(E_F-E_1)^2}{(E_4-E_1)(E_3-E_1)}, $
  $ C_2=\frac{V_T}{4V_G}\frac{(E_F-E_1)(E_F-E_2)(E_3-E_F)}{(E_4-E_1)(E_3-E_2)(E_3-E_1)}, $
  $ C_3=\frac{V_T}{4V_G}\frac{(E_F-E_2)^2(E_4-E_F)}{(E_4-E_2)(E_3-E_2)(E_4-E_1)}. $
 

 当$ E_3 $ < $ E_F $ < $ E_4 $时，
  $ \omega_1=\frac{V_T}{4V_G}-C_4\frac{E_4-E_F}{E_4-E_1}, $
  $ \omega_2=\frac{V_T}{4V_G}-C_4\frac{E_4-E_F}{E_4-E_2}, $
  $ \omega_3=\frac{V_T}{4V_G}-C_4\frac{E_4-E_F}{E_4-E_3}, $
  $ \omega_4=\frac{V_T}{4V_G}-C_4[4-(\frac{1}{E_4-E_1}+\frac{1}{E_4-E_2}+\frac{1}{E_4-E_3})(E_4-E_F)]. (48) $

 其中$ C_4 $为
 
  $ C_4=\frac{V_T}{4V_G}\frac{(E_4-E_F)^3}{(E_4-E_1)(E_4-E_2)(E_4-E_3)}. $
 

 当$ E_F $ > $ E_4 $时，
  $ \omega_1=\omega_2=\omega_3=\omega_4=\frac{V_T}{4V_G}. (49) $
 
 因此，对于给定的$ k_i $，可以找出其自身和等价点所属的四面体，通过分析各四面体内的能级分布情况由方程(45)-(49)计算出积分权重$ \omega_{ni} $。

 Blöchl修正

 在四面体中对于被积函数的线性插值近似可以获得简单的计算公式，同时不会过多地牺牲精度。但是对于费米面形状复杂且是部分填充的过渡金属体系而言，还是要考虑线性近似带来的误差及其相关的修正问题。

 首先讨论误差。真实的被积函数$ F(\vec{k}) $在四面体中总有曲率为正以及曲率为负的部分。而线性函数$ f(\vec{k}) $会高估正曲率区间的函数值，低估负曲率区间的函数值。对于绝缘体和半导体而言，因为积分区域覆盖整个四面体，所以这两部分在极大程度上可以抵消。但是对于金属体系，因为导带部分填充，也即积分区域只覆盖四面体的一部分而非全部，所以高估部分和低估部分只有一项占主要地位。这使得误差大大地增加。具体体现为总能随$ k $点数目收敛缓慢。

 再考虑对上述误差的修正。设一个普遍的二次函数$ X(\vec{k}) $，在以四面体三条边为坐标轴的坐标系内可以表示为
 
  $ X(\vec{k})=\tilde{a}+\sum_i\tilde{b}_ik_i+\frac{1}{2}\sum_{i,j}k_i\tilde{c}_{ij}k_j, i,j = 1,2,3. $
 
 相应的线性差值函数$ x(\vec{k}) $则为
 
  $ x(\vec{k})=\tilde{a}+\sum_i(\tilde{b}_i+\frac{1}{2}\tilde{c}_{ii})k_i $
 
 则四面体中的误差为

  $ \delta\langle X\rangle_T=\int_{V_T}d^3k\frac{1}{2}[\sum_{ij}k_i\tilde{c}_{ij}k_j-\sum_i\tilde{c}_{ii}k_i] =V_T\frac{1}{40}[\sum_{i\neq j}\tilde{c}_{ij}-3\sum_i\tilde{c}_{ii}] (50) $

 其中最后一步利用公式（参看文献[11]）

  $ \frac{1}{V_T}\int_{V_T}X(\vec{k})d^3k=\frac{1}{40}\sum F_v+\frac{9}{40}\sum F_f+\mathcal{O}(4). $

 其中$ F_v $和$ F_f $分别指四面体顶点处及面心处的函数值。利用方程(50)可以计算整个第一布里渊区（必须是全部区域，仅包含不可约部分的求和无效）中的误差总值$ \delta\langle X\rangle $。具体的计算过程可以在文献[5]中找到，这里不再详述，仅给出结果。利用高斯定理以及有限差分近似，最终可得

  $ \delta\langle X\rangle=\sum_TD_T(E_F)\frac{1}{40}\sum_{i=1}^4X_i\sum_{j=1}^4(E_j-E_i). (51) $

 则相应的修正权重$ d\omega_i $为

  $ d\omega_i=\frac{d\delta\langle X\rangle}{dX_i}=\sum_T\frac{1}{40}D_T(E_F)\sum_{j=1}^4(E_j-E_i). (52) $

 在利用方程(42)计算形如式(20)和(29)的积分时，对于金属体系，计入上述修正会有效地改善被积函数的收敛行为。

小结

 以上各节较为全面的介绍了现行的四面体积分方法。对于不同类型的积分，其计算表达式不尽相同。这里简要总结一下上文中比较重要的公式。

 1. 利用传统四面体法，对于形如方程(1)及(29)，采用公式(6)-(19)进行计算。
 2. 计算态密度，采用公式(23)、(24)计算；计算积分态密度，采用公式(25)。
 3 对于形如方程(20)，采用公式(27)、(28)进行计算。
 4 利用各$ k $点加权求和形式，对于形如方程(1)及(29)，采用公式(40)、(25)进行计算。
 5. 考虑Blöch修正的各$ k $点加权求和形式，对于形如方程(1)及(29)，采用公式(45)-(49)、(52)进行计算。
 

参考文献

[1] J. Molenaar, P. Coleridge, A. Lodder, J. Phys. C: Solid State Phys. 15, 6955 (1982).
[2] D. Zaharioudakis, Comput. Phys. Comm. 157, 17 (2004).
[3] G. Lehmann, M. Taut, Phys. Status Solidi B 54, 169 (1972).
[4] O. Jepsen, O.K. Andersen, Solid State Commun. 9, 1763 (1971).
[5] P.E. Blöchl. O. Jepsen, O.K. Andersen, Phys. Rev. B 49, 16223 (1994).
[6] P. Lambin, J.P. Vigneron, Phys. Rev. B 29, 3430 (1984).
[7] J. Rath, A.J. Freeman, Phys. Rev. B 11, 2109 (1975).
[8] L. Kleinman, Phys. Rev. B 28, 1139(R) (1983).
[9] O. Jepsen, O.K. Andersen, Phys. Rev. B 29, 5965 (1983).
[10] H.J. Monkhorst, J.D. Pack, Phys. Rev. B 13, 5188 (1976).
[11] M. Abranowitz, I.A. Stegun, Handbook of Mathematical Functions, Dover, p.895 (1972).