液態鉛的動力學性質研究

全文

(1)附件一封面格式樣本. 國立交通大學物理研究所碩士論文. 液態鉛的動力學性質研究 The Study of The Dynamic Properties of Liquid Pb. 研究生：黃楷智指導教授：吳天鳴. 教授. 中華民國九十四年六月.

(2) 附件二書名頁樣本. 液態鉛的動力學性質研究 The Study of The Dynamic Properties of Liquid Pb. 研究生：黃楷智. Student：Kai-Chih Huang. 指導教授：吳天鳴教授. Advisor：Dr. Ten-Ming Wu. 國立交通大學物理研究所碩士論文. A Thesis Submitted to Institute of Physics College of Science National Chiao Tung University in partial Fulfillment of the Requirements for the Degree of Master in Physics June 2005 Hsinchu, Taiwan, Republic of China. 中華民國九十四年六月.

(3) 附件五橫式中文摘要格式樣本. 液態鉛的動力學性質研究. 指導教授：吳天鳴教授. 學生：黃楷智. 國立交通大學物理研究所碩士班. 摘. 要. 本篇論文致力於利用分子動力模擬及瞬間正則模探討液態鉛的各項動力學性質，文中分別探討兩種不同理論基礎的虛擬液態鉛位能，並比較一系列的動力學性質結果以了解兩者的差異，包括粒子徑向分布函數、靜態結構因子、動態結構因子、粒子速度相干函數、瞬間正則模態密度以及參與比例，並試圖尋找瞬間共振模，最後的結論為如果在虛擬液態鉛位能中的核心排斥力部分有反曲點的話，就可以發現瞬間共振模的存在，且瞬間共振模的態密度分布會隨著溫度的升高而增加。. i.

(4) 附件五橫式中文摘要格式樣本. The Study of Dynamic Properties of Liquid Pb. Student： Kai-Chih Huang. Advisors： Dr. Ten-Ming Wu. Institute of Physics National Chiao Tung University. ABSTRACT. The main purpose of this thesis is to explore several dynamic properties of liquid Pb, using the molecule dynamic simulation and the instantaneous normal mode analysis. Two potentials of liquid Pb based on different theories are discussed, and the results of their dynamic properties, including radial distribution function, static structure factor, dynamic structure factor, velocity autocorrelation function, density of state and participation ratio, are compared in order to investigate the differences. Attempt is also made to find out the instantaneous resonant modes in liquid Pb in this study. It is concluded that the existence of instantaneous resonant modes can be found in the case that there is a reflection point in the repulsive core of the liquid Pb potential. It is found out that the density of states of instantaneous resonant mode increases with the increase of temperature.. ii.

(5) 謝誌回首這兩年，要感謝的人實在太多，如有疏漏，還望見諒。在研究的專業工作上，最感謝的當然是指導教授吳天鳴老師的耐心指導與諄諄教誨，讓我得以進入研究工作的領域，並了解到研究工作的態度與精神，也要感謝昆憲學長及同學凱驗在程式及物理上對我的幫助以及啟發，讓我能夠順利的在兩年之中完成這些工作。而在工作之餘，也很謝謝所上的所有同學，大家一起笑鬧，一起打球、聊天，讓這兩年除了研究工作之外，還有許多的歡笑，並結交許多知心的朋友。而國高中及大學的好友們也是給予我極大精神支柱的重要角色，感謝你們。而一直在背後默默給予我支持的，當然是最重要的家人，爸媽、老姐，還有乾爹乾媽總是在背後一直給我鼓勵和支持，讓我在遇到挫折與打擊時，依然能夠支持下去，繼續前進。而在我心中擁有極重要地位的玉婷，這兩年來我的喜怒哀樂，她是最清楚的，謝謝玉婷用心及耐心的陪伴，而我也終於完成了這個階段的工作，接下來就要邁向新的生活，這兩年的點點滴滴，我銘記在心。. iii.

(6) 目錄中文摘要 ……………………………………………………………………. i. 英文摘要 ……………………………………………………………………. ii. 謝誌 …………………………………………………………………………. iii. 目錄 …………………………………………………………………………. iv. 圖目錄 ………………………………………………………………………. vi. 符號說明 ……………………………………………………………………. viii. 第一章、緒論 ……………………………………………………………..... 1. 第二章、理論 ………………………………………………………………. 3. （一）Pair potential ………………………………………………. 3. （二）粒子徑向分佈函數(radial distribution function) g(r) 和. G 靜態結構因子（static structure factor） S k ………….... (). 5. （三）粒子速度相干函數 (velocity autocorrelation function) C(t) …………………….. 6. G （四）中間散射函數 (intermediate scattering function) F k , t 和. ( ). G 動態結構因子 (dynamic structure factor) S k , ω. ( ). ….……. 7. （五）瞬間正則模 INM(instantaneous normal mode) 與態密度 D（ω）（ density of state ） ………………...…. 8. （六）參與數 R Nα （ participation number ）、參與比例 PN (ω ) （ participation ratio ）、化約後的參與比例 Sα （ reduced participation ratio ）與瞬間共振模 IRM(instantaneous resonant mode)… iv. 10.

(7) 第三章、分子動力模擬 …………………………………………………….. 12. 第四章、結果與討論 ……………………………………………………….. 14. 第五章、結論 ……………………...………………………………………... 18. 參考文獻 …..……………………...………………………………………..... 19. v.

