第二章、
理論原理與計算方法2-1
密度泛函理論-(Density functional theory，DFT)此為一種用於研究多電子系統之電子結構的量子力學方法，近來多應用於多電子計算 表達出來，以下將針對 Kohn-Sham equation 及 Exchange-Correlation Functions 做一些簡單 的介紹。
2-1-1 The Kohn-Sham equation
Kohn 與 Sham 於 1965 年提出的 Kohn-Sham 系統為現今密度泛函理論中最為普遍之應
7
基態的電子密度及能量可藉由 Schroedinger-like non-interacting single-particle equations 的 計算而得到，另一方面 Hartree operator，VH 和 exchange-correlation operator，Vxc，並根據 電子密度ρ(r)解出上述方程式即為自洽場問題。欲求出方程式則據下圖中之反覆運算過程 為一必要過程。而所謂的自洽過程即是一開始先猜想一密度ρn-1 ，一 Hamiltonian，則得
到 HKSn-1，接著將求得之 eigenvalue 建立一新的密度，再由此新的密度創出新的 Hamiltonian，
HKSn，直至新舊密度彼此之間的值趨於相等。此想法主要在於如何將一多電子系統轉變成
一個獨立之大電子，如此將少了電子彼此之間的作用力，大大降低了計算的困難度，即是 利用唯一的決定因素-基態電子密度，來了解一個多電子系統的物理量。
圖 2-1-1 自洽過程示意圖
8
2-1-2
交換相關函數(Exchange-Correlation Functions)量子力學電子中彼此之間的 exchange energy 與 correlation energy 稱之為交換相關函數 (exchange-correlation functions)Vxc，因我們無法得出真正之數學型態，故用此近似法。而不 同之 DFT 方法，其差異即在於用的交換相關函數不同；其中若是只和密度函數本身有關之 方法稱之 Local (Spin) Density Approximation (LDA or LSDA)，而與密度函數梯度有關之近 似則稱之為 Generalized Gradient Approximation(GGA)。亦有將 exchange energy 和
correlation energy 先分開討論後再進行組合的方式如 BLYP-於 1988 年由 Beck’s exchange function(B88) 與 Lee，Yang and Parr’s correlation function and Parr’s correlation function (LYP) 組合得之；而 1991 年的 Perdew 和 Wang’s correlation (PW91)也是類似的例子。近期 對於 DFT 及 wave function 結合之研究亦有得到準確且有效率的理論方法。
2-1-3
基底(Basis set)於量子化學計算中唯一有限集合之函數組，一般用於建分子軌域或用來平均電子密度。
DFT 方法中使用的為適度基底函數組，相較於傳統的 ab initio 方法教不嚴謹苛求，如此較 適用於處理較大的系統。DFT 理論(或 ab initio 理論裡的分子軌域)中的準粒子(quasi particle) 波函數是以如下列所示之基底函數組來建構:
2-1-4
布洛赫定理與平面波基組 (Bloch theorem and plane wave basis sets)布洛赫定理(Bloch theorem)是研究完整晶體能帶理論及計算的基礎，用於描述含有週
9
φ(r)=∑Ckn，ke−ik∙r𝑢(r)
其中 k 是獨立的，每一個特徵函數皆有相同的 k 但卻有不同的 n，其皆被表示於此基底中，
假設特徵函數的 k 不同，則具不同 k 值之新基底，將用於尋找係數 Ckn，k。實際上，無限的 基底必須被限制於某處。如:波函數(plane wave)可輕易將 K 限制於 K ≦ Kmax，此可視為一
個球型在座標空間中離球心的半徑為 Kmax，所有的晶格向量在此球型內部皆必須考量此一
基底。現在大部分市售的軟體，內含的相關變數 cut-off energy(Ecut)自由能等同於 Kmax， 如下所示；
Εcut= ℏ2𝐾2𝑚𝑚𝑎𝑥2
2-1-5
有效核勢(Pseudopotential or Effective core potential，ECP) pseudo potential 若要計算一包含重原子或是過度金屬的粒子時，於硬體方面的計算需求較高，如 CPU 速度、硬體空間和記憶體等等，而有效核勢的加入企圖藉由取代原子中核的價電子、有效 核電位，以薛丁格方程式取代核電子中的庫倫力，以降低這些硬體的計算需求量。其中，偽電子及全電子價態須選一截止核半徑 rc (core cut-off radius)，且具有相同的能量態及振 幅。
圖 2-1-2 核之庫倫力波函數及 pseudo-wave function 彼此關係示意圖
核的庫倫力
pseudopotential
r
c= Cut off radius
10 中使用的伺服器平台為 IBM Cluster1350。IBM Cluster 1350 系統為 x86 架構之分散式平行 電腦系統，當中建置了 512 個計算節點，共 2048 個計算核心及 GPFS 檔案儲存空間，採 Intel Woodcrest 3.0 GHz 雙核心處理器，總計算效能(Rmax)為 19.91TF，理論值(Rpeak)為 24.86。
2-2-2
操作軟體 VASP(Vienna Ab-initio simulation Package)VASP 是一種由維也納大學 Hafner 小組研發出用於進行電子結構計算與量子力學-分子
Package(VASP)，以 Perdew-Wang 1991(PW91)的 Generalized Gradient Approximation (GGA) 作為本研究之交換關聯函數(Exchange correlation function)，同時亦使用了由 Blochl 等多位 作者發展出的投影綴加波(Projector Augmented Wave，PAW)的方法作為系統分子之理論計 算基礎。利用 VASP 進行計算流程時，須先建立一檔案資料夾其中主要包含四個主要 INPUT 檔案:INCAR、POTCAR、KPOINT、POSCAR，將所需的參數設定完成後即可進行計算程 序。以下將對這四個主要的 INPUT 檔案做設定及原理的介紹:
11
12
KPOINT:主要表示計算波函數時，取點做計算的數量，設定如下:
Automatic Mesh 0 Monkhorst-Pack 2 2 1 0 0 0
以上的參數為本實驗主要使用的設定，其中 2x2x1 分別表示 X、Y、Z 三軸方向所取 點的數目，此數目需搭配 POSCAR 中 3x3 的矩陣，舉例來說，當一軸的長度為 10Å，另一 軸長為 20Å，則取點的數目為 2 和 1；因為將這兩軸寫成反函數，分別為 1/10 及 1/20，搭 配取點為在第一條軸上每 1/10 處取兩點做計算，而另一軸則是每 1/20 處取一點做計算，
又於第一條軸上的兩點等距故此兩軸所取點的數目相同。下圖可簡單解釋各軸上取點與波 函數的關係。
圖 2-2-1 波函數取 1、2、4 點做計算之圖示
從圖中可知，當欲描述一波函數時，若僅取一點來計算或敘述事相當不準確的，而當 取點越多，其將波函數分割的越精細則其計算將越趨近於真實的波函數數值及圖形。
13 0.0000000000000000 11.2854799999999997 0.0000000000000000 0.0000000000000000 0.0000000000000000 20.7327499999999993
此方格中的矩陣為晶格常數，分別為 X、Y、Z 三軸；用於將 direct 座標轉換為 cartesian 座標 研究中建立的 slab model，是依據 crystal model 而得出的，crystal 是一個 X、Y、Z 軸 無限週期性延伸的金屬，於層層之間加入一層所謂的真空層(vacuum space)，即可免於受到 層與層之間的交互作用力，而加此真空層一是在做模擬時分子於金屬表面上活動及反應的 空間；在本次計算中，我們加入的真空層約厚為 10Å。
14
POTCAR:主要為提供計算工作中所包含的元素種類，其中包含了計算工作所需每個元 素的 pseudopotential(模擬勢能)，元素排列方式須與 POSCAR 中的相同。下表為 POTCAR 的設定範例:
PAW_GGA Pt 05Jan2001
--- RWIGS = 2.750; RWIGS = 1.455 wigner-seitz radius (au A)
ENMAX = 230.277; ENMIN = 172.708 eV ---
PAW_GGA Au 18Jul2000 11.0000000000000000
--- 以上為四個主要 INPUT(INCAR、KPOINT、POSCAR、POTCAR)檔案的設定介紹。將
設定完的 INPUT 檔上傳至國家高速網路中心，以 8-32CPU 做計算，隨著結構位置的不同、
吸附物種類的不同及欲求出的不同結果，其所需要的計算時間也不同，一般如果要求反應 物及產物最穩定的能量(Local minimum)，其計算量約需要 10-20 小時可完成；若是過度態 (Transition state)的運算，則其需要約 2-5 天的計算量。當計算工作完成後，一般會得到幾 個 OUTPUT 檔，包含 OUTCAR、CONTCAR、CHGCAR、WAVECAR、DOSCAR……等 等檔案；我們會從 OUTCAR 檔案中取得能量及原 direct 座標經矩陣轉換後得的 cartesian 座標，CONTCAR 中可得到最後結構收斂的 direct 座標，其中本次實驗計算中我們也使用 了 DOSCAR 這個檔案來分析繪製吸附物的 DOS(density of state)，藉此了解電子於分子軌 域中如何轉移。
15
OUTCAR:包含一開始設定的原始資料及每一大迴圈收斂後得到的能量值及結構的 cartesian 座標&其分子互相影響產生的斥力。以下為 OUTCAR 的範例介紹:
POSCAR:名稱 positions in direct lattice
velocities in cartesian coordinates
ion position nearest neighbor table Maximum number of real-space cells 4x 3x 2 Maximum number of reciprocal cells 2x 2x 4
FEWALD: VPU time 0.06: CPU time 0.06
之後會開始進行重複的迴圈運算，從大迴圈開始跑之後換小迴圈，直到最後達到所設定的 準確度，即會停止能量及結構上的收斂。
--- Iteration1( 1) --- Free energy of the ion-electron system (eV)
---
16
DOSCAR:包含各原子的密度態；能量與 s、p、d 軌域之間的的分布，此檔案用於分析 繪製 DOS 的圖。
DOSCAR 第一部分為所有參與計算的原子總和，再來是各個原子能量及 s、p、d 軌域的分布能量態；若有考慮磁性(即計算有磁性的金屬 Fe、Co、Ni，或是吸附物 O 則需加入 ISPIN=2 的參數)進行計算，則其會分裂成 s-up、s-down、p-up、p-down、d-up、
d-down 軌域的分布 88 88 1 0
0.2598632E+02 0.9773510E-09 0.1128548E-08 0.2073275E-08 0.5000000E-15 1.000000000000000E-004
CAR
名稱
7.93249810 -18.79229262 301 2.43130176 1.00000000 -18.792 0.0000E+00 0.0000E+00
-18.703 0.0000E+00 0.0000E+00 -18.614 0.0000E+00 0.0000E+00 -18.525 0.0000E+00 0.0000E+00 -18.436 0.4198E-02 0.3740E-03
-18.792 0.0000E+00 0.0000E+00 0.0000E+00 -18.703 0.0000E+00 0.0000E+00 0.0000E+00 -18.614 0.0000E+00 0.0000E+00 0.0000E+00
能量 s 軌域 p 軌域 d 軌域
17
2-3
計算流程圖 2-3-1 理論計算流程 model
Eads
ΔE & Ea
依據催化劑所使用的金屬種類，按其最穩定的排列方式建立金屬 表面(111、100…)，接著進行結構的優化，得到我們的 surface
將反應的反應物、中間物及生成物，置於表面模型上可能的吸附 位置，找出最穩定之吸附位置並得其吸附能
利用上個步驟中所求得反應物及生成物之穩定優化結構計算在 反應中中所需的反應能量(ΔE)和活化能(Ea)
18 型；平面波(plane wave)基礎設定中的 cutoff energy 為 600 eV；Brillouin-zone 積分設定，
利用 Monkhorst-Pack scheme 做取樣；KPOINT 設定為 2x2x1。而表面的選擇為相對較 穩定之 fcc 最密堆積(111)平面，並於 XYZ 軸無限週期性延伸的 super cell 中，單位晶 格(unit cell)為五層的 4x4 層板(slab)表面，並於 super cell 中加入一約 10 Å 的真空層，
使吸附的分子能於此空間反應；此五層單位晶格中，我們將固定最下面兩層(freeze)，
將上三層金屬設為可自由移動(relaxed)，隨反應收斂至最穩定的結構。本篇實驗計算 結構優化製能量最低點所使用的方式為 quasi-Newton method，能量之收斂極限設為 1x10-4 eV，而計算過度態所使用的方法為 Nudged Elastic Band(NEB) method，其能量 收斂極限設為 1x10-2 eV，此法主要是從反應物及產物之間的局部最小值中(local 另一種 surface(0-12-8)則是模擬實驗上 core-shell 的結構，亦是呈現一個以 Pt 為主體，並於 二、三層中混合其他金屬的表面結構，本篇計算所模擬的為一共 80 顆金屬的結構，其中 取代的金屬共 20 顆，符合中央王冠文老師實驗室所合成的 25%取代 Pt 金屬的結構。如下 圖所示的表面結構圖:
19
二元合金(包括 PtAu(111)、PtAg(111)、PtPd(111)、PtRh(111)、PtRu(111)、PtSn(111))
二元合金以 PtAu(111)為例:
Surface (4-8-8
a. 第一層取代四個 Pt
俯視圖 b. 第二層取代八個 Pt
俯視圖
One layer
Two layer
Three layer
Four & Five layer
T
H B
20
Surface (0-12-8) a. 第一層未取代
b. 第二層取代十二個 Pt
俯視圖
One layer
Two layer
Three layer
Four & Five layer
H B
T
21
三元合金(PtAgSn(111))
Surface (4-8-8)
a. 第一層取代四個 Pt
俯視圖 b. 第二層取代八個 Pt
俯視圖 One layer
Two layer
Three layer
Four & Five layer
22
Surface (0-12-8) a. 第一層未取代
b. 第二層取代十二個 Pt 圖示
俯視圖 One layer
Two layer
Three layer
Four & Five layer
23