布里渊区的特殊k点采样问题研究

介绍

在各种周期性边界条件的第一原理计算方法中,需要涉及到在布里渊区的积分问题,例如总能、电荷密度分布,以及金属体系中费米面的确定等等。如果采用普通的在布里渊区内均匀选取k点的方法,那么为了得到精确的结果$ k $点的密度必须很大,从而导致非常大的计算量。这使得计算的效率非常低下。因此,需要寻找一种高效的积分方法,可以通过较少的$ k $点运算取得较高的精度。而这些k点被称之为“平均值点”(Baldereschi)或者“特殊点”(Chadi, Cohen)。

基本思想

Chadi和Cohen最早提出了这种特殊点的数学基础[1]。考虑一个光滑函数$ g(\vec{k}) $,我们可以将其展为傅立叶级数:

$ g(\vec{k})=g_0+\sum_{m=1} g_m e^{i\vec{k}\cdot\vec{R_m}}.(1) $

假设另有一个拥有体系全部对称性(对称性用对称群$ G $表示)的函数$ f(\vec{k}) $,满足条件$ f(T\vec{k})=f(\vec{k})\forall T\in G $ ,则 我们可以将$ f(\vec{k}) $$ g(\vec{k}) $展开如下:

$ f(\vec{k})=\frac{1}{n_G}\sum_ig(T_i\vec{k})=g_0+\sum_{m=1}^{\infty}\sum_i\frac{1}{n_G}g_me^{iT_i\vec{k}\cdot\vec{R_m}}.(2) $

其中$ n_G $是对称群$ G $的阶数。设$ f_0=g_0 $,将上式的求和顺序重新组合可以得到
$ f(\vec{k})=f_0+\sum_{m=1}\frac{g_m}{n_G}\sum_{T_i\in G}e^{i\vec{k}T_i^{-1}\cdot\vec{R_m}}=f_0+\sum_{m=1}f_m+\sum{|\vec{R}|=C_m}e^{i\vec{k}\cdot\vec{Rl}}=f_0+\sum_{m=1}f_mA_m(\vec{k}).(3) $

其中$ C_m $是距离原点第$ m $近邻的球半径,按升序排列,且$ C_m $ $ \leq $ $ C_{m+1} $。需要注意的是限制条件$ C_m $ $ \leq $ $ |R| $ $ \leq $ $ C_{m+1} $具有球对称性,也即高于$ G $的对称性,所以满足限制条件的格点集合$ \{R\} $并不一定都是等价的——或说可以通过$ G $中的操作$ T $联系起来的——格点。方程(3)中的函数$ A_m $满足下列条件:
$ \Omega/(2\pi)^3\int_{BZ}A_m(\vec{k})d\vec{k}=0,\forall m\{m>0, m\in \mathbb{I}\}; $

$ \Omega/(2\pi)^3\int_{BZ}A_m(\vec{k})A_n(\vec{k})=N_n\delta_{nm}; $

$ A_m(\vec{k}+\vec{G})=A_m(\vec{K}); $

$ A_m(T_i\vec{k})=A_m(\vec{k}); $
$ A_m(\vec{k})A_n(\vec{k})=\sum_ja(j,m,n)A_j(\vec{k}).(4) $
上式中$ \vec{G} $是倒格矢,$ N_n $是满足条件$ |R|=C_n $的格点数。五个方程分别表明函数$ A_m(\vec{k}) $在第一布里渊区内成奇函数、具有正交性、周期性、体系对称性和完备性。对于特殊$ k $点法而言,前两条更为重要。

注意到上面公式中的求和$ m $从1开始,因此需要对$ m=0 $的情况进行单独定义。我们定义$ A_0(\vec{k})=1 $,则函数$ f(\vec{k}) $的平均值为:

$ \overline{f}=\Omega/(2\pi)^3\int_{BZ}f(\vec{k})d\vec{k}=f_0.(5) $

那么该如何得到$ f_0 $呢?注意方程(3),如果存在这样的特殊$ k_0 $点,使其满足:

$ A_m(\vec{k_0})=0,\quad \forall \{m $ > $ 0,m\in \mathbb{I}\}.(6) $

那么立刻可以得到$ f=f_0=f(k_0) $,这样的$ k_0 $点被称为“平均值点”。但是普遍的讲,满足上述条件的$ k $点并不存在。对这个问题的解决办法就是不用单个$ k $点,而采用满足一定条件的k点的集合$ \{k\} $,利用这些$ k $点上函数值的加权平均计算$ f_0 $。也即:

$ \sum_{i=1}^{n}\alpha_iA_m(\vec{k_i})=0, m=1,\ldots,N; $

$ \sum_i\alpha_i=1.(7) $

其中$ N $可以取有限值。

利用方程(3)左右端的两个式子,可以得到:

$ \sum_{i=1}^n\alpha_if(\vec{k_i})=f_0\sum_{i=1}^n\alpha_i+\sum_{m=1}^Nf_m\sum_{i=1}^n\alpha_iA_m(\vec{k_i})+\sum_{m=N+1}^{\infty}\alpha_iA_m(\vec{k_i})f_m.(8) $

根据方程(7),可以得到

$ f_0=\sum_{i=1}^nf(\vec{k_i})-\sum_{m=N+1}^{\infty}\sum_{i=1}^n\alpha_iA_m(\vec{k_i})f_m.(9) $

考虑到$ f_m $$ m $的增大迅速减小的性质,我们可以近似的得到$ f(\vec{k}) $的平均值:

$ f(\vec{k})\approx f_0 = \sum_{i=1}^{n} \alpha_i f(\vec{k_i}).(10) $

而将方程(9)的第二项作为可控误差。因此,如果我们可以找到一组$ k $点,使得(1)集合中的$ k $点尽量少;(2)这些$ k $点在$ N $尽量大的情况下满足方程(7),则我们进行布里渊区积分的时候可以尽可能快的得到精度较高的结果。这正是特殊$ k $点方法的要点所在。反过来讲,这也表明进行具体计算的时候我们需要对计算精度进行测试,也即保证所取$ k $点使得上式的第二项足够小。

Chadi-Cohen方法

上一节证明了$ k $点的可行性。Chadi和Cohen首先提出了一套可以得出这些特殊$ k $点的方法[1]。

首先找出两个特殊$ k $点——$ \vec{k_1}, \vec{k_2} $,二者分别在$ \{N_1\} $$ \{N_2\} $的情况下满足
     $ A_m(\vec{k})=0. $
然后通过这两个$ k $点构造新的$ k $点集合:
     $ \vec{k_i}=\vec{k_1}+T_i\vec{k_2}. $
且权重为$ \alpha_i=\frac{1}{n_G} $。下面证明$ \vec{k_i} $$ \{N_1\}\bigcup\{N_2\} $的情况下仍然满足方程(7):

根据$ k_1 $$ k_2 $的定义,可知对于$ m\in\{N_1\} $$ m\in\{N_2\} $

$ A_m(\vec{k_1})A_m(\vec{k_2})=0. $

也即

$ (\sum_{|R|=C_m}e^{i\vec{k_1}\cdot\vec{R}})(\sum_{|R|=C_m}e^{i\vec{k_2}\cdot\vec{R}})=0 $

上式等价于

$ (\sum_{|R|=C_m}e^{i\vec{k_1}\cdot\vec{R}})(\sum_i e^{i\vec{k_2}\cdot T_i\vec{R}})\Rightarrow(\sum_{|R|=C_m}e^{i\vec{k_1}\cdot\vec{R}})(\sum_le^{iT_l\vec{k_2}\cdot\vec{R}}) $
$ \Rightarrow\sum_l\sum_{|R|=C_m}e^{i\vec{k_1}+T_l\vec{k_2}}=0\Rightarrow \sum_lA_m(\vec{k_1}+T_l\vec{k_2})=0\Rightarrow \sum_lA_m(\vec{k_l})=0. $

因此可以用这种方法产生一系列$ k $点,用以计算布里渊区内的积分。如果此时的精度不够,则利用同样的方法继续生成新的$ k $点集合:

$ \vec{k_{ii}}=\vec{k_i}+T_i\vec{k_3} $

其中$ \vec{k_3} $为在$ m\in\{N_3\} $情况下满足$ Am(\vec{k})=0 $的特殊$ k $点。从而改进精度。

事实上,如果考虑体系的对称性,则$ \{k\} $中的$ k $点数目可以极大的减小。也就是说,对于给定的点$ \vec{k_2} $,可以找出其波矢群$ {\vec{k_2}^{*}} $,阶数为$ n^* $,那么实际上按上述方法构造出来的$ \vec{k_i} $只有$ n_G/n^{*} $个不同的点,此时各点上的权重为$ \alpha_i=n^{*}/n_G $。更进一步,通过点群的全部对称操作,可以将全部的$ \vec{k_i} $点转入第一布里渊区的不可约部分。如果$ \vec{k_i} $点的重叠度(即第一布里渊区不可约部分中占有同样位置的$ k $点个数)是$ n_{\vec{k_i}} $,则在最后的计算中,这个点的权重为

