第一性原理计算方法讲义

第一性原理计算方法引言前面讲述的有限元和有限差分等数值计算方法中，求解的过程中需要知道一些物理参量，如温度场方程中的热传导系数和浓度场方程中的扩散系数等，这些参量随着材料的不同而改变，需要通过实验或经验来确定，所以这些方法也叫做经验或者半经验方法。而第一性原理计算方法只需要知道几个基本的物理参量如电子质量、电子的电量、原子的质量、原子的核电荷数、布朗克常数、波尔半径等，而不需要知道那些经验或半经验的参数。第一性原理计算方法的理论基础是量子力学，即对体系薛定额方程的求解。量子力学是反映微观粒子运动规律的理论。量子力学的出现，使得人们对于物质微观结构的认识日益深入。原则上，量子力学完全可以解释原子之间是如何相互作用从而构成固体的。量子力学在物理、化学、材料、生物以及许多现代技术中得到了广泛的应用。以量子力学为基础而发展起来的固体物理学，使人们搞清了“为什么物质有半导体、导体、绝缘体的区别”等一系列基本问题，引发了通讯技术和计算机技术的重大变革。目前，结合高速发展的计算机技术建立起来的计算材料科学已经在材料设计、物性研究方面发挥着越来越重要的作用。但是固体是具有～1023数量级粒子的多粒子系统，具体应用量子理论时会导致物理方程过于复杂以至于无法求解，所以将量子理论应用于固体系统必须采用一些近似和简化。绝热近似（Born-Oppenheimei近似）将电子的运动和原子核的运动分开，从而将多粒子系统简化为多电子系统。Hartree-Fock近似将多电子问题简化为仅与以单电子波函数（分子轨道）为基本变量的单粒子问题。但是其中波函数的行列式表示使得求解需要非常大的计算量；对于研究分子体系，他可以作为一个很好的出发点，但是不适于研究固态体系。1964年，Hohenberg和Kohn提出了严格的密度泛函理论（DensityFunctionalTheory,DFT）。它建立在非均匀电子气理论基础之上，以粒子数密度()r作为基本变量。1965年，Kohn和Sham提出Kohn-Sham方程将复杂的多电子问题及其对应的薛定谔方程转化为相对简单的单电子问题及单电子Kohn-Sham方程。将精确的密度泛函理论应用到实际，需要对电子间的交换关联作用进行近似。局域密度近似（LDA）、广义梯度近似（GGA）等的提出，以及以密度泛函理论为基础的计算方法（赝势方法、全电子线形缀加平面波方法（FLAPW）等）的提出，使得密度泛函理论在化学和固体物理中的电子结构计算取得了广泛的应用，从而使得固体材料的研究取得长足的进步。第一性原理计算方法的应用1、体系的能量。进行第一性原理计算前，首先需要确定体系模型，即模型的晶胞和晶胞中原子的坐标。对于晶体具有周期对称性，具有三个基矢方向和基矢大小（晶格常数）。由于理论计算确定的平衡晶格常数和实验值有所差别，建立模型前需要确定平衡晶格常数。晶格常数的确定采用如下步骤：通过改变三个基矢的大小，改变单胞的体积（81-119%）。通过第一性原理计算可以得到具有不同体积的模型的能量。通过拟合Murnaghan方程，得到晶体的晶格常数以及单胞的能量：'0'000000''00()()11(1)BBVVVEVEVBBBVV其中，0V为基态平衡体积，00()EV为基态下体系的结合能（相对于对应孤立原子能量）。V为原胞体积，0B为体模量，'0B为体模量对压强的导数。如课件中图形所示，可以确定在一定体积下体系的能量达到极小值，即体系的基态能量，所对应的体积为体系的平衡体积，进而可以得到模型三个基矢的大小确定晶体的平衡晶格常数。这里需要指出的是不同的第一性原理计算方法给出的能量，代表的物理意义不同，但是本质上都可以反应体系的稳定性。如总能指构成体系的原子孤立时的能量减去原子成键放出的能量；结合能是以孤立原子的能量为零点，体系具有的总能，即原子构成晶体时放出的能量。在上面求得的晶格常数的基础上，根据要研究的物理问题，确定体系中包含原子数目的多少，建立第一性原理计算模型。第一性原理计算的模型通常选取一个或几个单胞（超单胞）作为模型，选取的模型具有三个基矢方向，应保证沿着三个基矢方向平移可以构成无限大的晶体。第一性原理计算输入的原子坐标有两种坐标形式，一种是笛卡尔坐标（Cartesiancoordinates），一种是分数坐标（fractionalcoordinates）。如对于Ni3Al高温合金，具有如图所示的晶体结构，铝原子位于立方体的顶点，镍原子位于立方体的面心位置。如果取一个单胞作为研究模型，则三个基矢1a，2a，3a分别为（a,0,0）（0,a,0）（0,0,a），其中a为体系的晶格常数。单胞中包含四个不等价原子：三个Ni和一个Al。采用笛卡尔坐标四个原子的坐标可表示为（0，0，0），（a/2,a/2,0），（a/2,0,a/2），（0,a/2,a/2）。如果采用分数坐标表示为（0，0，0），（1/2,1/2,0），（1/2,0,1/2），（0,1/2,1/2）。迪卡尔坐标(x,y,z)和分数坐标(a,b,c)之间关系为1a·a+1a·b+1a·c=(x,y,z)，其中a,b,c为一个原子的三个分数坐标，x,y,z为该原子的笛卡尔坐标。图中所示各点表示将晶格常数的大小a取不同值时得到的单胞体积作为横轴，而纵轴表示对应体积下将原子坐标输入进行第一性原理计算求得的体系能量。拟合后得到Ni3Al的平衡基矢大小以及体系能量。对于研究合金中的掺杂问题，由于掺杂元素的量很少，所以建立的模型需要取多个单胞（超晶胞）。随着模型中原子数目的增加，第一性原理计算方法的计算量指数增加，对于掺杂量很低的情况，如0.1%，需要模型中至少取1000个原子来和实际相符合，这超出了第一性原理计算在目前的计算机上的计算能力（100原子左右），所以建立模型时需要考虑能够反映要研究的实际问题就可以。假设一个超单胞中只存在一个掺杂原子，这样相邻两个超晶胞中掺杂原子的间距为超晶胞的基矢大小。一般两个原子之间相隔三到四个原子层，原子之间的相互作用就可以认为非常小了。所以选取八个单胞构成的超单胞就可以基本反应掺杂量很低的掺杂问题了。将建立好的模型，进行第一性原理计算可以得到体系的总量，对总能进行变换可以定义体系的内聚能、形成能以及择优占位能，进而可以对掺杂是否有利于形成，形成掺杂后对体系稳定的影响而进行分析，如课件中所列。内聚能：体系的总能减去所有原子孤立时的能量，即由于原子之间的相互作用而放出能量，从而内聚能为负值，越小表示形成的体系越稳定。通过和没有掺杂体系的内聚能相比较可以看出掺杂元素对体系稳定性的影响。形成能：体系的总能减去体系中各自元素对应的晶体中原子的能量。形成能可表示各种金属组成合金的能力。另外通过比较掺杂原子替代合金中不同元素原子时体系的形成能可以得到掺杂原子倾向位于合金的什么位置，这个差值就可以定义为择优占据能。这里需要强调，各种能量是根据要研究的具体问题来定义的，比如我们要就掺杂原子倾向于位于合金的什么位置，使用总能是不能得到，因为超单胞模型中各种原子的数目不相同，而每种原子的能量是不一样的，没有可比性，所以定义了择优占位能。2、电荷密度电荷密度就是晶体中电子密度的分布。通过电荷密度可以知道晶体中原子间的成键状态，如金属键、共价键、离子健、vanderWaals键和氢键。为了更好的表示原子形成晶体后原子间的电荷转移和成键情况，引入差分电荷密度，即两个体系中电荷密度的差值。这两个体系应该具有相同的超单胞，超单胞中原子类型可以不一样，而原子位置要基本一致，如课件中所示，用Ni3Al中一个Ni被掺杂元素替代时体系的电荷密度减去没有掺杂的Ni3Al的电荷密度，而得到差分电荷密度，通过图形可以清楚地看出由于掺杂元素的存在导致的电子分布状态的改变。再如，我们可以将Ni3Al晶体的电荷密度减去Ni和Al原子放在超单胞相同位置时孤立Ni和Al原子的电荷密度，可以得到Ni3Al中Ni和Al原子成键过程中电子密度分布的变化，从而更好的观察原子之间的成键情况。所以差分电荷密度相当于一个体系的电荷密度相对于另一个体系的电荷密度改变，目的是为了更好的研究体系中的成键状态。3、能带能带理论是目前研究固体中电子运动的一个主要理论基础。是在用量子力学研究金属电导理论的过程中开始发展起来的。最初的成就在于定性地阐明了晶体中电子运动的普遍性的特点，例如固体为什么会有导体、非导体的区别，晶体中电子的平均自由程为什么会远大于原子的间距等。在半导体技术上，能带论提供了分析半导体理论问题的基础，有利的推动了半导体技术的发展。能带理论是一个近似的理论。在固体中存在大量的电子，他们的运动是相互关联的，每个电子的运动都要受其他电子运动的牵连。价电子是人们最关心的。在原子结合成固体的过程中价电子的运动状态发生了很大的变化，而内层电子的变化比较小，可以把原子核和内层电子近似看成是一个离子实，价电子可以看成在一个等效势场中运动。能带理论的出发点是固体中的电子不再束缚于个别原子，而是在整个固体内运动，称为共有化电子。电子的能量状态从处于一个电子能级变到在一个能量范围内都会存在。我们下面从对自由电子的能量讨论，得出能带的表示方法。自由电子的能量（动能）2/2kEEpm==，根据德布罗依波长与动量关系：/phkl==h，其中h为普朗克常量，h为/2hp，2/kpl=成为波数或波矢。所以电子能量可表示为22/2Ekm=h，即电子能量为波矢k的函数，如课件中图形所示。我们讨论固体中电子的能量通常是在k空间（倒空间）内进行的。对于晶体具有周期对称性，其对应的倒空间也具有周期对称性，对于一维的情况，实空间的周期为a，则倒空间的周期为2/ap，定义第一个周期（/ap-，/ap）为第一布里渊区。则电子能量在实空间分布的随晶体的周期对称性变化；转化到k空间，电子能量随倒空间的周期对称性发生变化，即电子能量在第一布里渊区随k的变化在整个k空间中周期性重复。因此可以得到课件中能量随k变化的图形，相当于自由电子的能带经过周期性势场调制后的结果。所以我们讨论能带只需要考虑电子能量在第一布里渊区随k的变化关系就可以。由于计算所取晶体超单胞形状的不同会导致第一布里渊形状的变化，以及晶体中原子化合状态的不同导致电子所受到的周期性势场不同，所以电子能量随k的变化关系在不同的晶体中是不同的，即不同的晶体具有不同的能带结构，从而反映出不同的物理性质。4、电子状态密度由前面的讨论可知晶体中电子能量状态可以取一定的能量范围，在此能量范围内在不同的能量区间电子能量状态的多少或填充这些能量状态的电子数目是不一样的。电子状态密度反映了这个不同，即在一能量区间内电子状态（数目）的多少，/ddEh，dh为在dE能量区间内电子能量状态的数目。通过对体系电子状态密度的分析可以得到晶体中原子间的电子杂化情况。附：第一性原理计算方法的基本原理第一性原理计算方法的建立是基于对量子力学薛定鄂方程的求解，定态薛定鄂方程为：22[](,,)(,,)2VxyzExyzm-?Y=Yh，其中(,,)xyzY为粒子的波函数，2(,,)xyzY为波的强度，反映了粒子出现在（x,y,z）的概率，即概率密度。V为晶体中粒子受到的势场，E为体系的能量。22ˆH=2Vm-?h为体系的哈密顿量，即能量算符。例：氢原子的薛定鄂方程及其基态解。考虑到氢原子核外电子受核吸引的势能为：204eVrpe=-则氢原子体系的哈密顿算符是：222024eHmrpe=-?h氢原子的薛定鄂方程的直角坐标系的表达式为：222222220()(,,)(,,)(,,)24exyzxyzExyzmxyzrjjjpe抖?-++-=抖?h将薛定鄂方程从直角坐标变换到球坐标，其结果为：22222222201112()(sin)()0sinsin4merErxrrrrjjjqjqqqqfpe抖抖?-++++=抖抖?h由于氢原子核电场是求对称的，故方程解是球对称的，即：0jq¶=¶；0jf¶=¶在这种条件下，氢原子的薛定鄂方程简化为：2222012()()04merErxrrjjpe