CASTEP獲得電荷密度等高線的Matlab作圖法

1、資料來自小木蟲論壇;H.J.Zhao匯總CASTEP獲得電荷密度等高線的Matlab作圖法MATLAB程序用來處理MS電荷密度的等高線做法，分四步：1找到電荷密度圖的存儲文件，一般是在自己所建的計算project中，隱藏的，castep模塊是.charg_frm格式的，可以用txt打開后另存為txt格式。(這個地方注意把文件單獨拷貝出來，以免破壞原始文件)2讀取數(shù)據(jù)，在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數(shù)據(jù)的，其中18*18*144是可以改的，分別對應(yīng)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是可以改的，分別對應(yīng)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.)可以修改作圖的層數(shù)，分別是從前面到后面。最后的8表示顯示線條數(shù)目，可以自己修改成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);這是沿側(cè)面做等高線(即xz平面)，其它方向作圖類似。只要修改其中的參數(shù)就可以了，如把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計算設(shè)置在最新發(fā)布的MS版本中增加了LDA+U的計算，這是MS計算性能特別是CASTEP模塊的一次巨大提高，終于解決了已久的d,f軌道帶隙分裂問題，也不會再被審稿人輕視認為是簡單DFT計算結(jié)果了。這里講解一下MS.43里

6、面element設(shè)置方面的問題，我們選擇了Fe3O4體系，該體系是鐵磁性的，計算該體系主要涉及到兩個參數(shù)設(shè)置，一個是自旋設(shè)置（SPinparameters），其次是FeHubbard參數(shù)設(shè)置。首先介紹需要用到的菜單：打開MS之后看到主菜單，打開如下圖所示的菜單：（點擊此圖，可放大觀看）由于計算體系位Fe3O4,計算采用USPP完成,Fe的USPP價電子參數(shù)為3d64s2,O是2s22p6,結(jié)構(gòu)中的磁性主要來自于Fe原子磁矩的定向排列，計算中自旋參數(shù)和Hubbard參數(shù)主要是針對Fe調(diào)整。再建造好Fe3O4結(jié)構(gòu)之后（注MS結(jié)構(gòu)庫中包含此結(jié)構(gòu)），為一個面心立方晶格，一個慣用晶胞包含8個Fe3O4單

7、位,因此整個結(jié)構(gòu)有24個Fe原子，32個O。計算可以采用FCC結(jié)構(gòu)的初級晶胞結(jié)構(gòu)，也可以采用慣用晶胞。我們計算采用后者。打開Modifyelectronicconfiguration菜單，看到Spin和HubbardU兩個子菜單如圖：SPin主要是自旋控制參數(shù)，由于再鐵磁性結(jié)構(gòu)中Fe自旋是高度極化的，因此這個地方設(shè)置為Highspinstate,即3d6結(jié)構(gòu)為4個自旋平行同向電子。自旋方向不管是Up還是Down當然不會改變Fe原子自旋極化的性質(zhì)，但會改變總的磁矩（自旋電子數(shù)目，這點會再CASTEP計算參數(shù)Spin設(shè)置中體現(xiàn)出來，可以嘗試把部分Fe原子自旋改成High,SPinup,部分為Hig

8、h,Down,然后查看CASTEP計算控制面板的SPin數(shù)值）。下面按照Fe3O4具體情況設(shè)置，F(xiàn)e自旋為Highspin,方向Up。自旋電子數(shù)目+4,O原子不予考慮。HUbbard帶參數(shù)設(shè)置，主要針對d軌道，F(xiàn)e的Hubbard參數(shù)再CASTEP里面默認數(shù)值是2.5eV，根據(jù)一些文獻也可以選取其他參數(shù)，比如再FeO計算文獻中選用了6eV,也有用45eV的，根據(jù)具體計算情況和結(jié)果可以調(diào)整，這可能會要求計算幾次，比較結(jié)果。但為了保證結(jié)果準確性，是必要的。我們計算選Fe的dU=5eV。如下所示：元素的電荷這里不需要考慮，實際上一般QM計算原子電荷狀態(tài)是中性原子狀態(tài)，再MD計算中，可能對這個參數(shù)有要

9、求，或者計算原子必須要帶電荷。一般默認即可！接下來進入MSCASTEP計算控制菜單:MS4.3中CASTEP模塊功能增加了泛函類型，如增加了HF交換，HF交換+LDA關(guān)聯(lián)等，增加了LDA+U計算。在上面的計算中對Fe3O4Fe的參數(shù)做了修正，在CASTEP計算控制面板中可以看到默認的自旋電子數(shù)目是96個，因為我們采用了慣用晶胞，24個Fe原子，每個自旋電子數(shù)目是4個，因此是96。計算中無需再次對此進行設(shè)置，默認即可。Hubbard計算參數(shù)已經(jīng)設(shè)置，此時只要選中該計算項目即可，需要注意的是LDA+U計算只支持自旋極化計算模型(spinpolarizedscheme)，并且不能用于應(yīng)力，NMR，聲