(8) 圖目錄圖一. 位能 B 在溫度為 623K 的靜態結構因子 S(k)和 k 的關係圖……... 20 圖二. 位能 B 在溫度為 1170K 的靜態結構因子 S(k)和 k 的關係圖……. 21. 圖三. 位能 A 在三個不同溫度下和距離 r 的關係圖………….…………. 22. 圖四. 位能 A 在 1173K 和位能 B 和距離 r 的關係圖……………………. 23. 圖五. 粒子徑向分佈函數 g(r)和距離 r 的關係圖…………………………. 24 圖六. 靜態結構因子 S(k)和 k 的關係圖………………..…………………. 圖七. 速度相干函數 C(t)和時間 t 的關係圖………………………………. 25. 26. (G ). 圖八. 動態結構因子 S k , ω 示意圖 ………………………………………. 圖九. 鉛假位能函數化約單位取值示意圖……………………...…………. 27. 28 位能 A 的粒子徑向分佈函數 g(r)和距離 r 的關係……...………….. 29. 圖十一位能 B 的粒子徑向分佈函數 g(r)和距離 r 的關係圖……..……..…. 30. 圖十. 圖十二位能 A 在 1173K 和位能 B 在 1170K 的粒子徑向分佈函數 g(r)和距離 r 的關係圖……………………………………………………. 31. 圖十三位能 A 和位能 B 分別在 1173K 和 1170K 的靜態結構因子 S(k)和實驗值的比較………………………………………………………. 32. 圖十四位能 B 在 623K 的靜態結構因子 S(k)和實驗值的比較...…………. 33. G 圖十五位能 A 在 1173K 和位能 B 在 1170 的動態結構因子 S k , ω 和實驗. ( ). 值的比較……………………………………………………………. 34. 圖十六位能 A 的速度相干函數 C(t)和時間 t 的關係圖………..…………. 37. 圖十七位能 B 的速度相干函數 C(t)和時間 t 的關係圖…………..………. 38. vi.

(9) 圖十八位能 A 在 1173K 和位能 B 在 1170K 的速度相干函數 C(t)和時間 t 的關係圖…………………………………………………….……. 39. 圖十九位能 A 的 INM 狀態密度與頻率ω的關係圖…………………..…. 40. 圖二十位能 B 的 INM 狀態密度與頻率ω的關係圖………..……………. 41. 圖二十一位能 A 的參與比例 PN(ω)與頻率ω的關係圖……………………. 42. 圖二十二位能 A 的參與比例 P2N(ω)/ PN(ω)與頻率ω的關係圖……………. 43. 圖二十三位能 B 的參與比例 PN(ω)與頻率ω的關係圖…………………...... 44. 圖二十四位能 B 的參與比例 P2N(ω)/ PN(ω)與頻率ω的關係圖…………....... 45. 圖二十五位能 B 的 IRM 狀態密度與頻率ω的關係圖……………..……….. 46. 圖二十六位能 A 在 1173K 和位能 B 在 1170K 的位能及位能二次微分和距離 r 的關係圖…………………………………….…………….... vii. 47.

(10) 符號說明 ρ：系統密度 N：系統粒子總數 k：波向量 r：相對距離 kB：波茲曼常數 T：系統溫度 t：單位時間 ε：最穩定之能量 σ：粒子直徑 cs：系統聲速 ω：頻率 D： Hessian matrix INM：瞬間正則模 IRM：瞬間共振模. viii.

(11) 第一章緒論在液體、非晶態固體 (armophus) 等系統中，因為系統的無序性，不具有週期性晶格結構，所以很難長時間的追蹤其運動的軌跡，然而在短時間我們可以利用一些技巧，像是利用簡諧近似 (Harmonic. approximation) 等方法，可以描述系統的瞬間動力性質，上述的方法稱為瞬間正則模 (instantaneous normal mode 以下簡稱 INM) 分析. [6][7] 。而在過去的研究中 [5] ，發現在非晶態固體或是超冷液體 (supercoolede liquid) 中，有一種低頻的準侷限 (qusailocalized) 正則模 (NM)，稱之為共振模 (resonant mode, 以下簡稱 RM)，在系統中扮演著相當重要的角色，近年來，研究的觸角慢慢延伸到液體系統 [2][3][4]，因此，瞬間共振模 (instantaneous resonant mode 以下簡稱 IRM) 在液體中所扮演的角色漸漸成為一個很值得探討的主題，而 INM 的方法就非常適合用在對液體的研究上，而瞬間一詞旨在強調無序系統的極短時間有效性。要得到 INM 的結果需求解一系列瞬間系統組態的 Hessian. matrix( 或 dynamic matrix)，理論模擬上求得一系列系統組態的方法有二，其一為蒙第卡羅 (Monte Carlo ，以下簡稱 MC) ，另一個方法為分子動力模擬 (Molecule Dynamic simulation ，以下簡稱 MD. simulation) ，在本篇論文中採用後者。在理論模擬上要研究液體中的各項動力學性質， MD simulation 是相當有效的工具，而在實驗上，大都採用非彈性中子散射實驗 (Inelastic Neutron Scattering, 以下簡稱. INS) 或是非彈性 X 光散射實驗 (Inelastic X-ray Scattering, 以下簡稱 IXS) 來探討液體的動力學性質，為了更加了解模擬實驗結果和真實世界的液態系統之間的差異，文中會將模擬的結果和實驗的結果作一個 -1-.

(12) 比較分析。在之前的研究中顯示 [2]，在某些液態金屬 (liquid metal) 中，如鎵金屬 (Gallium) ，由於在系統位能的排斥力部分 (repulsive core) 有一個反曲點，因此造成 IRM 的發生，本論文將對象放在液態鉛 (liquid lead) 上，想知道在真實世界的液態鉛是否會存在 IRM，而液態鉛的虛擬位能也是一個探討已久的課題 [1] ，本論文中我們分別探討兩個不同的虛擬鉛位能系統，從結果可以發現，如果系統位能在排斥力部分具有反曲點的話，則我們就可以找到 IRM，這點和前面有關液態鎵的研究結論是一致的。下面第二章中解釋我們所使用的兩種不同的虛擬鉛位能系統和一些我們要探討的動力學性質理論基礎。第三章闡述我們所使用的模擬實驗方法和一些詳細數據。第四章則將模擬的結果配合圖表表示並分析比較。第五章則是最後的結果討論與結論。. -2-.

