首页     直问学姐学长     专家     文章          软件     小组     视频    

注册  |  登录
首页 > 文章 > DFT+U中U参数的确定(转)

DFT+U中U参数的确定(转)

作者:朝圣   (离线)   [DFT+U ]   时间:2016-05-23 08:46:59  向他请教

11.28 更新:加入插图和说明;更改对Marzari自洽U方法的介绍;
11.29 更新:修改个别表述和链接;

最近在做DFT+U计算,有些收获分享一下。本人数学功底先天不足,难免有理解错误的地方,欢迎指正。

首先是对+U必要性的感悟。DFT(LDA和GGA)对一般体系的结构、能带计算还是很令人满意的。这些一般体系是指前3周期的体系和纯金属体系。但是,对于包含d电子和f电子的体系,特别是过渡金属氧化物或氮化物,在“金属/绝缘体”的定性判断上常常出错。LDA和GGA往往把一些绝缘体/半导体的gap算的太小,甚至最高占据轨道/能带(HOMO)在Fermi面之上,即成金属了。这些最高占据轨道往往是来自这些金属原子的d电子和f电子。那么为什么DFT处理d/f电子会出现较大误差呢?我认为这些d/f电子的特殊性在于它们不仅离原子核较远,轨道空间伸展形状奇特(比如穿透效应,即穿透更加外层的s,p轨道)而且多数会发生自旋极化,即原应自旋相反的电子对发生自旋跃迁,导致金属原子神奇地出现了净的磁矩。如果单胞内所有金属原子的净磁矩方向相同,则单胞总体产生净磁矩,材料呈现铁磁性;如果单胞内相邻金属原子的净磁矩相反,则单胞没有总体净磁矩,但存在一上一下的磁矩排列,材料呈现为反铁磁性。正是d/f电子这些神奇的特性,使得DFT计算的过渡金属氧化物或氮化物能带结构频频失误。为什么呢?原因在于DFT计算电子-电子之间的排斥作用所用的交换-相关泛函,是基于单粒子近似发展起来的。单粒子近似是把自旋配对的电子对看成单个粒子,这样Kohn-Sham轨道数量为总电子数的一半。但是当处理d/f电子时,不得不把电子对分开为自旋向上和自旋向下的两个粒子,即开壳层计算open shell。对开壳层体系使用单粒子近似下的交换-相关泛函时,其中的相关泛函就不能完全考虑d/f电子的相关效应。d/f电子的相关效应,就是自旋电子与自旋电子之间的排斥,和自旋电子与内层电子之间的排斥。这些排斥作用使得轨道/能带之间分割较远,轨道/能带较窄。但是DFT的相关泛函对这些排斥考虑的不够,结果轨道/能带过宽,轨道与轨道相互接近甚至重叠。这就是为什么DFT把绝缘体计算成了金属。一般把这些过度金属半导体/绝缘体体系称为强相关体系strongly correlated system。

怎么办呢?可能你已经直觉到了,直接加上一项,把d/f轨道上的自旋电子间的排斥作用作为(u/2*int(n_up*n_down))一个能量项加进去。这个就是Hubbard U model:

E_DFT+U = E_DFT + U/2*int(n_up*n_down)
其中第一项E_DFT仍为正常的DFT计算。n_up, n_down分别代表自旋向上和自旋向下的电荷密度。二者相乘表示了它们之间的排斥作用。完整的Hubbard model还要考虑其中的重复计算,减去一个项,称为J项。+J项的计算刚刚在QE中实现。如果你追求结果的定量准确性,可以尝试玩玩+J。

这个+U项排斥作用前面有个系数U/2,也就是说这个排斥能并非直接的库仑排斥作用,而是有个系数。这个系数并非一个经验参数,而是代表了能量对d/f轨道上电子密度变化的二级导系数【1】。你可能会困惑,这有点像GGA?是的,不过GGA是通用的交换-相关能对密度的二阶展开的泛函。而这里的U是纯DFT基态总能量对特定体系d/f电子密度变化的二阶导。这有点像是个二次方项的收缩函数,即你们这些自旋作用的电子在轨道上的占据必须往特定方向收缩,否则会有能量惩罚。


