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,'r');head=fscanf(fp,&

2、#39;%s',4);data=fscanf(fp,'%f',4,18*18*144);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;c

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

4、(d(:,layer,:),18,144),8);>> hold on>> layer=1;contour(reshape(a(:,layer,:),18,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,14

5、4),reshape(c(layer,l:,:),18,144),reshape(d(layer,:,:),18,144),8);MS 4.3磁性体系LDAU计算设置在最新发布的MS版本中增加了LDAU的计算，这是MS计算性能特别是CASTEP模块的一次巨大提高，终于解决了已久的d，f轨道带隙分裂问题，也不会再被审稿人轻视认为是简单DFT计算结果了。这里讲解一下MS.43里面element设置方面的问题，我们选择了Fe3O4体系，该体系是铁磁性的，计算该体系主要涉及到两个参数设置，一个是自旋设置（SPin parameters），其次是Fe Hubbard参数设置。首先介绍需要用到的菜单：打开

6、MS之后看到主菜单，打开如下图所示的菜单：（点击此图，可放大观看） 由于计算体系位Fe3O4，计算采用USPP完成，Fe的USPP价电子参数为3d64s2，O是2s22p6，结构中的磁性主要来自于Fe原子磁矩的定向排列，计算中自旋参数和Hubbard参数主要是针对Fe调整。再建造好Fe3O4结构之后（注MS结构库中包含此结构），为一个面心立方晶格，一个惯用晶胞包含8个Fe3O4单位，因此整个结构有24个Fe原子，32个O。计算可以采用FCC结构的初级晶胞结构，也可以采用惯用晶胞。我们计算采用后者。打开Modifyelectronic configuration菜单，看到Spin和Hubbard

7、 U两个子菜单如图：SPin主要是自旋控制参数，由于再铁磁性结构中Fe自旋是高度极化的，因此这个地方设置为High spin state，即3d6结构为4个自旋平行同向电子。自旋方向不管是Up还是Down当然不会改变Fe原子自旋极化的性质，但会改变总的磁矩（自旋电子数目，这点会再CASTEP计算参数Spin设置中体现出来，可以尝试把部分Fe原子自旋改成High，SPin up，部分为High，Down，然后查看CASTEP计算控制面板的SPin数值）。 下面按照Fe3O4具体情况设置，Fe自旋为High spin， 方向Up。自旋电子数目4，O原子不予考虑。HUbbard带参数设置，主要针对d

8、轨道，Fe的Hubbard参数再CASTEP里面默认数值是2.5eV，根据一些文献也可以选取其他参数，比如再FeO计算文献中选用了6eV，也有用45eV的，根据具体计算情况和结果可以调整，这可能会要求计算几次，比较结果。但为了保证结果准确性，是必要的。我们计算选Fe的d U5eV。如下所示：元素的电荷这里不需要考虑，实际上一般QM计算原子电荷状态是中性原子状态，再MD计算中，可能对这个参数有要求，或者计算原子必须要带电荷。一般默认即可！ 接下来进入MS CASTEP计算控制菜单：MS4.3中CASTEP模块功能增加了泛函类型，如增加了HF交换，HF交换LDA关联等，增加了LDAU计算。在上面的

9、计算中对Fe3O4Fe的参数做了修正，在CASTEP计算控制面板中可以看到默认的自旋电子数目是96个，因为我们采用了惯用晶胞，24个Fe原子，每个自旋电子数目是4个，因此是96。计算中无需再次对此进行设置，默认即可。Hubbard计算参数已经设置，此时只要选中该计算项目即可，需要注意的是LDAU计算只支持自旋极化计算模型(spin polarized scheme )，并且不能用于应力，NMR，声子谱等的计算，只能用于能量计算任务，而不是结构优化任务。最后设置如下图： LDA+U是对平均场计算结果的一个微扰，只不过是对特定的轨道做了一个Coulomb相互作用的修正（Onsite），因此这种微扰

10、不会对LDA计算得到的几何结构有太大的改变，只是对电子结构有影响，而且特别是那些d或者f轨道能级分布。因此采用LDA计算结构，得到优化结构，然后引入U对d轨道能级的修正应该来说是不会对平衡结构有太大的影响Hubbard 参数U和Stoner 自旋极化参量在平均场近似或者说是一般得LDA计算中，能带得自旋分裂是由一个称为Stoner 参数I来主导得，而且平均场近似中认为这个交换分裂能是小于带宽得，一般而言I得数值在LSDA中大概在1eV左右，这样即使对于某些氧化物如NiO，FeO，MnO等即使能带发生交换分裂，整个体系依然是金属性得，试验观察到得这些氧化物实际上是绝缘体，主要原因在于在这些氧化物

11、中d轨道得能级位置不是由I来决定得，而是Hubbard参数U决定，U也称为Onsite Coulomb作用能，相当于把两个电子放在空间同一个位置需要得能量，U数值一般在10eV左右，如此大得分裂能足以将Fermi面附近连续分布得d能级分开，从而得到正确得基态性质，目前广泛采用得LDAU的算法就是针对某些定域轨道，如d或者f，这些轨道占据存在强烈的在位Coulomb排斥。正是U使得带隙分裂，而不是轨道极化参量I。 U计算目前仍然是一个研究热点，最近由人采用线性响应理论，同时结合轨道束缚的DFT计算自洽的求解了部分体系的U参数，在大部分情况下U数值要根据计算结果和试验参数的符合程度而定。 LDAU

12、计算核心思路是：首先将研究体系的轨道分隔成两个子体系（subsystem），其中一部分是一般的DFT算法（如LSDA，GGA）等可以比较准确描述的体系，另外是定域在原子周围的轨道如d或者f轨道，这些轨道在标准的DFT计算下不能获得正确的能量与占据数之间的关系（如DFT总是认为分数占据是能量最小的，而不是整数占据）；对于d或者f轨道，能带模型采用Hubbard模型，而其他轨道仍然是按照KohnSHam方程求解；d以及f轨道电子之间的关联能采用一个和轨道占据以及自旋相关的有效U表示；整体计算的时候需要将原来DFT计算过程中已经包含的部分关联能扣除，这部分一般叫Double Counting par

13、t，并且用一个新的U来表示，最终的结果是在DFT计算的基础上新增加一个和d或者f轨道直接相关的分裂势的微扰项，这部分能量可以采用一般微扰理论计算。 在CASTEP最新的版本中增加的LDAU的计算，U参数的设置一般主要是针对过渡金属氧化物（Charge transfer type insulator），包含非满层f轨道的元素等，高温超导体强关联体系。在参数设置方面主要是需要注意d和f轨道，至于s以及p轨道一般不需要设置，当然由文献也报道p轨道的这种关联性。过渡金属氧化物的有效U如下：3d轨道U和J计算如下所示：主要原理是改变d轨道的占据，在自旋极化的前提下计算不同自旋轨道能量的差值提取U和J，U

14、微Coulomb排斥能，J是交换能，U在所有电子中都存在，不管自旋是否相同，J只存在于自旋相同的电子上。 下面给出过渡金属（不包括稀土元素）U和J参数的选取(uint in Ry,1Ry=13.6eV)：对于其他过渡金属化合物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 limit in the local density approxim