锗的热力学性能第一性原理研究--计算材料学课程设计.doc

湖南工业大学

课 程 设 计

资 料 袋
理 学院（系、部） 2011 ~ 2012 学年第 一 学期
课程名称 计算材料学 指导教师 雷军辉 职称 讲师
学生姓名 刘进萍 专业班级 应用物理081班 学号 08411200201
题 目 锗的热力学性能第一性原理研究
成 绩 起止日期 2011年 12月 4日 ～ 2011年 12 月 12 日

目 录 清 单
序号
材 料 名 称
资料数量
备 注
1
课程设计任务书

2
课程设计说明书

3
课程设计图

张
4

5

6

湖南工业大学
课程设计任务书
2011—2012 学年第 1 学期

理学院 学院（系、部） 应用物理学 专业 081 班级
课程 名称： 计算材料学
一、 设计题目： 锗的热力学性能第一性原理研究
二、 完成期限：自 2011 年 12 月 4 日至 2011 年 12 月 12 日共 2 周

内
容
及
任
务

1. DFT基本理论，CASTEP使用方法
2. 晶体模型的建立与几何优化，相关性质的计算。
3. 计算热学性质
4. 结果分析
5. 报告写作与修改
进
度
安
排
起止日期
工作内容
11-12-4-6
熟悉DTT理论，软件安装，认识界面，熟悉基本操作
11-12-7
建立NaCl模型，进行结构优化，计算物理性质
11-12-8
计算Si的物理性质，力学性质
11-12-9
计算Ge的热学性质
11-12-10-12
写出课程设计的总结实验报告.，修改成文
主
要
参
考
资
料
[1] Kohn W, Sham L J, Self-consistent equations including exchange and correlation effects [J]. Physical review, 1965, 140(4):A1133-A1338.
[2] Hohenberg P, Kohn W. Inhomogeneous electron gas [J]. Physical review, 1964,136(3):B864- B871.
[3] 谢希德, 陆栋.固体能带理论[M].上海：复旦大学出版，1998.
[4] Perdew J P, Chevary J A, Vosko S H. Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation[J]. Physical review B, 1992, 46(11): 6671-6687.

指导教师（签字）： 年 月 日
系（教研室）主任（签字）： 年 月 日

23

（计算材料）
设计说明书

锗的热学性能第一性原理研究

起止日期： 2011 年 12月 4日 至 2011 年 12月 12日

学生姓名
刘进萍
班级
081
学号
08411200201
成绩

指导教师(签字)

理学院（部）
2011年 月 日

锗的热学性能第一性原理研究
背景
近年来，随着材料、物理、计算机和数学等学科的发展，应用计算的方法研究材料的结构、能量和性能已成为一门迅速发展的新兴学科-计算材料学。这种方法不仅能进行材料的计算模拟，而且能进行材料的计算机设计和相关性能的预测。随着计算机技术的飞速发展，第一性原理计算的方法在材料的结构和性能等方面的研究已取得了巨大的成功，第一性原理的方法是基于量子力学理论，从电子运动的层次研究材料的结构和相关性能。目前，CASTEP软件的主要功能是对半导体、非线性光学材料、金属氧化物、玻璃、陶瓷等固体材料，对电子工业、航空航天以及石化、化工等工业领域有着非常重要的战略意义。对这些材料而言，其电子的结构与性质，以及表面和界面的性质与行为都非常重要。CASTEP的量子力学方法，为深入了解固体材料的这些性质并进而设计新的材料，提供了强有力的工具。
基于密度泛函平面波赝势方法的CASTEP软件可以对许多体系包括像半导体、陶瓷、金属、矿石、沸石等进行第一性原理量子力学计算。典型的功能包括研究表面化学、能带结构、态密度、热学性质和光学性质。它也能够研究体系电荷密度的空间分布和体系波函数。CASTEP还可以用来计算晶体的弹性模量和相关的机械性能，如泊松系数等。半导体和其他固体材料的许多性能由电子性质决定，而电子性质又由原子结构决定，特别是缺陷在改变电子结构上的作用对半导体性质尤为重要。分子模拟，特别是量子物理技术，可用来预测原子和电子结构及分析缺陷对材料性能的影响。CASTEP能有效的研究存在点缺陷、空位、替代杂质、位错等的半导体和其它材料中的的性能。除此以外，它还可以被用来计算固体的振动性质，如声子色散关系、声子态密度等。这些计算结果可以用来分析表面吸附的振动性质，可以解释实验中的振动谱，可以研究在高温高压下的相稳定性等等。总的来说，它可以实现如下的功能：
1.计算体系的总能。
2.进行结构优化。
3.执行动力学任务：在设置的温度和关联参数下，研究体系中原子的运动行为。
4.计算周期体系的弹性常数。
5.化学反应的过度态搜索。
除此之外，计算一些晶体的性质，如能带结构、态密度、声子色散关系、声子态密度、光学性质、应力等。
下面介绍一下密度泛函理论、交换关联泛函近似、赝势方法和K-S方程迭代解法。