$ \alpha_{\vec{k_i}}=\frac{n_{\vec{k_i}}}{\sum_j n_{\vec{k_j}}}.(11) $
 

Monkhorst-Pack方法

上述Chadi-Cohen方法非常巧妙,但是在具体的应用上必须首先确定2$ \sim $3个性能比较好的$ k $点,由此构建出的$ k $点集合才拥有比较高的效率和精度。因此,对于每一个具体问题,在计算之前都必须经过相当的对称性上的分析。对于编写程序而言,这是一件很麻烦的事情。那么,是否存在一种比较简易的产生$ k $点网格的方法,同时又满足方程(7)呢?答案是肯定的,这就是通常所说的Monkhorst和Pack方法[2]。

晶体中的格点$ \vec{R} $总可以表述为$ R=R_1\vec{a_1}+R_2\vec{a_2}+R_3\vec{a_3} $,其中$ \vec{a_i} $是实空间三个方向上的基矢。Monkhorst和Pack建议按如下方法划分布里渊区

$ u_r=(2r-q+1)/2q, 1\leq r\leq q.(12) $

$ k $点写为分量形式,则可得到如下表达式

$ \vec{k}_{prs}=u_p\vec{b_1}+u_r\vec{b_2}+u_s\vec{b_3}.(13) $

其中$ \vec{b_i} $是倒空间的基矢。与Chadi-Cohen方法相似,Monkhorst和Pack定义函数$ A_m $为:

$ A_m(\vec{k})=\sqrt{N_m}\sum_{|R|=C_m} e^{i\vec{k}\cdot\vec{R}},\quad m $ > $ 1 $
$ A_1(\vec{k})=1.(14) $

则相应于Chadi-Cohen方法中的$ \Omega/(2\pi)^3\int_{BZ}A_m(\vec{k})A_n(\vec{k})d\vec{k} $,我们可以计算在如方程(12)所生成的离散化的网格点的相同的量:

$ S_{mn}(q)=1/q^3\sum_{p,r,s=1}^{q}A_m(\vec{k}_{prs})A_n(\vec{k}_{prs})=\frac{1}{\sqrt{N_mN_n}}\sum_{a=1}^{N_m}\sum_{b=1}^{N_n}\Pi_{j=1}^3W_j^{ab}(q).(15) $

其中

$ W_j^{ab}(q)=\frac{1}{q}\sum_{r=1}^{q}e^{i\pi(2r-q+1)/q\cdot(R_j^b-R_j^a)}.(16) $

注意到$ R_j^a, R_j^b,\quad (j=1,2,3) $都是整数,因此可以算出$ W_j^{ab}(q) $

$ W_j^{ab}(q)=1, |R_j^b-R_j^a=0,2q,4q,\ldots; $
                     $ =(-1)^{q+1}, |R_j^b-R_j^a|=q,3q,5q,\ldots; $
                     $ =0, {\rm otherwise}.(17) $

其中第三种情况是因为$ W_j^{ab}(q) $是奇函数。引入限制条件:$ |R_j^a| $ <$ q/2 $;$ |R_j^b| $ < $ q/2 $则可得:

$ S_{mn}(q)=\delta_{mn}. $

也即$ A_m $$ k $点网格上是正交的。与Chadi-Cohen方法类似,将函数$ f(\vec{k}) $$ A_m $展开:

$ f(\vec{k})=\sum_{m=1}f_mA_m(\vec{k}).(18) $

同时左乘$ A_m^{*}(\vec{k}) $并在布里渊区内积分,可得

$ f_m=\Omega/(2\pi)^3\int_{BZ}A_m^{\*}(\vec{k})f(\vec{k})d\vec{k}.(19) $

因为$ A_1(\vec{k})=1 $,所以从方程(19)可得

$ f=\int_{BZ}f(\vec{k})d\vec{k}=(8\pi^3/\Omega)f_1.(20) $

忽略前面的常数因子,可以看到Monkhorst-Pack方法中的表达式与Chadi-Cohen方法完全一样。现在将$ f_1 $的表达式代入上述方程,则

$ \overline{f}=(8\pi^3/\Omega)(\Omega/8\pi^3)\int_{BZ}A_1^{*}(\vec{k})f(\vec{k})d\vec{k}&nbsp;=\int_{BZ}A_1^{*}(\vec{k})\sum_{m=1}^{\infty}f_mA_m(\vec{k})=\sum_m \sqrt{N_m}f_m S_{m1}(q).(21) $

