密度泛函理论.docx

密度泛函理论.docx

《密度泛函理论.docx》由会员分享，可在线阅读，更多相关《密度泛函理论.docx（16页珍藏版）》请在冰点文库上搜索。

密度泛函理论

密度泛函理论

摘要：

介绍了密度泛函理论的发展与完善，运用密度泛函理论研究了钒（Vanadium）在高压下的结构相变。

通过计算体心立方结构的钒在不同压强下剪切弹性系数C44,发现当压强约95GPa时C44<0,说明体心立方结构的钒在此条件下是不稳定的。

进一步计算分析得到钒在高压下发生了从体心立方到菱面体的结构相变,相变压强约70GPa,这一结果与实验结果符合。

还首次发现当压强约380GPa时,将会发生菱面体到体心立方的结构相变,这有待实验的验证。

引言：

相变的研究受到广泛重视,通过相变研究可以认识物质的内部结构,可以了解原子核的内部性质。

尤其是极端条件下—高温、高压下相变的研究一直是人们关注的热点,能量很高的重离子反应能形成高温、高密的区域,在这种条件下会出现许多奇异现象[1]。

原子在高压下也会出现许多新的特征,如发生结构相变。

过渡金属钒由于有较高的超导转变温度Tc,最近成为实验和理论研究的主题[2—8]。

Ishizuka等[2]对钒的实验研究发现:

常压下钒的转变温度Tc为5.3K,并随压强成线性增长的关系,当压强为120GPa时Tc=17.2K（迄今是金属中最大的Tc）,但压强大于120GPa,Tc出现了反常,即不再随压强成线性增长而保持不变。

Takemura等[8]对高压下的钒进行了X射线衍射实验,结果显示状态方程并没有奇异性,体心立方结构的钒在压强达到154GPa时仍是稳定的。

Suzuki和Ostani利用第一性原理对进行了计算,发现横向声子模在加压下有明显的软化,当压强约130GPa时变成虚的,能说明可能发生了结构相变,但并未给出相变细节[3]。

Nirmal等[4]理论计算表明,压强约140GPa时会发生体心立方到简立方（sc）的结构相变。

Landa等[5,6]计算了体心立方结构的钒在加压下剪切弹性系数C44的大小,发现压强约200GPa时会出现力学不稳定,并用费米面嵌套解释了不稳定的原因,但并没有给出相变后的结构。

最近Ding等[7]在常温下首次从实验上得到当准静压约63GPa时钒会发生从体心立方到菱面体的结构相变,并分析了产生结构相变的原因。

他们认为,排除传统的s-d电子跃迁的驱动,相变可能与来自于费米面嵌套、带的Jahn-Teller扭曲以及电子拓扑跃迁等因素有关。

基于如上原因,本文运用密度泛函理论研究钒在高压下的结构相变,即通过计算体心立方结构的

电子坐标用

表示;离子实坐标用

表示。

若不考虑其他外场的作用,哈密顿量表示如下:

其中

和

分别表示电子体系和离子实的哈密顿量,为电子与离子实的相互作用能。

电子体系哈密顿量为:

其中

表示电子动能项;

表示库企相互作用势能项。

对

除外的所有电子进行求和。

离子实的哈密顿量为:

第一项表示离子实的动能;第二项表示离子实间的相互作用能。

对除

外的所有的离子实求和,

表示第

个离子实的质量。

这里假定离子实间的作用能与离子实之间的位矢差

有关。

电子与离子实的相互作用能形式上为:

式（2.1.1）~（2.1.5）中每立方米对/和/求和的数量级为1029,所以需要进一步的近似。

在电子体系的哈密顿量

中,只出现电子坐标

;在离子实的哈密顿量

中,只出现离子实坐标

;电子坐标和离子实坐标及同时出现在离子实与电子相互作用项中。

由于他们的作用都是一个数量级的,所以不可以简单略去。

离子实的运动幅度和运动速度都比电子要小得多,而且离子实的质量约为电子质量的1000倍以上,可认为离子实处于某一平衡位置做轻微的振动,电子则在晶体内大范围运动。