一、 基础理论
1. Hohenberg-Kohn 定理和密度泛函理论
密度泛函理论(DFT)是用量子力学的理论求解多电子体系基态能量方法，其核心是用电子密度函数取代波函数作为研究的基本量，由Hohenberg 和Kohn 在1964 年创建[1,2]。根据量子力量的相关知识，大量电子和原子核相互作用的多粒子体系，在非相对论前提下，系统粒子运动的波函数可以由以下定态薛定谔方程来描述：
(1-1)
哈密顿量 仅考虑电子-电子作用、电子-原子核作用、原子核-原子核作用以及各个粒子的动能，对其它外场的情况可忽略。因此其哈密顿量可以写成如下形式：
(1-2)
其中，
(1-3)

(1-4)

(1-5)

对于上述方程，是无法直接求解的，必须对多粒子系统的电子能级计算采用一些简化和近似。在实际的多粒子体系中，原子核的质量远远大约电子， 但是运动速度比电子小的多。因此考虑粒子运动时，将原子核的运动和电子的运动分开，考虑核的运动时忽略其电子分布，考虑电子运动时假定原子核处于相对静止的状态，这就是绝热近似[3]。通过近似，可以独立的处理原子核运动和电子的运动，因此可以将薛定谔方程写成电子运动方程和原子核运动方程。其电子运动方程是：
(1-6)

原子核的运动方程：
(1-7)
通过绝热近似，得到了多电子的薛定谔方程，但不能实际求解，要求解上述方程， 必须将多电子问题简化为单电子问题。单电子近似理论的源于H.Thomas 和E.Fermi 在1927 年的工作，就是用粒子数密度表示多粒子的基态系统的能量。P.Hohenberg 和W.Kohn 根据的均匀电子气的理论提出著名的Hohenberg-Kohn 定理[1]，这个定理包含如下内容：不计自旋的情况下，将粒子数密度函数表示成全同费米子系统的基态能量的唯一泛函；在粒子数不变的情况下，能量泛函对正确的粒子数密度取等于基态能量的极小值。因此，对于基态非间并多粒子系统，不考虑自旋的条件下，其哈密顿算符为
(1-8)
式(1-8)中，外场作用看成原子核-电子作用，相同的局域势对外场的作用用表示。对于给定的外场，多电子系统的能量表示成电子数密度的泛函为：
(1-9)
(1-10)
(1-11)
(1-12)
式中，包括体系中电子之间的相互作用能和电子的动能， 是外场对电子的作用能， 是系统中原子核间的排斥能。在式(1-10)中，前两项表示无相互作用粒子模型的动能和库仑排斥能，复杂的电子相互作用用交换关联能表示。根据Hohenberg-Kohn 定理，假设能得到能量泛函E(r )，然后就能将电子数密度r 变分，就能确定系统的基态和基态所有的性质，因此确定E(r )成为问题的关键所在，而要确定能量泛函E(r )，必须要确定动能泛函T[r ]、电子数密度 以及交换关联泛函。为了解决上述问题，W.Kohn 和L.J.Sham 提出了如下假设：假定已知无相互作用的电子系统和未知的有相互作用的电子系统密度函数相同，未知的相互作用电子系统的动能泛函T[r ]可用已知的无相互作用电子系统的动能泛函 来代替；假定密度函数 用N 个单电子波函数构成，于是有：
(1-13)

则
(1-14)
对能量泛函进行变分得到
(1-15)

(1-16)

