主页 > 案例大全 > 论文案例实践-基于alps软件的强关联电子计算研究

论文案例实践-基于alps软件的强关联电子计算研究

2021-03-26 11:18:52

  凝聚态物理学是当前物理学中一个重要的分支,本文主要探究的是凝聚态物理中的一个主流方向——固体物理学。众所周知,固体物理是一门研究固体物质的原子排列、电荷分布等特性的综合性学科,其涉猎范围包括力学、热学、光学、电磁学等等,其研究方法也被分为理论计算和实验两部分,而在本文中,我们主要关注的是理论计算部分,尤其是能带结构的计算。然而,有关固体问题的计算十分复杂,因为通常的固体物质中含有的原子数目和电子数目十分巨大,其数量级达到了1023,因此,我们引入了一系列的近似条件,不断降低计算量,进而提高计算精度,最终发展出一套成熟的能带计算理论——密度泛函理论。对于绝大多数体系而言,该理论能够在保证计算量不大的同时,得到一个相对精确的计算结果,但是对于过渡金属或稀土元素氧化物等d或f轨道上被电子填充的物质而言,由于电子间具有很强的关联的作用,密度泛函理论的计算结果与实验严重不符,因此,我们需要引出一种新的计算方法——LDA+DMFT方法,该方法在密度泛函理论的基础之上,保证了计算量的大小,同时又引入了电子关联项,确保了计算结果的准确性。最后,本文介绍了一款对于计算强关联体系十分方便的软件——alps的使用方法及其参数文件的设置方法,从而使得LDA+DMFT方法能够彻底地应用到实践。

  1凝聚态物理导读

  凝聚态物理学是一门研究凝聚状物质(固体或液体)的各种物理性质及其应用的综合性学科,包含实验和理论模拟两部分。众所周知,凝聚状物质的使用在人类历史上占据非常重要的地位。从原始人类使用简单的木制或石质工具开始,人类就已经算是开始了对凝聚状物质的研究。历史学中讲到,工具的进步标志着人类文明的进步,换句话讲,人类社会的不断发展的过程,其实也就是人类不断加深对凝聚状物质的理解的过程。然而,虽然凝聚态物质的研究最早可追溯到史前时代,但是很长的时间里,人类对这些物质的理解一直处于宏观视角,进步缓慢。

  20世纪是至今为止科学思想、技术进步最大的一个世纪,短短100年间,人类获得了远远超越过去数千年的知识储备,同样地,凝聚态物理这门学科,也在这个时期,抓住了历史的机遇,成为了一门主流学科。

  19世纪末,热力学、流体力学和弹性力学依然是人们研究凝聚状物质的主要方向,但是原子和分子等模型的建立,已经让当时的人开始思考如何从更微观的角度来描述这些物质,从而能得到一些更加一般的结论,因此,人们需要一些新思想来打破这个僵局。幸运的是,不久之后,量子力学的问世,为研究这个领域的科学家们提供了一种非常有力的工具。与此同时,一些新的实验方法,例如散射和光谱,也被引入到凝聚态物理的研究中来,这极大地增强了人们对于凝聚态物质在微观水平上的理解。然而,对于由大量原子或分子有机结合成的凝聚态物质而言,其物理化学属性与单个原子或分子间又存在着较大大差异。以价电子数为1的Na原子为例,根据量子力学,Na原子中存在着分立能级,每个能级上分布有一定数目的核外电子,其出现概率与该能级波函数的模的平方成正比。但是当两个Na原子距离很小的时候,成键和反键轨道就会建立[1],成键轨道的能量要低于反键轨道。又因为每个Na原子的3s轨道上都只存在一个电子,因此,根据洪德规则和泡利不相容原理,两个电子会成对出现在成键轨道上,并且自选方向相反。据此,我们可知,当原子结合在一起时,其能级会发生很大的变化,不再和单个原子一样是固定的,而是原子间距离的函数。所以,对于凝聚态物质而言,一般的量子力学方法显然是适用的。

  由于固体具有一定的晶格结构,也具有更加稳定的物理性质,因此相比较于液体而言,我们更倾向于研究固体来得到一些一般性结论,或者说建立一个广泛适用的理论模型。固体的量子特性的研究起源于20世纪20年代,直到今天仍然在不断发展,我们将其独立成一个学科,称为“固体物理学”(Solid State Physics)。这个学科的主要成就包括电子能带结构理论,超导理论和量子霍尔效应。其中,电子能带结构理论可以用来解释金属、半导体和绝缘体的导电性的差异。但是,至今为止,仍然有一个很重要的基本问题尚未完全解决,就是如何应对多电子体系中电子之间的强库珀作用。

  29世纪后半年叶,X射线、中子和光子散射图谱被引用到物理学中来观察一块固体物质的表面的微观形貌,最终能得到一副宏观的图像,使人们能够对其表面加以分析。这使得凝聚态物理往前迈入了一大步,因为从此之后,人们不仅可以在脑海里面想象某种固体的表面形貌,更是可以实实在在地进行观测,从而得到更加准确的判断。而既然人们想要知道固体的微观结构,那么我们必然需要一种模型来归化各种不同的固体的结构,在这个基础之上,布拉菲父子做出了卓越的贡献,最终他们也荣获诺贝尔物理奖,关于这一点我们会在下一章节中详细阐述。

  最后,走进21世纪,在这个科学与日俱新的时代,凝聚态物理迎来了它的第二春。石墨烯、碳纳米管等材料的发现,让人们不再满足于普通的三维固体材料,而是将眼光放在了二维、一维甚至零维材料上。除此之外,液晶,这种处于液体和固体之间的物质,也涌现出新的一些性质。总而言之,在21世纪的今后几十年间,我们相信,凝聚态物理学依然会是一门充满可能,充满希望的学科。

  2固体物理学

  2.1.固体物理介绍

  科学地讲,整个物理学只需要处理两类问题,一类是单体问题,一类是多体问题,其余的一切,例如二体问题或者三体问题,都可以通过简化模型,合并到这两个大类中来。之所以要进行这样的分类,是因为随着研究的不断深入和经验的不断积累,我们意识到对于数量高达1023这样的系统来讲,它的研究方法跟研究单体问题一样简单。由此可知,对于多体问题,我们拥有一套独立的处理手段。

  多体问题中最常见的就是固体问题,固体物理学也是凝聚态物理的最大分支,是一门从微观的角度来探究固体的宏观性质的学科。固体材料通常都是由大量原子构成,其数量级达到1019~1023,因此想要从微观上来解释电导、热导等宏观性质,我们需要运用大量理论知识,包括热力学与统计物理、理论力学、非相对论量子力学和电动力学等。除此之外,晶体结构学、固体比热理论、金属导电的自由电子论和铁磁理论,也为固体物理学科的发展提供了坚实的基础。在现代固体物理中,为了描述固体的性质,我们通常最关注的就是电子分布问题,这直接决定了固体的导电性能,对半导体工业具有十分重大的意义。根据量子力学,要描述电子的状态,就必须要建立相应的波函数,但是,任何固体材料中,电子的数目都是极其巨大的,因此,要建立具有几乎无穷个变量的波函数,并使用薛定谔方程进行求解,这是不现实的,所以,必须引入一些近似条件,才能够得到这些电子的运动规律:

  (1)绝热近似(adiabatic approximation):

  由于晶格上的原子中,原子核的质量远远大于核外电子的质量,因此价电子对外部环境的响应速度远远大于原子核或离子实,所以,电子能迅速调整以适应原子核的状态。据此,我们近似认为,某一时刻价电子的分布由该时刻离子在晶格上的位置决定,即电荷分布是离子位移的函数。

  (2)平均场近似(mean-field approximation):

  平均场近似又称单电子近似(one-electron approximation),当我们不考虑电子之间的关联作用时,我们仅仅只需要计算一个电子在有效场U()下运动时的波函数即可。该有效场由所有离子与除该电子之外的其余所有电子共同提供,因此,薛定谔方程简化为:

  (2-1)

  (3)周期场近似(periodic field approximation)

  根据晶体结构学,固体物质的晶格结构具有周期性,因此,我们可以合理假设由离子和电子组成的有效势场U()也具有周期性。

  U()=U(+)=n11+n22+n33(2-2)

  晶体学中,通常将三维固体的晶格结构归化成7种晶系、14种布拉菲格子、32种点群和230种空间群。结合物理学理论和晶体学知识,可以引出固体物理中最重要的一个概念——能带(energy band)。通过作出固体材料的能带图,我们可以精确地描述固体材料的电荷输运性质,从而判断固体材料的种类,即是属于导体,半导体还是绝缘体,这对现代半导体工业来讲意义重大,因为我们可以不必像过去一样,花大量时间来制备样品,然后通过测电导率来判断材料类型,而是可以直接通过能带理论来预言该材料的导电性,极大地提高了生产效率。

  2.2能带理论

  2.1.1能带模型的建立

  对于由大量原子或分子组成的固体而言,我们本能地会去想要计算每一个原子或分子的能级,这也确实能够给出一个关于电荷输运性质的正确的结果,但是基于我们现有的计算模型和计算能力,却是不现实的,尤其是针对一些晶格对称性较低的晶体。因此,我们急需一种能够在大量原子结合,彼此影响的情况下,仍能够准确描述电子运动状态的模型,而能带理论恰好提供了分析此类问题的基础。上世纪50年代,计算机的诞生,到如今的大型计算机群的建立,使得能带理论从定性的研究普遍规律发展到了能针对具体材料的极其复杂的能带结构进行精确计算,极大地促进了半导体工业的发展。

  然而能带理论也只是一个近似理论,因为固体中含有大量电子,且电子与电子、电子与离子之间都存在相互作用,所以我们无法针对这种多电子系统的波函数进行严格求解。因此,能带理论实际上是建立在单电子近似的基础之上,将每个单独的电子都看成是在等效势场中进行运动,从而降低所需计算量,换句话说,能带理论将固体的价电子近似为不再受束缚于个别原子,而是在整个固体材料中运动,这就是电子共有化(electronic sharing)。

  前面我们已经提到了,当两个Na原子相结合时,每个Na原子的价电子轨道能级会分裂成一个成键轨道能级和一个反键轨道能级(见第一章)。现在,我们将其推广到,当一个Na原子周围有N个Na原子时,Na原子的轨道能级将会分裂成N个非简并能级。由于每个Na原子3s轨道上又只拥有一个电子,又因为根据泡利不相容原理,每个非简并能级能够被两个自旋方向相反的电子所占据,因此,每个Na原子只有个非简并能级被电子所占据。当N的数目非常大时,每个非简并能级之间的能量差值非常小,从而可以在最低能级和最高能级之间形成一个准连续带(a quasi-continuum band of states),也就是所谓的能带,如图(2-1)[1]。在金属Na中,由于只有个非简并能级被电子填充,因此能带处于半满状态,又根据洪德规则,这些电子所占据的必然是从低能量往高能量不断填充。这个结论对金属的导电性具有十分重大的意义。

  图2-1

  图(a)表示当两个Na原子相互结合时,每个Na原子的3s将分裂为一个成键轨道和一个反键轨道,成键轨道被两个电子所占据,黑色箭头表示电子自旋方向相反

  图(b)表示当两个Na原子距离为a时,成键轨道和反键轨道能量差值最大

  图(c)表示当很多Na原子相互结合,且距离为a时,最低能级和最高能级之间会分裂出很多条非简并能级,并且能量较低的一半能级被电子所占据

  图(d)表示当Na原子数目逐渐增多时,最低能级和最高能级之间出现了一条准连续带,即能带,且下半部分被电子所占据

  依然以金属Na为例,当在其两端加上一定的直流电压时,金属中的电子会受到一个与电压方向相反的电场力,从而做定向运动,然而,在运动过程中,由于外部电场对电子做工,导致电子的能量不断增加,因此,他们必然要进入一个能量更高的状态。换句话讲,电子从外部获得能量时,会跃迁到未被电子占据的能量更高的状态。

  接下来,让我们来了解一下Si的能带结构,看看与金属有何不同。

  Si的最外层轨道包含2个3s电子和2个3p电子。若是根据金属能带的讨论,由于3p轨道上电子未把能带填满,因此Si应当表现出金属导体的性质,但实际上当Si原子之间相互结合时,会发生轨道杂化(orbital hybridize)现象,3s轨道和3p轨道会结合成2个能量较低的等效sp3轨道,即当电子数目很大时,会形成一条能量较低的能带,称为价带(valence band),电子会完全填充该能带。因此,当加入外部电场时,由于价带中的较高能级没有未被电子填充的空状态,导致较低能级中的电子无法跃迁到高能级中。而另一条能量较高的能带称为导带(conduction band),未被电子填充。导带和价带之间存在一个较大的能量间隙,被称为带隙(band gap)。除非获得极大的能量,否则电子无法越过带隙进入到导带中,因此Si表现出绝缘体的性质,如图(2-2)[1]。

  图2-2

  综上所述,固体材料分为金属(metal)和非金属(nonmetal),而非金属中又进一步划分为半导体(semiconductor)和绝缘体(insulator)。金属的能带结构是能带中只有能量较低的一半能级被电子完全填充,是半满状态,无带隙;半导体则是价带被电子完全填充,而导带则为全空状态,存在带隙,但是带隙较小,一半低于3eV;绝缘体与半导体类似,价带为全满状态,导带为全空状态,也存在带隙,但带隙较大,高于3eV。因此,根据能带理论,我们能够很容易地预测固体的导电性。

  2.2.2布洛赫定理

  根据单电子近似和周期场近似,能带理论的核心就是求解晶体周期性势场中的单电子薛定谔方程,得到电子的本征波函数和能量本征值。因此,首先应该考虑在该条件下,电子波函数的特点。据此,首先我们应该了解布洛赫定理(Bloch’s theorem)。

  所谓布洛赫定理,指的是当晶格势场具有周期性时,薛定谔方程的解应该为下列形式:

  )(2-3)

  其中,表示一矢量。(2-3)式表明当平移晶格矢量时,电子的波函数仅仅只是增加了一个相位因子。根据布洛赫定理,波函数可以表示为

  (2-4)

  其中,u()也具有与晶格一致的周期性

  (2-5)

  (2-5)式表示的为布洛赫函数(Bloch’s wave function),其实质为平面波波函数与周期性函数的乘积。

  2.2.3布里渊区

  由于能带的计算相比于实空间而言,更多地依赖于k空间,因此,我们需要一个在k空间中等效于实空间中的原胞的概念,被称为倒格子原胞。但是,倒格子原胞由于对称性不高,在实际计算中非常复杂,因此,我们提出一个等效于倒格子原胞,且具有高对称性的模型,即布里渊区(Brillouin zone)。此外,任意布里渊区的体积相等,且等于倒格子中原胞的体积。

  布里渊区的绘制方法是:在k空间中,以任意格点为原点,从原点出发,与周围所有最近邻格点相连,再作连线的垂直平分线,垂直平分线相交,所围成的区域,就被称为第一布里渊区,也被称为简约布里渊区;原点与所有次近邻格点的连线的垂直平分线所围成的区域,被称为第二布里渊区,以此类推出第三布里渊区、第四布里渊区……如图(2-3)

  图2-3

  以一维情况为例,由于倒格子基矢的基本性质为

  (2-5)

  所以,b的取值为

  (2-6)

  由此可得,第一布里渊区中,的取值为.

  2.2.4近自由电子模型

  在金属中,单个电子与周围的离子和其它电子之间的相互作用非常微弱,因此,我们可以忽略掉该效应,将金属中的电子假设为自由电子(金属因此被称为自由电子气)。因此,我们可以使用量子力学中对自由电子的讨论方式,设势场U()=0,然后确定一个特定的边界条件,从而求解出电子波函数。

  由于真实世界中的固体都具有边界,因此,在边界处,波函数会消失,从而形成一个驻波。但是,由于组成固体的原子数目非常巨大,高达1023,因此,边界处的原子数目相比较于固体内部的原子而言,可以忽略不计。据此,我们引入周期性边界条件(periodic boundary conditions),也被称为波恩-卡曼条件,其具体形式为

  (2-7)

  N1、N2、N3分别为沿、、方向的原胞数,因此总原胞数为N=N1*N2*N3。

  当晶体为立方晶格时,也可以表示为

  (2-8)

  L为晶体的边长,因此晶体的体积为V=L3.

  而在周期性边界条件中,某一方向上的k的取值范围为

  (h为任意整数)(2-9)

  引入边界条件之后,求解单电子定态薛定谔方程,可以得到本征波函数和能量本征值:

  ,(2-10)

  其中,根据周期性边界条件,波矢为

  (2-11)

  其中,、、为任意整数。

  然而这只是一种理想情况,在实际金属中,虽然U()很小,但是还是不能被完全忽略,而为了解决这个矛盾,我们采用了量子力学中的微扰理论,将势场对电子的作用作为一个微扰项,引入到薛定谔方程中,从而求解出更加准确的本征波函数和能量本征值,这便是近自由电子模型。由此,用平均场替代U(),我们写出波函数和本征值的零级近似为

  ,(2-12)

  能量本征值的一级近似为

  (2-13)

  其中,微扰项

  (2-14)

  据此,能量本征值的一级近似可表示为

  (2-15)

  电子波函数的一级近似为

  (2-16)

  其中,为倒格子格式,表示为

  (2-17)

  而表示周期性势场()的第n个傅里叶系数。

  能量本征值的二级近似为

  (2-18)

  综上所述,在近自由电子模型中,电子在波矢的波函数为

  (2-19)

  由(2-19式)可知,当为的整数倍时,括号内的函数大小不变,因此,该波函数符合布洛赫函数的形式,即为近自由电子模型的电子本征函数。当电子的k值接近布里渊区边界时,晶格作用会很强,因此,近自由电子模型不再适用。

  2.2.5紧束缚近似模型

  在前面我们提到,当两个原子结合时,原子的最外层能级会形成一个成键能级和反键能级,或者是原子的最外层轨道之间会发生杂化,总之,能级会分裂,单个原子的能级数目会上升,且分裂出的能级间的能量差值很小。因此,如果要做出能带图,就必须要考虑晶体内所有原子对某一单个原子作用,也就表明晶格势场的作用是较强的。据此,为了作出多原子体系的连续能带图,必须引入新的模型,即紧束缚近似模型(tight-model model),也被称为原子轨道线性组合法(linear combination of atomic obital,LCAO)。其实质为当电子在一个原子格点上时,主要受到该原子形成的势场的作用,而其它原子势场的作用视为微扰,而且一般而言,只考虑最近邻原子和次近邻原子势场的围绕作用。

  紧束缚近似模型适用的是共价键晶体,例如金属氧化物,而非和近自由电子模型一样适用于一般金属。但是对于过渡金属原子的d轨道电子,是一种局域化电子,也可以使用紧束缚近似模型进行分析。下列,我们将讨论其中一种最简单的情况。

  首先,我们需要建立一个组成固体的原子的哈密顿算符(仅考虑单原子固体),可以被表示为

  (2-20)

  其中,U()表示在各格点处的原子形成的势场。

  我们知道,对于理想晶体而言,组成晶体的每个原子都位于空间点阵的格点上,每个原子都能够对单一电子贡献一定的势能,因此,哈密顿算符可以进一步表示为

  (2-21)

  而某一格点上的原子周围所束缚的电子的波函数可以表示为该原子能级En上束缚的电子的波函数,则波动方程为

  (2-22)

  将方程两边同时乘以,再对整个体系进行积分,可得:

  (2-23)

  其中,-β表示某一格点上的原子的能级,因周围原子的作用而产生的小转变,即因周围原子的作用而产生的微扰项。因为原子的能级为负,所以为了保证β为正,将负号提出。

  由于晶体中一共有N个原子,因此对原子能级En对应的的哈密顿量求解,即计算(2-21)式,可以求解出N个形式相似的电子波函数,其具有相同的能量值,换句话讲,该能级是N重简并。因此,该能级上的电子的波函数可以表示为对原子轨道哈密顿量求解之后得到了N个波函数的线性组合,故此,该方法也被称为原子轨道线性组合法。同时,根据布洛赫定理,该波函数还应该具有布洛赫波的形式,因此,其波函数的具体形式为

  (2-24)

  其中,是归一化系数,没有其它作用。

  利用该波函数求解体系的能带结构E(),可得:

  (2-25)

  该式考虑了有限原子组成的晶体中的所有原子,但是由于周期性边界条件的限制,导致晶体中无论表面还是内部的原子都是等效的,因此,为了简便运算,令R’=0,得:

  (2-26)

  根据(2-22)式,又可将(2-26)改写为

  (2-27)

  同样根据式(2-22),可得:

  (2-28)

  由于和分别表示两个不同格点的波函数,根据正交归一化条件,(2-28)的前一项为0。令后一项为,则

  (2-29)

  因此,可将(2-27)式改写成

  (2-30)

  据此,得到紧束缚近似模型中的能带计算公式。

  2.2.6费米能级与能态密度

  原子的能级是分立的,每个能级中含有一定数目的量子态,其数量与其简并度相等。但是当N个原子结合在一起形成固体形成能带时,会消除能级的简并度,也就表明每一个能级上仅能存在一个电子(考虑电子自选),因此,能带中所蕴含的能级数目是巨大的,能带是一个准连续带。由此可见,讨论单个能级已经没有意义,针对这种情况,我们引入能态密度(density of states,DOS)的概念。

  考虑在区间内的能态数目,设其为,则其中能态密度可表示为

  (2-31)

  假设k空间中,状态的分布是均匀的,其分布密度为

  (2-32)

  其中,V为晶体的体积.

  能带中,能级间的体积可用体积元在k空间中的积分来表示,因此,能态数目可表示为

  (2-33)

  其中,为两能级之间的垂直距离,则

  (2-34)

  则能态数目可以进一步表示为

  (2-35)

  从而得到能态密度N(E)的一般式为

  (2-36)

  若考虑电子自选,则能态密度加倍

  (2-37)

  进一步讨论能带的性质,我们想要知道如何计算电子所能填充的最高能级。依然以金属Na为例,由于Na原子只具有一个价电子,因此由N个Na原子组成的体积为V的金属块材中,价电子的数量为N,因此价电子密度为

  (2-38)

  根据洪德规则和泡利不相容原理,电子从能带中最低能级向高能级填充,且填充的状态数为,进一步假设将所有价电子状态都填充进一个半径为的球中,则可以得到电子状态数为

  (2-39)(2-40)

  其中,为费米波矢(Fermi wave vector),也可用表示,为电子在能带中所能填充的最高能级中的电子本征态对应的波矢。假设电子为自由电子,则根据(2-10),可得

  (2-41)

  该能级为电子在能带中所能填充的最高能级,被称为费米能级(Fermi energy)

  2.3密度泛函理论

  2.3.1密度泛函理论介绍

  密度泛函理论(Density functional theory,DFT)又被称为非均匀电子气理论(theory of the inhomogeneous electron gas),是当前针对固体材料的电子结构计算的一个标准工具,目前,密度泛函理论已经发展出一套完备的体系来探索晶体的物理化学性质。但是,市面上大量质量参差不齐的介绍密度泛函理论的书籍与文章,使得初学者在学习过程中仍有较大的困难。因此,在这一节中,我们将要系统地梳理密度泛函理论的知识,使得读者能够掌握其核心思想,也为我们在下一章中将要开展的关于强关联电子问题的讨论与计算打下基础。

  密度泛函这个词,其实就是针对以电子云密度为基本变量的函数进行求解,因此,密度泛函理论其实质就是找到在一系列实际问题中描述原子或分子行为的一系列薛定谔方程的解。因此,我们可以知道,要想得到晶体的电学性质,能带理论是一切问题的出发点,我们的基本思想永远是构建电子的哈密顿量并进行求解,从而得到电子波函数和能量本征值。但是,在复杂的实际问题中,构建电子的哈密顿量却并不是一件容易的事情,需要进行大量复杂的数学物理处理,而密度泛函理论就是在众多数学工具中,最为成功的一个。这个方法,从最开始只有少部分处于量子力学前沿的物理学家和化学家的研究开始,发展到现在已经成为了物理、化学、材料科学、化学工程等学科的必备工具之一。根据对科学检索指数的研究,我们发现,在1986年,有关密度泛函理论的词条少于50条,但是1996年至2006年间却从1100条增长至5600条[2]。

  2.3.2第一性原理

  在前面的论述中我们已经知道,要想得到晶体的电子结构,就必须针对多粒子薛定谔方程求解,得到本征函数和本征值。而在这个求解过程中,为了得到准确的结果,我们将最大限度地进行非经验处理,即不涉及到任何的经验参数,只需要输入原子的核电荷数和一些环境参量之后便开始计算,从而得到一个与真实情况相符的结果。这种方法,被称之为第一性原理模拟(first principle simulation)。

  第一性原理是目前在计算物理领域非常流行的一个术语,基于该思想,我们能够对算法和模拟结果进行优化,从而使得理论模拟值能够更加符合真实情况。但是这种思维方式,不光是在计算物理学界建树颇深,在提升个人知识水平上也是有一定帮助的。例如,当看见一棵果树上结下一个果实之后,人们立刻就能明白这棵树是一棵果树,但是,对于擅长从第一性原理分析的人来讲,他会抛却掉这种“经验”,将整个系统进行解刨,从果实本身、树枝、树干不断延伸到根部,最终得到这棵树会结这种果实的原因。对于这种思维方式,我们不再是递进的,填鸭式地学习,而是从原始情况出发,构建出我们自己的逻辑思维体系,从而得到一个最真实的结果,这也是第一性原理名称的由来。

  2.3.3平面波展开

  从实用角度来看,平面波展开的方法并不是计算晶体能带结构的好办法,但是它的概念简单,便于明晰电子结构计算过程[3]。

  我们已经知道,计算电子结构,其实质就是求解晶体中电子的薛定谔方程。但由于晶体中的电子数目巨大,且电子与电子,电子与离子实均存在关联作用,因此我们无法得出一个精确解,只能通过一些近似条件,得到一些近似的解。而对于一个未知的波函数的求解,通常采用的方法是使用已知的波函数,根据波函数的正交完备性,将其作为基函数(basic wave function),构建一个表象空间,将目标波函数进行展开,从而求解出薛定谔方程。其中,最简单的基函数就是平面波波函数。这种方法的的优点是计算量小,因此,当参与计算的原子数目较大时,平面波展开是一种优秀的计算方法。但是,其计算精度却不是那么令人满意,当涉及到局域性的精确计算时,如紧束缚模型,一般选用分子轨道法,也就是原子轨道线性组合法(LCAO)。

  2.3.4赝势

  “赝”在现代汉语中意思是伪造的,因此,赝势(pseudo potential)表示的不是真正的势。目前在能带结构的计算中,赝势的概念已被广泛采用。

  在晶格原胞中心,接近原子核时,由于库伦相互作用的影响,使得势场偏离平均值很远,因此波函数会发生显著变化,从而导致使用平面波展开进行计算时,能量收敛很慢。但是,固体的形成的根本原因是原子间成键,可能是化学键,也可能是范德瓦尔斯键、氢键等,但其本质就是电子云,所以其与波函数的特性密不可分。因此,在接近原子核的部分,即便是波函数发生异常,但是其仍然对于电子能量给予贡献。据此,为了简化计算,同时又能给出一个关于外层电子能量的较为真实的结果,我们将内层,不参与导电的电子与原子核并入到一起,构成一个离子实,将该离子实的原子核和电子的势场也看作一个统一的,虚拟的势,即赝势。最初是由Hellman在1934年提出,后在1959年,经Phillips和Kleinman进行完善[3]。由赝势求出的波函数被称为赝波函数,在离子实和价电子之间的区域,真实的势场和赝势会给出同样的波函数。

  2.3.5瓦尼尔表象

  据此,我们仅仅讨论了理想晶体的能带理论,但是在实际情况中,还会存在外部场强或杂志和缺陷的影响。因此,在求解单电子薛定谔方程时,应该要加入外部势场U的作用,从而将薛定谔方程转变为

  (2-41)

  求解上类方程是,可根据现有的,可根据理想晶体的H所求解出的具有正交完备性的波函数为基函数展开,表示实际情况下的波函数Ψ。如果是利用布洛赫函数为基函数的话,则可以构建出布洛赫表象空间,但是,我们还可以通过布洛赫函数,定义出另外一套描述局域电子状态的正交完备的函数组——瓦尼尔函数,从而形成瓦尼尔表象空间。

  根据式(2-4)和(2-5),布洛赫函数是k空间中的,以一个格点为周期的周期函数

  (2-42)

  根据倒点阵周期函数按正格矢展开的点阵傅里叶级数

  (2-43)

  可将按正格矢展开为傅里叶级数

  (2-44)

  利用

  (2-45)

  求得其逆变换为

  (2-46)

  其中,与只与实空间中的位置有关,与无关。根据布洛赫定理可得

  (2-47)

  因此,可将(2-46)改写为

  (2-48)

  式(2-48)所定义的即为瓦尼尔函数。不同的能带和不同的格点有不同的瓦尼尔函数。

  瓦尼尔表象在计算完整晶体能带时并不是十分方便,但是在讨论局域杂质的电子能谱时却十分有效[4],尤其当局域杂质势能U()

  较强时,量子力学的微扰理论失效,利用瓦尼尔函数展开(2-41)中的波函数Ψ,可以得到一个有效哈密顿量,其求解非常方便。

  2.3.6 Hartree–Fock方程

  对于多体体系而言,要想求解出电子结构是非常复杂的,必须知道离子实对电子的势能,但是即便是由两个具有相互作用的电子构成的系统,其薛定谔方程为

  (2-49)

  由于相互作用项的存在,无法使用分离变量法求解该方程,因此需要进行近似处理。而之前提过的单电子近似就是处理此类问题的,不合理程度最低,结果最接近真实值的近似模型。接下来,我们将要讨论单电子近似下,离子实对电子的势能的具体形式。

  单电子近似首要考虑的是如何利用N个只含有关于一个电子的变量的单电子波函数构建以电子变量为宗量的N电子系统总的波函数Ψ()。Hartree最早采用的是连乘积形式:

  (2-50)

  根据量子力学变分原理,即微观体系的能量本征值可通过在一定的边界条件下求解薛定谔方程,并要求所得的波函数满足正交归一化条件而实现,对N个单电子波函数的选取,应当使N电子系统的能量本征值达到极小,如此一来才能形成一个稳定的多体系统。而N电子系统的能量本征值可表示为

  (2-51)

  对变分,并引用正交归一化条件

  (2-52)

  使得

  (2-53)

  其中是拉格朗日乘子(lagrangian multiplier)。据此得到单电子波函数须满足的有效薛定谔方程为

  (2-54)

  式(2-54)就是Hartree方程。式中左边第一项为动能,第二项为原子核对j电子形成的势能,第三项为其他N-1个电子对j电子的库伦作用能。

  但是,Hatree利用连乘积方式构建出的N电子体系总的波函数不符合泡利不相容原理,因此,在1929年,Slater考虑了交换电子的反对称性(电子是费米子),要求利用N各正交归一化的单电子波函数以行列式的方式来构建N电子体系总体波函数[3]:

  (2-55)

  当固体问题中可近似地不考虑自选—轨道耦合问题时,可将波函数利用分离变量法,写成空间函数与自选函数的乘积,即

  (2-56)

  将自选函数积分后,其便不再出现,如此一来,N电子的能量本征值为

  (2-57)

  其中,是固体离子对电子形成的势能,通过变分法,即

  (2-58)

  可得

  (2-59)

  通过幺正变换

  (2-60)

  使得波函数集满足(2-55)等号左边的形式,同时约化右边为,则可以得到

  (2-61)

  上式即为Hartree-Fock方程,相比于式(2-54)的Hartree方程,左边出现了没有经典对应的第四项:交换作用项。其中,求和部分只涉及与j电子态自选平行的j’态,是电子服从费米-狄拉克分布的体现。

  2.3.6密度泛函理论

  当利用Hartree-Fock方程计算能带时,有两个重要的缺陷[4],第一是当实际系统中,一个电子的状态发生改变时,难以保证其它(N-1)个电子的状态不发生变化。此外,Hartree-Fock方程仅仅只考虑了电子间的交换作用,而完全忽略了自旋反平行电子之间的相关能。因此,Hartree-Fock方程不能作为单电子近似(即能带论)的理论基础,而密度泛函理论则很好地弥补了这两个缺陷。

  密度泛函理论是建立Hohenberg-Kohn定理基础之上的,该定理最开始Hohenberg和Kohn两人提出,但在上世纪60年代中期,由Kohn和Sham进一步将其完善[2]。

  最初的Hohenberg-Kohn定理指的是薛定谔方程所计算出的基态能量是电子密度的唯一的泛函。这个定理将基态电子波函数和基态电子密度联系起来,而为了理解这个定理,我们需要理解什么是“泛函”。

  “函数”一词表示用一个或者一组数值变量来表示一个单一数值变化的概念,而与之相似,“泛函”一词则表示用一个或一组函数为变量,来表示另一组变量的变化。换句话讲,如果函数体现了坐标空间中的(x,y,z)变化趋势,那么泛函则表示了在函数空间中的函数f(x)的变化趋。那么现在,我们就可以说明,Hohenberg-Kohn定理实质就是将基态能E表示为,其中表示电子密度或电子云密度。而这也说明电子云密度可以决定体系的所有微观性质,包括基态能量和基态波函数。同时,这意味着,当我们想要求解由N个电子组成的系统时,我们只需要建立以x,y,z和电子云密度为变量的函数关系即可,而不需要和Hartree-Fock方程一样,建立N个表示不同电子的波函数。如此一来,求解薛定谔方程的过程将会被极大地简化,同时结果也更加精确。例如,针对一个由100个原子组成个纳米团簇,我们可以将其有一个23000维的问题简化成一个三维问题,也就是在一个三维空间中找到其电子云的分布即可[2]。然而,Hohenberg-Kohn定理却并没有告诉我们基态能量和电子云密度之间具体的关系式,因此,改进后的hohenberg-Kohn定理为在对应于薛定谔方程求解出波函数的电子云密度中,能够在泛函中使能量达到极小值的电子云密度是真实的电子云密度。如果说泛函空间的形式已被知晓,那么可通过一些对基函数的近似处理来变换电子云密度,使其成为一种相对电子云密度,使得能量极小化。

  众所周知,想要知道多体中的基态能量是非常困难的,但是,根据Hartree,Kohn和Sham等人所建立的理论,我们可以通过极小化能量泛函中的能量来实现,具体可以通过找到一系列的单粒子波函数的自洽计算的解来实现。而为了便于计算,我们假设系统中的电子密度是均匀的,也就表明电子密度为常数。虽然这种假设在真实材料中是不可能的,因为物质的组成是依赖于化学键,而化学键又是因为系统中电子密度不一致才形成的,但是这种思想对于我们求解泛函中的能量极小值,却是有很大帮助的。因为基于该假设,我们可以确定系统中某一格点周围的电子云密度,从而确定该格点的能量函数,进一步求解出薛定谔方程,这种方法被称为局域密度近似(local density approximation,LDA)。但是,这种方法得到的薛定谔方程的解是不精确的,因此,发展更加真实地描述自然的近似方法,一直是量子化学领域的热门课题。现在,在LDA之后,一种基于局域梯度电子云密度的方法应用最为广泛,这种方法被称为广义梯度近似(generalized gradient approximation,GGA)。因为该方法相比较于LDA而言,考虑了更多实在的物理信息,所以其运算结果更加精确。

  然而,基于单电子近似的密度泛函理论在很多情况下也是不适用的,尤其是对于强关联电子体系,主要是因为这种材料的d轨道和f轨道的电子之间会发生很强的,不可忽略的关联作用,因此无法做近似化处理,这也是我们在下一章中要进行论述的核心内容。

  2.3.7局域密度近似

  根据上一节,由Hohenberg和Kohn建立的密度泛函理论可知,电子密度的泛函是基态能量。该泛函可以通过在电子数为N,电子密度为处极小化与所有的多体波函数相关的能量期望值来确定:

  (2-62)

  其中,哈密顿量为

  (2-63)

  其中,和都是场算符,可以创造和湮灭在处,自旋为的电子,

  ,(2-64)

  代表在处的带着电荷数的所有离子和电子之间的势场以及所有电子与电子之间的势场。而式(2-62)还可以利用Hartree能和离子势能来表示:

  (2-65)

  (2-66)(2-67)

  其中,代表动能,代表未知的交换关联项,该项包含处Hartree项之外的电子相互作用能。因此,多体问题中的所有困难都可以转化为求解。尽管动能无法利用电子密度得到一个精确的表达,但是可以通过一个技巧来完成该项工作——不是极小化与密度相关的,而是极小化与相关的一系列单粒子波函数:

  (2-68)

  为了保证的正交归一性,拉格朗日参数被引入,通过下式(2-69)

  得到下列Kohn-Sham方程:

  (2-70)

  上述方程均构成了同一个单粒子薛定谔方程,然后,该方程通过单粒子波函数拟设的方法计算动能。如果是式(2-70)和(2-68)的基态的自洽解,那么这种基态密度所拥有的单粒子拟设的动能是符合下式

  (2-71)

  可是,要注意式(2-70)的单粒子势,例如下列式子

  (2-72)

  是在极小化方法中唯一的附加势能。因此,波函数和拉格朗日参数在该极小值点上不具有物理上的平均值。同时,这些方程也可以进行DFT/LDA计算,如图2-4。

  因为多体问题可以转化为解决未知的问题,因此利用一些近似条件是必须的。针对这样的问题,局域电子密度近似LDA被建立:

  (2-73)

  其中,通常是由胶体模型的微扰解或定量解计算而来,而这种模型定义。

  第一性原理信息:

  原子数;晶格结构(晶格类型和原子位置)

  选定初始电子密度

  反

  复

  进

  行

  自

  洽

  计

  算

  利用(2-72)计算有效势

  求解Kohn-Sham方程(2-70)

  计算电子密度(2-68)

  计算能带结构,部分和完整地DOS,LDA下自洽的哈密顿量LDA+DMFT,全部能量…

  2.4 DFT/LDA计算流程图

  原则上讲,DFT/LDA只能计算一些稳定的性质,例如基态能或其衍生物,但是LDA的一个主要应用就是计算能带结构,因此可以被认为是系统的单粒子能级。因为实际的基态并不是一个简单的单电子波函数,因此LDA是在DFT之上的一个近似。事实上,这种近似取代哈密顿量(2-63)为

  (2-74)

  对于实际计算而言,拓展与基函数相关的场算符是必要的(i代表晶格格点,l和m代表轨道指数),例如线性松饼盒轨道基函数(linearized muffin-tin orbital,LMTO)。在这个基函数中,

  (2-75)

  这样式(2-74)的哈密顿量可变换为

  (2-76)

  其中,

  (2-77)

  (2-78)

  其中,,代表对应对角部分。

  对于稳定的性质,LDA方法基于式(2-76),也可以非常成功地求解出能带结构,但这仅仅只针对于弱相互作用材料,对于一些强关联材料,如过渡金属或稀土元素氧化物来讲,LDA方法具有较大的误差,因此需要运用接下来将要深入讨论的LDA+DMFT方法。

  3强关联电子体系

  3.1强关联作用简介

  固体物理学解释了很多材料的性质,例如简单金属和一些半导体和绝缘体,但是对于一些具有d和f轨道电子的材料而言,这些电子会占据窄能带,因此这些材料会表现出一些独特的,难以解释的特性[5]。例如过渡金属中的钒、铁等元素以及其氧化物,由于他们的外层电子在空间中的局限性,导致其会出现强的库仑排斥作用。这些电子间的相互作用影响太过显著,因此无法被当成是独立的,换句话讲,电子彼此之间存在着相互影响的关联作用,而这些电子,则被称为强关联电子(strongly correlated electrons)。

  强关联电子问题是当今凝聚态物理学界的非常重要,也非常困难的问题之一。依据单电子近似建立的能带理论其实忽略了电子间的关联作用,而实际上,在窄能带中,这些作用力又十分重要,因为当电子从一个局域原子轨道跃迁到另一个局域原子轨道时,如果该轨道已被电子占据,则必须要考虑电子间的库伦作用,这会导致能带结构发生显著的变化[4]。在强关联体系中,不仅包含了一些经典的系统,如电子化合物,也包含了量子相变、重费米系统、非费米液体等蕴涵一些奇妙的物理现象的新兴系统,因此强关联电子相关的理论研究,是当前物理理论的发展和新材料的研发与制备的重要基石之一。

  目前,人们针对强关联电子体系已经做了很多针对性的研究,我们建立了一系列的有关强关联电子系统的物理模型,其中最经典的莫过于赫伯德模型,并且将清晰的物理思想和各种有效的数学工具相结合,这些工具包括量子蒙特卡洛方法、格林函数方法等等,在于第一性原理模拟相结合,从而得到强关联电子系统的各种晶体学性质。

  蒙特卡洛方法是应用于求解统计物理学问题的最重要的一种方法。对于一个多体问题而言,像经典力学一样,对单个物体构建微分方程,从而求解出其运动规律的意义并不大,因此我们会应用统计物理学的知识,讨论一系列物体的整体性质[6]。因此,针对这个全新的领域,我们也就将要使用一些独特的数学方法,也就是蒙特卡洛方法。另一方面,从上世纪50年代开始,量子场论中的格林函数方法已经被用于研究统计物理学中的问题。到了60年代中后期,格林函数理论已经在许多方面取得了重要成就。目前,格林函数方法已经成为了研究固体物理中具有各种复杂的相互作用的多粒子系统的有力工具[7]。对于强关联电子体系的具体计算理论,主要有解析解和数值解两种,这里我们主要讨论数值解的计算。而数值解的计算,最重要也最广泛的一种理论就是动力学平均场理论。上述这些方法都可以成为计算该理论的有效工具。

  综上所述,强关联电子究竟可以应用于哪些方面是一个我们非常关注的问题,而实验表明,用单电子近似基础之上能带理论来解释过渡金属化合物的物理特性是不正确的。例如,MnO晶体中,每个原胞中有5个3d电子构成未填满的3d能带,而氧离子的2p能带是满带,且不与3d重叠,因此不会发生分子轨道杂化,这么说来,根据能带理论,MnO的3d能带将具有金属导电性。但实际上,实验测得MnO的电导率很低,在室温下只有10-15,说明MnO实际上是绝缘体。除此之外,还有些过渡金属氧化物在温度变化时,会发生金属-绝缘体转变[4],这些都是普通的能带理论所不能解释的。而这些现象产生的主要原因,就是因为基于单电子近似的能带理论忽略了窄能带中的电子强关联效应。

  总而言之,强关联电子问题是当代凝聚态物理中的一个热点问题,有关其的理论知识浩瀚深远,接下来,我们会逐步深入地了解。

  3.2二次量子化

  在正式引入强关联电子模型之前,我们要介绍一个在非相对论量子力学中非常重要的概念——二次量子化。我们已经知道,在自然界中,存在着两种全同粒子——玻色子和费米子,名称来源于他们分别符合波色-爱因斯坦分布和费米-狄拉克分布,同时,他们又各自具有对称性和反对称性。而对于由N个全同粒子组成的系统的量子态的描述,一般的在坐标表象和动量表象的波函数难以胜任,因此,通过引入粒子数表象和交换、产生和消灭算符,我们能够更好地表示全同粒子系的状态,这种方法被称为二次量子化方法(the method of second quantization)。一次量子化是用波函数对量子系统的一种描述,而二次量子化是用算符对量子系统的另一种描述。

  3.3赫伯德模型

  赫伯德模型(Hubbard model)一开始是为了讨论过渡金属中的d轨道形成的窄能带中的电子而建立的,如今已经成为了解决电子强关联问题的最简单有效的模型之一[8]。

  首先,我们需要建立赫伯德模型的哈密顿量。考虑由N个原子组成的简单晶体,假设为单电子在周期性势场中的哈密顿量,则多体系统组成的总的哈密顿量为

  (3-1)

  其中

  (3-2)

  为简单期间,我们只考虑单个未填满电子的能带,例如孤立的s带。然后,我们引入一个二次量子化的概念,即使用产生和消灭算符在对称化的Hilbert空间中处理全同粒子的方法。然后,我们表示出在布洛赫表象中的电子相互作用哈密顿量的二次量子化表达式为

  (3-3)

  其中,和代表布洛赫表象中,自旋电子的产生和消灭算符。但是,为了更好地讨论窄能带中的电子强关联问题,最好采用的是瓦尼尔表象,因此利用变换关系(2-44)和(2-48)以及

  (3-4)

  将式(3-3)利用瓦尼尔表象表示为

  (3-5)

  其中,和是瓦尼尔表象空间中自旋电子的产生和消灭算符,单粒子项中

  (3-6)

  相互作用项中

  (3-7)

  (3-7)式一般为多中心积分,但可简化为单中心、双中心积分。赫伯德指出[4],对于窄带单中心积分

  (3-8)

  最重要,其代表同一原子周围能带电子之间的库伦作用,是近距离作用,其数量级可达10eV,比双中心、三中心等积分的贡献都要大得多。赫伯德模型是一个简单的模型,因此,可以略去相互作用项中其它所有的多中心积分,只留下一个单中心积分,这也就是赫伯德模型。因此,多中心积分近似表示为

  (3-9)

  其中

  (3-10)

  (3-11)

  式(3-11)代表了在i格点上的自旋的粒子数。根据泡利不相容原理,处于同一格点的两个电子自旋方向一定相反,因此

  (3-12)

  所以,式(3-5)可近似表示为

  (3-13)

  注意,在赫伯德模型中的关联自旋相反的电子的排斥势U对金属-绝缘体转变和窄带的磁性而言具有重要的影响。根据式(3-6),H的单粒子项中国i=j项的系数为

  (3-14)

  该式代表能带电子的平均能量。

  3.4动力学平均场理论

  3.4.1背景介绍

  动力学平均场理论(dynamic mean-field theory,DMFT)是近十几年来发展起来的,对处理强关联相关问题的一种强有力的方法。动力学平均场理论是在空间维数趋于无穷时精确的理论[9],其基本思想使用杂质模型来代替一个格点,来描述电子之间的局域相互作用。该杂质所处的没接被称为Weiss分子场,是由自洽方程所决定的媒介。是由此可看出,动力学平均场理论的核心在于用一个自洽方程决定的平均场,替代周围格点的影响,将它作用在一个格点上,从而将标准的格点模型简化成可以求解的杂质模型[9]。因此,动力学平均场理论又被称为局域杂质自洽近似(local impurity self-consistent approximation)。

  3.4.2平均场理论

  平均场理论(mean-field theory)是当某一系统的空间涨落可忽略时,对系统的热力学性质的一种近似,当粒子之间的相互作用为无穷时,该近似结果会变得非常精确。在空间维度很高的相变的过程中(例如临界指数),平均场理论可以针对这一进程做出定量的描述和预言,此外,还具有数学形式简单明了等优势。总而言之,平均场理论是第一个预言相图和新的实验过程的方法[10]。

  而平均场理论之所以被称之为平均场,是因为其使用了一个平均值来代替空间中的局域变量(例如自旋),因此忽略了其波动范围。但是,在某些情况中,这种波动是不可忽略的。

  当大量自旋能够与一检测自旋(test spin)具有相互作用时,若存在检测自旋,则会产生一个有效的平均场,但是如果自旋仅与两个检测自旋相互作用时,产生的平均场会很小,这时空间中的涨落现象是相对来讲很大的,因此不可忽略。由此可见,自旋的数量对于产生有效场而言是很重要的,平均场理论适用于相互作用范围广和维度高的系统,对于低维局域系统来说,无法提供一个正确的描述。

  3.4.3强关联体系计算方法---LDA+U

  LDA+DMFT是一种全新的针对强关联电子体系下的材料,如过渡金属及其氧化物等的第一性原理计算方法。这种方法将传统的局域密度近似下的能带理论,或者说密度泛函理论与动力学平均场理论很好地结合起来,能够在保证计算精度的同时,极大地提高考虑了电子将强关联作用的系统的计算速度。近几年来,LDA+DMFT方法已经发展成为了一种针对强关联体系计算的非常有力的工具。在接下来的讨论中,我们将主要以强关联体系中最简单赫伯德模型为例。

  尽管能带理论中的局域密度近似(LDA)在很多方面都取得了巨大成功,但是仍然存在着很多不足,例如s-d轨道结合的晶体,其结合能结果较真实值而言偏大;碱金属禁带宽度的计算结果也偏大;而对于一些具有3d或4f电子的体系,其计算结果离实验值相差甚远[8],因此,对于这些电子之间具有强库伯作用或者交换作用的体系,我们通常采用一种被称为LDA+U的方法[11]。LDA+U是一种在LDA基础之上,考虑了电子库伦作用和交换作用的轨道独立性,而这种方法其实也就是第一性原理的具体应用,因为其初始条件中没有考虑任何参数。当计算过渡金属或稀土元素氧化物时,无论是针对激发态,如带隙,还是基态,如磁矩的计算,LDA+U法都在LDA基础上做了一个定性的提升。

  在LDA下的DFT是一种最成功的第一性原理计算方法。多体问题都被划分进一个无相互作用的,只包含一个由电子云分布决定的单电子交换关联势中。LAD已经被证明在大规模原子组成的固体中,是一种极为有效的计算方法。但是,LDA却没有办法精确地得到材料的基态性质,强关联体系也无法计算。这些体系通常包含过渡金属和稀土元素,如上文所讲的那般,d轨道或f轨道都被电子填充。当应用LDA这种,基于轨道无耦合势的单电子计算时,岂能带结果是一个部分填充金属型电子和巡游电子的d带[11],但这是完全错误的,因为对于过渡金属和稀土元素化合物而言,d轨道上的电子是局域化的,而非公有化,在被占据和占据的轨道之间,能级差也较大,无法形成准连续带。因此,我们尝试了很多种提升LDA精度的方法。

  最广为人知的一种对LDA修正的方法就是自身反应修正(self-interaction correction,SIC),它重申了过渡金属或稀土元素氧化物中d轨道电子的局域性,但是这种方法所计算出的电子能量与光谱学数据不一致。Hartree-Fock方法是一种描述莫特绝缘体(Mott insulator)的近似方法[12],这种方法没有考虑自身反应,因为这导致计算结果与光谱学数据的不一致。但是Hartree-Fock近似却没有考虑库伦屏蔽作用,而库伦作用的裸值通常来讲比较大,高达15-20 eV,但是在固体中,由于屏蔽作用,则要小得多,通常是8 eV及其以下[13,14],因此Hartree-Fock近似的计算结果要比实验值大上许多。

  为了处理上述的问题,我们采用了一种被称为GW近似(GWA)的方法[15],其可以被看做是考虑了库伦屏蔽作用下的Hartree-Fock近似,其应用范围很广,从一般的单质金属到过渡元素金属都可以,但是对于一些更复杂的化合物,由于其计算量太大,因此并不是一个可行的方案。此外,利用GWA进行计算时,我们需要以LDA为基础的能带和波函数计算结果作为一个响应函数,来得到相应的库伦屏蔽作用能的大小。但是,尽管这些方法已经被证明在弱相互作用系统(例如半导体)是正确的,对于强关联系统来讲,比起LDA或者上述对LDA的修正方法而言,仍然需要一个更好的初始哈密顿量。

  在这样的背景之下,一种名为LDA+U的计算方法油然而生,在强关联系统中,可以提供给我们一个更加精确的结果。

  LDA+U方法是一种基于Hubbard模型的,对LDA进行纠正的方法[16],是一种描述强关联系统中的基态能级的最简单的方法。由于其表达形式简单,计算量也相对较低,仅略大于标准的LDA计算,因此,LDA+U(如果不特别指定的话,我们还可以将赫伯德模型中的U指定到其它近似模型中,例如局域轨道密度近似LSDA和广义梯度近似GGA)很快就成为了一种在第一性原理计算中计算强关联系统的通用方法。下面我们就来详细说明LDA+U方法的理论依据。

  根据安德森模型[17],我们可将体系中的电子划分到两个亚系统中:(1)以d和f轨道为代表的的局域化电子,其电子之间具有很强的库伦作用力,因此在考虑平均场近似时,要引入赫伯德模型中的(是被占据的d轨道);(2)以s和p为代表的公有化电子,这些电子使用轨道无耦合的单电子近似进行描述。下面,让我们来考虑一个具有d电子的离子,并且其d电子数目不固定。首先,根据LDA,d电子之间库伦作用的能量值与d电子的数目之间存在一定的函数关系,我们假定其库伦作用能量值为一个近似(尤其注意不是轨道能量值,其为d轨道上电子的本征值),然后结合上文中的两个表达式,将该能量修正为:

  (3-13)

  现在,让我们从LDA总能量中,减去该修正能,同时引入一个类赫伯德模型项,其结果如下:

  (3-14)

  而轨道电子的本征能量可由(3-14)中关于被占据的轨道电子数衍生出:

  (3-15)

  该变型的本质是当电子被占据()时,将LDA轨道能量降低;当其未被电子占据()时,将LDA轨道能量增加。而针对于轨道无耦合势则可应用一种相似的形式:

  (3-16)

  其中,该变量并不包含总的电子密度,而是一个特定的i轨道的电子密度。上式的轨道耦合势给定了在较高和较低的赫伯德能带中,能量的差值等于库伦参数U,从而可以定量地再现莫特-赫伯德绝缘体的正确物理性质。为了定性建立一个更合理的计算方案,我们需要以更一般的形式建立一个轨道基集,同时也需要考虑在部分填充的d或f轨道中的直接库伦作用和交换库伦作用。

  从物理上讲,上述需要都被建立在空间中一个特定区域的识别中,而电子本征态的特征很大程度上就保留在这些区域中,因此该区域被称为原子球(atomic spheres)。对d或f的电子来讲,这些讨论没有任何问题。在这些原子球中,我们可以衍生出一组局域正交基(i表示轨道位置,n为主量子数,l为轨道角动量,m为磁量子数,为自旋量子数)。虽然这组局域正交基不是严格必要的,但是接下来,还是以此为基础,详细阐述只有nl壳层被填满的一般情况。密度矩阵可以定义为:

  (3-17)

  其中

  (3-18)

  (3-18)为格林函数矩阵在局域化表征中的一个元素,而关于,我们会在之后进行定义。根据密度矩阵,广义的LDA+U方程可以被定义为[18]:

  (3-19)

  其中,表示考虑自旋的电子密度,表示一个标准的局域自旋密度近似(local spin-density approximation,LSDA)泛函。该式表明LSDA满足轨道无极化条件,而后一项则可表示为:

  (3-20)

  其中,表示在nl轨道上电子之间被屏蔽的库伦作用。而式(3-19)的最后一项纠正了双计数(当系统无轨道极化时,(3-19)可简化为),具体形式为:

  (3-21)

  其中,U和J代表被屏蔽的库伦和交换作用参数。

  在一般的LDA势中,我们发现一个有效单电子势,可以应用到有效单电子哈密顿量中:

  (3-22)

  其中

  (3-23)

  (3-21)和(3-23)中,

  ,(3-24)

  根据(3-23),我们可以得到。通过假设在原子球中,这些相互作用依然维持着他们原子性质,我们再次遵循着LDA+U的核心思想。此外,LSDA方法自己也可以通过被严格检测的超胞LSDA方法程序来确定自身的值,也就是中的元素。

  式(3-22)中,新的哈密顿量以算符的形式包含了式(3-23)中的轨道无耦合势,这意味着LDA+U方法本质上是依赖于在哈密顿算符下局域化轨道集的选择,这也是LDA+U方法的安德森模型式基本思想的结果:整个变化的空间分散进局域亚轨道,如d或f轨道,彼此之间的相互作用伴随着赫伯德型的的哈密顿量,以及库伦相互作用的LDA被认为是充分填充所有电子状态的亚空间。但是,局域化轨道的选择并没有像我们期望的那样重要。d或者f轨道,这些轨道需要考虑库伦关联作用,事实上他们也是很好地局域化空间,能够很好地保留固体中原子的性质,而使用LDA+U方法进行电子结构模拟的经验表明计算结果对于特定的局域化轨道而言并不灵敏。

  由于在式(3-22)中,哈密顿算符定义中缺少射影算符,最直接的计算方法应该是原子轨道式基集,例如LMTOs。但是,只要局域化d或f轨道被定义,式(3-22)中的哈密顿量就可以使用平面波波函数作为基集,比如赝势方法。

  3.4.4强关联体系计算方法---LDA+DMFT

  在上一节中,LDA方法在多体问题中拓展成LDA+U方法已经被我们进行了深入讨论,但是在LDA+U方法中,电子的库伦作用项是在Hartree-Fock近似中被描述的,因此LDA+U方法并不能代表一个真实的多体物理图像。虽然这种方法成功地描述了强关联体系中长程有序的,孤立的电子态,但是它未能给出一个强关联电子态的相关参数。因此,进一步优化LDA+U方法,从而获得多体问题中关于电子与电子间的相互作用的数值是一个急需解决的问题。而解决该问题的途径很多,但其中最著名,最先由Anisimov开始应用的方法,就是动力学平均场近似方法(DMFT),也被称为LDA+DMFT方法。在以LDA为基础进行拓展的所有方法中,只有LDA+DMFT方法可以描述一个考虑强关联效应的顺磁金属的物理性质,这种金属通常还具有较高或较低的赫伯德能带以及费米能级上的准粒子狭峰。而这种特征性的三峰结构是多体效应的一个重要标志。

  过去的几十年的时间中,DMFT方法被证明是一个探究带着局域库伦作用的强关联电子体系的物理性质的有效的方法,其计算结果在高的晶格配位数中十分精确以及保留局域近似的的动力学特征。因此,其才被称为动力学平均场近似。在这种方法中,有关晶格的问题是被描述在一个有效的单格点问题中(如3-1),其是由对k积分的与系统自身能量和在频率的格林函数有关的Dyson方程的自洽计算决定:

  (3-25)

  (3-26)

  其中,1表示单位矩阵,表示化学势,也被称为费米能级,表示原胞中相互作用的原子,表示k空间中的一个矩阵元,表示只有在具有相互作用的轨道上才是0的自能矩阵,[…]-1表示具有矩阵元和的矩阵的转置,而(3-25)式的积分是覆盖整个体积为的布里渊区。

  图3-1

  该图说明当临近格点的数目趋向于无穷时,中心极限定理和格点与格点间位置的波动可忽略。这意味着邻近格点的影响可被一个有效的平均场所取代,该动力学平均场可由自能进行描述。因此,该动力学平均场问题可等效于对k积分的Dyson方程的自洽计算和下列多能带Anderson求解器模型。

  DMFT单格点问题依赖于

  (3-27)

  该问题在杂化满足于

  (3-28)

  等效于Anderson求解器模型[19-21]。松原频率

  (3-29)

  (其中表示温度的倒数)上和磁量子数为m()以及自旋量子数为处的局域单粒子格林函数由下列基于Grassmann变量和的积分函数给出:

  (3-30)

  其中,

  (3-31)

  表示配分函数,单格点行为A具有下列形式(A的交互部分由虚时表示,例如的傅里叶变换)

  (3-32)

  式(3-30)中的单格点问题必须与以k为积分变量的Dyson方程一起进行自洽求解,来得到一个给定问题的DMFT解,如图3-2。

  选择初始自身能量

  从中计算G通过式(3-25)

  从g中计算G通过DMFT单格点问题,即式(3-30)

  整个过程利用进行迭代直到收敛,例如

  图3-2

  该图表示DMFT自洽周期的流程

  由于DMFT单格点问题和Anderson杂化问题的等价性,一系列的近似方法可以被应用到求解DMFT方程,例如迭代扰动理论(the iterated perturbation theory,IPT)[22,23]和非交叉近似(non-crossing approximation,NCA)[24],此外,也被应用到量子蒙特卡洛方法(quantum Monte Carlo method,QMC)的定量分析,精确对角化(exact diagonalization,EQ)以及定量的重整化群(numerical renormalization group,NRG)[23,25-26]上。QMC和NCA方法将要在下一章节中被描述,而IPT是表示非自洽的二次扰动理论,针对的是在半满过程中的Anderson模型中的U值的求解。它代表了一种拟设,即持正确的的扰动U2项和正确的非半满状态下的原子自洽能的限制[27],而ED方法可以直接在有限格点和轨道数目中对角化Anderson杂化问题。NRG首先利用一系列离散的在(D为带宽,n=0,1,2…,NS)处量子态来取代导带,然后反复地将在能量较低处,利用增加的精确度将该问题进行对角化,例如,增加NS的取值。原则上来讲,QMC和ED都是一种精确的方法,但是他们都必须附带外推法,例如在QMC中,虚时的离散度和在ED中的杂化模型的晶格格点数目。

  在LDA+DMFT背景下,我们倾向于利用上述的计算思路来求解之前讨论过的DMFT方程,具体表达为LDA+DMFT(X),其中X=IPT,NCA,QMC,最典型的例子就是La1-xSrxTiO3。同样的策略也被Lichtenstein和Katsnelson所使用来作为他们的LDA++方法之一。Lichtenstein和Katsnelson应用LDA+DMFT(IPT),他们也第一次运用了LDA+DMFT(QMC)来研究铁的光谱性质。近年来,V2O3,Ca2-xSrxRuO4,Ni,Ru和Ce等物质的相关性质均已被LDA+DMFT方法所研究。而巡游电磁体的实体调查,例如Ni,通过将密度泛函理论和Gutzwiller波函数结合,最近也成为了可能。

  3.4.5 QMC方法求解DMFT

  DMFT的自洽循环要求一种能够求解DMFT中的单格点动力学方程的方法,如式(3-30)。由Hirsch和Fye建立的QMC算法是一种很好的定量精确求解Anderson杂化模型的方法,此外,它也允许我们求解在处的杂化格林函数G和强关联方程。事实上,QMC是在一系列无相互作用问题之和上描绘电子相互作用问题(式(3-30)),其中单粒子始终在一个波动的,与时间关联的场中运动,具体参见图3-3。最后,泛函积分式(3-30)的虚时间隔被分散进大小为

  (3-33)

  中的步,得到了一些支撑点:

  (3-34)

  其中,。式(3-30)的指数项可以通过关于算符和算符的Trotter-Suzuki公式而被分开表[28]:

  (3-35)

  上式可在中被精确表达。式(3-32)中的单格点表达式A也可以利用离散的虚时写成下列形式:

  (3-36)

  其中,第一项是从松原频率到虚时的傅里叶变换。第二步,在单格点表达式(3-32)中,项可以通过引入一个经典附加场来实现:

  (3-37)

  其中,双曲余弦

  (3-38)

  以及M代表相互作用的轨道数量,这被称为离散Hirsch-Fye-Hubbard-Stratonovich转变,可以被应用到库伯排斥作用和洪德规则耦合的Z轴分量上。现在,这种泛函积分可以被简单的高斯积分所求解,因为费米算符仅仅是二次化的,例如,由于一个被给定的附加场结构,整个系统是无相互作用的,然后量子力学问题就可以简化成一个矩阵问题:

  (3-39)

  其中,Z为配分函数,矩阵

  (3-40)

  以及矩阵的元素

  (3-41)

  其中,

  (3-42)

  因为式(3-39)的和是由因数组成,关于大Λ的计算是不可能的,因此在高维求和和积分是有效的蒙特卡洛方法在式(3-39)中起着重要作用。在这个方法中,积分被分裂进归一化的概率分布P和余项O:

  (3-43)

  其中,

  ,(3-44)

  在统计物理学中,玻尔兹曼分布可以作为P的一种选择:

  (3-45)

  对于式(3-39),这种概率分布可转化为:

  (3-46)

  而余项也可表示成:

  (3-47)

  蒙特卡洛模拟生成了与概率分布有关的结构,以及在利用来平均化可观察量,而不是计算所有可能的结构。因此,具有很高的玻尔兹曼权重的相空间中的相关的部分相比于权重较小的相空间而言是更广泛的,创造了在蒙特卡洛方法中的重要的样本。在中心极限理论中,该理论给予了N个统计学上独立的参数,可以利用下列近似:

  (3-48)

  选择任意附加场结构

  利用式(3-47)计算当前格林函数:

  利用式(3-40)中的M以及引入类似于(3-27)形式的下列表达式:

  预热

  蒙特卡洛扫描:MC-sweep

  测量

  扫描

  蒙特卡洛扫描:MC-sweep

  图3-3

  QMC库通过使用MC-sweep程序来计算格林函数矩阵G的流程图

  选择乘以一个集合,定义为s,除了这些元素,其具有相反的符号

  利用下列式计算翻转频率

  然后计算式(3-40)中的M

  随机数?

  是否

  ;根据式(3-47)计算保持s

  图3-4

  其中,误差与所需因数的数量几乎是是与积分次数无关的,因此蒙特卡洛方法的计算仅仅是多项式阶数随积分次数上升,而不是作为积分的指数。利用Markvo过程以及在附加场中的单自旋翻转,该算法的在Λ上的主要阶数的消耗是

  (3-49)

  其中a为单自旋翻转的接受速率。

  QMC方法的优势是他在数值上非常精确,允许我们计算单粒子格林函数以及双粒子(甚至更高)格林函数。在现有的工作站中QMC方法能够在温度高于室温的情况下,处理高达7个的相互作用轨道,但是在低温下却不可以,因为数值会随着温度发生变化,其关系为。因为QMC计算或时会带有一定的误差,因此它也要求熵的极大值化来得到在真实频率下的格林函数。

  3.4.6 NCA方法求解DMFT

  NCA方法是一种分散扰动理论,主要应用于有效Anderson杂化问题中的杂化参数,因此,如果库伯相互作用参数U相比于带宽而言是很大的话,这种方法是很可靠的,而且也提供一种在计算上很方便有效的方法来检查其他情况下的普通光谱特征。为了了解NCA是如何利用DMFT进行计算,我们重新定义式(3-25)

  (3-50)

  其中

  (3-51)

  再次说明,,和之后的和都是轨道空间中的矩阵。注意仅仅对于计算轨道来讲有非零通道。在基态,式(3-50)可以被改写为

  (3-52)

  其中

  (3-53)

  K点的数目为以及

  4alps软件

  4.1 alps软件介绍

  Alps软件,全称为Algorithms and Library for Physics Simulations,即物理模拟计算中所用到的算法和库,是一款为强关联量子体系提供高端模拟代码的开源软件,同时也提供了C++库来优化该类代码。此外,alps所依托的编程语言是Python,其具有内容简洁,运算效率高等特点。

  4.2 alps软件的使用指导

  4.2.1 alps程序运行的前提准备

  在alps软件中,模拟计算需要调动一系列的库,但是使用这些库进行计算的前提是使用者能够提详尽的参数,即使这些参数具有多重含义,例如一个物理系统中的几个温度参量。然后在程序运行过程中,从单个参量开始,而不是多个参量并行计算,此外,该程序还使用一系列的检查点来防止在计算超时时的数据丢失。因此,这些被调用的库要求一个文件来作为参数文件,来明确地指出每一个参与将要运行的蒙特卡洛模拟的参数。而在alps中,一个基本的参数文件的格式如下:

  <SIMULATION>

  <PARAMETERS>

  <PARAMETER name="L"

  <PARAMETER name="SWEEPS"

  <PARAMETER name="T"

  <PARAMETER name="THERMALIZATION"

  </PARAMETERS>

  </SIMULATION>

  另一方面,在模拟计算之前,除了要给出相应的参数文件,还应该要给出工作文件,来调用alps的各类库函数进行计算,并生成最终的输出文件。而工作文件的准备,也是有一个明确的书写格式,并且符合Python的语法规则,下面也给出一个基本例子:

  <JOB>

  <OUTPUT file="parm.xml"/>

  <TASK status="new">

  <INPUT file="parm.task1.in.xml"/>

  <OUTPUT file="parm.task1.xml"/>

  </TASK>

  <TASK status="new">

  <INPUT file="parm.task2.in.xml"/>

  <OUTPUT file="parm.task2.xml"/>

  </TASK>

  <TASK status="new">

  <INPUT file="parm.task3.in.xml"/>

  <OUTPUT file="parm.task3.xml"/>

  </TASK>

  </JOB>

  综上所述,利用alps进行强关联量子模拟的过程可分为三步:

  (1)准备输入文件;

  (2)运行模拟的程序;

  (3)分析输出的计算结果。

  4.2.2 alps中如何应用DMFT

  动力学平均场理论(DMFT)针对量子多体问题提供了一种近似的求解方法。在这种方法中,系统局域态的物理性质可以被很精确地讨论,但是空间中的相关性却会被忽略。DMFT被应用的最广泛的一个方面就是强关联材料的有关计算,其典型的方法是LDF+DMF。在这种近似之中,晶格问题是在包含时间的杂化模型上被描述的,同时该模型也能够进行自洽计算。含时的杂化模型可以被杂化求解器求解,而alps中,通常提供了两种这样的杂化求解器的相关代码:the hybridization expansion algorithm和the interaction expansion algorithm。

  4.2.3 The hybridization expansion algorithm

  在这一节中,我们将运行一个连续时间蒙特卡洛模拟的例子,来说明the hybridization expansion algorithm。下面,我们将要计算贝斯晶格中的赫伯德模型计算,其相互作用势为:

  (4-1)

  其中,D表示半带宽宽度。该模型的电子填充状态为半满,冷却后进入反铁磁相。计算该模型的完整代码为:

  import pyalps

  import numpy as np

  import matplotlib.pyplot as plt

  import pyalps.plot

  #prepare the input parameters

  parms=[]

  for b in[6.,12.]:

  parms.append(

  {

  'ANTIFERROMAGNET':1,

  'CONVERGED':0.003,

  'FLAVORS':2,

  'H':0,

  'H_INIT':0.03*b/8.,

  'MAX_IT':6,

  'MAX_TIME':300,

  'MU':0,

  'N':250,

  'NMATSUBARA':250,

  'N_MEAS':10000,

  'OMEGA_LOOP':1,

  'SEED':0,

  'SITES':1,

  'SOLVER':'hybridization',

  'SC_WRITE_DELTA':1,

  'SYMMETRIZATION':0,

  'U':3,

  't':0.707106781186547,

  'SWEEPS':int(10000*b/16.),

  'THERMALIZATION':1000,

  'BETA':b

  }

  )

  #write the input file and run the simulation

  for p in parms:

  input_file=pyalps.writeParameterFile('parm_beta_'+str(p['BETA']),p)

  res=pyalps.runDMFT(input_file)

  接下来,就让我们来详细分析每一个参数的物理意义:

  Antiferromagnet:表示该模型是否为反铁磁相,若是则设置成1,反之为0;

  Converged:表示迭代收敛准则,也就是代表一个判断是否发生收敛的依据;

  Flavors:表示电子的自选状态,若为0则表示只考虑自旋向上,为1则表示只考虑自旋向下,为2则表示两种状态均要考虑;

  H:表示在量子轴方向的磁场,一般取0;

  MAX_IT:表示迭代次数的上限(一旦达到收敛准则,自洽计算会提前结束);

  MAX_TIME:表示杂化求解器中每一个迭代的计算时间上限,单位为秒;

  MU:表示化学势或费米能级,即电子填充状态。若为0,则表示在该模型中电子填充状态为半满;

  N:表示虚时格林函数的辅助离散度;

  Nmatsubara:表示松原频率的截止频率;

  N_MEAS:表示测量间的更新次数;

  OMAGE_LOOP:设置为松原频率下运行的自洽计算;

  Seed:蒙特卡洛随机数种子设置;

  Sites:杂化格点数,对于单格点DMFT模拟而言,其值为1;

  Solver:表示应用何种杂化求解器;

  SC_WRITE_DELTA:表示在书写杂化求解器的hybridization函数时的自洽选项;

  Symmetrization:表示求解模型的对称性;

  U:表示相互作用强度,一般取3。其值若为0,则表示该材料为金属导体,反之则为绝缘体;

  T:表示跳跃参数;

  Sweeps:表示整个计算过程中的扫描总数,但如果达到了计算时间上限MAX_TIME,计算也会停止;

  Thermalization:表示热化扫描;

  Bata:表示温度的倒数;

  在准备好输入文件之后,在Python窗口中,执行命令:

  Alpspython tutorial.py

  等待一段时间之后,就能够得到对应的计算结果。而在这个例子中,我们可以得到下列图4-1。

  图4-1

  该图表示在和(从底部到顶部)处,通过L=16的量子蒙特卡洛方法得到的关于半填满赫伯德模型(包含一个半带宽为D的半圆态密度)格林函数中的自洽解。在这个计算中,当(无法被体现)时,用到的自洽条件允许反铁磁命令。