如何考察结构能带

1、如何考察一个能带（DOS）结构和复杂的相互作用Part 1 Electric conductivity and Band structures固体计算最终结果将以能带结构展示出来，关于能带结构，固体中化学键分析，轨道之间的相互作用的解释等是一个复杂的过程，这里只是简单的根据本人的经验对此作定性的描述.根据Fermi面附近能带的分布情况，固体分为绝缘体（insulator），半导体（semiconductor），导体（conductor），导体比较典型的是金属，能带在Fermi面附近是连续分布，主要由于金属d，s以及p轨道之间能级重叠导致了Fermi面能带的联系分布，金属电导的好坏不仅仅是看Fe

2、rmi附近是不是存在可供电子跃迁的能级，还要看这些能级是不是扩展态（extended  or delocalized states），如果是定域态（localized）那么及时Fermi附近呈现Metallic特性，电导不会比金属好，比如过渡金属化合物电导就要比金属本身差很多。过渡金属本身电导也会受到d轨道扩展程度的影响，比如3d系列Fe，Co，Ni等电导率不是很大，比起Cu，Ag等就差的远了，对于Fe等金属Fermi面主要陈分是3d轨道，而对于Cu和Ag，由于3d（4d）轨道已经成满层排列，因此Fermi面落在了扩展的s轨道上面，这些轨道上的电子类似于自由电子气，能带呈现

3、抛物线的形式，E(k)=h2k2/2me; 具有比较高的电导率，相反Fe等的3d轨道成分也可以分为巡游电子（自由电子，轨道为扩展性，能带呈现抛物线特点）和定域轨道两类，定于轨道能带在k空间是离域的，色散关系比较平直，但在晶体实空间内高度的定域，受到原子核的Coulomb吸引作用比较强烈，难以发生迁移，因此如果填充电子落在这些d轨道上面，电导性会大大降低。当然具体取决于DOS或者能带是如何分布的，这个和晶体结构有关系。在一些化合物中如TiC等结构中，Fermi面最后落在以C2p轨道为主要成分的能带上面，p轨道主要参与结构共价键形成，这些电子能级一般定域在Ti和C原子周围，电子处于紧束缚状态，难以

4、在外加电场下发生迁移，因此这时候化合物的电导会进一步下降。 Part 2 关于半导体能带的特点：半导体能带类似于绝缘体，区别在于带隙数值，一般认为宽带隙半导体的能带最大在4eV左右。如果比这个更大，可以认为是绝缘体。半导体能带主要分成三个部分：valence band， band gap， conduction band。Valence band：主要由电负性较大的原子组成，如InP，价带主要是P的3s，3p轨道，导带一般是金属原子组成，如In的s，p轨道等。从化学键角度考虑，价带一般是Bonding，当然也有部分结构表现出Antibonding状态。Band gap： 带隙，电子动

5、量空间的绝对真空地带。Conduction band： 电负性较小的原子轨道组成，一般是反键态为主，也可能由Nonbonding states。半导体能带需要注意the top of valence band和the bottom of conduction band，确定他们在Brillouin 区位置有助于了解带隙特点，如直接带隙（direct band gap） 还是间接带隙（indirect），不同的带隙特点在发光效率方面由比较大的影响，由于间接带隙涉及导声子驰豫过程，因此发光效率远远不如直接带隙半导体。 带隙比较大的固体结构，能带结构类似于半导体，能带的色散关系和固体结构中

6、化学键的类型有关系，比如共价键强烈的固体能带扩展能量范围很广，如果存在离子键，能带会比较平滑一些。在磁性材料中需要考虑交换分裂对能带的影响，一般的DFT计算，轨道电子占据如果是非自旋极化的，电子内禀属性spin是不考虑的，不同自旋的电子占据同一轨道的概率相同，但如果是考虑自旋极化，轨道的这种对称性将被打破，能量的在最小化过程中可以让某些轨道保留多于的自旋极化电子，电子的磁矩可以保留下来。自旋极化对应的交换能（原子之间）一般用Heisenberg模型：EJS1S2CosaJ为交换系数，可正可负，Si为自旋角动量，a是自旋角动量之间矢量夹角。一般轨道自旋极化产生的磁性用Stoner 模型。在Sto