(13) 第二章理論（一） Pair potential 本篇論文中採用兩種不同的液態鉛虛擬位能，分述如下：（ 1 ）假位能（ Pseudopotential ）在真實世界中，液態金屬的組成包括離子及價電子，故其位能應包含三項，分別為系統中離子與離子的作用力、離子與電子之間的作用力還有電子與電子之間的作用力。在電腦模擬時，需要系統的位能才能進行實驗，但我們很難精確的掌握上述三項的每一個作用力，特別是系統粒子數目越多時，模擬位能的工作就變得更加複雜，因此有理論學家將液態金屬中的電子看成均勻的電子海，而可將三項作用力簡化成只有離子與離子之間的作用力，這也就是我們所採用的第一個位能－鉛假位能（ Lead Pseudopotential ），應用 first-principles. generalized energy-independent nonlocal model-pseudopotential (GEINMP) theory ，液體金屬的假位能函數可由下式得到 [12]. V (r; Ω ) =. 2 Z eff. r.  2 1 −  π.  ( ) ( ) Ω dkF k ; sin rk / k  ∫0 . ∞. 其中 Ω 是原子體積， Zeff 是等效離子之原子量， k 是波向量， F (k ; Ω ) 是. exchange-correlation factor 。理論上，此為一個虛擬的位能，但是卻可以有效的將真實世界中液態金屬裡離子間行為充分的表現出來。以下簡稱假位能為位能 A 。. （２） Model potential 模擬作用位能的方法有很多種，本篇論文中所採用的第二個虛擬位能引述自 [8] ，利用一個固定的模型，但其中總共有７個可以調整 -3-.

(14) G 的參數，藉由和實驗的靜態結構因子 S k 反覆作比較進而調整出能夠. (). 表達真實世界中液態鉛行為的位能，但能夠和實驗資料吻合的的參數可能不只一組，並且在得到此位能後，利用計算系統擴散係數並和實驗值作比對來當作一個確認的條件，藉由上述兩個物理量來確認此位能的正確與否。此位能較為特別的是不受系統的溫度和密度影響，雖然是針對鉛在 623K 時的狀況模擬出來的，但用在 1170K 時卻可以和. G 實驗在 1170K 的靜態結構因子 S k 一樣有很好的吻合，從圖一和圖二. (). G 中可以看到此位能所模擬出的兩個不同溫度下之 S k 和實驗值的比. (). 對，除了在 k 值大約 1 A - 1 左右的地方有比較明顯的誤差以外，大體上都和實驗值很符合，而我們採用的各項參數都跟 [8] 裡的一樣，一個不同的地方在 [8] 中使用的粒子數目是 16384 顆而我們接下來的一系列的電腦模擬中採用的粒子數目是 500 顆。此虛擬位能的形式如下：. V (r ) = V1 (r ) + V2 (r ) + V3 (r )   b1 12 −1 a1   exp (r − c1 ) , V1 (r ) =   r  0 , . (. ). (. ).   b2  −1 a2  − c2  exp (r − b2 ) , V2 (r ) =   r  0 , . V3 (r ) = a3r −3 cos(2 K F r ) 其中各項參數如下 a1 = 102.5mev a2 = 90.0mev 3 a3 = 418.3mev × A b = 3.284A 1. -4-. for. r < c1. for. r ≥ c1. for. r < b2. for. r ≥ b2.

(15) b 2 = 4 . 83 A c = 5 . 72 A 1. c 2 = 0 .5 K F = 1 . 5417 A -1 以下簡稱為位能 B。首先我們比較上述兩個不同的虛擬鉛位能系統，從圖三中可以看出位能 A 會隨著溫度的變化而有些微的不同，當溫度越高時，能量的最低點會越低，但能量最低點之後的 Friedel part 則沒有太大的差異。而位能 B 的特性是不會隨著溫度的改變而不同，可以用在兩個不同的溫度上，一為接近熔點（約 600K）的 623K，另一個溫度為明顯高出許多的 1170K。從圖四中可以看到位能 A 和位能 B 有明顯的差別，首先，位能 B 在核心排斥力（ repulsive core）的部分較位能 A 弱，從位能零點截距σ的取值可知位能 B 系統中的虛擬液態鉛粒子直徑較位能 A 系統中虛擬液態鉛粒子直徑大，而兩個位能中之能量最低點的位置上，位能 B 較位能 A 明顯較遠離零點，顯示在位能 B 中的兩粒子間達到最穩定狀態的距離較遠，而過了能量最低值後的 Friedel part 的震盪幅度則明顯較位能 A 小。. （二）粒子徑向分佈函數 (radial distribution function) g (r ) 和靜態結構因子 (static structure factor) S(k) 粒子徑向分佈函數 (radial distribution function) g (r ) 是在描述系統中所有粒子在空間中的分佈情況，選擇任意一顆粒子，並以自身為座. G 標的中心，觀察其它粒子在空間中隨著 r 變化的分佈情形，空間中粒子的總數 N 可表示成： ∞. ∫ ρ ⋅ g (r ) ⋅ 4πr. 2. G ⋅ dr ≈ N. （１）. 0. -5-.

(16) 其中 ρ 代表系統的密度，從 (1) 式中的 ρ ⋅ g (r ) ⋅ 4πr 2 說明了以中心粒子的. G G G 眼光所觀察到由 r 到 r + dr 之間的粒子數目，從 ( 圖五 ) 中我們可以看出，在空間中，系統粒子是呈現層次性的分佈的。圖中的第一個波峰代表最接近中心粒子的粒子數目，第二個波峰代表次接近中心粒子的 G 粒子數目，以此類推，當 r → ∞ ，則 g (r ) → 1 ，是表示中心粒子所感受. 到的系統密度已無太大的變化，因此可視為一常數。. G 靜態結構因子 (static structure factor) S k 是經由 INS 和 IXS 實驗. (). 所觀察到的結果，而我們可藉由上述的粒子徑向分佈函數 g (r ) − 1 做傅利葉轉換後得到：. [ ( )]. ∞ G G G G S k = ρ ⋅ ∫ [g (r ) − 1]⋅ exp i ⋅ k ⋅ r ⋅ dr + 1. (). （２）. 0. G 2π ( 圖六 ) 中 S k 的最大波峰位置 k max ，值對應到 ( 圖一 ) 中 g (r ) 振盪波 k max. (). G 波長。從圖中可看出，當 k → ∞ ，則 S k → 1。. (). （三）粒子速度相干函數 (velocity autocorrelation function) C (t ) 歸一化的粒子速度相干函數 (velocity autocorrelation function). C (t ) 的定義為： C (t ) =. v(t ) ⋅ v (0 ). （３）. v(0 ) ⋅ v (0). G 其中， v (t ) 是粒子在時間 t 的速度。在系統中，隨著時間的變動而粒. 子的速度也會跟著改變，粒子速度相干函數就是觀察當我們以一顆粒子為觀察中心時，此粒子隨著時間所發生的速度變化與時間ｔ＝０時速度的相關程度，從（圖七）中我們可以看到隨著時間的變化，粒子速度相干函數 C (t ) 會漸漸趨近於零，因為系統中粒子會和其他的粒子有所碰撞，隨著時間的演化，此時的粒子速度 v(t ) 會和最初速度 v (0 ) 的 -6-.