那么就可以分两部分来处理问题:

处理电子的运动的时候认为离子实是于某处静止不动,而处理离子实运动的时候则不需考虑电子的空间分布,这就是绝热近似也称为玻恩-奥本海默近似[32]。

系统薛定愕方程的解即体系波函数可以通过绝热近似表示为离子实波函数和电子波函数的乘积:

其中

表示原子核运动的波函数,

为多电子系统的波函数。

在绝热近似方法的基础上,在考虑电子运动时,认为离子实是在平衡位置不动,忽略其动能,选取离子实所在位置为零势能点,可得多电子体系的薛定愕方程

2.1.2哈特利-福克（Hartree~Fock）近似

通过绝热近似得到可以得到多电子体系薛定愕方程（2.1.7）,式中电子之间的相互作用仍然无法精确求解,哈特利提出假设可以忽略电子之间的相互作用,该多电子体系的薛定愕方程表示为:

其电子波函数为柳的连乘积:

令

可得单电子方程为:

实际上电子与电子间的相互作用是不能忽略的,多电子薛定愕方程（2.1.7）的近似解仍然可以用单电子波函数乘积式（2.1.9）来表示。

可以根据变分原理通过总能量对每个单电子波函数进行变分使被忽略电子间相互作用的多电子系统的能量更接近有电子间相互作用的多电子系统的能量,继而得到哈特利方程:

该方程描写的是在晶格势

和其他全部电子平均势场中,位于坐标

处的单电子的运动。

此时体系总的波函数仍为:

根据泡利不相容原理,每个电子具有不同的量子态。

在满足电子交换反对称性的基础上,可将总的波函数看做是所有单粒子波函数的线性组合。

通过进一步的变分处理,可得哈特利-福克方程:

此式中;表示交换电荷密度,即考虑了电子的交换反对称性。

但是其中没有包括由于自旋反平行导致的电子间的排斥相互作用,也就是说电子的关联作用没有包含在内。

2.1.3霍恩伯格□科恩（Hohenberg-Kohn）定理

霍恩伯格-科恩定理是基于非均匀电子气理论提出的:

定理一为不计自旋的全同费米子系统的基态能量是粒子数密度函数的唯一泛函[34];定理二为能量泛函在粒子数不变条件下对正确的粒子数密度函数取极小值,并等于基态能量[34]。

如不考虑自旋,该系统的电子哈密顿量写为:

其中

为动能项,

为库企排斥项,

表示外势场影响。

定理一说明多粒子系统的基态波函数都是由密度函数决定,进而可以得到体系所有的电子性质。

定理二说明如果可以得到系统的基态密度函数,然后利用变分原理的计算,通过能量泛函的最小化就可以找到准确的系统基态能量。

基态的电子密度分布函数的泛函

是未知的。

为进一步进行说明,将无相互作用粒子的对应项从中分离出:

无相互作用的粒子模型其动能用

表示,

表示的是库合排斥项。

交换关联项用

来表示,代表所有没有包含在无相互作用的粒子模型中的、相互作用的项,包含相互作用的全部复杂性。

在这里仍然是未知,因为它是P的泛函。

上式（2.1.16）中还要进一步明确:

粒子数密度函数

、动能的泛函

和交换关联能的泛函

。

2.1.4Kohn-Sham方程

以Hohenberg-Kohn定理为基础,根据粒子数不变,通过能量泛函对密度函数进行变分,就能得出基态的能量和粒子数密度函数:

上式中的

具有化学势的意义,称作拉格朗日乘子。

（2.1.17）式描述的是粒子处于一有效势场中的形式,表示为:

但是,动能的泛函

还是没有解决。

Kohn和Sham于1965年提出有效方案用于解决有相互作用粒子的动能项:

定义一个已知粒子的无相互作用动能泛函,

用来替代复杂的有相互作用粒子的动能泛函

