超軟贗勢(Ultrasoft Pseudopotentials)

超軟贋勢使用的平面波基底函數更少，特色是讓贋化波函數變得更平滑， Norm-Conserving 條件成為廣義 Norm-Conserving 條件。

以下就是Vanderbilt 的三階段介紹：

28

我們把它當作投影算符一般，便可以定義出一個non-local potential operator 如下：

| |

| 2.2.25

29

其實這就是 Kleinman-Bylander 型的贋勢，我們可以利用剛才創造出來的 non-local pseudopotential 來證明在 0~ 的情況下 是

的本徵態。

在同樣的參考能量下，含贋勢之新 Hamiltonian 向外積分，獲得的贋化 波函數，與真實位勢向外積分所產生的真實波函數，在 以外是一模一樣 的，這正是我們所希望看到的結果，由於這是標準的Norm-Conserving 位勢，

其散射性質具有波函數對數的能量導數在本徵值處之性質在贋勢與真實位 勢都是一致的。

第二個步驟，我們引進每個角動量通道的投影算符，這樣可以使用一個 以上之參考能量的技術。所以我們讓新的參考態在同一個角動量 ，給予兩 個以上的參考能量 ，並且要讓這些新的參考態滿足廣義 Norm-Conserving 條件 0如下：

2.2.26 接著利用新的參考態，可以定義一組新的投影算符如下：

| 2.2.27

其中 ，然後寫下新的 non-local pseudpotential 如下：

,

2.2.28

30

而non-local pseudopotential 也必須重新定義成

,

31

0 | | 2.2.33

利用這個關係可以證明真實波函數與贋化波函數之對數導數的一致性仍然 是具備的。

Vanderbilt 所提出的超軟贋勢，由於捨棄了 Norm-Conserving 條件，讓 原本即無結點的波函數也能被改造的非常平滑，這可以讓我們可以使用很 少的平面波基底展開，也因此在計算量上比使用Norm-Conserving 型的贋勢 更迅速且更划算。

圖6 實線是真實的徑向波函數，虛線則是超軟贋勢，可明顯看見核心部 分少了一大塊，變得非常的平滑。

資料來源：Vanderbilt PRB 1990 paper 2.2.5 廣義梯度近似法(GGA)

由於在LDA中，我們將 寫成 之函數，加上電子電荷密度梯度的修正，

此即所謂的廣義梯度近似法(generalized gradient approximation)，廣義梯度

32 hyperfine fields,in Kg)、磁矩(magnetic moment, in )…等，我們就必須考 慮自旋極化的部份。此時電荷密度 必須分解為兩部份，上旋密度

(spin-up density) 、下旋密度(spin-down density) 部份。

2.2.35

33

, , 2.2.38 交換相干位勢為， _, ^, ，經過一些推導可知總能量為：

E , 1

2

, ,

2.2.39

如果我們想要得知基態電荷密度與基態能量，也可以透過非自旋極化系 統的流程，解上述的式子，其中需要將電荷密度改為上旋密度和下旋密度 的相加，然後分解兩個自旋極化的Kohn-Sham方程式，而在Kohn-Sham方程 式中之Hamitonian只有交換相干位勢與上旋密度和下旋密度有關，至於在總 能量的表示式中，交換相干能為上旋密度和下旋密度的函數，而動能項雖 與自旋相關，但卻為未知函數，所以其後推導消去。

34

第三章 計算方法與結果

在這章我會介紹我整個系統的架構以及如何去計算系統的熱導(

Thermal

conductance)

，在3.1 將介紹如何去架構整系統的結構與大小，3.2 節我將說明

如何利用 Müller-Plathe 的方法得到系統的熱流及溫度分佈圖，並進一步計 算系統的熱導，3.3 節先介紹如何利用另一個分子動力學軟體(DL_POLY_4) 建構分子動力學（Lammps）所需的 force field，然後將 force field 導入分子 動力學（Lammps）3.4 節詳細介紹如何將密度泛函理論軟體(VASP)與分子 動力學（Lammps）結合，最後在 3.5 節開始分子動力學的模擬，並計算求 得我們想要的分子接面的熱傳導性質。

3.1 系統的架構

先利用密度泛函理論軟體（VASP）來確認 nano wire 與電極的接點位置，

接下來做吸附能以確定系統所需要的層數，以及利用分子動力學（Lammps）

計算 Z 方向層數大小與熱導率的關係、X-Y 截面積大小與熱導率的關係，

這樣一來系統的大小就可以確定了。

35

36

3.2 Müller-Plathe 的方法

