DFT计算的基本要素¶

倒易空间和k点¶

所有的一切来自于bloch定理,为了简化起见,让我们从一维的情况开始考虑,对于一个一维空间中自由运动的电子,其动量算符的本征方程为:

\[ -\frac{\hbar}{i} \frac{d}{dr} \phi(r) = p \phi(r) \]

其解为:

\[ \phi(r)= e^{-\frac{ip}{\hbar}r} \]

这个解被称作平面波,波函数可以简写为\(\phi_k(r)=e^{-ikr}\),k被称作波矢,这也是自由运动的电子的薛定谔方程的解.

现在,考虑引入一维周期性势场来模拟晶体的情况,那么,薛定谔方程的解需要进行一定的修正,不难发现,电子密度的分布一定具有周期性,所以\(|\phi(r)|\),故我们考虑分离变量,使用一个周期性函数对平面波函数进行调制:

\[ \phi_k(r)=e^{-ikr}u_k(r) \]

其中\(u_k(r+R)=u_k(r)\),即满足空间周期性,周期性势场中薛定谔方程的解一定具有这种形式,这就是bloch定理,推广一下写成三维的形式:

\[ \phi_{\mathbf{k}}(\mathbf{r})=e^{-i\mathbf{k}\cdot \mathbf{r}}u_{\mathbf{k}}(\mathbf{r}) \]

布里渊区¶

bloch定理表明,我们可以根据每一个波矢单独求解薛定谔方程,而这就是核心所在,使用k求解DFT问题会便利很多,所以接下来的一切都要围绕k展开.

矢量r的空间被称作实空间,而k的空间被称作倒易空间,倒易空间中的点被称作k点.倒易空间的单位矢量的定义依赖实空间:

\[ \mathbf{a_i}\cdot \mathbf{b_j}=2\pi\delta_{ij} \]

对于简单立方晶体,可以化简关系为:\(|\mathbf{b}|=2\pi/a\),其中a是晶格常数.

就像实空间中存在原胞这种周期性最小重复单元一样,倒易空间中存在倒易原胞,这被我们称作布里渊区,不难发现一个规律,实空间的可重复单元体积越大,布里渊区体积越小.

k空间中的积分¶

对于实际的DFT计算,其大量工作可以归结为计算以下积分:

\[ \bar{g}=\frac{V_{cell}}{(2\pi)^3}\int_{\text{BZ}}g(\mathbf{k})d\mathbf{k} \]

其中\(\text{BZ}\)表示布里渊区,\(g(k)\)是\(k\)的函数,即在布里渊区内做积分运算.

由于g函数一次的计算成本是很高的,并且大多数情况是不解析的,所以我们必须选用收敛较快的数值积分方法,例如高斯数值积分法:

\[ \int_{-1}^{1}f(x)dx=\sum_{i=1}^{n}w_if(x_i) \]

其中\(w_i\)是权重,而\(x_i\)是积分点,通过合理选择权重和积分点,可以大大提高积分的收敛速度.

在布里渊区中选择k点¶

k点就是我们积分的采样点,其采样往往使用三维网格采样的方式,即确定每个方向上的采样点数,然后计算每个采样点上的函数值,最后求和.

根据上面的讨论我们知道,如果可重复单元是超胞,那么布里渊区会相应缩小,反之,使用原胞来进行计算,布里渊区会相应增大,此时需要的k点数目也会相应增多,k点的数目越多计算结果越精确,但是计算成本也会相应增加.

书上指出了一个奇特的现象,当采用[M,M,M]这样的采样方式的时候,M和M+1的计算结果相差不大,这可以用对称性来解释,倒易空间中的积分可以用部分区域即 不可约布里渊区(IBZ) 来计算.

另一个奇特的现象是,能量的差值比能量的绝对值收敛的更快,这具有重要的实际意义,因为我们往往关注的是从一种结构到另一种结构能量的变化,故可以用更少的k点抵达收敛精度以内.

当我们要同时计算多个体系的并且要将他们进行比较的时候,这些体系的超胞规格往往是不相同的,为了让结果具有可比性,我们计算的时候需要使用相同的倒易空间中的k点密度.

金属¶

根据能带理论,金属的电子占据区域被费米面所分割,形成间断的区域,这导致g函数往往是不连续的,出现阶跃的现象,这对数值积分是致命的,往往需要使用大量的k点才能达成较好的收敛精度.

为了解决这一问题,主流的解决方案有两个:

四面体插值
模糊化处理(Smearing)

截断能¶

为了简化计算,我们需要引入一些截断操作,下面我们来说明这个截断是怎么来的,首先,我们将\(u_k(r)\)进行无穷基组展开:

\[ u_{\mathbf{k}}(\mathbf{r})=\sum_{\mathbf{G}}c_{\mathbf{G}}e^{-i\mathbf{G}\cdot \mathbf{r}} \]

于是,原本的解就可以写成:

\[ \phi_{\mathbf{k}}(\mathbf{r})=\sum_{\mathbf{G}}c_{\mathbf{G}}e^{-i(\mathbf{G}+\mathbf{k})\cdot \mathbf{r}} \]

这样,原本的波函数就可以写成平面波的叠加,每个平面波的动能为:

\[ E_{\mathbf{G}}=\frac{\hbar^2}{2m}|\mathbf{G}+\mathbf{k}|^2 \]

较高的动能意味着离核较近的电子,这部分电子不参与化学性质的体现,低动能的电子具有较强的物理意义,于是我们可以截断基组:

\[ \phi_{\mathbf{k}}(\mathbf{r})=\sum_{|\mathbf{G}+\mathbf{k}|<G_{cut}}c_{\mathbf{G}}e^{-i(\mathbf{G}+\mathbf{k})\cdot \mathbf{r}} \]

其中\(G_{cut}\)是截断能,这个能的大小决定了基组的大小,基组越大,计算结果越精确,但是计算成本也会相应增加.

当我们不指定\(G_{cut}\)的时候,程序会自动选择一个合适的值,这个值往往是化合物具有最高截断能的元素的截断能.,值得注意的时候,如果要计算化学反应带来的构型变化,所有化合物的截断能必须统一.

赝势¶

进一步的,既然核附近的电子不参与化学性质的体现,那么我们就可以拟合一个核势,这个核势代表内层电子和原子核提供的势能,从而减少自由电子的数目,这就是赝势.

对于重原子,其电子个数非常多,赝势的引入可以大大减少计算成本.对于每一种赝势,都规定了使用最小的截断能,这个最小截断能的大小区分了赝势的软硬,软赝势的截断能小,硬赝势的截断能大,对应的截断半径则是相反.Vanderbilt等人提出了一种只需要更低截断能但却具有高精度的赝势,这种赝势被称为超软赝势(USPP),他平衡了精度与效率,适合复杂体系,但是,缺点是,对于每个原子,都需要确定一系列经验参数,只有那些经过严格标定的元素是可靠的,此外,还有一种被称作投影缀加平面波（PAW）的方法,是目前最精确的赝势类方法.