这样无相互作用粒子系统就与复杂的相互作用粒子系统的粒子密度数相同[34]。

定义一个交换关联能泛函,

将有相互作用粒子的动能泛函和无相互作用粒子的动能泛函相比较复杂的、无法转换的部分归入其中,然后用

个单粒子的波函数

构成整体的密度函数,完成单粒子图像

用对

的变分代替对

的变分,用五,代替拉格朗日乘子

可得:

得到对于单粒子波函数的单电子方程。

以上三式统称Kohn-Sham方程。

根据Hohenberg-Kohn定理,将该系统的波函数、基态的能量以及各物理算符的期待值等由求出的粒子数密度函数准确地得出,单粒子能量为得到的本征值。

交换关联近似

应用Kohn-Sham方程,能用有效的单电子问题描述有效的多电子系统基态特性问题,但是由于在交换关联能泛函仍然包含着多粒子系统相互作用的所有复杂性,因此解Kohn-Sham方程,必须找到交换关联能泛函。

目前比较常用的交换关联能泛函是局域密度近似（LDA）和广义梯度近似（GGA）,前者适用于均勾电子电子体系,后者适用于非均匀电子体系[34]。

局域密度近似是对交换能最简单的近似,其主要思想是对于一个密度函数分布平坦的体系,其电子密度的空间变化非常缓慢,其空间可以分割成许多无穷小的体积元,可以认为电子在每个体积元中分布是近似均匀的,这样,各小体积元对交换关联能的贡献可视为与体积相同的均匀电子气对交换关联能的贡献相同[38]。

对于电荷密度变化不是很明显的体系,局部被视为均勾电子气,局域密度近似可取得较好的效果,在凝聚态物理及量子化学中的应用已十分广泛。

但是,局域密度近似也存在很多缺点,例如在处理半导体材料带隙、有电子强关联效应的材料和超大体系研究时仍然存在问题。

为克服局域密度近似存在的缺点,在此基础上引入了对粒子数密度的梯度修正,大幅度减少了局域密度近似的误差,成为广义梯度近似,与局域密度近似相比,广义梯度近似全面的考虑到了每个区间内电子数密度对交换关联能的贡献和相邻区间的不同电子数密度对交换关联能的贡献,更适合处理具有非均匀密度的电子体系。

但广义梯度近似亦有缺点,如高估晶格常数,还有无论是局域密度近似还是广义梯度近似都会低估绝缘体和半导体的带隙。

布洛赫（Bloch）定理

实际计算中,在不能精确求解的无穷大体系中抽取具有代表性的一部分作为研究对象,或把无穷大的体系用有限体系模型模拟出来,这样既能达到使计算简化的目的,又能体现该无穷大体系的整体性质。

通常是创建一个有限大小的基本单元,使它在整个空间中进行周期性的繁衍来填充整个空间,形成具有周期性的宏观物质体系。

根据周期性边界条件,宏观物质体系具有与基本单元相同的周期性,可以通过久期方程在第一布里渊区进行求解。

在描述周期性体系中的电子行为时,布洛赫电子波是一个十分重要的概念。

布洛赫（Bloch）定理说明了周期性的体系其电子波函数具有“调幅波”的性质[33,34]。

晶体电子哈密顿算符的本征函数为布洛赫（Bloch）电子波函数。

在周期性晶格势场中,用有限大小的体系来描述无限大的体系,把体系中拥有相同矢量左的电子归于同一本征空间中的不同本征态,该体系本征空间的数目与简约波矢的数目相同（事实上本征空间的数目可能会略小于简约波矢的数目）。

所以凝聚态物理学中把i称作量子数,有时直接当成一个与体系的宏观周期性有关的“量子数”。

