第3章赝势平面波方法(I)汇总.pdf

《第3章赝势平面波方法(I)汇总.pdf》由会员分享，可在线阅读，更多相关《第3章赝势平面波方法(I)汇总.pdf（37页珍藏版）》请在三一文库上搜索。

1、第 3 章赝势平面波方法(I) Computational Materials Science ：From Basic Principle to Practical Design Methodology - 49 - 第3章赝势平面波方法 (I) 基于密度泛函理论的赝势平面波方法可以计算很大范围不同体系的基态属性，它采用了 平面波来展开晶体波函数，用赝势方法作有效的近似处理。由于平面波具有标准正交化和能 量单一性的特点，对任何原子都适用且等同对待空间中的任何区域，不需要修正重叠误差。 因此平面波函数基组适合许多体系，其简单性使之成为求解Kohn-Sham 方程的高效方案之 一。另外， 赝势的引

2、入可以保证计算中用较少的平面波数就可以获得较为可靠的结果。该方 法具有较高的计算效率，使之日益发展成为有效的计算方法。本章首先对赝势平面波方法进 行重点讨论，其次介绍了基于第一性原理计算软件一般步骤，最后结合Materials Studio 软 件包应用，对锐钛矿型TiO2(101)表面及其点缺陷结构进行建模和计算。 3.1基本原理 基于密度泛函理论的第一性原理计算实质是求解Kohn-Sham 方程。 实际求解 Kohn-Sham 方程时， 由于原子核产生的势场项在原子中心是发散的，波函数变化剧烈，需要采用大量的 平面波展开， 因而计算成本变得非常大，所以在计算中选取尽可能少的基函数。计算中选

3、择 的基函数与最终波函数较接近则收敛较快，当然包含的维度也应该尽量少。众所周知， 根据 研究对象不同， 选择基函数的方法也不同的，如原子轨道线性组合法(LCAO-TB) 、正交平面 波法 (OPW)、平面波赝势法(PW-PP)、缀加平面波法(APW) 、格林函数法(KKR) 、线性缀加 平面波法 (LAPW) 、Muffin-tin轨道线性组合法(LMTO) 等，选取典型代表方法在随后的章节 中重点展开讨论。与LAPW ， LMTO等精度较高的第一性原理计算方法比较，平面波赝势 法是计算量较少的方法，适用于计算精度要求不严格，因原胞较复杂而导致计算量陡增加的 体系。 为此，本章将重点学习赝势平

4、面波方法，先学习电子能带的平面波基底展开以及赝势 等相关基本概念，然后再讨论赝势引入原理。 3.1.1平面波展开与截断能 1.平面波展开 平面波是自由电子气的本征函数，由于金属中离子芯与类似的电子气有很小的作用，因 此很自然的选择是用它描述简单金属的电子波函数。众所周知， 最简单的正交、 完备的函数 集是平面波exp ()i kGr，这里G是原胞的倒格矢。根据晶体的空间平移对称性，布洛赫 (Bloch) 定理 (将在第 4.1.1 节中说明 )证明，能带电子的波函数( , )r k 总是能够写成 ( , )( )exp()r krik r (3.1) 式中k是电子波矢，( )r是具有晶体平移周

5、期性的周期函数。对于理想晶体的计算，这是很 自然的， 因为其哈密顿量本身具有平移对称性，只要取它的一个原胞就行了。对于无序系统 (如无定型结构的固体或液体)或表面、界面问题，只要把原胞取得足够大，以至于不影响系 统的动力学性质，还是可以采用周期性边界条件的。因此， 这种利用平移对称性来计算电子 结构的方法， 对有序和无序系统都是适用的。采用周期性边界条件后，单粒子轨道波函数可 第 3 章赝势平面波方法(I) Computational Materials Science ：From Basic Principle to Practical Design Methodology - 50 - 以

