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

资料来自小木虫论坛；汇总

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,'%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;

contour(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(

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,144]),reshape(c(layer,l:,:),[18,144]),reshape(d(layer,:,:),[18,144]),8);

MS 4.3磁性体系LDA＋U计算设置

在最新发布的MS版本中增加了LDA＋U的计算，这是MS计算性能特别是CASTEP模块的一次巨大提高，终于解决了已久的d，f轨道带隙分裂问题，也不会再被审稿人轻视认为是简单DFT计算结果了。

这里讲解一下MS.43里面element设置方面的问题，我们选择了Fe3O4体系，该体系是铁磁性的，计算该体系主要涉及到两个参数设置，一个是自旋设置（SPin

parameters），其次是Fe Hubbard参数设置。

首先介绍需要用到的菜单：

打开MS之后看到主菜单，打开如下图所示的菜单：（点击此图，可放大观看）

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

打开Modify－electronic configuration菜单，看到Spin和Hubbard 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轨道，Fe的Hubbard参数再CASTEP里面默认数值是2.5eV，根据一些文献也可以选取其他参数，比如再FeO计算文献中选用了6eV，也有用4－5eV的，根据具体计算情况和结果可以调整，这可能会要求计算几次，比较结果。但为了保证结果准确性，是必要的。我们计算选Fe的d U＝5eV。如下所示：

元素的电荷这里不需要考虑，实际上一般QM计算原子电荷状态是中性原子状态，再MD计算中，可能对这个参数有要求，或者计算原子必须要带电荷。一般默认即可！

接下来进入MS CASTEP计算控制菜单：

MS4.3中CASTEP模块功能增加了泛函类型，如增加了HF交换，HF交换＋LDA关联等，增加了LDA＋U计算。在上面的计算中对Fe3O4Fe的参数做了修正，在CASTEP计算控制面板中可以看到默认的自旋电子数目是96个，因为我们采用了惯用晶胞，24个Fe原子，每个自旋电子数目是4个，因此是96。计算中无需再次对此进行设置，默认即可。Hubbard计算参数已经设置，此时只要选中该计算项目即可，需要注意的是LDA＋U计算只支持自旋极化计算模型(spin

polarized scheme )，并且不能用于应力，NMR，声子谱等的计算，只能用于能量计算任务，而不是结构优化任务。最后设置如下图：

LDA+U是对平均场计算结果的一个微扰，只不过是对特定的轨道做了一个Coulomb相互作用的修正（On－site），因此这种微扰不会对LDA计算得到的几何结构有太大的改变，只是对电子结构有影响，而且特别是那些d或者f轨道能级分布。因此采用LDA计算结构，得到优化结构，然后引入U对d轨道能级的修正应该来说是不会对平衡结构有太大的影响

Hubbard 参数U和Stoner 自旋极化参量

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

U计算目前仍然是一个研究热点，最近由人采用线性响应理论，同时结合轨

道束缚的DFT计算自洽的求解了部分体系的U参数，在大部分情况下U数值要根据计算结果和试验参数的符合程度而定。

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

在CASTEP最新的版本中增加的LDA＋U的计算，U参数的设置一般主要是针对过渡金属氧化物（Charge transfer type insulator），包含非满层f轨道的元素等，高温超导体强关联体系。在参数设置方面主要是需要注意d和f轨道，至于s以及p轨道一般不需要设置，当然由文献也报道p轨道的这种关联性。过渡金属氧化物的有效U如下：

3d轨道U和J计算如下所示：主要原理是改变d轨道的占据，在自旋极化的前提下计算不同自旋轨道能量的差值提取U和J，U微Coulomb排斥能，J是交换能，U在所有电子中都存在，不管自旋是否相同，J只存在于自旋相同的电子上。

下面给出过渡金属（不包括稀土元素）U和J参数的选取(uint in

Ry,1Ry=13.6eV)：

对于其他过渡金属化合物U一般在5－10eV之间。如在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 approximation and the electronic

structure of d impurities in Rb, .B 50,23 (1994）；

Band theory and Mott insulators: Hubbard U instead of Stoner I, .B.

44 No.3 (1991);

Materials Studio 4.3版本中也给出了一些元素默认的U数值（实际是Ueff＝U－J）：