(17) 關聯性變得越來越低。另外我們會看到 C (t ) 會有負值，是因為當粒子在運動的時候，剛開始可能會受到週遭粒子分布的影響，而被暫時阻擋住，經過一段時間的運動後就可以突破此侷限。在 IXS 或 INS 的實驗中並無法直接得到粒子速度相干函數 C (t ) 的數值，大部分都是藉由電腦 MD simulation 得到。. G （四）中間散射函數 (intermediate scattering functions) F k , t. ( ). 與. G 動態結構因子 (dynamics structure factor) S k , ω. ( ). G S k , ω 為可以在 IXS 或 INS 實驗中觀察到的密度波動頻譜 (density. ( ). fluctuation spectrum) ，但在電腦 MD simulation 中，我們無法直接得 G 到動態結構因子 S k , ω ，必須要先求得中間散射函數 (intermediate. ( ). G scattering functions) F k , t ：. ( ). G 1 F k , t ＝ ⋅ ρ kG (t ) ⋅ ρ −kG (0) N. ( ). （４）. 其中 ρ kG (t ) 是密度波動 (density fluctuations) 的平面波： N. [. G G. ]. ρ kG (t ) = ∑ exp i ⋅ k ⋅ R j (t ) j =1. （５）. G 其波向量 (wavevector) 為 k ，將 (4) 式做時間的傅利葉轉換後就可得到. G 動態結構因子 (dynamics structure factor) S k , ω ：. ( ). ∞ G G 1 S k ,ω = ⋅ ∫ F k , t ⋅ exp(i ⋅ ω ⋅ t ) ⋅ dt 2π −∞. ( ). ( ). （６）. 根據理論和實驗結果，可在密度頻譜中得到三個勞倫茲分量（圖八），在 ω = 0 處的波峰稱為 Rayleigh 波峰， Rayleigh 波峰關係到系統的熱. -7-.

(18) 擴散，當系統溫度愈高，波峰就愈寬；而在 ω = ±c S k 處的波峰稱為. Brillouin 波峰，此分量關係到系統中的聲速 (c S ) ，此三個波峰的寬度都和 k 2 成正比，如果我們將 S (k , ω ) 的頻譜積分，範圍從頻率的負無限大到正無限大，則可以得到下列關係式：. G dω G G S k , = S k = F k ,0 ω ∫ 2π −∞ ∞. ( ). () ( ). （７）. G G 而經由理論計算所得到的 S k , ω 和實驗觀察所得到的 S exp k , ω 曲線，會. ( ). ( ). 因為兩項因素而有所差異：（１）平衡因子 (detail balance) 。（２）儀器解析度 (instrument resolution) 。兩者之間的關係我們用下列式子表示： G S exp k , ω = ∫. ( ). G =ω / k BT ⋅ S k , ω ⋅ R ω − ω ' ⋅ dω 1 − exp − =ω / k BT. ) (. (. ) (. ). （８）. (8) 式中積分內第一項即平衡因子 (detail balance) ，其中的符號 k B 為波茲曼常數， T 為系統溫度，第三項為儀器解析度 (instrument G resolution)，通常以高斯函數代替。而經由實驗得到的 Sexp k , ω 頻譜在. ( ). 呈現時經常會用對稱後的動態結構因子（ symmetrized dynamic. ~ G structure factor ） S k , ω 來表示，關係如下：. ( ). G ~ G S k , ω ＝ exp(− =ω / 2k BT ) × S exp k , ω. ( ). ( ). （９）. （五）瞬間正則模 INM(instantaneous normal mode)與態密度 D(ω ) （ density of state）有別於固體，在液體中，分子沒有固定的晶格結構，因此想要長時間追蹤液體分子的位置變化相當困難。因此，我們可以利用一個在固態物理中用到的技巧，對系統取極短時間的系統組態，然後利用簡 -8-.

(19) 協近似來探討其中的性質。就稱為瞬間正則模 (instantaneous normal mode)，以下稱 INM。 G G G G 首先定義系統組態為 R ≡ (r1 , r2, r3, ...) ，假設此系統在某一瞬間的總 G 位能為 V R0 ，經過很短時間的運動後，此時系統的組態發生變化，. ( ). 而我們可以利用簡協近似的方法得到此時的位能為：. G G G G G G G G G 1 G V Rt ≈ V R0 + ∇V R0 ⋅ Rt − R0 + Rt − R0 ⋅ ∇∇V R0 ⋅ Rt − R0 2 G G 其中 ∇ V R0 ≡ − F. ( ) ( ). ( )(. ) (. ). ( )(. ). (１０ ). ( ). G ∇∇V R0 ≡ D. ( ). G F 是瞬間作用力；. D. 稱作 Hessian matrix （ dynamic matrix ），為. 系統位能的二次微分矩陣，定義如下： Dij ;µν = −[Τij ] µν ; i ≠ j Dij ;µν = ∑ [Τ jl ] µν ; i = j l≠ j. V ' (rij ) G [Τij (rij )] µν = [ (Ι 3 − rîj rîj ) + V '' (rij )rîj rîj ] µν rij. 其中 i,j 為粒子的 index ， μ ,ν 表示 x 軸、 y 軸或是 z 軸， I 3 為 3×3 之單位矩陣。每一個瞬間系統組態，都有不同的作用力和 Hessian matrix ， G 2 求解此 Hessian matrix 是一個本徵問題，解後可得到本徵值 ω α ( R0 ). （ eigenvalue ）和本徵向量（ eigen vector ），求得的本徵值開平方後即為該模α的頻率，但頻譜並非全部都是實數，也有虛數的可能，表示該模並不穩定。經由電腦 MD simulation 求得一系列的瞬間系統組態後，求解 Hessian matrix ，將其本徵值作統計分析後，即可得到系統的歸一化態密度分布狀態 D(ω ) ，我們將虛數的頻率作圖在負數區段，歸一化態密度分布狀態 D(ω ) 定義如下： -9-.