熱導率是物質的一種基本傳輸性質，但熱導率已經被證明是最難計算的 傳輸系數之一，在這裡我們引用Müller-Plathe 的一篇文章[25]，文章中提到 使用非平衡分子動力學方法來計算的熱導率。在模擬這個實驗時，必須符 合以下這些條件：

一、 Homogeneity：

在模擬的系統中，不可以有外加的邊界，因為這樣將會讓這 個結構與bulk 不同。

二、 Boundary conditions：

必須選擇一個適合的邊界條件，通常要選擇具有週期性的邊 界條件。但是其實我們的系統不是週期性的，那是因為使用密度 泛函理論軟體(VASP)與分子動力學（Lammps），這些現有的軟 體在編寫程式的時候必須寫成週期性的邊界條件。

三、 Small temperature gradient：

如果系統中的 temperature 和 density 變化太過劇烈，那將會 使我們在不同的 temperature 和 density 所計算出來的 thermal conductivity 相差不大，如此系統將不再出現於線性區間。

37

四、 Conservation of energy and momentum：

如果能量守恆，整個系統就不會有定溫的問題。

五、 Convergence：

在波動沒有很大的情況下，有意義的統計必須合理的獲得。

六、 Hamiltonian：

必須在適當的Hamiltonian 下獲得運動方程式，如此一來與統 計力學做連接會比較容易。

在巨觀上，熱通量等於負熱導係數乘以溫度梯度，而且熱通量方向是垂 直於表面的。一般情況下，熱導係數是一個3x3 的張量。下面我們就以 z 方向為溫度梯度的方向，我們可以定義熱傳導公式如下：

3.2.1

⁄lim ∞lim

∞ ⁄ 3.2.2 (3.2.1)式中， 代表熱通量、 代表溫度梯度、 代表熱導係數(率)。

現在我們將系統視為一個模擬的盒子，並且將盒子重複，以便加強統計，

並且讓這個盒子是一個對稱的系統，然後將盒子分為 個厚度相等並且垂直 於 方向的板子，如圖10。

38

圖10 假設 Slab0 是冷端板，而 SlabN/2 是熱端板。動能從冷端板轉移至 熱端板，然後再透過熱傳導由熱端板回流至冷端板。而溫度分佈的計算方 法決定於Slab1 至 SlabN/2-1 以及 SlabN/2+1 至 SlabN-1，slabs 的厚度、體 積都是相等的。

資料來源： Florian Müller-Plathe,J. Chem. Phys. ,Vol. 106, No.14, 8 April 1997.

每一片slab k 的瞬時運動的溫度 可以用下式計算求得：

1

3 ² 3.2.3 (3.2.3)式中，每個 slab 裡包含 i 顆原子，並且每顆原子有相對的質量 與速 度 ，而 為波茲曼常數。然後計算平均時間就可以得到溫度分佈。

每經過t 步我們交換冷端板中最熱的原子與熱端板中最冷的原子一次，

如此一來便可以產生的熱通量。相對的，熱端板和冷端板之間將產生溫度 梯度。有了以上兩個模擬數據後，便可以計算出熱導係數，公式如下：

39

∑ 2

^{2 2}

2 ⁄

3.2.4 (3.2.4)式中， 是系統 x 方向的大小， 是系統 y 方向的大小。

在溫度梯度小的條件下，我們可以計算出有效的熱量轉換，並可以計算 出Slab0 與 SlabN/2 之間的熱傳導系數。(3.2.4)式分母中的參數 2，是因為 系統有週期性的邊界，所以除以2 便可以在系統中求得我們想要的部份。

3.3 DL_POLY_4 建構之 force field 對應至 Lammps

由於Lammps 所提供的勢能函數較少，所以利用 DL_POLY_4 來建構系 統所需要的force field，然後將這些 force field 導入 Lammps 來做分子動力 學的計算，在此我將DL_POLY_4 與 Lammps 之 force field 對應及公式詳細 列出[26][27]：

一、 Bond Potentials：

圖11 原子間的鍵結向量。

Bond potentials 是用來描述原子之間的化學鍵，在 DL_POLY_4 我們 使用Harmonic bond，如下：

1

2 3.2.5

40

對應到Lammps 我們使用 harmonic，如下：

3.2.6 r = distance (computed by LAMMPS)

coeff1 = K (energy/distance^2) (the usual 1/2 is included in the K) coeff2 = r0 (distance)

二、 Valence Angle Potentials：

圖12 The valence angle 和相關向量。

Valence Angle Potentials 用來描述鍵結之間的夾角，在 DL_POLY_4 我們使用Harmonic cosine，如下：

2 cos cos 3.2.7 對應到Lammps 我們使用 harmonic，如下：