上图中,不加U时左边的最高占据轨道和右边的空轨道之间的gap很窄,是不正确的;+U后变宽,更加接近实验结果。需要注意的是,+U是对电子结构的改善,对原子几何、晶体结构的影响可能很有限。此时只要不涉及电子性质,一般的GGA, LDA的结果是可以接受的。


那么这个+U能量惩罚项是人为的经验项,还是真实的物理存在呢?Cococcioni曾在【1】PRB71(http://prb.aps.org/pdf/PRB/v71/i3/e035105)中试图用图1来说明Hubbard U能量项的本质。他说真实的体系与环境交换电子,其能量变化应该是交换前后的两个状态的线性组合,这样的能量变化应该是线性的。但LDA和GGA对部分(非整数的)占据的Kohn-Sham轨道的自作用的不正确处理导致能量变化是非线性的。这二者之间的能量差就是Hubbard U能量项的来源。从这个解释可以看出,

(1) +U能量惩罚项是物理存在的,并非经验校正。参数U的确定也不应该是经验的,而是必须考察DFT与真实势能面的差。

(2) 这个二级导是体系环境变化而变化的。这就导致必须对每个体系,甚至体系的不同状态(比如键长变化)计算相应的参数U。注意如果体系内有对称性不相同的金属位点,需要对每个位点上的金属原子计算其相应的参数U,不管它们是否为同种金属!获得合适的参数U以后,才能使用DFT+U模型精确确定自旋极化和能带、带隙等问题。(当然,更精确的模型还要考虑旋轨耦合spin-orbital coupling。)

怎么确定某个体系参数U呢?因为U是纯DFT基态总能量对这个体系d/f电子密度变化的二阶导,那么很自然的一个想法,我们对密度加以扰动,求能量的变化(响应),可以得到数值二阶导(微扰法)。不过,密度并非一个理想的微扰项,因为一个金属原子(Hubbard位点)涉及到多个电子。怎么办呢?

Cococcioni提出了一种线性响应法。它向单个Hubbard位点施加一个势alpha,然后引起该位点的密度(占据数)重新分配,并得到纯DFT(即不含U,或者使U近似为0)能量对此alpha势变化的响应。这种总能响应包括2部分,1是能量对势本身的线性响应(占据数发生变化前)称为bare response(X0);2是能量对加势和占据数变化后总的响应,称为self consistent response(X)。二者的线性响应系数的倒数之差即为U:

U = 1/X0 - 1/X  
此法原始文献即PRB71. 




图2,
上图中,红色为bare response,其线性系数为X0=-0.32,数据从自洽前第一个循环计算得到的Hubbard site上(即带d/f电子的金属原子)的占据数得到;蓝色为self-consistent response, 线性系数为X=-0.15,数据从自洽完成后Hubbard site的占据数得到. U = 1/X0 - 1/X = 3.47eV
 


在这个教程里有很详细的例子: 点击打开并下翻到
July, 23 Tursday: LSDA, collinear and noncollinear magnetism; spin-orbit coupling; DFT+U
在其中的Labs 1.  Matteo Cococcioni,Examples of spin-polarized and DFT+U calculations

注意上面这种方法基于一个假定,即势alpha足够小,此时得到的U近似等于alpha=0时的U,即U_0。这是线性响应理论。即体系能量实际对势alpha的响应并不是线性的,但我们假定扰动足够小的情况下假定这种响应是线性的。上述方法对单个Hubard位点施加势alpha,但是在周期性体系中,任何扰动也都是周期性的。怎么办呢?

Cococcioni提出的方法是延伸单胞为超胞,做响应计算;然后再次延伸为更大超胞,做响应计算。。。这样随着盒子越来越大,势在其中的作用越来越小,直到盒子无限大,则响应系数收敛。这种方法并不需要计算实际超胞,而是对密度响应进行外推。

Marzari提供了另一种不通过超胞外推即可获得无扰动的U_0的自洽方法(self-consistent U)(原始文献PRL97,103001(2006).)注意到尽管我们需要的是纯DFT能量对alpha的响应,但是上面的所有计算其实都是按照DFT+U的模型计算,只不过,进行alpha势扰动的时候,设置U为一个很小的值(称为初始U值,U_in),以获得纯DFT基态能对alpha的响应,从而得到U(称为U_out)。

Marzari并不让初始U_in近似为0, 而是设置它等于一系列值(0.5, 1.0, 1.5, 2.0, 2.5, 3.0 ...),然后运用上述方法去计算相应的U_out,然后将上述(U_in, U_out)数据系列在线性相关区拟合直线,可外推到U_in=0时的U_scf。如下图所示。




图3,
上图中,对给定的U_in = 0.5, 1.0, 1.5, 2.0, 2.5, 3.0做DFT+U计算,结束后可根据Hubbard位点上占据数得到一个U_out = Uscf - U_in / m。然后外推U_in=0时,U_scf=3.56即为正确的参数U。其中的U_0是使用第一种方法在势alpha扰动近似为0得到的近似U_0值,并非alpha=0的U值。


第一种超胞法(称为linear response法)外推和U_in/U_out外推法得到的参数U应该是一致的,但第二种方法更加优雅。

这个方法在这里有很好的教程: Calculating the Hubbard U
其中对两种方法都有相应的例子和Python脚本。注意两种方法的切换要修改其中的Variable.py文件,设定hubbylist或者alphalist。其中所用的例子是6重态MnO的cluster,并非固体。

虽然参数U的确定很麻烦,但是还是值得的。从PRL 97, 103001 (2006)中可以看到,GGA+U后算FeO的平衡键长、简谐振动频率,非谐频率的结果,直接逼近CCSD(T),当时哥就震惊了。

后记:MS和Castep的朋友不要问我参数怎么算,我也不知道。我只会用PWScf即quantum espresso。欢迎补充。

在使用+U model的过程中,我一直在思考这个U参数的本质,即U大一点意味着什么,小一点意味着什么?我们想,既然它是DFT能量能量对电荷密度变化的二阶导,那么有点类似与振动频率。在分子中,如果一根键的强度很高,那么振动频率就很高;如果键强度低,振动就弱,比如C=C双键就比C-C单键的振动频率高。类比的话,应该说U代表了d/f电荷密度随外部势(比如其它原子核对它的作用)的变化而变化的难易程度。如果U大,那么不容易变化;如果U小,则容易变化。在不使用+U模型时,那就是说U=0,这时d/f电子最容易受外部势的影响而左右摇摆,从宏观上来看,就是分布比较均匀,能带很宽。一旦+U以后,立即受到束缚,从而能带变窄。这和+U的实际作用正好一致。

进一步,Hubbard原子上的U参数与材料的金属性-非金属性的有一种联系:材料的金属性越强,U越小;材料非金属性越强,U越大。这从Cococcioni博士论文中的数据可以看出来:在金属态中,铁Fe的U为2.2, 铁磁FeO中为4.29, 而绝缘体Fe2SiO4中为 4.9/4.6. 所以,如果你得到的U很小,很可能意味着你的体系的金属性很强。当然,这并不意味着你可以不使用U,U虽小仍然可能大大改善能带的分布。

"在不使用+U模型时,那就是说U=0,这时d/f电子最容易受外部势的影响而左右摇摆,从宏观上来看,就是分布比较均匀,能带很宽。一旦+U以后,立即受到束缚,从而能带变窄。"
这句话容易引起歧义,补充一下,U越大,能带越宽,对应带隙越小;U越小,能带变窄,对应带隙越大。不过这只是一个简单理论分析,是否真的对所有材料都是如此,由于本人研究过的体系有限,不敢下定结论。欢迎补充。

[ Last edited by franch on 2013-1-29 at 13:48 ]

帖子里的Stanford的教程: Calculating the Hubbard U
已经转移到了MIT: http://hjklol.mit.edu/content/calculating-hubbard-u

[ Last edited by ChemiAndy on 2014-10-14 at 05:08 ]

Calculating the Hubbard U的参考文献PRL91直接下载: http://pan.baidu.com/s/13x4x4

发表提问

您有0条新消息