國 立 交 通 大 學
機械工程學系
碩 士 論 文
比較兩種位能函數在矽沉積模擬系統之異同
Comparisons between two potential functions for MD simulation
on silicon depositions
研 究 生: 賴益璋
指導教授: 陳宗麟 博士
比較兩種位能函數在矽沉積模擬系統之異同
Comparisons between two potential functions for MD
simulation on silicon depositions
研究生 : 賴益璋 Student: Yi-Jhang Lai
指導教授 : 陳宗麟 博士 Advisor:Dr. Tsung-Lin Chen
國 立 交 通 大 學
機 械 工 程 學 系
碩 士 論 文
A Thesis
Submitted to Department of Mechanical Engineering College of Engineering
National Chiao Tung University in Partial Fulfillment of the Requirements
for the Degree of Master of Science
in
Mechanical Engineering September 2008
Hsinchu, Taiwan, Republic of China
比較兩種位能函數在矽沉積模擬系統之異同
學生:賴益璋 指導教授:陳宗麟博士國立交通大學 機械工程學系 碩士班
摘要
位能函數決定原子的運動,也決定系統的特性。本論文利用分子動力學(Molecular dynamics)模擬以分子束法(Beam epitaxy)進行的矽沉積系統。矽原子間的位能關係 分別採用 S.W.及 B.H.兩種位能函數,並藉由比較兩種位能函數和模擬所得的晶體結構 的差異,從而得知兩位能函數的優劣。 常見於分子動力學模擬中的空間分割(Cell subdivision)技術,目前僅適用於 two-body potential 的位能計算。本論文以原有的空間分割理論為基礎,提出一套適用 於 three-body potential 計算的空間分割方式。希望藉由本方法有系統的空間選取,減 少重覆計算和除去作用力有效範圍外的計算,進而大幅降低模擬程式的計算量。本論 文詳細說明此一新式的空間選取法則,並將其應用於本論文的程式模擬上。Comparisons between two potential functions for MD simulation
on silicon deposition
Student:Yi-Jhang Lai adviser:Dr. Tsung-Lin Chen
Department of mechanical Engineering National Chiao-Tung University
Abstract
Potential functions, which describe the energy relations between atoms in a material, can determine the movements of atoms and thus properties of a material. In the thesis, we use the molecular dynamics (MD) to simulate a beam expitaxy silicon deposition system. Two well-known potential functions for silicon atoms, S.W. and B.H., were used in these simulations respectively for comparisons.
“Cell subdivision” techniques are very popular in MD simulations because they can greatly reduce the computation load. However, they are applicable to the calculations of “two-body potential,” but not the “three-body potential.” In this thesis, we develop a systematic way of selecting spaces in cell subdivisions so that the method can be extended to the case of “three-body potential.” This method is discussed in details and used in the MD simulations in this thesis.
誌謝
二年來的耕耘,終於使本篇論文結纍了,能走完這段辛苦的旅程,我的指導教授 陳宗麟老師惠我良多。在研究陷入困境時,如良師般引導我如何解決;在苦悶的研究 過程中,如益友般開導我面對;令我印象最深的是,老師那開明的教學態度,以鼓勵 代替強迫,以責任代替規定,以引導代替教學,使我在研究過程中,享受到絕大部份 的自由,讓我在沒有研究進度的壓力下,培養對自己負責的態度,增加解決問題的經 驗。對於老師的諄諄教誨,銘記於心,對於老師無私的教導,在此獻上誠摯的敬意和 謝意 一起抱怨,讓我在研究過程中不覺得孤 單, 盾,適時的關心和安慰,溫暖了我寂寞的心,讓我能順利完成這段辛苦的旅程。 。 研究生活是苦閟的,在此感謝實驗室學長、同學和學弟,以及大學同學、室友和 好友的陪伴,一起討論學業,一起出遊玩樂, 和寂寞的研究生活交織成難忘的回憶。 最後要感謝我的父親、母親、弟弟和妹妹,有他們在我身後無條件支持,做我的 後目錄
中文摘要 英文摘要 誌謝 目錄 圖目錄 第一章 諸論 動機與目的 文獻回顧 論文大綱 第二章 分子動力學模擬 牛頓力學 位能函數 交互作用 無因次 初始條件 邊界條件 積分法 空間分割 第三章 矽模擬計算理論 設定參數 初始位置和速度 位能函數 ...i ...ii ...iii ...iv ...vi ... 1 1.1 ... 1 1.2 ... 2 1.3 ... 3 ... 4 2.1 (Newton’s law) ... 4 2.2 (Potential function) ... 5 2.3 (Interactions)... 6 2.4 (Dimensionless) ... 7 2.5 (Initial conditions) ... 8 2.6 (Boundary Conditions) ... 9 2.7 (Integration methods) ... 10 2.8 (Cell subdivision) ... 11 ... 13 3.1 (Set parameters)...143.2 (Initial coordinates and velocities) ...14
3.4 空間分割(Cell subdivision) ... 17
3.5 力的計算(Compute forces) ... 22
3.6 積分法-蛙跳法(Leapfrog methods) ... 24
3.7 邊界條件-週期性邊界(Periodic boundary condition) ... 25
3.8 溫度控制-定溫(Constant temperature) ... 25 第四 第五 考文獻 ...38 章 矽沉積模擬 ...26 4.1 使用相同單位 ... 26 4.2 位能函數分析 ... 27 4.3 參數設定 ... 29 4.4 模擬流程 ... 30 4.5 模擬結果 ... 30 章 結論與未來計劃 ...36 5.1 結論 ... 36 5.2 未來計劃 ... 37 參
圖目錄
圖 4-10 B.H.結晶上視圖(z=1 ~ z=5) ...35 圖 2-1 原子向量示意圖...5 圖 2-2 Lennard-Jones 位能函數圖 ...6 圖 2-3 週期性邊界圖... 9 圖 2-4 空間分割圖...12 圖 3-1 MD 模擬流程圖 ... 13 圖 3-2 鑽石晶格結構圖...15 圖 3-3 三原子結構圖...17 圖 3-4 two-body 空間分割圖 ...19 圖 3-5 three-body 三維空間分割圖 ...20 圖 3-6 three-body 三維空間分割圖(灰色區域)...21 圖 3-7 四排列組合...22 圖 3-8 three-body 力矩向量圖 ...24 圖 4-1 S.W.和 B.H.的 two-body potential ...28 圖 4-2 S.W.和 B.H.(已改變參數)的 two-body potential...28 圖 4-3 S.W.和 B.H.的 three-body potential ...29 圖 4-4 矽沉積模擬流程圖...30 圖 4-5 S.W.位能函數結晶圖 ...31 圖 4-6 動能變化圖(S.W.)...31 圖 4-7 完全結晶矽(鑽石結構) ...32 圖 4-8 鑽石和纖維鋅礦晶格結構上視圖...33 圖 4-9 B.H.(已改變參數)位能函數結晶圖...34第一章
諸論
1.1 動機與目的
在目前的製程中,CMOS(Complementary metal oxide semiconductor)和 BJT(Bipolar junction transistor)等半導體元件,大都是以單結晶矽(Mono-crystalline silicon)為原料加 以製成,而得到的單結晶矽經過一連串的加工過程後,將可做為半導體元件並使用之。 如果能把加工步驟的一部分和單結晶矽的形成製程合併,那麼將可節省不少的時間。 單結晶矽主要是利用 epitaxial growth[1]所製成,分為氣態沉積、液態沉積和分子 束長晶(Molecular beam epitaxy)三種方法,其中分子束長晶的速度雖然最慢(每小時不 超過 1000nm),但是可以藉由控制分子束的速度,得到不同的排列方式。 尋找合適的參數來達到想要的晶體結構,在實驗上是所費不貲的,若利用現今電 腦的快速運算,和科學家推導的原子位能函數,建構出真實情況的模型,藉由模擬模 型中系統的動態,來尋找合適的參數,將能有效地減少實驗上錯誤的嘗試,進而降低 高昂的研究費用。 位能函數決定了原子的特性,函數愈近似真實原子,模擬出來的參數就愈接近真 實,本論文以半導體工業最重要的元素-矽,利用分子動力學模擬(Molecular dynamics simulation)和兩種不同的矽位能函數,模擬原子沉積並結晶的過程,來計算出分原子 在微觀系統中的軌跡移動及當時的速度,從這兩個資訊推算出當時的晶體結構和平均 溫度。進而討論兩種矽位能函數的特性及優缺點。
除此之外還有另一個重點:three-body 空間分割法。在計算 three-body potential 的 過程中,以 two-body 空間分割法為理論基礎,配合有系統地且不重覆地選取,推得 three-body 空間分割法,可除去多餘的距離長度判斷(在作用力範圍之內或之外),大幅 簡化 three-body potential 的計算。
1.2 文獻回顧
本論文所採用的兩種位能函數為 S.W.(Stillinger and Weber)位能函數[2]和 B.H.(Biswas and Hamann)位能函數[3]。
S.W.位能函數是以統計力學[4][5][6]和 L-J(Lennard-Jones)位能函數為基礎,推導 出了由 two-body 和 three-body 所組成的位能函數,用來描述矽在液態和固態的交互作 用。雖然 S.W.位能函數有必須在真空的先決條件限制[2],以及無法清楚分辨鑽石晶格 結構和纖維鋅礦(wurtzite)晶格結構的缺點[7],但是因為 S.W.位能函數在模擬中的反應 很快(以快速收斂為目的設計而成)[2],所以還是常被用來模擬矽的交互作用。例如以 (100)和(111)方向排列的結晶矽和熔融速度的關係[8],在(100)或(111)表面模擬 epitaxial growth[7][9],模擬非結晶矽的形成[10]。
B.H.位能函數是以 LDA(Local density approximation)為計算方法,去推得符合真實 矽的參數,可分為 two-body 和 three-body potential,適用液態和固態矽的交互作用。 B.H.位能函數雖然比 S.W.來的精確,但是沒有作用力範圍的設定,使其在計算作用力 時,會需要大量的運算時間,因為如此,只有被少許論文用來模擬矽的交互作用。例 如沿(111)方向模擬 epitaxial growth[11]。 分子動力學模擬中的計算簡化方法-空間分割法,目前應用在 two-body potential 所推導的作用力計算上[13][14]。空間分割法參考位能函數的作用力有效範圍,把系統 切割成相同的區塊大小,輪流選取每個區塊為計算起點,以作用力與反作用力為原理 基礎,對稱簡化切割好的區塊,再對簡化後的區塊做計算。比起直接選取原子,空間 分割法除去了重覆和作用力有效範圍外的計算,減少無謂的模擬時間。 我們模擬的矽沉積系統,是由分子束長晶(Molecular-beam epitaxy)來射入原子。基 底表面的原子因為只受到基底下原子的作用力,所以會產生表面重建的現象[12],表 面原子不會在結晶位置上,而是在最低能量的位置上,而射入的原子能量必須大於熔 點,撞擊表面原子使其有足夠的能量遷移(可視為熔融狀態),遷移過程中因能量流失 而穩定,再受撞擊而遷移,直到有足夠原子在其之上並構成新表面,才可能穩定形成
結晶層。 液態矽主要是以簡單流體來運動,主要和 two-body 函數有關,而固態矽則是以傾 向固定結構來運動,和 three-body 函數有關,S.W.和 B.H 位能函數都含有 two-body 和 three-body 的形式,因此能成功的描述矽的熔融狀態和固態。 S.W.和 B.H.位能函數都有和 epitaxial growth 相關的論文[7][9][11],但是卻沒有完 整的結晶圖,因此本論文的模擬除了比較兩種位能函數的特性外,還可驗證其是否能 正確的描述矽的結晶。 空間分割法是利用作用力和反作用力的理論,減少作用力的計算量。本論文將利 用上述原理,對稱簡化需要計算的區塊後,以計算起點與區塊的距離和平面的異同為 分類依據,逐一討論區塊的排列並避免重覆選取,推得空間分割法在 three-body potential 的應用。
1.3 論文大綱
本篇論文有兩大重點,其中之一為模擬矽沉積系統,用分子束法進行 epitaxial growth,由分子動力學模擬和 C 語言寫成,探討不同位能函數所形成的晶體結構有何 不同;另一個為 three-body potential 空間分割法。第二章介紹模擬需要使用到的力學, 和本篇論文最主要的模擬方法-分子動力學模擬;第三章講解矽沉積和分子動力學模 擬的結合與模擬所需的理論和進行方法,以及改進 three-body potential 的計算方式- three-body potential 空間分割法;第四章為矽沉積模擬與兩位能函數和晶體結構的討 論;第五章總結結論和未來發展;最後為所參考的文獻。第二章
分子動力學模擬(molecular dynamics simulation)
為了計算出分子最後的位置,也就是最穩定、能量最低的分子結構,我們選擇了 分子動力學模擬,來做為數值模擬的方法。分子動力學模擬是用來計算固體、液體和 氣體單原子或單分子的運動,並用以時間為變數的函數-位置、速度和方向,來描述 所模擬原子和分子的運動。本章將說明分子動力學模擬的理論和方法。2.1 節介紹基本 力學;2.2 節和 2.3 節例舉描述惰性氣體的位能函數,並推得原子間的作用力公式;2.4 節說明無因次在分子動力學模擬的應用;2.5 節說明初始條件的設定;2.6 節介紹週期 性邊界條件;2.7 節說明藉由積分法得到原子的速度和位置;2.8 節講解空間分割法的 原理與用法。
2.1 牛頓力學(Newton's law)
牛頓第一運動定律: 物體不受外力作用,或所受外力合力為零時,原來為靜止的物體還是靜止, 原來為運動的物體則呈等速度直線運動,又稱慣性定律。 牛頓第二運動定律: 物體受力後的加速度,和物體所受的合力成正比,和物體質量成反比。 ri:原子位置向量 2 2 2 i i i d r F mr m d t = = (2.1)Z Y
r
i atom i X 圖 2-1 原子向量示意圖 牛頓第三運動定律: 當兩物體互相對彼此作用力時,此兩作用力大小相等,方向相反,並分別作 用在不同物體上,又稱為作用與反作用定律。 12 21 F = −F (2.2)2.2 位能函數(Potential function)
選擇一個簡單的位能函數來說明,Lennard-Jones 位能函數(Lennard-Jones potential function)是由數學家 John Lennard-Jones 於 1924 年所提出的,由於其方程式簡單而被 廣泛使用,特別是用來描述惰性氣體原子間相互作用尤為精確。 (2.3)式為 Lennard-Jones 位能函式。 12 6 ( )ij 4 [( ) ( ) ] ij ij u r r r σ σ ε = − , rij≤rc (2.3) ( )ij u r :兩原子的位能值 c r :兩原子間無作用力的最短距離 ε:作用力強度的單位 σ :距離長度的單位 ij r :兩原子間的距離
(2.3)式中的 12 次方項,代表的是兩原子間的斥力,是由尚未鍵結且重疊的電子雲 所產生的,而(2.3)式的 6 次方項,代表的是吸力,是由凡得瓦力(Van der Waals force) 所產生的。把(2.3)式繪圖後可得圖 2-2,在原子間距離為 2.5σ 時,能量變化不大且幾 乎為 0,故可設定rc =2.5σ 。 0 0.5 1 1.5 2 2.5 -1 0 1 2 3 4 5 兩原子間的距離(σ) 兩 原 子 的 位 能 值 (ε) 圖 2-2 Lennard-Jones 位能函數圖
2.3 交互作用(Interactions)
為了要模擬原子的運動,需取得原子的位置和速度,這二個參數可由加速度求得, 從(2.1)式了解到,加速度可由作用力和質量的關係求得,作用力是由所選擇的位能函 數所推導出來的。因此,位能函數的選取,影響到整個系統所有的狀態(溫度、壓力等)。 (2.4)式為作用力和位能的關係。 ( ) f = −∇u r (2.4) 把(2.3)式代入(2.4)式,可推得原子間的作用力。14 8 2 48 1 ( )[( ) ( ) ] 2 ij ij ij ij f r r r ε σ σ σ = − , rij ≤rc (2.5)
2.4 無因次(Dimensionless)
原子的大小只有奈米的等級,相關的參數大都不是簡單的整數或普通的分數,所 以原子的運算就會變的很繁雜,有鑑於此,我們利用單位的無因次化來簡化它。 使用無因次化有以下優點: 1. 有些繁雜的參數由於無因次化,而被吸收進單位中,簡化了運動方程式。例 如(2.8)式中的質量 m ,併入於時間的單位後,(2.1)式就簡化為作用力等於加 速度。 2. 使參數間的數值不會差異過大,降低數值運算錯誤的風險 對長度、質量、能量、時間和溫度無因次化: 長度: r→rσ (2.6) 質量: e→eε (2.7) 時間: 2 / t→t mσ ε (2.8) 溫度:K →Kε/kB (2.9) 把(2.7)、(2.8)和(2.9)式代入由(2.4)式組成且含 N 顆原子的 d 維系統,可推得簡化 後的式子。 每顆原子的作用力方程式: 14 8 ( ) 1 48 ( ) 2 ij ij ij ij j i f r− r− ≠ =∑
− r (2.10) 每顆原子的平均動能: 2 1 1 2 N k i i E v N = =∑
(2.11) 每顆原子的平均位能: 12 6 1 4 ( u ij i j N ) ij E r r N − − ≤ < ≤ =∑
− (2.12)每顆原子的平均溫度: 2 1 1 N i i T dN = =
∑
v (2.13)2.5 初始條件(Initial Conditions)
為了計算原子的運動,初始條件的設定是必要的,一般來說,初始條件的設定可 以不用太過準確,因為系統會以能量最低的形式,去進行系統內原子的運動,若設定 一個接近穩定狀態的初始條件,則可有效的降低從初始條件到穩定狀態所需要的電腦 運算時間。 在初始位置的設定上,我們常常有規則的排列,將原子平均的配置在所要模擬的 系統中。除了平均散佈外,原子的相對位置愈接近穩定狀態,平衡所需要的運算時間, 也將愈少。因此,我們所要模擬的原子,大多使用固體狀態的晶格結構排列,既符合 平均,也可減少運算,不過使用晶格結構的初始位置設定,比較適合固體和液體的原 子模擬,至於氣體的原子模擬初始位置設定,則符合平均的配置就好。 在初始速度和溫度的設定上,如果對原子沒有特別的要求(例如特定方向的粒子 束),速度方向都是視系統維度去給任意方向的,而速度的大小則是由初始系統溫度去 設定,根據(2.13)式,可以反推得到式子右邊的速度平方的總和,接下來有兩種方法去 決定各個原子速度的大小,一個是把速度平方的總和,平均後開根號,讓所有的原子 速度大小都一樣 v= dT (2.14) 另一個是利用馬克斯威爾-波茲曼速度分佈(Maxwell-Boltzmann distribution of speed) 2 3/ 2 2 ( ) 4 ( ) exp[ ] 2 2 m f v v kT kT π π − = mv (2.15) 由於(2.15)式是以氣體為假設前提下所推導而成的,所以適用於氣體[11]。而固體和液 體的速度大小,則是用前一個方法去設定。2.6 邊界條件(Boundary Conditions)
使用分子動力學模擬時,除非目的是為了模擬邊界附近的特性,否則為了減少計 算量,大都採用週期性邊界條件(Periodic boundary conditions)做為計算的方法。
C D B A 圖 2-3 週期性邊界圖 以二維系統為例,圖 2-3 正中央黑框部分為我們模擬的區塊,週期性邊界簡單地 說,就是複製模擬的區塊,並把它放在模擬區塊的四周,當要用邊界條件來計算時, 便使用剛剛複製的區塊做計算。也可以看做是,在一個廣大的系統中,取一個很小的 區塊來模擬,而其他區塊都用此小區塊來無限複製,而小區塊模擬出來的特性,可代 表整個區塊的特性。 注意事項: 1. 因為模擬區塊的四周皆為模擬區塊的複製,見圖 2-3,當原子 A 離開模擬的 區塊並進入右邊的複製區塊,同樣地在左邊的複製區塊也有原子 B 同時進入 模擬的區塊中,因此,在模擬原子運動時,原子越過邊界離開模擬的區塊後,
要立即地從另一邊的邊界,再進入模擬的區塊裡。 2. 計算原子和原子之間的作用力,雖然只需要算模擬區塊內的原子即可,但如 果某原子有效作用力的範圍,超過了模擬的邊界,則必須要考慮到相鄰複製 區塊中的原子,對模擬區塊內原子的影響。也就是要把邊界外有效力範圍內 的原子,納入模擬區塊內原子的計算。 3. 為了要使用週期性邊界,模擬區塊形狀的選擇是很重要的,以二維系統為例, 矩形是最簡單的區塊形狀,但區塊形狀的選擇不限於矩形,只要複製後能填 滿模擬區塊四周的空間,該區塊形狀都可以被拿來使用。例如菱形、六角形 等。 4. 如果要使用週期性邊界條件,則模擬區塊的大小的選取有所限制,過大會增 加過多的計算量,失去使用週期性邊界的意義;過小會有重覆計算的現象產 生。見圖 2-3,假設模擬區塊橫軸長為1.8r (c r 為兩原子間無作用力的最短距c 離),AC=0.85rc,推得AD=0.95rc,原子 C 和原子 D 都在原子 A 的作用力 範圍內,同一原子重覆計算是不合理的,所以區塊的最短距離(相對邊界的距 離),不得小於原子間作用力有效範圍的二倍。
2.7 積分法(Integration methods)
利用由位能函數推導得到的運動方程式,可以從原子的位置得到原子之間的作用 力,推得原子的加速度,積分後可得速度和位置,積分運算可由蛙跳法(Leapfrog methods)和預估修正法(Predictor-corrector methods)等數值運算方法來近似。在此簡單 地介紹本文所採用的積分法-蛙跳法。 蛙跳法為最簡單的數值積分方法,雖然相較於其他複雜的數值積分,蛙跳法求出 來的數值準確率比較低,但是在大量數值運算下,能確實縮短計算時間。 蛙跳法的計算公式如下: ( ) ( ) ( ) 2 2 i i h h v t+ =v t + ×a ti (2.16)( ) ( ) ( 2 i i i h r t+h =r t +hv t+ ) (2.17) h:取樣時間 由於取樣時間很小和加速度的變化不大,所以為了計算方便, ( 2 i h v t+ ) ) 可視為 。若是想更為精確,可以把以上兩式在計算作用力之前先執行,計算作用力之 後再執行(2.18)式,將可得到符合該時間的速度值,相對的,也要付出較多的時間去計 算之。 ( i v t+h ( ) ( ) ( 2 2 i i i h h v t+h =v t+ + ×a t+ )h (2.18)
2.8 空間分割(Cell subdivision)
因為計算量是隨著原子的數量呈現次方成長,所以如果模擬的原子數量很大,那 麼所需的時間就會非常長,為了有效的減少計算所需的時間,我們對模擬的原子進行 空間的分割,只需計算在作用力範圍內的空間即可。在參考作用力的有效範圍後,決 定分割的長度(通常是等於或略大於有效範圍),對整個系統做相等大小的分割,然後 把所有的原子依照它們各自的位置,分配到相對應的空間,計算力的時候,以空間為 單位,某空間只需要和相鄰的空間以及自身空間內,進行作用力的運算即可。 注意事項: 1. 使用空間分割時,各個空間大小要相等 2. 模擬系統不可太小,否則空間分割的效用不大 3. 當相鄰的空間超過邊界時,利用週期性邊界的觀念,取另一邊的空間做 運算。圖 2-4 為圖 2-3 中央模擬區塊的空間分割圖,空間 A 除了和空間 BCDEF 計算外,還需和空間 GHI 計算。I
H
G
E F
D
B C
A
圖 2-4 空間分割圖第三章
矽沉積之模擬計算理論
在了解分子動力學模擬的原理後,為了應用在矽原子的沉積上,本章將把理論和 矽原子的特性結合,選取適合的方法去模擬它。在本章的 3.4 節,有著本論文另一個 重點:three-body 空間分割法的使用。此方法以 two-body 空間分割法為理論基礎推導 而成,採用系統地選取空間組合,來減少 three-body potential 的計算量。 模擬的內容包括,相關的參數標準化、初始值設定、方程式設定、邊界條件和溫 度設定,流程請見圖 3-1。Set parameters
Initial coordinates
Initial velocities
for each atom
Cell subdivision
Compute forces
Leapfrog method
Periodic boundary condition
Constant temperature
Add a atom
Count < final count
Count ≧ final count
End
3.1 設定參數(Set parameters)
自然界的參數大部分都是很繁雜的大數或小數,極少有漂亮的整數或是位數少的 數,矽的各項參數也不例外,因此,利用無因次化來簡化常數,將會對展現數據變化 和減少數值運算的錯誤大有幫助。由 2.4 節可以知道無因次化的好處。我們所做的模 擬當然也使用了無因次化,在單位長度、單位能量、單位時間和單位溫度都有簡化繁 雜計算的好處,本文直接利用 S.W.位能函數的無因次化單位,其中單位溫度為了計算 方便,取整數運算,和正確的溫度有些微誤差(大約 50K)。 單位長度:σ =0.20951nm=2.0951Å (3.1) 單位能量: -19 (3.2) 3.4723 10 J ε = × 單位時間:t*=σ(
m/ε)
1/ 2 =7.6634 10× −14s (3.3) 4 K 單位溫度: * (3.4) / B 2.52 10 T =ε k = ×3.2 初始位置和速度(Initial coordinates and velocities)
由於要模擬的過程是單結晶矽的 epitaxial growth,而做為原子依附的基底,就必 須是單結晶矽的排列結構。矽的晶格結構是鑽石晶格(見圖 3-2),矽原子的四個鍵結分 別連接四個矽原子,這四個矽原子組成正四面體,正四面體的重心(也是外心和內心) 位置有個這四個矽原子所共同鍵結的矽原子,依照此正四面體的模型,可用數學的方 法計算出單結晶矽內各原子的相對位置,依此相對位置給予原子座標模擬。要被射入 基底的矽原子位置則是設定在固定範圍高度的隨機位置。 為了避免射入的原子直接被降溫(詳見 3.8 節),所以要在溫控層上加上適量的緩衝 層,避免入射原子和溫控層原子之間能量轉換過快,而使得原子彈飛基底,或是刺穿、 鑲入基底,增加模擬的誤差。 原子速度由溫度所推導而得,詳細內容參考 2.5 節,基底的矽原子速度用(2.14)式 推得,而要被射入基底的矽原子速度則是要利用(2.15)式推得,方向垂直基底往下。
圖 3-2 鑽石晶格結構圖
3.3 位能函數(Potential function)
本論文模擬所用的位能函數有兩個:S.W.和 B.H 位能函數,都是用來描述矽原子 間的交互作用。
1. S.W. 位能函數:由 two-body potential 和 three-body potential 所組成,當距離大於作 用力有效距離時,則位能為 0,藉此可以忽略掉間距過長的原子作用力,減少計 算所需的時間。two-body potential 是兩原子間距離的函數,能使原子們趨向保持在 適當的距離;three-body potential 是兩個兩原子間距離和其夾角的函數,能使原子 們趨向保持在固定的角度,形成某個結構,對於本論文的目的-矽結晶,是重點 所在。 Two-body potential: 1 2( ) ( 1) exp[( ) ] p ij ij ij f r = A Br− − r −a − , rij <a (3.5) Three-body potential,見圖 3-3:
3( , , )i j k ( ,ij ik, jik) ( ,ji jk, ijk) ( ,ki kj, ikj)
f r r r =h r r θ +h r r θ +h r r θ , r<a (3.6) 1 1 1 ( , , ) exp[ ( ) ( ) ] (cos ) 3 ij ik jik ij ik jik h r r θ =λ γ r −a − +γ r −a − × θ + 2 (3.7)
:兩原子間的距離 ij r :兩原子間無作用力的最短距離 a , , , , A B p λ γ :S.W.位能函數的參數,詳見 4.1 節
2. B.H. 位能函數:由 two-body potential 和 three-body potential 所組成,和 S.W.位能函 數不同的是,沒有作用力有效距離的設定,改用 cut-off 函數取代,如此作法可得 到更為精確的模擬結果,但也需要更大量的時間去計算,因此可視所需精確度的 多少,設定使計算有效率的作用力有效距離。B.H.位能函數的變數和設計,基本 上和 S.W.位能函數相同,不同的是 three-body potential 增加了 1 3 (cos ) 3 jik θ + 項,在 非最穩定角( 11 cos 3 jik θ = − )的情況下,對角度提供擾動,使結構不會停留在亞穩態。 Two-body potential: 2 2 2( )ij ( 1exp( 1 ij ) 2exp( 2 ij )) ( )c ij f r = A −λr +A −λ r f r (3.8) Three-body potential,見圖 3-3:
3( , , )i j k ( ,ij ik, jik) ( ,ji jk, ijk) ( ,ki kj, ikj)
f r r r =h r r θ +h r r θ +h r r θ (3.9) 2 3 1 1 1 2 2 2 1 1 ( , , ) [ ( ) ( )(cos ) ( ) ( )(cos ) ] 3 3
ij ik jik ij ik jik ij ik jik
h r r θ = Bψ r ψ r θ + +Bψ r ψ r θ + ( ) ( ) c ij c ik f r f r × (3.10) Cut-off 函數: 1 ( ) {1 exp[( ) / ]} c ij ij c f r = + r −r μ − (3.11) 收斂函數: 2 ( ) exp( ) i rij i ijr ψ = −α (3.12) :兩原子間的距離 ij r 1, 2, 1, 2, 1, 2, 1, 2, ,c A A B B λ λ α α r μ :B.H.位能函數的參數,詳見 4.1 節
θ
JIKr
IKr
JKr
IJθ
IKJθ
IJKI
K
J
圖 3-3 三原子結構圖3.4 空間分割(Cell subdivision)
因為我們要模擬的原子數量不少,所以要利用空間分割來減少計算量,詳細方法 請見 2.8 節,參考作用力的效範圍並分割原子到各個空間後,開始計算之間的作用力, 計算力的方法可以分成直接計算和對稱計算。 1. 直接計算:針對其中一顆原子,把其所在的空間以及相鄰的空間內的原子們,計算 它們對此原子的作用力並加總,得到此原子的作用力。依序地計算每一顆原子所受 的力,雖然計算費時,但程式的撰寫簡單。 2. 對稱計算:原理是利用(2.2)式,作用力和反作用力的關係,還有 2.6 節、2.8 節和 3.7 節講解的週期性邊界的影響。有別於直接計算,只需一半多的計算量,不過空 間的取法卻相當複雜且容易出錯。 3.4.1 對稱計算 以二維空間系統為例,見圖 3-4a,假設先以空間 1 為計算起點,對作用力可達的 範圍做運算,和直接計算不同的是,空間 1 只需和空間 1、2、3、4、5 做運算,但不代表空間 1 沒有和空間 6、7、8、9 做運算,而是當以空間 6 為計算起點時,已和空間 1 做運算了,同理,空間 7、8、9 亦然。 圖 3-4a two-body 二維空間分割圖 3.4.2 two-body potential 的空間分割 本論文所模擬的系統為三維空間系統,見圖 3-4b,因為我們是以作用力的有效範 圍作為切割的長度,所以如果以黑色方塊為起點作直接運算,要計算的空間就如圖 3-4 所示,共 3 3 3 個方塊空間。為了減少計算所需的時間,利用在上小節說明的內 容,在三維空間時,利用透明和灰色方塊對稱的性質,只需對黑色和透明方塊進行運 算即可,雖然增加了些許程式寫法的難度,不過可大幅地縮減計算所需的時間。 27 × × =
圖 3-4b two-body 三維空間分割圖 3.4.3 three-body potential 的空間分割 上小節所提及的 two-body potential,有效作用力範圍所涵蓋的空間,可包括 個方塊所組成的正立方體。本小節所要討論的 three-body potential,則要 多算一個有效作用力的距離,空間包括 5 5 5 3 3 3× × =27 125 × × = 個方塊所組成的正立方體,若是 對正立方體做直接計算,則共要取125 125 15625× = 次方塊位置,而有效的僅有數百種 取法,十分浪費運算時間,因此,我們利用對稱計算去簡化它。 由於三個原子可能在同一個方塊中,或是二個在同方塊,或是三個都不在同方塊, 以下針對同異方塊做討論: 1. 三個原子在同方塊:只要每個方塊都有各自計算到即可 2. 二個原子在同方塊:當計算起點在 two-body 所在的方塊時,會受到作用力的原 子只會在其周圍 26 個方塊中,所以只要依序計算 26 個方塊,並使每一個方塊 都當個起點,所有的排列組合都能計算到。對稱計算也可用在此,但是只是增 加複雜度,對於減少計算量並無貢獻。 3. 三個原子在異方塊:在此我們使用對稱的觀念來減少方塊的計算量,並有系統
地和不重覆地去取方塊來運算。
圖 3-5a three-body 三維空間分割圖(上層)
圖 3-5b three-body 三維空間分割圖(中層)
圖 3-5c three-body 三維空間分割圖(下層)
3.4.2 節的原理和用法,將正立方體對稱簡化後,可得圖 3-5a、3-5b 和 3-5c,上中下三 層共 63 個方塊,灰色區域在黑色方塊周圍(距起點一個作用力有效距離),透明區域在 灰色區域外圍(距起點二個作用力有效距離),以黑色方塊為計算起點,有系統地和不 重覆地去取該被計算的方塊(作用力有效範圍內),選取和黑色距離最遠的方塊來分開 討論: a. 最遠的方塊在透明區域:構成鍵結的三個原子中,有兩個原子分別在黑色和透 明方塊中,第三個原子如果要和另二個原子鍵結,勢必都在另二個原子的作用 有效範圍內,因此,對黑色周圍的方塊和透明周圍的方塊取交集,即為第三個 原子所在的方塊(不只有以上三圖的方塊列入考慮),數出在每個不同的透明方 塊下,有多少種選取方式。 b. 最遠的方塊在灰色區域:把灰色和黑色方塊拿出來,以三方塊是否在同一平面 分開來討論。 圖 3-6a three-body 三維空間分割圖(灰色區域上層) 圖 3-6b three-body 三維空間分割圖(灰色區域下層) 圖 3-6a 和圖 3-6b 為圖 3-5 的灰色和黑色方塊部分,共上下兩層。
z 三方塊不在同一平面上:代表著至少有一個方塊在上層,選取距黑色最遠的方塊 進行討論。見圖 3-6a,若最遠的方塊是四個角落的透明方塊之一,取透明和黑色 周圍方塊的交集處,即為第三個方塊的所在;若最遠的方塊是灰色方塊之一,第 三個方塊則是和灰色、黑色都有面接觸的方塊,而非周圍方塊的交集(若取交集將 重覆運算);淡色是最遠方塊的情形不需討論,已被上兩種選取方法所包括。 z 三方塊都在同一平面上:見圖 3-6b,要包含黑色方塊且不重覆選取,則可以得到 四個組合,見圖 3-7a、3-7b、3-7c 和 3-7d。 圖 3-7a 四排列組合之一 圖 3-7b 四排列組合之二 圖 3-7c 四排列組合之三 圖 3-7d 四排列組合之四 選取可能需要計算的方塊,過程雖然繁瑣且容易出錯,但是卻能有效的降低計算所需 的時間。
3.5 力的計算(compute forces)
原子間的作用力由函數推得,而為了減少計算量和誤差,使用了無因次化,詳細 內容請見 2.2 節、2.3 節和 2.4 節。 本篇論文所使用的函數是 S.W.位能函數和 B.H.位能函數,是為了矽在固液兩種狀態都能使用而推導出來的公式,詳見 3.3 節。 從 S.W.位能函數所推得作用力:
(
)
1 2 2 2 ( 1) exp[( ) ] p ij p ij ij ij ij A Br df F ABpr r a dr r a − 1 − − ⎛ − ⎞ ⎜ ⎟ = − = − − − − ⎜ − ⎟ ⎝ ⎠ − (3.13) 1 1 3 ( ,ij ik, jik) ( ) exp[ ( ij ) (ik ) ] ij ik jik F h r r r a r a r r r θ λ γ θ γ − − ∂ ∂ ∂ = −∇ = − + + − + − ∂ ∂ ∂ 2 1 (cos ) 3 jik θ × + (3.14) 從 B.H.位能函數所推得作用力: 2 2 2 2( )ij {( 1exp( 1 ) 2exp( 2 )) c( )}ij d d F f r A r A r f r dr dr λ λ = − = − − + − (3.15) 2 3 1 1 12 1 13 1 ( , , ) ( ) {[ ( ) ( )(cos ) 3 ij ik jik jik ij ik jik F h r r B r r r r r θ ψ ψ θ θ ∂ ∂ ∂ = −∇ = − + + × + ∂ ∂ ∂ 3 2 2 12 2 13 12 13 1 ( ) ( )(cos ) ] ( ) ( )} 3 jik c c Bψ r ψ r θ f r f r + + × (3.16) 2 ( ) 2 exp( ) i ij i i d r r drψ = − α − rα (3.17) 2 1 ( ) exp[( ) / ] {1 exp[( ) / ]} c ij c c d f r r r r r dr μ μ μ − = − − × + − (3.18) 力矩 M 等於位能 U 對角度θ的微分: dU M dθ = (3.19) 作用力 F 等於力矩 M 除以力臂 r: M F r = (3.20) 將(3.19)式代入(3.20)式,可得作用力 F 和角度θ的關係: M du F r rdθ = = (3.21) 由(3.21)式可推得(3.14)式和(3.16)式中對角度偏微分項為 jik r θ ∂ ∂ 。 以上幾式的作用力推導,只能算出力的大小,無法從數值中取得力的方向,因此, 必須由原子的相對位置去推得力的方向。由 two-body potential 所算出的力,其方向為 兩原子共線方向;由 three-body potential 所算出的力,受力方向可分為兩種,見圖 3-8,一種為兩原子共線方向 AB、AC(函數對長度的微分所得),另一種與 AB、AC 垂直且 在平面上的力矩 M(函數對角度的微分所得),後者可採用向量外積去計算出確切的受 力方向。
圖 3-8 three-body 力矩向量圖
見圖 3-8,首先,由向量 AB 和向量 AC 外積,可推得垂直於 ABC 平面的向量 AB×AC,再由向量 AC 和 AB×AC 外積,可推得垂直 AC 且在 ABC 平面上的向量 AC×(AB×AC),此向量為圖中右邊力矩所推算的力方向,同理,可得垂直 AB 且在 ABC 平面上的向量(AB×AC)×AB,為圖中左邊力矩所推算的力方向。
3.6 積分法-蛙跳法(Leapfrog methods)
蛙跳法為最簡單的積分法,雖然精確度不及其他積分法,但是計算量少,詳細方 法及理論請見 2.7 節,本文所採用的是計算較快的第一種方法,只採用(2.17)和(2.18) 式。A
B
C
A
J
B AC
×
M
M
(
)
AC
×
AB AC
×
JJJG
JJJG JJJG
JJG JJJG
(
J
A
JJG JJJG
B AC
×
)
×
JJJG
AB
3.7 邊界條件-週期性邊界(Periodic boundary condition)
模擬的模型是用分子束法沉積矽,因為基底的垂直方向(模擬中為 Z 方向)要進行 矽原子射入基底的模擬,所以沒有週期性邊界,在基底平面方向(模擬中為 XY 平面) 才有週期性邊界的應用,詳細用法及理論請見 2.6 節。3.8 溫度控制-定溫(Constant temperature)
如果只依照上述的程序進行模擬,射入的矽原子將使單結晶矽的基底融化,若沒 有對基底做溫度控制,根據能量守恆定律,熔融狀態的矽不會結晶,失去了原有的目 的,所以為了結晶成矽,必須對基底做溫度的控制,使得融化的矽在融點以下,慢慢 地凝固成固體。若由程式表現上述的行為,則是把部分原子(溫控層的原子)的平均速 度控制在由溫度換算而來的值(見(2.14)式),藉由修正溫度控制層的整體速度,來達成 最簡單的整體溫控。溫控的頻率,也就是修正原子速度的頻率不可太慢,否則會造成 溫控層原子速度劇變,這在物理上並非定溫,還會讓模擬出現錯誤;而修正的頻率也 不可太快,會增加多餘的計算量,本文以 10 個取樣時間(time-step)修正一次溫控層原 子的速度。第四章
矽沉積模擬
本章節以前兩章所敘述的理論為基礎,用 C 語言來寫程式碼,模擬用分子束沉積 形成單結晶矽的過程,模擬的模型是以單結晶矽為基底,給予要沉積的矽原子一個向 下的速度,每經過一段時間沉積一顆,模擬矽沉積的過程。矽原子間的交互作用,則 是用 S.W.和 B.H.位能函數來描述,此兩位能函數都是為了符合現實中矽的作用力,所 設計出來的,藉由比較兩者的結晶情形,討論兩函數的利弊。4.1 使用相同單位
本論文使用的是已無因次化後的單位(S.W.位能函數),詳細理論請見 2.4 節和 3.1 節,目的是為了方便運算和減少複雜度。由於另一個要比較的 B.H.位能函數,其單位 和 S.W.位能函數不同,所以要把這兩個位能函數的各項參數,都換算成同一單位,本 論文是以 S.W.位能函數為主,換算 B.H.位能函數的參數。位能函數詳見 3.3 節。 S.W.位能函數的參數: 7.049556277 A= 0.6022245584 B= 4 p= 1.8 a= 21 λ= 1.2 γ = B.H.位能函數的參數(換算後): 1 65.65600238 A = 2 49.38718653 A = − 1 6.012213915 B = 2 0.310106144 B =1 2.282877843 λ = 2 1.846608808 λ = 1 1.331921039 α = 2 1.40106795 α = 1.886657296 c r = 0.148946589 μ =
4.2 位能函數分析
把 S.W.和 B.H.位能函數換算成同單位後,圖解分析此兩函數。 圖 4-1 是 two-body potential 圖,S.W.的位能變化稍大,B.H.則較為平滑,從此可 以推斷出 S.W.的收斂速度較快。觀察位能的最低點,即兩原子最穩定的距離,在大約 1.1225σ時,S.W.位能值為最小;在大約 1.066σ時,B.H.位能值為最小。然而,矽原 子單鍵長度大約為 2.35Å,換算後大約是 1.1217σ,很明顯地,B.H.位能函數和實驗所 得的參數相差過大且短的多,在計算已結晶的基底時,僅有吸引力而無排斥力,因此, 本論文將對 B.H.位能函數稍作修改,使其能符合真實情況。 圖 4-2 是修改後的 B.H.和 S.W. two-body potential,改變的參數為 , 圖中兩函數穩定位置的原子距離差不多,B.H.略大於 S.W.,是為了使作用力距離較長 的 B.H.函數更穩定。 1 70.25600238 A = 圖 4-3 是 three-body potential 圖,S.W 圖形較 B.H.陡峭,位能變化較大,可推得 收斂速度較快。觀察位能的最低點,B.H.和 S.W.位能函數的角度平衡位置相同,大約 為 109.47 角度,剛好為鑽石結晶結構的夾角,所以不需要對函數多加變更。0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 -1 0 1 2 3 4 原子距離(σ) 位 能 (ε) S.W. B.H. 圖 4-1 S.W.和 B.H.的 two-body potential 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 -1 0 1 2 3 4 原子距離(σ) 位 能 (ε) S.W. B.H. 圖 4-2 S.W.和 B.H.(已改變參數)的 two-body potential
0 20 40 60 80 100 120 140 160 180 0 0.2 0.4 0.6 0.8 1 1.2 1.4 夾角(角度) 位 能 (ε) S.W. B.H. 圖 4-3 S.W.和 B.H.的 three-body potential
4.3 參數設定
基底的設定詳見 3.2 節,矽的單鍵鍵長為1.12166σ =2.35Å,初始溫度設定為 ,晶格方向為(1,1,1)向 z 軸方向,基底由六層矽原子組成,每層10 顆原子,共 480 顆原子。最底部兩層為固定層,雖然列入作用力的計算,但是不會移 動;固定層上兩層為溫控層,溫度控制在 ,詳見 3.8 節;溫控層上兩層為緩衝 層,使原子的能量轉換不會過快,詳見 3.2 節。 * 0.03T =756K 5 8 80 × = * 0.03T 入射原子的設定詳見 3.2 節,位置設定在固定範圍高度的隨機位置,溫度設定為 (熔點以上),由(2.16)式和駁回採樣法(Rejection sampling),給每個入射原子某個 範圍內的隨機速度,每顆入射原子的間隔為 200 個取樣時間,共 700 顆入射原子,取 樣時間為 。 * 0.1T * 1 0.02t =1.53268 10× − s 參考 3.3 節和 4-2 節,設定 S.W.的作用力範圍為1.8σ ,B.H.為 2.1σ 。4.4 模擬流程
Set parameters Initial coordinates Initial velocities for each atom
Cell subdivision
Compute forces
Leapfrog method Periodic boundary condition Constant temperature Add a atom every X time-step
End
≦ XX time-step
≧XX time-step
≦ XXX deposition atoms Count ≧ final count
≧ XXX deposition atoms 圖 4-4 矽沉積模擬流程圖 圖 4-4 為本論文的模擬流程圖,不同於圖 3-1,在入射原子前先計算基底至穩定的 狀態,模擬其為真實情況-表面重建[12];在入射全部的原子後,再計算系統至穩定 狀態,得到最穩定的晶體結構。
4.5 模擬結果
圖 4-5 是採用 S.W.位能函數的矽沉積模擬,除了最上層原子因為數量不足,而無 明顯的結晶層產生外,其他矽所形成的結構非常清楚,從最底層可以清楚地依序往上 辨識,固定層、溫控層和緩衝層,還有沉積原子形成的六層結晶層。由圖 4-6 動能變 化圖可知,系統已經穩定,所以圖 4-5 為最終的晶體結構。比較圖 4-5 和圖 4-7 後,可 看出 S.W.的缺點-無法辨別鑽石晶格結構和纖維鋅礦(wurtzite)晶格結構,(111)方向的 排列分別為 AaBbCcAa…和 AaBbAaBa…,詳見圖 4-8,最上層為灰色原子,在其之下 皆為透明原子,原子由下往上堆積,圖 4-8d 和圖 4-8e 為鑽石晶格結構和纖維鋅礦晶 格結構的主要差異。圖 4-5 S.W.位能函數結晶圖 0 500 1000 1500 2000 2500 3000 100 150 200 250 300 350 400 450 500 550 600 取樣時間(100t*) 系 統 動 能 (ε) 原子沉積結束 原子沉積開始 系統穩定過程 圖 4-6 動能變化圖(S.W.)
圖 4-7 完全結晶矽(鑽石結構)
圖 4-8a 鑽石和纖維鋅礦晶格結構上視圖(Aa 層)
圖 4-8c 鑽石和纖維鋅礦晶格結構上視圖(AaBb 層) 圖 4-8d 鑽石晶格結構上視圖(AaBbC 層) 圖 4-8e 纖維鋅礦晶格結構上視圖(AaBbA 層) 圖 4-9 是 B.H.位能函數(已改變參數)的結晶圖,很明顯地知道結晶並不完美,除了 不動的固定層外,溫控層和緩衝層難以辨別。圖 4-10 為圖 4-9 的系統中的一段區塊, 區塊中原子間的距離大都維持在定值,但是沒有排列成漂亮的鑽石晶體結構,由此可 知,two-body potential 確實穩定,而 three-body potential 沒有到達穩定值,使系統不是 最低的能量形式-單結晶鑽石結構,由此可見,目前的系統只是亞穩態,並非真正的 穩態,可能有幾點原因造成此種結構。
1. 入射原子不均勻:由於入射原子的位置是隨機的,所以位置不可能完全地平均 分佈,力的分佈也不會平均,也就不會形成清楚的晶格結構。 2. 平衡時間不足夠:系統尚未到達最穩定的狀態,原子還在緩慢地向最穩定的位 置移動中,所以無法呈現清楚的晶格結構。 3. 數據誤差:除了參數的錯誤外,經過無數次的計算所累積出的誤差值,可能對 最後的晶格結構造成影響。 圖 4-9 B.H.(已改變參數)位能函數結晶圖
第五章
結論與未來計劃
5.1 結論
論文的兩大重點,three-body 的空間分割法和兩種不同位能函數的結果比較,分 別可在 3.4 節和第四章讀到,前者以現有的 two-body 空間分割法為理論基礎推導而 成,後者從換算相同單位、比較位能函數圖,到觀察最後的矽晶體結構,比較兩者的 不同。在計算 three-body potential 時,該如何和 two-body potential 一樣,可藉由作用力 和反作用力的觀念,來簡化計算的步驟,現有的書籍資料未提及這方面的方法,因此, 本論文利用以上觀念,配合有系統的選取,來簡化 three-body potential 的計算,除去 大幅無用的步驟,第四章 S.W.位能函數的模擬結果,漂亮的晶體結構已證明此方法正 確無誤。 在第四章中,兩位能函數圖的比較,B.H.函數圖形收斂慢,作用力範圍較長,但 是能量的最低點卻不符合矽的參數,在修改參數後,可預期其將比 S.W.函數更接近真 實矽的描述,然而在兩者結晶圖的比較後,卻是 S.W.函數形成排列整齊的結晶,B.H. 函數則完全無結晶。比較三種可能的原因: 1. 入射原子不均勻:入射原子足夠多,上幾層可能不均勻,但是中間幾層則是均 勻分佈,不均勻的情況可以不考慮 2. 平衡時間不足:B.H.函數的作用力範圍廣,平衡時間比 S.W.函數需要的更多, 模擬的時間可能不足 3. 數據誤差:S.W.和 B.H.函數各有不到五顆的原子跑離系統,兩函數累積出的誤 差差不多,而 B.H.函數的能量最低點不符,故 B.H.函數參數可能有誤。
5.2 未來計劃
從 5.1 節可知,B.H.函數雖然較貼近真實的矽原子,但是收斂速度和晶體結構不 夠完美,期許未來能持續的研究參數關係,將其改寫成快速又正確的位能函數,期望 能藉由此兩函數模擬在不同環境參數下,矽晶體的變化,以及加入不同於矽的物質, 觀察晶體形成的過程與結果。
參考文獻
[1] J. W. Matthews, Epitaxial Growth, Academic, New York, 1975.
[2] Frank H. Stillinger and Thomas A. Weber, Computer Simulation of Local Order in Condensed Phases of Silicon, Physical Review B, Volume 31, 1985.
[3] R. Biswas and D. R. Hamann, New Classical Models for Silicon Structural Energies, Physical Review B, Volume 36, 1987.
[4] S. A. Rice and P. Gray, The Statistical Mechanics of Simple Liquids, Wiley-Interscience, New York, 1965.
[5] C. A. Croxton, Liquid State Physics-A Statistical Mechanical Introduction, Cambridge University Press, London, 1974.
[6] J. P. Hansen and J. R. McDonald, Theory of Simple Liquids, Academic, New York, 1976.
[7] M. Schneider, Ivan K. Schuller and A. Rahman, Epitaxial Growth of Silicon:A Molecular Dynamics Simulation, Physical Review B, Volume 36, 1987.
[8] Farid F. Abraham and Jeremy Q. Broughton, Pulsed Melting of Silicon (111) and (100) Surfaces Simulated by Molecular Dynamics, Physical Review Letters, Volume 56, 1986.
[9] E. T. Gawlinski and J. D. Gunton, Molecular Dynamics Simulation of Molecular Beam Epitaxial Growth of the Silicon (100) Surface, Physical Review B, Volume 36, 1987. [10] Mark D. Kluge and John R. Ray, Amorphous Silicon Formation by Rapid Quenching:
A Molecular Dynamics Study, Physical Review B, Volume 36,1987.
[11] R. Biswas, Gary S. Grest and C. M. Soukoulis, Molecular Dynamics Simulation of Cluster and Atom Deposition on Silicon (111), Physical Review B, Volume 38, 1988. [12] H. J. Gossmann and L. C. Feldman, Initial Stages of Silicon Molecular Beam
Epitaxy:Effects of Surface Reconstruction, Physical Review B, Volume 32, 1984. [13] D. C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University
press, 1995.
[14] D. C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University press, second edition, 2004.
[15] Peter Atkins, Julio de Paula, Atkins’ Physical Chemistry, 7th edition, Oxford, p.818 ~ p.819.