山东大学博士学位论文 铁电体的第一性原理研究
第二章 研究方法与程序介绍
§2.1 全电子法和赝势法
应用于铁电体的第一性原理计算方法和工具很多,根据对势函数及内层电子
的处理方法不同主要分为两大类,一种是波函数中包含了高能态和内层电子,而
势函数只是原子核的贡献,这称为全电子(all electron calculation)法,另一种处
理方法是势函数为原子核和内层电子联合产生的势,称为离子赝势,波函数只是
高能态电子的函数,这称为赝势(pseudo-potential)法。
因为内层电子对价电子的排斥作用部分地抵消了原子核对价电子的强吸引
作用,所以赝势是一种比较弱和比较平坦的势。引入赝势的要点在于,赝势对应
的薛定谔方程与真实势对应的薛定谔方程有相同的能量本征值。在这一前提下,
引入赝势的方法不是唯一的。在第一性原理计算中,用的是所谓模守恒赝势法。
这种赝势所对应的波函数有一个特点,在离开原子核一定距离的空间,它与真实
势对应的波函数不但形式相同,而且幅度相等,故称模守恒。这种方法从原子势
算起,不引入任何实验参数,所以又称为从头算起(ab initio)赝势方法。
一般来说,赝势法计算量较小,但其中消去了内层电子态,相对于全电子法
多引入了一个近似。该方法的优点是较便于计算离子受到的作用力,后者等于总
能量对原胞内离子位矢导数的负数,称为 Hellmann-Feynman 力。赝势法用平面
波展开来表示价电子态,如果晶体中的原子有 2p 未满壳层(如氧)或 3d 未满壳
层(如钛),则赝势将很“硬”,为满足模守恒,需要为数很多的平面波基函数,
计算量太大。为此发展了超软赝势(ultro-soft pseudo-potential)法。对波函数引
入一个重叠算符,使赝势变软,减少了平面波基函数。在铁电体研究中用的赝势
法通常是这种方法。
全电子法表示电子态时将空间分为两部分:一是原子核附近的球形区,称为
丸盒(muffin-tin)区,二是原子核间的其它区域。在球形区,基函数、电荷密度
16
第二章 研 究 方 法 与 程 序 介 绍
和势均用径向函数展开,在其它区域,这些量用平面波或球面波展开。若用平面
波展开,则称为线性缀加平面波(LAPW)法;若用球面波展开,则称为线性丸
盒轨道(LMTO)法。全电子法计算所得的总能量较赝势法的要高,因为它包含
了赝势法中未计入的内层电子的能量,但它能得出精确的能量差,后者是十分重
要的。这种方法的缺点是较难计算离子受到的作用力。
§2.2 线性缀加平面波方法
固体能带理论以确定固体电子能级为主要任务,是深入研究固体物理的理
论基础。固体的许多基本物理性质,如电导率、热导率、磁有序、光学介电函
数、振动谱等,原则上都可由固体的能带 理论加以阐明和解释,或需要对具体
材料能带结构的了解作为基础。近年来,由于计算机技术日新月异的发展,能
带理论在方法上的创新得以在计算中实现。实际上,能带理论已成为计算物理
的一个重要组成部分。
原则上,从量子力学的基本方程薛定谔方程出发,可以计算任何分子或晶
体的性质。但是,由于所有电子彼此之间存在相互作用,数学上严格求解这种
多电子问题是不可 能的。解决这一困难需 采用不同的近似法,单电子近似正是
其中的一种,在分子量子化学和晶体的能带理论中它属于首选方法。
单电子近似的概念,在量子力学的前身玻尔理论中早就存在。单电子模型
的基础是假设体系的所有核以及所有其余电子对某一给定电子的作用可近似
地用一平均“有效”场的作用来代替(哈特利—福克近似,更为严格精确的处
理是后来发展的密度泛函理论)。在这个场中,电子的势能即所谓“有效单电
子势能”V
=
r
)(rV
,而玻尔原子模型中的单电子态就是相应于量子力学中的单
电子波函数ψ=
)(rrψ
,也即单电子轨道。在求得体系单电子轨道序列{ψn}以及
相应的单电子能级序列{εn}或体系的单电子“能谱”之后,就获得了该体系的
电子结构和电子性质各个方面的信息。
17
山东大学博士学位论文 铁电体的第一性原理研究
在量子力学中,单电子轨道与能谱是由单电子的定态薛定谔方程确定的,
采用原子单位后,该方程具有如下的形式:
ψ+∇−
2
1
2
V
εψψ=
. (2.1)
实际上,在将单电子近似法用于原子、分子和固体时,就会依次引出著名的原
子壳层理论、分子量子化学中的分子轨道法以及固体的能带理论。
原子壳层理论是单电子模型最简单的例子,它常当作是原子结构的严格描
述。假设单电子所在的有效场是球对称的V )(r =v
)(rV
(即为辏(音 còu)力场),则
严 格 的 量 子 力 学 计 算 表 明 , 原 子 中 电 子 的 定 态 相 当 于 单 电 子 能 级 组
1s,2s,2p,3s...,每一个能级相应于适当的单电子波函数或称为原子轨道(AO)。
单电子能级对磁量子数是简并的(这是辏力场的性质),但对轨道量子数却是非
简并的(因为这里的有效场不是库仑场)。在壳层模型中,多电子原子的状态由
该原子的电子组态所描述,它表明了被占据的单电子能级和其上的电子数(即
能级布居数)。
与单原子不同的是,将单电子模型应用于多原子体系时,由于势能表达
式中有多个吸引中心(原子核或原子实)存在,方程(2.1)的求解将变得更为复
杂。可以设想,在任何给定原子的周围,整个体系的势能多少接近于该原子
单独产生的势能,而在这个范围内,体系薛定谔方程的解也接近于此单原子
的方程之解,即接近它的原子轨道 AO。从数学上考虑,此时单电子波函数
的一种简便形式就是体系中所有原子的 AO 的线性组合。通过将有效单电子
哈密顿算符的本征函数展开为 AO 之和来求解多原子体系薛定谔方程,这种
方法即为原子轨道线性组合法(LCAO)。在分子体系中,它被称为分子轨道的
原子轨道线性组合法(MO LCAO);在晶体的情况,此即紧束缚的单电子近似
法。
固体能带理论在单电子近似的基础上,进一步通过将固体抽象为具有平移
对称性的理想晶体,借助布洛赫定理,将能带问题归结为单电子在周期性势场
18
第二章 研 究 方 法 与 程 序 介 绍
中的运动。同时在密度泛 函理论之 下,单 电子薛定谔方程(2.1)演化为著名的
Kohn-Sham 方程。LCAO 法有可能将晶体的性质与原子的性质以及原子的相互
作用联系起来,但如果这不是所期望的目的,则还可以用平面波或是发散 波形
式的“基”函数与 AO 波函数一起作展开基,甚或完全替代后者。就精确度而
言,这比 LCAO 法更具优点,由此发展了另外两大类能带结构计算方法:平面
波法和散射法 [3,4](或称格林函数法,又常称为 KKR 方法)。此外,根据研究对
象的物理性质对晶体势作合理有效的近似处理,比如采用一个赝势替代原子核
与芯态电子对价电子的作用,而集中地把计算重点放在价态和类价态电子上,
这样的思路即为另一类能带计算方法赝势法。
I
I
II
图 2.1 Muffin-tin:将原胞分为原子球区(I)和球间区(II)
平面波法和格林函数方法有一个共同的特点,即采用 J.C. Slater 提出的所
谓的 Muffin-tin 势[5]。它是这种处理实际上是基于早先原胞法的思路,即以一
个 Wigner-Seitz 原胞中电子的能量和波函数为出发点,将晶体的波函数用原胞
中电子波函数为基函数展开,并通过晶体波函数在原胞边界上所满足的条件来
确定晶体波函数。Muffin-tin 处理进一步将晶体的原胞空间分为两类区域:以
每个原子为中心的 Muffin-tin 球形区 I 和球外的球间区 II,如图 2.1 所示。在
球形区取球对称势,在球间区外则取常数势,通常选取适当的能量零点,使此
常数为零。Muffin-tin 球半径的选取要求球落在原胞之内,各原子的 Muffin-tin
球半径可以互不相等,但彼此不能相交。显然,和原胞法中选取的势场相比,
它更接近实际情况,而且也避免了原胞法中要满足边界条件的困难。同时,可
以看到,Muffin-tin 势可以方便地推广到更加复杂的格子,即原胞中不仅含有
一个原子的情况。这时可以按各个原子为中心做各自的原子球,半径可以不等,
19
山东大学博士学位论文 铁电体的第一性原理研究
但要互不相交。球内都有球对称势,球外的势场为零。所以这个模型更灵活,
适用性强。即使对于晶体势场不能完全用 Muffin-tin 势描述的情况,也可以按
微扰的思想进行考虑。
基于 Muffin-tin 势是选取,可以建立起一套缀加平面波(augmented plane
wave APW)。在区域 I 中有球对称势,Kohn-Sham 方程的解应有如下形式:
(2.2)
=ρ
)(
)ˆ(
ρ
)
ρ
,
Y
lm
ER
l
(
ϕ
lm
式中 是以原子为中心的矢径 ρ 的角度部分;Y
ρˆ
)ˆ(ρlm
是球谐函数;R
( ρEl
)
,
是径
向波函数,它满足如下径向 Kohn-Sham 方程:
−
d
1
d
2
ρρ
2
ρ
dR
l
d
ρ
+
ll
(
+
2
ρ
)1
+
V
ν
)
(
ρ
ER
l
(
,'
)
ρ
=
RE
l
'
( ρE
)
,'
(2.3)
式中
)ρνV 是第ν个球内的球对称势,l 为角量子数,在第ν个球内,缀加平面波
(
函数可定义为各个分波函数的线性组合:
φν
=ρ
)(
l
+
∑ ∑∞
lm
−=
l
=
0
Ya
lm
lm
)ˆ(
R
ρ
l
( ρE
)
,'
(2.4)
在球外,势场为零,解应有平面波形式。设第ν个球心的位矢为 ,则νr
所以在第ν个球外,
=r
ρν +r
,
e
i
k
⋅
r
ν
e
i
k
ρ
⋅
e
⋅ =rki
而后一个因子可以按球谐函数展开:
∑ ∑∞
4
e
ik
+
l
⋅ =
πρ
l
=
0
i
lm
−=
l
kYkj
(
)(
l
)
ρ
∗
lm
(lmY ρ
)
(2.5)
其中 , 分别表示 k , ρ 的角度部分。
ρˆ
kˆ
( ρkjl
)
为l 级球贝塞耳函数。在球表面
上 ρρ=
ν
,根据波函数连续的条件,可以得到关于系数 a 的方程:
lm
a
lm
=
i
rk
e
4
⋅
π ν
kjkYi
l
(
)ˆ(
∗
lm
l
/)
ρ
R
l
(
,' νρE
)
(2.6)
于是,缀加平面波函数可以写为:
20
第二章 研 究 方 法 与 程 序 介 绍
l
kj
(
l
)ˆ(
YkYi
l
)
ρ
lm
∗
lm
)ˆ(
ρ
ER
l
(
,'
/)
ρ
R
l
(
E
,' νρ
)
(2.7)
∞
+
l
r
ν
⋅ ∑ ∑
i
lm
−=
=
0
l
ik
e
4
π
(
φ
rk,
)
=
其中
r =
e
r +ν
i
rk
⋅
ρ
。
可以看到,在原胞法中原胞边界上存在的难处,用缀加平面波基函数以后
就避免了,各个原子附近的球对称势场所句定的波函数是借助于球间的平面波
形式的解来相互连接的。这里并没有要求波函数在球面上连续,不过这一不连
续性只是在单个基函数中存在,由这些基函数线性组合构成的晶体波函数仍可
以是一个连续性并有连续导数的函数。缀加平面波函数是基于 Muffin-tin 势建
立起来的的一套函数,它将作为缀加平面波法的基函数使用。
在两个在球间区,由于势场为零,Kohn-Sham 方程的解应有平面波的形式。可
以在球间区采用平面波形式的基函数,它“缀加”在球形区其他形式的基函数
(径向函数与球谐函数的乘积)上一起组成了缀加平面波方法(APW)的展开基矢
集。随着计算机的发展,对于各种能带理论方法,原则上对一个给定的周期性
势 场 , 可 以 将 单 电 子 问 题 求 解 到 任 何 需 要 的 精 度 , 但 实 际 上 , 对 同 样 的
Muffin-tin 势,用不同的能带方法,计算结果会有一定的差别。随着局域交换
关联势理论的发展和应用,通过自恰能带运算,可以得到晶体基态的性质,可
以得到晶体基态的性质,如结合能、原子间距,电荷转移、磁矩等,还可以得
到单电子的激发谱。但是,APW 方法存在着以下三个计算及概念上的困难:
(i)在建立 Muffin-tin 球的径向解时隐含了能量,它们的矩阵元都是能量的函数,
使久期方程成为能量的超越方程,需要另外用自洽计算来求解能量本征值。(ii)
当径向解的节点落在 Muffin-tin 球面上时,久期方程出现奇异性。(iii)基函数
的导数在球面上不连续。引入线性化的概念[6,7]可以解决上述各问题。目前有
两种十分有效的线性化方法:一是选取一套与能量无关的 Muffin-tin 轨道,使
得矩阵元中不含能量,此即线性化的 Muffin-tin 轨道方法(LMTO);一是在
Muffin-tin 球内给 APW 基函数增加一对能量求导的项,使得球内径向函数的解
21
山东大学博士学位论文 铁电体的第一性原理研究
不再是能量本征值的函数,而代之以某个待选定的能量参数值,这就是线性缀
加平面波方法(LAPW)。
在密度泛函理论内,LAPW 方法是迄今为止已经发展起来的许多固体材料
能带计算方法中最为有效和最为精确的方法之一。LAPW 方法可以用来计算固
体材料的体相、表面和界面的电子态,归纳起来,它能够给出以下有关电子态
的基本物理信息:(1)能带结构:电子态能量与波函数;(2)电荷密度:总电子
密度和各电子态电荷密度;(3)态密度:总态密度、分波态密度及分区态密度;
(4)总能量。
用 LAPW 方法可以更加方便地引入非 Muffin-tin 效应,假定在非 Muffin-tin
势 VMT(r)中增加了两项修正:
V(r)=VMT(r)+V1(r)+VNS(r) (2.8)
其中 V1(r)为球间区域内的修正,它在球内为零;VNS(r)是球内势的非对称部分,
在哈密顿量中的能量矩阵元中多了两项:
H
ij
=
H
TM
ij
+
V
+1
ij
ijV
NS
(2.9)
这里V 只在球间区域内才不为零,此时基函数有平面波形式,平面波实际上
NS
ij
足够多,完全可以用适当的变分自由度来求解球间区内的修正。球内的修正常
可以用球谐函数展开,计算时更加方便,
V
NS
)r(
∑=
LM
υ
LM
r
Y
)(
)ˆ(
LM r
(2.10)
采用 Muffin-tin 势场的 LAPW 方法对密堆积结构的金属是很好的近似,但
对于非密堆积的开结构,如半导体、清洁的或有吸附物的表面、薄膜、空洞材
料、准低维材料等,用这样是势是不合适的。E.Wimmer 等人[8]基于多极势的
概念,提出了一个新的方法,改进了对势的形状限制,对 LAPW 方法进行了
修正和发展,称之为全电势 LAPW 方法或 FLAPW 方法。FLAPW 方法的基本
思想是:一个局域电荷外面的势场仅通过多极矩而依赖于电荷,为了求得晶体
中球间区域的势场,只需要知道缓变的球间区内的电荷密度(r)和各个球内
22
第二章 研 究 方 法 与 程 序 介 绍
电荷的多极矩(快速收敛的)傅立叶表示。既然多极矩并不和唯一的电荷密度
对应,可以把球内的真实电荷用一个赝电荷来代替,使之和真实电荷分布有相
同的多极矩,但有快速收敛的傅立叶表示。这样将 Muffin-tin 球内(尤其是近
核处)快速变化的电荷分布用另一个平滑的赝电荷分布来代替,这由要求多极
矩相等来确定。然后求解泊松方程,得出球外空间及球面上的势场分布。最后,
用球内的真实电荷密度,用上面求得的球边界上的势场 作为边界条件,通过格
林函数方法可以求得球内势场的简谐展式。采用这一势场以后,在求能量矩阵
元时波函数中也考虑 l≠0 的非对称项。详细推导见文献[8~10]。在改进了势的
形状限制后,计算量会有所增加,但原则上说,可以用 FLAPW 方 法 处 理
各类体系,包括表面、界面、非密堆积的开结构晶体等。
总 之 , 在众 多 的 第 一 性 原 理 计 算 方 法 中 , 全 电 势 线 性 缀 加平 面 波 法
(FLAPW)是一种计算精确度高、性能可 靠的运算方法,它的可靠性已被理
论领域的科研工作者所认同。它特别适合晶体材料研究。铁电体是一类成 分和
结构都比较复杂的晶体材料,所以采 用这一方法比较合适。
§2.3 几种常用的程序
针对第一性原理计算编写的程序包很多。在铁电体及有关功能材料的研究
中,常用的主要有以下几种:
1. 针对 FLAPW 方法编写的各种计算程序中,由维也纳工业大学量子理论
计算研究小组开发的基于 LINUX 操作系统的 WIEN 系列程序最为著名。此程
序的主要优点是:(1)它是基于 LINUX 平台的软件,我们都知道 LINUX 的
高稳定性是众所周知的,在科研领域的许多软件都是基于 Unix 或 Linux 平台。
(2)它采用图形化友好界面,可以大大简化操作者的工作量,并降低了难度。
(3)它是采用 FORTRAN 语言编写的。FORTRAN 语言特别适合科研方面的
需求。我们在论文中采用的就是 WIEN97。
23