(20) 1 3N. D (ω ) =. ∫ D (ω ) d ω. 3N. K. ∑ δ [ω − ω α ( R 0 )]. α =1. ，. =1. G. （六）參與數 R Nα （ participation number）、參與比例 PN (ω ) （ participation ratio）、化約後的參與比例 Sα （ reduced participation ratio）與瞬間共振模 IRM(instantaneous resonant mode) 解本徵方程式，我們可以得到本徵解（ eigenvalue ）還有本徵向量（ eigen vector ），本徵向量可以用來計算參與數 R Nα （ participation. number ）：  N G R N =  ∑ e αj  j =1  α. 4.    . −1. (１１ ). G 其中 e αj 為對第 α 個模的歸一化（ normalized ）後本徵向量之投影分. 量。參與數 R Nα 是表達平均粒子數的參與程度。參與比例 PN (ω ) 的定義為. R Nα ，其中 N 表示系統中粒子總數。 N. 化約後的參與比例 Sα （ reduced participation ratio ）：. Sα ＝ R Nα. Q Nα. (１２ ). 其中 Q Nα 的定義與參與數 R Nα 的定義很相似，唯一不同的是累加的計算中摒棄掉本徵分量中最大的一個分量。如果 Sα 趨近於１，則我們稱此模為延展模（ extented mode ），反之，若趨近於０，則稱為侷限模（ localized mode ）。而在 INM 中，有一些在低頻（ low frequency ）、介於延展模與侷限模之間的模，我們稱為瞬間共振模 (instantaneous resonant mode)，. - 10 -.

(21) 以下稱 IRM， IRM 是一種延展模與侷限模的混合（準侷限（ quasilocalized ）），為了和在低頻正數區段的 INM 區分，我們利用化約後的參與比例 Sα 來加以辨識。當 Sα ＜ 0.5 ，且在低頻的 INM ，我們就稱之為瞬間共振模 IRM。. - 11 -.

(22) 第三章方法分子動力模擬（ molecular dynamics simulation ）我們利用電腦虛擬一個立方體的盒子，裡面有 500 顆鉛粒子，並採用週期性邊界條件（ periodic boundary condition ）及 Verlet 演算法。系統是第二章（一）中所陳述的鉛假位能（ pesudopotential ），位能. A 。系統溫度有三個，分別是 823 Ｋ、 1023K 和 1173K ，然後在運算的過程中，需採用化約單位，如圖九，位能的最低能量深度為 ε，零點位能的距離截距為 σ，ε 和 σ 的值會隨著系統而有所不同，所有的 ε 和 σ 取值以及系統的密度 ρ 在表一中，粒子質量取值為１，而系統中所採用的各項化約單位定義如下：溫度 T* =. k BT. ε. ，密度 ρ * = ρσ 3 ，. 1.  mσ 2  2 1  ，頻率 ω 0 = ，每一步的時間單位 dt * = 0.001t 0 ，時間單位 t 0 =  t0  ε  其中密度的定義為每單位體積內的粒子數目（ particle/V），所有化約單位的最後取值在表二中。待系統到達我們所要的狀態之後，經過 10 ６步以確定系統達到平衡後讀取此時的系統組態，計算此時系統的. G 粒子徑向分佈函數 g (r )、靜態結構因子 S k 、粒子速度相干函數 C (t ) 以. (). G 及動態結構因子 S k , ω ，波向量的範圍為 1.8A - 1 到 6.5A - 1 之間的六個. ( ). K k 值。接下來每 3500 步收集一個系統組態，共收集 200 個系統組態，用來求解 Hessian matrix。另外也收集 1000 顆鉛粒子時的系統組態，化約單位都和上面的定義一樣，不同的是經過每 5000 步收集系統組態，總共收集 100 個系統組態。後來我們另外選定 [8] 中的虛擬鉛位能，即第二章（一）中所陳述的位能 B ，用一樣的方法進行分子動力模擬，系統的溫度為 623K - 12 -.

(23) 和 1170K ，一樣 10 ６步讀取此時的系統組態，然後再收集一系列的系統組態。 500 顆鉛粒子每 4000 步收集一個系統組態，總共收集 200 個，而 1000 顆鉛粒子每 5000 步收集一個系統組態，共收集 100 個。接下來呈現兩個不一樣的系統位能所求得的實驗結果。. 表一位能. 溫度. ε (×10 - 4 au). σ (A). ρ ( A-3). A. 823K. 1.083102. 4.278960388. 0.03375. A. 1023K. 1.153826. 4.279191616. 0.03320. A. 1173K. 1.241366. 4.279670306. 0.03276. B. 623K. 1.134769. 4.648919542. 0.03094. B. 1170K. 1.134769. 4.648919542. 0.02890. 表二位能. 溫度\化約單位. ρ*. A. 823K. A. t 0 (ps). dt *. 2.644111. T* 24.077815. 11.553888. 0.001. 1023K. 2.601334. 28.094539. 11.194953. 0.001. A. 1173K. 2.567807. 29.942273. 10.794063. 0.001. B. 623K. 3.108682. 17.396706. 12.263909. 0.001. B. 1170K. 2.903714. 32.671181. 12.263909. 0.001. - 13 -.

