CASTEP获得电荷密度等高线的Matlab作图法

1、资料来自小木虫论坛;H.J.Zhao汇总CASTEP获得电荷密度等高线的Matlab作图法MATLAB程序用来处理MS电荷密度的等高线做法，分四步：1找到电荷密度图的存储文件，一般是在自己所建的计算project中，隐藏的，castep模块是.charg_frm格式的，可以用txt打开后另存为txt格式。(这个地方注意把文件单独拷贝出来，以免破坏原始文件)2读取数据，在matlab命令里输入一下命令：path,fn=uigetfile(*txt,Open);fp=fopen(fn,path,l);head=fscanf(fp,%s,4);data=fscanf(fp,%f,4,18*18*14

2、4);fclose(fp);这是读取刚才保存的txt数据的，其中18*18*144是可以改的，分别对应a,b,c的重复单元。3作矩阵。a=reshape(data(1,:),18,18,144);b=reshape(data(2,:),18,18,144);c=reshape(data(3,:),18,18,144);d=reshape(data(4,:),18,18,144);其中18*18*144是可以改的，分别对应a,b,c的重复单元，即长宽高的值。4作等高线图layer=1;contour(reshape(a(:,layer,:),18,144),reshape(c(:,layer,:

3、),18,144),reshape(d(:,layer,:),18,144),8);layer=1(2，3，4.)可以修改作图的层数，分别是从前面到后面。最后的8表示显示线条数目，可以自己修改成9，10或其它。holdon表示累积作图，holdoff表示清空前面的作图举例；做表面层和第九层的等高线图layer=9;contour(reshape(a(:,layer,:),18,144),reshape(c(:,layer,:),18,144),reshape(d(:,layer,:),18,144),8);holdonlayer=1;contour(reshape(a(:,layer,:),1

4、8,144),reshape(c(:,layer,:),18,144),reshape(d(:,layer,:),18,144),8);这是沿侧面做等高线(即xz平面)，其它方向作图类似。只要修改其中的参数就可以了，如把a(:,layer,:),18,144),reshape(c(:,layer,:),18,144),reshape(d(:,layer,:),18,144),8);改为a(layer,:,:),18,144),reshape(c(layer,l:,:),18,144),reshape(d(layer,:,:),18,144),8);1.00.0-0.50.55-1o-D5000

5、5101502.0-1.0-1.5-2：.O2.01.51.00.50.0QJfl匚負An-ixx-O.7-5OO-Oj6s250-O.5DOO0.3750ZPI00.125002SM1.00.80.60.40.20.0-0.2-0.4-0.6-0.8-1.0-1.0-0.8-0.6-0.4-020.00.20.40.60.81.0AxisTitleMS4.3磁性体系LDA+U计算设置在最新发布的MS版本中增加了LDA+U的计算，这是MS计算性能特别是CASTEP模块的一次巨大提高，终于解决了已久的d,f轨道带隙分裂问题，也不会再被审稿人轻视认为是简单DFT计算结果了。这里讲解一下MS.43里

6、面element设置方面的问题，我们选择了Fe3O4体系，该体系是铁磁性的，计算该体系主要涉及到两个参数设置，一个是自旋设置（SPinparameters），其次是FeHubbard参数设置。首先介绍需要用到的菜单：打开MS之后看到主菜单，打开如下图所示的菜单：（点击此图，可放大观看）由于计算体系位Fe3O4,计算采用USPP完成,Fe的USPP价电子参数为3d64s2,O是2s22p6,结构中的磁性主要来自于Fe原子磁矩的定向排列，计算中自旋参数和Hubbard参数主要是针对Fe调整。再建造好Fe3O4结构之后（注MS结构库中包含此结构），为一个面心立方晶格，一个惯用晶胞包含8个Fe3O4单

7、位,因此整个结构有24个Fe原子，32个O。计算可以采用FCC结构的初级晶胞结构，也可以采用惯用晶胞。我们计算采用后者。打开Modifyelectronicconfiguration菜单，看到Spin和HubbardU两个子菜单如图：SPin主要是自旋控制参数，由于再铁磁性结构中Fe自旋是高度极化的，因此这个地方设置为Highspinstate,即3d6结构为4个自旋平行同向电子。自旋方向不管是Up还是Down当然不会改变Fe原子自旋极化的性质，但会改变总的磁矩（自旋电子数目，这点会再CASTEP计算参数Spin设置中体现出来，可以尝试把部分Fe原子自旋改成High,SPinup,部分为Hig

8、h,Down,然后查看CASTEP计算控制面板的SPin数值）。下面按照Fe3O4具体情况设置，Fe自旋为Highspin,方向Up。自旋电子数目+4,O原子不予考虑。HUbbard带参数设置，主要针对d轨道，Fe的Hubbard参数再CASTEP里面默认数值是2.5eV，根据一些文献也可以选取其他参数，比如再FeO计算文献中选用了6eV,也有用45eV的，根据具体计算情况和结果可以调整，这可能会要求计算几次，比较结果。但为了保证结果准确性，是必要的。我们计算选Fe的dU=5eV。如下所示：元素的电荷这里不需要考虑，实际上一般QM计算原子电荷状态是中性原子状态，再MD计算中，可能对这个参数有要

