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

.第3章赝势平面波方法(I)基于密度泛函理论的赝势平面波方法可以计算很大范围不同体系的基态属性，它采用了平面波来展开晶体波函数，用赝势方法作有效的近似处理。由于平面波具有正交化和能量单一性的特点，对任何原子都适用且等同对待空间中的任何区域，不需要修正重叠误差。因此平面波函数基组适合许多体系，其简单性使之成为求解Kohn-Sham方程的高效之一。另外，赝势的引入可以保证计算中用较少的平面波数就可以获得较为可靠的结果。该方法具有较高的计算效率，使之日益发展成为有效的计算方法。本章首先对赝势平面波方法进行重点讨论，其次介绍了基于第一性原理计算软件一般步骤，最后结合MaterialsStudio软件包应用，对锐钛矿型TiO2(101)表面及其点缺陷结构进行建模和计算。3.1基本原理基于密度泛函理论的第一性原理计算实质是求解Kohn-Sham方程。实际求解Kohn-Sham方程时，由于原子核产生的势场项在原子中心是发散的，波函数变化剧烈，需要采用大量的平面波展开，因而计算成本变得非常大，所以在计算中选取尽可能少的基函数。计算中选择的基函数与最终波函数较接近则收敛较快，当然包含的维度也应该尽量少。众所周知，根据研究对象不同，选择基函数的方法也不同的，如原子轨道线性组合法(LCAO-TB)、正交平面波法(OPW)、平面波赝势法(PW-PP)、缀加平面波法(APW)、格林函数法(KKR)、线性缀加平面波法(LAPW)、Muffin-tin轨道线性组合法(LMTO)等，选取典型代表方法在随后的章节中重点展开讨论。与LAPW，LMTO等精度较高的第一性原理计算方法比较，平面波赝势法是计算量较少的方法，适用于计算精度要求不严格，因原胞较复杂而导致计算量陡增加的体系。为此，本章将重点学习赝势平面波方法，先学习电子能带的平面波基底展开以及赝势等相关基本概念，然后再讨论赝势引入原理。3.1.1平面波展开与截断能1.平面波展开平面波是自由电子气的本征函数，由于金属中离子芯与类似的电子气有很小的作用，因此很自然的选择是用它描述简单金属的电子波函数。众所周知，最简单的正交、完备的函数集是平面波，这里是原胞的倒格矢。根据晶体的空间平移对称性，布洛赫(Bloch)定理(将在第4.1.1节中说明)证明，能带电子的波函数总是能够写成(3.1)式中是电子波矢，是具有晶体平移周期性的周期函数。对于理想晶体的计算，这是很自然的，因为其哈密顿量本身具有平移对称性，只要取它的一个原胞就行了。对于无序系统(如无定型结构的固体或液体)或表面、界面问题，只要把原胞取得足够大，以至于不影响系统的动力学性质，还是可以采用周期性边界条件的。因此，这种利用平移对称性来计算电子结构的方法，对有序和无序系统都是适用的。采用周期性边界条件后，单粒子轨道波函数可以用平面波基展开为(3.2)式中是归一化因子，其中是原胞体积；这里是原胞的倒格矢，是第一Brillouin区的波矢，是展开系数。Bloch定理表明，在对真实系统的模拟中，由于电子数目的无限性，矢量的个数从原则上讲是无限的，每个矢量处的电子波函数都可以展开成离散的平面波基组形式，这种展开形式包含的平面波数量是无限多的。基于计算成本的考虑，实际计算中只能取有限个平面波数。采用的具体办法是一方面由于随点的变化在点附近是可以忽略的，因此我们可以使用点取样通过有限个点进行计算。另一方面，为了得到对波函数的准确表示，矢量的个数也应该是无限的，但由于对有限个数的矢量求和已经能够达到足够的准确性，因此对的求和可以截断成有限的。给定一个截断能(3.3)对的求和可以限制在的范围内，即要求用于展开的波函数的能量小于。当时，即在点，有很大的计算优势，因为这时波函数的相因子是任意的，就可以取实的单粒子轨道波函数。这样，对Fourier系数满足关系式，利用这一点，就可以节约不少的计算时间。2.截断能选取原则为了取有限个的平面波数，通常的做法是确定一个截断能量(Energycutoff)，如图3-1所示，此时函数基组并不完备，总能量计算会产生相应误差，通过增加截断能量可以减小误差幅度。为了使计算出的体系总能量达到设定精度，一般截断能量必须选取到足够高。有限平面波基组的误差可以加以校正，较好的解决方法是引入一个校正因子(correctionfactor)，由此可以在一个恒定数量基组下进行计算，即使采用了恒定的截止能量这个强制条件也可以校正相应的计算结果。进行这种校正所需要的唯一的参数就是，Etot是体系总能量，Ecut是截止能量。例如，当它的数值小于0.01eV/atom时，计算就达到了良好的收敛精度，对于大多数计算0.1eV/atom就已足够。3.平面波基展开特征用平面波基来展开电子波函数是因为用平面波基来计算有很多优点。平面波基能很方便地采用快速傅里叶变换(FFT)技术，使能量、力等的计算在实空间和倒空间快速转换，这样计算尽可能在方便的空间中进行。如前面讲到的哈密顿量中的动能项的矩阵元，在倒空间中只有对角元非零，就比实空间减少了工作量。第二，平面波基函数的具体形式并不依赖于核的坐标。这样，一方面，价电子对离子的作用力可以直接应用Hellman-Feymann定理(将在3.1.5节中进行说明)得到解析的表达式，计算显得非常方便。另一方面也使总能量的计算在不同的原子构型下有基本相同的精度。此外平面波计算的收敛性和精确性比较容易控制，因为通过截断能Ecut的选择可以方便地改变平面波基的多少。当然平面波基也有缺点，一般电子轨道具有一定的局域性，而平面波是空间均匀的，因此电子轨道展开时与原子轨道基相比，平面波基的个数要多得多。为了尽量减少平面波基的个数，一般在平面波的计算中都采用赝势(pseudopotentials)来描述离子实与价电子之间的相互作用，使电子轨道波函数在离子实内部的分布尽量平缓些。下面将讨论赝势概念及其引入思路。3.1.2赝势4.赝势引入平面波函数作为展开基组具有很多优点，然而截断能的选取与具体材料体系密切相关。由于原子核与电子的库仑相互作用在靠近原子核附近具有奇异性，导致在原子核附近电子波函数将剧烈振荡。因此，需要选取较大的截断能量才能正确反映电子波函数在原子核附近的行为，这势必大大地增加计算量。另一方面，在真正反映分子或固体性质的原子间成键区域，其电子波函数较为平坦。基于这些特点，将固体看作价电子和离子实的集合体，离子实部分由原子核和紧密结合的芯电子组成，价电子波函数与离子实波函数满足正交化条件，由此发展出所谓的赝势方法。1959年，基于正交化平面波方法，Phillips和Kleinman提出了赝势的概念。基本思路是适当选取一平滑赝势，波函数用少数平面波展开，使计算出的能带结构与真实的接近。换句话说，使电子波函数在原子核附近表现更为平滑，而在一定范围以外又能正确反映真实波函数的特征，如图3-2所示。所谓赝势，即在离子实内部用假想的势取代真实的势，求解波动方程时，能够保持能量本征值和离子实之间的区域的波函数的不变。原子周围的所有电子中，基本上仅有价电子具有化学活性，而相邻原子的存在和作用对芯电子状态影响不大。这样，对一个由许多原子组成的固体，坐标空间根据波函数的不同特点可分成两部分(假设存在某个截断距离)。(1)以内的核区域，所谓的芯区。波函数由紧束缚的芯电子波函数组成，对周围其它原子是否存在不敏感，即与近邻的原子的波函数相互作用很小；(2)以外的电子波函数(称为价电子波函数)承担周围其它原子的作用而变化明显。5.原子赝势全电子DFT理论处理价电子和芯电子时采取等同对待，而在赝势中离子芯电子是被冻结的，因此采用赝势计算固体或分子性质时认为芯电子是不参与化学成键的，在体系结构进行调整时也不涉及到离子的芯电子。在赝势近似中用较弱的赝势替代芯电子所受的强烈库仑势，得到较平缓的赝波函数，此时只需考虑价电子，在不影响计算精度情况下，可以大大降低体系相应的平面波截断能Ecut，从而降低计算量。图3-3为Si原子赝势示意图。赝原子用于描述真实原子自身性质时是不正确的，但是它对原子-原子之间相互作用的描述是近似正确的。近似程度的好坏，取决于截断距离的大小。越大，赝波函数越平缓，与真实波函数的差别越大，近似带来的误差越大；反之，越小，与真实波函数相等的部分就越多，近似引入的误差就越小。可将真实价波函数看作是由赝势波函数和内层波函数线性组合，即(3.4)其中系数可由正交条件确定，即图3-2赝波函数与势图3-3Si原子赝势示意图联合真实波函数所满足的薛定谔方程：可得到赝波函数满足如下方程(3.5)其中称为原子赝势。根据密度泛函理论，原子赝势包括离子赝势和价电子库仑势和交换关联势：，其中后两项和可以从真实电荷密度计算，此时等于对应的全电子势和。从上面可知，赝势应具有以下特征：(1)赝波函数和真实波函数具有完全相同的能量本征值，这是赝势方法的重要特点；(2)赝势第二项是排斥势，与真实的吸引势有相消趋势，因此比真实势弱；(3)赝势包括局域项，其中非局域项同时与和处的赝波函数和有关，而且依赖于能量本征值。6.赝势分类从上面的推导可以看出，赝势实际上是一种算子，但可以近似的将它处理成原子间距离的简单函数，叫做局部赝势。其结果是只要给定周期性排列的原子的位置，系统的能量就可以算出来。确定赝势的方法不是唯一的，主要有以下几种：(1)经验赝势方法，利用实验数据拟合有限几个V(K)的值，主要用途是在现代从头算原子赝势自洽迭代计算中作初始值使用。(2)模型赝势是半经验的，能用于自洽计算的原子赝势。(3)模守恒赝势是第一性原理从头算原子赝势，是核与芯电子联合产生的有效势，是从原子的薛定谔方程从头计算得到的，它可以给出价电子或类价电子的正确电荷分布，适于作自洽计算。(4)模非守恒赝势(超软赝势)，其优点是容易选择芯区的赝势波函数，减少了必须的平面波函数的数目，较大地减轻了计算工作量。目前在第一性原理计算中应用较多的为模守恒赝势(Norm-conservingpseudopotentials，NCPP)和超软赝势(UltraSoftpseudopotentials，USPP)两种方案。下面将分别讨论。模守恒赝势1.模守恒赝势构造在第一性原理计算中，常采用由HamannDR等提出的模守恒赝势。这种赝势所对应的波函数不仅与真实势对应的波函数具有相同的能量本征值，而且在(原子芯半径)以外与真实波函数的形状和振幅都相同(即模守恒)，且在以内比较平缓。采用赝势计算关键在于可以有效的对化学键的价电子进行可再现的近似，赝势与全势在超过离子实半径以外具有完全相同的函数形式。这种赝势能生成正确的电荷密度，适合作自洽计算。构造模守恒赝势基本思想是选择某个特定的电子排部状态(不一定就是基态)全部电子计算在一个孤立的原子中进行，从而得到原子价电子能量本征值和价电子波函数。选择一个离子赝势或赝波函数参数形式，通过对参数的调节，使得赝原子计算和全电子原子赝势计算采用相同的交换-相关势，在超过截止半径后与价电子波函数形式相同，赝势的本征值等于价电子的本征值。如果电子波函数和赝势波函数满足正交归一，两者在截止半径以外的匹配性决定了模守恒条件自动成立。模守恒赝势要求赝势波函数满足：(1)本征值与真实本征值相等；(2)没有节点；(3)在原子核区以外()与真实波函数相等；(4)在内层区()内的赝电荷与真实电荷相等，将赝波函数插入到薛定谔方程中即得对应的赝势。一般说来，小的移植性好，可用于不同环境，但平面波收敛慢。第一性原理模守恒赝势可分为局域和非局域两部分，(3.6)其中是对离子势求和。考虑到原子球对称性，得用球谐函数将赝势的非局域部分写成：如果将取成半局域形式，即径向是局域的，只有角部分是非局域的：，并定义角动量的投影算符。则半局域的原子赝势可以写成如下的形式：，(3.7)为了简化计算，Kleinman和Bylander(KB)将上面半局域赝势部分用一个非局域赝势来近似：(3.8)2.模守恒赝势应用特征模守恒赝势方法可以在局域密度近似下，采用平面波基精确有效地计算固态性质，且可移植性好，但在描述局域价轨道平面波基仍然很大，因而在第Ⅰ族元素和过渡族金属中的应用受到了限制。通过优化光滑的赝波函数或赝势和增大截断半径的改进有一定效果，但模守恒条件的限制使得在一些情况下，如O原子2p或Ni原子3d轨道，很难构造出比全电子波函数更光滑的赝波函数，收敛仍然很慢。模守恒赝势最早由HamannDR等提出，后来建立一组涵盖整个周期表的参数，之后发展出的Kerker、TM赝势和Optimised赝势，都是在朝着兼顾准确性的情况下，尽可能使必须使用的平面波基底数目越少越好，平面波基底数是直接影响着所需计算量大小的量。一个赝势所需的基底数多少，可由对的收敛性来判断，即平面波截断动能用到多大时则固态计算所求得系统总能就不再改变，所需越小，也就是所谓的赝势越“软”。使用Optimised或TM赝势虽然能够把模守恒型赝势变的很“软”，但模守恒条件对于原本就已经没有节点价电子云分布的改造及最佳化的程度，与现今日渐普遍的超软赝势(它不必遵守模守恒条件)来比，节省计算的程度仍是有限的。总之，计算量的大小是取决于原子的种类这一点，是十分明确而普遍的认识，也就是说不同种类元素其势的“软硬”差异会令人明显感觉到。超软赝势1.超软赝势构造对于过渡族元素和第一周期元素，模守恒赝势不能明显降低所需平面波截断能Ecut。Vanderbilt提出了超软赝势，其赝波函数在核心范围是被作成尽可能平滑，可以大幅度地减少截断能，即可使计算所需的平面波函数基组更少。就技术上而言，这是靠放宽模守恒的要求，采用广义的正交条件来达成的。为了重建整个总的电子密度，波函数平方所得到电荷密度必须在核心范围再附加额外的密度进去。这个电子云密度由此就被分成两个部分，第一部分是一个延伸在整个单位晶胞平滑部分，第二部分是一个局域化在核心区域的自旋部分。前面所提的附加部分是只出现在电子密度，并不在波函数。超软赝势中总能量与采用其它赝势平面波方法时相同，非定域势表达如下(3.9)式中投影算符和系数分别表征赝势和原子种类的差别，指数对应于一个原子位置。总能量用电子密度可以表示为：(3.10)式中是波函数，是严格位于芯区的附加函数。超软赝势完全由定域部分，和系数，，确定。赝势是通过引进一系列正交条件来建立的，，是哈密顿重叠算符，可以表示为，系数是通过对积分得到。从而，超软赝势的Kohn-Sham方程可以写为：，可以表示为动能和定域势能之和(3.11)(3.12)2.应用特点超软赝势产生算法保证了在预先选择的能量范围内会有良好的散射性质，这导致了赝势更好的转换性与精确性。超软赝势通常也借着把多套每个角动量通道当作价电子来处理浅的内层电子态。这也会使精确度跟转换性更加提升，虽然计算代价会比较高。与模守恒赝势对比，不同之处在于在超软赝势中存在重叠算符，波函数与有关。而且投影算符函数数量要比模守恒赝势中大两倍多。与附加电荷相关的一系列计算可以在实空间中进行，这与函数中定域势的性质有关，而多余的步骤不会对计算效率产生较大的影响。在Laasonen文献中提供了超软赝势计算的详细方法以及总能量微分表达式。3.1.3Hellmann-Feynman力Hellmann和Feynman在量子力学框架下给出了作用于离子实上(位置坐标为)的力。离子受的力为总能对离子位置的偏导，(3.13)E作为系统哈密顿量的能量本征值，满足Kohn-Sham方程，可以得到：(3.14)将式(3.14)代入(3.13)得(3.15)由于是一个归一化常数，上式的第一项等于零。最终得到作用在离子上的力(3.16)这就是著名的Hellmann-Feynman定理Hellmann-Feynman定理计算出的力是和电子波函数相联系的，它的误差与波函数误差的一级修正量成正比，只有波函数非常接近真实的本征态时这个力才是精确的。所以在计算时需要同时考虑到离子弛豫和电荷密度自洽，即，离子在受力后到达一个新的位置，此时电子也需要接近瞬间基态，然后在新的离子位置和新的电子密度下进行计算，直至总能到达局部极小值。在得出离子受的受力后，需要对离子进行弛豫，即需要知道离子弛豫的方向和大小。【练习与思考】3-1.查阅有关文献和籍，找出3种以上不同的经验赝势，试分析他们的适用对象及特点？3-2.查阅MartinRM《ElectronicStructure》一书，写出模守恒赝势方法中模守恒的条件。3-3.查阅谢希德《固体能带理论》和有关书籍，试说明赝势平面波法采用了那些近似处理？3.2数值处理方法与技巧3.2.1超原胞方法对于固体体系，目前主要用三种模型来模拟材料特性，一是团簇模型(ClusterModel)，二是嵌入模型(EmbeddedClusterModel)，三是层状模型(SlabModel)。在以往的量子化学计算中，研究人员往往使用团簇模型来模拟固体材料，这种简化在化学上的确存在一定依据，因为根据化学家的直觉，一个分子体系作用是受局部相互作用支配的，从这一观点出发，可以用分子与原子簇的作用来反映分子与固体相互作用的性质。由于计算方法和计算条件的限制，没有考虑表面结构弛豫的影响，从而使计算的模型体系和实验体系存在较大的误差。第二种嵌入模型的提出主要是为了克服团簇模型在模拟材料表面时存在边界性问题，但该方法在计算过程中涉及到大量近似，需要针对不同的体系使用不同的计算方法，比如由KorringaJ、KohnW和RostokerN提出的格林函数(GreenFunction)方法和戴逊方程(DysonEquation)，所以该模型并不为广泛接受。由此，人们提出了超原胞模型，该模型在模拟体系时采用了周期性边界条件，特别适合研究金属、半导体这类具有周期性的凝聚态体系。在后面的实例中基本上都是选取了超原胞模型。层状模型(SlabModel)中的超原胞模型将体系看作沿晶格矢周期性排列的体系，在计算中所研究的原子都放在超原胞中，原子坐标或者其对应的周期性位置可用下式表示：(3.17)其中为原子在超原胞中的坐标，为原子种类的序号，为多个同类原子之间的序号。为格矢。目前许多第一性原理计算软件采取了超原胞模型，来构造周期性结构，包括三维或低维周期性结构。对于特殊体系如掺杂、缺陷、表面等，采取多倍原胞进行平移扩展，以保证物理上相邻原胞中的原子或分子没有相互作用。例如研究表面的分子吸附可假设它们在一个“盒子”里面成为周期体系，层与层之间用足够厚度的真空层隔离以忽略在盒子间原子的相互作用，如图3-4所示；再如，在研究中使用超原胞，认为它是可以在三维方向无限拓展，超原胞是没有外形的限制，假如这个晶体具有高点群的对称性，则它也可以用来加速计算。锐钛矿型TiO2(101)的超原胞如图3-5所示。 图3-4研究表面的分子吸附原胞 图3-5锐钛矿型TiO2(101)表面原子层超原胞3.2.2自洽电子弛豫方法表示电子结构松弛有多种方法，其中密度混合方法是最有效的。该方法是使用在一个在固定位势之下将电子能量本征值的总和极小化，而不是将总能做自洽式的极小化。在步骤的最后新的电荷密度就会与初始的电荷密度来混合产生新的电荷密度再计算体系能量以重复叠代直到收敛为止。密度混合有以下几种形式：(1)线性混合，新的密度是上一次输入与输出密度的线性组合。(2)Kerker混合，可用下式表达：(3.18)其中为截断波矢，为混合幅度。(3)Pulay形式，是最有效的一种方式，在Pulay形式中新的输入密度上所有先前迭代步中密度的组合，它不仅与和有关，而且和所有的先前迭代步有关。(4)Broyden混合形式，它与Pulay形式有类似之处。将本征值之和来做最小化既可以使用共扼梯度方法也可以使用加权余量方法。电子波函数是以平面波基底来表示，并且展开系数逐渐会被变化以便达到最小的总能。此极小化可利用每个波函数取独立的最佳化的band-by-band的技术，或允许同时更新所有波函数的all-band方法来达成。此方式用了Payne等人所提出的预先调节式的共轭梯度技术。传统总能极小化方法可能在具有晶胞一个方向上拉长的金属系统中计算会不穏定。而这是在表面做超晶格计算的典型设置中无法避免的。密度泛函方法对于绝缘体跟金属的状况都一样收敛良好。密度混合方法对于中等大小的绝缘体系统甚至都还提供3倍快的加速。密度混合方式主要的优势是当处在金属系统时可在相当少的次数很可靠的收敛。3.2.3几何结构优化技术对于给定各原子位置、元素种类的研究体系，通过密度泛函理论自洽求解Kohn-Sham方程可以得到整个系统处于多电子基态时的总能。总能量对系统虚拟微位移的导数就是各原子的受力，即Hellmann-Feynman力。这为我们理论预言物质的结构提供了一种行之有效的方法。因为自然界稳定的结构应该具有最低的总能，只要根据原子受力来调整原子的位置，直到整个体系的总能达到最低，所有原子受力为零。当然在实际的计算过程中，给出希望达到而且有限的计算精度，可找到能量面的(全局)最小值，这时所对应的物质结构就是自然界最稳定的结构，该过程被称为几何结构优化(GeometryOptimization)。为了确保搜索能量面的最小值时能找到全局最小而不是局域最小，并提高整个搜索过程的效率，需要一些强有力的搜索算法以使原子最快地运动到最稳定结构的位置。能量最小化算法一般分为两大类：全局极小和局部极小。全局极小算法可以得到基态构型，如模拟退火和遗传算法；局部极小算法找寻的是亚稳态结构，最常用的方法有直接能量最小化、最陡下降法、共轭梯度法、Newton-Rapshon方法、阻尼动力学方法等等。这里主要介绍材料模拟常用的三种优化方法。1.最陡下降法(SteepestDescentMethod)最陡下降法沿着局部净受力方向行走，以进行能量极小。从初始点开始，沿着局部梯度的反方向，并通过在此方向上的一维极小化，移动到该方向的极小点，再从这个点，开始重复以上过程，直到达到所要求的精度。最速下降法在远离极小点效率很好，在接近极小时效率不高，而且沿梯度方向每前进一步将对接下来一步都引入一个正比于它梯度的误差，常常只在优化的最初几步使用这种方法。2.共轭梯度法(Conjugale-gradientMethod)共轭梯度方法克服了最陡下降法的困难。在此方法中，每相邻两次优化的起始点的仍是正交的，但优化方向由当前梯度结合前一次优化方向和梯度共同决定，和互为共轭：(3.19)其中是一个标量数，由前一次优化起点的梯度和优化终点的梯度(即当前时刻梯度)共同决定，不同的算法给出各自的确定公式：Flrtcher-Reeves算法：Polak-Ribiete算法：对于能量函数，可以按如下方式进行优化：若初始点在，令，沿方向运用一维极小化方法到达该方向的一极小点，则。由和可得到，则，再沿找极小，重复以上过程，如果函数是含个变量的二次型，则通过次一维极小化就可以找到极小。上面两种方法在优化中只用了势能函数的一阶导数，即梯度。3.Newton-Rapshon法任何给定点的能量都可以展成泰勒展式：(3.20)是在x处的一阶导数矢量，为二阶导数矩阵，称为Hessian矩阵。对泰勒展式只取到二阶导数，而忽略高阶的，则位移矢量为：(3.21)其中一般情况下，可以重复以上过程直到能量最小，这就是所谓的Newton-Raphson方法。Newton-Raphson方法收敛速度很快，但是不能保证收敛的方向，而且需要计算Hessian矩阵，计算量很大。为了避免函数值上升，采用以下修正公式：(3.22)其中通过一维搜索算法确定，保证函数值向最小方向收敛。为了避免频繁计算Hessian矩阵，常采用更新校正方案(updatingscheme)，常见算法有：Davidon—Fletcher—Powell(DFP)和Broyden—Fletcher—Goldfarb—Shanno(BFGS)。具体表达式如下：(3.23)(3.24)其中。3.2.4快速傅立叶变换(FFT)平面波函数中使用的赝势常常是非局域化的。这意味着，为了要模拟全电子原子的散射性质，必须要在每个角动量信道里采用不同的散射位势。然而这些位势只在原子核心区域之内有所不同，采用实空间的表示方法是有效的。目前一般的计算软件采用的是选择在倒易空间中使用赝势，也就是说在倒易空间中进行对波函数和势操作的加总。这对于以平面波基底系数来表达的波函数是一种很常用的做法。然而对于大晶胞的Hamiltonian波函数和势能项进行实空间计算，会变得比在倒易空间中求值更快。这是因为赝势的非局域部份只有在核心内的区域才不是零，并且在远比各个单一原子核心区域还要大的晶胞中几乎到处都是零。因此，在实空间求值波函数和势的积，并且使用傅立叶转换来获得在倒易空间的值就变的更为有效率。这就意味着波函数和各种势在计算中需要一个有效的途径完成从实空间和倒易空间的相互转换，傅立叶变换为此提供了有力的工具。平面波赝势方法中采用一组有限个正交的平面波作为基矢，这就将久期方程中的矩阵元计算：，归结为傅立叶变换。快速傅立叶变换在计算中频繁使用。傅立叶变换网格不仅影响计算精度，也影响计算速度，因此也存在计算结果对傅立叶变换网格的依赖关系，严格地说，这也是一个收敛问题。傅立叶变换网格空间稀疏，将导致体系不易收敛或收敛结果粗糙，网格越密，则将大量消耗计算资源。3.2.5k空间取样规则在通常的能带图中，经常会出现比如Γ、Λ、Κ、Χ……等等的符号。这些符号表示的是布里渊区内具有高对称性的一些特殊的k点，这些k点有特殊重要的意义。在进行体系总能计算时，通常要对布里渊区内的波函数或本征值进行积分，在实际计算过程中，积分是通过对部分特殊选取的k点求和完成的。比较常见的k点网格方法有Monkhorst-Pack方法和四面体网格。根据Bloch定理，周期体系中的电子波函数可以表述为调幅平面波的形式，即。其本征能量和本征矢量为和。不同的电子状态按照量子数进行分类，而量子数n则表征能态的分立性。研究多体体系的价电子问题，归根结底是计算出不同类的电子状态的本征值和本征矢量，体系处于基态情况下，哪些不同的低能量状态被电子占据。因体系具有周期性，所以，第一布里渊区的所有可以代表所有的。但是，由于周期边界条件确定的k有无穷多个，而且计及相互作用势的实际体系中许可的k点在倒易空间内是不均匀的，实际计算过程中只能选取有限个点。第四章还将对k点的选取方法进行深入分析。在赝势平面波计算工作中，有限个k点在第一布里渊区内等权重均匀选取，这种选取k点的方法称之为Monkhorst-Pack方法。实际操作中考虑体系的对称性，将第一布里渊区依据点对称性划分为几个等价的“不可约空间”，自洽计算只在这个较小的不可约空间内进行。换言之，将研究那些互不等价的k量子数的集合，然后再用以描述整个电子体系的状态。如果体系的第一布里渊区的不可约空间大小为布里渊区体积的1/n，则总体性质取不可约体积内的计算结果n倍即可。但要注意的是，由于不可约布里渊区之间相交，第一布里渊区和第二布里渊区之间也有相交点，所以总有一些点为两个或几个不可约空间共有或为相邻布里渊区共有，这时如果进行占据态总能量和其它物理性质计算时采用简单倍乘就将导致完全错误的结果。关于k点网格选取方法的进一步讨论，请参考第4章相关内容。3.2.6基于密度泛函理论的第一性原理计算框架及步骤结合第2和3章相关计算理论和近似方法，本节总结出第一性原理计算框架中的关键步骤及典型实现方法，以及对应的波函数变化形式，如表3-1所示。表3-1基于密度泛函理论的第一性原理计算框架及步骤 关键步骤 典型方法 波函数形式 理论基础 Hohenberg-Kohn 唯一性，变分原理 多粒子波函数 Kohn-Sham 平均场近似 单粒子波函数 交换关联能 LDA，GGA 周期性 原胞(超胞)选择 Bloch定理 原胞周期函数 倒空间K网格 Monkhorst-Pack Smearing Fermi-DiracGaussMethfessel-Paxton 数值计算 基矢展开 PWAPW，MTOAO 展开基矢 电子、离子作用 赝势muffin-tin势紧束缚 实空间展开格点 FFT 离散数组 效率与精度 对角化哈密顿量优化电子波函数 steepestdescentdavidsonconjugate-gradient 电子密度自洽优化电荷密度 线性混合Kerker混合Pulay混合Broyden混合 几何优化优化原子位置 steepestdescentBFGSDIIS 收敛参数 电荷密度体系总能原子受力 【练习与思考】3-4.快速傅立叶变换(FFT)对体系哈密顿量进行数值化计算是如何实现的？3-5.写出超晶胞的周期性边界条件。并指出采用周期性结构计算的有优点。3.3基于MaterialsStudio的实践3.3.1软件介绍及计算步骤MaterialsStudio是Accelrys公司专为材料科学领域开发的可运行于PC机上的新一代材料计算软件，可帮助研究人员解决当今化学及材料研究中的许多重要问题。MaterialsStudio软件采用Client/Server结构(图3-6所示)，客户端可以是Windows2000、XP或NT系统，计算服务器可以是本机的Windows2000、XP或NT，也可以是网络上的WindowsXP、WindowsNT、Linux或UNIX系统。因此使得任何的材料研究人员可以轻易获得与世界一流研究机构相一致的材料模拟能力。图3-6MaterialsStudio软件采用Client/Server结构MaterialsStudio系列软件有一个核心模块，称为Visualizer。它是一个搭建分子、晶体及高分子结构模型的图形界面，同时提供各个计算模块的基本环境和分析工具，如图3-7所示。在Visualizer中提供了一个完整且容易使用的接口来操作CASTEP等计算程序，使得Visualizer所提供的交互式图形接口与仿真环境有效发挥CASTEP等程序的量子力学功能。下面TiO2表面模型构建都是在该模块下面完成。图3-7Visualizer图形界面MaterialsStudio软件中计算的模块很多，这里以MaterialsStudio系列软件中4.0版的模块CASTEP和Dmol3为例进行介绍。CASTEP模块是Materialsstudio系列软件中运用赝势平面波方法的模块之一。在Castep中，模型的建立基于超晶胞结构(supercell)的方法，并且使用周期性边界条件。由于周期边界条件的使用，晶体的表面用一个有限长度的模型(slab)就可以建立起来。周期性边界条件和Bloch原理联系在一起，即在一个周期性结构中，每个电子波函数都可以表示成和单胞中的波函数类似的形式，，而，这样每个电子波函数都可以写成平面波叠加的形式，即。有关更为详细的说明请看3.3.3节中的操作方法。Dmol3使用原子中心网格的数值函数作为其原子基，原子的基函数通过解不同原子的DFT方程得到，并将其储存为一系列的三次样条函数。这种基组是十分精确的，高精度的基组减少了重叠效应，于是体系可以得到准确的描述，能表示正确的电荷分布，因而可以提高对分子极化描述的精度。在Dmol3中，电子密度依据以原子为中心的多极部分密度展开，这提供了一种简洁而精确的表示密度的方法，通过求解Posson方程，可以用电子密度的多极表示来估算库仑势，从而将库仑势计算这一较快的步骤所代替，这一操作使计算所用的时间与体系的大小成线性正比。Hamilton矩阵元素用复杂的数值积分算法来计算，其计算量也接近与体系大小呈线性相关。此外，Dmol3的计算法允许对数值积分过程进行高效的并行处理。无论是对于分子体系或周期性结构，几何结构和过渡态的优化都采用非定域的内坐标。这包括在进行内坐标优化时，对其施加笛卡尔几何限制的能力。Dmol3采用了一种新的过渡态搜索方法，将传统的LST/QST算法与共轭梯度方法(conjugategradientmethods)结合使用。这种新的快速的方法允许在进行过渡态优化时不计第二阶导数矩阵，而这是传统过渡态搜索算法所必需的。计算步骤可以概括为三步：首先建立目标物质的几何结构；其次对建立的结构进行优化，这包括体系电子能量的最小化和几何结构稳定化；最后是计算要求的性质，如电子密度分布(Electrondensitydistribution)，能带结构(Bandstructure)、状态密度分布(Densityofstates)、声子能谱(Phononspectrum)、声子状态密度分布(DOSofphonon)，轨道群分布(Orbitalpopulations)以及光学性质(Opticalproperties)等。第一性原理计算的一般步骤如图3-8所示，下面将依照该步骤分别说明其细节及相关选项设置。图3-8MaterialsStudio软件计算步骤方框图3.3.2几何结构的构建本节以锐钛矿型TiO2的解理表面(101)面的构建为例，来说明构建材料的计算模型结构的具体方法。1.构建体相晶体结构采用MaterialsVisualizer模块来建立晶体表面结构。先建立材料的体相晶体结构，在此基础上沿某一特定方向进行剪切，从而获得所需要的理想表面。获得体相晶体结构有两种途径：一是直接从软件的结构库(structure)中导入。具体操作：点击File|import…命令，显示出ImportDocument(导入文档)对话框，注意，这个对话框也可以通过使用标准工具条中的导入按钮打开，在导入文档对话框中选择structure文件夹中需要的材料的晶体结构。在该文件夹包含了常见的半导体、陶瓷、金属、矿物以及聚合物等结构。另一种方式是已知晶体结构数据，如空间群，内坐标等，通过MaterialsVisualizer模块中的建模工具建立晶体结构，即创建3D文档。具体操作：选择FileNew…命令，在NewDocument对话框中选择3DAtomistic，在工程管理器中显示创建了一个3DAtomisticDocument.xsd的文件。在文件名上右击，选择改名命令，可以键入新的文件名。选择File|Save命令保存文件。选择BuildCrystalsBuildCrystal…命令(如图3-9)，打开建立晶体(BuildCrystal)对话框(如图3-10)。图3-9Build菜单图3-10BuildCrystal对话框在空间群(SpaceGroup)选项卡，输入空间群(Entergroup)文本框内输入空间群，也可以在下拉列表中选择空间群或者直接输入空间群序号。在晶格参数(LatticeParameters)中输入晶体的a，b，c轴的晶格参数。例如，锐钛矿型TiO2为四方晶格，空间群为I41/amd，晶格参数a=3.782Å，c=9.502Å。一旦空间群输入后，其α，β，γ的晶格参数自动设置。按Build按钮，一个定义了晶格参数的空晶格就出现在三维窗口中。此后，由于对称性已经产生，只需要在相应位置加入原子，对称位置的原子自动产生。选择“添加原子(BuildAddAtoms)”命令，打开“增加原子(AddAtoms)”对话框(如图3-11)。也可以通过单击“增加原子(AddAtoms)”按钮打开这个对话框。在“选项(Options)”部分，不要选择“原子成键测试(Testforbondsasatomsarecreated)”。在此选项选定时，MaterialsStudio会在建立晶体过程中自动加入键。MaterialsStudio也有一个灵活的“键计算(BondCalculation)”工具，可以选择、编辑、定义键参数。一般使用自动选项就行了。仍然是在“选项(Options)”部分，把“坐标系统(Coordinatesystem)”设为”分数(Fractional)”。回到“原子(Atoms)”选项页，元素列表中选择相应的Ti，键入以下a，b和c的值a=0，b=0，c=0。按添加(Add)按钮，一个相应的原子和它的对称位置原子就加入晶胞中了。在“原子Atoms”选项页，选择O元素，键入下列数值a=0，b=0，c=0.208。单击对话框底部“添加(Add)”按钮，一个氧原子和它对称位置的原子就被加入了。关闭“添加原子(AddAtoms)”对话框。为了提高结构的可视化，可在DisplayStyle对话框中的原子(Atom)选项中将“棍(Stick)”改为“球棍(Ballandstick)”选项(如图3-12)。最简单的操作过程为：右键视图的空白区域，从下拉菜单中选择DisplayStyle即可。图3-11Addatom对话框图3-12DisplayStyle对话框2.构建晶体表面构建晶体表面主要由以下四个步骤组成：步骤一：获得体相晶体结构。前面已经进行了具体介绍。步骤二：在上面构建的体相晶体结构的基础上，选择BuildCleavesurface，然后在Cleaveplane对话框里设定需要解理的晶面指数，如图3-13所示。例如锐钛矿型TiO2的解理表面是(101)面。步骤三：选择layer的厚度，就在Depth里选择Fraction，也就是重复的晶面数，注意这里不是周期数。在Cleave之后，如果需要修改这个Depth值，可以根据厚度来Recleave。Top的设置也很重要，它直接决定了剪切的位置，即暴露在表面的是那类原子。这个Top值究竟设置多少需要研究人员十分谨慎的选择，否则无法获得所需要的真正的实际表面。步骤四：同样在BuildCrystal选项里选择建立真空层(BuildVacuumSlab)对话框。如图图3-14，在该对话框中设定真空厚度。在这里，可以设定的厚度不能太小。一般情况下最好大于4Å。点击建立(Build)按钮，修改一下显示风格(Displaystyle)，获得图3-15显示的结构。 图3-13CleaveSurface对话框 图3-14BuildVacuumSlabCrystal对话框 图3-15锐钛矿型TiO2(101)面3.数据的导出事实上前面的晶体结构坐标数据可以方便导出，用于其它软件，如Wien2K，VASP等的输入文件。利用FileExport菜单命令将晶体结构数据导成cif或者pdb文件。这里以pdb文件导出为例举例说明。以下是导出的锐钛矿型TiO2(101)数据。从该文件中我们可以看到晶体表面的晶格常数以及原子坐标数据。|-----------------------------------------------------------------------------|REMARKMaterialsStudioPDBfileREMARKCreated:SatJun0906:45:50中国标准时间2008CRYST15.4433.77614.35390.0090.00110.30P1ATOM1TIMOL2-0.000-0.000-0.0001.000.00Ti...ATOM6TIMOL21.888-0.000-2.3721.000.00TiATOM7OMOL2-5.6641.888-2.7701.000.00O...ATOM18OMOL2-0.000-0.000-1.9731.000.00OTER|-----------------------------------------------------------------------------|3.3.3计算操作方法模型构建好之后，然后启动CASTEP程序进行几何优化和性质的计算，1.计算任务的设置在CASTEP软件中执行任务设置，主要是通过Visualizer应用窗口中的工具条之一“计算(Calculation)”来进行。可以更改工具框中的相应选项，来配置诸如：“电子选项(ElectronicOptions)”、“结构优化选项(GeometryOptimization)”、“属性选项(Properties)”等。首先建立周期性晶体结构(或者导入一个结构文件)。关于如何建立晶体结构上面已经进行了详细的说明。打开其中CASTEP中的Calculation对话框，在Setup中选择要计算的任务，单点能计算、结构优化、弹性常数、动力学过渡态等计算的设置。其次在“电子选项”进行精度设置。在CASTEP中在程序运行之前，从研究的具体问题出发，将软件中关键的一些任务参数设置成符合计算需要的值，我们才能得到所期望的运算结果。下面结合应用实例具体展开论述。(1)设置Setup选项选择需要计算的任务，如单点能计算、结构优化、动力学过渡态等，然后是精度设置和交换-关联函数的设置和电子自旋设置。·精度设置如表3-2所示，主要分为差、中等、好、超好四个等级。在涉及SCF收敛精度、K点取样精度、截断能等的设置时都要进行适当的选取。对材料的表面性质研究中，由于计算的精度要求较高，在Setup中Quality选项中，单点能和几何优化中最好选择Fine或Ultra-Fine，这个设置是对后面各种精度的总设置，高的精度设置直接导致完成计算收敛需要更多的时间。如果在电子选项中重新进行精度设置，这个对话框自动变成了Customed。表3-2电子选项中的精度设置对照表 Value Coarse Medium Fine Ultra-Fine SCFtolerance(eV/atom) 1.0e-5 2.0e-6 1.0e-6 5.0e-7 k-pointsseparation(1/Å) 0.07 0.05 0.04 0.04 Energycutoff Coarse Medium Fine Ultra-Fine·交换-关联函数的设置CASTEP提供了两种交换-关联函数的设置，一种是局域密度近似(LDA)，它使用的是CA-PZ形式的赝势，是引用CeperleyDM，AlderB和PerdewJP，Zunger采用的赝势形式。另外一种是广义梯度近似(GGA)，它有三种形式可供选择，分别是PBE(Perdew-Burke-Ernzerhoffunctional)、RPBE(RevisedPerdew-Burke-Ernzerhoffunctional)、PW91(Perdew-Wanggeneralized-gradientapproximation)，对应的是三种不同的广义梯度近似产生形式。以上两种交换-关联函数LDA和GGA及其对应的可选形式是通过CASTEP“Calculation”中的“Setup”选项来进行统一配置，如图3-16和3-17所示，改变其中“Functional”的类型和每一类型对应的函数产生形式，即可得到相应的设置。交换关联函数设置对计算结果的影响后面将加以测试说明。图3-16LDA设置对话框图3-17GGA设置对话框在范例中，在计算六种表面原子终止结构的能量时分别采用了局域密度近似(LDA)中CA-PZ方案和广义梯度近似(GGA)中PW-91方案，在研究表面原子弛豫和电子结构时采用了GGA-PW91方案。·电子自旋选项设置该设置一般是在研究电子自旋的磁性材料中才考虑，如果是不考虑材料的磁学性能，可以不做选择(默认设置是不考虑电子自旋)。如果要对某些原子的自旋属性进行设置，有两种方法。方法一：选中要设置的原子，在Modify的下拉菜单中选择“ElectronicConfiguration”，在出现的对话框中可以修改原子的“FormalCharge”和“Formalspin”等。HubbardU选项中，可以对LDA+U中的“U”进行设置，这仅仅在考虑电子自旋情况的时候有效。方法二：选中要设置的原子，在“属性窗口”(Properties)中直接选择“FormalCharge”、“Formalspin”和“Formalspindirection”等，双击鼠标后，可以在对话框中设置与电子自旋有关的相应参数。(2)设置电子选项在利用CASTEP做有关单点能量、结构优化、动力学、弹性常数、过渡态等计算时，必须对电子选项进行设置。在电子选项中主要有以下方面的设置。·截断能的设置CASTEP中分子轨道是通过平面波基来扩展。平面波基的数目是通过截断能的高低来控制，选择的截断能过低会影响计算结果的正确性，而选择的截断能过高会使计算量过大。因此，在计算中要选取合适的截断能。一般情况下，在计算前可以选取几个截断能来测试，哪个更合适，在保证计算精度的前提下选择尽可能低的截断能。因此，截断能是CASTEP计算中最重要的参数之一。设置截断能的最简单方式是在“CASTEPCalculationElectronic”选项中选择“More”，然后在对话框中选择“基本”(Basic)选项，在“能量截断”(Energycut-off)中填入设置的值(图3-18)。计算中截断能(Ecut)取为380eV。图3-18截断能的设置对话框·赝势的设置在CASTEP中有Ultrasoft赝势和Norm-conserving赝势这两种赝势，它们有各自的优势，Norm-conserving赝势一般适用于金属体系，对于纳米碳管的研究，采用低截断能的快速计算，使用Ultrasoft赝势更适合。设置的方法是在“CASTEPCalculationElectronic”选项中，在“Pseudopotential”的下拉框来进行选择。如图3-19所示。在描述离子实与价电子之间的相互作用时采用超软赝势(ultrasoftpseudopotentials，USPP)来，尽量减少平面波基的个数。为了减少计算需要的时间，采用了超软赝势研究锐钛矿型TiO2(101)表面结构。图3-19赝势的设置对话框·K点的设置布里渊区的设置是通过K点的设置来确定的。在K点的设置中，使用的是Monkhorst-Pack法表示。在倒格矢空间的划分适当的选择k点对于达成精确度与效率的平衡是很重要的。预设的Monkhorst-Pack点是在给绝缘体的0.1e-1到给金属的0.05e-1之间，这是因为金属系统需要更好的取样。如此通常就能产生足够的点数。例如传统硅晶胞所需的222，应该检查增加一个，直到建议出来的奇数值Monkhorst-Pack参数能更为有利。我们必须推荐用k点取样的增加来减低有限基底集的修正并促使在一个固定能量下晶体松弛更加精确。设置方法是在“CASTEPCalculationElectronic”选项中选择“更多”(More)，然后在对话框中选择“k-points”。在“k-points”的选项中又有几种方式。第一种是“只取Γ点”(Gammapointonly)，这对于计算体系具有较大的原子数目并且对称性低的情况下可以考虑。第二种是按“精度”(Quality)来选择，它有三个等级，分别是“粗糙”(Course)、“中等”(Medium)、“精细”(Fine)。这三个等级对应着不同的Monkhorst-Pack点。比如在做结构优化任务时一般精度都要选择“精细”(Fine)这一等级。第三种是可以给定k点的间隔，这样也就定下了在布里渊区中k点的设置。第四种是直接给定沿着超原胞倒格矢空间三个基矢a、b、c的k点取值，如在图3-20中显示的666。图3-20K点的设置对话框·自洽精度设置自洽循环计算中上次计算总能的结果与下次计算的结果的差值小于选择的自洽精度设置时，自洽计算终止，完成计算任务。2.结构优化任务的设置结构优化是CASTEP计算中经常要进行的计算任务，特别是计算所关注体系的各种性质时，必须首先进行结构优化的计算，在得到结构优化结果文件以后，才能进行性质的计算。所以，正确设置结构优化参数是非常重要的。在CASTEP软件中，有四个收敛参数来控制结构优化的精度。如表3-3所示。表3-3结构优化任务中的精度 Value Coarse Medium Fine Ultra-Fine Energytolerance(eV/atom) 5.0e-5 2.0e-5 1.0e-5 5.0e-6 Max.forcetolerance(eV/Å) 0.1 0.05 0.03 0.01 Max.stresstolerance(GPA) 0.2 0.1 0.05 0.02 Max.displacementtolerance(Å) 5.0e-3 2.0e-3 1.0e-3 5.0e-4第一个是能量收敛精度(Eergytolerance)，单位为eV/atom，是体系中每个原子的能量值；第二个是作用在每个原子上的最大力收敛精度(Max.forcetolerance)，单位为eV/Å。第三个是最大应变收敛精度(Max.stresstolerance)，单位为GPa；第四个是最大位移收敛精度(Max.displacementtolerance)，单位为Å。这些收敛精度指的是两次迭代求解之间的差，只有当某次计算的值与上一次计算的值相比小于设置的值时，计算才停止。设置方法是在“计算”(CASTEPCalculation)中选择“装载”(Setup)选项，再在“任务”(Task)选项中选择“几何优化”(GeometryOptimization)。在“更多”(More)选项中可以进行收敛精度的设置(见图3-21)。这里自洽精度不能设置过高，否则出现无法收敛的情况。图3-21结构优化任务的设置对话框3.计算体系性质的设置在CASTEP中可以计算体系的性质，如能带结构、态密度、布居数分析、声子色散关系、声子态密度、光学性质、应力等。在这里介绍一下在涉及到的能带和态密度这两项的计算设置。在计算这两项性质之前，必需先进行自洽计算得到基态能量，而结构优化能够做到这一点，这就是之所以要在计算能带和态密度之前对体系进行结构优化的原因。图3-22能带结构计算设置对话框设置能带计算时如图3-22所示。在“计算”(CASTEPCalculation)中选择“属性”(Properties)选项，然后选中“能带结构”(Bandstructure)，在下面的“空带”(Emptybands)中设置好研究体系导带中选取空带数目。同时在“k点设置”(k-pointset)项里选择计算能带时k点的设置，它也分为“粗糙”(Course)、“中等”(Medium)、“精细”(Fine)三种。点开“更多”(More)，即可看见详细的沿各高对称路径的k点设置。4.0版本中还允许能带计算中的交换-关联函数与自洽计算中的不同，所以在CASTEP的PropertiesBandStrutureMore中设置了交换-关联函数的选取项。能带中能量的收敛精度在“能带收敛精度”(Bandenergytolerance)进行设置。态密度的设置与能带类似，它在“属性”(Properties)选项下有部分态密度(PDOS)的选择项，其k点的设置和在“电子选项”中介绍的完全相同，并且也可以进行与自洽计算中的不同的交换-关联函数，如图3-23所示。图3-23态密度设置对话框4.计算结果的分析操作计算结果从计算服务器上返回以后，在Visualizer界面中进行分析。打开“分析”(CASTEPAnalysis)对话框，对于单机计算的结果可以直接看到。计算时把计算模型取名为“锐钛矿型TiO2”，那么当锐钛矿型TiO2的能带计算完成后，在文件夹锐钛矿型TiO2下有3DAtomistic_BandStr.castep文件生成。首先在Visualizer界面中把该文件打开，接着点击Visualizer应用窗口中的工具条“Analysis”，“”就会有如图3-24的对话框出现。该对话框中的“Scissors”选项是剪刀工具，可以把能带作一个微调，使之显示的结果符合实验值。“Graphstyle”选项是选择图形显示的方式，分为“点”(Points)、“线”(Line)两种，一般选择“线”(Line)。在计算了能带以后，可以同时把总的态密度显示出来。所以可以看到，在3-24对话框中进行了一个作总的态密度图的设置。然后选中“查看”(View)按纽，则在Visualizer界面中会显示能带和对应的总的态密度图得到的能带和总的态密度图还可以导出到图形数据处理专业软件，如OriginLab软件中进行处理，以利于更直观的分析。3.3.4结果可靠性测试为了分析该软件的在材料方面的可靠性，对金红石型TiO2结构的单点能进行了计算。在计算前，对相关参数进行了设置，分别考虑平面波截断能Ecut，SCF以及K点设置等对计算结果的影响，并采用LDA和GGA以及Ultrasoft赝势和Norm-conserving赝势的方法对结果的影响进行了分析。使用的计算资源配置如下：CPU(IntelP43.0GHz)、内存(kingkong，1G)，硬盘(80G，串口)。TiO2结构为四方晶系P42/mnm空间群，对称性D4h，采用CromerDT和HerringtonK在1955年测得的实验数据。晶胞参数a=b=0.45936nm，c=0.29587nm，Z=2，内坐标u=0.305，其结构如图3-25所示。在设置参数时采用保持其中一个参数变化而其它参数不变，依次判断各参数对计算结果的影响。为了计算的可靠性，各参量的初始设置均较高，分别为：平面波截断能Ecut=380eV。SCF设置为“极精细”(Ultra-fine)，即5×10-7；k点取为6×6×6；FFTgrid设置为“精细”(Fine)。对SCF的收敛性进行测试，发现收敛性设置对总能的影响很小，总能差值在10-5以下。且发现SCF设置越高，获的总能越精确，但是将花费计算机的大量时间。一般认为在体相结构，例如掺杂，缺陷等，SCF设置为中等质量就可以了。在我们后面对锐钛矿TiO2(101)面的几何结构和表面研究中，SCF设置应该设置得更高。否则将影响几何结构和电子结构的准确预测。如图3-26所示，图中的结果是在部分参量变化，而其它参数保持上面的设置不变。 (a)超软赝势下截断能与总能量的关系 (b)超软赝势下截断能与计算时间的关系 (c)模守恒赝势下截断能与总能量的关系 (d)模守恒赝势下截断能与计算时间的关系图3-26平面波截断能Ecut对单点能的影响 n×n×n(a)超软赝势下k点设置与总能量的关系 n×n×n(b)超软赝势下k点设置与计算时间的关系 n×n×n(c)模守恒赝势下k点设置与总能量的关系 n×n×n(d)模守恒赝势下k点设置与计算时间的关系图3-27k点设置对单点能的影响正确设置优化的参数，是结构优化非常重要的步骤。对于精度的设置，可以根据前面的表3-3中的四个参数来控制结构优化的收敛参数，也可以进行任意设置四个不同参数的。设置的原则是：既要达到好的收敛精度，获得可靠的原胞晶体结构，又要尽量减少优化所花费的时间。为了讨论TiO2在高压下的晶型转变，结合对金红石型TiO2不同静水压下的测试，分别采用了GGA和LDA两种方法进行几何优化，其结果如图3-28所示。在同样的晶胞体积下，用GGA的方法计算的单点能比LDA方法计算的能量要高12eV/cell。事实上在计算其他结构也呈现同样的现象。从图3-28(a)可以看出，用LDA方法计算的体积比实验值要小2.6%，而用GGA方法计算的能量最低点的单胞体积偏大1.9%。在图3-28(b)中，用GGA方法计算的晶格常数a，c与实验的误差在0.94%之内，c/a的比值误差在0.83%。图3-28(c)是可以看出，用GGA方法明显要优于LDA方法。综上研究结果分析表明，在研究表面弛豫或结构时，当精度设置较低，可以很快收敛，能量达到最低，但是获得的弛豫和表面结构方面的信息很不可靠。因此在诸如此类研究中应该适当提高精度设置。结构优化精度设置是个技术问题，考虑到精度与时间之间的矛盾，很多时候无法完成优化。一般的体相结构或少量的掺杂，在结构优化方面的要求比较低，而在研究弛豫和表面结构方面却要求很高。关键是需要提高到多少为佳。关于这方面的问题的深入讨论可以参见后面的CASTEP的实际应用。 (a)能量与体积的关系 (b)能量与c/a的关系 (c)体积与c/a的关系(虚线交点代表实验值)实验值来自文献Shang-DiMoandChingWY.Phys.Rev.B199551(19)13023～13031图3-28GGA和LDA方法两种计算的结构常数结果的对比【练习与思考】3-6.利用MaterialsStudio或者CrystalStudio重复表面模型的建立过程，并建立NaCl晶体结构和bccFe(001)表面。3-7.简述一个完整的计算流程，结合计算成本，你认为其中哪些参数的设置对计算结果的精度以及计算所需要的时间影响比较大？3-8.利用MaterialsStudio软件包中Visualizer模块，建立聚合物乙烯结构，利用Dmol3模块计算材料的电子结构。3-9.建立了硅的一个(310)的孪晶界(twinboundary)，主要步骤有：晶面切割，选择合适的切片厚度，按照建立表面模型的方法建立(310)和(-310)孪晶界结构，并且采用分子力场进行优化。3.4锐钛矿型TiO2(101)表面分析计算锐钛矿型TiO2是一种重要的宽禁带半导体材料，具有独特的物理化学性质。尤其是具有大比表面积的纳米锐钛矿型TiO2材料的低指数表面备受人们关注。其中(101)面是纳米锐钛矿型TiO2微晶材料中主要存在的表面，该表面具有表面能低、结构稳定等特点而备受重视。在实验上，因很难获得符合实验要求的锐钛矿型TiO2清洁表面而使研究变得困难。然而运用可靠的理论计算研究晶体表面，则可以提供实验分析技术所不能提供的数据，并指导TiO2的实验研究。为了进一步理解第一性原理在实践中的具体运用，这里选取TiO2材料为研究对象，运用MaterialsStudio软件中的CASTEP模块对该材料表面原子几何及弛豫结构和电子结构的开展研究。对于计算结果进行计算与分析，首先要明确研究的是材料的电学、磁学还是光学性质，或者是材料化学吸附、催化、化学反应等问题。要有针对性的设置计算任务和结果分析。根据前面的参数设置，对锐钛矿型TiO2(101)面计算的结果进行分析。重点研究表面结构建立、表面弛豫、电子布居分析及表面的电子结构等方面。3.4.1表面结构的建立锐钛矿型TiO2(I41/amd)属四方晶系，每个晶胞包含四个TiO2单元。对建立的晶胞进行优化后，晶胞参数分别为a=b=0.3784nm；c=0.9712nm；V=0.1391nm3，与Burdett等用中子衍射实验在15K的低温下获得的值相比，误差分别仅为+0.057%，+2.214%，+2.324%。对该结构在垂直[101]方向进行剪切，由于断键情况和原子位置的不同，存在六种不同的表面原子终止结构。即A类结构：最表层终止于Ti原子，为二配位，次层是O原子，为二配位；B类结构：最表层终止于O原子，为二配位，次层是Ti原子，为四配位，第三层是O原子，为二配位；C类结构：最表层终止于O原子，为一配位，次层是O原子，为二配位，第三层是Ti原子，为五配位，第四层是O原子，为二配位；D类结构：最表层终止于Ti原子，为四配位，次层是O原子，为三配位，第三层是O原子，为三配位，第四层是Ti原子，为五配位；E类结构：最表层是O原子，为二配位，次层是Ti原子，为五配位；F类结构：最表层是O原子，为一配位，次层是O原子，为二配位。其剪切位置如图3-29所示。根据上面六种表面原子终止结构，分别建立12个原子的晶体表面原胞，其化学式均为Ti4O8，将它们置入一个方形空腔中，选取真空厚度为0.8nm，初始晶胞参数a0=0.5545nm，b0=0.3784nm，r=109.95o，分别计算体系的总能量。哪种结构是最可能的(101)面需要从能量上比较验证。3.4.2参数设置在计算六种不同的表面原子终止结构的体系能量时分别采用了局域密度近似(LDA)中的CA-PZ方案和广义梯度近似(GGA)中的PW91方案。在研究表面原子弛豫时采用了PW91方案，从表3-4中可以看出，该方案计算出的TiO2几何结构参数比PBE和RPBE方案更接近实验值。计算中平面波截断能Ecut取为380eV。自洽场运算中，应用了Pulay密度混合法，自洽精度设为每个原子能量收敛至2.0×10-6eV，原子受力不超过0.5eV/nm。在对模型的结构优化中，采用了BFGS算法，每个原子能量收敛至2.0×10-5eV以内。第一布里渊区按5×5×3进行分格。表3-4TiO2几何结构优化后结果与实验值的比较 方案 a/nm c/nm c/a PBE 0.3788 0.9770 2.579 RPBE 0.3802 0.9910 2.607 PW91 0.3784 0.9712 2.566 实验值 0.3782 0.9502 2.512实验值来自文献BurdettJK,HughbanksT,MillerGJ,RichardsonJWandSmithJV1987Chem.Soc.10936393.4.3表面原子几何结构固体表面几个原子层中原子的排列情况，包括表面区各层中原子的种类、彼此间的相对位置(键长和键角等)、各层间的间距、单元网格相对体内网格的取向等。稳定的固体表面其总能量相对较低。前面六种结构模型体系总能量的计算结果如表3-5所示，计算精度都较好地达到了设定的收敛度。可以看出用LDA算法和GGA算法计算出的E种表面结构体系能量比其它五种结构体系能量分别低2.6—10.2eV和1.9—9.4eV。根据总能量最小原理，E种表面原子终止结构是TiO2(101)面最可能出现的表面，即锐钛矿型TiO2(101)最表层应终止于两配位的O原子，次层为五配位的Ti原子。从[010]方向看，TiO2(101)表面明显为波浪形的弓形几何结构，图3-5显示了2×2×3超胞表面。表3-5不同的表面原子终止结构的体系总能量和收敛度 模型 总能量(LDA) 收敛度 总能量(GGA) 收敛度 类型 /eV ×10-7(eV/atom) /eV ×10-7(eV/atom) A -9908.09119 2.75 -9940.39862 4.07 B -9915.74661 4.28 -9947.86295 4.99 C -9908.09329 6.54 -9940.40392 1.10 D -9911.62948 1.68 -9943.32798 3.25 E -9918.31973 1.37 -9949.80965 1.59 F -9911.62833 6.21 -9943.33177 7.593.4.4表面能计算晶体中原子排列的周期性在垂直于表面的方向上于表面处突然中断，存在大量的断键，形成具有一定表面能的结构。在建立的表面模型中，真空层厚度和原子层数改变对表面能和表面弛豫的计算结果有很大关系。为了获得较为准确的计算数据，首先考虑模型的表面能与真空层厚度和原子层数的关系，以确定合适的模型参数。表面能定义为：Esurf=(Eslab-NEbulk)/2A(3.25)式中N为表面单胞所含的TiO2分子数；Esurf为含有N个分子的表面所具有的表面能；Eslab为含有N个TiO2分子的表面单胞能量；Ebulk为体相中每个TiO2分子的能量；2A为单胞上下总表面积。图3-30显示了原子层数为12层时表面能与真空厚度的关系。从图3-31中可以看到，表面能与建立表面层数有很大关系。当原子层数达到18层时，表面能收敛于0.462J/m2，收敛度小于0.02J/m2，这与文献中用其它方法计算的表面能很接近，但是没有发现像金红石型TiO2(110)面出现的表面能与原子层数关系的那样的奇偶行为。图3-30原子层数为12层时表面能与图3-31真空厚度为8Å时表面能与真空厚度的关系原子层数的关系3.4.5表面原子弛豫由于体相的三维周期性在表面处突然和断，表面上原子的配位情况发生变化，相应地表面原子附近的电荷分布将有所改变，表面原子所处的力场与体相内原子也不相同。为使体系能量尽可能降低，表面上的原子常常会产生相对于正常位置的上、下位移，结果表面相中原于层的间距偏离体相内原子层的间距，产生压缩或膨胀。表面上原子的这种位移称为表面弛豫。表面弛豫往往不限于表面上第一层原子，还会波及到下面几层原子，但愈深入体相，弛豫效应愈弱。为了定量描述表面原子的几何结构弛豫，取X，Y，Z轴分别沿[01]，[010]，[101]晶向，表面原子层数为i，设ΔX=Xi-X0i；ΔY=Yi-Y0i；ΔZ=Zi-Z0i为原子偏离理想晶格位移，X0i，Y0i，Z0i为理想表面晶格第i层原子位置坐标分量，Xi，Yi，Zi为弛豫后表面第i层原子坐标分量。优化结构后表面原子发生的位移如表3-6所示(以表面以下第七原子层的O原子为坐标原点)。由表3-6可以看出TiO2(101)面有明显的结构弛豫现象，弛豫后体系的总能量降低了约1.9eV。对比表3-6中18层和24层模型计算的结果，可以看出建立模型原子层数对表面原子的位移有一定的影响。可以知道：当模型中原子层数取得越多，表面模型越接近实际的表面。从24层模型计算的结果来看，最明显向表面外弛豫的原子是O3c和Ti6c，弛豫距离分别为0.246Å和0.181Å，同时Ti5c原子向里移动0.155Å，在[01]方向上移动0.139Å，表面最外层的O2c有向里的微小移动，而在[01]反方向上移动0.102Å，使在[01]方向上的O2c和Ti5c相互靠近，距离减小，O2c—Ti5c键连作用增强。由于Ti5c原子向里移动，O3c原子向表面外移动，使原来位于Ti5c原子的下方O3c移动到Ti5c原子上方。表3-6锐钛矿型TiO2(101)表面原子位移理论计算结果与文献报道值的比较(正值表示弛豫向[01]，[010]，[101]方向，负值表示弛豫为反方向)。有关原子符号请参见图3-5。 Label Nlay 本文工作 Lazzeri[a] Nilsing[b] ΔX/Å ΔY/Å ΔZ/Å ΔX/Å ΔZ/Å ΔX/Å ΔZ/Å O2c 18 -0.094 -0.081 -0.016 -0.29 -0.02 -0.072 -0.040 24 -0.102 -0.072 -0.012 O3c 18 -0.004 -0.010 0.251 -0.16 0.19 -0.003 0.197 24 -0.010 -0.009 0.246 O3c′ 18 0.001 -0.031 -0.010 -0.17 0.06 0.016 0.032 24 -0.003 -0.024 0.009 O3c〞 18 -0.056 -0.062 -0.064 -0.15 -0.07 -0.054 -0.050 24 -0.058 -0.053 -0.053 Ti5c 18 0.144 0.009 -0.169 -0.02 -0.18 0.091 -0.245 24 0.139 0.017 -0.155 Ti6c 18 -0.059 -0.070 0.183 -0.17 0.20 -0.020 0.280 24 -0.051 -0.056 0.181 [a]LazzeriM,VittadiniA,SelloniA.Phys.Rev.B,2001,63:155409[b]NilsingM,LunellS,PerssonP,OjamaeL.Surf.Sci.,2005,582:493.4.6电子布居分析通过电子布居数的分析可以了解固体单胞中的电荷分布、转移和化学键性质。根据密里根(Mulliken)布居数的概念，分布于原子周围的键的电子电荷，称为原子布居数，分布在两原子之间的重叠电子电荷，称为键布居数，它表征了两原子之间形成化学键的强度以及成键特征，如共价键或是离子键。键布居数越大，表明两原子形成的化学键越强；反之，表示成键越弱。零值表明一个完美的离子键，小于零值表明离子键的程度，大于零值表明原子间共价键的程度。然而Mulliken方法也有一些严重的缺陷使得Mulliken布居分析在某些情况下得到的数值不合理，影响了其应用。为克服Mulliken布居的这些缺陷，有人提出将重叠区电荷在两原子间加权分配，Pollak等则提出将重叠电荷按Löwdin提出的保持重叠电荷偶极矩不变的原则分配，Davidson、Politzer和Harris，Pople等也提出过许多种其它方案。MS的Castep模块可以直接分析Mulliken布居。表3-7列出了理想表面和弛豫表面原子间的键长及键布居数。从表中可以看出，弛豫后表面几层原子的键长发生了较大的变化，其中Ti5c—O2c键长变化最为明显，收缩了约0.171Å，同时Ti5c—O3c〞键长也有较大的收缩，使其成为表面上最短的键长，仅有1.78245Å由于O2c离子半径比Ti离子半径大的多，而且弛豫后O2c离子中心位于TiO2(101)最表面处。因而锐钛矿型TiO2(101)暴露在最表面是未满配位的O原子，易在表面形成O空位。表3-7原子间的键长及键布居数 理想表面 弛豫表面 Ti—O键 键长/Å 布居数/eV 键长/Å 布居数/eV Ti5c—O2c 2.00071 0.40 1.82937 0.36 Ti5c—O3c 1.93974 0.69 1.96267 0.68 Ti5c—O3c′ 2.00071 0.27 2.06119 0.19 Ti5c—O3c〞 1.93974 0.43 1.78245 0.54 Ti6c—O2c 1.93974 0.47 1.83965 0.46 Ti6c—O3c〞 2.00071 0.19 2.07603 0.16表3-8Mulliken净电荷布居数(单位电量) 表面原子 O2c O3c O3c′ O3c〞 Ti5c Ti6c 理想表面 -0.64 -0.71 -0.69 -0.64 1.34 1.33 弛豫表面 -0.59 -0.73 -0.70 -0.64 1.37 1.29TiO2(101)表面原子发生弛豫后，其表面原子核外电荷分布发生改变。表3-8为理想表面和弛豫表面Mulliken布居分析后的结果。可以看出弛豫后Ti5c和O2c核外电荷电量分别减少了0.03、0.05，而Ti6c增加了0.04，其它原子核外电荷的转移相对较小。总体上表现出表面悬挂键上多余的核外电荷向表面内层原子转移。3.4.7表面原子弛豫对电子结构的影响电子态密度(DOS)描述了电子的能量分布，设在之间的能态数目为，则电子态密度就为(3.26)电子态密度表示了在能量空间中电子态的分布，可以给出局域分波轨道间相互作用(轨道杂化)以及电子态能级移动和能级弥散信息，可以直接用于分析掺杂效应等。通过分解的总态密度，获得的表面的O2c，O3c，O3c′，O3c〞，Ti5c和Ti6c原子分别对态密度贡献，如图3-32未弛豫表面的态密度在费米能级非常高，这主要是由于二配位的O2c原子的贡献。费米能级附近的态密度越高，其活性越强。弛豫后的表面能带结构并没有出现上述明显的现象。对于弛豫后的表面，二配位的O2c主要贡献价带的中部，弛豫后费米能级上的态密度要比理想表面的态密度要低很多，且禁带宽度增大，这表明弛豫后的锐钛矿型TiO2(101)表面趋于稳定。图3-32表面的O2c，O3c，O3c′，O3c〞，Ti5c和Ti6c原子分态密度虚线表示弛豫前的分态密度，实线表示弛豫后的分态密度【练习与思考】3-10.选取金红石型SnO2晶体，构建其(110)表面结构，采用CASTEP进行参数设置，要求计算能带结构和态密度，并对结果进行分析。3-11.对锐钛矿型TiO2进行非金属掺杂，如B，N，S元素，分析掺杂对TiO2能带结构的影响，分析其光学性质。3.5*锐钛矿型TiO2(101)表面缺陷结构的计算表面原子的周围环境不同于体相结构，表面上的原子受到的静电力发生改变，使表面原子发生弛豫或者重构，以降低表面能，表面结构趋于稳定。在前一节中主要讨论的是锐钛矿型TiO2(101)无缺陷的表面弛豫等几何结构问题。而实际的材料表面往往存在各种表面缺陷，如O空位或Ti空位等。早期对于TiO2表面空位的研究主要集中在金红石型TiO2(110)面，而对于锐钛矿的表面的缺陷研究相对较少。由于锐钛矿TiO2(101)表面原子与体相原子环境有很大差异，这将导致表面原子的电子结构与体相有很大不同。本章将系统研究锐钛矿型TiO2体相电子结构以及清洁表面(101)面的电子结构。3.5.1表面缺陷模型本节在前面锐钛矿型TiO2(101)表面模型的基础上考虑表面点缺陷，即采用Slab模型模拟半无穷大晶体(表面)。选取真空厚度为0.8nm，以保证上下两层表面之间的作用可以忽略，表面厚度选取18层原子，其完整表面能收敛度小于0.02J/m2.锐钛矿型TiO2(101)面最表层应终止于两配位的O原子，次层为五配位的Ti原子。考虑到实际TiO2表面存在的空位和间隙原子的浓度可能比较低，为此沿[010]方向建立1×2表面超胞，降低点缺陷的浓度，以减小相邻点缺陷间的相互作用。用该表面超胞模拟有缺陷的表面是足够的。本文模型中的点缺陷是通过在表面上移去或增加原子(Ti或O)的方式获得，分别考虑了O空位()，Ti空位()，O间隙(Oi)，Ti间隙(Tii)以及不同位置的缺陷结构，如图3-33所示。为了节约计算时间，采用了固定超原胞常数的方法，仅仅让表面上的12层原子层进行弛豫优化。SHAPE\*MERGEFORMAT为O原子为Ti原子图3-33锐钛矿型TiO2(101)面18层原子模型图，其中Tii1和Tii2为Ti间隙原子，Oi为O间隙原子3.5.2表面缺陷原子几何结构不同类型的表面点缺陷对周围原子间的作用是不同的，这些缺陷会影响周围原子位置及周围电子的分布。表3-9列出了表面空位缺陷导致周围原子的弛豫结果。对于空位缺陷，弛豫后的总能比未弛豫时要低0.57eV；空位周围的Ti1和Ti2原子向空位外移动的位移较大，向表面里分别弛豫了0.0142和0.0255nm，这主要源于空位上有效的正电荷与邻近的阳离子之间的静电排斥作用，或者说O1空位的产生使两个Ti原子之间失去了屏蔽它们之间排斥作用的原子。根据原子球内总电荷数和有效离子半径，即所谓的ShannonPrewitt半径，Ti4+离子和O2-离子的半径分别取为0.061和0.140nm，计算出缺陷周围原子的Mulliken净电荷布居数。从表3-9可以看到围绕着O1空位的两个Ti原子净布居数分别为1.23和1.15e。与产生缺陷前的布居数相比较，均增加了0.14e的电荷，其它Ti原子的电荷改变小于0.02e，这说明空位的产生导致了Ti离子周围局部电子密度的增加，Ti与O的键联作用有所增强。表3-9表面空位缺陷导致的周围原子的弛豫(沿[101]方向的位移，向上弛豫为正，向下弛豫为负；单位:nm)及空位周围各原子Mulliken净电荷布居数(单位电量) 表面类型 O1 O2 O3 O4 Ti1 Ti2 Perfect 0 -0.59 0 -0.73 0 -0.70 0 -0.64 0 1.37 0 1.29 — — -0.0009 -0.74 0.0227 -0.69 0.0140 -0.68 -0.0142 1.23 -0.0255 1.15 -0.1330 -0.08 0.0021 -0.65 0.0080 -0.63 -0.0021 -0.32 — — -0.0170 1.41与空位缺陷周围的Ti原子移动情况相似，空位缺陷上有效的负电荷与周围O原子间的静电排斥作用导致周围O原子向空位外移动，其中O2，O3原子分别向[010]方向移动了0.042和0.026nm，而O1原子由于失去Ti1的键连作用，使Ti2—O4达到0.2979nm，而O1向表面里弛豫了近0.1nm，O1有填补Ti空位的趋势.此时O1与O4的距离缩短至0.1334nm(氧气分子的键长0.1214nm)，这说明了还原气氛下失去Ti1原子，导致表面O原子有自发地团聚成氧气分子的趋势。考察这些原子的Mulliken电荷布居数可以发现，为了达到静电平衡，原子之间的电荷重新分布。O1、O4原子周围的电子显著减少，O1—O4之间的共价特征明显，而Ti2原子周围的电子减少，显然与它形成低配位的几何构形有关。而对空位结构进行优化，发现其结构同空位结构相似。其中O1原子向表面里有较大的弛豫，致使O1与O3的距离仅仅有0.1273nm。这些说明表面Ti原子空位缺陷的出现也会导致表面O原子团聚成氧气分子。对O间隙缺陷结构进行研究，发现O间隙原子与周围的O原子形成了双键(键长为0.1423nm)，即O间隙原子容易吸附在晶格上的O位上，并与Ti成键，形成Ti(O(O二聚物。特别值得关注的是间隙原子的出现对结构的影响比较小。结果显示了间隙原子与周围的最近邻的五个O原子之间的距离在0.1953～0.2138nm之间，主晶格结构没有显著变化，因而该缺陷形成能很低。然而间隙原子存在晶格间隙中，由于其周围最近邻原子均是Ti，其排斥作用相当大，导致了表面结构显著的扭曲，其总能比缺陷结构高1.35eV，因此该位置的Ti间隙缺陷结构是不稳定的。3.5.3点缺陷形成能与体相结构中点缺陷不同的是研究表面点缺陷不仅要考虑点缺陷类型，还要考虑点缺陷在表面上出现的位置。为此，分别计算了表面上各种可能出现的点缺陷类型和不同位置点缺陷的形成能。缺陷形成能定义为：(3.27)其中为有缺陷的表面总能，为化学计量比的表面总能，、分别为O和Ti原子的化学势，、分别为移去或增加原子产生空位和间隙缺陷的数目。例如，表示一个Ti原子间隙缺陷；，表示一个Ti空位缺陷。为了获得有意义的结果，对和的计算均采用了相同的计算参数。晶体生长时原子的化学势与生长环境中的气氛有关。因不能高于Ti金属原子的化学势；在还原性气氛下生长的TiO2满足上限。同样因不能高于氧气分子中原子的化学势；在氧化性气氛下生长的TiO2满足上限。由可知，的上限对应的下限，的上限对应的下限。为了得到一个O2(g)的总能量，将两个氧原子置入一个1nm3的立方真空晶胞中，优化后键长为0.1214nm，其单点能为-872.04eV，该结果与其它理论计算的结果误差小于0.2%。本计算的Ti单晶中原子的单点能为-1606.6eV。根据方程(3.27)和极端生长条件下的化学势，表3-10列出了计算出的各种点缺陷形成能。表3-10TiO2(101)面各种缺陷形成能(eV)(采用克勒格尔-温克的点缺陷符号) 缺陷类型 Oi 氧化性气氛 5.61 6.08 5.82 5.78 12.41 10.90 11.65 2.80 3.58 10.33 8.98 1.70 还原性气氛 0.86 1.33 1.07 1.02 2.91 1.40 2.14 12.30 13.08 0.83 -0.52 6.45图3-34TiO2晶体表面生长的和的关系，其中水平虚线对应Ti金属原子的化学势，竖直虚线对应O2(g)的化学势为了进一步研究点缺陷结构的稳定性，分别计算了上面讨论的几种点缺陷的形成焓。缺陷形成焓定义为：(3.28)其中为存在点缺陷的表面结构总能，、分别为表面结构中Ti和O原子个数.在还原气氛下的Tii缺陷形成焓为-9.49eV/TiO2(或3.16eV/atom)。图3-34中的和关系说明了不同缺陷结构下允许的生长条件，这里定义，。其中实线为化学计量比TiO2平衡生长时、满足的条件：，在此条件下TiO2晶体缓慢平衡生长。我们计算出化学计量比TiO2分子的形成焓为-10.3eV/TiO2，它表示为保持TiO2的热力学稳定性，要求满足-10.3eV。如果<-10.3eV，则表面结构出现分解，晶体无法生长，对应图3-34中实线下方区域。如果>-10.3eV，则晶体快速生长，但生长不均匀，对应图3-34中由水平虚线和竖直虚线与实线构成的三角形区域。可以看到当表面上存在点缺陷时，和关系直线全部处于这个三角形区域内，说明了在适当的生长条件下，晶体的快速生长必然伴随这些点缺陷的出现。通过表达式(3.28)计算出的各种点缺陷形成焓均为负，表明了TiO2形成点缺陷结构均为放热反应，形成的缺陷结构较为稳定。在氧化性气氛下，偏离化学计量比的最稳定结构为Oi和缺陷结构；在还原性气氛中，最稳定的结构却是Tii和缺陷结构。由此可知，采用第一性原理的赝势平面波法方法分析点缺陷是切实可行，其研究结果对实验研究具有一定的参考价值和指导意义。【练习与思考】3-12.利用MaterialsStudio软件建立一个两端闭合的(6，6)的碳纳米管，计算其能带结构和光学性质？并将其中的部分碳原子用硼和氮原子取代，计算其能带结构和光学性质，并分析之。3-13.建立金红石型TiO2(110)和锐钛矿型TiO2(101)面，采用第一性原理的方法研究表面几何结构及稳定性。3.6习题3-1.当前流行的赝势平面波计算软件主要有哪些，试比较其共同点和各自应用特点。3-2.建立水分子吸附在锐钛矿型TiO2(101)表面的模型，计算其吸附能并分析结构稳定性。3-3.查阅有关文献，建立金属型碳纳米管结构模型，将其中一个碳原子用钛原子取代，讨论该钛原子吸附氢的能力。3-4.氧化物半导体SnO2(110)表面氧缺陷的几何结构研究。3-5.计算纳米晶GaAs的电子结构和光学性质，分析其光学性能。3-6.纳米金属Fe(100)表面磁学特性研究。3.7推荐阅读材料[1]SegallMD,LindanPLD,ProbertMJ,etal.First-principlessimulation:ideas,illustrationsandtheCASTEPcode.JPhys:CondensMatter,2002,14:2717−2744.[2]Al-KhatatbehY,LeeKKM,KieferB.High-pressurebehaviorofTiO2asdeterminedbyexperimentandtheory.PhysRevB,2009,79:134114-134122[3]HeYB,DulubO,ChengHZ,SelloniAandDieboldU,EvidenceforthePredominanceofSubsurfaceDefectsonReducedAnataseTiO2(101).PhysRevLett.2009,102:106105-106108[4]GaiYQ,LiJB,LiSS,XiaJB,WeiSH.DesignofNarrow-GapTiO2:APassivatedCodopingApproachforEnhancedPhotoelectrochemicalActivity.PhysRevLett.2009,102:036402-036405[5]WendtS,SprungerPT,LiraE.TheRoleofInterstitialSitesintheTi3dDefectStateintheBandGapofTitania,Science,2008,320:1755-1759图3-1截断能示意图[010][�EMBEDEquation.3���01][101]Ti5cO3cTi6cO3c′O2cO3c〞图3-24计算结果的分析对话框图3-25金红石型TiO2结构，图中黑色大球为O原子，灰色小球为Ti原子。图3-29锐钛矿TiO2体相晶胞及垂直[101]方向的六种剪切位置，其中(黑色大球为O原子，灰色小球为Ti原子）O2cO3cO3c’Ti5cTi6c�EMBEDEquation.DSMT4\*MERGEFORMAT����EMBEDEquation.DSMT4\*MERGEFORMAT����EMBEDEquation.DSMT4\*MERGEFORMAT����EMBEDEquation.DSMT4\*MERGEFORMAT����EMBEDEquation.DSMT4\*MERGEFORMAT����EMBEDEquation.DSMT4\*MERGEFORMAT����EMBEDEquation.DSMT4\*MERGEFORMAT����EMBEDEquation.DSMT4\*MERGEFORMAT����EMBEDEquation.DSMT4\*MERGEFORMAT����EMBEDEquation.DSMT4\*MERGEFORMAT����EMBEDEquation.DSMT4\*MERGEFORMAT����EMBEDEquation.DSMT4\*MERGEFORMAT����EMBEDEquation.DSMT4\*MERGEFORMAT���精选范本_1234568018.unknown_1234568082.unknown_1234568116.unknown_1234568132.unknown_1234568148.unknown_1234568159.unknown_1234568167.unknown_1234568171.unknown_1234568173.unknown_1234568175.unknown_1234568176.unknown_1234568174.unknown_1234568172.unknown_1234568169.unknown_1234568170.unknown_1234568168.unknown_1234568163.unknown_1234568165.unknown_1234568166.unknown_1234568164.unknown_1234568161.unknown_1234568162.unknown_1234568160.unknown_1234568152.unknown_1234568157.unknown_1234568158.unknown_1234568154.unknown_1234568156.bin_1234568155.unknown_1234568153.unknown_1234568150.unknown_1234568151.unknown_1234568149.unknown_1234568140.unknown_1234568144.unknown_1234568146.unknown_1234568147.unknown_1234568145.unknown_1234568142.unknown_1234568143.unknown_1234568141.unknown_1234568136.unknown_1234568138.unknown_1234568139.unknown_1234568137.unknown_1234568134.unknown_1234568135.unknown_1234568133.unknown_1234568124.unknown_1234568128.unknown_1234568130.unknown_1234568131.unknown_1234568129.unknown_1234568126.unknown_1234568127.unknown_1234568125.unknown_1234568120.unknown_1234568122.unknown_1234568123.unknown_1234568121.unknown_1234568118.unknown_1234568119.unknown_1234568117.unknown_1234568098.bin_1234568106.unknown_1234568110.unknown_1234568114.unknown_1234568115.unknown_1234568112.unknown_1234568113.unknown_1234568111.unknown_1234568108.unknown_1234568109.unknown_1234568107.bin_1234568102.unknown_1234568104.unknown_1234568105.unknown_1234568103.unknown_1234568100.unknown_1234568101.unknown_1234568099.unknown_1234568090.bin_1234568094.bin_1234568096.bin_1234568097.bin_1234568095.bin_1234568092.bin_1234568093.bin_1234568091.bin_1234568086.bin_1234568088.bin_1234568089.bin_1234568087.bin_1234568084.unknown_1234568085.unknown_1234568083.unknown_1234568050.unknown_1234568066.unknown_1234568074.unknown_1234568078.unknown_1234568080.unknown_1234568081.unknown_1234568079.unknown_1234568076.unknown_1234568077.unknown_1234568075.unknown_1234568070.unknown_1234568072.unknown_1234568073.unknown_1234568071.unknown_1234568068.unknown_1234568069.vsd�������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������IR���_1234568067.unknown_1234568058.unknown_1234568062.unknown_1234568064.unknown_1234568063.unknown_1234568060.unknown_1234568061.unknown_1234568059.unknown_1234568054.unknown_1234568056.unknown_1234568057.unknown_1234568055.unknown_1234568052.unknown_1234568053.unknown_1234568051.unknown_1234568034.unknown_1234568042.unknown_1234568046.unknown_1234568048.unknown_1234568049.unknown_1234568047.unknown_1234568044.unknown_1234568045.unknown_1234568043.unknown_1234568038.unknown_1234568040.unknown_1234568041.unknown_1234568039.unknown_1234568036.unknown_1234568037.unknown_1234568035.unknown_1234568026.unknown_1234568030.unknown_1234568032.unknown_1234568033.unknown_1234568031.unknown_1234568028.unknown_1234568029.unknown_1234568027.unknown_1234568022.unknown_1234568024.unknown_1234568025.unknown_1234568023.unknown_1234568020.unknown_1234568021.unknown_1234568019.unknown_1234567953.unknown_1234567985.unknown_1234568001.unknown_1234568010.unknown_1234568014.unknown_1234568016.unknown_1234568017.unknown_1234568015.unknown_1234568012.unknown_1234568013.unknown_1234568011.unknown_1234568005.unknown_1234568007.unknown_1234568009.unknown_1234568008.unknown_1234568006.unknown_1234568003.unknown_1234568004.unknown_1234568002.unknown_1234567993.unknown_1234567997.unknown_1234567999.unknown_1234568000.unknown_1234567998.unknown_1234567995.unknown_1234567996.unknown_1234567994.unknown_1234567989.unknown_1234567991.unknown_1234567992.unknown_1234567990.unknown_1234567987.unknown_1234567988.unknown_1234567986.unknown_1234567969.unknown_1234567977.unknown_1234567981.unknown_1234567983.unknown_1234567984.unknown_1234567982.unknown_1234567979.unknown_1234567980.unknown_1234567978.unknown_1234567973.unknown_1234567975.unknown_1234567976.unknown_1234567974.unknown_1234567971.unknown_1234567972.unknown_1234567970.unknown_1234567961.unknown_1234567965.unknown_1234567967.unknown_1234567968.unknown_1234567966.unknown_1234567963.unknown_1234567964.unknown_1234567962.unknown_1234567957.unknown_1234567959.unknown_1234567960.unknown_1234567958.unknown_1234567955.unknown_1234567956.unknown_1234567954.unknown_1234567921.unknown_1234567937.unknown_1234567945.unknown_1234567949.unknown_1234567951.unknown_1234567952.unknown_1234567950.unknown_1234567947.unknown_1234567948.unknown_1234567946.unknown_1234567941.unknown_1234567943.unknown_1234567944.unknown_1234567942.unknown_1234567939.unknown_1234567940.unknown_1234567938.unknown_1234567929.unknown_1234567933.unknown_1234567935.unknown_1234567936.unknown_1234567934.unknown_1234567931.unknown_1234567932.unknown_1234567930.unknown_1234567925.unknown_1234567927.unknown_1234567928.unknown_1234567926.unknown_1234567923.unknown_1234567924.unknown_1234567922.unknown_1234567905.unknown_1234567913.unknown_1234567917.unknown_1234567919.unknown_1234567920.unknown_1234567918.unknown_1234567915.unknown_1234567916.unknown_1234567914.unknown_1234567909.unknown_1234567911.unknown_1234567912.unknown_1234567910.unknown_1234567907.unknown_1234567908.unknown_1234567906.unknown_1234567897.unknown_1234567901.unknown_1234567903.unknown_1234567904.unknown_1234567902.unknown_1234567899.unknown_1234567900.unknown_1234567898.unknown_1234567893.unknown_1234567895.unknown_1234567896.unknown_1234567894.unknown_1234567891.unknown_1234567892.unknown_1234567890.unknown