第3章赝势平面波方法(I)

第3章赝势平面波方法(I)-49-基于密度泛函理论的赝势平面波方法可以计算很大范围不同体系的基态属性，它采用了平面波来展开晶体波函数，用赝势方法作有效的近似处理。由于平面波具有标准正交化和能量单一性的特点，对任何原子都适用且等同对待空间中的任何区域，不需要修正重叠误差。因此平面波函数基组适合许多体系，其简单性使之成为求解Kohn-Sham方程的高效方案之一另外，赝势的引入可以保证计算中用较少的平面波数就可以获得较为可靠的结果。该方法具有较高的计算效率，使之日益发展成为有效的计算方法。本章首先对赝势平面波方法进行重点讨论，其次介绍了基于第一性原理计算软件一般步骤，最后结合MaterialStudio软件包应用，对锐钛矿型TiO2(101)表面及其点缺陷结构进行建模和计算。基本原理基于密度泛函理论的第一性原理计算实质是求解Kohn-Sham方程。实际求解Kohn-Sham方程时，由于原子核产生的势场项在原子中心是发散的波函数变化剧烈，需要采用大量的平面波展开，因而计算成本变得非常大，所以在计算中选取尽可能少的基函数。计算中选择的基函数与最终波函数较接近则收敛较快，当然包含的维度也应该尽量少。众所周知，根据研究对象不同，选择基函数的方法也不同的，如原子轨道线性组合法(LCA0-TB)、正交平面波法(0PW)、平面波贋势法(PW-PP)、缀加平面波法(APW)、格林函数法(KKR)、线性缀加平面波法(LAPW)、Muffin-tin轨道线性组合法(LMTO)等，选取典型代表方法在随后的章节中重点展开讨论。与LAPW,LMTO等精度较高的第一性原理计算方法比较，平面波贋势法是计算量较少的方法，适用于计算精度要求不严格，因原胞较复杂而导致计算量陡增加的体系。为此，本章将重点学习赝势平面波方法，先学习电子能带的平面波基底展开以及赝势等相关基本概念，然后再讨论赝势引入原理。平面波展开与截断能1.平面波展开平面波是自由电子气的本征函数，由于金属中离子芯与类似的电子气有很小的作用，因此很自然的选择是用它描述简单金属的电子波函数。众所周知，最简单的正交、完备的函数集是平面波e某p[i(kG)r)，这里G是原胞的倒格矢。根据晶体的空间平移对称性，布洛赫(Bloch)定理(将在第4.1.1节中说明)证明，能带电子的波函数(r,k)总是能够写成(r,k)(r)e某p(ikr)(3.1)式中k是电子波矢，(r)是具有晶体平移周期性的周期函数。对于理想晶体的计算，这是很自然的，因为其哈密顿量本身具有平移对称性，只要取它的一个原胞就行了。对于无序系统(如无定型结构的固体或液体)或表面、界面问题，只要把原胞取得足够大，以至于不影响系统的动力学性质，还是可以采用周期性边界条件的。因此，这种利用平移对称性来计算电子结构的方法，对有序和无序系统都是适用的。采用周期性边界条件后单粒子轨道波函数可-50-以用平面波基展开为(r)式中11N(G)e某p(i(KG)r)(3.2)GN是归一化因子，其中是原胞体积；这里G是原胞的倒格矢，K是第一Brillouin区的波矢，(G)是展开系数。Bloch定理表明，在对真实系统的模拟中，由于电子数目的无限性，K矢量的个数从原则上讲是无限的，每个K矢量处的电子波函数都可以展开成离散的平面波基组形式，这种展开形式包含的平面波数量是无限多的。基于计算成本的考虑，实际计算中只能取有限个平面波数。采用的具体办法是一方面由于(r)随K点的变化在K点附近是可以忽略的，因此我们可以使用K点取样通过有限个K点进行计算另一方面，为了得到对波函数的准确表示，G矢量的个数也应该是无限的,但由于对有限个数的G矢量求和已经能够达到足够的准确性，因此对G的求和可以截断成有限的。给定一个截断能2Ecut(GK)2(3.3)2m对G的求和可以限制在(GK)2/2Ecut的范围内，即要求用于展开的波函数的能量小于Ecut。当K0时，即在点，有很大的计算优势，因为这时波函数的相因子是任意的，就可以取实的单粒子轨道波函数。这样，对Fourier系数满足关系式l(G)l某(G)，利用这一点，就可以节约不少的计算时间。2.截断能选取原则为了取有限个的平面波数，通常的做法是确定一个截断能量(Energycutoff)，如图3T所示，此时函数基组并不完备，总能量计算会产生相应误差，通过增加截断能量可以减小误差幅度。为了使计算出的体系总能量达到设定精度，一般截断能量必须选取到足够高。有限平面波基组的误差可以加以校正，较好的解决方法是引入一个校正因子(correctionfactor)，由此可以在一个恒定数量基组下进行计算，即使采用了恒定的截止能量这个强制条件也可以校正相应的计算结果。进行这种校正所需要的唯一的参数就是dEtot，Etot是体系总能量，Ecut是截止能量。例如，当它的数dlnEcut图3-1截断能示意图值小于0.01eV/atom时，计算就达到了良好的收敛精度，对于大多数计算0.1eV/atom就已足够。3.平面波基展开特征用平面波基来展开电子波函数是因为用平面波基来计算有很多优点。平面波基能很方便地采用快速傅里叶变换(FFT)技术，使能量、力等的计算在实空间和倒空间快速转换，这样计算尽可能在方便的空间中进行。如前面讲到的哈密顿量中的动能项的矩阵元，在倒空间中只有对角元非零，就比实空间减少了工作量。第二，平面波基函数的具体形式并不依赖于核的坐标。这样，一方面，价电子对离子的作用力可以直接应用Hellman-Feymann定理(将在3.1.5节中进行说明)得到解析的表达式，计算显得非常方便。另一方面也使总能量的计算在不同的原子构型下有基本相同的精度。此外平面波计算的收敛性和精确性比较容易控制，因-51-为通过截断能Ecut的选择可以方便地改变平面波基的多少。当然平面波基也有缺点，一般电子轨道具有一定的局域性，而平面波是空间均匀的，因此电子轨道展开时与原子轨道基相比，平面波基的个数要多得多。为了尽量减少平面波基的个数，一般在平面波的计算中都采用赝势(peudopotential)来描述离子实与价电子之间的相互作用，使电子轨道波函数在离子实内部的分布尽量平缓些。下面将讨论赝势概念及其引入思路赝势1.赝势引入平面波函数作为展开基组具有很多优点，然而截断能的选取与具体材料体系密切相关。由于原子核与电子的库仑相互作用在靠近原子核附近具有奇异性，导致在原子核附近电子波函数将剧烈振荡。因此，需要选取较大的截断能量才能正确反映电子波函数在原子核附近的行为，这势必大大地增加计算量。另一方面，在真正反映分子或固体性质的原子间成键区域其电子波函数较为平坦。基于这些特点，将固体看作价电子和离子实的集合体，离子实部分由原子核和紧密结合的芯电子组成，价电子波函数与离子实波函数满足正交化条件，由此发展出所谓的赝势方法。1959年，基于正交化平面波方法，Phillip和Kleinman提出了贋势的概念。基本思路是适当选取一平滑赝势，波函数用少数平面波展开，使计算出的能带结构与真实的接近。换句话说，使电子波函数在原子核附近表现更为平滑，而在一定范围以外又能正确反映真实波函数的特征，如图3-2所示。所谓贋势，即在离子实内部用假想的势取代真实的势，求解波动方程时，能够保持能量本征值和离子实之间的区域的波函数的不变。原子周围的所有电子中，基本上仅有价电子具有化学活性，而相邻原子的存在和作用对芯电子状态影响不大。这样，对一个由许多原子组成的固体，坐标空间根据波函数的不同特点可分成两部分(假设存在某个截断距离rc)。rc以内的核区域，所谓的芯区。波函数由紧束缚的芯电子波函数组成对周围其它原子是否存在不敏感，即与近邻的原子的波函数相互作用很小rc以外的电子波函数(称为价电子波函数)承担周围其它原子的作用而变化明显。2.原子贋势全电子DFT理论处理价电子和芯电子时采取等同对待，而在贋势中离子芯电子是被冻结的，因此采用贋势计算固体或分子性质时认为芯电子是不参与化学成键的，在体系结构进行调整时也不涉及到离子的芯电子。在贋势近似中用较弱的贋势替代芯电子所受的强烈库仑势，得到较平缓的贋波函数，此时只需考虑价电子，在不影响计算精度情况下，可以大大降低体系相应的平面波截断能Ecut,从而降低计算量。图3-3为Si原子贋势示意图。赝原子用于描述真实原子自身性质时是不正确的，但是它对原子-原子之间相互作用的描述是近似正确的。近似程度的好坏，取决于截断距离rc的大小。rc越大，贋波函数越平缓，与真实波函数的差别越大，近似带来的误差越大；反之，rc越小，与真实波函数相等的部分就越多近似引入的误差就越小。可将真实价波函数n(r,k)看作是由贋势波函数n(r,k)和内层波函数J(r,k)线性组合，即-52-n(r,k)n(r,k)nJ(k)J(r,k)(3.4)J某其中系数nJ(k)可由正交条件dr'J(r,k)n(r,k)0确定，即nJ(k)dr'J某(r,k)n(r,k)图3-2贋波函数与势图3-3Si原子贋势示意图联合真实波函数n(r,k)所满足的薛定谔方程：[TV(r)]n(r,k)En(k)n(r,k)可得到贋波函数满足如下方程[TUp]n(r,k)En(k)n(r,k)Upn(r,k)V(r)n(r,k)dr'VR(r,r')n(r',k)(3.5)其中VR(r,r')J某(r',k)[En(k)EJ]J(r',k)JUp称为原子贋势。根据密度泛函理论，原子贋势包括离子贋势Uionp和价电子库仑势和交换pp关联势：UpUionpVH(r)V某c(r)，其中后两项VH(r序叮某c(r)可以从真实电荷密度计算，此时等于对应的全电子势V(r)和V某c。从上面可知，赝势应具有以下特征：(1)赝波函数和真实波函数具有完全相同的能量本征值En(k)，这是贋势方法的重要特点；(2)贋势第二项是排斥势，与真实的吸引势有相消趋势，因此比真实势弱；(3)赝势包括局域项，其中非局域项同时与r和r'处的贋波函数n(r,k)和n(r',k)有关，而且依赖于能量本征值En(k)。3.贋势分类从上面的推导可以看出，贋势实际上是一种算子，但可以近似的将它处理成原子间距离的简单函数，叫做局部贋势。其结果是只要给定周期性排列的原子的位置，系统的能量就可以算出来。确定贋势的方法不是唯一的，主要有以下几种：⑴经验贋势方法，利用实验数据拟合有限几个V(K)的值，主要用途是在现代从头算原子贋势自洽迭代计算中作初始值使用。模型贋势是半经验的，能用于自洽计算的原子贋势。-53-模守恒贋势是第一性原理从头算原子贋势，是核与芯电子联合产生的有效势，是从原子的薛定谔方程从头计算得到的，它可以给出价电子或类价电子的正确电荷分布，适于作自洽计算。模非守恒贋势(超软贋势)，其优点是容易选择芯区的贋势波函数减少了必须的平面波函数的数目，较大地减轻了计算工作量。目前在第一性原理计算中应用较多的为模守恒贋势(Norm-conervingpeudopotential，NCPP)和超软贋势(UltraSoftpeudopotential,USPP)两种方案。下面将分别讨论。3.1.3模守恒贋势模守恒贋势构造在第一性原理计算中，常采用由HamannDR等提出的模守恒贋势。这种贋势所对应的波函数不仅与真实势对应的波函数具有相同的能量本征值而且在rc(原子芯半径)以外与真实波函数的形状和振幅都相同(即模守恒)，且在rc以内比较平缓。采用贋势计算关键在于可以有效的对化学键的价电子进行可再现的近似，贋势与全势在超过离子实半径以外具有完全相同的函数形式。这种贋势能生成正确的电荷密度，适合作自洽计算。构造模守恒贋势基本思想是选择某个特定的电子排部状态(不一定就是基态)全部电子计算在一个孤立的原子中进行，从而得到原子价电子能量本征值和价电子波函数。选择一个离子贋势或贋波函数参数形式，通过对参数的调节，使得贋原子计算和全电子原子贋势计算采用相同的交换-相关势，在超过截止半径rc后与价电子波函数形式相同，贋势的本征值等于价电子的本征值。如果电子波函数和贋势波函数满足正交归一，两者在截止半径以外的匹配性决定了模守恒条件自动成立。模守恒贋势要求贋势波函数满足：(1)本征值与真实本征值相等；(2)没有节点；(3)在原子核区以外(rrc)与真实波函数相等；(4)在内层区(rrc)内的贋电荷与真实电荷相等，将贋波函数插入到薛定谔方程中即得对应的贋势。一般说来，小的rc移植性好，可用于不同环境，但平面波收敛慢。第一性原理模守恒贋势可分为局域和非局域两部分，UpVloc(rRv)UNL(rRv,r'Rv)(3.6)其中v是对离子势求和。考虑到原子球对称性，得用球谐函数将