式(1-13)、(1-15)和(1-16)就是Kohn-Sham 方程。这个方程的核心就是有相互作用动能泛函能否用未知的无相互作用的动能泛函来代替。而将所有复杂问题都归入中，所以求解Kohn-Sham 方程的关键是找到准确的，这样密度泛函理论精确求解量子多体问题的中心是构造交换关联泛函。
2．交换关联能近似
根据密度泛函理论，能将多电子的基态特性问题转化成等效的单电子问题， 而其它所有复杂问题都归结到交换关联能泛函，但是交换关联泛函是未知的。因此得到可靠并准确的交换关联，成为求解Kohn-Sham 方程的关键。W.kohn 和L.J.Sham 提出了交换关联泛函局域密度近似(LDA，Local Density Approximation)，其基本思想是：在局域密度近似中，利用均匀电子气密度函数来获得非均匀电子气密度泛函。对变化平缓的密度函数，非均匀交换关联能密度用均匀电子气代替，则可表示为：
(1-17)
相应的局域交换关联势可以表示为：
(1-18)

局域密度近似虽然在大多数的材料计算中显示出巨大的成功，但是由于点r处的交换关联作用仅依赖于点r处的近邻和近邻的电荷密度，因此，对于与均匀电子气或空间变化缓慢的电子气相差太远的体统，LDA 不适用。因此，人们对局域密度近似应用多种方法进行修正，应用较广的是广义梯度近似(GGA)，其泛函与局域密度和密度梯度都有关[4]，因此能更好的描述真实体系的电子密度的不匀性，其交换关联能密度泛函 可表示为
(1-19)
目前，在交换关联泛函GGA 的构建上有两个方向，一个是Becke 为首的， 这类泛函包含若干个实验参数，这些参数通过计算和实验数据来获得，这种形式的好坏由实践的工作所决定的。另一个是Perdew-Wang 91 的，这类泛函以物理规律为基础，不包含实验参数。随着研究的不断深入，不仅出现了非局域的相互作用交换关联泛函，还有密度高阶梯度的近似交换关联泛函，如Vaner Wals 和Meta-GGA 等。GGA 和LDA 相比在能量精确度和开放体系方面更有优势。
3. 赝势法
在晶体的近自由电子能带计算中，计算量大而且收敛速度慢。对于固体而言，价电子的化学性质活泼，对于结构和性质的影响较大，而内层电子的能带较窄，较稳定，而且相邻原子的作用对内层电子的状态影响较小。因此，人们关注的是价电子，将原子核和内层电子近似看出粒子实。对于固体中的价电子波函数而言，在离子实的内部区域，变化剧烈，存在若干个节点；而在离子实之间的区域，变
化平缓。离子实内部的这一特点要求价电子波函数与内层电子波函数正交，而价
电子与内层电子波函数正交起了一种排斥势的作用，在很大程度抵消了离子实内
部V(r )的吸引作用。据此，离子实内部的势函数用假想势代替，在离子实之间
的区域波函数和电子的能量本征值保持不变的条件下求解固体单电子波函数方
程，假想的势叫赝势，用赝势求出的波函数叫赝势波函数。对于多原子固体而言，
根据波函数的不同特征坐标空间被分成c r 以内的原子核区域(芯区)和以外的其它区域两部分(假定存在某个截断距离c r )。芯区(r< c r )其波函数与紧邻原子波函数相互作用很小，赝势和赝势波函数变化缓慢，比较平坦；芯区外(r> c r )价电子波函数相互交叠作用，和真实的势和波函数相比，其形状和幅度都一样。目前，除了经验赝势、半经验的模型赝势外，还有没有附加经验参数的第一性原理从头算原子赝势，包括模守恒赝势(norm-conserving)、超软赝势方法(ultra-soft pseudo-potentials，USPP)[5]
4. 分子轨道的自洽求解
(1)分子轨道的自洽场方程
密度泛函理论是基于Hohenberg--Kohn定理，该定理表明体系基态的性质由电荷密度决定，体系的总能量是电荷密度ρ的函数。总能Et可以表达为：
(1-20)
T[ρ]是密度为ρ的电子的动能，U[ρ]是经典的库仑相互作用静电能，Exc[ρ]包括了多体相互作用对总能的贡献，其中交换-关联能是主要的部分。我们从波函数Ψ来构造电荷密度。对于波函数Ψ可以写成具有反对称性单?粒子波函数（分子轨道）的Slater行列式：
(1-21)
当分子轨道是正交时，即
(1-22)

