ρn-1 = ρ
n
14
交換相關函數(Exchange-Correlation Functions)
討論量子力學中的電子間的 exchange 及 correlation energy 稱為交換相關函數 ，由 於我們無法得到真正的數學型態因此使用近似法，不同的 DFT 方法差異就在使用不同的交 換相關函數，如果只跟密度函數本身有關稱為 Local (Spin) Density Approximation (LDA or LSDA)，另一與密度函數梯度有關的近似方法為 Generalized Gradient Approximation(GGA)。
[127]也有將 exchange 及 correlation energy 兩者分開討論再進行組合的，1988 年的 Beck’s
exchange function(B88)與 Lee, Yang and Parr’s correlation function (LYP)組合稱為 BLYP 以 及 1991 年的 Perdew 和 Wang’s correlation(PW91)。近期也研究出 DFT 與 wave function 結 合獲得準確並有效率的理論方法。
Vˆ
xc15
iC
i2-2 基底(Basis set)
基底(Basis set)為一個有限集合的函數組，用來建構分子軌域或是來帄均電子的密度。
DFT方法中使用的適度基底函數組相較於傳統的ab initio方法不苛求，適用於處理較大的系 統。DFT理論(或ab initio理論裡的分子軌域)中的準粒子(quasi particle)波函數是用下列這種 基底函數組來建構，表示方式如下：
(6) 其中的 Ci 是線性組合(linear combination )的係數，各種不同的基底是為了滿足實際上的需 求，主要包含兩種分類：第一種是以原子為中心(atom-centered basis set)的基底，此類含是 主要運用在局部性的系統中；第二種是帄面波基底 (plane-wave basis set)，主要運用在有週 期性的系統。本次使用的軟體為 VASP(Vienna Ab-initio simulation Package，於下節中詳細 介紹)，利用 VASP code 來作帄板模型計算都是利用帄面波基底，因此，以下針對帄面波基 組作探討。
布洛赫定理與帄面波基組 (Bloch theorem and plane wave basis sets)
Bloch theorem 是 完 整 晶 體 能 帶 理 論 與 計 算 的 基 礎 ， 其 應 用 在 含 有 週 期 能 量 的 Hamiltonian 式子中(V(r+R) = V(r) ) 例如：固態晶體。R 代表周期原子核的座標， r 代表電 子的座標，下列公式中，Hamiltonian 的 eigenvalue 也是周期性且可寫成 ug(r)，且含有晶格 的周期性和 planewave，此公式為熟知的 Bloch thermo:
16
因此，Bloch’s therom 可被表示為較普遍的形式：任何 eigenfunction 可被表示為一個方
程式，uk(r)此含有晶格的周期性和 plane wave eikr此處 k 為 Brillouin zone 的任何向量。
(10) 週期 Hamiltonian 的任何 eigenfunction 可藉由係數的無限設定帄均表示：
(11)
17
於 Kmax ，例子如下： (12)
有效核勢( Effective core potential，ECP )
波函數最重要擺盪部分為從外部接近內核區域的末端，無論如何，固態的區域被很多 原子外層區域所保護住，且電子從自由原子電子中不會不一致的活動。此外，計算一個包 含重原子或是過渡金屬的例子時，硬體來源方面有著相當高的計算需求，例如：CPU 速度、
硬體空間和記憶體。可利用 effect core potential(ECP)取代接近內核區域的勢能，可大量減 低這些需求，此為被物理學家所熟知的 pseudopotential。ECP 被設計為產生較帄坦的波函 數末端，在此 core 區域中沒有結點。ECP 連續性的包含著真實的勢能，此勢能在價區域中 行動如所有真實的電子波。
在特殊的元素沒有特定的 ECP 建構的方式。可能含有無限的選擇，ECP 用於計算幾何 結構和包含重元素的分子能量的工作已經可做得很好且以文件化。ECP 基底連結 DFT 可以 操作出相當正確的結果，所消耗的計算花費也減少，因此被應用於大部分的計算研究。
m E
cutK
2
2 max
2
18
2-3 系統與軟體
國家高速網路與計算中心
在國家高速網路與計算中心創建一個帳號，將已設定好的工作在電腦上上傳至伺服器計 算，使用的伺服器帄台為 IBM cluster1350(圖 2-2)，IBM Cluster 1350，屬於分散式帄行電 腦系統，採 Intel Woodcrest 3.0 GHz 雙核心處理器，擁有 512 個節點(node)、共 2048 個 CPUs，
總計算效能（Rmax）為 19.91TF，理論值（Rpeak）為 24.86。
圖 2-2 國家高速網路與計算中心所提供的 IBM cluster 1350。
操作軟體 VASP(Vienna Ab-initio simulation Package)
本實驗的計算是使用 Vienna Ab-initio simulation Package(VASP) [128]計算軟體，以 Perdew-Wang 1991 (PW91)的 Generalized Gradient Approximation (GGA) [127]作為本實驗系 統的交換關聯函數(Exchange correlation function)，同時也使用 Blochl 等作者發展的投影綴 加波(Projector Augmented Wave；PAW)[129]的方法，作為系統分子的計算理論基礎。VASP 廣泛運用在下列的應用中：
19
1. 採用週期性邊界條件(或超原胞模型)處理原子、分子、表面體系及固體。
2. 計算材料的結構參數(鍵長、鍵角、晶格常數及原子位置等)和構型。
3. 計算材料的電子結構(能級、電荷密度分佈、能帶、電子態密度)。
利用 VASP 計算得先建立一個檔案資料夾其內部需包含四個主要檔案： INCAR POTCAR、KPOINTS 和 POSCAR，將所需參數設定好後即可進行計算工作，以下將針對 此四個檔案的設定及原理做介紹。
(一) INCAR
SYSTEM = Bulk CeO2
ISMEAR = -1 ! -5 for insulator JCP vol(114) p5816 SIGMA = 0.1
EDIFF = 1E-4 NSW = 300 ENCUT = 600 ALGO = FAST ISPIN = 2 ISIF = 2 GGA = 91 IBRION = 2 IALGO = 48
!NBAND = 1000 NPAR = 16 NELM = 300 LREAL = Auto
ISYM = 0 ! default for PAW LWAVE = F ! write WAVECAR
20
INCAR 為計算時的操作設定，為 DFT 裡所提及的基底設定，下列為常用參數的介紹：
NSW = 300 代表在限度內整個外迴圈的結構收斂次數。
NELM =300 代表內迴圈電子自洽場(self-consistency)收斂次數。
EDIFF = 1E-4 電子自洽場迴圈總能量的允許誤差範圍。
GGA = 91 表示使用 Perdew -Wang 91 此種近似方式計算法則。
ENCUT = 600 為 Cutoff Energy，越大值計算準確度高但耗時久。
ISIF = 2 or 3 此設定為 2 時，則固定晶格常數計算，稱為 fix；相反為 3 時 則自由優化，稱為 relax。
ISPIN = 2 計算工作中的自旋極化(spin-polarized)，對 Co 及 Ni 含有鐵磁 性性質金屬所需要具備的參數，其他金屬可不用此設定。
IBRION = -1 or 2
離子如何移動及收斂，設定為 2 表示 conjugate-gradient，用於 結構收斂。優化結構進行 DOS 計算時，設定為-1 表示原子不 動，也可設定 NSW=0。
RWIGS= 1.201 0.802
輸入元素的 wigner-seitz radius。當需要做 DOS 計算時，會加 入此參數， RWIGS 值可取自 POTCAR，並依照元素在 POTCAR 的順序填入。
NPAR =1 or 16 帄行計算的數量。做 DOS 計算時，將此參數設定為 1。
21
(二) KPOINT Automatic mesh 0
Monkhorst-pack 2 2 1
0 0 0
KPOINTS 主要代表計算波函數時取點做計算的數量，研究所使用的為 2 x 2 x 1 的 KPOINTS，
分別代表 X、Y、Z 軸方向分別所取的點數，頇搭配 POSCAR 中 3 X 3 的矩陣來觀察，舉 例來說，一軸的長度為 10Å ，另一軸為 20Å ，所取的點分別為 2 跟 1，因為將這兩軸寫成 反函數，分別為 1/10 及 1/20，搭配取點為在第一條軸上每 1/10 份取兩點做計算，另一條 軸則是每 1/20 份取一點做計算，但第一條軸因為在每 1/10 處的兩點是等距，所以跟另一 條軸的點數為相同，而各軸上的取點跟波函數的關係為何？利用圖 2-3 波函數圖加以解 釋。
圖 2-3 波函數從上而下分別取 1、2 和 4 點去計算波函數的值。
22
10.8575653100000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 9.4410598310000005 0.0000000000000000 0.0000000000000000 0.0000000000000000 25.9663481300000001 96 2 2
23
原子種類的排序頇和 POTCAR 中的順序相同，在這邊的順序為 96 個 Ru、2 個 C 及 4 個氫，
第 8 行中，direct 為下列座標的種類，若為 cartesian 則表示下列座標為已經過矩陣轉換。
第 9 行以下則是 modle 中所有原子的座標，座標順序和第 6 行相同，前 96 個為 Ru 的座標，
接著 2 個 C 及 4 個氫的座標。座標中，第 1、2 和 3 列分別為 X、Y 和 Z 軸上的座標，第 4、
5 和 6 列中的 F 和 T 分別代表 X、Y 和 Z 軸上的座標在計算時固定或可以移動。
研究中所建立的 slab model，是根據 crystal model 而來的，crystal 是一個 X、Y、Z 軸 無限週期性延伸的金屬，在層與層間加入了真空層(vacuum space)，就可將層與層間相隔，
而真空層也是在做模擬工作時，分子在金屬表面上活動及反應的空間(圖 2-4，真空層厚度 約為 10Å ，為避免吸附分子與上一層表面有作用，吸附物至上一層表面至少要有 5 Å 的真 空層。
圖 2-4 金屬晶體中加入真空層(vacuum space)形成 slab 金屬表面。
24
(四) POTCAR
PAW_GGA Ru 03Mar1998 8.00000000000000000 parameters from PSCTR are:
VRHFIN =Ru: s1 d7 parameters from PSCTR are:
VRHFIN =C: s2p2
25
POTCAR 內容主要是提供計算的工作中所包含的元素種類，並提供工作中所需每個元
素的 pseudopotential(模擬勢能)，元素編排順序需和 POSCAR 中相同。而 POTCAR 中 pseudopotential 的模擬方式可參考圖 2-5：
圖 2-5 全電子波函數以及 pseudopotential 的關係圖。
一個原子中的內核電子並不會受化學環境的改變而受影響，所以當鍵結時的貢獻並不會有 很大的改變，主要的能量差異來自於外層電子部分，從圖 2-5 可得知，rc 左半邊為內核電 子區，黑色實線部分為真實的函數部分，在模擬(紅色虛線部分)的過程中，由於內核電子 的差異對於能量上貢獻度不大，所以可以用較簡單的方式模擬以減少工作量，到達 rc 後為
-Z/r
r
cpseudo po te n ti al
Ψ
pseudor
(a.u.)Vps
Ψr
26
結構
外層電子的部分，就可模擬的與真實函數一樣的帄滑曲線。
綜合以上四個檔案(INCAR、KPOINTS、POSCAR、POTCAR)的設定後，丟入高速電 腦中以 8~16CPU 去計算，以反應物及產物的 Local minimum 的計算量約 10 小時到 20 小時 間完成，若是計算過渡態(transition state)則約兩天到一個禮拜間。計算完成，所得到的檔 案包含 OUTCAR、OSZICAR、CONTCAR、CHGCAR、WAVECAR、DOSCAR、CHG 等 的檔案，可從 OUTCAR 中的資料中取得能量及經過轉換 cartesian 座標；CONTCAR 中可 得到最後結構的 direct 收斂座標；另外本次工作也利用了 DOSCAR 來繪製吸附物的 density
of state，以及結合 CHGCAR 和 bader charge 來得到分子與表面間的感應電荷以及繪製圖形，
以下簡略介紹：
(五) DOSCAR
100 100 1 0
0.2661730E+02 0.1085757E-08 0.9441060E-09 0.2596635E-08 0.5000000E-15 1.000000000000000E-004
CAR
Bulk CeO2
5.20225594 -14.25335317 301 2.00872117 1.00000000 -14.253 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
-14.189 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 -14.124 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 ………..
………..
-14.253 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 -14.189 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
-14.124 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
27
DOSCAR 可來繪製吸附物的 density of state(簡稱 DOS)，DOS 計算時，INCAR 的 設定修改為固定原子的位置將 IBRION 設為-1 或是 NSW 設為 0，且將 NPAR 設為 1，
只利用一個 node 去 run，並使用 RWIGS 半徑去侷限範圍，讓已優化的結構進行電荷 分布的計算。
DOSCAR 中第一部分為所有參與計算中 100 個原子總和，共有五列，第一列為能 量，第二列為 DOS 的累計值，第三、四和五列分別為 s、p 和 d 軌域的分布。其他以 下的部分為 100 個原子的各別原子密度態，共有四列，第一列為能量，第二、三和四 列分別為 s、p 和 d 軌域的分布；若考慮磁性(計算有磁性的金屬 Fe、Co 或 Ni，或是吸 附物 O 需加入 ISPIN=2 的參數)進行計算，共有七列，第一列為能量，第二-七列分別 為 s-up、s-down、p-up、p-down、d-up 和 d-down 軌域的分布。
利用 VASP 指令將 DOSCAR 分割成 101 個檔案，其中包括整體 DOS 以及各個原 子的 DOS，在選擇相關原子的 DOS，將軌域分布對能量作圖(如圖 2-6)進行分析。
圖 2-6 在 Ru 表面上吸附物 CHCH 上 C 的 DOS。
0.0 0.1 0.2 0.3 0.4
-20 -10 0 10
Density
Energy
28
(六) OUTCAR
vasp.4.6.28 25Jul05 complex
executed on LinuxIFC date 2012.10.22 11:34:45
29
30
--- Iteration 1( 1) ---
……….
--- Iteration 2( 1) ---
……….
reached required accuracy - stopping structural energy minimisation