7、ner模型中整个体系的电子总数是N0，不考虑自旋的时候up spin 和down spin电子数目都是No/2，如果考虑自旋极化，自旋极化率为s，则up spin（1s）×No/2，对应自旋down spin（1s）×No/2。交换能存在于自旋相同电子之间，假设为U，则E（up spin）E（k）U×N（up spin）；同理E（down spin）E（k）U×N（down spin）；显然U<0时候，E（Up） Hubbard 参数U和Stoner 自旋极化参量I 在平均场近似或者说是一般得LDA计算中，能带得自旋分裂是由一个

8、称为Stoner 参数I来主导得，而且平均场近似中认为这个交换分裂能是小于带宽得，一般而言I得数值在LSDA中大概在1eV左右，这样即使对于某些氧化物如NiO，FeO，MnO等即使能带发生交换分裂，整个体系依然是金属性得，试验观察到得这些氧化物实际上是绝缘体，主要原因在于在这些氧化物中d轨道得能级位置不是由I来决定得，而是Hubbard参数U决定，U也称为Onsite Coulomb作用能，相当于把两个电子放在空间同一个位置需要得能量，U数值一般在10eV左右，如此大得分裂能足以将Fermi面附近连续分布得d能级分开，从而得到正确得基态性质，目前广泛采用得LDAU的算法就是针对某些定域轨道，如

9、d或者f，这些轨道占据存在强烈的在位Coulomb排斥。正是U使得带隙分裂，而不是轨道极化参量I。U计算目前仍然是一个研究热点，最近由人采用线性响应理论，同时结合轨道束缚的DFT计算自洽的求解了部分体系的U参数，在大部分情况下U数值要根据计算结果和试验参数的符合程度而定。LDAU计算核心思路是：首先将研究体系的轨道分隔成两个子体系（subsystem），其中一部分是一般的DFT算法（如LSDA，GGA）等可以比较准确描述的体系，另外是定域在原子周围的轨道如d或者f轨道，这些轨道在标准的DFT计算下不能获得正确的能量与占据数之间的关系（如DFT总是认为分数占据是能量最小的，而不是整数占据）；对于

10、d或者f轨道，能带模型采用Hubbard模型，而其他轨道仍然是按照KohnSHam方程求解；d以及f轨道电子之间的关联能采用一个和轨道占据以及自旋相关的有效U表示；整体计算的时候需要将原来DFT计算过程中已经包含的部分关联能扣除，这部分一般叫Double Counting part，并且用一个新的U来表示，最终的结果是在DFT计算的基础上新增加一个和d或者f轨道直接相关的分裂势的微扰项，这部分能量可以采用一般微扰理论计算。在CASTEP最新的版本中增加的LDAU的计算，U参数的设置一般主要是针对过渡金属氧化物（Charge transfer type insulator），包含非满层f轨道的元

11、素等，高温超导体强关联体系。在参数设置方面主要是需要注意d和f轨道，至于s以及p轨道一般不需要设置，当然由文献也报道p轨道的这种关联性。过渡金属氧化物的有效U如下：Species         U              J                    UJ（Ueff）NiO 

12、;             8.0           0.95                  7.1CoO            7.8          

13、;  0.92                 6.9FeO             6.8            0.89                 5.9MnO 

14、60;         6.9            0.86                 10.3VO               6.7           &#

15、160;0.81                 5.9TiO              6.6            0.78                  5.8Refer

16、ence: Band theory and Mott insulators: Hubbard U insteat of Stoner I, PRB Vol44 No 3 (1991);3d轨道U和J计算如下所示：主要原理是改变d轨道的占据，在自旋极化的前提下计算不同自旋轨道能量的差值提取U和J，U微Coulomb排斥能，J是交换能，U在所有电子中都存在，不管自旋是否相同，J只存在于自旋相同的电子上。下面给出过渡金属（不包括稀土元素）U和J参数的选取(uint in Ry,1Ry=13.6eV)：3d series:Elements       

17、   U                  JV                   0.25               0.05Cr             

18、;      0.26               0.053Mn                  0.28               0.055Fe           &#

19、160;       0.3                 0.058Co                  0.31                0.059Ni      &

20、#160;             0.31                0.064d series:Elements          U                    JNb  

21、0;                 0.19               0.04Mo                    0.2            &

22、#160;    0.04Tc                    0.21                0.042Ru                    0.22  &#

23、160;             0.042Rh                    0.25                0.044Pd               

24、0;    0.29                0.0445d series:Elements           U                     JTa         &