9、求，或者计算原子必须要带电荷。一般默认即可！接下来进入MSCASTEP计算控制菜单:MS4.3中CASTEP模块功能增加了泛函类型，如增加了HF交换，HF交换+LDA关联等，增加了LDA+U计算。在上面的计算中对Fe3O4Fe的参数做了修正，在CASTEP计算控制面板中可以看到默认的自旋电子数目是96个，因为我们采用了惯用晶胞，24个Fe原子，每个自旋电子数目是4个，因此是96。计算中无需再次对此进行设置，默认即可。Hubbard计算参数已经设置，此时只要选中该计算项目即可，需要注意的是LDA+U计算只支持自旋极化计算模型(spinpolarizedscheme)，并且不能用于应力，NMR，声

10、子谱等的计算，只能用于能量计算任务，而不是结构优化任务。最后设置如下图：LDA+U是对平均场计算结果的一个微扰，只不过是对特定的轨道做了一个Coulomb相互作用的修正（Onsite）,因此这种微扰不会对LDA计算得到的几何结构有太大的改变，只是对电子结构有影响，而且特别是那些d或者f轨道能级分布。因此采用LDA计算结构，得到优化结构，然后引入U对d轨道能级的修正应该来说是不会对平衡结构有太大的影响Hubbard参数U和Stoner自旋极化参量在平均场近似或者说是一般得LDA计算中，能带得自旋分裂是由一个称为Stoner参数I来主导得，而且平均场近似中认为这个交换分裂能是小于带宽得,一般而言I

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

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

13、；d以及f轨道电子之间的关联能采用一个和轨道占据以及自旋相关的有效U表示；整体计算的时候需要将原来DFT计算过程中已经包含的部分关联能扣除，这部分一般叫DoubleCountingpart，并且用一个新的U来表示，最终的结果是在DFT计算的基础上新增加一个和d或者f轨道直接相关的分裂势的微扰项，这部分能量可以采用一般微扰理论计算。在CASTEP最新的版本中增加的LDA+U的计算，U参数的设置一般主要是针对过渡金属氧化物(Chargetransfertypeinsulator)，包含非满层f轨道的元素等，高温超导体强关联体系。在参数设置方面主要是需要注意d和f轨道，至于s以及p轨道一般不需要设置

14、，当然由文献也报道p轨道的这种关联性。过渡金属氧化物的有效U如下：SpeciesUJU-J(UEfDNiO3.00.957.1CoO7.30.926.9FeO6.S0.S95.9MnO6.90.8610.3VO6.70.815.9TiO6.60.735.3Reference:BandtheoryandMottinsulators:HubbardUinsteatofStonerI,PRBVol44No3(1991);3d轨道U和J计算如下所示：主要原理是改变d轨道的占据，在自旋极化的前提下计算不同自旋轨道能量的差值提取U和J，U微Coulomb排斥能，J是交换能，U在所有电子中都存在，不管自旋是

15、否相同，J只存在于自旋相同的电子上。下面给出过渡金属(不包括稀土元素)U和J参数的选取(uintinRy,1Ry=13.6eV)：3dseries:ElementsUJV0.250.05Or0.260.053Mn0.2S0.055Fe0.30.050Co0.310.059Ni0.310.064dseries:ElementsUJNb0.190.04Mo0.20.04Tc0.210.042Ru0.220.042Rh0.250.044Pd0.290.0445dseries:ElementsUJTa0.190.039W0.200.038Re0.2050.039Os0.20.039Ir0.210.03

16、3Pt0.2150.038对于其他过渡金属化合物U般在510eV之间。如在PRB73,134418(2006)这个文献中作者在计算Co掺杂的ZnO时采用的U是6和8eV。过渡金属的U数值和d电子排列以及价态有关系，因此上面给出的数值只是一个大概的估算数值，具体文献见PhysicalReviewBVol50,No23,1994.LDA+U算法主要原创作者是俄罗斯金属研究所的V.Anisimov，重要文献有：CorrectedatomiclimitinthelocaldensityapproximatenandtheelectronicstructureofdimpuritiesinRb,Phys

17、.Rev.B50,23(1994)；BandtheoryandMottinsulators:HubbardUinsteadofStonerI,Phys.Rev.B.44No.3(1991);MaterialsStudio4.3版本中也给出了一些元素默认的U数值(实际是Ueff=UJ）：ElementNameAtomicnumberAngularMomentumHubbardUSc21d2.5eVTi22d2.5eVV23d2.5eVCr24d2.5贺Mn25d2.5eVFe26d2.5eVCo27d2.5eVNi2Sd2.5eVCu29d2.5eVY39d2.0eVZr40d2.0贺Nb41d2.0eVMo42d2.0eVTc43d2.0eVRu44d2.0eVRh45d2.0eVPd46d2.0eVAg47d2.0贺Cd48d2.0eVLa57f6.0贺Ce58f6.0eVPr59f6.0