VASP的全稱Vienna Ab-initio Simulation Package,是維也納大學Hafner小組開發的進行電子結構計算和量子力學-分子動力學模擬軟體包。它是目前材料模擬和計算物質科學研究中最流行的商用軟體之一。與Material Studio軟體包中的CASTEP功能類似,但是VASP的精度相對要高一點。不同於CASTEP的圖形界面,VASP是一套沒有界面的計算軟體,建模、可視化、數據分析都需要依賴第三方工具如P4VASP、ASE、Pymatgen、VESTA軟體等。VESTA、P4VASP主要是用來建模、可視化和分析部分數據。而ASE、Pymatgen這些軟體擅長於數據處理,但是安裝比較麻煩,同時入門門檻比較高,需要使用者有一定的編程水平。VASP用戶的學科分布很廣,有做催化的,有做光學的,有做材料的,各個領域的數據後處理方式大相逕庭。很多用戶開發並貢獻了自己所在領域用到的的腳本或者小程序,本人就開發了一款用來處理結構文件的POSCARtookit腳本。但是對於新用戶來說,找到並成功使用這些腳本是不太容易的。因此一款容易上手、功能強大的預-後數據處理軟體VASPKIT應運而生。
最新版的VASPKIT是王偉老師、許楠、劉錦程,唐剛和李強共同努力的成果。我作為貢獻者近期參與了VASPKIT項目的維護與開發工作,後期將會和諸多開發者一起增加催化相關的功能。VASPKIT release 版本是一款用FORTRAN編寫,在LINUX環境下運行的二進位軟體。它幾乎不依賴於其他庫,軟體體積僅僅4.5M,無需安裝即可使用,同時EXAMPLES目錄下面有主要功能的測試例子,方便用戶學習使用。
主要功能有:1.生成不同計算任務的INCAR文件 2.結構對稱性查找 3.生成倒空間K點網格 4.生成晶體的能帶路徑(包括雜化泛函),並處理能帶數據 5.處理態密度DOS和投影態密度PDOS 6.處理電荷密度、靜電勢,繪製是空間波函數 7.其他功能,比如熱力學量校正(新增),光學、分子動力學、導電率和半導體方面的小工具。
2. VASPKIT的配置和使用vaspkit是一款運行在LINUX環境下的軟體,為了確保能夠使用vaspkit的完整功能。用戶可以配置vaspkit的環境目錄,在下一次運行時生效。通過在bash終端下運行以下命令將環境變量文件複製到用戶目錄下。
cp how_to_set_environment_variable ~/.vaspkit
並編輯.vaspkit文件
vi ~/.vaspkit
該配置文件主要用來設置vaspkit的環境變量,包括VASP版本信息,贗勢庫的目錄,泛函方法的選擇,並選擇是否按照VASP官方的推薦生成元素的贗勢文件,設置生成的INCAR模板是覆蓋、追加還是備份原有的INCAR。
# cp how_to_set_environment_variable ~/.vaspkit and modify the ~/.vaspkit file based on your settings!
VASP5 .TRUE. # .TRUE. or .FALSE.; Set .FALSE. if you are using vasp.4.x
GGA_PATH
PBE_PATH
LDA_PATH
POTCAR_TYPE PBE # PBE, GGA or LDA;
GW_POTCAR .FALSE. # .TRUE. or .FALSE.;
RECOMMENDED_POTCAR .TRUE. # .TRUE. or .FALSE.;
SET_FERMI_ENERGY_ZERO .TRUE. # .TRUE. or .FALSE.;
SET_INCAR_WRITE_MODE OVERRIDE # OVERRIDE, APPEND,BACK-UP-OLD,BACK-UP-NEW;
設置好POTCAR的目錄和完成一些其他的設置後,就可以啟動vaspkit了。
為了方便,可以將vaspkit的絕對路徑加入到環境變量裡,如果是LINUX的用戶,可以這樣操作:
echo 'export PATH=$PATH:/home/vaspkit.0.70/bin/vaspkit' >> ~/.bashrc
source ~/.bashrc
其中/home/vaspkit.0.70/bin/vaspkit用自己的vaspkit絕對路徑替代。
在終端直接輸入vaspkit或者/home/vaspkit.0.70/bin/vaspkit開始運行vaspkit程序。不出意外,你將會得到一個與我展示的一致的界面:
+---+
| VASPKIT Version: 0.70 (06 Nov. 2018) |
| A Pre- and Post-Processing Program for VASP Code |
| Official Website: http://vaspkit.sourceforge.net |
| Developed By Vei WANG (wangvei@me.com) |
| Contributor: Nan XU (tamas@zju.edu.cn) |
+---+
===================== Structural Options ========================
1) Customize INCAR File 2) Elastic-Constant Calculator
3) Poscar2other Toolkit 4) Build Supercell
5) Equation-of-State Fitting 6) Structure Toolkit
7) K-Mesh Generator 8) Band-Path Generator
===================== Electronic Options ========================
11) Total Density-of-States 12) Projected Density-of-States
21) Band-Structure 22) Projected Band-Structure
23) 3D Band-Structure 25) Hybrid-DFT Band-Structure
26) Fermi-Surface Calculator
=========== Charge & Potential & Wavefunction Options ===========
31) Charge Density 32) Spin Density
33) Spin-Up & -Down Density 34) Charge-Density Difference
41) Planar-Average Charge-Density 42) Planar-Average Potential
51) Real-Space WaveFunction Plot
====================== Misc Utilities ===========================
71) Linear Optics 72) Molecular-Dynamics Toolkit
73) VASP2BoltzTraP Interface 74) Energy-Related Toolkit
91) Semiconductor Calculator 92) 2D-Materials Toolkit
0) Quit
-->>
如果出現以下問題,說明你的LINUX運行依賴庫版本太低,需要升級(不建議),可以聯繫開發者獲得在低版本LINUX環境下編譯的版本。
vaspkit: /lib64/libc.so.6: version `GLIBC_2.14' not found (required by vaspkit)
如果出現-bash: line 7: ./vaspkit: Permission denied權限問題,只需賦予vaspkit執行權限即可:
chmod u+x /home/vaspkit.0.70/bin/vaspkit
本教程將介紹使用vaspkit生成VASP的輸入文件和很多用戶關心的能帶計算及熱力學量校正功能,能帶計算包含了使用普通泛函和雜化泛函兩個例子,而熱力學量校正主要計算零點振動能及溫度對自由能和焓的貢獻。
3.1 VASP輸入文件的生成為了成功運行VASP計算任務,我們至少需要4個文件:INCAR、POSCAR、POTCAR及KPOINTS,INCAR是告訴 VASP算什麼任務,怎麼算的控制文件;POSCAR是包含晶格信息,原子坐標信息和原子速度信息(MD用)的文件;POTCAR是贗勢文件,也就是將內層電子用勢函數表示;KPOINTS(可包含在INCAR內,不推薦省略)包含了倒易空間內K點信息,波函數會在這些點上積分得到電荷密度。POSCAR一般由軟體生成或者從資料庫中獲得,簡單體系可自己搭建。比如從資料庫(COD)中獲得磷化鎵的CIF文件後,可以通過VTST腳本cif2pos.pl或者VESTA轉化成POSCAR文件(原子位置分數佔據的問題需要注意)。接下來優化的原子坐標得到穩定的結構,我們採用vaspkit預設的INCAR組合生成所需的INCAR文件。在有POSCAR的目錄下運行vaspkit,輸入1選擇功能Customize INCAR File,會得到以下的顯示信息:
+
You MUST Know What You Are Doing
Some Parameters in INCAR File Neet To Be Set/Adjusted Manually
+
======================== INCAR Options ==========================
ST) Static-Calculation SR) Standard Relaxation
MG) Magnetic Properties SO) Spin-Orbit Coupling
D3) DFT-D3 no-damping Correction H6) HSE06 Calculation
PU) DFT+U Calculation MD) Molecular Dynamics
GW) GW0 Calculation BS) BSE Calculation
DC) Elastic Constant EL) ELF Calculation
BD) Bader Charge Analysis OP) Optical Properties
EC) Static Dielectric Constant PC) Decomposed Charge Density
FD) Phonon-Finite-Displacement DT) Phonon-DFPT
NE) Nudged Elastic Band (NEB) DM) The Dimer Method
FQ) Frequence Calculations LR) Lattice Relaxation
0) Quit
9) Back
Input Key-Parameters (STH6D2 means HSE06-D2 Static-Calcualtion)
由於進行能量弛豫任務,故輸入SR,就會得到一個預設好的用於做弛豫任務的INCAR(有些模板需要手動修改。比如DFT+U的U值設定,NEB的IMAGES數目等)。如果已經有INCAR文件,則原來的INCAR文件會被覆蓋。你可以編輯~/.vaspkit更改INCAR的寫出設置。只需將最後一行的SET_INCAR_WRITE_MODE由默認的override更改為 append,back-up-old,back-up-new中的一個,分別對應於新的內容增加到原有的INCAR後面,備份原有的INCAR再寫入新的INCAR和寫入到新的INCAR.new裡面。
接下來生成KPOINTS文件。對於非能帶計算,只需用程序自動撒點的方式,但是需要用戶選擇撒點方式和K點密度。具體內容可以參考李強的教程Learn VASP The Hard Way (Ex1): VASP基本輸入文件的準備。啟動vaspkit,輸入7選擇功能K-Mesh Generator,接下來輸入2選擇Gamma Scheme撒點方式(穩妥的選擇),會得到以下的顯示信息:
-->> (1) Reading Lattices & Atomic-Positions from POSCAR File...
+- Warm Tips -+
* Accuracy Levels: (1) Low: 0.06~0.04;
(2) Medium: 0.04~0.03;
(3) Fine: 0.02-0.01.
(4) Gamma-Only: 0.
* 0.04 is Generally Precise Enough!
+---+
Input KP-Resolved Value (unit: 2*PI/Angstrom):
-->>
根據提示0.04已經足夠精確了,因此輸入0.04,將會得在當前目錄下得到KPOINTS和POTCAR文件。
KPOINTS內容如下:
K-Mesh Generated with KP-Resolved Value : 0.040
0
Gamma
8 8 8
0.0 0.0 0.0
Learn VASP The Hard Way Ex19 誰偷走的我的機時?(三)提到了一個簡單的K點數目判斷準則,對於半導體k點數目乘實空間對應晶矢量大於20Å(參考下圖),本例中ka=21.802Å已經符合經驗規律。劉錦程提到對於非正交體系 ,倒格矢長度和晶常數不滿足反比關係,所以採用ka≈kb≈kc的經驗準則並不能保證 K點密度在各個方向長相等。而vaspkit嚴格計算了倒空間晶格矢量比例,因此使用vaspkit就能根據所選的K點密度自動生成各個方向的K點數。
grep ENMAX POTCAR
可以看到默認的INCAR參數中ENCUT=400eV被注釋掉了,但是保留了PREC = Normal,程序會自動將ENCUT設為max(ENMAX)。當然也可以自行設置ENCUT參數,只要參數大於所有元素的ENMAX,這時候以自己設定的ENCUT參數為準。注意在優化晶胞常數時,需要用較高的ENCUT(Learn VASP The Hard Way (Ex36):直接優化晶格常數),因此LQ任務(優化晶格常數)模板生成的INCAR中默認設置了PREC=High,程序會自動將ENCUT設為max(ENMAX)*1.3。
3.2 能帶計算能帶計算的KPOINTS與普通計算的KPOINTS不一樣,通常需要第一布裡淵區內的一條或幾條高對稱點路徑來計算能帶性質。 傳統的做法是通過SeeK-Path網站(https://www.materialscloud.org/tools/seekpath/)或者Material Studio軟體獲得晶體倒易空間第一布裡淵區內的高對稱點,再通過腳本插值生成高對稱點路徑上的K點,得到滿足要求的KPOINTS。好消息是新版的VASPKIT集成了與SeeK-path一致的算法分析晶體的高對稱點,可以方便地生成PBE泛函和HSE06雜化泛函所需的KPOINTS。
在vaspkit.0.70/examples/hybrid_DFT_band目錄下有一個使用HSE06雜化泛函計算磷化鎵的能帶結構的例子。同樣也可以使用PBE泛函計算該磷化鎵結構的能帶,只是普通的PBE泛函會低估帶隙。為了計算能帶,首先得獲得晶體的第一布裡淵區內的一條或幾條高對稱點路徑。在有POSCAR的目錄下運行vaspkit,輸入8選擇功能Band-Path Generator,在下一個界面輸入3選擇3D bulk structure (Experimental),你會得到以下信息:
+
See An Example in vaspkit/examples/seek_kpath.
This Feature Is Experimental & Check Your System with SeeK-Path.
For More details See [www.materialscloud.org/work/tools/seekpath].
+
+
Prototype: AB
Total Atoms in Input Cell: 2
Lattice Constants in Input Cell: 3.854 3.854 3.854
Lattice Angles in Input Cell: 60.000 60.000 60.000
Total Atoms in Primitive Cell: 2
Lattice Constants in Primitive Cell: 3.854 3.854 3.854
Lattice Angles in Primitive Cell: 60.000 60.000 60.000
Crystal System: Cubic
Crystal Class: -43m
Bravais Lattice: cF
Extended Bravais Lattice: cF2
Space Group: 216
Point Group: 31 [ Td ]
International: F-43m
Symmetry Operations: 24
Suggested K-Path: (shown in the next line)
[ GAMMA-X-U|K-GAMMA-L-W-X ]
+
+
| * DISCLAIMER * |
| CHECK Your Results for Consistency if Necessary |
| Bug Reports and Suggestions for Improvements Are Welcome |
| Citation of VASPKIT Is Not Mandatory BUT Would Be Appreciated |
+
VASPKIT會分析晶體的對稱性並得到兩條建議的能帶路徑[ GAMMA-X-U|K-GAMMA-L-W-X ],同時生成了晶體的單胞結構PRIMCELL.vasp,並生成了KPATH.in用於能帶結構計算。KPATH.in的第二行數字表示了每小段路徑中插值的K點的數目,如果默認的數值都算不動的話,可以考慮將其設小。能帶路徑只針對於primitive cell,因此需要執行下面的命令,用生成的primitive cell作為計算的結構文件。
\cp -f PRIMCELL.vasp POSCAR
下圖展示的是磷化鎵的primitive cell和其第一布裡淵區的高對稱點位置,由SeeK-path網站生成。
PBE泛函計算能帶分為兩步,第一步使用普通K點網格(主功能7)進行自洽計算 ,然後重啟vaspkit,輸入1選擇Customize INCAR File功能,輸入ST,生成靜態自洽的INCAR。本例中ISMEAR=0,即Gaussian Smearing方法,如果是金屬體系可以選擇換成ISMEAR=1。接著調用VASP計算。第二步:使用KPATH.in裡的高對稱點信息作為新的KPOINTS,然後讀入電荷CHGCAR進行能帶非自洽計算,即:
\cp -f KPATH.in KPOINTS
echo "ICHARG=11" >> INCAR
也就是讀入上一步生成的CHGCAR並保持不變,調用VASP計算。
待第二步完成後,通過功能21從能量本徵值文件中讀入能帶結構。
雜化泛函相比於PBE泛函和DFT+U方法,在計算帶隙方面很有優勢,但是HSE06的雜化泛函需要KPOINTS裡既有權重不為0的K點進行自洽計算,又要求有權重為0的高對稱點計算能帶性質。因此操作流程頗為繁瑣。
重啟vaspkit,輸入25選擇功能Hybrid-DFT Band-Structure,在下一個界面輸入251選擇Generate KPOINTS File for Hybrid Band-Structure Calculation。再輸入1選擇Monkhorst-Pack Scheme用MP方法生成自洽用的K點網格並根據建議輸入0.04選擇較密的K點密度(權重不為0的K點)即可生成HSE06計算所需要的KPOINTS。
因為HSE06計算非常耗時,因此本例中採用兩步法加速收斂,當然也可以跳過第一步直接用HSE06進行自洽。第一步:使用PBE泛函產生波函數和電子密度,第二步:保持KPOINTS不變,讀入波函數進行HSE06計算。
利用功能Customize INCAR File生成第一步PBE自洽需要的INCR。重啟vaspkit,輸入1選擇Customize INCAR File功能,輸入ST,生成靜態自洽的INCAR。本例中ISMEAR=0,即Gaussian Smearing方法,如果是金屬體系可以選擇換成ISMEAR=1。調用VASP計算,待自洽完成後執行第二步HSE06計算。重啟vaspkit,輸入1選擇Customize INCAR File功能,輸入STH6,生成HSE06計算所需要的INCAR。再次調用VASP計算後,就完成了HSE06的自洽計算。
接下來使用vaspkit提取能帶數據,並輸出高對稱點在能帶圖中的坐標信息。輸入25選擇功能Hybrid-DFT Band-Structure,在下一個界面輸入252選擇Read Band-Structure for Hybrid-DFT Calculation處理能帶數據。處理結果如下所示:
-->> (1) Reading Input Parameters From INCAR File...
-->> (2) Reading Fermi-Level From DOSCAR File...
-->> (3) Reading Energy-Levels From EIGENVAL File...
-->> (4) Reading Lattices & Atomic-Positions from POSCAR File...
-->> (5) Reading K-Paths From KPATH.in File...
-->> (6) Written BAND.dat File!
-->> (7) Written BAND-REFORMATTED.dat File!
-->> (8) Written KLINES.dat File!
-->> (9) Written KLABELS File!
BAND-REFORMATTED.dat是我修改後生成的能帶信息,格式如下所示,第一列是K-PATH,相鄰的距離為高對稱點在倒空間的距離,第二列,第三列等即為所需的不同的能帶線。如果開了自旋,第二列和第三列為能帶1的spin up和spin down,以此類推。將dat文件拖入Origin軟體後,即可得到能帶圖。
# Kpath Energy-Level (in eV)
0.00000 -20.67711 -19.12020 -19.03327
0.07227 -20.66779 -18.93770 -19.03277
0.14454 -20.63874 -19.02325 -19.03142
0.21680 -20.59133 -19.66623 -19.02948
KLABELS是通過我參與修改的讀取能帶標識的子模塊生成的。通過讀取KPATH.in(PBE能帶計算時為KPOINTS)高對稱點後標識,寫入KLABELS文件。
KPATH.in文件內容如下:
Generated by VASPKIT based on Hinuma et al.
10
Line-Mode
Reciprocal
0.0000000000 0.0000000000 0.0000000000 GAMMA
0.5000000000 0.0000000000 0.5000000000 X
0.5000000000 0.0000000000 0.5000000000 X
0.6250000000 0.2500000000 0.6250000000 U
0.3750000000 0.3750000000 0.7500000000 K
0.0000000000 0.0000000000 0.0000000000 GAMMA
0.0000000000 0.0000000000 0.0000000000 GAMMA
0.5000000000 0.5000000000 0.5000000000 L
0.5000000000 0.5000000000 0.5000000000 L
0.5000000000 0.2500000000 0.7500000000 W
0.5000000000 0.2500000000 0.7500000000 W
0.5000000000 0.0000000000 0.5000000000 X
KLABELS文件如下:
K-Label K-Coordinate in band-structure plots
GAMMA 0.000
X 1.153
U|K 1.560
GAMMA 2.783
L 3.781
W 4.597
X 5.173
因為路徑分了兩條,所以能帶圖中第二個高對稱點位置的標識為U|K。
用戶在Origin軟體中處理了BAND-REFORMATTED.dat後,需要按照KLABELS文件中的位置在能帶圖中標識高對稱點位置。下一版本的版本將會加入生成繪製能帶圖的python腳本的功能,需要用戶安裝了Matplotlib庫。下一版本的彩蛋:例子裡有一個band_plot.py腳本,只要執行這個python腳本,就能在文件夾下生成一個band.png,即為程序自動繪製的能帶圖,如下圖所示。
眾所周知,VASP計算的是在體系在0K下的電子能量,沒有考慮溫度的貢獻。吉布斯自由能的計算可以參考教程Learn VASP The Hard Way (Ex68) 頻率,零點能,吉布斯自由能的計算和Learn VASP The Hard Way (Ex69) 表面吸附物種熵的計算,系列教程可以在群217821116下載到PDF版本。本模塊的原理正是基於此教程,由許楠、李強和劉錦程貢獻。以/vaspkit.0.70/examples/thermo_correction為例。下圖展示的是在氧氣在磷摻雜的石墨烯上的解離吸附的過渡態(Carbon, 2016, 105:214-223.),圖片用VESTA軟體可視化。為了求得溫度的貢獻,需要進行振動分析(IBRION=5),為了節省計算資源,只有P、O原子放開,而C原子固定住,在計算Energy Profile時,需要對其他體系做振動分析,保持放開的原子一致。
+- Warm Tips -+
See An Example in vaspkit/examples/thermo_correction.
This Feature Was Contributed by Nan XU, Qiang LI and Jincheng LIU.
Only vibrations! No Translation & Rotation & Electron contribitions.
+---+
Please Enter The Temperature (K):
-->>
298.15
-->> (1) Reading OUTCAR File...
+- Summary ---+
H = E_DFT + E_ZPE + E_H + (nRT = n* 0.0256926 if in gas phase)
G = H - TS = E_DFT + E_ZPE + E_H - T*S
Temperature (K): 298.1
Entropy (eV/K): 0.0005
Entropy contribution T*S (eV): 0.1541
Enthalpy contribution E_H (eV): 0.0858
Zero-point energy E_ZPE (eV): 0.1944
Total energy at Zero K (eV) : -159.3511
Gibbs free energy at 298.1 K (eV) : -159.2250
+---+
+---+
| * DISCLAIMER * |
| CHECK Your Results for Consistency if Necessary |
| Bug Reports and Suggestions for Improvements Are Welcome |
| Citation of VASPKIT Is Not Mandatory BUT Would Be Appreciated |
+---+
298.15K下體系的吉布斯自由能由以下公式給定:
H = E_DFT + E_ZPE + E_H + (nRT = n* 0.0256926 如果體系處於氣相狀態還需加上nRT)
G = H - TS = E_DFT + E_ZPE + E_H - T*S
E_DFT為體系在0K下的靜態電子能量。本例中體系在298.15下的自由能為-159.2250 eV。要記住的是DFT計算中絕對能量沒有意義,只有相對能量才make sense。
現階段VASPKIT的手冊不是太完善,開發者正在努力開發。但是可以在以下QQ群裡(364586948,217821116)找到開發者交流,歡迎大家使用VASPKIT,並及時反饋BUG。VASPKIT的下載地址:http://vaspkit.sourceforge.net/?page_id=7 ,最新版會發布在以上兩個群裡和我的Github Repository (https://github.com/tamaswells/VASP_script)。引用參考: Wang, V. VASPKIT. "a Post-Processing Program for the VASP Code." (2013).