ABINIT實(shí)戰(zhàn)手冊(cè).doc

ABINIT-5.8.4軟件的Windows版本實(shí)戰(zhàn)守則（by H.J.Zhao）前言 第一性原理計(jì)算發(fā)表文章的要點(diǎn)1. 你算的是什么，2. 你為什么要算這個(gè)東西，3. 你用了怎樣的算法，4. 你算出來了什么，5. 你怎樣解釋你算出來的東西，有什么結(jié)論。第一章 ABINIT軟件的運(yùn)行/user1/11542/archives/2005/363031.shtml以BaTiO3為例，需建立輸入文件*.in，在abinit網(wǎng)站下載Ba、Ti、O的贗勢(shì)文件，之后執(zhí)行abinis程序，依次輸入：1. bto.in2. bto.out3. btoi4. btoo5. bto.temp6. 把三個(gè)贗勢(shì)文件拖入程序執(zhí)行框中運(yùn)行即可！Abinit使用篇簡(jiǎn)介abinis的輸入文件分為三類：定義了所有輸入文件名的文件（比如in.files，這個(gè)文件就是用來告訴abinit哪些文件是abinit讀入?yún)?shù)的文件名，以及輸出主要結(jié)果的文件的名稱），定義了計(jì)算的控制參數(shù)的文件（比如取名為INP)，贗勢(shì)的文件。下面舉例如何做好輸入文件計(jì)算金剛石結(jié)構(gòu)的Si的狀態(tài)方程，由此也得到Si的晶格常數(shù)和體彈性模量。大致的步驟是在Si的晶格常數(shù)的實(shí)驗(yàn)值附近取11個(gè)數(shù)據(jù)點(diǎn)，也就是說取11個(gè)晶格常數(shù)或體積，然后計(jì)算在這些晶格常數(shù)下的總能。在計(jì)算得到總能后，采用狀態(tài)方程擬合得到狀態(tài)方程、平衡態(tài)時(shí)的體積（或晶格常數(shù)）和體彈性模量。本例子中采用的是LDA-HGH贗勢(shì)。贗勢(shì)的文件名為：14si.4.hgh。in.files的內(nèi)容為（紫色標(biāo)示）：#BeginINP #設(shè)置關(guān)鍵詞的文件名為INPOUT #主要的輸出文件為OUT，該文件將被寫入計(jì)算最重要的結(jié)果sii sisi14si.4.hgh #贗勢(shì)的文件名#ENDINP文件的內(nèi)容為：# Crystalline cubic Si #ndtset 11 #說明下面將有11組數(shù)據(jù)acell: 3*9.8112 # 晶格常數(shù)a=b=c，將從9.8112. a.u.開始增加acell+ 3*0.09 #晶格常數(shù)將以0.09 a.u.的間隔進(jìn)行增加#Ground state calculationkptopt 1 #在 k點(diǎn)網(wǎng)格取樣時(shí)根據(jù)對(duì)稱性來取樣，并由下面的 # ngkpt和kptrlatt, 或者nshiftk和shiftk來確定k點(diǎn)的數(shù)目iscf 5 #采用CG方法對(duì)能量進(jìn)行優(yōu)化，用在基態(tài)計(jì)算中。#Definition of the unit cell rprim 0.0 0.5 0.5 #下面三行定義了原胞的基矢，本例子中Si是fcc結(jié)構(gòu) 0.5 0.0 0.5 0.5 0.5 0.0 #Definition of the atom typesntypat 1 #定義原胞中原子的類別的數(shù)目，本例子中只有1類原子znucl 14 #定義原胞中原子的核電荷數(shù)#Definition of the atomsnatom 2 #定義原胞中原子的總個(gè)數(shù)，本例子中有2個(gè)原子typat 2*1 #定義每類原子的個(gè)數(shù)，本例子中第一類原子有2個(gè)xred #下面定義了原胞中原子的坐標(biāo) 0.0 0.0 0.0 0.25 0.25 0.25 #Gives the number of band, explicitely (do not take the default)nband 16 #定義了要計(jì)算的能帶的數(shù)目，最好按這樣來設(shè)置: # nband 原胞中總的價(jià)電子數(shù)目/2 + 10#Exchange-correlation functionalixc 1 #定義交換關(guān)聯(lián)函數(shù)，本例子中，采用的是Teter Pade 參數(shù)化的LDA形式#Definition of the planewave basis setecut 20.0 #定義了平面波的切斷動(dòng)能 #Definition of the k-point gridngkpt 8 8 8 #下面定義了k點(diǎn)網(wǎng)格取樣的大小nshiftk 4 shiftk 0.5 0.5 0.5 0.5 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.5#Definition of the SCF procedurenstep 60 #電子自洽迭代的最大步數(shù)diemac 12.0 #介電常數(shù)設(shè)置tolvrs 1.0d-20 #電子自洽收斂的標(biāo)準(zhǔn)- END-計(jì)算完后，得到OUT文件，用下面的命令：grep volume OUT得到如下的內(nèi)容：Unit cell volume ucvol= 2.3610688E+02 bohr3Unit cell volume ucvol= 2.4266422E+02 bohr3Unit cell volume ucvol= 2.4934185E+02 bohr3Unit cell volume ucvol= 2.5614088E+02 bohr3Unit cell volume ucvol= 2.6306239E+02 bohr3Unit cell volume ucvol= 2.7010748E+02 bohr3Unit cell volume ucvol= 2.7727725E+02 bohr3Unit cell volume ucvol= 2.8457279E+02 bohr3Unit cell volume ucvol= 2.9199518E+02 bohr3Unit cell volume ucvol= 2.9954553E+02 bohr3Unit cell volume ucvol= 3.0722493E+02 bohr3然后用下面的命令：grep Etotal OUT得到如下的內(nèi)容： Etotal= -7.92750029752797E+00 Etotal= -7.92997465524506E+00 Etotal= -7.93167675973445E+00 Etotal= -7.93266612552653E+00 Etotal= -7.93299797094926E+00 Etotal= -7.93272412167288E+00 Etotal= -7.93189304516315E+00 Etotal= -7.93055314505404E+00 Etotal= -7.92874706916830E+00 Etotal= -7.92651686675655E+00 Etotal= -7.92390242137627E+00因此，Volume 和Etotal對(duì)應(yīng)的關(guān)系為：2.3610688E+02 -7.92750029752797E+002.4266422E+02 -7.92997465524506E+002.4934185E+02 -7.93167675973445E+002.5614088E+02 -7.93266612552653E+002.6306239E+02 -7.93299797094926E+002.7010748E+02 -7.93272412167288E+002.7727725E+02 -7.93189304516315E+002.8457279E+02 -7.93055314505404E+002.9199518E+02 -7.92874706916830E+002.9954553E+02 -7.92651686675655E+003.0722493E+02 -7.92390242137627E+00下面就可以用 Birch-Murnaghan 3階狀態(tài)方程進(jìn)行 (Birch F, Phys. Rev. 71, p809 (1947)擬合得到體彈性模量和平衡狀態(tài)下的體積：V0 = 263.276940709097 a.u.3B0 = 95.497(GPa)第二章 ABINIT參數(shù)設(shè)定與收斂測(cè)試/u2/85704/showart_1935274.html ABINIT計(jì)算晶體時(shí)主要參數(shù)的確定(切斷動(dòng)能和k點(diǎn)網(wǎng)格)。 在采用平面波贗勢(shì)法進(jìn)行固體的電子結(jié)構(gòu)計(jì)算時(shí)，為了確保計(jì)算的精度和計(jì)算結(jié)果的可靠性，在計(jì)算晶體的物理性質(zhì)之前，要進(jìn)行幾個(gè)重要參數(shù)的測(cè)試，以保證這些參數(shù)的選取使得計(jì)算結(jié)果有很好的收斂性，這些參數(shù)是平面波的切斷動(dòng)能和k點(diǎn)網(wǎng)格取樣的大小。 下面以采用ABINIT計(jì)算立方的ZrO2晶體為例： 采用的贗勢(shì)是40zr.psp_mod和8o.psp_mod(它們均是LDA的TM贗勢(shì))。 輸入文件，in.files的內(nèi)容為： INP OUT zroi zroo zro ./40zr.psp_mod ./8o.psp_mod 在測(cè)試平面波切斷動(dòng)能的收斂性時(shí)，我們通過設(shè)置平面波切斷動(dòng)能從20 Ha開始，以2 Ha遞增，直到58 Ha，其他參數(shù)不變計(jì)算ZrO2的總能。在ABINIT的輸入文件中很方便的通過ndtset來設(shè)置。 輸入文件INP如下（紫色標(biāo)示）： # Crystalline ZrO2-cubic ndtset 20 #表示有20組數(shù)據(jù) ecut: 20.0 #平面波切斷動(dòng)能從20 Ha開始 ecut+ 2 #以2 Ha遞增，也就是Ecut為20.0 + i*2.0, (i從1到20) #Definition of the unit cell acell 3*9.65285 #設(shè)置晶格常數(shù)a=b=c為9.65285 a.u. rprim 0.0 0.5 0.5 #同上面的acell確定了原胞的基矢 0.5 0.0 0.5 0.5 0.5 0.0 #Definition of the atom types ntypat 2 #設(shè)置原胞的原子種類數(shù)，這里有2類原子 znucl 40 8 #每類原子的核電荷數(shù)Z natom 3 #原胞中總的原子數(shù)目，這里原胞總共有3個(gè)原子 typat 1 2*2 #第一類原子有1個(gè)，第二類原子有2個(gè) xred #下面的按分?jǐn)?shù)坐標(biāo)給原胞中原子的坐標(biāo)位置 0.00 0.000 0.000 0.25 0.25 0.25 0.75 0.75 0.75 #Definition of the k-point grid kptopt 1 #設(shè)置生成k點(diǎn)的方法，這里表明有ngkpt和nshiftk來確定k點(diǎn)網(wǎng)格的大小 ngkpt 8 8 8 #設(shè)置對(duì)布里淵區(qū)進(jìn)行8x8x8網(wǎng)格的劃分 nshiftk 4 #對(duì)劃分得到的k點(diǎn)按下面的偏移量進(jìn)行平移 shiftk 0.5 0.5 0.5 0.5 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.5 #Definition of the SCF procedure #ecut 40.0 iscf 5 #自洽迭代中的算法，這里采用cg方法優(yōu)化。 toldfe 1.0d-10 #自洽迭代迭過程總能收斂的標(biāo)準(zhǔn) diemac 3.0 #體系宏觀的介電常數(shù) nstep 60 #最大自洽迭代步數(shù) #Definition of the outfile information prtwf 0 #不輸出波函數(shù)文件。 運(yùn)行下面的命令 abins OUT 對(duì)OUT采用進(jìn)行分析，用下面的命令取出一些數(shù)據(jù): grep ecut OUT summary grep Total energy(eV) OUT summary 得到summary類似如下數(shù)據(jù): 2.00000000E+01 -2.16713704994688E+03 2.20000000E+01 -2.17244914205186E+03 2.40000000E+01 -2.17629539439455E+03 2.60000000E+01 -2.17888687489711E+03 2.80000000E+01 -2.18045824630388E+03 3.00000000E+01 -2.18139622512928E+03 3.20000000E+01 -2.18192168325676E+03 3.40000000E+01 -2.18219516090481E+03 3.60000000E+01 -2.18232706572898E+03 3.80000000E+01 -2.18237668482667E+03 4.00000000E+01 -2.18239923089573E+03 4.20000000E+01 -2.18241955706857E+03 4.40000000E+01 -2.18244775901834E+03 4.60000000E+01 -2.18248714583389E+03 4.80000000E+01 -2.18253519502611E+03 5.00000000E+01 -2.18258901301065E+03 5.20000000E+01 -2.18264245486986E+03 5.40000000E+01 -2.18269485143783E+03 5.60000000E+01 -2.18274313138312E+03 5.80000000E+01 -2.18278585826504E+03 畫圖可以看到當(dāng)Ecut取40 Ha時(shí)，體系的總能有很好的收斂了。 下面測(cè)試對(duì)k點(diǎn)網(wǎng)格的，in.files的輸入文件同上。我們這里把k點(diǎn)網(wǎng)格從6x6x6開始增加，以2x2x2遞增。INP的輸入文件（紫色標(biāo)示）： # Crystalline ZrO2-cubic ndtset 10 ngkpt: 6 6 6 ngkpt+ 2 2 2 #Definition of the unit cell acell 3*9.65285 rprim 0.0 0.5 0.5 0.5 0.0 0.5 0.5 0.5 0.0 #Definition of the atom types ntypat 2 znucl 40 8 natom 3 typat 1 2*2 xred 0.00 0.000 0.000 0.25 0.25 0.25 0.75 0.75 0.75 #Definition of the k-point grid kptopt 1 #ngkpt 8 8 8 nshiftk 4 shiftk 0.5 0.5 0.5 0.5 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.5 #Definition of the SCF procedure ecut 40.0 iscf 5 toldfe 1.0d-10 diemac 3.0 nstep 60 #Definition of the outfile information prtwf 0 運(yùn)行 abinis OUT& 計(jì)算完后對(duì)OUT進(jìn)行分析， grep ngkpt OUT comment grep Total energy(eV) OUT comment 取出得到的數(shù)據(jù)如下： 28.0000 -2182.416564 60.0000 -2182.416514 110.0000 -2182.416506 182.0000 -2182.416526 280.0000 -2182.416499 408.0000 -2182.416512 570.0000 -2182.416506 770.0000 -2182.416514 1012.0000 -2182.416515 1300.0000 -2182.416504 畫圖可以看到8x8x8的k點(diǎn)網(wǎng)格就能保證體系的總能很好的收斂了。 第三章 ABINIT計(jì)算晶體的能帶結(jié)構(gòu)/user1/11542/archives/2006/1273410.shtml采用第一原理的電子結(jié)構(gòu)計(jì)算方法來計(jì)算晶體的能帶結(jié)構(gòu)一般來說，要進(jìn)行兩個(gè)步驟。這不論是采用VASP、PWSCF還是ABINIT這些程序。步驟為：先進(jìn)行自洽的電子結(jié)構(gòu)迭代得到自洽計(jì)算得到的電荷密度，然后讀入這個(gè)自洽得到的電荷密度，進(jìn)行非自洽的計(jì)算得到體系的本征值。 下面以立方的ZrO2晶體為例，采用ABINIT來計(jì)算。下面采用ndset這個(gè)關(guān)鍵詞，在輸入文件中輸入這兩步的控制參數(shù)，一次性計(jì)算完成得到能帶結(jié)構(gòu)。并對(duì)相關(guān)的相關(guān)的輸入?yún)?shù)進(jìn)行解釋： 采用的贗勢(shì)是:40zr.psp_mod和8o.psp_mod，輸入文件 in.files的內(nèi)容為 INP OUT zroi zroo zro ./40zr.psp_mod ./8o.psp_mod 此in.files文件的內(nèi)容設(shè)置了主要輸入文件的名稱和贗勢(shì)文件所在的目錄。 主要輸入文件INP的內(nèi)容為（紫色標(biāo)示）： # Crystalline ZrO2-cubic ndtset 2 #表示有兩組控制參數(shù)：第一組控制參數(shù)用來設(shè)置自洽計(jì)算，第二組是非自洽的本征值計(jì)算 #Dataset 1 : usual self-consistent calculation kptopt1 1 #自洽計(jì)算中設(shè)置k點(diǎn)網(wǎng)格取樣的方法，表明采用ngkpt和shfitk來設(shè)置k點(diǎn)網(wǎng)格。 nshiftk1 4 #使生成的k點(diǎn)進(jìn)行平移。 shiftk1 0.5 0.5 0.5 0.5 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.5 ngkpt1 8 8 8 #K點(diǎn)網(wǎng)格取樣，網(wǎng)格劃分的分割數(shù) prtden1 1 #表明輸出電荷密度文件。 tolvrs1 1.0d-20 # 自洽計(jì)算收斂的標(biāo)準(zhǔn) iscf1 5 #自洽迭代計(jì)算時(shí)，采用 CG方法來優(yōu)化有效勢(shì) #Dataset 2 : #用來設(shè)置本征值計(jì)算時(shí)的參數(shù) iscf2 -2 #表明非自洽計(jì)算 getden2 -1 #讀入上一組數(shù)據(jù)進(jìn)行自洽計(jì)算得到的電荷密度文件 kptopt2 -5 #負(fù)數(shù)，表示下面計(jì)算能帶計(jì)算時(shí)，有5段特殊線(由6個(gè)特殊k點(diǎn)來確定) ndivk2 10 12 18 8 8 #每段特殊線上分幾等份 kptbounds2 #特殊k點(diǎn)的坐標(biāo) 0.5 0.25 0.75 # W point 0.5 0.0 0.0 # L point 0.0 0.0 0.0 # Gamma point 0.5 0.0 0.5 # X point 0.5 0.25 0.75 # W point 0.0 0.0 0.0 # Gamma enunit2 0 #輸出本征值時(shí)，本征值的單位，這里為0表示是以Hatree為單位給出 prteig2 1 #表明輸出本征值到文件中 #Definition of the unit cell acell 3*9.48196 #設(shè)置晶常數(shù) rprim 0.0 0.5 0.5 #設(shè)置計(jì)算原胞的基矢(同上面的acell一起構(gòu)成原胞的基矢) 0.5 0.0 0.5 0.5 0.5 0.0 #Definition of the atom types ntypat 2 #原胞的原子種類數(shù)目 znucl 40 8 #每類原子的核電荷數(shù) natom 3 #原胞中總的原子數(shù)目 typat 1 2*2 #表示第一個(gè)原子是第一類的，后面兩個(gè)原子是第二類的。 xred #以分?jǐn)?shù)坐標(biāo)給出原子的位置 0.00 0.000 0.000 0.25 0.25 0.25 0.75 0.75 0.75 #Definition of the SCF procedure nband 30 #在計(jì)算中考慮多少條能帶 ecut 40.0 #平面波切斷動(dòng)能 diemac 3.0 #體系的宏觀介電常數(shù)，給一個(gè)近似值就可以了。 nstep 60 #自洽迭代時(shí)的最大步數(shù) 計(jì)算得到的本征值文件zroo_DS2_EIG內(nèi)容為： Eigenvalues (hartree) for nkpt= 57 k points: kpt# 1, nband= 30, wtk= 1.00000, kpt= 0.5000 0.2500 0.7500 (reduced coord) -1.41786 -0.62938 -0.62938 -0.61620 -0.26378 -0.26378 0.14567 0.14721 0.23910 0.23910 0.27711 0.32275 0.49670 0.50302 0.57458 0.57458 0.61760 0.78769 0.90930 0.92430 0.92430 1.05159 1.05159 1.10749 1.20089 1.28891 1.28891 1.35618 1.42279 1.46219 kpt# 2, nband= 30, wtk= 1.00000, kpt= 0.5000 0.2250 0.6750 (reduced coord) -1.41788 -0.63312 -0.62545 -0.61596 -0.26959 -0.25851 0.13111 0.16454 0.22812 0.24784 0.27846 0.32146 0.49126 0.50549 0.56954 0.58345 0.61790 0.78670 0.89599 0.91521 0.94726 1.02556 1.08306 1.08822 1.21602 1.28607 1.29698 1.37003 1.41244 1.48297 。 把它們轉(zhuǎn)換為k點(diǎn)的絕對(duì)值和本征值后，畫圖如下：附1 處理abinit計(jì)算得到的本征值文件/user1/11542/archives/2006/1275248.shtml在采用abinit計(jì)算固體的能帶結(jié)構(gòu)時(shí)，abinit輸出的本征值文件格式為： Eigenvalues (hartree) for nkpt= 27 k points: kpt# 1, nband= 20, wtk= 1.00000, kpt= 0.0000 0.0000 0.0000 (reduced coord) -0.37264 -0.18371 -0.15298 -0.15298 -0.11442 -0.01609 0.01563 0.02497 0.03905 0.03905 0.04309 0.05183 0.05183 0.07649 0.07991 0.07991 0.08411 0.08411 0.08768 0.11715 kpt# 2, nband= 20, wtk= 1.00000, kpt= 0.0000 0.0000 0.0192 (reduced coord) -0.37255 -0.18459 -0.15289 -0.15289 -0.11343 -0.01601 0.01572 0.02505 0.03913 0.03913 0.04317 0.04982 0.04982 0.07658 0.08000 0.08000 0.08627 0.08627 0.08776 0.11724 。 第一行告訴你輸出的本征值的單位是Hatree，共計(jì)算了27個(gè)特殊k點(diǎn)。 下面是分別給出每個(gè)k點(diǎn)的的坐標(biāo)，以及所計(jì)算的本征值。 但是要把這些本征值和k點(diǎn)的坐標(biāo)對(duì)應(yīng)起來畫圖，怎么辦呢？這里提供一個(gè)我根據(jù)pwscf程序包中的一個(gè)代碼改寫后的小工具plotband_abinit.f90，它用來把每個(gè)k點(diǎn)與前一個(gè)k點(diǎn)之間的間隔一一加起來得到的距離和本征值對(duì)應(yīng)起來。 plotband_abinit.f90的內(nèi)容如下： ! for ABINIT ! program prog real, allocatable : e(:,:) real, allocatable : k(:,:) real, dimension(3) :k0,a character(len=32): input, output character(len=32): xx, yy write(6,*) number of calculated bands read(5,*) nbands write(6,*) number of bands to be plotted read(5,*) nbands2 write(6,*) number of k-points read(5,*) nk write(6,*) fermi level (eV) read(5,*) ef write(6,*) name bande.in e bande.out read(5,*) input,output allocate(e(nk,nbands) allocate(k(nk,3) open(10,file=input, status=old) open(11,file=output, status=new) read(10,*) xx, xx,xx,xx, nkpt, xx, xx do i=1,nk read(10,*)yy, ik, yy, nb, yy, wtk, yy, (k(i,j),j=1,3),yy write(6,*)(k(i,j),j=1,3) read(10,*) (e(i,n),n=1,nbands) write(13,9030) (e(i,n),n=1,nbands) 9030 format (8f9.5) enddo do j=1,nbands2 dk=0 do i=1,nk if (i.eq.1) then k0=k(i,:) endif a=k(i,:)-k0 dk=dk+sqrt(dot_product(a,a) write(11,*)dk,(e(i,j)-ef)*27.2116 k0=k(i,:) enddo write(11,*) enddo stop end program prog 采用任何一種fortran90的編譯器(比如pgf90和ifc，ifort)編譯后。 pgf90 plotband_abinit.f90 它的輸入文件為計(jì)算完后的本征值。注意這里是針對(duì)abinit計(jì)算輸出的本征值單位為Hatree。要補(bǔ)充的是：abinit再輸出本征值文件時(shí)有個(gè)小小的bug，在以eV單位輸出本征值時(shí)，輸出的本征值文件中會(huì)漏掉第51個(gè)的本征值。因此建議在采用abinit計(jì)算輸出本征值時(shí)以Hatree為單位。 編譯后得到a.out，運(yùn)行a.out按提示輸入相應(yīng)的值； 第一次提示要輸入的量是你計(jì)算了多少條能帶；第二次提示要輸入的量是你要畫出多少條能帶。 第三次提示要輸入的量是你計(jì)算了多少個(gè)k點(diǎn)；第四次提示要輸入的是體系的費(fèi)米能級(jí)。(注意不要從自洽計(jì)算得到的 OUT中取出費(fèi)米能級(jí)，abinit給出的有錯(cuò)誤。而是要從計(jì)算態(tài)密度時(shí)的得到的態(tài)密度文件中取出。) ；第五此提示要輸入的是包含了計(jì)算得到的本征值的文件名，以及轉(zhuǎn)換后輸出的文件名。 然后就可以對(duì)轉(zhuǎn)換后的文件進(jìn)行畫圖了。附2 不同晶格的簡(jiǎn)約布里淵區(qū)形狀及特殊K點(diǎn)坐標(biāo)1. FCC晶格的簡(jiǎn)約布里淵區(qū)形狀及特殊K點(diǎn)坐標(biāo)2. BCC晶格的簡(jiǎn)約布里淵區(qū)形狀及特殊K點(diǎn)坐標(biāo)3. SC晶格的簡(jiǎn)約布里淵區(qū)形狀及特殊K點(diǎn)坐標(biāo)4. 六角晶格的簡(jiǎn)約布里淵區(qū)形狀及特殊K點(diǎn)坐標(biāo)5. 正交晶格的簡(jiǎn)約布里淵區(qū)行狀及特殊k點(diǎn)坐標(biāo)以上是采用xcrysden軟件畫出了各種晶格的簡(jiǎn)約布里淵區(qū)形狀，并標(biāo)出了幾個(gè)特殊K點(diǎn)的位置、標(biāo)記及其坐標(biāo)（以倒格子基矢為坐標(biāo)系）。第四章 ABINIT計(jì)算電子態(tài)密度/user1/11542/archives/2007/1664470.shtml電子態(tài)密度(Density of States, DOS)定義為單位能量范圍內(nèi)所允許的能級(jí)數(shù)。在電子結(jié)構(gòu)計(jì)算程序中進(jìn)行計(jì)算時(shí)，一般是按如下的步驟來進(jìn)行的：先采用小的k點(diǎn)網(wǎng)格進(jìn)行自洽計(jì)算，然后采用較密的k點(diǎn)網(wǎng)格進(jìn)行非自洽計(jì)算。下面介紹采用ABINIT程序電子態(tài)密度時(shí)的步驟，解釋需要輸入的關(guān)鍵詞，以及計(jì)算結(jié)果的處理說明。這里以采用ABINIT5.2.4版本計(jì)算MgO電子態(tài)密度（包括總態(tài)密度和spd分波態(tài)密度等)為例子。在例子中MgO的晶格常數(shù)(它是面心立方晶體，NaCl結(jié)構(gòu)，a=7.921 a.u.)1取為實(shí)驗(yàn)值，平面波切斷動(dòng)能和k點(diǎn)網(wǎng)格大小分別取為24 Ha，10x10x10（自洽計(jì)算中）和16x16x16（態(tài)密度計(jì)算中）。1 O.L. Anderson and P. Andreatch, J. Am. Ceram. Soc. 49, 404 (1966).1. 計(jì)算總態(tài)密度 通過abinit中的ndtset設(shè)置兩組數(shù)據(jù)進(jìn)行計(jì)算，先是自洽計(jì)算，后是態(tài)密度計(jì)算的設(shè)置。在計(jì)算態(tài)密度時(shí)，需將prtdos設(shè)置為1、2或3。當(dāng)prtdos設(shè)置1時(shí)，采用離散(smearing)的方法來確定態(tài)密度，即由本征值、k點(diǎn)權(quán)重(參數(shù)wtk的值)、每個(gè)本征值的占有數(shù)（由occopt設(shè)置的參數(shù)所選擇的方法來確定占有數(shù))以及離散展寬參數(shù)(tsmear參數(shù)設(shè)置的值)來計(jì)算，此時(shí)需要注意或設(shè)置wtk、occopt和tsmear參數(shù)。一般對(duì)一維體系或者在k點(diǎn)較少的情況采用prtdos=1的方法來計(jì)算態(tài)密度，以得到較為光滑的DOS曲線。當(dāng)設(shè)置prtdos=2時(shí)，表示采用四面體方法來計(jì)算總態(tài)密度，此時(shí)不需要設(shè)置occopt和tsmear了，所采用的k點(diǎn)需要由參數(shù)ngkpt或kptrlatt參數(shù)來確定。當(dāng)prtdos=3時(shí)，表示采用四面體方法計(jì)算分波態(tài)密度(angular-momentum projected)，同時(shí)需設(shè)置原子球半徑。以計(jì)算面心立方MgO的總態(tài)密度為例子。i). in.files輸入文件的內(nèi)容:dos.outmgoimgoomgo12mg.pspnc8o.pspncii)主要輸入文件tdos.in的內(nèi)容（紫色標(biāo)注）：ndtset 2kptopt 1 nshiftk1 4shiftk1 0.5 0.5 0.50.5 0.0 0.00.0 0.5 0.00.0 0.0 0.5ngkpt1 8 8 8prtden1 1 toldfe1 1.0d-6iscf1 5 #Dataset 2 : DOS calculationiscf2 -3getden2 1getwfk2 1prtdos2 2dosdeltae 0.00005 ngkpt2 16 16 16tolwfr2 1.0d-16shiftk20.0 0.0 0.00.0 0.5 0.50.5 0.0 0.50.5 0.5 0.0#Definition of the unit cellacell 3*7.921rprim0.0 0.5 0.50.5 0.0 0.50.5 0.5 0.0ntypat 2znucl 12 8natom 2typat 1 2 xred0.0 0.0 0.00.5 0.5 0.5nband 12ixc 1 ecut 24.0nstep 45diemac 9.0 這里利用ndtset設(shè)置兩組數(shù)據(jù)，第一組是進(jìn)行自洽計(jì)算以得到電荷密度，第二組是非自洽計(jì)算以得到總態(tài)密度。在第一組數(shù)據(jù)中:kptopt 1 #設(shè)置通過ngkpt由程序自動(dòng)生成k點(diǎn)nshiftk1 4 #由于是fcc結(jié)構(gòu)，為了使得的產(chǎn)生的k點(diǎn)對(duì)計(jì)算更有效，這里將所產(chǎn)生的k點(diǎn)進(jìn)行移動(dòng)shiftk1 #設(shè)置了移動(dòng)的量0.5 0.5 0.50.5 0.0 0.00.0 0.5 0.00.0 0.0 0.5ngkpt1 8 8 8 #設(shè)置k點(diǎn)網(wǎng)格的大小。由于是在自洽計(jì)算中，此時(shí)設(shè)置的數(shù)較小。prtden1 1 #表示要將電荷密度輸出到文件中，以供第二組中的計(jì)算用到。toldfe 1.0d-6 #自洽迭代計(jì)算中總能變化收斂的標(biāo)準(zhǔn)iscf1 5 #表示采用CG方法來優(yōu)化總能以得到基態(tài)在第二組數(shù)據(jù)中：iscf2 -3 #表示進(jìn)行非自洽計(jì)算，一般用在態(tài)密度和STM計(jì)算以及線性響應(yīng)計(jì)算中對(duì)ddk微擾的計(jì)算。getden2 1 #表示從上一步中得到的電荷密度作為初始的電荷密度讀入getwfk2 1 #表示讀入上一步的波函數(shù)prtdos2 2 #表示要輸出總態(tài)密度（或進(jìn)行總態(tài)密度計(jì)算)dosdeltae 0.00005 #態(tài)密度計(jì)算時(shí)能量軸上的能量刻度ngkpt2 16 16 16 #在態(tài)密度計(jì)算時(shí)，k點(diǎn)網(wǎng)格應(yīng)該設(shè)置的密一些。tolwfr2 1.0d-16shiftk20.0 0.0 0.00.0 0.5 0.50.5 0.0 0.50.5 0.5 0.0其他的參數(shù)是用來設(shè)置晶格結(jié)構(gòu)以及原子的種類和坐標(biāo)，平面波切斷動(dòng)能，自洽收斂參數(shù)。在計(jì)算出來的_DOS文件中給出了總態(tài)密度的值，其格式如下：ABINIT package : DOS filensppol = 1, nkpt = 145, nband(1)= 12Tetrahedron methodFor identification : eigen(1:3)= -0.429 0.203 0.203Fermi energy : 0.20340927The DOS (in electrons/Hartree/cell) and integrated DOS (in electrons/cell) are computed,at 38001 energies (in Hartree) covering the intervalbetween -0.5000 and 1.4000 Hartree by steps of 0.00005 Hartree.index energy(Ha) DOS integrated DOS0 -0.50000 0.0000 0.00001 -0.49995 0.0000 0.00002 -0.49990 0.0000 0.0000.index給出的是能量點(diǎn)數(shù)的序號(hào)，第二列給出的是能量的數(shù)值（單位是Hatree），在畫圖是一般以eV，并減去前面給出的Fermi energy。第三列是總態(tài)密度值，第四列是電子數(shù)。2. 計(jì)算各個(gè)原子的分波態(tài)密度（紫色標(biāo)注）ndtset 2kptopt 1nshiftk1 4shiftk10.5 0.5 0.50.5 0.0 0.00.0 0.5 0.00.0 0.0 0.5ngkpt1 8 8 8prtden1 1toldfe1 1.0d-6iscf1 5#Dataset 2 : DOS c