因此

$ \overline{f}=f_1+\sum_{m=2}f_m\sqrt{N_m}S_{m1}(q).(22) $

其中

$ S_{m1}(q)=(-1)^{(q+1)(R_1+R_2+R_3)/q},R_j=nq,j=1,2,3, n\in \mathbb{I}; $
                 $ =0, {\rm otherwise}.(23) $

与我们在Chadi-Cohen方法中看到的一样,$ f(\vec{k}) $在第一布里渊区的平均值可以用$ f_1 $近似(在Chadi-Cohen方法中是$ f_0 $)。而且误差(方程右边第二项)可控,即可以通过增加$ k $点密度$ q $的方法提高精度。这是因为$ q $增大,根据上面所述$ S_{m1}(q) $的取值可知,在$ R_j $更大的时候仍能保证方程(7)成立。

但是根据方程(3)可得

$ f_m=(1/q)\sum_{p,r,s=1}^qf(\vec{k}_{prs})A_m(\vec{k}_{prs})\Rightarrow f_1=(1/q^3)\sum_{p,r,s=1}^qf(\vec{k}_{prs}).(24) $

如果$ q $值取得比较大,那么所需计算的$ k $点数目就会非常大,如何提高Monkhorst-Pack方法的效率呢?考虑到体系的对称性,则$ k $点的数目会大大的减少。重新写出$ f_1 $如下:

$ f_1=(1/q^3)\sum_{j=1}^{P(q)}\omega_jf(\vec{k_j}).(25) $

其中$ \omega_j $是体系所属点群阶数与$ k_j $点的波矢群阶数的比值:$ \omega_j=n_G/n^{*} $$ P(q) $是对所有$ k $点进行对称及平移操作后第一布里渊区中所有不重合的$ k $点数。进一步考虑不可约部分,那么通过改变$ \omega_j $(变为$ \omega_j\times j $,其中$ j $定义见上节)可以进一步减少$ P(q) $。因为处于高对称位置上的$ k $点其波矢群阶数也比较高,因此相应的这些高对称$ k $点的权重就比较小。这也是为什么在VASP的OUTCAR文件中高对称的$ k $点权重比较小的理论根本,也是特殊点法尽量避开高对称点的原因所在。与Chadi-Cohen方法一样,$ P(q) $的大小是Monkhorst-Pack方法效率高低的重要标志。文献中给出了bcc以及fcc两种格子中的$ P(q) $

BCC:

$ P(q)=q(q+4)(q+8)/192, {\rm mod}(q/2)=0; $
      $ =(q+2)(q+4)(q+6)/192, {\rm mod}(q/2)\ne 0. $

FCC:

$ P(q)=q(q+2)(q+4)/96, {\rm mod}(q/2)=0; $
      $ =(q+2)(q^2+4q+2), {\rm mod}(q/2)\ne 0. $

可以看出,即使对于较大的$ q $值,$ P(q) $也是比较小的,因此Monkhorst-Pack方法效率是比较高的。

应该注意的是,Monkhorst-Pack方法的关键一点是将三维空间的问题转化为三个独立的一维问题。因此,对于六角格子或者单斜格子,基矢之间不正交,上述Monkhorst-Pack方法并不适用,而必须加以修改[3]。以六角格子为例,Pack指出$ k $点网格应按下述方法生成[4]:

$ u_p=u_r=(p-1)/q_a, p,r\in [1,q_a]; u_s=(2s-q_c-1)/2q_c, s\in [1,q_c]. (28) $

也即$ a $轴和$ c $轴分别设置。相应的,$ P(q) $的大小可计算如下:

$ P_a(q_a)=(\alpha+1)(3\alpha+\beta)+\delta_{\beta0},\quad\beta={\rm mod}(q_a,6),\alpha=(q_a-\beta)/6.(29) $

$ P_c(q_c)=q_c/2,\quad {\rm mod}(q_c,2)=0; $
             $ =(q_c+1)/2,\quad {\rm mod}(q_c,2)\ne0.(30) $

上述生成$ k $点的方法对应于VASP手册中对于$ k $点设置的建议“对于六方体系应该将$ k $点置于原点处”。

需要强调的是,我们在以上所讨论的所有对称性均指纯旋转操作,也即点群对称性。因此,对于同属一种晶系而属于不同空间群的两种体系而言,其操作可能并不一致。

Chadi-Cohen方法的实例

