第三章 計算方法與結果
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
一、 一個苯環系統在不同時間步與模擬圈數下之熱導與誤差：
Temperature(K) Thermal Conductance(W/K) Error
30 1.02E-11 3.08E-13
50 1.35E-11 4.73E-13
70 1.11E-11 3.95E-13
90 1.13E-11 6.14E-13
110 1.21E-11 4.48E-13
130 1.16E-11 5.44E-13
150 1.30E-11 6.11E-13
170 1.37E-11 4.25E-13
190 1.38E-11 6.87E-13
210 1.81E-11 8.45E-13
230 2.25E-11 8.03E-13
250 2.32E-11 1.15E-12
270 2.74E-11 9.76E-13
300 2.78E-11 1.69E-12
表2 系統在時間步為0.001 皮秒，模擬圈數為 10000 之熱導與誤差。
Temperature(K) Thermal Conductance(W/K) Error
30 9.49E-12 3.39E-13
50 9.70E-12 3.35E-13
70 1.02E-11 4.36E-13
90 1.09E-11 5.20E-13
110 1.26E-11 5.27E-13
130 1.14E-11 5.57E-13
150 1.25E-11 3.86E-13
170 1.34E-11 5.19E-13
190 1.38E-11 5.31E-13
210 1.72E-11 7.21E-13
230 1.79E-11 1.06E-12
250 2.40E-11 9.37E-13
270 2.49E-11 8.93E-13
300 2.84E-11 1.31E-12
表3 系統在時間步為0.0005 皮秒，模擬圈數為 20000 之熱導與誤差。
47
我們畫出一個苯環系統，在0.001 皮秒與 0.0005 皮秒之兩個不同時間步 下所計算出來的熱導-系統溫度圖，如下：
圖16 一個苯環系統，在兩個 0.001 和 0.0005 之時間步下所計算出來的熱 導-系統溫度圖。
48
二、 兩個苯環在不同時間步與模擬圈數下之熱導與誤差：
Temperature(K) Thermal Conductance(W/K) Error 30 9.00E-12 4.18E-13
Temperature(K) Thermal Conductance(W/K) Error 30 8.95E-12 3.29E-13
49
我們畫出兩個苯環系統，在0.001 皮秒與 0.0005 皮秒之兩個不同時間步 下所計算出來的熱導-系統溫度圖，如下：
圖17 兩個苯環系統，在兩個 0.001 和 0.0005 之時間步下所計算出來的熱 導-系統溫度圖。
50
最後，我們將一個苯環系統與兩個苯環系統，在時間步為0.0005 皮秒，
模擬圈數為20000 的參數設定下，比較計算出來的熱導-系統溫度圖，如下：
圖18 在 0.0005 之時間步下對一個苯環系統與兩個苯環系統模擬計算出 來的熱導-系統溫度圖。
51
第四章 結論
本篇論文主要以分子動力學方式研究分子接面的熱傳導性質，並且探討 在不同溫度下的相關性質。我們將苯環與金電極連接，透過結合分子動力 學（Lammps）與密度泛函理論軟體(VASP)來計算出分子接面的熱傳導性質，
當分子放置在兩個電極中，在不同的溫度下將會有不同的熱導，並找出溫 度與其對應之熱導的關係，然後更進一步計算一個苯環與兩個苯環系統之 間的差異。
由於密度泛函理論軟體(VASP)是從第一原理出發，其計算量非常的龐 大，所以我們採取先執行分子動力學（Lammps）的模擬。首先，我們對一 個苯環的系統，做了相同的總時間(1 10 秒)下，模擬不同時間步與模擬 圈數，其結果差異不大，代表當系統收斂達穩定以後，熱導將不受到時間 步與模擬圈數的影響，在兩個苯環的系統也有相同的結果。
再來我們比較一個苯環系統與兩個苯環系統的差異，發現跟巨觀下熱 導與長度成正比不相同，在這兩個系統裡，兩個苯環的熱導反而優於一個 苯環的熱導，這與前面圖8 ，Bulk 在 Z 方向大小與熱導率的關係，有一樣 的結果，原因為 Z 方向長度太小時，系統將不會考慮到長波長的影響，所 以兩個苯環的系統之熱導反而優於一個苯環的系統。
52
最後，我們透過分子動力學模擬的計算，讓我們更了解在原子尺度下 的分子接面的熱傳導性質，希望藉由這些分析對未來在製作奈米尺度的元 件上有所幫助。
53
第五章 未來展望
一、 目前估算溫度的方法仍然是古典的統計力學中，由理想氣體方程式 推導計算而來，這只適用於高溫的系統，當溫度降低到德拜溫度以下，
就必須考慮到量子效應，因此必須引進phonon density of states(PDOS)。
二、 可將系統推廣至多個苯環，如此一來當系統長度足夠考慮長波長時，
就有機會消除兩個苯環之熱導優於一個苯環的情況。
三、 為了從第一原理的觀點出發，必須將密度泛函理論軟體(VASP)引進 來計算原子間的力，然後將上述之結果繼續模擬計算。
四、 將密度泛函理論軟體(VASP)引進計算當中時，將可以使用贋勢來模 擬計算，便可不受限於分子動力學的force field。
五、 目前的計算，是使用CPU 來計算，但是由於第一原理的密度泛函理 論軟體(VASP)，其計算量非常的龐大，如果改用 GPU 來計算，便可以 有效加快運算的速度。
六、 可以將整個系統推廣到不同的分子接面，希望藉由這些分析對熱傳 導性質有更多的了解，更希望對未來在製作奈米尺度的元件上有所幫 助。
54
參考文獻
1. Zhixin Guo, Dier Zhang, and Xin-Gao Gong, Thermal conductivity of grapheme nanoribbons. Appl. Phys. Lett. 95, 163103(2009)
2. Nuo Yang, Gang Zhang, and Baowen Li, Ultralow Thermal Conductivity of Isotope-Doped Silicon Nanowires. Nano Lett., 8 (1), pp 276–280(2008) 3. G. Mie, Ann. Physik 11, 657 (1903)
4. J. E. Lennard-Jones,Proc. Roy. Soc. ,106A, 441(1924) 5. J. E. Lennard-Jones,Proc. Roy. Soc. ,106A, 463(1924)
6. P. M. Morse, Diatomic molecules according to the wave mechanics. II.
Vibrational levels. Phys. Rev. 1929, 34, 57-64 7. M. Born and J. E. Mayer, Z. Physik 75, 1 (1932).
8. M. S. Daw and M. I. Baskes, "Embedded-atom method: Derivation and application to impurities and other defects in metals" ,Phys Rev B29, 6443 (1984)
9. Finnis, M. W. , Sinclair J. E. ,"A simple empirical N-body potential for transition metals". Phil. Mag. A 50 (1): 45. (1984)
10. J. H. Irving, and J. G. Kirkwood,“The Statistical Mechanics of Transport Process. IV. The Equation of Hydrodynamics”, J. Chem. Phys., 18, pp.817-820(1950)
11. Alder, B. J. and Wainwright, T. E., J. Chem. Phys. 27, 1208 (1957)
12. Gibson, J. B., Goland, A. N., Milgram, M., and Vineyard, G. H., Phys. Rev.
120, 1229 (1960)
13. A. Rahman, “Correlations in the Motion of Atoms in Liquid Argon”,Phys.
Rev. 136A, 405(1964)
14. Girifalco, L. A. and Weizer, V. G., “Application of the Morse Potential Function to Cubic Metails”,Phys. Rev., 113, No.3 ,p687-690(1959) 15. L. Verlet, “Computer ‘Experiments’on Classical Fluids. I.
Thermodynamical Propertier of Lennard-Jones Fluid,”Mol. Phys. Rev., 159, 98(1967)
16. Reif, F., Fundamentals of Statistical and Thermal Physics, McGraw Hill, pp.48-49(1985)
17. Haile, J. M., Molecular Dynamic Simulation: Elementary Methods, John Wiely & Sons, Inc., USA(1992)
18. Hohenberg, P. and W. Kohn, Inhomogeneous electron gas. Phys. Rev, 1964.136(3B):p. B864-B871.
55
19. Kohn, W. and L. Sham, Self-consistent equations including exchange and correlation effects. Phys. Rev, 1965. 140(4A):p. A1133-A1138.
20. 江進福，波函數與密度泛函，物理雙月刊，二十三卷五期，2001 年 10 月。
21. J.D. Pack H.J. Monkhorst, Phys. Rev. B16, 1748 (1977) 22. D.R.Hamann,M.Schluter and
C.Chiang,Phys.Rev.Lett.43,pp.1494-1497(1979)
23. Xavier Gonze, Peter Kackell, and Matthias Scheffler, Phys. Rev. B41, p.12264 (1990)
24. David Vanderbilt, Phys. Rev. B41 p.7892 (1990)
25. Florian Müller-Plathe,J. Chem. Phys. ,Vol. 106, No.14, 8 April 1997.
26. I.T. Todorov & W. Smith,THE DL_POLY_4 USER MANUAL, Version 4.01.0, October 2010
27. Steve Plimpton, Aidan Thompson, and Paul Crozier ,LAMMPS Users Manual