(24) 第四章結果與討論圖十為位能 A 在三個不同溫度下的粒子徑向分佈函數 g (r )，當系統的溫度越高時，系統的密度也隨著減小，分子間的動能會變強，其他粒子可以更靠近中心粒子，因此我們可以看到 g (r ) 的第一峰值的位置會往距離較小的一端靠近，且第一峰值的高度也會越來越低。溫度越高時，可以看到第一峰值後的震盪幅度也會變小，表示分子間的分佈更為均勻。圖十一為位能 B 的粒子徑向分佈函數 g (r )，位能 B 一樣具有當系統的溫度越高時， g (r ) 的第一峰值的位置會往距離較小的一端靠近，且第一峰值的高度也會越來越低，第一峰值後的震盪幅度也會變小，因為位能 B 的系統溫度差異較大，因此這些現象也更為明顯。如果把位能 A 在 1173K 和位能 B 在 1170K 時的模擬結果做比較時（圖十二），可以發現，位能 B 的第一峰值明顯較高，且位置也比較偏離中心粒子，顯示兩不同系統中粒子的分布狀況有所不同，在位能 B 中，第一層的粒子數目比較多，但過了第一峰值之後，外層粒子的分布是比較均勻的。接下來將兩位能的 g (r ) 經過傅利葉轉換後得到的靜態結構因子 S (k ) ，並將得到的 S (k ) 和 1170K 時的實驗值 [8] 作比對，如圖十三，圖中的位能 A 是粒子總數 N 為 3000 顆的結果，可以明顯的看到不管是那一個位能在那一個溫度，在 k 值較小的時候的吻合程度都不好，因為週期性邊界條件的關係，在 k 值較小時，此時的波向量波長接近系統長度，會造成波和波有重疊的現象，使得模擬的結果和實驗值會有明顯的差距，但如果我們將系統粒子數目增加的話 [8][9]，就可以盡量減少這樣的誤差（從圖六中可以明顯看出改善的情況）。隨著 k 到大約１ A - 1 左右，位能 B 已經和實驗值越來越一致了，且第一峰值和實驗值幾乎完全吻合，但是位能 A 卻明顯往 k 值大的一 - 14 -.

(25) 邊偏離且峰值大了很多，但後面的第二峰值和第三峰值兩位能就都能和實驗值有很好的吻合。另外從圖十四中也可以看到位能 B 在接近熔. G 點（約 600K）的 623K 之 S k ，除了 k 值較小的部分以外，整體的結. (). 果也都和實驗有很好的吻合。除了靜態結構因子 S (k ) 外，接下來我們探討兩位能動態結構因子 S (k , ω ) 的結果，我們所選擇的波向量 k 的範圍為 1.8 A - 1 -6.5 - 1 A - 1 之間，. 共取 6 個不同的 k 值，圖十五顯示兩位能在溫度接近 1173K 時的動態結構因子 S (k , ω ) 和實驗值在 1173K[10] 的比對，圖中的實心圓圈符號為實驗值，實線（位能 A）和虛線（位能 B）為我們電腦模擬的結果，但由於我們所探討的波向量 k 值較大，因此無法明顯看出 Brillouin 波峰和系統的聲速 (c S ) ，而且可以發現兩位能在 k 值較小時和實驗的結果都有相當程度的差異，雖然位能 B 在靜態結構因子 S (k ) 上比位能 A 來得接近實驗，但是在動態結構因子 S (k , ω ) 上，和實驗的結果卻也一樣有著滿大的誤差，波向量 k 為 1.8 A - 1 時兩位能的結果明顯不同，位能 B 的 Rayleigh 波峰明顯比位能 A 和實驗值高出很多，且在波向量 k 為 2.2 A - 1 時的誤差更為明顯，位能 B 的 Rayleigh 波峰比實驗值和位能 A 要低太多了，從第二章中的（ 7）式中看出我們可以從 S (k , ω ) 頻譜去計算得到 S (k ) 的值，而比對圖十三中兩位能系統在 k 為 1.8 A - 1 和. 2.2 A - 1 時的 S (k ) 值，可以發現確實在 1.8 A - 1 時的位能 B 之 S (k ) 值比位能 A 大（位能 A 約為 0.2337 位能 B 約為 0.5863），但是到了 2.2 A - 1 時卻相反，位能 A 之 S (k ) 值比位能 B 大（位能 A 約為 2.4315 位能 B 約為 2.0431），因此在 1.8 A - 1 和 2.2 A - 1 這兩個不同的 k 值時，兩位能的 S (k , ω ) 頻譜會有這樣的差異是合理的。但隨著 k 值越大，兩位能所模擬出的結果是差不多的，但在和實驗值的吻合度上都不是很好，顯示兩虛擬位能似乎都無法準確表現出液態鉛在此溫度的動力學性質。. - 15 -.

(26) 圖十六為位能 A 在三個不同溫度下的粒子速度相干函數 C (t )，當溫度越高時，分子間的動能越強，則系統的擴散性也會更好，故我們可以看到隨著溫度越高， C (t ) 的最低值會越大，隨著時間的演進， C (t ) 值會漸漸趨近於零，表示此時粒子間就不在有速度上的關聯性了。位能 B （圖十七）的 C (t ) 也有一樣的趨勢，隨著溫度越高， C (t ) 的最低值會越大，但是當我們把兩位能在接近 1170K 的結果做比較時（圖十八），會發現位能 B 的 C (t ) 的最低值明顯高於位能 A ，顯示在位能. B 系統中粒子較容易擴散。接下來我們應用 INM 分析，求解一系列的瞬間系統組態，經過統計後可以得到歸一化的狀態密度 D(ω ) ，圖十九為為位能 A 在三個不同溫度下的 D(ω ) ，可以看到隨著溫度越高， D(ω ) 在實數區段 (real lobe)的頻率分佈範圍越大，但峰值高度也隨之降低，但在虛數區段 (image lobe)則是隨著溫度增加，不僅頻率分佈範圍越大，峰值高度也同時增加，這和之前液態鎵 [2]的研究結果是一樣的，是大多數的液體都會具有的現象。接下來，從圖二十中來看位能 B 的在兩個不同溫度下之歸一化的 INM 狀態密度 D(ω ) ，隨著溫度越高， D(ω ) 在實數區段 (real lobe)的頻率分佈範圍越大，但峰值高度也隨之降低，虛數區段 (image lobe)隨著溫度增加，不僅頻率分佈範圍越大，峰值高度也同時增加，這一點和位能 A 的趨勢是一致的。圖二十一為位能 A 的參與比例 PN (ω ) 和頻率 ω 的關係圖，系統溫度越高， PN (ω ) 在頻率 ω ≈ 0 的值沒有太大的差異，另外配合圖二十二，可以看到. P2 N (ω ) （ N=500 ）隨著頻率 ω 變化的關係，隨著溫度增加， PN (ω ). 在實數的低頻區段（小於 5mev ）的參與比例幾乎為 1，表示在位能 A 中此低頻區段大都為延展模，這結果有點類似一些簡單液體，像是鋰、鈉 [3] 等，在頻率 ω 接近零的附近都是屬於延展模，但是像鈉的. - 16 -.