事实上i并不是实际的动量访2,但又具有动量的属性,所以在描述周期单元中电子运动时{被称为准动量,或晶体动量。

2.3赝势（Pseudopotential）方法

在固体材料的原子核附近,周期性势场由于电子受到库企相互作用远离平均值,.电子波函数因受原子核的势场的作用,波函数发生明显的变化,这个周期性势场不可以被看作是微扰势场而被忽略,这就需要用势（Pseu（k）potential）方法来解决。

求解薛定愕方程时,不改变离子实间区域的波函数和能量本征值,离子实内部真实的势能用假想势能来代替,这个假想的势能就叫做赝势。

对于晶体来说,它是由大量原子构成的,可以根据波函数的特点将其坐标空间分为以下两个部分:

1、芯区处于截断距离r。

以内,是靠近原子核的区域,其波函数主要是由紧束缚的芯电子波函数构成,它与邻近的原子波函数的相互作用很微弱。

2、其余区域处于截断距离r。

以外,是除芯区以外区域,价电子的波函数是相互交叠和相互作用的。

价电子被芯区的势强烈吸引,然而在正交化平面波方法中需要价态和芯态正交因此产生了极大动能,这个动能对价电子起着一种排斥的作用,在很大程度上它抵消了在芯区周期性势场的作用。

赝势思想就是在原子的芯区取截断距离r。

为了使离子实周围的势场变得平缓,需要部分地抵消原子核的吸引势,也就是在半径之内球体中构造短程排斥势,从而有利于平面波的展开。

在构造赝势的时候,用一个趋于平滑的波函数来描述振动剧烈的芯区,并且不改变芯区内电荷密度的分布。

赝势是一种势函数,它给出的波函数在其余区域（即芯区以外的区域,具有和真

实波函数一样的形状。

在展势的应用中赝势的硬度是一个十分重要的概念。

当一个赝势用很少的傅立叶变换组成就可以准确描述时,这个赝势就称为软赝势,反之成为硬赝势。

超软展势是在上世纪90年代提出的软，芯区的赝平面波函数尽可能的“软”,从而截断能就可大幅度的减少。

模守恒赝势[4i】有以下的特点:

其相对应的波函数和真实势对应的波函数的能量本征值相同;截断距离r。

以外的其余区域和真实波函数的幅度和形状相同,也就是“模守恒”;截断距离r。

以内电子波函数和赝势波函数具有相同的微分函数,变化缓慢,没有大的动能。

其主要缺点是模守恒赝势需要的截断动能比较高。

基于超软赝势结合线性缀加平面波思想提出的PAW方法,经推导验证得到了广泛的应用。

第一性原理计算方法

第一性原理计算又称为从头计算，是从电子结构出发,只需使用某些基本的常量和一些合理的近似来进行计算,从而分析物质的结构及性质的一种方法,它是计算材料科学中十分重要的研究方法。

第一性原理的基本思想是把固体看成是由原子核和电子共同组成的多粒子系统,应用量子力学的方法对系统内的粒子的运动进行“非经验”计算,也就是在计算.时只需要确定构成微观体系各元素的原子类型及原子位置,进行理论模拟计算时不需要任何其他参数（包括经验、拟合等）,计算结果可以得到微观体系的总能量、电子结构等性质来预测宏观物质的性质[31]。

第一性原理计算的优势在于,一方面它可以和真实的实验相互补充、相互促进;另一方面与真实的实验相比,它能更快地设计出符合要求的实验及材料。

现如今,计算机软件和硬件的飞速发展,在研究领域中计算材料科学与理论和实验具有同等重要的地位。

第一性原理计算,特别是基于密度泛函理论的计算,近年来己成为计算材料科学的基础和核心,在新材料的开发和设计等诸多领域有广泛的应用。

运用密度泛函理论研究了钒（Vanadium）在高压下的结构相变。

相变的研究受到广泛重视,通过相变研究可以认识物质的内部结构,可以了解原子核的内部性

质。

尤其是极端条件下———高温、高密（或高压）下相变的研究一直是人们关注的热点,对能量很高的重离子反应能形成高温、高密的区域,在这种条件下会出现许多奇异现象[1]。

原子在高压下也会出现许多新的特征,如发生结构相变。

过渡金属钒由于有较高的超导转变温度Tc,最近成为实验和理论研究的主题[2—8]。

Ishizuka等[2]对钒的实验研究发现:

常压下钒的转变温度Tc为5.3K,并随压强成线性增长的关系,当压强为120GPa时Tc=17.2K（迄今是金属中最大的Tc）,但压强大于120GPa,Tc出现了反常,即不再随压强成线性增长而保持不变。

Takemura等[8]对高压下的钒进行了X射线衍射实验,结果显示状态方程并没有奇异性,体心立方结构的钒在压强达到154GPa时仍是稳定的。

Suzuki和Ostani利用第一性原理对钒进行了计算,发现横向声子模在加压下有明显的软化,当压强约130GPa时变成虚的,能说明可能发生了结构相变,但并未给出相变细节[3]。

Nirmal等[4]理论计算表明,压强约140GPa时会发生体心立方到简立方（sc）的结构相变。

Landa等[5,6]计算了体心立方结构的钒在加压下剪切弹性系数C44的大小,发现压强约200GPa时会出现力学不稳定,并用费米面嵌套解释了不稳定的原因,但并没有给出相变后的结构。

最近Ding等[7]在常温下首次从实验上得到当准静压约63GPa时钒会发生从体心立方到菱面体的结构相变,并分析了产生结构相变的原因。

他们认为,排除传统的s-d电子跃迁的驱动,相变可能与来自于费米面嵌套、带的Jahn-Teller扭曲以及电子拓扑跃迁等因素有关。

基于如上原因,本文运用密度泛函理论研究钒在高压下的结构相变,即通过计算体心立方结构的钒在不同压强下的剪切弹性系数C44得到结构相变的信息。

研究得到钒在高压下发生从体心立方到菱面体的结构相变,相变压强为70GPa;当压强约380GPa时,又会发生菱面体到体心立方的结构相变,这需要实验的论证。

计算前对几何结构进行优化,常压下钒为体心立方结构,优化得到的体心立方结构钒的平衡晶格常数a0,体模量B0列在表1中。

从表中可看出,计算得到的a0和B0与实验理论符合得较好。

说明我们的优化方法是合适的。

表1　常压下体心立方结构的钒的晶格常数a0和体模量B0

图1给出了体心立方结构钒约化体积V/V0随着压强P变化的关系图,V0=13.923为常压下实验测得的单胞的体积。

从图中可以看出,随着压强的增大,约化体积是逐渐减小的,减小的程度随压强的增大而变小。

图1约化体积V/V0随压强P的变化

对钒进行几何结构进行优化后,我们将讨论加压下体心立方结构钒的结构相变。

首先计算剪切弹性系数C44=ω2ρ/K2,其中ω,ρ,K分别为声子频率、密度和声子波矢。

在长波极限下,与体心立方结构的三角形扭曲相联系。

由压强导致钒的结构相变可通过声子软化引起C44的不稳定性来分析[7]。

常压下C44=18,随着压强的增大而增大,压强约为60GPa时,C44达到最大值,之后随着压强的增大而减小;压强约为95GPa时,C44=0;压强再增大时,C44<0（如图2所示）。

图2体心立方结构的钒的剪切弹性系数C44随压强P的变化

图2给出了体心立方结图2体心立方结构的钒的剪切弹性系数C44随压强P的变化构的钒在不同压强下的剪切弹性系数C44,所得结果与Landa的C44结果定性的一致,定量的差别可能是选的交换势不同,本文选的交换势为GGA而Landa选的局域密度近似LDA[12]。

一个晶体的力学稳定性要求它的弹性能必须是正的,对于一个立方晶体要求C44>0[13]。

因此,可以得出95GPa是C44不稳定的上限,超过这一压强,体心立方结构的钒中一个有限的剪切会导致无限大的张力而造成结构的崩塌,在到达这一极限以前将发生结构相变[7]。

为了验证体心立方结构的钒在高压下的结构相变,我们通过计算来讨论是否存在体心立方到简立方的结构相变。

体心立方结构和简立方结构的钒的总能量Etol随体积的变化如图3所示。

体积相同时,钒的体心立方结构的总能量始终小于简立方结构的总能量,即体心立方结构始终比简立方结构稳定,说明钒在高压下不会发生从体心立方到简立方的结构相变。

简立方结构计算得到的C44=-121GPa（压强P=0）,C44=-276GPa（压强P=100GPa）,这也说明简立方结构的钒在常压和加压下都是不稳定的。

图3bcc和sc结构的钒在不同体积下的总能量

　　钒菱面体（Rhombohedral）结构空间群为R3-m,是一个三角Bravais晶格,用三角晶格单胞的棱长和两个三角轴的夹角αr表征。

当αr=109.47°时,菱面体结构变为特殊的体心立方结构;αr=90°时,为简立方结构;αr=60°时,为面心立方结构（fcc）。

吉布斯自由能为G=Etol+PV+TS。

因为理论计算中T=0K,所以吉布斯自由能等G于焓H=Etol+PV[14]。

为了验证钒是否发生体心立方（bcc）到菱面体（Rhombohedral）的结构相变,我们计算了菱面体结构（αr=109.65°）和体心立方结构的焓差ΔH=HR-Hbcc,如图4所示。

ΔH=0时,压强P≈90

GPa;当P<90GPa时,ΔH>0;而P>90GPa时,ΔH<0。

所以,体心立方到菱面体（αr=109.

65°）的结构相变的相变压强约为90GPa,与实验测得的80GPa时钒从体心立方结构转变为αr=109.65°的菱面体结构较符合[7]。

图4不同压强下菱面体（αr=109.65°）结构和体心立方结构的ΔH焓差ΔH

另外,还详细计算了不同约化体积（或压强）下,钒的总能量最低时对应的αr的变化,来研究结构相变,如图5所示。

图5中右纵坐标为钒在不同约化体积V/V0下,所有的αr中对应的最低总能量Etol,其中V0=13.923为常压下实验测得的单胞的体积,V是计算所取的体积。

左纵坐标为不同约化体积V/V0下最低的总能量Etol对应的αr。

总能量随着图5钒在不同约化体积V/V0下最低的总能量Etol和对应的αr约化体积的减少先减少后增大,即有一极小值对应着基态,此时约化体积V/V0=0.95,αr=109.47°,与实验测得的基态较符合。

αr随着约化体积的变化分为4个部分:

（1）V/V0（1.2※0.95）,αr从112°减少到109.47°,这一部分并不对应常压或加压下的情形,因为计算的基态约化体积V/V0为0.95;

（2）V/V0（0.95※0.75）对应压强约（0※70GPa）,αr保持109.47°不变,即为体心立方结构;（3）V/V0（0.75※0.5）对应压强约（70※380GPa）,这一范围的αr都偏离109.47°,即不再是体心立方结构而是菱面体结构;（4）V/V0（0.5※0.4）对应压强约

（380※940GPa）,αr保持109.47°不变,即为体心立方结构。

钒在加压下发生两次结构相变:

在V/V0=0.75对应压强约为70GPa时,发生体心立方到菱面体的结构相变,与实验结果较符合[7];当V/V0=0.5对应压强约为380GPa时,发生了菱面体到体心立方的结构相变,与Te发生菱面体到体心立方的结构相变相似[14,15],该结果尚没有实验对比,有

待实验的进一步检验。

本文运用密度泛函理论研究了钒在高压下的结构相变,体心立方结构的钒在不同压强下计算得到的剪切弹性系数C44表明,外界压强约为95GPa时C44<0,即体心立方结构的钒是不稳定的。

通过计算发现,钒在高压下不会发生体心立方到简立方的结构相变,而是发生体心立方到菱面体的结构相变,相变压强约为70GPa,与实验得到的结果符合较好。

这将对钒在高压下的超导转变温度T0的研究提供重要信息。

当压强约为380GPa时,将会发生菱面体到体心立方结构的相变,这有待实验的进一步检验。