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

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

基于密度泛函理论的赝势平面波方法可以计算很大范围不同体系的基态属性，它采用了平面波来展开晶体波函数，用赝势方法作有效的近似处理。由于平面波具有标准正交化和能量单一性的特点，对任何原子都适用且等同对待空间中的任何区域，不需要修正重叠误差。因此平面波函数基组适合许多体系，其简单性使之成为求解Kohn-Sham 方程的高效方案之一。另外，赝势的引入可以保证计算中用较少的平面波数就可以获得较为可靠的结果。该方法具有较高的计算效率，使之日益发展成为有效的计算方法。本章首先对赝势平面波方法进行重点讨论，其次介绍了基于第一性原理计算软件一般步骤，最后结合Materials Studio 软件包应用，对锐钛矿型TiO 2(101)表面及其点缺陷结构进行建模和计算。

3.1 基本原理

基于密度泛函理论的第一性原理计算实质是求解Kohn-Sham 方程。实际求解Kohn-Sham 方程时，由于原子核产生的势场项在原子中心是发散的，波函数变化剧烈，需要采用大量的平面波展开，因而计算成本变得非常大，所以在计算中选取尽可能少的基函数。计算中选择的基函数与最终波函数较接近则收敛较快，当然包含的维度也应该尽量少。众所周知，根据研究对象不同，选择基函数的方法也不同的，如原子轨道线性组合法(LCAO-TB)、正交平面波法(OPW)、平面波赝势法(PW-PP)、缀加平面波法(APW)、格林函数法(KKR)、线性缀加平面波法(LAPW)、Muffin-tin 轨道线性组合法(LMTO)等，选取典型代表方法在随后的章节中重点展开讨论。与LAPW ，LMTO 等精度较高的第一性原理计算方法比较，平面波赝势法是计算量较少的方法，适用于计算精度要求不严格，因原胞较复杂而导致计算量陡增加的体系。为此，本章将重点学习赝势平面波方法，先学习电子能带的平面波基底展开以及赝势等相关基本概念，然后再讨论赝势引入原理。

3.1.1 平面波展开与截断能

1. 平面波展开

平面波是自由电子气的本征函数，由于金属中离子芯与类似的电子气有很小的作用，因此很自然的选择是用它描述简单金属的电子波函数。众所周知，最简单的正交、完备的函数集是平面波exp[())i k G r +?，这里G 是原胞的倒格矢。根据晶体的空间平移对称性，布洛赫(Bloch)定理(将在第4.1.1节中说明)证明，能带电子的波函数(,)r k ψ总是能够写成

(,)()exp()r k r ik r ψμ=? (3.1)

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

以用平面波基展开为

()()exp(())G

r G i K G r N ψμ=

+?Ω

(3.2)

式中1

N Ω其中Ω是原胞体积；这里G 是原胞的倒格矢，K 是第一Brillouin

区的波矢，()G μ是展开系数。Bloch 定理表明，在对真实系统的模拟中，由于电子数目的无限性，K 矢量的个数从原则上讲是无限的，每个K 矢量处的电子波函数都可以展开成离散的平面波基组形式，这种展开形式包含的平面波数量是无限多的。基于计算成本的考虑，实际计算中只能取有限个平面波数。采用的具体办法是一方面由于()r ψ随K 点的变化在K 点附近是可以忽略的，因此我们可以使用K 点取样通过有限个K 点进行计算。另一方面，为了得到对波函数的准确表示，G 矢量的个数也应该是无限的，但由于对有限个数的G 矢量求和已经能够达到足够的准确性，因此对G 的求和可以截断成有限的。给定一个截断能

22

()2cut

G K E m

+=

h (3.3) 对G 的求和可以限制在2()/2cut G K E +≤的范围内，即要求用于展开的波函数的能量小于cut E 。当0K =时，即在Γ点，有很大的计算优势，因为这时波函数的相因子是任意的，就可以取实的单粒子轨道波函数。这样，对Fourier 系数满足关系式*()()l l G G μμ-=，利用这一点，就可以节约不少的计算时间。 2. 截断能选取原则

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

ln tot

cut

dE d E ，E tot 是体系总能量，E cut 是截止能量。例如，当它的数

值小于0.01 eV/atom 时，计算就达到了良好的收敛精度，对于大多数计算0.1 eV/atom 就已足够。

3. 平面波基展开特征

用平面波基来展开电子波函数是因为用平面波基来计算有很多优点。平面波基能很方便地采用快速傅里叶变换(FFT)技术，使能量、力等的计算在实空间和倒空间快速转换，这样计算尽可能在方便的空间中进行。如前面讲到的哈密顿量中的动能项的矩阵元，在倒空间中只有对角元非零，就比实空间减少了工作量。第二，平面波基函数的具体形式并不依赖于核的坐标。这样，一方面，价电子对离子的作用力可以直接应用Hellman-Feymann 定理(将在3.1.5节中进行说明)得到解析的表达式，计算显得非常方便。另一方面也使总能量的计算在不同的原子构型下有基本相同的精度。此外平面波计算的收敛性和精确性比较容易控制，因

图3-1 截断能示意图

为通过截断能E cut 的选择可以方便地改变平面波基的多少。当然平面波基也有缺点，一般电子轨道具有一定的局域性，而平面波是空间均匀的，因此电子轨道展开时与原子轨道基相比，平面波基的个数要多得多。为了尽量减少平面波基的个数，一般在平面波的计算中都采用赝势(pseudopotentials)来描述离子实与价电子之间的相互作用，使电子轨道波函数在离子实内部的分布尽量平缓些。下面将讨论赝势概念及其引入思路。