25、#160;          0.19                0.039W                      0.20               0.038Re 

26、                   0.205              0.039Os                    0.2           &

27、#160;      0.039Ir                       0.21                0.038Pt                  

28、0;  0.215              0.038对于其他过渡金属化合物U一般在510eV之间。如在PRB73，134418（2006）这个文献中作者在计算Co掺杂的ZnO时采用的U是6和8eV。过渡金属的U数值和d电子排列以及价态有关系，因此上面给出的数值只是一个大概的估算数值，具体文献见Physical Review B Vol50,No23,1994.LDA+U 算法主要原创作者是俄罗斯金属研究所的V.I. Anisimov，重要文献有：Corrected atomic l

29、imit in the local density approximation and the electronic structure of d impurities in Rb, Phys.Rev.B 50,23 (1994）；Band theory and Mott insulators: Hubbard U instead of Stoner I, Phys.Rev.B. 44 No.3 (1991);Materials Studio 4.3版本中也给出了一些元素默认的U数值（实际是UeffUJ）：Element Name        Atom

30、ic number        Angular Momentum        Hubbard USc        21        d        2.5 eVTi        22        d        2.5 eVV     

31、;   23        d        2.5 eVCr        24        d        2.5 eVMn        25        d        2.5 eVFe        26

32、0;       d        2.5 eVCo        27        d        2.5 eVNi        28        d        2.5 eVCu        29     

33、;   d        2.5 eVY        39        d        2.0 eVZr        40        d        2.0 eVNb        41        d 

34、       2.0 eVMo        42        d        2.0 eVTc        43        d        2.0 eVRu        44        d     

35、  2.0 eVRh        45        d        2.0 eVPd        46        d        2.0 eVAg        47        d        2.0 eVCd

36、        48        d        2.0 eVLa        57        f        6.0 eVCe        58        f        6.0 eVPr    &

37、#160;   59        f        6.0 eVNd        60        f        6.0 eVPm        61        f        6.0 eVSm        62

38、        f        6.0 eVEu        63        f        6.0 eVGd        64        f        6.0 eVTb        65    &

39、#160;   f        6.0 eVDy        66        f        6.0 eVHo        67        f        6.0 eVEr        68        f&

40、#160;       6.0 eVTm        69        f        6.0 eVYb        70        f        6.0 eVFr        87        f    &#

41、160;   2.0 eVRa        88        f        2.0 eVAc        89        f        2.0 eVTh        90        f        2.0

42、 eVPa        91        f        2.0 eVU        92        f        2.0 eVNp        93        f        2.0 eVPu  

43、0;     94        f        2.0 eVAm        95        f        2.0 eVCm        96        f        2.0 eVBk       

44、; 97        f        2.0 eVCf        98        f        2.0 eVEs        99        f        2.0 eVFm        100 

45、60;     f        2.0 eVMd        101        f        2.0 eVNo        102        f        2.0 eV1)3d轨道的带宽很窄即局域性很强，而对于这种局域性很强的带其就与其他能带的交叠较少，从而体现出的库伦作用很强

46、从，所以此时必须进行U的修正才能体现出强的关联作用来。至于2p轨道现在也有文献说其有较强的关联作用，但是这种体现在绝缘体和宽带隙半导体中，对于导体和带隙较小的效果不明显。至于U的引入其实就是使我们理论要和实验符合这样我们的结果能更好的解释于实验，无论预测还是对实验的支持都需要我们对各种方面考虑的很周全，所以都算一算看看哪个有个好的效果。2）至于d轨道加U的影响我没有做过，但是同组人有做过的，一般是这样d轨道U值的引入后我们都是固定费米面对比d带的变化，一般的情况都是p提供的价带，d提供导带，深层价带由d提供，所以DOS图中会发现，p带的宽和位置没有变化，而是深层d带明显下移。dx2-y2轨道对