10、子譜等的計算，只能用于能量計算任務(wù)，而不是結(jié)構(gòu)優(yōu)化任務(wù)。最后設(shè)置如下圖：LDA+U是對平均場計算結(jié)果的一個微擾，只不過是對特定的軌道做了一個Coulomb相互作用的修正（Onsite）,因此這種微擾不會對LDA計算得到的幾何結(jié)構(gòu)有太大的改變，只是對電子結(jié)構(gòu)有影響，而且特別是那些d或者f軌道能級分布。因此采用LDA計算結(jié)構(gòu)，得到優(yōu)化結(jié)構(gòu)，然后引入U對d軌道能級的修正應(yīng)該來說是不會對平衡結(jié)構(gòu)有太大的影響Hubbard參數(shù)U和Stoner自旋極化參量在平均場近似或者說是一般得LDA計算中，能帶得自旋分裂是由一個稱為Stoner參數(shù)I來主導得，而且平均場近似中認為這個交換分裂能是小于帶寬得,一般而言I

11、得數(shù)值在LSDA中大概在1eV左右,這樣即使對于某些氧化物如NiO,FeO,MnO等即使能帶發(fā)生交換分裂，整個體系依然是金屬性得，試驗觀察到得這些氧化物實際上是絕緣體，主要原因在于在這些氧化物中d軌道得能級位置不是由I來決定得，而是Hubbard參數(shù)U決定，U也稱為OnsiteCoulomb作用能，相當于把兩個電子放在空間同一個位置需要得能量，U數(shù)值一般在10eV左右，如此大得分裂能足以將Fermi面附近連續(xù)分布得d能級分開，從而得到正確得基態(tài)性質(zhì)，目前廣泛采用得LDA+U的算法就是針對某些定域軌道，如d或者f,這些軌道占據(jù)存在強烈的在位Coulomb排斥。正是U使得帶隙分裂，而不是軌道極化參

12、量I。U計算目前仍然是一個研究熱點，最近由人采用線性響應(yīng)理論，同時結(jié)合軌道束縛的DFT計算自洽的求解了部分體系的U參數(shù)，在大部分情況下U數(shù)值要根據(jù)計算結(jié)果和試驗參數(shù)的符合程度而定。LDA+U計算核心思路是：首先將研究體系的軌道分隔成兩個子體系(subsystem),其中一部分是一般的DFT算法(如LSDA,GGA)等可以比較準確描述的體系，另外是定域在原子周圍的軌道如d或者f軌道，這些軌道在標準的DFT計算下不能獲得正確的能量與占據(jù)數(shù)之間的關(guān)系(如DFT總是認為分數(shù)占據(jù)是能量最小的，而不是整數(shù)占據(jù))；對于d或者f軌道，能帶模型采用Hubbard模型，而其他軌道仍然是按照KohnSHam方程求解

13、；d以及f軌道電子之間的關(guān)聯(lián)能采用一個和軌道占據(jù)以及自旋相關(guān)的有效U表示；整體計算的時候需要將原來DFT計算過程中已經(jīng)包含的部分關(guān)聯(lián)能扣除，這部分一般叫DoubleCountingpart，并且用一個新的U來表示，最終的結(jié)果是在DFT計算的基礎(chǔ)上新增加一個和d或者f軌道直接相關(guān)的分裂勢的微擾項，這部分能量可以采用一般微擾理論計算。在CASTEP最新的版本中增加的LDA+U的計算，U參數(shù)的設(shè)置一般主要是針對過渡金屬氧化物(Chargetransfertypeinsulator)，包含非滿層f軌道的元素等，高溫超導體強關(guān)聯(lián)體系。在參數(shù)設(shè)置方面主要是需要注意d和f軌道，至于s以及p軌道一般不需要設(shè)置

14、，當然由文獻也報道p軌道的這種關(guān)聯(lián)性。過渡金屬氧化物的有效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軌道的占據(jù)，在自旋極化的前提下計算不同自旋軌道能量的差值提取U和J，U微Coulomb排斥能，J是交換能，U在所有電子中都存在，不管自旋是

15、否相同，J只存在于自旋相同的電子上。下面給出過渡金屬(不包括稀土元素)U和J參數(shù)的選取(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數(shù)值和d電子排列以及價態(tài)有關(guān)系，因此上面給出的數(shù)值只是一個大概的估算數(shù)值，具體文獻見PhysicalReviewBVol50,No23,1994.LDA+U算法主要原創(chuàng)作者是俄羅斯金屬研究所的V.Anisimov，重要文獻有：CorrectedatomiclimitinthelocaldensityapproximatenandtheelectronicstructureofdimpuritiesinRb,Phys

17、.Rev.B50,23(1994)；BandtheoryandMottinsulators:HubbardUinsteadofStonerI,Phys.Rev.B.44No.3(1991);MaterialsStudio4.3版本中也給出了一些元素默認的U數(shù)值(實際是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