3.1.2 赝势

1. 赝势引入

平面波函数作为展开基组具有很多优点，然而截断能的选取与具体材料体系密切相关。由于原子核与电子的库仑相互作用在靠近原子核附近具有奇异性，导致在原子核附近电子波函数将剧烈振荡。因此，需要选取较大的截断能量才能正确反映电子波函数在原子核附近的行为，这势必大大地增加计算量。另一方面，在真正反映分子或固体性质的原子间成键区域，其电子波函数较为平坦。基于这些特点，将固体看作价电子和离子实的集合体，离子实部分由原子核和紧密结合的芯电子组成，价电子波函数与离子实波函数满足正交化条件，由此发展出所谓的赝势方法。1959年，基于正交化平面波方法，Phillips 和Kleinman 提出了赝势的概念。基本思路是适当选取一平滑赝势，波函数用少数平面波展开，使计算出的能带结构与真实的接近。换句话说，使电子波函数在原子核附近表现更为平滑，而在一定范围以外又能正确反映真实波函数的特征，如图3-2所示。

所谓赝势，即在离子实内部用假想的势取代真实的势，求解波动方程时，能够保持能量本征值和离子实之间的区域的波函数的不变。原子周围的所有电子中，基本上仅有价电子具有化学活性，而相邻原子的存在和作用对芯电子状态影响不大。这样，对一个由许多原子组成的固体，坐标空间根据波函数的不同特点可分成两部分(假设存在某个截断距离c r )。(1)c r 以内的核区域，所谓的芯区。波函数由紧束缚的芯电子波函数组成，对周围其它原子是否存在不敏感，即与近邻的原子的波函数相互作用很小；(2)c r 以外的电子波函数(称为价电子波函数)承担周围其它原子的作用而变化明显。 2. 原子赝势

全电子DFT 理论处理价电子和芯电子时采取等同对待，而在赝势中离子芯电子是被冻结的，因此采用赝势计算固体或分子性质时认为芯电子是不参与化学成键的，在体系结构进行调整时也不涉及到离子的芯电子。在赝势近似中用较弱的赝势替代芯电子所受的强烈库仑势，得到较平缓的赝波函数，此时只需考虑价电子，在不影响计算精度情况下，可以大大降低体系相应的平面波截断能E cut ，从而降低计算量。图3-3为Si 原子赝势示意图。赝原子用于描述真实原子自身性质时是不正确的，但是它对原子-原子之间相互作用的描述是近似正确的。近似程度的好坏，取决于截断距离c r 的大小。c r 越大，赝波函数越平缓，与真实波函数的差别越大，近似带来的误差越大；反之，c r 越小，与真实波函数相等的部分就越多，近似引入的误差就越小。

可将真实价波函数(,)n r k ψ看作是由赝势波函数(,)n r k λ和内层波函数(,)J r k φ线性组合，即

(,)(,)()(,)n n nJ J J

r k r k k r k ψλσφ=-∑ (3.4)

其中系数()nJ k σ可由正交条件*

'(,)(,)0J n dr r k r k φψ=?确定，即

*()'(,)(,)nJ J n k dr r k r k σφλ=?

图3-2 赝波函数与势 图3-3 Si 原子赝势示意图

联合真实波函数(,)n r k ψ所满足的薛定谔方程：

[()](,)()(,)n n n T V r r k E k r k ψψ+=

可得到赝波函数满足如下方程

[](,)()(,)ps n n n T U r k E k r k λλ+=