47、应2dxz轨道对应1dz2轨道对应0dyz轨道对应-1dxy轨道对应-2最后提示一下大多数参数的首字母代表了参数的性质I 初始化 L 逻辑开关 E 能量 N 数目 T 温度便于记忆spin polarize是指在这个计算过程中要考虑到电子的自旋作用，因为加上U参数本身就是对于能量的一种修正，这个修正项就来自于电子之间由于自旋和库仑力产生的相互排斥作用。U定义为两个电子放在同一位置时花费的库仑作用能可以通过这样的方法计算：两个同样价态n的相同离子，一个离子失去一个电子n+1，另一个离子得到一个电子n-1。然后根据两者的能量差，(E(n)+E(n)-(E(n+1)+E(n-1)，可以得到U值。不过

48、算出的结果不是太准确。还会有大约0.02Ry的能量差。总结：U值很难精确的确定。而且相同原子在不同体系中U值不同。现在一般用经验值。【原创】MS电荷密度在Origin中作图方法1、Total density map现在MS中save为文件，当然你需要将你的隐藏文件和扩展名都能显示，在结果文件夹中找到*.chrgm(好像是).charg_frm文件从origin中导入这个文件,得到4列数据,分别为a, b, c, 和charge.一般的排序是先将a=1,对应的b,c不同,然后才是a=2的,以此类推.如果你想做(100)面, 将b, c, charge列分别设为x, y, z轴，那么就只选中a=1

49、对应的所有的b, c, charge数据，然后直接可以做出(100)面得电荷密度等高线图，通过调节，可以得到很好的图片。（？怎么插图片啊？）如果想得到（010）面的，将原来的数据的b列推至首列，然后选中首列，右键单击选择sort worksheet，这样就会对b列的数据按从小到大的顺序排列，而与之对应的a，c，charge随之改变。按同样的方法可得到（010）面得电荷密度。2、差分电荷密度MS中只能给出显示，存为文件后却无法打开。昨天在一个会上问M C Payne这个问题，他解释说不是他的错，是公司的问题。（具一个牛人说，可用命令行运行后，可得到数据文件。还没试）。如果就现在的差分密度图片如何

50、在Origin中做出来呢？你需要做的也很简单。将MS的差分密度图片显示后，隐藏其它（包括原子），存为.bmp文件。然后在origin中新建一个空matrix表格。从File的import中选择image，按照默认方式导入刚才的.bmp文件。然后选择你需要的区域，点击右键形成新的Matrix表（其实现实的还是图片）。再到view中的view mode选择数据，你会得到一个很奇怪的matrix表，别管他，直接对这个matrix做contour图，你会得到差分密度图了，要想获得好的效果，需要很仔细地调节。具体办法，实在太难写出来，等以后再说吧。本想贴几张我的结果，但不知如何贴。请斑竹给我指示！好，知

51、道了。上面一张是根据下面MS给出的差分密度map用Origin做出来的。由于为了满足附件的大小限制，对origin的输出做了压缩和色度限制，这只是今天初略做的，没有做详细的调节。其实关键在于，国外期刊一般倾向你提交.eps格式的图片文件，如果你将.bmp通过画图来转为.jpg，然后再用jpeg2ps程序转为.eps，效果实在不敢恭维。所以利用origin直接输出.eps格式的图片，效果很好，而且不失真，属于矢量图。 最近在做DFT+U计算，有些收获分享一下。本人数学功底先天不足，难免有理解错误的地方，欢迎指正。首先是对+U必要性的感悟。DFT（LDA和GGA）对一般体系的结构、能带计算还是很令

52、人满意的。这些一般体系是指前3周期的体系和纯金属体系。但是，对于包含d电子和f电子的体系，特别是过渡金属氧化物或氮化物，在“金属/绝缘体”的定性判断上常常出错。LDA和GGA往往把一些绝缘体/半导体的gap算的太小，甚至最高占据轨道/能带(HOMO)在Fermi面之上，即成金属了。这些最高占据轨道往往是来自这些金属原子的d电子和f电子。那么为什么DFT处理d/f电子会出现较大误差呢？我认为这些d/f电子的特殊性在于它们不仅离原子核较远，轨道空间伸展形状奇特（比如穿透效应，即穿透更加外层的s，p轨道）而且多数会发生自旋极化，即原应自旋相反的电子对发生自旋跃迁，导致金属原子神奇地出现了净的磁矩。如

53、果单胞内所有金属原子的净磁矩方向相同，则单胞总体产生净磁矩，材料呈现铁磁性；如果单胞内相邻金属原子的净磁矩相反，则单胞没有总体净磁矩，但存在一上一下的磁矩排列，材料呈现为反铁磁性。正是d/f电子这些神奇的特性，使得DFT计算的过渡金属氧化物或氮化物能带结构频频失误。为什么呢？原因在于DFT计算电子-电子之间的排斥作用所用的交换-相关泛函，是基于单粒子近似发展起来的。单粒子近似是把自旋配对的电子对看成单个粒子，这样Kohn-Sham轨道数量为总电子数的一半。但是当处理d/f电子时，不得不把电子对分开为自旋向上和自旋向下的两个粒子，即开壳层计算open shell。对开壳层体系使用单粒子近似下的交

54、换-相关泛函时，其中的相关泛函就不能完全考虑d/f电子的相关效应。d/f电子的相关效应，就是自旋电子与自旋电子之间的排斥，和自旋电子与内层电子之间的排斥。这些排斥作用使得轨道/能带之间分割较远，轨道/能带较窄。但是DFT的相关泛函对这些排斥考虑的不够，结果轨道/能带过宽，轨道与轨道相互接近甚至重叠。这就是为什么DFT把绝缘体计算成了金属。一般把这些过度金属半导体/绝缘体体系称为强相关体系strongly correlated system。怎么办呢？可能你已经直觉到了，直接加上一项，把d/f轨道上的自旋电子间的排斥作用作为（u/2*int(n_up*n_down)）一个能量项加进去。这个就是H

55、ubbard 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】。你可能会困惑，这有点像

56、GGA？是的，不过GGA是通用的交换-相关能对密度的二阶展开的泛函。而这里的U是纯DFT基态总能量对特定体系d/f电子密度变化的二阶导。这有点像是个二次方项的收缩函数，即你们这些自旋作用的电子在轨道上的占据必须往特定方向收缩，否则会有能量惩罚。 上图中，不加U时左边的最高占据轨道和右边的空轨道之间的gap很窄，是不正确的；+U后变宽，更加接近实验结果。需要注意的是，+U是对电子结构的改善，对原子几何、晶体结构的影响可能很有限。此时只要不涉及电子性质，一般的GGA, LDA的结果是可以接受的。那么这个+U能量惩罚项是人为的经验项，还是真实的物理存在呢？Cococcioni曾在【1】PRB71(h

57、ttp://pdf/PRB/v71/i3/e035105)中试图用图1来说明Hubbard U能量项的本质。他说真实的体系与环境交换电子，其能量变化应该是交换前后的两个状态的线性组合，这样的能量变化应该是线性的。但LDA和GGA对部分（非整数的）占据的Kohn-Sham轨道的自作用的不正确处理导致能量变化是非线性的。这二者之间的能量差就是Hubbard U能量项的来源。从这个解释可以看出，(1) +U能量惩罚项是物理存在的，并非经验校正。参数U的确定也不应该是经验的，而是必须考察DFT与真实势能面的差。(2) 这个二级导是体系环境变化而变化的。这就导致必须对每个体系，甚至

58、体系的不同状态（比如键长变化）计算相应的参数U。注意如果体系内有对称性不相同的金属位点，需要对每个位点上的金属原子计算其相应的参数U，不管它们是否为同种金属！获得合适的参数U以后，才能使用DFT+U模型精确确定自旋极化和能带、带隙等问题。（当然，更精确的模型还要考虑旋轨耦合spin-orbital coupling。）怎么确定某个体系参数U呢？因为U是纯DFT基态总能量对这个体系d/f电子密度变化的二阶导，那么很自然的一个想法，我们对密度加以扰动，求能量的变化（响应），可以得到数值二阶导（微扰法）。不过，密度并非一个理想的微扰项，因为一个金属原子（Hubbard位点）涉及到多个电子。怎么办呢？

59、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，数据从自洽前第

60、一个循环计算得到的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-polar

61、ized and DFT+U calculations注意上面这种方法基于一个假定，即势alpha足够小，此时得到的U近似等于alpha=0时的U，即U_0。这是线性响应理论。即体系能量实际对势alpha的响应并不是线性的，但我们假定扰动足够小的情况下假定这种响应是线性的。上述方法对单个Hubard位点施加势alpha，但是在周期性体系中，任何扰动也都是周期性的。怎么办呢？Cococcioni提出的方法是延伸单胞为超胞，做响应计算；然后再次延伸为更大超胞，做响应计算。这样随着盒子越来越大，势在其中的作用越来越小，直到盒子无限大，则响应系数收敛。这种方法并不需要计算实际超胞，而是对密度响应进行外

62、推。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，然后

63、将上述(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。其中所用的例子