电荷密度可表示为：
(1-23)
由总能的表达式和电荷密度的表达式，动能项（原子单位）可表示为：
(1-24)
库仑相互作用项为：
(1-25)
方程中表示原子核的带电量，表示电子与核的吸引作用, 表示电子与电子的排斥作用，表示核与核的排斥作用。总能表达式的最后一项交换-关联能需要作一些近似，比如局域自旋密度近似（LSDA），广义梯度近似（GGA）等。这样，总能的表达式可写成：
(1-26)
利用分子轨道的正交归一性，基态的能量由上式对密度的变分得
到： (1-27)
化简上式，得到Kohn-Sham方程：
(1-28)
式子中是与交换-关联能对应的交换-关联势。事实上，分子轨道可以通过原子轨道来展开，也就是说分子轨道是原子轨道的线性组合，可以表为：
(1-29)
在这里原子轨道称为原子轨道基函数，为扩展系数。也可以使用其它的基函数，而在CASTEP中分子轨道用平面波基来展开。不象分子轨道，原子轨道是非正交的，在使用原子轨道基函数展开时，Kohn-Sham化为下列形式：
(1-30)

其中
(1-31)
(1-32)
它是分子轨道的自洽场方程，是非线性方程，只能用迭代方法求解。
二． CASTEP软件的使用方法
1、模型的建立方法

点击file，选择new，则出现下图：

图2.1
其中有多个选项，可以选择3D Atomistic，点击确定，打开一个工作窗口。

图2.2

确定空间点群，选择lattice Parameters，就可以建立晶包结构。
最后选择加入到原子：选择built---Add Atoms

图2.3
在原子相应位置上添加原子，就可以建立计算模型。
以NaCl为例，重复上述过程就可以建立如图所示模型：

图2.4

2.计算任务的设置
在materials studio软件中行任务设置，主要是通过CASTEP应用窗口中的工具条之一“Calculation”来进行。我们可以更改工具框中的相应选项，来配置诸如：“电子选项”、“结构优化选项”、和“电子和结构性质选项”等。这几个选项是我们在运用CASTEP进行性质计算研究中，非常重要的几个技术参数。其中，“电子选项”是很多其它计算任务也要涉及的。在CASTEP中还有如动力学、结构优化、弹性常数、过渡态等计算的设置。在程序运行之前，从研究的问题出发，要将软件中关键的一些任务参数设置成符合计算需要的值，我们才能得到所期望的运算结果。
图2-5

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

表2.1
如表2.1所示，主要分为差、中等、好、超好四个等级。在涉及SCF收敛精度、K点取样精度、截断能等的设置时都要进行适当的选取。比如，在我们做掺杂纳米碳管结构优化时，选取的SCF收敛精度为“好”K点间隔为0.02/A,截断能为350eV。

b、交换-关联函数的设置CASTEP提供了两种交换-关联函数的设置，一种是局域密度近似(LDA),它使用的是CA-PZ形式的赝势，是引用Ceperley,D.M.;Alder,B
[10]和Perdew,J.P.;Zunger[11]采用的赝势形式。另外一种是广义梯度近似(GGA)，它有三种形式可供选择，分别是PBE[9]、RPBE[8]、PW91[7],对应的是三种不同的广义梯度近似形式。以上两种交换-关联函数LDA和GGA及其对应的可选形式是通过”CASTEP Calculation”中的“Setup”选项来进行配置，如图2-6和2-7所示,改变其中”Functional”的类型和每一类型对应的函数产生形式，即可得到相应的设置。在我们做纳米碳管结构优化任务时,选择的是广义梯度近似(GGA)下的PBE形式的关联函数，如图2-7所示。
图2-6 LDA设置 图2-7 GGA设置

c、赝势的设置
在这里我们介绍在CASTEP中设置Ultrasoft赝势和
Norm-conserving赝势这两种赝势，它们有各自的优势，Norm-conserving赝势一般适用于金属体系，对于我们研究的纳米碳管，更适合使用的是Ultrasoft赝势。设置的方法是在”CASTEP Calculation”中的“Electronic”选项中，在“Pseudopotential”的下拉框来进行选择。如图2-8所示。