6、用平面波基展开为 1 ( )()exp( () G rGi KGr N (3.2) 式中1N是归一化因子， 其中是原胞体积； 这里G是原胞的倒格矢，K是第一 Brillouin 区的波矢，()G 是展开系数。 Bloch 定理表明，在对真实系统的模拟中，由于电子数目的无 限性，K矢量的个数从原则上讲是无限的，每个K矢量处的电子波函数都可以展开成离散 的平面波基组形式，这种展开形式包含的平面波数量是无限多的。基于计算成本的考虑，实 际计算中只能取有限个平面波数。采用的具体办法是一方面由于( )r随K点的变化在K点 附近是可以忽略的，因此我们可以使用K点取样通过有限个K点进行计算。另一方面，为 了

7、得到对波函数的准确表示，G矢量的个数也应该是无限的，但由于对有限个数的G矢量 求和已经能够达到足够的准确性，因此对G的求和可以截断成有限的。给定一个截断能 22 () 2 cut GK E m (3.3) 对G的求和可以限制在 2 () / 2 cut GKE的范围内，即要求用于展开的波函数的能量小 于 cut E。当0K时，即在点，有很大的计算优势，因为这时波函数的相因子是任意的， 就可以取实的单粒子轨道波函数。这样， 对 Fourier 系数满足关系式 * ()() ll GG，利用这 一点，就可以节约不少的计算时间。 2.截断能选取原则 为了取有限个的平面波数，通常的做法是确定一个截断能

8、量(Energy cutoff) ，如图3-1 所示，此时函数基组并不完备，总能量计算会产生相应误差， 通过增加截断能量可以减小误差幅度。为了使计算出的体系总 能量达到设定精度，一般截断能量必须选取到足够高。有限平 面波基组的误差可以加以校正，较好的解决方法是引入一个校 正因子 (correction factor)，由此可以在一个恒定数量基组下进行 计算，即使采用了恒定的截止能量这个强制条件也可以校正相 应 的 计 算 结 果 。 进 行 这 种 校 正 所 需 要 的 唯 一 的 参 数 就 是 ln tot cut dE dE ，Etot是体系总能量，Ecut是截止能量。例如，当它的数

9、值小于 0.01 eV/atom 时，计算就达到了良好的收敛精度，对于大多数计算0.1 eV/atom 就已 足够。 3.平面波基展开特征 用平面波基来展开电子波函数是因为用平面波基来计算有很多优点。平面波基能很方便 地采用快速傅里叶变换(FFT)技术，使能量、力等的计算在实空间和倒空间快速转换，这样 计算尽可能在方便的空间中进行。如前面讲到的哈密顿量中的动能项的矩阵元，在倒空间中 只有对角元非零，就比实空间减少了工作量。第二， 平面波基函数的具体形式并不依赖于核 的坐标。这样，一方面，价电子对离子的作用力可以直接应用Hellman-Feymann 定理 (将在 3.1.5 节中进行说明)得到

10、解析的表达式，计算显得非常方便。另一方面也使总能量的计算在 不同的原子构型下有基本相同的精度。此外平面波计算的收敛性和精确性比较容易控制，因 图 3-1 截断能示意图 第 3 章赝势平面波方法(I) Computational Materials Science ：From Basic Principle to Practical Design Methodology - 51 - 为通过截断能Ecut的选择可以方便地改变平面波基的多少。当然平面波基也有缺点，一般电 子轨道具有一定的局域性，而平面波是空间均匀的，因此电子轨道展开时与原子轨道基相比， 平面波基的个数要多得多。为了尽量减少平面波基

11、的个数，一般在平面波的计算中都采用赝 势(pseudopotentials)来描述离子实与价电子之间的相互作用，使电子轨道波函数在离子实内 部的分布尽量平缓些。下面将讨论赝势概念及其引入思路。 3.1.2赝势 1.赝势引入 平面波函数作为展开基组具有很多优点，然而截断能的选取与具体材料体系密切相关。 由于原子核与电子的库仑相互作用在靠近原子核附近具有奇异性，导致在原子核附近电子波 函数将剧烈振荡。 因此，需要选取较大的截断能量才能正确反映电子波函数在原子核附近的 行为，这势必大大地增加计算量。另一方面， 在真正反映分子或固体性质的原子间成键区域， 其电子波函数较为平坦。基于这些特点， 将固体看

12、作价电子和离子实的集合体，离子实部分 由原子核和紧密结合的芯电子组成，价电子波函数与离子实波函数满足正交化条件，由此发 展出所谓的赝势方法。1959 年，基于正交化平面波方法，Phillips 和 Kleinman 提出了赝势的 概念。 基本思路是适当选取一平滑赝势，波函数用少数平面波展开，使计算出的能带结构与 真实的接近。 换句话说， 使电子波函数在原子核附近表现更为平滑，而在一定范围以外又能 正确反映真实波函数的特征，如图3-2 所示。 所谓赝势， 即在离子实内部用假想的势取代真实的势，求解波动方程时，能够保持能量 本征值和离子实之间的区域的波函数的不变。原子周围的所有电子中，基本上仅有价

13、电子具 有化学活性， 而相邻原子的存在和作用对芯电子状态影响不大。这样， 对一个由许多原子组 成的固体， 坐标空间根据波函数的不同特点可分成两部分(假设存在某个截断距离 c r)。(1) c r 以内的核区域， 所谓的芯区。 波函数由紧束缚的芯电子波函数组成，对周围其它原子是否存 在不敏感，即与近邻的原子的波函数相互作用很小；(2) c r以外的电子波函数(称为价电子波 函数 )承担周围其它原子的作用而变化明显。 2.原子赝势 全电子 DFT 理论处理价电子和芯电子时采取等同对待，而在赝势中离子芯电子是被冻 结的，因此采用赝势计算固体或分子性质时认为芯电子是不参与化学成键的，在体系结构进 行调

14、整时也不涉及到离子的芯电子。在赝势近似中用较弱的赝势替代芯电子所受的强烈库仑 势，得到较平缓的赝波函数，此时只需考虑价电子，在不影响计算精度情况下，可以大大降 低体系相应的平面波截断能Ecut，从而降低计算量。图3-3 为 Si 原子赝势示意图。赝原子用 于描述真实原子自身性质时是不正确的，但是它对原子-原子之间相互作用的描述是近似正 确的。近似程度的好坏，取决于截断距离 c r的大小。 c r越大，赝波函数越平缓，与真实波 函数的差别越大，近似带来的误差越大；反之， c r越小，与真实波函数相等的部分就越多， 近似引入的误差就越小。 可将真实价波函数( ,) n r k看作是由赝势波函数(

15、,) n r k和内层波函数( ,) J r k线性组 合，即 第 3 章赝势平面波方法(I) Computational Materials Science ：From Basic Principle to Practical Design Methodology - 52 - ( ,)( ,)()( ,) nnnJJ J r kr kkr k(3.4) 其中系数() nJ k可由正交条件 * ( ,)( ,)0 Jn drr kr k 确定，即 * ( )( , )( , ) nJJn kdrr kr k 图 3-2 赝波函数与势图 3-3 Si 原子赝势示意图 联合真实波函数( ,) n

16、 r k 所满足的薛定谔方程： ( )( , )( )( , ) nnnTV rr kEkr k 可得到赝波函数满足如下方程 ( , )( )( , ) psnnnTUr kEkr k ( ,)( )( , )( ,)(, ) psnnRn Ur kV rr kdr Vr rrk(3.5) 其中 * ( ,)(,)()(,) RJnJJ J Vr rrkEkErk ps U称为原子赝势。根据密度泛函理论，原子赝势包括离子赝势 ion ps U和价电子库仑势和交换 关联势：( )( ) ionps pspsHxc UUVrVr，其中后两项( ) ps H Vr和( ) xc Vr可以从真实电荷密

17、度计算， 此时等于对应的全电子势( )V r和 xc V 。 从上面可知，赝势应具有以下特征：(1)赝波函数和真实波函数具有完全相同的能量本 征值( ) n Ek，这是赝势方法的重要特点；(2)赝势第二项是排斥势，与真实的吸引势有相消趋 势，因此比真实势弱； (3)赝势包括局域项， 其中非局域项同时与r和r处的赝波函数( ,) n r k 和(,) n rk有关，而且依赖于能量本征值() n Ek。 3.赝势分类 从上面的推导可以看出，赝势实际上是一种算子，但可以近似的将它处理成原子间距离 的简单函数， 叫做局部赝势。 其结果是只要给定周期性排列的原子的位置，系统的能量就可 以算出来。确定赝势

18、的方法不是唯一的，主要有以下几种： (1) 经验赝势方法， 利用实验数据拟合有限几个V(K) 的值，主要用途是在现代从头算原 子赝势自洽迭代计算中作初始值使用。 (2) 模型赝势是半经验的，能用于自洽计算的原子赝势。 第 3 章赝势平面波方法(I) Computational Materials Science ：From Basic Principle to Practical Design Methodology - 53 - (3) 模守恒赝势是第一性原理从头算原子赝势，是核与芯电子联合产生的有效势，是从 原子的薛定谔方程从头计算得到的，它可以给出价电子或类价电子的正确电荷分布，适于作

19、自洽计算。 (4) 模非守恒赝势(超软赝势 )，其优点是容易选择芯区的赝势波函数，减少了必须的平 面波函数的数目，较大地减轻了计算工作量。 目前在第一性原理计算中应用较多的为模守恒赝势(Norm-conserving pseudopotentials， NCPP)和超软赝势 (Ultra Soft pseudopotentials ，USPP)两种方案。下面将分别讨论。 3.1.3模守恒赝势 1.模守恒赝势构造 在第一性原理计算中，常采用由Hamann D R 等提出的模守恒赝势。这种赝势所对应的 波函数不仅与真实势对应的波函数具有相同的能量本征值，而且在 c r (原子芯半径 )以外与真 实

20、波函数的形状和振幅都相同(即模守恒 )，且在 cr 以内比较平缓。采用赝势计算关键在于可 以有效的对化学键的价电子进行可再现的近似，赝势与全势在超过离子实半径以外具有完全 相同的函数形式。这种赝势能生成正确的电荷密度，适合作自洽计算。 构造模守恒赝势基本思想是选择某个特定的电子排部状态(不一定就是基态)全部电子 计算在一个孤立的原子中进行，从而得到原子价电子能量本征值和价电子波函数。选择一个 离子赝势或赝波函数参数形式，通过对参数的调节，使得赝原子计算和全电子原子赝势计算 采用相同的交换-相关势，在超过截止半径 cr 后与价电子波函数形式相同，赝势的本征值等 于价电子的本征值。如果电子波函数和

21、赝势波函数满足正交归一，两者在截止半径以外的匹 配性决定了模守恒条件自动成立。模守恒赝势要求赝势波函数满足：(1)本征值与真实本征 值相等； (2)没有节点； (3)在原子核区以外( crr )与真实波函数相等； (4)在内层区 ( crr )内 的赝电荷与真实电荷相等，将赝波函数插入到薛定谔方程中即得对应的赝势。一般说来， 小 的 c r 移植性好，可用于不同环境，但平面波收敛慢。 第一性原理模守恒赝势可分为局域和非局域两部分， ()(,) pslocvNLvv vv UVrRUrR rR(3.6) 其中 v 是对离子势求和。考虑到原子球对称性，得用球谐函数将赝势的非局域 部分写成： * ,

22、 ( ,)(,)( ,)( , )( ,) NLlmlmll l ml m Ur rYYV r rlmlm Vr r 如果将( ,) l V r r取成半局域形式，即径向是局域的，只有角部分是非局域的： ( ,)( ) () ll Vr rV rrr，并定义角动量l的投影算符 l m Plmlm 。则半局域的原子赝势可 以写成如下的形式： ( )()() pslocvlvl vvl UrVrRV rR P，(3.7) 为了简化计算，Kleinman 和Bylander(KB) 将上面半局域赝势部分用一个非局域赝势来近 似： , ( ,) lmlm NL l m mlm VV Ur r V (3

23、.8) 2.模守恒赝势应用特征 第 3 章赝势平面波方法(I) Computational Materials Science ：From Basic Principle to Practical Design Methodology - 54 - 模守恒赝势方法可以在局域密度近似下，采用平面波基精确有效地计算固态性质，且可 移植性好， 但在描述局域价轨道平面波基仍然很大，因而在第族元素和过渡族金属中的应 用受到了限制。 通过优化光滑的赝波函数或赝势和增大截断半径 c r 的改进有一定效果，但模 守恒条件的限制使得在一些情况下，如O 原子 2p 或 Ni 原子 3d 轨道，很难构造出比全电子

24、波函数更光滑的赝波函数，收敛仍然很慢。 模守恒赝势最早由Hamann D R 等提出，后来建立一组涵盖整个周期表的参数，之后发 展出的 Kerker、TM 赝势和 Optimised 赝势，都是在朝着兼顾准确性的情况下，尽可能使必 须使用的平面波基底数目越少越好，平面波基底数是直接影响着所需计算量大小的量。一个 赝势所需的基底数多少，可由 tot E对 cut E的收敛性来判断， 即平面波截断动能 cut E用到多大时 则固态计算所求得系统总能就不再改变，所需 cut E越小，也就是所谓的赝势越“ 软 ” 。 使用 Optimised 或 TM 赝势虽然能够把模守恒型赝势变的很“ 软” ，但模

25、守恒条件对于原 本就已经没有节点价电子云分布的改造及最佳化的程度，与现今日渐普遍的超软赝势(它不 必遵守模守恒条件)来比，节省计算的程度仍是有限的。总之，计算量的大小是取决于原子 的种类这一点， 是十分明确而普遍的认识，也就是说不同种类元素其势的“ 软硬 ” 差异会令人 明显感觉到。 3.1.4超软赝势 1.超软赝势构造 对于过渡族元素和第一周期元素，模守恒赝势不能明显降低所需平面波截断能Ecut。 Vanderbilt 提出了超软赝势，其赝波函数在核心范围是被作成尽可能平滑，可以大幅度地减 少截断能， 即可使计算所需的平面波函数基组更少。就技术上而言， 这是靠放宽模守恒的要 求，采用广义的正

26、交条件来达成的。为了重建整个总的电子密度，波函数平方所得到电荷密 度必须在核心范围再附加额外的密度进去。这个电子云密度由此就被分成两个部分，第一部 分是一个延伸在整个单位晶胞平滑部分，第二部分是一个局域化在核心区域的自旋部分。前 面所提的附加部分是只出现在电子密度，并不在波函数。 超软赝势中总能量与采用其它赝势平面波方法时相同，非定域势 NLV 表达如下 (0) , II NLnmnm nm I VD(3.9) 式中投影算符和系数 (0) D分别表征赝势和原子种类的差别，指数I 对应于一个原子位置。 总能量用电子密度可以表示为： 2 ( ) , ( )( )( ) III inminmi in

27、m I n rrQr(3.10) 式中是波函数，( )Q r 是严格位于芯区的附加函数。超软赝势完全由定域部分，( ) ion loc Vr 和系 数 (0) D，( )Q r，确定。赝势是通过引进一系列正交条件来建立的， ijij S，S是 哈密顿重叠算符，可以表示为1 II nmnm Sq，系数q是通过对 ( )Q r积分得到。 从而， 超 软赝势的Kohn-Sham 方程可以写为： iii HS，H可以表示为动能和定域势能之和 () , III effnmnm nm I HTVD(3.11) ( )(0) ( )( ) II nmnmeffnm DDdrVr Qr(3.12) 第 3 章

28、赝势平面波方法(I) Computational Materials Science ：From Basic Principle to Practical Design Methodology - 55 - 2.应用特点 超软赝势产生算法保证了在预先选择的能量范围内会有良好的散射性质，这导致了赝势 更好的转换性与精确性。超软赝势通常也借着把多套每个角动量通道当作价电子来处理浅的 内层电子态。 这也会使精确度跟转换性更加提升，虽然计算代价会比较高。与模守恒赝势对 比，不同之处在于在超软赝势中存在重叠算符S，波函数与D有关。而且投影算符函数数 量要比模守恒赝势中大两倍多。与附加电荷相关的一系列计算

29、可以在实空间中进行，这与函 数中定域势的性质有关，而多余的步骤不会对计算效率产生较大的影响。在Laasonen 文献 中提供了超软赝势计算的详细方法以及总能量微分表达式。 3.1.5Hellmann-Feynman力 Hellmann 和 Feynman 在量子力学框架下给出了作用于离子实上(位置坐标为 I R )的力 I F 。离子受的力为总能对离子位置的偏导， I I E F R (3.13) E 作为系统哈密顿量的能量本征值，满足Kohn-Sham 方程， HE 可以得到： EH(3.14) 将式 (3.14)代入 (3.13)得 I II H FE RR (3.15) 由于是一个归一化

30、常数，上式的第一项等于零。最终得到作用在离子上的力 I H F R (3.16) 这就是著名的Hellmann-Feynman 定理 Hellmann-Feynman 定理计算出的力是和电子波函数相联系的，它的误差与波函数误差的 一级修正量成正比，只有波函数非常接近真实的本征态时这个力才是精确的。所以在计算时 需要同时考虑到离子弛豫和电荷密度自洽，即，离子在受力后到达一个新的位置，此时电子 也需要接近瞬间基态，然后在新的离子位置和新的电子密度下进行计算，直至总能到达局部 极小值。 在得出离子受的受力后，需要对离子进行弛豫，即需要知道离子弛豫的方向和大小。 【练习与思考】 3-1.查阅有关文献和

31、书籍， 找出 3 种以上不同的经验赝势，试分析他们的适用对象及特点？ 3-2.查阅 Martin R M Electronic Structure 一书，写出模守恒赝势方法中模守恒的条件。 3-3.查阅谢希德 固体能带理论 和有关书籍， 试说明赝势平面波法采用了那些近似处理？ 第 3 章赝势平面波方法(I) Computational Materials Science ：From Basic Principle to Practical Design Methodology - 56 - 3.2数值处理方法与技巧 3.2.1超原胞方法 对于固体体系，目前主要用三种模型来模拟材料特性，一是团簇

32、模型(Cluster Model) ， 二是嵌入模型 (Embedded Cluster Model) ，三是层状模型(Slab Model) 。在以往的量子化学计 算中，研究人员往往使用团簇模型来模拟固体材料，这种简化在化学上的确存在一定依据， 因为根据化学家的直觉，一个分子体系作用是受局部相互作用支配的，从这一观点出发，可 以用分子与原子簇的作用来反映分子与固体相互作用的性质。由于计算方法和计算条件的限 制，没有考虑表面结构弛豫的影响，从而使计算的模型体系和实验体系存在较大的误差。第 二种嵌入模型的提出主要是为了克服团簇模型在模拟材料表面时存在边界性问题，但该方法 在计算过程中涉及到大量近

33、似，需要针对不同的体系使用不同的计算方法，比如由Korringa J、Kohn W 和Rostoker N提出的格林函数(Green Function) 方法和戴逊方程(Dyson Equation) ， 所以该模型并不为广泛接受。由此， 人们提出了超原胞模型，该模型在模拟体系时采用了周 期性边界条件， 特别适合研究金属、半导体这类具有周期性的凝聚态体系。在后面的实例中 基本上都是选取了超原胞模型。 层状模型 (Slab Model) 中的超原胞模型将体系看作沿晶格矢周期性排列的体系，在计算 中所研究的原子都放在超原胞中，原子坐标或者其对应的周期性位置可用下式表示： s a x a l l l

34、 l R RR(3.17) 其中 a lsl 为原子在超原胞中的坐标， s l为原子种类的序号， a l为多个同类原子之间的序号。R 为格矢。 目前许多第一性原理计算软件采取了超原胞模型，来构造周期性结构，包括三维或低维 周期性结构。对于特殊体系如掺杂、缺陷、表面等，采取多倍原胞进行平移扩展，以保证物 理上相邻原胞中的原子或分子没有相互作用。例如研究表面的分子吸附可假设它们在一个 “ 盒子 ” 里面成为周期体系，层与层之间用足够厚度的真空层隔离以忽略在盒子间原子的相互 作用，如图 3-4所示；再如，在研究中使用超原胞，认为它是可以在三维方向无限拓展，超 原胞是没有外形的限制，假如这个晶体具有高

35、点群的对称性，则它也可以用来加速计算。锐 钛矿型 TiO2 (101)的超原胞如图 3-5所示。 图 3-4 研究表面的分子吸附原胞 图 3-5 锐钛矿型 TiO2(101)表面原子层超原胞 010 101 101 Ti5c O3c Ti6 O3c O2 O3c 第 3 章赝势平面波方法(I) Computational Materials Science ：From Basic Principle to Practical Design Methodology - 57 - 3.2.2自洽电子弛豫方法 表示电子结构松弛有多种方法，其中密度混合方法是最有效的。该方法是使用在一个在 固定位势之下

36、将电子能量本征值的总和极小化，而不是将总能做自洽式的极小化。在步骤的 最后新的电荷密度就会与初始的电荷密度来混合产生新的电荷密度再计算体系能量以重复 叠代直到收敛为止。密度混合有以下几种形式： (1) 线性混合，新的密度是上一次输入与输出密度的线性组合。 (2) Kerker 混合，可用下式表达： 2 22 max ()()()() newinoutin G GGAGG GG (3.18) 其中G为截断波矢，A为混合幅度。 (3) Pulay 形式，是最有效的一种方式，在Pulay 形式中新的输入密度上所有先前迭代步 中密度的组合，它不仅与A和 max G有关，而且和所有的先前迭代步有关。 (

37、4) Broyden 混合形式，它与Pulay 形式有类似之处。 将本征值之和来做最小化既可以使用共扼梯度方法也可以使用加权余量方法。电子波函 数是以平面波基底来表示，并且展开系数逐渐会被变化以便达到最小的总能。此极小化可利 用每个波函数取独立的最佳化的band-by-band 的技术，或允许同时更新所有波函数的 all-band 方法来达成。此方式用了Payne 等人所提出的预先调节式的共轭梯度技术。传统总 能极小化方法可能在具有晶胞一个方向上拉长的金属系统中计算会不穏定。而这是在表面做 超晶格计算的典型设置中无法避免的。密度泛函方法对于绝缘体跟金属的状况都一样收敛良 好。密度混合方法对于中

38、等大小的绝缘体系统甚至都还提供3 倍快的加速。 密度混合方式主 要的优势是当处在金属系统时可在相当少的次数很可靠的收敛。 3.2.3几何结构优化技术 对于给定各原子位置、元素种类的研究体系，通过密度泛函理论自洽求解Kohn-Sham 方程可以得到整个系统处于多电子基态时的总能。总能量对系统虚拟微位移的导数就是各原 子的受力，即Hellmann-Feynman 力。这为我们理论预言物质的结构提供了一种行之有效的 方法。因为自然界稳定的结构应该具有最低的总能，只要根据原子受力来调整原子的位置， 直到整个体系的总能达到最低，所有原子受力为零。当然在实际的计算过程中，给出希望达 到而且有限的计算精度，

39、可找到能量面的(全局 )最小值，这时所对应的物质结构就是自然界 最稳定的结构，该过程被称为几何结构优化(Geometry Optimization) 。为了确保搜索能量面 的最小值时能找到全局最小而不是局域最小，并提高整个搜索过程的效率，需要一些强有力 的搜索算法以使原子最快地运动到最稳定结构的位置。 能量最小化算法一般分为两大类：全局极小和局部极小。全局极小算法可以得到基态构 型，如模拟退火和遗传算法；局部极小算法找寻的是亚稳态结构，最常用的方法有直接能量 最小化、最陡下降法、共轭梯度法、Newton-Rapshon 方法、阻尼动力学方法等等。这里主 要介绍材料模拟常用的三种优化方法。 1.

40、最陡下降法 (Steepest Descent Method) 第 3 章赝势平面波方法(I) Computational Materials Science ：From Basic Principle to Practical Design Methodology - 58 - 最陡下降法沿着局部净受力方向行走，以进行能量极小。从初始点开始， 沿着局部梯度 的反方向 l g ，并通过在此方向上的一维极小化，移动到该方向的极小点，再从这个点，开 始重复以上过程，直到达到所要求的精度。最速下降法在远离极小点效率很好，在接近极小 时效率不高，而且沿梯度方向每前进一步将对接下来一步都引入一个正比于它

41、梯度的误差， 常常只在优化的最初几步使用这种方法。 2.共轭梯度法 (Conjugale-gradient Method) 共轭梯度方法克服了最陡下降法的困难。在此方法中， 每相邻两次优化的起始点的 lg 仍 是正交的，但优化方向 l v 由当前梯度 l g 结合前一次优化方向 1l v 和梯度 1l g 共同决定， l v 和 1lv 互为共轭： 1llll vgd(3.19) 其中是一个标量数，由前一次优化起点的梯度 1l g 和优化终点的梯度(即当前时刻梯 度) l g 共同决定，不同的算法给出各自的确定公式： Flrtcher-Reeves 算法： 11 ll l ll gg gg P

42、olak-Ribiete 算法： 11 11 () lll l ll ggg gg 对于能量函数()f P ，可以按如下方式进行优化：若初始点在 0P ，令000 ()vgf P， 沿 0 v 方向运用一维极小化方法到达该方向的一极小点 1 P ，则 11 ()gf P。由 0 g 和 1 g 可得到 1，则1100 vgv，再沿 1v 找极小，重复以上过程，如果函数是含 N个变量的二次型，则 通过N次一维极小化就可以找到极小。上面两种方法在优化中只用了势能函数的一阶导数， 即梯度。 3.Newton-Rapshon 法 任何给定点的能量都可以展成泰勒展式： 2 1 ()( )( )( )()

43、. 2 U xxU xUxxUxx(3.20) ( )Ux 是在 x 处的一阶导数矢量，( )Ux 为二阶导数矩阵，称为Hessian 矩阵。对泰勒展式只 取到二阶导数，而忽略高阶的，则位移矢量x为： 1 xHg(3.21) 其中( ),( )HUxgUx 一般情况下，可以重复以上过程直到能量最小，这就是所谓的Newton-Raphson 方法。 Newton-Raphson 方法收敛速度很快， 但是不能保证收敛的方向，而且需要计算Hessian 矩阵， 计算量很大。为了避免函数值上升，采用以下修正公式： 1 xHg (3.22) 其中通过一维搜索算法确定，保证函数值向最小方向收敛。 为了避免

44、频繁计算Hessian 矩阵，常采用更新校正方案(updating scheme)，常见算法有： DavidonFletcherPowell(DFP) 和 BroydenFletcherGoldfarb Shanno(BFGS)。具体表达 式如下： 第 3 章赝势平面波方法(I) Computational Materials Science ：From Basic Principle to Practical Design Methodology - 59 - 1 ()() DFPDFP DFPDFP ii ii DFP i HgHgxx HH xggHg (3.23) 1 ()() BFG

45、SBFGS BFGSBFGSBFGS ii iii BFGS i HgHgxx HHg Hg vv xgg Hg (3.24) 其中 BFGS i BFGS i Hgx v xggHg 。 3.2.4快速傅立叶变换 (FFT) 平面波函数中使用的赝势常常是非局域化的。这意味着， 为了要模拟全电子原子的散射 性质，必须要在每个角动量信道里采用不同的散射位势。然而这些位势只在原子核心区域之 内有所不同， 采用实空间的表示方法是有效的。目前一般的计算软件采用的是选择在倒易空 间中使用赝势， 也就是说在倒易空间中进行对波函数和势操作的加总。这对于以平面波基底 系数来表达的波函数是一种很常用的做法。 然

46、而对于大晶胞的Hamiltonian 波函数和势能项进行实空间计算，会变得比在倒易空间 中求值更快。 这是因为赝势的非局域部份只有在核心内的区域才不是零，并且在远比各个单 一原子核心区域还要大的晶胞中几乎到处都是零。因此， 在实空间求值波函数和势的积，并 且使用傅立叶转换来获得在倒易空间的值就变的更为有效率。这就意味着波函数和各种势在 计算中需要一个有效的途径完成从实空间和倒易空间的相互转换，傅立叶变换为此提供了有 力的工具 。 平面波赝势方法中采用一组有限个正交的平面波作为基矢，这就将久期方程中的矩阵元 计算： ijii HkGH KG，归结为傅立叶变换。 快速傅立叶变换在计算中频繁使用。傅

47、立叶变换网格不仅影响计算精度，也影响计算速 度，因此也存在计算结果对傅立叶变换网格的依赖关系，严格地说，这也是一个收敛问题。 傅立叶变换网格空间稀疏，将导致体系不易收敛或收敛结果粗糙，网格越密， 则将大量消耗 计算资源。 3.2.5k 空间取样规则 在通常的能带图中，经常会出现比如 、 、 、 等等的符号。这些符号表示的 是布里渊区内具有高对称性的一些特殊的k 点，这些 k 点有特殊重要的意义。在进行体系总 能计算时， 通常要对布里渊区内的波函数或本征值进行积分，在实际计算过程中，积分是通 过对部分特殊选取的k 点求和完成的。比较常见的k 点网格方法有Monkhorst-Pack 方法和 四面

48、体网格。 根据Bloch定理，周期体系中的电子波函数可以表述为调幅平面波的形式，即 ( , )()exp()r krik r 。 其本征能量和本征矢量为 ( ) n Ek 和 () n k 。不同的电子状态按照量子 数k进行分类，而量子数n 则表征能态的分立性。研究多体体系的价电子问题，归根结底是 计算出不同类的电子状态的本征值和本征矢量，体系处于基态情况下，哪些不同k的低能量 状态被电子占据。因体系具有周期性，所以，第一布里渊区的所有k可以代表所有的k。但 第 3 章赝势平面波方法(I) Computational Materials Science ：From Basic Principl

49、e to Practical Design Methodology - 60 - 是，由于周期边界条件确定的k 有无穷多个， 而且计及相互作用势的实际体系中许可的k 点 在倒易空间内是不均匀的，实际计算过程中只能选取有限个点。第四章还将对k 点的选取方 法进行深入分析。 在赝势平面波计算工作中，有限个k 点在第一布里渊区内等权重均匀选取，这种选取k 点的方法称之为Monkhorst-Pack 方法。实际操作中考虑体系的对称性，将第一布里渊区依 据点对称性划分为几个等价的“不可约空间” ， 自洽计算只在这个较小的不可约空间内进行。 换言之， 将研究那些互不等价的k 量子数的集合， 然后再用以描述整个电子体系的状态。如 果体系的第一布里渊区的不可约空间大小为布里渊区体积的1/n，则总体性质取不可约体积 内的计算结果n 倍即可。 但要注意的是， 由于不可约布里渊区之间相交，第一布里渊区和第 二布里渊区之间也有相交点，所以总有一些点为两个或几个不可约空间共有或为相邻布里渊 区共有，这时如果进行占据态总能量和其它物理性质计算时采用简单倍乘就将导致完全错误 的结果。关于k 点网格选取方法的进一步讨论，请参考第4 章相关内容。 3.2.6基于密度泛函