计算材料计算BN的弹性常数

湖南工业大学
课程设计
资料袋
理学院（系、部）2011 ~ 2012 学年第一学期
课程名称计算材料学指导教师雷军辉职称讲师
学生姓名余晓燕专业班级应用物理081班学号***********
题目计算BN的弹性常数
成绩起止日期2011年12月4日～2011年12 月12 日
目录清单
1
湖南工业大学
课程设计任务书
2011—2012 学年第1 学期
理学院学院（系、部）应用物理学专业081 班级课程名称：计算材料学
一、设计题目：计算BN的弹性常数
指导教师（签字）：年月日
系（教研室）主任（签字）：年月日
2
（计算材料）
设计说明书
计算BN的弹性常数
起止日期：2011 年12月4日至2011 年12月12日
学生姓名余晓燕
班级081
学号***********
成绩
指导教师(签字)
理学院（部）
2011年12 月12 日
3
计算BN的弹性常数
背景：
近年来，随着材料、物理、计算机和数学等学科的发展，应用计算的方法研究材料的结构、能量和性能已成为一门迅速发展的新兴学科-计算材料学。

这种方法不仅能进行材料的计算模拟，而且能进行材料的计算机设计和相关性能的预测。

随着计算机技术的飞速发展，第一性原理计算的方法在材料的结构和性能等方面的研究已取得了巨大的成功，第一性原理的方法是基于量子力学理论，从电子运动的层次研究材料的结构和相关性能。

目前，CASTEP软件的主要功能是对半导体、非线性光学材料、金属氧化物、玻璃、陶瓷等固体材料，对电子工业、航空航天以及石化、化工等工业领域有着非常重要的战略意义。

对这些材料而言，其电子的结构与性质，以及表面和界面的性质与行为都非常重要。

CASTEP的量子力学方法，为深入了解固体材料的这些性质并进而设计新的材料，提供了强有力的工具。

基于密度泛函平面波赝势方法的CASTEP软件可以对许多体系包括像半导体、陶瓷、金属、矿石、沸石等进行第一性原理量子力学计算。

典型的功能包括研究表面化学、能带结构、态密度、热学性质和光学性质。

它也能够研究体系电荷密度的空间分布和体系波函数。

CASTEP还可以用来计算晶体的弹性模量和相关的机械性能，如泊松系数等。

半导体和其他固体材料的许多性能由电子性质决定，而电子性质又由原子结构决定，特别是缺陷在改变电子结构上的作用对半导体性质尤为重要。

分子模拟，特别是量子物理技术，可用来预测原子和电子结构及分析缺陷对材料性能的影响。

CASTEP能有效的研究存在点缺陷、空位、替代杂质、位错等的半导体和其它材料中的的性能。

除此以外，它还可以被用来计算固体的振动性质，如声子色散关系、声子态密度等。

这些计算结果可以用来分析表面吸附的振动性质，可以解释实验中的振动谱，可以研究在高温高压下的相稳定性等等。

总的来说，它可以实现如下的功能：
1.计算体系的总能；
2.进行结构优化；
3.执行动力学任务：在设置的温度和关联参数下，研究体系中原子的运动行为；
4.计算周期体系的弹性常数；
5.化学反应的过度态搜索。

除此之外，计算一些晶体的性质，如能带结构、态密度、声子色散关系、声子态密度、光学性质、应力等。

下面介绍一下密度泛函理论、交换关联泛函近似、赝势方法和K-S方程迭代解法。

一、基础理论：
4
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]，这个定理包含如下内容：不
计自旋的情况下，将粒子数密度函数表示成全同费米子系统的基态能量的唯一泛函；
在粒子数不变的情况下，能量泛函对正确的粒子数密度取等于基态能量的极小值。

因此，对于基态非间并多粒子系统，不考虑自旋的条件下，其哈密顿算符为
5
(1-8)
式(1-8)中，外场作用看成原子核-电子作用，相同的局域势对外场的作用用表示。

对于给定的外场，多电子系统的能量表示成电子数密度的泛函为：
(1-9)
(1-10)
(1-11)
(1-12)式中，包括体系中电子之间的相互作用能和电子的动能，是外场对电子的作用能，是系统中原子核间的排斥能。

在式(1-10)中，前两项表示无相互作用粒子模型的动能和库仑排斥能，复杂的电子相互作用用交换关联能表示。

根据Hohenberg-Kohn 定理，假设能得到能量泛函E(ρ )，然后就能将电子数密度ρ 变分，就能确定系统的基态和基态所有的性质，因此确定E(ρ )成为问题的关键所在，而要确定能量泛函E(ρ )，必须要确定动能泛函T[ρ ]、电子数
密度 以及交换关联泛函。

