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

资料来自小木虫论坛；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,%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，SPi