3.2.8 theta = radians (computed by LAMMPS)

41

coeff1 = K (energy/radian^2) (the usual 1/2 is included in the K) coeff2 = theta0 (degrees) (converted to radians within LAMMPS) 三、 Dihedral Angle Potentials：

圖13 The dihedral angle 和相關向量。

Dihedral Angle Potentials 用來描述二面角潛力相互作用而產生的扭 轉力量，在DL_POLY_4 我們使用 Cosine potential，如下：

1 cos 3.2.9

對應到Lammps 我們使用 harmonic，如下：

1 3.2.10 phi = radians (computed by LAMMPS)

coeff1 = K (energy)

coeff2 = d (always +1 or -1) coeff3 = n (1,2,3,4,6)

42

四、 Inversion Angle Potentials：

圖14 The inversion angle and 和相關向量。

Inversion Angle Potentials 用來描述三個原子圍繞一個中心原子的潛 力互動所產生的特定的幾何。

在DL_POLY_4 我們使用 Planar potential，如下：

1 cos 3.2.11

對應到Lammps 我們使用 cvff，如下：

1 cos 3.2.12 coeff1 =K (energy)

coeff2 =d (+1 or -1) coeff3 =n (0,1,2,3,4,6)

五、 Short Ranged (van der Waals) Potentials：

在DL_POLY_4 我們使用 Lennard-Jones potential，如下：

43

4 3.2.13

對應到Lammps 我們使用 lj/cutoff，如下：

4 3.2.14 r = distance (computed by LAMMPS)

coeff1 = epsilon (energy) coeff2 = sigma (distance) 六、 Metal Potentials

在Lammps 有提供金的勢能函數 EAM，所以我們無需使用

DL_POLY_4 來提供 force field。EAM 勢能函數如下：

1

2 3.2.15 3.4 結合密度泛函理論軟體(VASP)與分子動力學（Lammps）

我們編寫了程式用於計算系統的熱傳導性質。首先，密度泛函理論軟體 (VASP)是從第一原理出發，在此我們是用於計算原子間的力。其次，將力 量導入分子動力學（Lammps）的模擬，此模擬採用積分牛頓方程。最後，

使用 Müller-Plathe 的方法去計算獲得分子接面的熱傳導性質。該流程圖的 計算過程描述如下：

44

圖15 計算獲得分子接面的熱傳導性質之流程圖。

在圖15 中，分成五個步驟第一個步驟是先執行分子動力學（Lammps）

的模擬，在此必須先輸入系統的座標，並在這一步指定系統的溫度及各項 參數的設定，當中會啟動第二個步驟及第三個步驟，將座標檔轉換為密度 泛函理論軟體(VASP)的格式，並且執行密度泛函理論軟體(VASP)，第四個 步驟及第五個步驟為讀取密度泛函理論軟體(VASP)計算完單點能(single point energy)後的結果，結果將會寫在 OUTCAR 這個檔案裡面，利用我們 寫的程式會自動去讀取OUTCAR 中每顆原子的力並且導入分子動力學

（Lammps）的繼續模擬下一個迴圈，如此一來便可以計算出分子接面的熱 傳導性質。

45

3.5 以分子動力學方式計算分子接面的熱傳導性質之模擬結果

由於從第一原理出發的密度泛函理論軟體(VASP)，其計算量非常的龐 大，所以我們採取先執行分子動力學（Lammps）的模擬，達到平衡後再將 密度泛函理論軟體(VASP)引進來計算原子間的力。

我們有兩個系統，分別是金電極接一個苯環與金電極接兩個苯環，每 個系統分別模擬14 個溫度(30K、50 K、70 K、90 K、110 K、130 K、150 K、

170 K、190 K、210 K、230 K、250 K、270 K、300 K)，使用 Müller-Plathe 的方法去計算獲得分子接面的熱傳導性質時，有幾個參數需要設定，我們 將系統分成六個等厚的slab，並且每次交換 20 個粒子，而模擬的總時間固 定為1 10 秒，但分別計算兩個不同的時間步與模擬圈數，如下表 1：

總時間 時間步 模擬圈數

一個苯環 1 10 秒 0.001 皮秒 10000 0.0005 皮秒 20000 兩個苯環 1 10 秒 0.001 皮秒 10000 0.0005 皮秒 20000 表1 系統在相同的總時間下，模擬不同時間步與模擬圈數之參數。

46

一、 一個苯環系統在不同時間步與模擬圈數下之熱導與誤差：

一、 一個苯環系統在不同時間步與模擬圈數下之熱導與誤差：

在文檔中 以分子動力學方式研究分子接面的熱傳導性質 (頁 37-0)