(,)()(,)'(,')(',)ps n n R n U r k V r r k dr V r r r k λλλ=+? (3.5)

其中*(,')(',)[()](',)R J n J J J

V r r r k E k E r k φφ=-∑

ps U 称为原子赝势。根据密度泛函理论，原子赝势包括离子赝势ion ps U 和价电子库仑势和交换

关联势：()()ion ps ps ps H xc U U V r V r =++，其中后两项()ps H V r 和()xc V r 可以从真实电荷密度计算，

此时等于对应的全电子势()V r 和xc V 。

从上面可知，赝势应具有以下特征：(1)赝波函数和真实波函数具有完全相同的能量本征值()n E k ，这是赝势方法的重要特点；(2)赝势第二项是排斥势，与真实的吸引势有相消趋势，因此比真实势弱；(3)赝势包括局域项，其中非局域项同时与r 和'r 处的赝波函数(,)n r k λ和(',)n r k λ有关，而且依赖于能量本征值()n E k 。 3. 赝势分类

从上面的推导可以看出，赝势实际上是一种算子，但可以近似的将它处理成原子间距离的简单函数，叫做局部赝势。其结果是只要给定周期性排列的原子的位置，系统的能量就可以算出来。确定赝势的方法不是唯一的，主要有以下几种：

(1) 经验赝势方法，利用实验数据拟合有限几个V(K)的值，主要用途是在现代从头算原子赝势自洽迭代计算中作初始值使用。

(2) 模型赝势是半经验的，能用于自洽计算的原子赝势。

(3) 模守恒赝势是第一性原理从头算原子赝势，是核与芯电子联合产生的有效势，是从原子的薛定谔方程从头计算得到的，它可以给出价电子或类价电子的正确电荷分布，适于作自洽计算。

(4) 模非守恒赝势(超软赝势)，其优点是容易选择芯区的赝势波函数，减少了必须的平面波函数的数目，较大地减轻了计算工作量。

目前在第一性原理计算中应用较多的为模守恒赝势(Norm-conserving pseudopotentials ， NCPP)和超软赝势(Ultra Soft pseudopotentials ，USPP)两种方案。下面将分别讨论。

3.1.3 模守恒赝势

1. 模守恒赝势构造

在第一性原理计算中，常采用由Hamann D R 等提出的模守恒赝势。这种赝势所对应的波函数不仅与真实势对应的波函数具有相同的能量本征值，而且在c r (原子芯半径)以外与真实波函数的形状和振幅都相同(即模守恒)，且在c r 以内比较平缓。采用赝势计算关键在于可以有效的对化学键的价电子进行可再现的近似，赝势与全势在超过离子实半径以外具有完全相同的函数形式。这种赝势能生成正确的电荷密度，适合作自洽计算。

构造模守恒赝势基本思想是选择某个特定的电子排部状态(不一定就是基态)全部电子计算在一个孤立的原子中进行，从而得到原子价电子能量本征值和价电子波函数。选择一个离子赝势或赝波函数参数形式，通过对参数的调节，使得赝原子计算和全电子原子赝势计算采用相同的交换-相关势，在超过截止半径c r 后与价电子波函数形式相同，赝势的本征值等于价电子的本征值。如果电子波函数和赝势波函数满足正交归一，两者在截止半径以外的匹配性决定了模守恒条件自动成立。模守恒赝势要求赝势波函数满足：(1)本征值与真实本征值相等；(2)没有节点；(3)在原子核区以外(c r r >)与真实波函数相等；(4)在内层区(c r r <)内的赝电荷与真实电荷相等，将赝波函数插入到薛定谔方程中即得对应的赝势。一般说来，小的c r 移植性好，可用于不同环境，但平面波收敛慢。

第一性原理模守恒赝势可分为局域和非局域两部分，

()(,')ps loc v NL v v v v

U V r R U r R r R =-+--∑∑ (3.6)

其中v 是对离子势求和。考虑到原子球对称性，得用球谐函数将赝势的非局域 部分写成：

*

,,(,')(',')(,)(,')(,')NL lm lm l l l m l m

U r r Y Y V r r lm lm V r r θ?θ?==∑∑

如果将(,')l V r r 取成半局域形式，即径向是局域的，只有角部分是非局域的：

(,')()(')l l V r r V r r r δ=-，并定义角动量l 的投影算符l m

P lm =∑。则半局域的原子赝势可

以写成如下的形式：

()()()ps loc v l v l v

vl

U r V r R V r R P =-+-∑∑， (3.7)

为了简化计算，Kleinman 和Bylander(KB)将上面半局域赝势部分用一个非局域赝势来近似：

,(,')l m l m

NL l m

m l m

V V U r r V λλλλ=∑

(3.8)

2. 模守恒赝势应用特征

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

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

使用Optimised 或TM 赝势虽然能够把模守恒型赝势变的很“软”，但模守恒条件对于原本就已经没有节点价电子云分布的改造及最佳化的程度，与现今日渐普遍的超软赝势(它不必遵守模守恒条件)来比，节省计算的程度仍是有限的。总之，计算量的大小是取决于原子的种类这一点，是十分明确而普遍的认识，也就是说不同种类元素其势的“软硬”差异会令人明显感觉到。

3.1.4 超软赝势

1. 超软赝势构造

对于过渡族元素和第一周期元素，模守恒赝势不能明显降低所需平面波截断能E cut 。Vanderbilt 提出了超软赝势，其赝波函数在核心范围是被作成尽可能平滑，可以大幅度地减少截断能，即可使计算所需的平面波函数基组更少。就技术上而言，这是靠放宽模守恒的要求，采用广义的正交条件来达成的。为了重建整个总的电子密度，波函数平方所得到电荷密度必须在核心范围再附加额外的密度进去。这个电子云密度由此就被分成两个部分，第一部分是一个延伸在整个单位晶胞平滑部分，第二部分是一个局域化在核心区域的自旋部分。前面所提的附加部分是只出现在电子密度，并不在波函数。

超软赝势中总能量与采用其它赝势平面波方法时相同，非定域势NL V 表达如下

(0)

,I

I

NL nm n m

nm I

V D ββ=∑ (3.9) 式中投影算符β和系数(0)D 分别表征赝势和原子种类的差别，指数I 对应于一个原子位置。总能量用电子密度可以表示为：

2

()

,()[()()]I I

I

i nm i n m i i

nm I n r r Q r φββφ=+∑∑ (3.10)

式中φ是波函数，()Q r 是严格位于芯区的附加函数。超软赝势完全由定域部分，()ion

loc V r 和系

数(0)D ， ()Q r ，β确定。赝势是通过引进一系列正交条件来建立的，i j ij S φδ=，S 是哈密顿重叠算符，可以表示为1I

I

nm n m

S q ββ=+，系数q 是通过对()Q r 积分得到。从而，超软赝势的Kohn-Sham 方程可以写为：i i i H S φεφ=，H 可以表示为动能和定域势能之和

()

,I I I

eff nm

n m nm I

H T V D ββ=++∑ (3.11)

()(0)()()I I

nm nm eff nm D D drV r Q r =+? (3.12)

2. 应用特点

超软赝势产生算法保证了在预先选择的能量范围内会有良好的散射性质，这导致了赝势更好的转换性与精确性。超软赝势通常也借着把多套每个角动量通道当作价电子来处理浅的内层电子态。这也会使精确度跟转换性更加提升，虽然计算代价会比较高。与模守恒赝势对比，不同之处在于在超软赝势中存在重叠算符S ，波函数与D 有关。而且投影算符函数β数量要比模守恒赝势中大两倍多。与附加电荷相关的一系列计算可以在实空间中进行，这与函数中定域势的性质有关，而多余的步骤不会对计算效率产生较大的影响。在Laasonen 文献中提供了超软赝势计算的详细方法以及总能量微分表达式。

3.1.5 Hellmann-Feynman 力

Hellmann 和Feynman 在量子力学框架下给出了作用于离子实上(位置坐标为I R )的力I F 。离子受的力为总能对离子位置的偏导，

I I

E

F R ?=-

? (3.13) E 作为系统哈密顿量的能量本征值，满足Kohn-Sham 方程，

H E ??=

可以得到：

E H ?= (3.14)

将式(3.14)代入(3.13)得

I I I

H F E

R R ??????=--?? (3.15) 由于??是一个归一化常数，上式的第一项等于零。最终得到作用在离子上的力

I

H

F R ?

??=-? (3.16) 这就是著名的Hellmann-Feynman 定理

Hellmann-Feynman 定理计算出的力是和电子波函数相联系的，它的误差与波函数误差的一级修正量成正比，只有波函数非常接近真实的本征态时这个力才是精确的。所以在计算时需要同时考虑到离子弛豫和电荷密度自洽，即，离子在受力后到达一个新的位置，此时电子也需要接近瞬间基态，然后在新的离子位置和新的电子密度下进行计算，直至总能到达局部极小值。在得出离子受的受力后，需要对离子进行弛豫，即需要知道离子弛豫的方向和大小。 【练习与思考】

3-1. 查阅有关文献和书籍，找出3种以上不同的经验赝势，试分析他们的适用对象及特点？ 3-2. 查阅Martin R M 《Electronic Structure 》一书，写出模守恒赝势方法中模守恒的条件。 3-3. 查阅谢希德《固体能带理论》和有关书籍，试说明赝势平面波法采用了那些近似处理？

3.2 数值处理方法与技巧

3.2.1 超原胞方法

对于固体体系，目前主要用三种模型来模拟材料特性，一是团簇模型(Cluster Model)，二是嵌入模型(Embedded Cluster Model)，三是层状模型(Slab Model)。在以往的量子化学计算中，研究人员往往使用团簇模型来模拟固体材料，这种简化在化学上的确存在一定依据，因为根据化学家的直觉，一个分子体系作用是受局部相互作用支配的，从这一观点出发，可以用分子与原子簇的作用来反映分子与固体相互作用的性质。由于计算方法和计算条件的限制，没有考虑表面结构弛豫的影响，从而使计算的模型体系和实验体系存在较大的误差。第二种嵌入模型的提出主要是为了克服团簇模型在模拟材料表面时存在边界性问题，但该方法在计算过程中涉及到大量近似，需要针对不同的体系使用不同的计算方法，比如由Korringa J 、Kohn W 和Rostoker N 提出的格林函数(Green Function)方法和戴逊方程(Dyson Equation)，所以该模型并不为广泛接受。由此，人们提出了超原胞模型，该模型在模拟体系时采用了周期性边界条件，特别适合研究金属、半导体这类具有周期性的凝聚态体系。在后面的实例中基本上都是选取了超原胞模型。

层状模型(Slab Model)中的超原胞模型将体系看作沿晶格矢周期性排列的体系，在计算中所研究的原子都放在超原胞中，原子坐标或者其对应的周期性位置可用下式表示：

s a

x a

l l l l R R R τ=+r r r (3.17)

其中a

lsl τ为原子在超原胞中的坐标，s l 为原子种类的序号，a l 为多个同类原子之间的序号。R

r

为格矢。

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

图3-4 研究表面的分子吸附原胞

图3-5锐钛矿型TiO 2(101)表面原子层超原胞

[010]

[101]

[101]

Ti 5c O 3c Ti 6

O 3c ′ O 2

O 3c 〞

3.2.2 自洽电子弛豫方法

表示电子结构松弛有多种方法，其中密度混合方法是最有效的。该方法是使用在一个在固定位势之下将电子能量本征值的总和极小化，而不是将总能做自洽式的极小化。在步骤的最后新的电荷密度就会与初始的电荷密度来混合产生新的电荷密度再计算体系能量以重复叠代直到收敛为止。密度混合有以下几种形式：

(1) 线性混合，新的密度是上一次输入与输出密度的线性组合。 (2) Kerker 混合，可用下式表达：

2

22

max

()()[()()]new in out in G G G A G G G G ρρρρ=+-+ (3.18) 其中G 为截断波矢，A 为混合幅度。

(3) Pulay 形式，是最有效的一种方式，在Pulay 形式中新的输入密度上所有先前迭代步中密度的组合，它不仅与A 和max G 有关，而且和所有的先前迭代步有关。

(4) Broyden 混合形式，它与Pulay 形式有类似之处。

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

3.2.3 几何结构优化技术

对于给定各原子位置、元素种类的研究体系，通过密度泛函理论自洽求解Kohn-Sham 方程可以得到整个系统处于多电子基态时的总能。总能量对系统虚拟微位移的导数就是各原子的受力，即Hellmann-Feynman 力。这为我们理论预言物质的结构提供了一种行之有效的方法。因为自然界稳定的结构应该具有最低的总能，只要根据原子受力来调整原子的位置，直到整个体系的总能达到最低，所有原子受力为零。当然在实际的计算过程中，给出希望达到而且有限的计算精度，可找到能量面的(全局)最小值，这时所对应的物质结构就是自然界最稳定的结构，该过程被称为几何结构优化(Geometry Optimization)。为了确保搜索能量面的最小值时能找到全局最小而不是局域最小，并提高整个搜索过程的效率，需要一些强有力的搜索算法以使原子最快地运动到最稳定结构的位置。

能量最小化算法一般分为两大类：全局极小和局部极小。全局极小算法可以得到基态构型，如模拟退火和遗传算法；局部极小算法找寻的是亚稳态结构，最常用的方法有直接能量最小化、最陡下降法、共轭梯度法、Newton-Rapshon 方法、阻尼动力学方法等等。这里主要介绍材料模拟常用的三种优化方法。 1. 最陡下降法(Steepest Descent Method)

最陡下降法沿着局部净受力方向行走，以进行能量极小。从初始点开始，沿着局部梯度的反方向l g ，并通过在此方向上的一维极小化，移动到该方向的极小点，再从这个点，开始重复以上过程，直到达到所要求的精度。最速下降法在远离极小点效率很好，在接近极小时效率不高，而且沿梯度方向每前进一步将对接下来一步都引入一个正比于它梯度的误差，常常只在优化的最初几步使用这种方法。 2. 共轭梯度法(Conjugale-gradient Method)

共轭梯度方法克服了最陡下降法的困难。在此方法中，每相邻两次优化的起始点的l g 仍是正交的，但优化方向l v 由当前梯度l g 结合前一次优化方向1l v -和梯度1l g -共同决定，l v 和1l v -互为共轭：

1l l l l v g d γ-=+ (3.19)

其中γ是一个标量数，由前一次优化起点的梯度1l g -和优化终点的梯度(即当前时刻梯度)l g 共同决定，不同的算法给出各自的确定公式：

Flrtcher-Reeves 算法： 11

l l

l l l g g g g γ--?=?

Polak-Ribiete 算法： 11

11

()l l l l l l g g g g g γ-----?=

?

对于能量函数()f P ，可以按如下方式进行优化：若初始点在0P ，令000()v g f P =-=-?，沿0v 方向运用一维极小化方法到达该方向的一极小点1P ，则11()g f P

=?。由0g 和1g 可得到1γ，则1100v g v γ=-+，再沿1v 找极小，重复以上过程，如果函数是含N 个变量的二次型，则通过N 次一维极小化就可以找到极小。上面两种方法在优化中只用了势能函数的一阶导数，即梯度。

3. Newton-Rapshon 法

任何给定点的能量都可以展成泰勒展式：

21

()()'()''()() (2)

U x x U x U x x U x x δδδ+=+++ (3.20)

'()U x 是在x 处的一阶导数矢量，''()U x 为二阶导数矩阵，称为Hessian 矩阵。对泰勒展式只

取到二阶导数，而忽略高阶的，则位移矢量x ?为：

1x H g -?=- (3.21)

其中''(),'()H U x g U x ==

一般情况下，可以重复以上过程直到能量最小，这就是所谓的Newton-Raphson 方法。Newton-Raphson 方法收敛速度很快，但是不能保证收敛的方向，而且需要计算Hessian 矩阵，计算量很大。为了避免函数值上升，采用以下修正公式：

1x H g α-?=- (3.22)

其中α通过一维搜索算法确定，保证函数值向最小方向收敛。

为了避免频繁计算Hessian 矩阵，常采用更新校正方案(updating scheme)，常见算法有：Davidon —Fletcher —Powell(DFP)和Broyden —Fletcher —Goldfarb —Shanno(BFGS)。具体表达式如下：

1

()()

DFP DFP DFP DFP i i i i

DFP i H g H g x x H H

x g g H g +????????=+-

??????? (3.23) 1

()()[]BFGS BFGS BFGS

BFGS BFGS i i i i

i BFGS

i H g H g x x H

H

g H g v v x g g H g

+????????=+-+???????????? (3.24) 其中BFGS i BFGS i H g x

v x g g H g

???=-

???????。 3.2.4 快速傅立叶变换(FFT)

平面波函数中使用的赝势常常是非局域化的。这意味着，为了要模拟全电子原子的散射性质，必须要在每个角动量信道里采用不同的散射位势。然而这些位势只在原子核心区域之内有所不同，采用实空间的表示方法是有效的。目前一般的计算软件采用的是选择在倒易空间中使用赝势，也就是说在倒易空间中进行对波函数和势操作的加总。这对于以平面波基底系数来表达的波函数是一种很常用的做法。

然而对于大晶胞的Hamiltonian 波函数和势能项进行实空间计算，会变得比在倒易空间中求值更快。这是因为赝势的非局域部份只有在核心内的区域才不是零，并且在远比各个单一原子核心区域还要大的晶胞中几乎到处都是零。因此，在实空间求值波函数和势的积，并且使用傅立叶转换来获得在倒易空间的值就变的更为有效率。这就意味着波函数和各种势在计算中需要一个有效的途径完成从实空间和倒易空间的相互转换，傅立叶变换为此提供了有力的工具。

平面波赝势方法中采用一组有限个正交的平面波作为基矢，这就将久期方程中的矩阵元计算：ij i i H k G H K G =++，归结为傅立叶变换。

快速傅立叶变换在计算中频繁使用。傅立叶变换网格不仅影响计算精度，也影响计算速度，因此也存在计算结果对傅立叶变换网格的依赖关系，严格地说，这也是一个收敛问题。傅立叶变换网格空间稀疏，将导致体系不易收敛或收敛结果粗糙，网格越密，则将大量消耗计算资源。

3.2.5 k 空间取样规则

在通常的能带图中，经常会出现比如Γ、Λ、Κ、Χ ……等等的符号。这些符号表示的是布里渊区内具有高对称性的一些特殊的k 点，这些k 点有特殊重要的意义。在进行体系总能计算时，通常要对布里渊区内的波函数或本征值进行积分，在实际计算过程中，积分是通过对部分特殊选取的k 点求和完成的。比较常见的k 点网格方法有Monkhorst-Pack 方法和四面体网格。

根据Bloch 定理，周期体系中的电子波函数可以表述为调幅平面波的形式，即(,)()exp()r k r ik r ψμ=?。其本征能量和本征矢量为()n E k 和()n k ?。不同的电子状态按照量子

数k 进行分类，而量子数n 则表征能态的分立性。研究多体体系的价电子问题，归根结底是计算出不同类的电子状态的本征值和本征矢量，体系处于基态情况下，哪些不同k 的低能量状态被电子占据。因体系具有周期性，所以，第一布里渊区的所有k 可以代表所有的k 。但

是，由于周期边界条件确定的k有无穷多个，而且计及相互作用势的实际体系中许可的k点在倒易空间内是不均匀的，实际计算过程中只能选取有限个点。第四章还将对k点的选取方法进行深入分析。

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

3.2.6基于密度泛函理论的第一性原理计算框架及步骤

结合第2和3章相关计算理论和近似方法，本节总结出第一性原理计算框架中的关键步骤及典型实现方法，以及对应的波函数变化形式，如表3-1所示。

表3-1 基于密度泛函理论的第一性原理计算框架及步骤

几何优化优化原子位置steepest descent

BFGS

DIIS

收敛参数电荷密度体系总能原子受力

【练习与思考】

3-4.快速傅立叶变换(FFT)对体系哈密顿量进行数值化计算是如何实现的？

3-5.写出超晶胞的周期性边界条件。并指出采用周期性结构计算的有优点。

3.3基于Materials Studio的实践流程

3.3.1软件介绍及计算步骤

Materials Studio是Accelrys公司专为材料科学领域开发的可运行于PC机上的新一代材料计算软件，可帮助研究人员解决当今化学及材料研究中的许多重要问题。Materials Studio 软件采用Client/Server结构(图3-6所示)，客户端可以是Windows 2000、XP或NT系统，计算服务器可以是本机的Windows 2000、XP或NT，也可以是网络上的Windows XP、Windows NT、Linux或UNIX系统。因此使得任何的材料研究人员可以轻易获得与世界一流研究机构相一致的材料模拟能力。

图3-6 Materials Studio软件采用Client/Server结构

Materials Studio系列软件有一个核心模块，称为Visualizer。它是一个搭建分子、晶体及高分子结构模型的图形界面，同时提供各个计算模块的基本环境和分析工具，如图3-7所示。在Visualizer中提供了一个完整且容易使用的接口来操作CASTEP等计算程序，使得Visualizer所提供的交互式图形接口与仿真环境有效发挥CASTEP等程序的量子力学功能。下面TiO2表面模型构建都是在该模块下面完成。

图3-7 Visualizer 图形界面

Materials Studio 软件中计算的模块很多，这里以Materials Studio 系列软件中4.0版的模块CASTEP 和Dmol 3为例进行介绍。

CASTEP 模块是Materials studio 系列软件中运用赝势平面波方法的模块之一。在Castep 中，模型的建立基于超晶胞结构(supercell)的方法，并且使用周期性边界条件。由于周期边界条件的使用，晶体的表面用一个有限长度的模型(slab) 就可以建立起来。周期性边界条件和Bloch 原理联系在一起，即在一个周期性结构中，每个电子波函数都可以表示成和单胞中的波函数类似的形式，()()ik R n n r e r ψψ?=，而,()iG R n n G n

r C e ψ?=∑，这样每个电子波函数都可以写成平面波叠加的形式，即exp[()]i k G R +?。有关更为详细的说明请看3.3.3节中的操作方法。

Dmol 3使用原子中心网格的数值函数作为其原子基，原子的基函数通过解不同原子的DFT 方程得到，并将其储存为一系列的三次样条函数。这种基组是十分精确的，高精度的基组减少了重叠效应，于是体系可以得到准确的描述，能表示正确的电荷分布，因而可以提高对分子极化描述的精度。在Dmol 3中，电子密度依据以原子为中心的多极部分密度展开，这提供了一种简洁而精确的表示密度的方法，通过求解Posson 方程，可以用电子密度的多极表示来估算库仑势，从而将库仑势计算这一较快的步骤所代替，这一操作使计算所用的时间与体系的大小成线性正比。Hamilton 矩阵元素用复杂的数值积分算法来计算，其计算量也接近与体系大小呈线性相关。此外，Dmol 3的计算法允许对数值积分过程进行高效的并行处理。无论是对于分子体系或周期性结构，几何结构和过渡态的优化都采用非定域的内坐标。这包括在进行内坐标优化时，对其施加笛卡尔几何限制的能力。Dmol 3采用了一种新的过渡态搜索方法，将传统的LST/QST 算法与共轭梯度方法(conjugate gradient methods)结合使用。这种新的快速的方法允许在进行过渡态优化时不计第二阶导数矩阵，而这是传统过渡态搜索算法所必需的。

计算步骤可以概括为三步：首先建立目标物质的几何结构；其次对建立的结构进行优化，

这包括体系电子能量的最小化和几何结构稳定化；最后是计算要求的性质，如电子密度分布(Electron density distribution)，能带结构(Band structure)、状态密度分布(Density of states)、声子能谱(Phonon spectrum)、声子状态密度分布(DOS of phonon)，轨道群分布(Orbital populations)以及光学性质(Optical properties)等。第一性原理计算的一般步骤如图3-8所示，下面将依照该步骤分别说明其细节及相关选项设置。

图3-8 Materials Studio软件计算步骤方框图

3.3.2几何结构的构建

本节以锐钛矿型TiO2的解理表面(101)面的构建为例，来说明构建材料的计算模型结构的具体方法。

1.构建体相晶体结构

采用Materials Visualizer模块来建立晶体表面结构。先建立材料的体相晶体结构，在此基础上沿某一特定方向进行剪切，从而获得所需要的理想表面。获得体相晶体结构有两种途径：一是直接从软件的结构库(structure)中导入。具体操作：点击File | import…命令，显示出Import Document(导入文档)对话框，注意，这个对话框也可以通过使用标准工具条中的导入按钮打开，在导入文档对话框中选择structure文件夹中需要的材料的晶体结构。在该文件夹包含了常见的半导体、陶瓷、金属、矿物以及聚合物等结构。

另一种方式是已知晶体结构数据，如空间群，内坐标等，通过Materials Visualizer模块中的建模工具建立晶体结构，即创建3D文档。具体操作：选择→…命令，在New Document 对话框中选择3D Atomistic，在工程管理器中显示创建了一个3D Atomistic Document.xsd的文件。在文件名上右击，选择改名命令，可以键入新的文件名。选择命令保存文件。

选择Build→Crystals→Bu ild Crystal…命令(如图3-9)，打开建立晶体(Build Crystal)对话框(如图3-10)。

图3-9 Build菜单图3-10 Build Crystal对话框

在空间群(Space Group)选项卡，输入空间群(Enter group)文本框内输入空间群，也可以在下拉列表中选择空间群或者直接输入空间群序号。在晶格参数(Lattice Parameters)中输入晶体的a，b，c轴的晶格参数。例如，锐钛矿型TiO2为四方晶格，空间群为I41/amd，晶格参数a=3.782?，c=9.502 ?。一旦空间群输入后，其α，β，γ的晶格参数自动设置。按Build 按钮，一个定义了晶格参数的空晶格就出现在三维窗口中。

此后，由于对称性已经产生，只需要在相应位置加入原子，对称位置的原子自动产生。选择“添加原子(Build Add Atoms)”命令，打开“增加原子(Add Atoms)”对话框(如图3-11)。也可以通过单击“增加原子(Add Atoms)”按钮打开这个对话框。在“选项(O ptions)”部分，不要选择“原子成键测试(Test for bonds as atoms are created)”。在此选项选定时，Materials Studio 会在建立晶体过程中自动加入键。Materials Studio也有一个灵活的“键计算(Bond Calculation)”工具，可以选择、编辑、定义键参数。一般使用自动选项就行了。仍然是在“选项(Options)”部分，把“坐标系统(Coordinate system)”设为”分数(Fractional)”。回到“原子(Atoms)”选项页，元素列表中选择相应的Ti，键入以下a，b和c的值a =0，b = 0，c = 0。按添加(Add)按钮，一个相应的原子和它的对称位置原子就加入晶胞中了。在“原子Atoms”选项页，选择O元素，键入下列数值a = 0，b = 0，c = 0.208。单击对话框底部“添加(Add)”按钮，一个氧原子和它对称位置的原子就被加入了。关闭“添加原子(Add Atoms)”对话框。

为了提高结构的可视化，可在Display Style对话框中的原子(Atom)选项中将“棍(Stick)”改为“球棍(Ball an d stick)”选项(如图3-12)。最简单的操作过程为：右键视图的空白区域，从下拉菜单中选择Display Style即可。

图3-11 Add atom对话框图3-12 Display Style对话框

2.构建晶体表面

构建晶体表面主要由以下四个步骤组成：

步骤一：获得体相晶体结构。前面已经进行了具体介绍。

步骤二：在上面构建的体相晶体结构的基础上，选择Build→Cleave surface，然后在Cleave plane对话框里设定需要解理的晶面指数，如图3-13所示。例如锐钛矿型TiO2的解理表面是(101)面。

步骤三：选择layer的厚度，就在Depth里选择Fraction，也就是重复的晶面数，注意这里不是周期数。在Cleave之后，如果需要修改这个Depth值，可以根据厚度来Recleave。

Top的设置也很重要，它直接决定了剪切的位置，即暴露在表面的是那类原子。这

个Top值究竟设置多少需要研究人员十分谨慎的选择，否则无法获得所需要的真

正的实际表面。

步骤四：同样在Build→Crystal选项里选择建立真空层(Build Vacuum Slab)对话框。如图图3-14，在该对话框中设定真空厚度。在这里，可以设定的厚度不能太小。一般情况下最好大于4 ?。点击建立(Build)按钮，修改一下显示风格(Display style)，获得图3-15显示的结构。

图3-13 Cleave Surface对话框图3-14 Build Vacuum Slab Crystal对

话框图3-15 锐钛矿型

TiO2 (101)面

3.数据的导出

事实上前面的晶体结构坐标数据可以方便导出，用于其它软件，如Wien2K，VASP等的输入文件。利用 菜单命令将晶体结构数据导成cif或者pdb文件。这里以pdb文件导出为例举例说明。以下是导出的锐钛矿型TiO2(101)数据。从该文件中我们可以看到晶体表面的晶格常数以及原子坐标数据。

|-----------------------------------------------------------------------------| REMARK Materials Studio PDB file

REMARK Created: Sat Jun 09 06:45:50 中国标准时间2008

CRYST1 5.443 3.776 14.353 90.00 90.00 110.30 P1

ATOM 1 TI MOL 2 -0.000 -0.000 -0.000 1.00 0.00 Ti

...

ATOM 6 TI MOL 2 1.888 -0.000 -2.372 1.00 0.00 Ti

ATOM 7 O MOL 2 -5.664 1.888 -2.770 1.00 0.00 O

...

ATOM 18 O MOL 2 -0.000 -0.000 -1.973 1.00 0.00 O

TER

|-----------------------------------------------------------------------------|

3.3.3计算操作方法

模型构建好之后，然后启动CASTEP程序进行几何优化和性质的计算，

1.计算任务的设置

在CASTEP软件中执行任务设置，主要是通过Visualizer应用窗口中的工具条之一“计算(Calculation)”来进行。可以更改工具框中的相应选项，来配置诸如：“电子选项(Electronic Options)”、“结构优化选项(Geometry Optimization)”、“属性选项(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 SCF tolerance (eV/atom) 1.0e-5 2.0e-6 1.0e-6 5.0e-7

k-points separation (1/?) 0.07 0.05 0.04 0.04

Energy cutoff Coarse Medium Fine Ultra-Fine

?交换-关联函数的设置

CASTEP 提供了两种交换-关联函数的设置，一种是局域密度近似(LDA)，它使用的是

CA-PZ形式的赝势，是引用Ceperley D M，Alder B和Perdew J P，Zunger采用的赝势形式。

另外一种是广义梯度近似(GGA)，它有三种形式可供选择，分别是PBE(Perdew-Burke-Ernzerhof functional)、RPBE(Revised Perdew-Burke-Ernzerhof functional)、

PW91 (Perdew-Wang generalized-gradient approximation)，对应的是三种不同的广义梯度近似

产生形式。以上两种交换-关联函数LDA和GGA及其对应的可选形式是通过CASTEP “Calculation”中的“Setup”选项来进行统一配置，如图3-16和3-17所示，改变其中“Functional”

的类型和每一类型对应的函数产生形式，即可得到相应的设置。交换关联函数设置对计算结果的影响后面将加以测试说明。

图3-16 LDA设置对话框图3-17 GGA设置对话框在范例中，在计算六种表面原子终止结构的能量时分别采用了局域密度近似(LDA)中

CA-PZ方案和广义梯度近似(GGA)中PW-91方案，在研究表面原子弛豫和电子结构时采用了GGA-PW91方案。

?电子自旋选项设置

该设置一般是在研究电子自旋的磁性材料中才考虑，如果是不考虑材料的磁学性能，可

以不做选择(默认设置是不考虑电子自旋)。如果要对某些原子的自旋属性进行设置，有两种

方法。方法一：选中要设置的原子，在Modify的下拉菜单中选择“Electronic Configuration”，

在出现的对话框中可以修改原子的“Formal Charge”和“Formal spin”等。Hubbard U选项

中，可以对LDA+U中的“U”进行设置，这仅仅在考虑电子自旋情况的时候有效。方法二：

选中要设置的原子，在“属性窗口”(Properties)中直接选择“Formal Charge”、“Formal spin”

和“Formal spin direction”等，双击鼠标后，可以在对话框中设置与电子自旋有关的相应参数。

(2)设置电子选项

在利用CASTEP 做有关单点能量、结构优化、动力学、弹性常数、过渡态等计算时，必须对电子选项进行设置。在电子选项中主要有以下方面的设置。

?截断能的设置

CASTEP中分子轨道是通过平面波基来扩展。平面波基的数目是通过截断能的高低来控制，选择的截断能过低会影响计算结果的正确性，而选择的截断能过高会使计算量过大。因此，在计算中要选取合适的截断能。一般情况下，在计算前可以选取几个截断能来测试，哪个更合适，在保证计算精度的前提下选择尽可能低的截断能。因此，截断能是CASTEP计算中最重要的参数之一。

设置截断能的最简单方式是在“CASTEP Calculation→Electronic”选项中选择“More”，然后在对话框中选择“基本”(Basi c)选项，在“能量截断”(Energy cut-off)中填入设置的值(图3-18)。计算中截断能(E cut)取为380 eV。

图3-18 截断能的设置对话框

?赝势的设置

在CASTEP中有Ultrasoft赝势和Norm-conserving赝势这两种赝势，它们有各自的优势，Norm-conserving赝势一般适用于金属体系，对于纳米碳管的研究，采用低截断能的快速计算，使用Ultrasoft赝势更适合。设置的方法是在“CASTEP Calculation→Electronic”选项中，在“Pseudopotential”的下拉框来进行选择。如图3-19所示。在描述离子实与价电子之间的相互作用时采用超软赝势(ultrasoft pseudopotentials，USPP)来，尽量减少平面波基的个数。为了减少计算需要的时间，采用了超软赝势研究锐钛矿型TiO2(101)表面结构。