图2-8赝势的设置
d、截断能的设置
CASTEP中分子轨道是通过平面波基来扩展。平面波基的数目是通过截断能的高低来控制，选择的截断能过低会影响计算的结果的正确性，而选择的过高会影响计算量。因此，在计算中要选取合适的截断能。一般情况下，在计算前可以选取几个截断能来试验，哪个更合适，在保证计算精度的前提下选择尽可能低的截断能。因此，截断能是CASTEP计算中最重要的参数之一。
设置截断能的最简单方式是在”CASTEP Calculation”中的

图2-9 截断能的设置
“Electronic”选项中选择”more”,然后在对话框中选择”Basic”,在“Use custom energy cut-off”中填入欲设置的值（图2-9）。在计算超细纳米碳管的计算中我们首先尝试了350eV和400eV的截断能，结果显示，两者计算的总能差别不到0.1eV,差别非常小，于是我们采用350eV进行计算。

e、 K点设置
（1）布里渊区的设置
布里渊区的设置是通过K点的设置来反映的。在K点的设置中，使用的是按Monkhorst-Pack表格在倒格矢空间的划分。

图2-10 K点的设置

适当的选择k点对于达成精确度与效率的平衡是很重要的。预设Monkhorst-Pack点是在给绝缘体的0.1E-1到给金属的0.05E-1之间，这是因为金属系统需要更好的取样。如此通常就能产生足够的点数。例如传统硅晶胞所需的2X2X2，应该检查增加一个，直到建议出来的奇数值Monkhorst-Pack参数能更为有利。我们必须推荐用k点取样的增加来减低有限基底集的修正并促使在一个固定能量下晶体松弛更加精确。设置的方法是在”CASTEP Calculation”中的“Electronic”选项中选择”more”,然后在对话框中选择”k-point”。在”k-point”的选项中又有几种方式。第一种是只取Gamma点，这对于计算体系具有较大的原子数目并且对称性低的情况下可以考虑。第二种是按精度(Quality)来选择，它有三个等级，分别是“course”、”medium”、“fine”。这三个等级对应着不同的Monkhorst-Pack点。比如在做结构优化任务时一般精度都要选择“fine”一等级。第三种是可以给定k点的间隔，这样也就定下了在布里渊区中k点的设置。第四种是直接给定沿着超原胞倒格矢空间三个基矢a、b、c的k点取值，如在图2-10中显示的3X3X10。

(2)结构优化任务的设置
结构优化是CASTEP计算中经常要进行的计算任务，特别是想要计算所关注体系的各种性质的时候，必须首先进行结构优化的计算，在得到结构优化结果文件以后，才能进行性质的计算。所以，正确的设置结构优化的参数是非常重要的。在CASTEP软件中，有四个参数来控制结构优化的收敛参数，它们如表2.2所示：

表2.2
第一个是能量的收敛精度，单位为eV/atom，是体系中每个原子的能量值；第二个是作用在每个原子上的最大力收敛精度，单位为eV/?第三个是最大应变收敛精度，单位为GPa；第四个是最大位移收敛精度，单位为?。这些收敛精度指的是两次迭代求解之间的差，只有当某次计算的值与上一次计算的值相比小于
设置的值时，计算才停止。设置的方法是在”CASTEP Calculation”中选择“Setup”选项，再在”Task”选项中选择“Geometry Optimization”。在”more” 选项中可以进行收敛精度的设置(图2-11)。如下图，我们设置的四个方面的精度分别为5.0e-4eV/atom、0.01eV/?、0.05Gpa、0.001?,如图2-11所示。

图2-11

（3）计算体系性质的设置
在CASTEP中可以计算体系的性质，如能带结构、态密度、聚居数分析、声子色散关系、声子太密度、光学性质、应力等。
在这里我们介绍一下如何进行能带和态密度的计算设置。在计算这两项性质之前，必需先进行自洽计算得到基态能量，而结构优化能够做到这一点，这就是之所以要在计算能带和态密度之前对体系进行结构优化的原因。

