1.1 引言
在理论上最具有诱惑力，且将来最有可能开展真正意义上的材料设计的计算 就是求解体系的 Schrődinger 方程，也就是计算材料学中的第一原理计算。因为从 物理上讲第一原理计算有着完善的理论基础， 在求解体系 Schrődinger 方程的过程 中不涉及任何经验参数，需要输入的只是原子的核电荷数和一些模拟环境参量。 计算所求得的结果是体系 Schrődinger 方程的本征值和本征函数（波函数），有了 这两项结果从理论上讲就可以推导出体系的所有特性，如力学、电学、光学及磁 学等性质。随着计算机速度和精度的逐渐提高，再加上计算模型的不断改进，计 算模拟的结果就会和实验结果之间有着更好的可比性，从而增加计算模拟工作的 可靠性。1.2 多粒子体系的 Schrődinger 方程
[1] 多粒子体系的 Schrődinger 方程表达式为： 2 2 1 2 1 ( , , ) 2 N i N i i i U r r r t y y y m = ¶ = - Ñ + ¶å
h h L （1.1） 当体系的势场 U 与时间无关时，上面的 Schrődinger 方程的解可以用分离变量法进行简化。同时它的解也满足简化的 Schrődinger 方程，即定态 Schrődinger 方程： 2 2 1 2 1 ( , , ) 2 N i N i i i U r r r E Y y y m = é ù - Ñ + = ê ú ë
å
û h L （1.2） 对于多粒子体系，上述方程从数学上仍不能求解。为了求上述多粒子体系的 定态 Schrődinger 方程的解，必须在物理模型上作一系列的简化。分子轨道理论在 这方面作了三个近似处理 [2] ，即引入了非相对论近似、BornOppenheimer 近似和 轨道近似。 1.2.1 非相对论近似 电子在原子核附近运动但又不被原子核俘获，必须保持很高的运动速度。根 据相对论，此时电子的质量 μ 不是一个常数，而是由电子运动速度 ν，光速 c 和 电子静止质量 μ0 共同来决定： 0 2 1 v c m m = - （1.3） 用原子单位表示的多粒子体系定态 Schrődinger 方程为： 2 2 , 1 1 1 2 2 p q p p i i i i p pq ik pi p i p q i k p i Z Z Z E M < R < r r y y ì ü ï ï - Ñ - Ñ + + - = í ý ï ï þ îå
å
å
å
å
（1.4） 上述方程把电子的质量视为其静止质量 μ0，这仅在非相对论条件下才成立。 式（1.4）中，p 和 q 标记原子核，Rpq 为核 p 和核 q 间的距离，Zp 和 Zq 分别为核 p和核 q 所带的电荷，Mp 为核 p 的质量，用 i 和 k 标记电子，rik 为电子 i 和电子 k
间的距离，rip 为核 p 和电子 i 间的距离。
1.2.2 BornOppenheimer 近似
由于体系中原子核的质量比电子大 10 3 到 10 5 倍， 因而电子运动速度比原子核
起与变化后核力场相应的运动状态。这意味着，在任一确定的核的排布下，电子 都有相应的运动状态；同时，核间的相对运动可视为电子运动的平均作用结果。 据此， Born 和 Oppenheimer 处理了体系的定态 Schrődinger 方程，使核运动和电子 运动分离开来，这就是 BornOppenheimer 近似。用 ( , ) V r R 代表方程式（1.5）中 的势能项： , 1 ( , ) p q p pq ik pi p q i k p i Z Z Z V r R R r r < < =
å
+å
-å
（1.5） 分离变量后得到的电子运动方程为： 2 1 ( , ) ( , ) ( , ) ( ) ( , ) 2 i iy r R V r Ry r R E Ri y r R -å
Ñ + = （1.6） 原子核的运动方程为： 2 1 ( ) ( ) ( ) ( ) ( ) 2 p iF R E Ri F R E Ri F R -å
Ñ + = （1.7） 1.2.3 轨道近似 对于多电子体系， 上述简化后的定态 Schrődinger 方程仍然不可能严格求解， 原因是多电子势函数中包含了如 1 ik r - 形式的电子间排斥作用算符，不能分离变 量。因而近似求解多电子的 Schrődinger 方程还需要引入分子轨道法的第三个基 本近似——轨道近似，就是把 N 个电子体系的总波函数写成 N 个单电子波函数 的乘积： 1 2 1 2 2 ( ,x x , xN) ( )xi (x ) N(xN ) Y L = y y L y （1.8） 其中每一个单电子波函数 ψi(xi)只与一个电子的坐标 xi 有关。 这个近似隐含的 物理模型是一种“独立电子模型” ，有时又称为“单电子近似” 。用上式乘积波函 数描述多电子体系状态时，须使其反对称化，写成 Slater 行列式，以满足电子的 费米子特性，即：1 1 2 1 1 1 1 2 2 2 2 2 1 2 1 2 ( ) ( ) ( ) ( ) ( ) ( ) ( , , ) ( ) ( ) ( ) ( ) N N N N N N N x x x x x x x x x Nl x x x y y y y y y Y y y y = L L L M M M L （1.9） 根据数学完备集理论，体系状态波函数 ψ 应该是无限个 Slater 行列式波函数 的线性组合，即把式（1.9）中的单个行列式波函数记为 Dp，则有： p p p c D Y =
å
（1.10） 理论上，只要 Slater 行列式波函数个数取得足够多，通过变分处理就一定能 得到 BornOppenheimer 近似下的任意精确的能级和波函数。这个方法最大的优点 就是计算结果的精确性，它是严格意义上的从头算（abinitio）方法。但它也存在 着现在还难以克服的困难， 即此计算方法的计算量随着电子数的增多呈指数增加。 因此，这种计算对计算机的内存大小和 CPU 的运算速度有非常苛刻的要求，它使 得对较多电子数的计算变得不可能，如含有过渡元素或重金属元素体系的计算。 一般此方法多用于轻元素的计算，如 C、H、O、N 等，因而多用于量子化学计算。 这在很大程度上也是导致密度泛函理论产生的原因。1.3 密度泛函理论
[3-5] 对于一个实际体系来说， 一个多电子体系的定态 Schrődinger 方程如式 （1.1）。 在解这个多电子体系的 Schrődinger 方程的过程中遇到的最大困难就是求解电子 间重迭积分，即考察单个电子与其他电子间的相互作用，它的计算量令现在的计 算机难以承受。而密度泛函理论则视体系的多电子为电子气，在不同的区域具有 不同的密度，并可用一个公式来表述。进而上述多电子体系的 Schrődinger 方程被 简化为单电子的 Schrődinger 方程，从而使得求解变成可能。上述公式的引入也就 是第一原理计算的近似所在，因而形成相对意义上的从头算（abinitio）方法。它用密度泛函近似代替了先前的轨道近似，这种近似在计算中不能改善和消除，而 严格意义上的从头算方法中的轨道近似可以通过增加行列式的个数得到改善。在 HohenbergKohn 密度泛函理论（DFT） [6] 的基础上产生了著名的 KohnSham 方程 [7] （文献[8]，[9]对密度泛函理论作了开创性工作） ，即在密度泛函理论框架下的定 态 Schrődinger 方程： 2 2 , , , ( ) ( ) ( ) 2 m vext r vH r vxc r jks eksjk s æ Ñ ö - + + + = ç ÷ è ø h （1.11） 式中 vext 是原子核产生的外势，vH 是来源于电子密度的 Hartree 势， 2 , , ( ) ( ) occ k k r s r s r j =- ¯ =
å å
（1.12） 剩下依赖于密度的电子的多体效应可以通过 vxc 势或一函数的导数来表示， [ ] ( ) xc xc E v r d r dr = （1.13） 其中交换相关能 Exc 是电子密度 ρ(r)的泛函，但是其表达式不能被确定。Exc 的物 理意义在原理上是很简单的，即：当一个电子在一个多电子系统之间有交互关联 作用，如两个电子在同一时刻不能占据同一位置，从而产生交换相关能 Exc，最简 单的近似就是认为像在 ThomasFermi 理论中一样，Exc 仅与系统的体积有关，而 与电子密度的波动无关。事实证明它是一个比较粗糙的模型，它仅对电子密度波 动不大的系统适用。一个较好的近似就是构造一个函数，使得 Exc 不仅依赖于电 子密度 ρ(r)， 而且依赖于它的导数。 于是就出现了下面几种最为常见的泛函表达式。 1.3.1 局域密度近似（LDA） [1015] 在局域密度近似下的交换关联能的表达式为： [ ] d ( )[ ( ( )) ( ( ))] d ( ) [ ( )] LDA xc x c xc E r =ò
rr r e r r +e r r =ò
rr r e r r （1.14） 此式假定体系是无极化的电子气系统，式中 ex 和 ec 代表每个电子的交换能和 相关能。它们有明确的表达式。1.3.2 自旋极化局域密度近似（LSD） [15] 在自旋极化局域密度近似下的交换关联能的表达式为： 3 [ , ] d ( )[ ( ( )) ( ( )) ( ( ), ( ))] LSD xc x c s E r r- ¯ =
ò
rr r e r r f V r + e r r V r （1.15） 式（1.15）中，ex 和 ec 与 LDA 中的意义相同。其中 2 , ( ) occ k k r r- =å
j - （1.16） 2 , ( ) occ k k r r¯ =å
j ¯ （1.17） ( )r ( )r ( ) r r =r- + r¯ （1.18） ( ) ( ) ( ) ( ) r r r r r r V r r - ¯ - ¯ - = + （1.19） 4 / 3 4 / 3 1 ( ) [(1 ) (1 ) ] 2 f V = +V + - V （1.20） 1.3.3 广义梯度近似（GGA） [1627] 在 GGA 近似中，相关积分函数采用的表达式为： 3 [ , ] d [ ( , ) ( , , )] GGA c c s s E r r- ¯ =ò
rre r z + H r z t （1.21） 式中： 1/ 3 2 4 / 3 3 4 2 9π 2(3π ) B s s= Ñ r = æç ö ÷ a Ñ r è ø （1.22） 1/ 2 1/ 6 1/ 2 π 9π 4 4 s s t r f æ ö æ ö = ç ÷ ç ÷ è ø è ø （1.23） 2 / 3 2 / 3 1 ( ) (1 ) (1 ) 2 f z = é +z + - z ù ë û （1.24） 2 3 2 2 2 4 1 ln 1 1 At H r t At A t b f g ì é + ù ü ï ï = í + ê ú ý + + ï ë û ï î þ （1.25） 3 1 exp[ c( , ) /s ] 1 A e r r b g z f æ ö = ç ÷ - - è ø （1.26）式中的常数 b = 0.066725， g = 0.03191。 交换函数关联能采用的表达式为： 3 [ ] d ( ) (s) x x x E r =
ò
r re r F （1.27） 式中： 2 1 1 / x F s k k m k = + - + （1.28） 其中，常数 k = 0.804和 m = 0.21951。1.4 CASTEP 理论背景及基本原理方法
CASTEP（Cambridge Sequential Total Energy Package）为一商业软件包，始 于剑桥大学凝聚态理论研究组开发的一系列程序。它是一个基于密度泛函方法的 从头算量子力学程序，计算时仅需输入的是所研究体系的初始几何构型和组成原 子的数目。在这种方法中，离子势被赝势（即只作用于系统价电子的有效势）所 替代。电子波函数通过一平面波基组扩展，电子－电子相互作用通过密度泛函理 论得以充分的考虑。结合赝势和平面波基组的应用，体系中所有原子上的作用力 的计算变得极为容易。这些计算给出系统的基态能量和电子密度，允许计算与总 能量相关的任何物理量，如：晶格常数、弹性常数、几何结构、能带结构、合金 形成热、结合能等。它可以用于模拟固体界面和表面的性质，适用于多种材料体 系，包括陶瓷、半导体、金属及其合金。 1.4.1 总能量的计算 1.4.1.1 总能表达式 在 CASTEP 对各量的计算中，其中总能的计算是最为关键的量之一，在密度
泛函理论框架下，CASTEP 程序计算的总能分为三个部分 [5,6] ，由下式表示
[ ]
[ ]
[ ]
[ ]
t xc E r =T r +U r + E r （1.29） 式中，r 是体系中的电荷密度， [ ] T r 和 [ ] U r 是分别用电荷密度表征的体系的动 能和静电能， E r xc [ ] 是包含体系中所有多体效应的交换关联能。因为电荷密度是 一个可测量，它的值可以通过构造体系波函数的方法来获得，采用轨道近似（详 见 1.2.3 节） 把体系中 N 个电子体系的总波函数写成 N 个单电子波函数的乘积 （见 式（1.8） ），且这些单电子波函数满足正交关系： i j ij f f = d （1.30） 于是利用式（1.30）可求得电荷密度的值： 2 ( ) i ( ) i r r r =å
f （1.31） 由式（1.29）和式（1.31）可得到式（1.29）中各项的表达式： 2 2 n i i i T =å
f -Ñ f （1.32） 1 2 1 2 1 2 1 1 2 1 1 2 1 1 1 1 1 ( ) ( ) ( ) ( ) ( ) ( ) 2 1 1 ( ) ( ) ( ) 2 ( ) ( ) ( ) 2 n N i i i j i j i i j N N N e N NN Z U r r r r r r R r r r Z Z R R Z Z Z r r r R r r r R R V r r V r V a a a b a b a a b a b a a a a b a a b f f f f f f r r r r r < < - = + - - + - = - + + - - - = - + +åå
åå
å å
å
å å
（1.33） 式（1.33）中，Zα 表示的是 N 个原子的体系中核 α 的电荷。第一项，ρVN 表示电 子－核之间的吸引作用；第二项，ρVe/2 表示的是电子－电子之间的排斥作用；最 后一项， V NN 表示的是原子核－原子核之间的排斥作用。对于式（1.29）中最后 一项交换关联部分 Exc r [ ] ，需要作一些近似，通常所用到的是局域密度近似 LDA （详见 1.3.1 节）和广义梯度近似 GGA（详见 1.3.3 节）。在 CASTEP 软件中，对 局域密度近似（LDA）只提供了一种函数，即：CAPZ [10,12,26] 。而在广义梯度近似（GGA）中提供了三种函数：PW91 [22] （PerdewWang GGA 近似） 、PBE [23]
（PerdewBurkeErnzerhof 函数）和 RPBE [27] （修正的 PBE 函数） 。
综合式（1.14）、式（1.32）和式（1.33）可得体系中总能量的表达式为：
[ ]
[
]
2 1 1 1 ( ) ( ) ( ) 2 2 e t i i xc N NN V r E r = f + -Ñ f + r r éêe r r + -V ù ú + V ë ûå
（1.34） 要得到准确的体系总能量，Et 中各变量应满足正交归一化条件： 0 t ij i j i j E d e f f dr -åå
= （1.35） 其中 e ij是 Lagrangian 乘子。在求解式（1.34）的过程中，出现了一系列的方程式， 如 Hartree 方程、HartreeFock（HF）方程、KohnSham（KS）方程等。其中应 首推在 HobenbergKohn（HK）定理下建立起来的 KS 方程。 1.4.1.2 KS 方程 Hobenberg 和 Kohn [6] 在 1964 年指出： 对于一个处于外势场 V(r)中的 N 电子体 系，外势场 V(r)和体系基态电荷密度 r (r ) 之间是一一对应的（可相差常数）。由于 体系电子数 N 与电荷密度 r (r ) 直接联系，这样就决定了多电子 Schrődinger 方程 解的 N 和外势场 V(r)都由电荷密度 r (r ) 唯一确定。这样体系的基态总能量也是密 度 r (r ) 的唯一确定的泛函。随后他们又进一步证明了：在粒子总数 N 保持不变的 约束条件下，所得到的基态密度就对应于系统的总能量泛函取最小值。这样，如 果知道了体系能量的密度泛函表示，就可以从它的变分条件求出基态密度，进而 得到体系基态的所有物理性质。 随后，Kohn 与 Sham 又发展了此理论，于 1965 年 [6] 推出了著名的 KS 方程：[ ]
2 2 Vn Ve mxc r fi e fj i ì-Ñ ü ï ï - + + = í ý ï ï î þ （1.36） 为了区分输出结果 Exc，在式（1.36）中用 μxc 来表示交换关联势。在局域自 旋密度近似下，交换关联势 μxc 可表示为：( ) xc xc m re r ¶ = ¶ （1.37） 求解式（1.36）中的能量的本征值，用它来表示体系的总能为：
[ ]
[ ]
1 1 ( ) ( ) 2 e t i xc xc NN i V r E = e + r r êée r -m r - ù ú + V ë ûå
（1.38） 1.4.1.3 赝势 至于如何处理式（1.38）中的原子核－电子的相互作用势（即： 2 e V r ）这一 项，CASTEP 引入了赝势理论。赝势是一种有效的虚拟的原子核－电子的相互作 用势 [1] 。它最初是从能带计算中的正交化平面波方法所引出的一个概念，其普遍 的定义如下：在原子球内部用一个假设的势代替真正的原子势，求解原子间空间 的薛定谔方程式，若所得能量本征值不变（即与真实原子势的结果一样） ，则此虚 拟的势便是赝势。由于在正交化平面波方法中出现的波函数正交项趋向于抵消真 实势能项，从而可以给出一个有效势，此有效势比真实的晶体势弱很多，当时就 将它称为赝势。用这个较弱的势代替薛定谔方程中的真实的势，得到一个有效波 函数，此波函数称为赝波函数。由于赝势比真实势弱，所以赝波函数比真实波函 数平滑，因而展开式收敛较快，但本征值不变，因此可得到正确的能谱。虽然正 交平面波方法引出的赝势仅是赝势的一种具体形式，但通过理论证明：任何一个 假想（或虚拟）的势场，只要它能给出原子半径 r 处径向波函数正确的对数微商 0 Ll(E)，它就是一个赝势。由于波函数总是对应于某一确定的能量 E，所以赝势也 是对应于某一能量 E。此外，还存在赝势的另一种表述法：由于一个势的散射特 性完全决定于 r 处的对数微商 L0 l(E)，而并不依赖于 r< 区域中势或波函数的具 r 0 体细节。因此可将赝势选择为方便的形式，含有可调节参数，调节这些参数已给 出正确的散射。具体确定赝势的方法包括：①模型势法；②经验势法；③自洽势 法。常见的模型势有空模型、AbarankovHeine 模型。研究发现，根据一组试验数据所确定的半经验赝势用于预测其他不同试验，时常可以得到很相符的结果。赝 势理论已成功地用于除固体能带计算外的许多物质结构研究，甚至包括某些生物 物质。在 CASTEP 软件中，提供了两种赝势，即：模守恒（Normconserving）和 超软（Ultrasoft）赝势。 1.4.1.4 自洽计算（SCF） 在实际计算中，对于式（1.36）中的分子轨道（MOs）波函数 f i通常可以用 扩展原子轨道（AOs）波函数的方法来获得： i Ci m m m f =
å
c （1.39） 式（1.39）中原子轨道波函数c m 通常被称为原子基函数，C i m 是 MO 扩展系数。 目前有很多种基函数可以选择，如高斯（Gaussian）函数 [28] 、Slater 函数 [29] 、平面 波基函数 [30,31] 、数字轨道（Numerical Orbitals）函数等。在 CASTEP 程序包中采 用的是平面波基函数。与 MOs 不同的是，AOs 波函数不是正交的，于是 KS 方 程（式（1.36） ）变为： HC= eSC （1.40） 其中，哈密顿矩阵算符 H 和 S 算子由式（1.41）和式（1.42）确定 2 1 1 ( ) ( ) 2 v N e xc v Hm = xm r -Ñ -V +V + m rr x r （1.41） 1 1 ( ) | ( ) v v Sm = cm r c r （1.42） 因为哈密顿矩阵算符 H 依赖波函数 C，式（1.40）需要通过下面的自洽计算 （SCF）程序来求解： （1）设定一组初始 Ci m 。 （2）构造一组初始分子轨道（MOs）波函数 f i。 （3）由式（1.30）构造体系的电荷密度r 。 （4）用r 构造 V 和 e m xc。（5）构造 Hm v 。
（6）解方程（1.40） ，得到一组新的 Ci m 。
（7）由此可以求得新的 f i和r 。
（8）判断：如果 rnew = rold，由式（1.38）计算总能 E ，自洽计算程序终止； t
如果 rnew ¹ rold，返回（4） 。 一般来说，对于有机分子体系，大约只需要十多次自洽积分计算便能获得 6 10 new old r -r < - 的收敛值；对于金属体系，则需要更多的自洽积分计算。 1.4.1.5 超晶胞方法及周期边界条件 在 CASTEP 中运用的基本方法是在周期性重复的单位晶胞上执行的， 这种 单位晶胞通常被称作“超晶胞” 。假定晶胞中的晶格矢用 ai( i = 1, 2, 3 )来表示， 例如在立方晶胞中， a = 1 (0,1, 0) ， a = 2 (0,1, 0) ， a = 3 (0, 0,1) 。则基函数必须满足 晶 体 平 移对称性，因此有 ( )xm r =xm(r+ R) ，式中 r 表示晶体中任意一点， 1 1 2 2 3 3 R=n a +n a + n a ，n 为整数。为了满足这个条件，CASTEP 采用了具有晶格 周期性的一组平面波作为周期性基组： , ( , )k R eik R x r( ) Y = （1.43） 原则上讲，要准确地描述晶胞空间，需要无穷多的 k 矢和 R 矢。而实际上， 平面波的数目由截断动能来决定，一般只需要少量的 k 空间中的点（大约 1～10 2 个）便能准确地描述体系的能带结构。 1.4.2 CASTEP 计算输出结果 程序计算输出的结果主要包括以下几项： 1.4.2.1 电子密度数据 利用输出的电子分布密度，通过 CASTEP 程序自带的相应数据处理和绘图软 件， 可以绘制所计算的超晶胞内任意面的电子分布密度的等高线图或三维立体图，
从而可以了解各原子间的成键、电荷转移等情况。 1.4.2.2 超晶胞总能 通过此项数据和各孤立原子的总能即可获得晶体的结合能，用此项数据和形 成晶体各单位原素的固有能量可获得合金的形成热，分析合金结构的稳定性。通 过不同晶格常数下总能曲线 可以得到晶胞的平衡 晶格常数。也可直接通过 CASTEP 几何优化计算来得到体系的平衡晶格常数。 1.4.2.3 晶体的能带结构和态密度 通过计算晶体超晶胞的总能，可以进一步计算晶体的能带结构、晶体总态密 度，或各组成原子的态密度和分波态密度，从而可以了解晶体的电学性质。 1.4.2.4 弹性常数 CASTEP 可以计算任意晶体结构的弹性常数 Cij、体弹性模量、剪切模量、泊 松比、压缩比等，可以进一步用这些数据来分析晶体结构的稳定性和晶体的力学 性质等。 【参考文献】 [1] 周世勋．量子力学教程[M]．北京：高等教育出版社，1992，2132． [2] 廖沐真，吴国是，刘洪霖．量子化学从头计算方法[M]．北京：清华大学出 版社，1984，110． [3] Kohn W. Density Functional Theory[M]. New York, 1995, 120.
[4] Paul Z, Stefan K, John P. Perew. Density functional from LDA to GGA[J]. Comput Mater Sci, 1998, 11(1): 122127.
[5] 黄美纯．密度泛函理论的若干进展[J]．物理学进展，2000，20(3):199219.
[6] Hohenberg P, Kohn W. Inhomogeneous electron gas[J]. Phys. Rev. B, 1964, 136:
[7] Kohn W, Sham L J. Selfconsistent equations including exchange and correlation effects[J]. Phys. Rev. A, 1965, 140(4A): 11331138.
[8] Perdew J P. Accurate density functional for the energy: realspace cutoff of the gradient expansion for the exchange hole[J]. Phys. Rev. Lett., 1985, 55(16): 16651668.
[9] Perdew J P, Wang Y. Accurate and simple density fuctional for the electronic exchange energy: generalized gradient approximation[J]. Phys. Rev. B, 1986, 33(12): 88008802.
[10] Hedin L, Lundqvist B I. Explicit local exchangecorrelation potentials[J]. J. Phys. C: Solid State Phys, 1971, 4(14): 20642083.
[11] Barth UV, Hendin L. A local exchangecorrelation potential for the spin polarizied[J]. J. Phys. C: Solid State Phys, 1972, 5(13): 16291642.
[12] Ceperley D M, Alder B J. Ground state of the electron gas by a stochastic method[J]. Phys. Rev. Lett., 1980, 45(7): 566569.
[13] Nusair M, Wilk L, Vosko S H. A comparison of spindensity functional calculations for the Knight shift in Mg[J]. J. Phys. F: Met. Phys, 1981, 11(8): 16831690.
[14] Janak J F. Itinerant ferromagnetism in fcc cobalt[J]. Solid State Commun.,1978, 25(2): 5355.
[15] Perdew J P, Wang Y. Accurate and simple analytic representation of the electrongas correlation energy[J]. Phys. Rev. B, 1992, 45(23): 1324413249. [16] Langreth D C, Mehl M J. Beyond the localdensity approximation in calculations
of groundstate electronic properties[J]. Phys. Rev. B, 1983, 28(4): 18091834. [17] Perdew J P. Density functional approximation for correlation energy of the
[18] Becke A D. Density functional calculations of molecular bond energies[J]. J. Chem. Phys., 1986, 84(8): 45244529.
[19] Becke A D. Density functional exchange energy approximation with correct asymptotic behavior[J]. Phys. Rev. A, 1988, 38(6): 30983100.
[20] Lee C, Yang W, Parr R G. Development of the ColleSalvetti correlation energy and correlation[J]. Phys. Rev. B, 1988, 37(2): 785789.
[21] Perdew J P, Ziesche I P, Eschrig H. Electronic structure of solids[M]. Berlin: Akademie Verlag, 1991, 11.
[22] Perdew J P, Chevary J A, Vosko S H, et al. Fiolhais C, Atoms, molecules, solids and surfaces: applications of the GGA for exchange and correlation[J]. Phys. Rev. B, 1992, 46(11): 66716687.
[23] Perdew J P, Burke K, Ernzerhof M. Generalized gradient approximation made simple[J]. Phys. Rev. Lett., 1996, 77(18): 38653868.
[24] 程大勇. 金属材料弹性性质与电子结构的第一原理计算[D]，博士论文， 2001，中国科学院金属研究所.
[25] Hedin L, Lundqvist B I. Explicit local exchange correlation potentials[J]. J. Phys. C, 1971, 4: 20642083.
[26] Lundqvist S, March N. Eds. Theory of the Inhomogeneous Electron Gas[M]. New York: Plenum, 1983.
[27] Hammer B, Hansen L B, Norskov J K. Improved adsorption energetics within densityfunctional theory using revised PerdewBurkeErnzerhof functionals[J]. Phys. Rev. B, 1999, 59: 74137421.
[28] Andzelm J, Wimmer E, Salahub D R. Spin density functional approach to the chemistry of transition metal clusters: Gaussiantype orbital implementation, in The Challenge of d and fElectrons: Theory and Computation[M]. ACS Symp. ser., 1989, 394.
[29] Versluis L, Ziegler T. The determination of molecular structures by density functional theory: The evaluation of analytical energy gradients by numerical integration[J]. J. Chem. Phys., 1988, 88: 322328.
[30] Ashcroft N W, Mermin N D. Solid State Physics[M]. Holt Saunders, Philadelphia, 1976, 113.
[31] Payne M C, Teter M P, Allan D C, et al. Iterative Minimization Techniques for Ab Initio Total Energy Calculations: Molecular Dynamics and Conjugate Gradients[J]. Rev. Mod. Phys., 1992, 64:10451097.