(27) 延展模分佈範圍明顯較廣。如果我們在實數的低頻區段中（小於. 5mev ），加入一個辨別條件，即化約後的參與比例 Sα ，當此區段的 INM 之 Sα ＜ 0.5 時，我們就將之稱為瞬間共振模 IRM (instantaneous resonant mode)，藉以和此區段的 INM 區分開來，而實際研究的結果顯示在位能系統 A 中的低頻實數區段找不到 IRM。接下來看到的圖二十三，是位能 B 的參與比例 PN (ω ) 和頻率 ω 的關係圖，和位能 A 的趨勢不同的是隨著溫度越高， PN (ω ) 在頻率 ω ≈ 0 的峰值明顯下降，配合圖二十四，. P2 N (ω ) （ N=500 ）和頻率 ω 的關係 PN (ω ). 圖中顯示，位能 B 隨著溫度的顯著增加，在低頻（小於 5mev ）區段的參與比例不全為 1，表示此低頻區段中並非完全都是延展模，且隨著溫度升高，可以看出有越來越高的比例並非延展模，為了在此實數低頻區分出 INM 和 IRM，我們一樣加入化約後的參與比例 Sα ＜ 0.5 的條件，從統計的結果顯示，我們可以明顯的發現 IRM 的存在，並在圖二十四中呈現 IRM 隨著頻率 ω 變化的關係圖，圖中的圓圈（ 623K ）與叉叉（ 1170K ）符號為原始的資料點，其中的實線為將數據平滑（ smoothed ）後的結果，可以看到隨著溫度的增加， IRM 的密度也明顯的增加。接下來我們進一步探究原因，從圖二十六之中，兩個不同的系統位能一個明顯的差別在於，對照位能和位能二次微分的關係圖來看，位能 B 在核心排斥力的區段裡面，有一個相當明顯的反曲點存在，而位能 A 的反曲點則是出現在後面的 Friedel part，在核心排斥力的區段中並沒有反曲點的存在，而從過去有關鎵金屬 [2] 的研究結果中顯示，我們可以進一步確認 IRM 的產生是由於系統位能在核心排斥力的部分有反曲點。. - 17 -.

(28) 第五章結論經過上述一系列的電腦模擬實驗結果，我們可以看到，兩種不一樣的液態鉛虛擬位能，位能 A 會隨著溫度的不同而有所差異，而位能 B 卻不受系統的溫度和密度不同的影響，而可以同時應用在兩個差異甚大的溫度上。當我們把兩位能的各項動力學性質之模擬實驗結果做一個比對分析時，可以發現兩位能系統的結構狀態明顯不同，但是在. G 溫度在 1173K 附近的靜態結構因子 S k 分析上，位能 B 明顯的較位能. (). G A 接近實驗的結果，而且位能 B 在接近熔點的 623K 之 S k 也和實驗. (). 結果有很好的吻合，表示位能 B 較能夠呈現出真實世界中液態鉛的結構狀態。且在位能 B 系統中，粒子明顯較容易擴散。但是在動態結構因子上，兩個位能系統卻都和實驗結果有相當程度的誤差，顯示兩者似乎都無法準確的表現出液態鉛在真實世界中的動力學性質。而在位能 B 中我們可以找到低頻的（小於 5mev）、準侷限（ quasilocalized） INM，稱為瞬間共振模（ IRM），但在位能 A 中卻找不到 IRM 的存在，因此我們可以進一步驗證如果液態鉛的系統位能在核心排斥力的部分存在有反曲點的話，我們就可以在實數的低頻區段發現瞬間共振模（ IRM），而且可以發現瞬間共振模（ IRM）的態密度會隨著溫度的升高而增加，以上兩點和之前對液態鎵的研究結果是一致的。. - 18 -.

(29) 參考文獻 [1] Z. A. Rycerz, and P. W. M. Jacobs, Molecular Physics 74, 425-441 (1991) [2] T. M. Wu, S. F. Tsay, S. L. Chang, and W. J. Ma, Phys. Rev. B 64, 064204 (2001) [3] T. M. Wu, and W. J. Ma, J. Chem. Phys. 105, 9281 (1996) [4] T. M. Wu, W. J. Ma, S. L. Chang, and S. F. Tsay, Physica B 316-317,606-609 (2002) [5] B. B. Laird, and H. R. Schober, Phys. Rev. E 66, 636 (1991) [6] Micheal Buchner, Branka M. Ladanyi, and Richard M. Stratt, J. Chem. Phys. 97, 8522 (1992) [7] Richard M. Stratt, Acc. Chem. Res. 28, 201 (1995) [8] Mikhail Dzugutov, and Karl-Erik Larsson, Phys. Rev. A 38, 3609 (1988) [9] R. Mentz-Stern, and C. Hoheisel, Phys. Rev. A 40, 4558 (1989) [10] O. Soderstrom, U. Dahlborg, and M. Daridovic, Phys. Rev. A 27, 470 (1983) [11] M. Canales, and J. A. Padro, Phys. Rev. E 56, 1759 (1997) [12] S. F. Tsay, Phys. Rev. B 48, 5945 (1993) [13] March Norman Henry and Tosi M. P., Atomic dynamics in liquids (1976) [14] March Norman H. Street Robert A. and Tosi Mario ed., Amorphous solids and the liquid state (1985) [15] March, Norman H. q(Norman Henry) and Tosi, M. P., Introduction to liquid state physics (2002). - 19 -.