Cunningham[5]对于二维情况依照Chadi-Cohen方法分别生成了$ k $点集合。我们选择四方格子和正方格子两种情况进行具体的分析。

四方格子

实空间和倒空间的基矢及格点坐标分别为:

$ \vec{a_1}=a(1,0); \vec{a_2}=a(0,\beta) \beta &lt; 1, \vec{R}=a(l,n\beta); $
$ \vec{b_1}=(2\pi/a)(1,0); \vec{b_2}=(2\pi/a)(0,1/\beta); $
$ \vec{K}=(2\pi/a)(k,n/\beta).(32) $

选择$ \vec{k_1^0}=(\pi/a)[1/2,1/2\beta] $$ l $$ n $为奇数时$ A_m(\vec{k})=0 $)以及$ \vec{k_2^0}=(\pi/a)[1/4,1/(4\beta)] $$ l/2 $$ n/2 $为奇数时$ A_m(\vec{k})=0 $)这个格子的对称操作为$ \{\hat{E},\hat{c_2},\hat{\sigma_v^1},\hat{\sigma_v^2}\} $,按照Chadi-Cohen方法,可以构建$ k_i $点如下:

$ \vec{k_1}=\vec{k_1^0}+\hat{E}\vec{k_2^0}=[1/2,1/(2\beta)]+[1/4,1/(4\beta)]=[3/4,3/(4\beta)]; $
$ \vec{k_2}=\vec{k_1^0}+\hat{c_2}\vec{k_2^0}=[1/2,1/(2\beta)]+[-1/4,-1/(4\beta)]=[1/4,1/(4\beta)]; $
$ \vec{k_3}=\vec{k_1^0}+\hat{\sigma_v^1}\vec{k_2^0}=[1/2,1/(2\beta)]+[-1/4,1/(4\beta)]=[1/4,3/(4\beta)]; $
$ \vec{k_4}=\vec{k_1^0}+\hat{\sigma_v^2}\vec{k_2^0}=[1/2,1/(2\beta)]+[1/4,-1/(4\beta)]=[3/4,1/(4\beta)].(32) $

每个$ k $点的权重$ \alpha_i=1/4 $

正方格子

将上例中的$ \beta=1 $,则四方格子转变为正方格子。两种情况最主要的不同是布里渊区不可约部分有了变化,从上式可以看出,在正方格子下,$ \beta=1 $,$ \vec{k_3} $$ \vec{k_4} $重合。因此只有三个不同的$ k $点,每个$ k $点的权重为$ \alpha_3=1/2;\alpha_1=\alpha_2=1/4 $,而且$ \sum_{i=1}^3\alpha_i=1 $

利用特殊$ k $点计算电荷密度

将Bloch函数用Wannier函数展开,有[6]:

$ \Psi_{\vec{k}}(\vec{r})=\frac{1}{\sqrt{N}}\sum_me^{i\vec{k}\cdot\vec{R_m}}a(\vec{r}-\vec{R_m}).(34) $

则在给定$ k $点的电荷密度为:

$ \rho_{\vec{k}}(\vec{r})=\Psi_{\vec{k}}^{*}(\vec{r})\Psi_{\vec{k}}(\vec{r})=\frac{1}{N}\sum_{mn}e^{i\vec{k}\cdot(\vec{R_m}-\vec{R_n})}a(\vec{r}-\vec{R_m})a^{*}(\vec{r}-\vec{R_n}).(35) $

$ \rho(\vec{r})=\int_{BZ}\rho_{\vec{k}}(\vec{r})d\vec{k}.(36) $

我们重写$ \rho_{\vec{k}}(\vec{r}) $如下:

