http://hi.baidu.com/xijunw/item/d8bcdc3f1902fdb9623aff0e
从波的角度理解平面波, K点采样，动能截断
关于平面波PW —— 从波的角度理解平面波密度泛函
1. PW是与原子无关的基函数
量子化学中有2大计算流派，基于分子轨道的的计算(HF理论)，和基于密度泛函的计算(DFT)。HF使用“轨道-轨道”和“轨道-核”的作用来计算能量，而密度泛函使用“密度-密度”，和“密度-外势(原子核)“以及“密度->交换相关能”这样的方式来计算能量。因此，对于密DFT来说，只要密度计算正确，你就别管我用了什么基函数来对密度进行展开表示。因此，有几种不同的基函数类型。其中由固体电子理论发展而来的平面波表示应用非常广泛，因为其中平面波基函数的使用数量主要由体积决定，而与电子数目无关，因而可以应用到较大体系。
我们来看看Bloch平面波是什么样子的：
[color=green][size=2]上图就是一个平面波沿空间某个方向伸展的样子，是的，它是个正弦波！[/size][/color]
我们来看看它的数学形式。通常的教科书上会这样写：
[latex]{f(r)={\frac 1 \sqrt{\Omega}}* \exp[i*G*r]}[/latex]
其中[latex]{\Omega}[/latex]为单胞体积。
即平面波的复指数形式。
而实际它可以写成更容易理解的正弦波的形式：
[latex]f(r)={\frac 1 \sqrt{V}}* [\sin(G*r)+i*cos(G*r)][/latex]
上面仍然是复数形式，不要被吓到，后面虚数部分不过是表示正弦波的相位的一种数学形式，它是一个彻头彻尾的实的正弦波，不过其振幅f(r)的变量并非时间，而是实空间坐标r(x, y, z)。其中的G，代表平面波的频率，其值等于倒空间晶格矢量(1/a, 1/b, 1/z)乘以一个数组[i,j,k]，不同的[i,j,k]代表不同的频率值。这样，就把晶格矢量和平面波的频率通过某种形式联系了起来。[i,j,k]取值越大，对应的频率越大。在一般的教科书上，你看到的是(G+k)，而不是G，下面会谈到，G和k都是频率空间上的取值。
因此，周期空间上的体系波函数就是用不同频率的正弦波的线性组合。这个与“分子轨道是原子轨道基函数的线性组合”是类似的。平面波理论，就是要用这些个正弦波的组合来构建电子在三维晶胞内的复杂分布和形状。
周期体系的波函数与原子核无关，而是是空间坐标r的函数。这个r是什么呢？ r=h*s，其中h=[a,b,c]为晶格矢量，而s为scale系数。换句话说，这个晶格空间分成了很多均匀区间，而波函数在每个区间r上展开。对于每个区间r，都需要使用一系列频率的平面波基组。
从上面可以看到，平面波DFT并不是像Gaussian中那样，使用原子轨道基函数来展开波函数和密度。这些不同频率的正弦波的系数由原子核对它们的吸引、电子自作用和交换相关泛函计算来确定。换句话说，单胞的体积形状和原子核的位置（外势）决定了平面波的系数。有时候我们的确想得到轨道的图像，这时候就要对平面波构成的波函数进行变换，得到围绕原子的s, p, d轨道，称为局域化Localization，这是另外一个话题，将专题补充。
2. PW DFT中密度在频率空间的展开：”倒空间“概念的引入与傅立叶变换
使用平面波基函数一个神奇的好处，就是波函数的积分可以通过傅立叶变换，把实空间坐标变换到频域空间即G空间(或G+k)进行更加方便的处理。
下面是对傅立叶变换的描述的示意图(Fourier transform)，图都是从网上随便找的，是个大概意思，并不精确。
[color=green][size=2]上图中，左边是实空间中的一个一维波函数示意图，它由一系列频率不同的正弦波组成，这个波函数被一个原子核吸引。右边是这些波的频率(x轴)及每个频率上波的数量，对应于平面波的线性组合系数。右图称为频域图，其中的频率由G的取值范围决定。实空间的傅里叶变换就是倒空间，倒空间的傅里叶变换就是实空间，它们组成一个傅里叶变换对。[/size][/color]
如上图注解所述，每个格点r上的波函数在倒空间中沿G展开，这就是说，我们这里有多少G值，那么真实空间中就有多少个对应频率的正弦波基函数，而这些正弦波基函数再线性组合而成体系波函数，从而密度泛函问题转化成正弦波的系数问题，而正弦波的系数问题，又转变为倒空间里G值上的强度问题。
3. 平面波展开的高频区和低频区：G采样与(G+k)采样
从前面可知，每个格点(r,r=h*s)上的波函数是由不同频率的平面波来拟合的。所以，怎么取频率值很关键。也就是说，如何确定G？ 前面说过，G是矢量，等于倒空间晶格矢量(1/a, 1/b, 1/z)乘以一个数组[i,j,k]，不同的[i,j,k]代表不同的频率值。如果令G的[i,j,k]均为整数，以这样的频率的平面波来构建体系波函数称为在gamma点展开。gamma点采样不包含低于[0,0,0]和[1,1,1]之间的低频平面波（布里渊区）。
那么，gamma点采样，充分吗？当然不行！我们知道，对于C,H,N,O等紧束缚体系，gamma点采样通常我们认为它是足够的，而对于金属等近自由电子体系，gamma点采样不行。为什么呢？因为，你想，在紧束缚体系中，波函数紧紧包围一个电子，其起伏是快速的，相当于都是高频波。而近自由电子呢，由于离域的特性，它的波函数是缓慢变化跨越多个单胞的，相当于低频波。因此，对于金属半导体体系，我们需要低频波，[0,0,0]~[1,1,1]之间的取分数值的频率的平面波！
[0,0,0]~[1,1,1]之间对应的倒空间称为第一布里渊区，波长在这个区间的波都是频率很低的波。我们现在要在这个区间上增加采样点，换句话说，就是在基函数里增加一些频率更低的平面波。所谓K空间，其实就是第一布里渊区对应的频率区间。所谓k点采样，就是要把这些低频波增加到基函数当中去，以更好地逼近近自由电子体系里平缓变化的电子波！这个称为(G+k)采样。很多资料直接把整个G+k空间成为k空间，其中的每个点称为k波矢。所谓“采样”，就是确定用来组合波函数的平面波包含哪些频率的波，包含多少。这些平面波用于在每个实空间r上展开成波函数。
怎么在k空间设置采样点呢？首先选取的点的数量越多越好，这就是说，基函数越完备，组合出来的波函数必定越精准。其次，分割方式，可以均匀分割，成为MK方式。或者，可以选取高对称性的特殊k点。
我们来看最常用的MK方式，假如要进行2*2*2的采样方式，那就是说在每个维度把0~1均匀分割成2份并组合，
[0,0,0] [0,0,1/2] [0,1/2,0] [1/2,0,0]
[0,1/2,1/2] [1/2,0,1/2] [1/2,1/2,0] [1/2,1/2,1/2] [1,1,1]
即在第一布里渊区增加处[0,0,0][1,1,1]外的7个平面波。可以看到，k点采样的数目是N^3个，会极大地增加计算量。
而高对称点比较复杂，它与盒子的类型有关，其k点多为一些特殊的分数的组合。就采样(选取平面波基组)而言，特殊k点并不比MK方式更好。那为什么要用特殊k点呢？这是因为特殊k点所对应的实空间的r点往往是波函数变号或极值点，这样当你要画能带或态密度DOS图时，就能看到波函数在实空间某些点间的起伏变化的范围。这些点就对应特殊k点。
4. 原子核附近区域采样：截断与赝势的使用
对于G采样，我们知道[i,j,k]中的值为整数。而且值越大，对应的频率越高。对于真实体系，如果波函数的起伏变化比较大，例如对靠近原子核的区域r进行展开，那么我们就需要很多高频的波函数。实际上，随着越来越靠近原子核，所需要的高频平面波的数量急剧增加，这会大大增加计算量。
为了减少计算量，我们的策略呢，首先是只对价层电子使用平面波。内层电子呢？不使用平面波表示，相当于把它从平面波表示中剔除了。然后使用一组(l=s,p,d,...)等效的原子中心波函数表示，称为赝势。如果同时要求赝势的密度积分等于内层电子的密度积分，则称为正则守恒赝势。赝势是是以原子核的坐标为中心的，采用高斯型函数即原子轨道基函数的形式。通过拟合单原子波函数确定不同元素的赝势参数，并假定它不随化合环境的变化而变化，称为“冻核近似”。
其次，我们还要对平面波的频率进行截断。虽然只对价层电子使用平面波，但[i,j,k]的取值我们不会取到无穷，尽管取的越多，波函数越逼近于真实。关于截断，一般资料用一句话来说，叫“Kohn-Sham势V_KS随|G|增加而快速收敛”，意思是说高频平面波对核-电子作用的贡献随频率增高而迅速下降，可以截断。对于波函数，[latex](h^2/2m_e)*|G|^2[/latex]代表动能，因此，这个截断又称平面波的动能截断：
[latex]{\frac {h^2} {2m_e}}*|G|^2 \leq E_{cut}[/latex]
这样，平面波的数量就由|G|^2值的上限给出了。假如知道单胞晶格矢量[a,b,c]，即可计算得到最大G，假如是[l,m,n]，你可以统计从[1,1,1]到[l,m,n]之间有多少个取值，那就代表你需要多少个平面波。
下式可以从单胞体积估算平面波数量：
[latex]N_{PW}={\frac 1 {2\pi^2}} *\Omega*{E_{cut}}^{\frac 3 2}[/latex]
5. “变盒子”问题与平面波数量问题
在VASP， CPMD等使用周期平面波作为基函数的DFT软件中，优化盒子是个比较难实施的问题，因为伴随体积变化，若给定能量截断cutoff不变，则基函数的数目必定变化，这样，不同基函数数目计算得到的能量不具有可比性，产生了所谓基组大小不一致的问题。
在量子力学计算中，基函数数目不相等则能量不能相加减，除非基组已经收敛，即基组数目已经足够多，其多少已经不再影响能量计算精确度。不过很少有计算是达到基组极限的。在Gaussian等以原子为中心的球对称基函数情形下，这种基函数数目不匹配而产生的波函数能量误差称为基组叠加误差(BSSE, basis set superposition error)。比如，计算两个水分子的结合能binding energy，需要用2个水分子dimer的总能量减去2倍的单个水分子能量。但如果假设你计算dimer用了100个基函数，计算单水分子用了50个基函数，那么不能直接相减，而需要使用100个基函数去计算单个水分子的能量后才能相减。当然，Gaussian中处理这个问题的方案是在计算dimer的时候直接计算一个误差值，除去此误差以后的dimer能量即可用来减去单水分子的能量。
在平面波DFT中，对于盒子优化以及依赖于体积的性质比如Bulk modulus，则也存在潜在的基函数数量变化问题。如果由于cutoff选取和k点采样不充分，这样的体系能量比较就可能有问题。这个问题怎么办呢？一般要验证一下能量是否对基组平面波收敛(cutoff / k sampling)；如果由于计算量限制不能使用收敛的基组，那么CPMD和Quantum Espresso提供一种constant number of plane wave方法，即强制保持基函数的数量不变。在CPMD中，可以用如下关键词来指定所需的平面波的数量：
DENSITY CUTOFF NUMBER
192000
（上面192000仅为举例）
那么，为什么使用DENSITY CUTOFF而不使用KINETIC ENERGY CUTOFF即WAVEFUNCTION CUTOFF呢？这是因为二者其实是关联的。在CPMD中，DENSITY CUTOFF默认是WAVEFUNCTION CUTOFF的4倍。可以使用DUAL关键词修改这个倍数，在使用超软赝势时 DULA应至少8～10。更重要的是，DENSITY CUTOFF所需要的平面波数量远大于WAVEFUNCTION CUTOFF所需数量，因此，可以看出，平面波code实施当中一般是有动能截断和密度截断两个标准的，并且一般是以后者DENSITY CUTOFF为标准来确定所需平面波数量的。
1. PW是与原子无关的基函数
量子化学中有2大计算流派，基于分子轨道的的计算(HF理论)，和基于密度泛函的计算(DFT)。HF使用“轨道-轨道”和“轨道-核”的作用来计算能量，而密度泛函使用“密度-密度”，和“密度-外势(原子核)“以及“密度->交换相关能”这样的方式来计算能量。因此，对于密DFT来说，只要密度计算正确，你就别管我用了什么基函数来对密度进行展开表示。因此，有几种不同的基函数类型。其中由固体电子理论发展而来的平面波表示应用非常广泛，因为其中平面波基函数的使用数量主要由体积决定，而与电子数目无关，因而可以应用到较大体系。
我们来看看Bloch平面波是什么样子的：
[color=green][size=2]上图就是一个平面波沿空间某个方向伸展的样子，是的，它是个正弦波！[/size][/color]
我们来看看它的数学形式。通常的教科书上会这样写：
[latex]{f(r)={\frac 1 \sqrt{\Omega}}* \exp[i*G*r]}[/latex]
其中[latex]{\Omega}[/latex]为单胞体积。
即平面波的复指数形式。
而实际它可以写成更容易理解的正弦波的形式：
[latex]f(r)={\frac 1 \sqrt{V}}* [\sin(G*r)+i*cos(G*r)][/latex]
上面仍然是复数形式，不要被吓到，后面虚数部分不过是表示正弦波的相位的一种数学形式，它是一个彻头彻尾的实的正弦波，不过其振幅f(r)的变量并非时间，而是实空间坐标r(x, y, z)。其中的G，代表平面波的频率，其值等于倒空间晶格矢量(1/a, 1/b, 1/z)乘以一个数组[i,j,k]，不同的[i,j,k]代表不同的频率值。这样，就把晶格矢量和平面波的频率通过某种形式联系了起来。[i,j,k]取值越大，对应的频率越大。在一般的教科书上，你看到的是(G+k)，而不是G，下面会谈到，G和k都是频率空间上的取值。
因此，周期空间上的体系波函数就是用不同频率的正弦波的线性组合。这个与“分子轨道是原子轨道基函数的线性组合”是类似的。平面波理论，就是要用这些个正弦波的组合来构建电子在三维晶胞内的复杂分布和形状。
周期体系的波函数与原子核无关，而是是空间坐标r的函数。这个r是什么呢？ r=h*s，其中h=[a,b,c]为晶格矢量，而s为scale系数。换句话说，这个晶格空间分成了很多均匀区间，而波函数在每个区间r上展开。对于每个区间r，都需要使用一系列频率的平面波基组。
从上面可以看到，平面波DFT并不是像Gaussian中那样，使用原子轨道基函数来展开波函数和密度。这些不同频率的正弦波的系数由原子核对它们的吸引、电子自作用和交换相关泛函计算来确定。换句话说，单胞的体积形状和原子核的位置（外势）决定了平面波的系数。有时候我们的确想得到轨道的图像，这时候就要对平面波构成的波函数进行变换，得到围绕原子的s, p, d轨道，称为局域化Localization，这是另外一个话题，将专题补充。
2. PW DFT中密度在频率空间的展开：”倒空间“概念的引入与傅立叶变换
使用平面波基函数一个神奇的好处，就是波函数的积分可以通过傅立叶变换，把实空间坐标变换到频域空间即G空间(或G+k)进行更加方便的处理。
下面是对傅立叶变换的描述的示意图(Fourier transform)，图都是从网上随便找的，是个大概意思，并不精确。
[color=green][size=2]上图中，左边是实空间中的一个一维波函数示意图，它由一系列频率不同的正弦波组成，这个波函数被一个原子核吸引。右边是这些波的频率(x轴)及每个频率上波的数量，对应于平面波的线性组合系数。右图称为频域图，其中的频率由G的取值范围决定。实空间的傅里叶变换就是倒空间，倒空间的傅里叶变换就是实空间，它们组成一个傅里叶变换对。[/size][/color]
如上图注解所述，每个格点r上的波函数在倒空间中沿G展开，这就是说，我们这里有多少G值，那么真实空间中就有多少个对应频率的正弦波基函数，而这些正弦波基函数再线性组合而成体系波函数，从而密度泛函问题转化成正弦波的系数问题，而正弦波的系数问题，又转变为倒空间里G值上的强度问题。
3. 平面波展开的高频区和低频区：G采样与(G+k)采样
从前面可知，每个格点(r,r=h*s)上的波函数是由不同频率的平面波来拟合的。所以，怎么取频率值很关键。也就是说，如何确定G？ 前面说过，G是矢量，等于倒空间晶格矢量(1/a, 1/b, 1/z)乘以一个数组[i,j,k]，不同的[i,j,k]代表不同的频率值。如果令G的[i,j,k]均为整数，以这样的频率的平面波来构建体系波函数称为在gamma点展开。gamma点采样不包含低于[0,0,0]和[1,1,1]之间的低频平面波（布里渊区）。
那么，gamma点采样，充分吗？当然不行！我们知道，对于C,H,N,O等紧束缚体系，gamma点采样通常我们认为它是足够的，而对于金属等近自由电子体系，gamma点采样不行。为什么呢？因为，你想，在紧束缚体系中，波函数紧紧包围一个电子，其起伏是快速的，相当于都是高频波。而近自由电子呢，由于离域的特性，它的波函数是缓慢变化跨越多个单胞的，相当于低频波。因此，对于金属半导体体系，我们需要低频波，[0,0,0]~[1,1,1]之间的取分数值的频率的平面波！
[0,0,0]~[1,1,1]之间对应的倒空间称为第一布里渊区，波长在这个区间的波都是频率很低的波。我们现在要在这个区间上增加采样点，换句话说，就是在基函数里增加一些频率更低的平面波。所谓K空间，其实就是第一布里渊区对应的频率区间。所谓k点采样，就是要把这些低频波增加到基函数当中去，以更好地逼近近自由电子体系里平缓变化的电子波！这个称为(G+k)采样。很多资料直接把整个G+k空间成为k空间，其中的每个点称为k波矢。所谓“采样”，就是确定用来组合波函数的平面波包含哪些频率的波，包含多少。这些平面波用于在每个实空间r上展开成波函数。
怎么在k空间设置采样点呢？首先选取的点的数量越多越好，这就是说，基函数越完备，组合出来的波函数必定越精准。其次，分割方式，可以均匀分割，成为MK方式。或者，可以选取高对称性的特殊k点。
我们来看最常用的MK方式，假如要进行2*2*2的采样方式，那就是说在每个维度把0~1均匀分割成2份并组合，
[0,0,0] [0,0,1/2] [0,1/2,0] [1/2,0,0]
[0,1/2,1/2] [1/2,0,1/2] [1/2,1/2,0] [1/2,1/2,1/2] [1,1,1]
即在第一布里渊区增加处[0,0,0][1,1,1]外的7个平面波。可以看到，k点采样的数目是N^3个，会极大地增加计算量。
而高对称点比较复杂，它与盒子的类型有关，其k点多为一些特殊的分数的组合。就采样(选取平面波基组)而言，特殊k点并不比MK方式更好。那为什么要用特殊k点呢？这是因为特殊k点所对应的实空间的r点往往是波函数变号或极值点，这样当你要画能带或态密度DOS图时，就能看到波函数在实空间某些点间的起伏变化的范围。这些点就对应特殊k点。
4. 原子核附近区域采样：截断与赝势的使用
对于G采样，我们知道[i,j,k]中的值为整数。而且值越大，对应的频率越高。对于真实体系，如果波函数的起伏变化比较大，例如对靠近原子核的区域r进行展开，那么我们就需要很多高频的波函数。实际上，随着越来越靠近原子核，所需要的高频平面波的数量急剧增加，这会大大增加计算量。
为了减少计算量，我们的策略呢，首先是只对价层电子使用平面波。内层电子呢？不使用平面波表示，相当于把它从平面波表示中剔除了。然后使用一组(l=s,p,d,...)等效的原子中心波函数表示，称为赝势。如果同时要求赝势的密度积分等于内层电子的密度积分，则称为正则守恒赝势。赝势是是以原子核的坐标为中心的，采用高斯型函数即原子轨道基函数的形式。通过拟合单原子波函数确定不同元素的赝势参数，并假定它不随化合环境的变化而变化，称为“冻核近似”。
其次，我们还要对平面波的频率进行截断。虽然只对价层电子使用平面波，但[i,j,k]的取值我们不会取到无穷，尽管取的越多，波函数越逼近于真实。关于截断，一般资料用一句话来说，叫“Kohn-Sham势V_KS随|G|增加而快速收敛”，意思是说高频平面波对核-电子作用的贡献随频率增高而迅速下降，可以截断。对于波函数，[latex](h^2/2m_e)*|G|^2[/latex]代表动能，因此，这个截断又称平面波的动能截断：
[latex]{\frac {h^2} {2m_e}}*|G|^2 \leq E_{cut}[/latex]
这样，平面波的数量就由|G|^2值的上限给出了。假如知道单胞晶格矢量[a,b,c]，即可计算得到最大G，假如是[l,m,n]，你可以统计从[1,1,1]到[l,m,n]之间有多少个取值，那就代表你需要多少个平面波。
下式可以从单胞体积估算平面波数量：
[latex]N_{PW}={\frac 1 {2\pi^2}} *\Omega*{E_{cut}}^{\frac 3 2}[/latex]
5. “变盒子”问题与平面波数量问题
在VASP， CPMD等使用周期平面波作为基函数的DFT软件中，优化盒子是个比较难实施的问题，因为伴随体积变化，若给定能量截断cutoff不变，则基函数的数目必定变化，这样，不同基函数数目计算得到的能量不具有可比性，产生了所谓基组大小不一致的问题。
在量子力学计算中，基函数数目不相等则能量不能相加减，除非基组已经收敛，即基组数目已经足够多，其多少已经不再影响能量计算精确度。不过很少有计算是达到基组极限的。在Gaussian等以原子为中心的球对称基函数情形下，这种基函数数目不匹配而产生的波函数能量误差称为基组叠加误差(BSSE, basis set superposition error)。比如，计算两个水分子的结合能binding energy，需要用2个水分子dimer的总能量减去2倍的单个水分子能量。但如果假设你计算dimer用了100个基函数，计算单水分子用了50个基函数，那么不能直接相减，而需要使用100个基函数去计算单个水分子的能量后才能相减。当然，Gaussian中处理这个问题的方案是在计算dimer的时候直接计算一个误差值，除去此误差以后的dimer能量即可用来减去单水分子的能量。
在平面波DFT中，对于盒子优化以及依赖于体积的性质比如Bulk modulus，则也存在潜在的基函数数量变化问题。如果由于cutoff选取和k点采样不充分，这样的体系能量比较就可能有问题。这个问题怎么办呢？一般要验证一下能量是否对基组平面波收敛(cutoff / k sampling)；如果由于计算量限制不能使用收敛的基组，那么CPMD和Quantum Espresso提供一种constant number of plane wave方法，即强制保持基函数的数量不变。在CPMD中，可以用如下关键词来指定所需的平面波的数量：
DENSITY CUTOFF NUMBER
192000
（上面192000仅为举例）
那么，为什么使用DENSITY CUTOFF而不使用KINETIC ENERGY CUTOFF即WAVEFUNCTION CUTOFF呢？这是因为二者其实是关联的。在CPMD中，DENSITY CUTOFF默认是WAVEFUNCTION CUTOFF的4倍。可以使用DUAL关键词修改这个倍数，在使用超软赝势时 DULA应至少8～10。更重要的是，DENSITY CUTOFF所需要的平面波数量远大于WAVEFUNCTION CUTOFF所需数量，因此，可以看出，平面波code实施当中一般是有动能截断和密度截断两个标准的，并且一般是以后者DENSITY CUTOFF为标准来确定所需平面波数量的。
No comments:
Post a Comment