图2-12
1）设置能带计算时如图2-12所示。在”CASTEP Calculation”中选择“properties”选项,然后选中”Band structure”,在“Emptybands”中设置好假设的体系导带中取的空带的数目。同时在”K-point set”项里选择计算能带时k点的设置，它也分为“course”、”medium”、“fine”三种。点开“more”,即可看见详细的沿各高对称路径的k点设置。3.0版本中还允许能带计算中的交换-关联函数与自洽计算中的不同，所以在”CASTEP BandStruture Options”中设置了交换-关联函数的选取项。能带中能量的收敛精度在“Band energy tolerance”进行设置。如图2-12中显示的是我们计算能带的设置，空带的数目为16条，k点取了”fine”,交换-关联函数与自洽计算中的相同，带的能量收敛值为1.0e-4eV。

2)态密度的设置与能带类似，它在“properties”选项下有部分态密度(PDOS)的选择项,其k点的设置和在”电子选项“中介绍的完全相同，并且也可以进行与自洽计算中的不同的交换-关联函数，如图2-13所示。
图2-13
3、计算结果的分析
在CASTEP软件中，计算结果从计算服务器上返回以后，在Visualizer界面中就可以进行分析
三．模拟过程与结果
1.模型建立:
点击file,import,导出库中Ge,如图

图3-1
为了减少计算时间,点击build—symmetry—primary cell

图3-2
2.几何优化。

图3-3
点击run，得到优化结构如图

图3-4

比较原来建模时的晶格常数如图3-5

图3-5

比较可知，晶格常数由原来的4.000457A优化为3.947015A。

3.计算热学性质
（1）计算设置
其中交换关联能选择LDA

图3-6

（2）.电子选项设置
K–point 设置为，截断能e-cut为340eV,SCF自洽精度为1.0e-6eV/atom,赝势选择为Norm-conserving模守恒赝势。

图3-7

图3-8

（3).特性设置
在图3-9的左边框图中，点击more，就会出现右边框图，设置如图3-9.

图3-9

（4）.运算设置

图3-10

然后点击run,运行约两个小时,可得到结果，如图3-11.

图3-11
4.分析结果
将目标文件定在 Ge-PhonDOS.cestep.点击Modules—Castep—Ananysis

图3-12
点击view可以等到声子态密度，色散关系，焓，自由能，熵（TS），比热，德拜温度。结果如下图

图3-13

图3-14

图3-15

图3-16

5．结果分析与比较。
分析能量图像可以得知，随着温度T的升高，自由能是下降的，熵是增加的，这就表示，该结构的稳定性随温度增加而提高。由比热容可知，在低温时，与温度呈一次关系，随着温度升高，与温度呈立方关系，高温部分呈水平，与固体物理教材[12]吻合。

参考文献
[1] Kohn W, Sham L J, Self-consistent equations including exchange and correlation effects [J]. Physical review, 1965, 140(4):A1133-A1338.
[2] Hohenberg P, Kohn W. Inhomogeneous electron gas [J]. Physical review, 1964,136(3):B864- B871.
[3] 谢希德, 陆栋.固体能带理论[M].上海：复旦大学出版，1998.
[4] Perdew J P, Chevary J A, Vosko S H. Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchangeandcorrelation[J]. Physical review B, 1992, 46(11): 6671-6687.
[5] Vanderbilt D, Soft self-consistent pseudopotentials in a generalized eigenvalue formalism[J]. Physical review B, 1990,41(11):7892-7895.
[6] Blöchl P E, Projector augmented-wave method[J]. Physical review B, 1994,50(24): 17953-17979. 第87页，共92页
[7]J.P.Perdew and Y.Wang,“Accurate and simple analyticrepresentation of the electron-gas correlation energy”,Phys.Rev.B 45,(1992)13244~13249.
[8]Hammer,B.;Hansen,L.B.;Norskov,J.K."Improvedadsorption energetics within density-functional theory usingrevised Perdew-Burke-Ernzerhof functionals",Phys.Rev.B59,7413-7421(1999).
[9]Vanderbilt,D."Soft self-consistent pseudopotentials in ageneralized eigenvalue formalism",Phys.Rev.B 41,7892-7895(1990).
[10]Ceperley,D.M.;Alder,B.J."Ground State of the ElectronGas by a Stochastic Method",Phys.Rev.Lett.,45,566-569(1980).
[11]Perdew,J.P.;Zunger,A."Self-interaction correction todensity-functional approximations for many-electron systems",Phys.Rev.,B23,5048(1981).
[12]黄昆，韩汝琦.固体物理学. 高等教育出版社.(1985)130-131.