$ \rho_{\vec{k}}(\vec{r})=\frac{1}{N}\sum_m|a(\vec{r}-\vec{R_m})^2|+\frac{1}{N}\sum_j^{'}\sum_me^{i\vec{k}\cdot\vec{R_j}}a(\vec{r}-\vec{R_m})a^{*}(\vec{r}+\vec{R_j}-\vec{R_m}).(37) $

其中求和号中的$ ' $表明$ \vec{R_j}\ne0 $而且$ \vec{R_j}=\vec{R_m}-\vec{R_n} $。因此,考虑到对称性,$ \rho_{\vec{k}}(\vec{r}) $又可写为:

$ \rho_{\vec{k}}(\vec{r})=\frac{1}{n_G}\sum_{T_i}\rho_{T_i\vec{k}}(\vec{r}) $
$ \Rightarrow\rho_{\vec{k}}(\vec{r})=\frac{1}{Nn_G}\sum_{T_i}\sum_m|a(\vec{r}-\vec{R_m})^2|&amp;+\frac{1}{Nn_G}\sum_j^{'}\sum_m\sum_{T_i}e^{i\vec{k}\cdot\vec{R_j}}a(\vec{r}-\vec{R_m})a^{*}(\vec{r}+\vec{R_j}-\vec{R_m}).(38) $

上式中,第一项与$ T_i $$ \vec{k} $无关,相当于Chadi-Cohen方法中的$ f_0 $。而第二项中因为对所有的$ j $求和,因此可以将这一项写为如下形式:

$ F(\vec{r})=\frac{1}{Nn_G}\sum_j^{'}e^{i\vec{k}\cdot\vec{R_j}}\sum_m\sum_{T_l}a(\vec{r}-\vec{R_m})a^{*}(\vec{r}-T_l\vec{R_j}-\vec{R_m})=\frac{1}{Nn_G}\sum_j^{'}e^{i\vec{k}\cdot\vec{R_j}}\sum_mS_m(\vec{r}).(39) $

上式中$ S_m(\vec{r}) $$ R_j $无关,且随$ |T_l\vec{R_j}-\vec{R_m}| $增大而递减,相当于$ f_m $。因此$ \rho_{\vec{k}}(\vec{r}) $可写为

$ \rho_{\vec{k}}(\vec{r})=f_0+\sum_m\sum_{|R_j|=C_m}e^{i\vec{k}\cdot\vec{R_j}}f_m=f_0+\sum_mA_m(\vec{k})f_m.(40) $

如果存在$ \vec{k_0} $,满足$ A_m(\vec{k_0})=0,m=1,2,\ldots $,则立即可以得到

$ \rho(\vec{r})=f_0=\frac{1}{N}\sum_m|a(\vec{r}-\vec{R_m})^2|=\rho_{\vec{k_0}}(\vec{r}).(41) $

但是普遍的讲,这样的$ \vec{k_0} $并不存在。例如,在fcc格子中考虑第一、二、三近邻,写出$ A_m(\vec{k}) $

$ \cos k_1\cos k_2+\cos k_1 \cos k_3+\cos k_2\cos k_3=0; $
$ \cos2k_1+\cos2k_2+\cos2k_3=0; $
$ \cos2k_1\cos k_2\cos k_3+\cos k_1\cos2k_2\cos k_3+\cos k_1\cosk_2\cos2k_3=0.(42) $

不存在单独的$ \vec{k_0} $点同时满足上述三个方程。因此,需要寻找一系列特殊的$ k $点,满足

$ \sum_{i=1}^n\sum_{|R_j|=C_m}\alpha_ie^{i\vec{k}\cdot\vec{R_j}}=\sum_{i=1}^n\alpha_iA_m(\vec{k_i})=0; $
$ \sum_i\alpha_i=1.(43) $

$ \rho(\vec{r})=\sum_i\alpha_i\rho_{\vec{k_i}}(\vec{r}) $

Chadi和Cohen[6]采用$ \vec{k_1}=(0.5,0,0) $$ \vec{k_2}=(1.0,0.5,0) $$ \vec{k_3}=(0.5,0.5,0) $三个$ k $点计算$ \rho(\vec{r}) $的值:$ \rho(\vec{r})=1/4\rho_{\vec{k_1}}(\vec{r})+1/2\rho_{\vec{k_2}}(\vec{r})+1/4\rho_{\vec{k_3}}(\vec{r}) $,取得了较好的结果。而在文献1中,他们利用$ \vec{k_1}=(3/4,1/4,1/4) $$ \vec{k_2}=(1/4,1/4,1/4) $改进了计算结果:$ \rho(\vec{r})=3/4\rho_{\vec{k_1}}(\vec{r})+1/4\rho_{\vec{k_2}}(\vec{r}) $
 

 

参考文献

[1] D.J. Chadi, M.L. Cohen. Phys. Rev. B 8, 5747 (1973).
[2] H.J Monkhorst, J.D. Pack. Phys. Rev. B 13, 5188 (1976).
[3] D.J. Chadi. Phys. Rev. B 16, 1746 (1977).
[4] J.D. Pack, H.J. Monkhorst, Phys. Rev. B 16, 1748 (1977).
[5] S.L. Cunningham. Phys. Rev. B 10, 4988 (1974).
[6] D.J. Chadi, and M.L. Cohen. Phys. Rev. B 7, 692 (1973).