(30) S(k). k (A-1) 圖一、位能 B 在溫度為 623K 的靜態結構因子 S(k) 和 k 的關係圖，其中＋號為 [8] 中位能 B 的模擬結果，方塊為實驗值，粒子數目為 16384. - 20 -.

(31) S(k). k (A-1) 圖二、位能 B 在溫度為 1170K 的靜態結構因子 S(k) 和 k 的關係圖，其中＋號為 [8] 中位能 B 的模擬結果，方塊符號為實驗值，三角符號是利用 [8] 中 (1) 式計算得到的結果，粒子數目為 16384. - 21 -.

(32) V(mev). r (A) 圖三、位能 A 在三個不同溫度下和距離 r 的關係圖. - 22 -.

(33) V(mev). r (A) 圖四、位能 A （實線）在 1173K 和位能 B （虛線）和距離 r 的關係圖. - 23 -.

(34) g(r). r (A) 圖五、位能 A 在溫度為 1173K 的的粒子徑向分佈函數 g(r)和距離 r 的關係圖. - 24 -.

(35) S(k). k (A-1) 圖六、位能 A 在溫度為 1173K 的靜態結構因子 S(k) 和 k 的關係圖，系統粒子數分別為 500 顆（實線）與 3000 顆（虛線）. - 25 -.

(36) C(t). t(ps) 圖七、位能 A 在溫度為 1173k 的速度相干函數 C(t) 和時間 t 的關係圖. - 26 -.

(37) ω G 圖八、在 IXS 或是 INS 實驗中所觀察到的動態結構因子 S k , ω 。其中. ( ). Cs 是系統聲速， D T 是熱擴散係數。. - 27 -.

(38) V(r). σ. ε. r (A) 圖九、鉛假位能函數化約單位取值示意圖. - 28 -.

(39) g(r). r (A) 圖十、位能 A 在三個不同溫度下的粒子徑向分佈函數 g(r)和距離 r 的關係. - 29 -.

(40) g(r). r (A) 圖十一、位能 B 在溫度分別為 623K（實線）和 1170K（虛線）的粒子徑向分佈函數 g(r)和距離 r 的關係圖. - 30 -.

(41) g(r). r (A) 圖十二、位能 A（實線）在 1173K 和位能 B（虛線）在 1170K 的粒子徑向分佈函數 g(r) 和距離 r 的關係圖. - 31 -.

(42) S(k). -1. k (A ) 圖十三、位能 A 和位能 B 分別在 1173K（實線，系統粒子數為 3000 顆）和 1170K（虛線，系統粒子數為 500 顆）的靜態結構因子 S(k)和實驗值 [8]的比較，其中實心圓點為實驗值. - 32 -.

(43) S(k). k (A-1) 圖十四、位能 B 在 623K（實線）的靜態結構因子 S(k)和實驗值 [8]的比較，其中實心圓點為實驗值. - 33 -.

(44) S (k , ω ) (mev-1). 1.8A-1. 2.2 A-1. ω (mev) 圖十五之一、位能 A（實線）在 1173K 和位能 B（虛線）在 1170K 的. G 動態結構因子 S k , ω 和實驗值 [10]（實心圓圈）的比較. ( ). - 34 -.

(45) S (k , ω ) (mev-1). 3.0A-1. 4.1 A-1. ω (mev) 圖十五之二、位能 A（實線）在 1173K 和位能 B（虛線）在 1170K 的. G 動態結構因子 S k , ω 和實驗值 [10]（實心圓圈）的比較. ( ). - 35 -.

(46) S (k , ω ) (mev-1). 5.6A-1. 6.5 A-1. ω (mev) 圖十五之三、位能 A（實線）在 1173K 和位能 B（虛線）在 1170K 的. G 動態結構因子 S k , ω 和實驗值 [10]（實心圓圈）的比較. ( ). - 36 -.

(47) C(t). t(ps) 圖十六、位能 A 在三個不同溫度下的速度相干函數 C(t)和時間 t 的關係圖. - 37 -.

(48) C(t). t(ps) 圖十七、位能 B 在溫度分別為 623K（實線）和 1170K（虛線）的速度相干函數 C(t)和時間 t 的關係圖. - 38 -.

(49) C(t). t(ps) 圖十八、位能 A（實線）在 1173K 和位能 B（虛線）在 1170K 的速度相干函數 C(t)和時間 t 的關係圖. - 39 -.

(50) D(ω). ω(mev) 圖十九、位能 A 在三個不同溫度下的 INM 狀態密度與頻率 ω 的關係圖. - 40 -.

(51) D(ω). ω(mev) 圖二十、位能 B 在溫度為 623K（實線）和 1170K（虛線）的 INM 狀態密度與頻率ω的關係圖. - 41 -.

(52) PN(ω). ω(mev) 圖二十一、位能 A 在三個不同溫度下的參與比例與頻率ω的關係圖. - 42 -.

(53) P2 N (ω ) PN (ω ). ω(mev) 圖二十二、位能 A 在三個不同溫度下的參與比例與頻率ω的關係圖. - 43 -.

(54) PN(ω). ω(mev) 圖二十三、位能 B 在溫度分別為 623K（實線）和 1170K（虛線）下的參與比例與頻率ω的關係圖. - 44 -.

(55) P2 N (ω ) PN (ω ). ω(mev) 圖二十四、位能 B 在溫度分別為 623K（實線）和 1170K（虛線）的參與比例與頻率ω的關係圖. - 45 -.

(56) D irm (ω). ω(mev) 圖二十五、位能 B 在溫度 623K(圓圈符號 )和 1170K(叉叉符號 )下的 IRM 狀態密度與頻率ω的關係圖，實線為將所有資料點平滑後的結果。. - 46 -.

(57) V(mev). r (A) 圖二十六、位能 A 在 1173K 和位能 B 在 1170K 的位能（實線）及位能二次微分（虛線）和距離 r 的關係圖. - 47 -.

(58)