为了解决上述问题，W.Kohn 和L.J.Sham 提出了如下假设：假定已知无相互作用的电子系统和未知的有相互作用的电子系统密度函数相同，未知的相互作用电子系统的动能泛函T[ρ ]可用已知的无相互作用电子系统的动能泛函来代替；假定密度函数 用N 个单电子波函数构成，于是有：
(1-13)则
(1-14)对能量泛函进行变分得到
6
(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 为首的，这类
7
泛函包含若干个实验参数，这些参数通过计算和实验数据来获得，这种形式的好坏由实践的工作所决定的。

另一个是Perdew-Wang 91 的，这类泛函以物理规律为基础，不包含实验参数。

随着研究的不断深入，不仅出现了非局域的相互作用交换关联泛函，还有密度高阶梯度的近似交换关联泛函，如Vaner Wals 和Meta-GGA 等。

GGA 和LDA 相比在能量精确度和开放体系方面更有优势。

3. 赝势法：
在晶体的近自由电子能带计算中，计算量大而且收敛速度慢。

对于固体而言，价电子的化学性质活泼，对于结构和性质的影响较大，而内层电子的能带较窄，较稳定，而且相邻原子的作用对内层电子的状态影响较小。

因此，人们关注的是价电子，将原子核和内层电子近似看出粒子实。

对于固体中的价电子波函数而言，在离子实的内部区域，变化剧烈，存在若干个节点；而在离子实之间的区域，变化平缓。

离子实内部的这一特点要求价电子波函数与内层电子波函数正交，而价电子与内层电子波函数正交起了一种排斥势的作用，在很大程度抵消了离子实内部V(r )的吸引作用。

据此，离子实内部的势函数用假想势代替，在离子实之间的区域波函数和电子的能量本征值保持不变的条件下求解固体单电子波函数方程，假想的势叫赝势，用赝势求出的波函数叫赝势波函数。

对于多原子固体而言，根据波函数的不同特征坐标空间被分成c r 以内的原子核区域(芯区)和以外的其它区域两部分(假定存在某个截断距离c r )。

芯区(r< c r )其波函数与紧邻原子波函数相互作用很小，赝势和赝势波函数变化缓慢，比较平坦；芯区外(r> c r )价电子波函数相互交叠作用，和真实的势和波函数相比，其形状和幅度都一样。

目前，除了经验赝势、半经验的模型赝势外，还有没有附加经验参数的第一性原理从头算原子赝势，包括模守恒赝势(norm-conserving potentials NCP)、超软赝势方法(ultra-soft pseudo-potentials，USPP)[5]
5.分子轨道的自洽求解：
(1)分子轨道的自洽场方程：
密度泛函理论是基于Hohenberg--Kohn定理，该定理表明体系基态的性质由电荷密度决定，体系的总能量是电荷密度ρ的函数。

总能Et可以表达为：
(1-20)
T[ρ]是密度为ρ的电子的动能，U[ρ]是经典的库仑相互作用静电能，Exc[ρ]包括了多体相互作用对总能的贡献，其中交换-关联能是主要的部分。

我们从波函数Ψ来构造电荷密度。

对于波函数Ψ可以写成具有反对称性单?粒子波函数（分子轨道）的Slater行列式：
(1-21) 当分子轨道是正交时，即
(1-22) 电荷密度可表示为：
(1-23) 8
由总能的表达式和电荷密度的表达式，动能项（原子单位）可表示为：
(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，则出现下图：
9
10 图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点取样精度、截断能等的设置时都要进行适当的选取。

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赝势一般适用于金属体系，列如我们设置Ultra-soft赝势。

设置的方法是在”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）。

e、K点设置
布里渊区抽样的设置：
布里渊区的设置是通过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/A第三个是最大应变收敛精度，单位为GPa；第四个是最大位移收敛精度，单位为A。

这些收敛精度指的是两次迭代求解之间的差，只有当某次计算的值与上一次计算的值相比小于设置的值时，计算才停止。

设置的方法是在”CASTEP Calculation”中选择“Setup”选项，再在”Task”选项中选择“Geometry Optimization”。

在”more”选项中可以进行收敛精度的设置(图2-11)。

如下图，我们设置的四个方面的精度分别为5.0e-4eV/atom、0.01eV/A、0.05Gpa、0.001A,如图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. 优化BN立方晶体的结构
在计算弹性常数之前并不一定要进行几何优化，可以由实验观测到的结构计算出Cij数据。

尽管如此，如果我们完成晶胞的几何优化，可以获得更多相容的结果，进而得到理论基态对应的弹性常数。

弹性常数的精确度，尤其是切变常数的精确度，主要取决于SCF计算的品质，布里渊区取样和波函数收敛程度的品质。

所以我们设置SCF、k点取样和FFT格子的精度为Fine。

BN结构的建立
在菜单栏中选择File/ Import，从structures/semiconductors中选中
BN.msi，按Import按钮，输入BN的晶体结构，见下图。

为了节省计算时间，由Build / Symmetry / Primitive Cell将此conventional representation 转化为primitive representation.
几何优化的设置
从工具栏中选择CASTEP工具，然后从下拉列中选择Calculation（或从菜单栏中选择Modules / CASTEP / Calculation）。

CASTEP Calculation对话框见下图：
在Setup标签中，把Task设置为Geometry Optimization。

按下more按钮，选中Optimize cell。

关闭CASTEP Geometry Optimization 对话框。

选择Electronic标签，按下More...按钮以得到CASTEP Electronic Options 对话框。

把Derived grid的设置从Standard改为Fine。

关闭CASTEP Electronic Options对话框。

选择Job Control标签，设定本地机运算。

按下CASTEP Calculation对话框中的Run按钮。

结果如下图。

2. 计算BN的弹性常数
BN CASTEP GeomOpt/BN.xsd处于激活状态。

选择CASTEP Calculation对话框中的Setup标
签，从Task的下拉清单中选择Elastic
Constants。

按下More...按钮，CASTEP Elastic
Constants对话框见右图。

将Number of steps
for each strain由4增加为7，按Run运行。

CASTEP的弹性常数计算任务的结果以一
批.castep输出文件的形式给出。

这些文件中的
每一个文件都代表确定的晶胞在假设的应变模
式和应变振幅下的几何优化运行结果。

这些文件
的命名约定为：seedname_cij__m__n。

对于给定
的模式来说，m代表当前的应变模式，n代表当前
的应变振幅。

CASTEP可以使用这些结果来分析每一个运行计算出来的压力张量，产生一个有关弹性性质的文件。

从工具栏中选择CASTEP工具，然后选择Analysis或者从菜单栏中选择Modules | CASTEP | Analysis。

从属性清单中选择Elastic constants，从BN的弹性常数计算工作中得到的结果文件BN.castep应自动显示在Results file选框中。

按下Calculate按钮。

计算结束后产生一个新的文档BN Elastic Constants.txt。

此文档中的信息包括：
*输入的应变和计算出的应力的总结
*每一种应变模式线性拟合和拟合质量的计算结果
*给定对称性下计算出的应力与弹性常数之间的对应
*弹性常数Cij和弹性柔量Sij的表格
*导出量：体积模量和其倒数、压缩系数、杨氏模量、Poisson比、 Lame 常数(用于模拟各向同性介质) 。

3 结果分析
1）.弹性刚度常数：
2）.弹性柔度常数
3）.体弹模量：
4）.杨氏模量和泊松比：
四．心得体会：
一学期的学习结束了，通过这一学期的学习我了解到，计算材料学是材料科学与计算机
科学的交叉学科，是一门正在快速发展的新兴学科，是关于材料组成、结构、性能、服役性能的计算机模拟与设计的学科，是材料科学研究里的“计算机实验”。

它涉及材料、物理、计算机、数学、化学等多门学科。

正是由于它涉及到的学科之广和意义之重大让我们在拿到课题时真的有些无从下手，真正的感觉到了理论和实践的差距。

还好我们知道团队的力量是伟大的，我们很幸运的拥有了一个优秀的团队。

团员们各取所长，分工明确，我主要负责相关资料的搜集和最后的文档整理。

面对复杂的课题，有了这只强大的团队，我们很快就有了头绪，也迈出了课程设计的第一步。

千里之行始于足下，有了一个很好的开始我们接下来的设计进行的也算顺利。

但是，面对第一次接触的软件Materials studio，对各种功能的熟悉，我们还是遇到了不少的麻烦也走了不少歪路。

这时，雷老师总能来到我们身边进行细心指导和说明，为我们解决并讲解一个个问题。

真的很感谢雷军辉老师的帮助让我们能按时并完美的完成课程设计。

路漫漫其修远兮，吾将上下而求索。

未来的路还很长，我相信以后的路上我还会遇到很多材料方面的问题，但我更相信我们定能学以致用，一定不辜负雷老师的期望和谆谆教诲！
参考文献
[1] Kohn W, Sham L J, Self-consistent equations including exchange and
correlation effects [J]. Physical review, 1965, 140(4):A1133-A1338.
21。