水及鹽類水溶液之熱力學性質、鍵結行為及系統結構之探討:分子動力學模擬的研究

國

立

交

通

大

學

應用化學系

碩

士

論

文

水及鹽類水溶液之熱力學性質、鍵結行為及系統結

構之探討：分子動力學模擬的研究

Studies on the Thermal Dynamic Properties, Bonding Behavior and

Structure of Water and Ionic Aqueous Systems: A Molecular

Dynamics Simulation Study

研 究 生：周東和

指導教授：吳建興 教授

水及鹽類水溶液之熱力學性質、鍵結行為及系統之探討：分子動

力學模擬

Studies on the Thermal Dynamic Properties, Bonding Behavior

and Structure of Water and Ionic Aqueous Systems: A Molecular

Dynamics Simulation Study

研 究 生：周東和 Student：Dong-Han Chou

指導教授：吳建興 Advisor：Jiann-Shing Wu

國 立 交 通 大 學

應 用 化 學 系

碩 士 論 文

A Thesis

Submitted to Department of Applied Chemistry

College of Science

National Chiao Tung University

In partial Fulfillment of the Requirements

For the Degree of

Master

In

Applied Chemistry

June 2007

Hsinchu, Taiwan, Republic of China

中文摘要

本論文之目的，是找出一個適合模擬高分子溶液(Polymer solution)的 水分子模型(Water model)並且以此模型建立一水溶液(包含離子)模擬系 統。相較於其他水分子模型，Fexibile three center (F3C) model 已經用於大 尺度的巨分子(Macromolecules)溶液的模擬，並且相當符合實驗結果。由 於其勢能的短距離截斷(Short range truncation)的補償計算相當好，因此它 適合於需要高計算速度的模擬系統。在統計性質上，如水溶液的結構、動 態與熱力學性質在不同溫度、壓力以及系統尺度，都與實驗結果相當符 合。本研究在相同模擬系統下，統計性質也與文獻相當符合。

研究分為三部分：帶電的水溶液(Charged aqueous solutions)的結構與 動態性質；改變溫度、加入不同的離子種類與數量對奈米液滴的影響；給 予帶電水溶液的液滴不同的合併速度，對於液滴合併過程的影響。目前對 於帶電荷的奈米液滴的結構性質、碰撞以及合併過程的瞭解相當少。研究 結果顯示，水合作用(Hydration)對水溶液性質影響甚劇，主要是因為庫倫 靜電力影響水分子的空間排列與方位，直接影響氫鍵的形成。但是，由於 Cl...H 鍵結相當的強，因此能加強水分子在溶液中的結構性質。而模擬帶 電水溶液的奈米液滴(液滴半徑約 1nm)在真空系統中，是模擬液滴擁有不

同種類與數量的離子，如 Sodium ions、Chloride ions 或兩種離子等量存

在。純水溶液或是僅含有 Sodium ions 的液滴皆無法在 330-360K 的真空 環境下，維持奈米液滴的形狀。給予兩液滴額外的速度在特定方向，可觀 察兩液滴合併的過程。結果顯示，給予不同的速度，5 m/s、50 m/s、150 m/s 與200 m/s，兩液滴會從僅僅靠近乃至於完全合併。在 150 m/s 時，兩液 滴間會產生數個水分子所形成的架橋結構(Bridge structure)。而液滴的表 面積與體積是相當重要的參考依據，本研究藉由 Molecular modeling

積，得到了單一個液滴、兩液滴合併以及與合併液滴是相同狀態的液滴， 後二者的液滴的表面積對體積比率，確定此合併液滴是一穩定結構。

Abstract

This study investigates on obtaining a water model suited to long simulation time of macromolecules in solution and constructing a simulation system of aqueous solutions. Comparison between other models revealed that flexible three-center model has been already used in many large-scale simulation and it’s provided with experimental data. Because the model works well with short-range truncation suited to high-speed computation. It’s tested by comparing the structural, dynamic and thermodynamic properties of water in aqueous solution, at several temperature and density, with the other models and experimental data. Our program also was tested by calculation such properties and fitted these literature very well.

Therefore, the aims of research divide into three parts: First is tested for many water properties comparing with literature；Second investigates on the structural and dynamic properties of brine solutions；the other simulates charged aqueous nanodroplets for different condition, which were different temperatures, the number of ions, the type of ion and the size of droplets, in vacuum and these nanodroplets were given some additional velocity ranging from 1 m/s to 200 m/s to observe two nanodroplets bumped into or merged each other or merged . Studies show that ionic solvation shell effects strongly on the water structure in aqueous solution, like Cl- anion makes water more slow meaning ionic solvation shells are rigid. Computation of Bond time correlation functions shows that Cl--water pair can hold longer than water-water pair. The rigidity can play an important role in charged aqueous nanodroplets. At several conditions, the nanodroplet including Cl- ions were stable. Giving two nanodroplets a velocity in a direction to overcome the

surface energy of the droplet made a formation of bridge structure and giving more kinetic energy performed the merged process. To calculate the surface area and volume of a nanodroplet, there are the merged nanodroplet and a nanodroplet had the same number water molecules and ions of the merged nanodroplet, we use the molecular modeling software TINKER implemented the algorithms of solvent accessible surface. The result of computation can prove the merged nanodroplet is the stable structure.

誌謝

光陰冉冉似浮光，在即將畢業的前夕，回顧碩士生涯感觸良多。首先， 感謝指導教授，吳建興博士的指導與關懷，其學業上的諄諄教誨與事必躬 親的態度，使學生獲益良多；特別在學術研究上、其細心、仔細與專注的 程度，更令學生印象深刻，且深深體會學海無涯，為勤是岸的道理，在往 後的研究或工作生涯中，必定時時警惕自己。另外，更要感謝清大化工系 張榮語教授，在 CAE 研究室中，張教授提供最完善的資源與良好的研究 環境。每週撰寫研究與成果報告的訓練，更讓學生收穫不少，在口語表達 以及報告撰寫的能力大幅提升。另外，其學者風範與工程師的反應能力， 也深植學生腦海，期許自我能夠成為兼具工程師與科學家能力的人才。十 分有緣能夠與各位學長姐、同學以及學弟妹渡過這些共同研究、成長與玩 樂的時光，藉由謝誌一一表達對你們的感謝與祝福；煥錩、阮沈、文彥、 鎮杰以及啟夫學長祝你們的研究能順順利利；穎玫、明熙以及丘宏同學祝 你們在未來的人生旅途上能夠平平安安；建超、長福、煜智以及玉蕙祝你 們寫程式能夠順心如意，跑程式時，不當機、不發散以及結果符合物理意 義，順利地畢業。感謝我的父母親以及女友，幫助我跨越許多的挑戰與難 關。最後，由衷感謝各位！儘管光陰似箭，歲月如梭，但與大家回憶及情 誼深深烙印在心田中，難以忘懷！ 周東和 民國九十六年五月二十八日星期一

目錄

中文摘要...III Abstract ... V 誌謝... VII 目錄...VIII 圖目錄...IX 表目錄...XIII 符號說明...XIV 第1 章、緒論...1 第2 章、文獻回顧...5 2.1 分子動力學模擬...5 2.2 勢能模型...8 2.3 鹽類水溶液的性質... 11 2.4 水溶液的奈米液滴(Aqueous Nanodroplets) ...14 第3 章、研究方法...16 3.1 理論基礎與模擬流程圖...16 3.2 勢能模型...18 3.3 計算與分析的方法...21 第4 章、結果與討論...35 4.1 水分子之徑向分佈函數的模擬結果...35 4.2 模擬水溶液的超臨界狀態...41 4.3 水溶液的擴散係數...46

4.4 鹽類水溶液(Ionic Aqueous Solutions)的性質...48

4.5 帶電水溶液的奈米液滴(Charged Aqueous Nanodroplets)...62

4.5.1 液滴性質 ...62

4.5.2 液滴的運動與合併(Merged)過程...73

第5 章、結論與未來展望...87

圖目錄

Fig. 1-1 噴流系統示意圖。 ...3

Fig. 1-2 左邊圖為日本Sony公司 2007 年所發表的電子紙。...4

Fig. 2-1 最簡化的勢能形式。 ...7

Fig. 2-2 Engineering Village 2 資料庫搜尋以Molecular dynamics simulation 為關鍵字的文獻數量與發表年代之長條圖。...7 Fig. 2-3 氫鍵形成時的幾何限制。 ...12 Fig. 3-1 本研究的模擬步驟流程圖。 ...17 Fig. 3-2 水溶液系統初始化結果圖 (216 個水分子及 4 個氯離子)。 ...21 Fig. 3-3 二維週期性邊界作用示意圖。 ...23 Fig. 3-4 二維週期性邊界作用示意圖。 ...24 Fig. 3-5 二維週期性邊界作用示意圖。 ...24 Fig. 3-6 鄰近列表法示意圖。 ...27 Fig. 3-7 徑向分佈函數物理示意圖。 ...31 Fig. 3-8 方均根位移統計的方法。 ...32 Fig. 3-9 氫鍵形成的幾何限制條件。 ...34 Fig. 4-1 依據文獻[2]的系統設定。 ...36 Fig. 4-2 文獻[2]的O…O徑向分佈函數結果。 ...37 Fig. 4-3 本研究的O…O徑向分佈函數結果。 ...37 Fig. 4-4 文獻[2]的O…H徑向分佈函數結果。 ...38 Fig. 4-5 本研究的O…H徑向分佈函數結果。 ...38 Fig. 4-6 文獻[2]的O…H徑向分佈函數結果。 ...39 Fig. 4-7 本研究的O…H徑向分佈函數結果。 ...39 Fig. 4-8 依據Liew[63]對液氣界面之模擬系統設定。...42 Fig. 4-9 計算不同溫度下，氣液介面的密度分佈情形。...43 Fig. 4-10 Liew[63]以可撓性勢能(SPC-MTR)模擬氣液界面的密度分佈結 果。...43

Fig. 4-11 將各溫度液氣界面厚度(Thickness between liquid and vapor interface)繪製成X-Y圖。 ...44

Fig. 4-12 將統計之結果與文獻模擬的結果[63]以及實驗數據(NBS/NRC Steam Tables, Hemisphere, WA, 1984)進行比較。 ...45

Fig. 4-13 各溫度之方均根位移的收斂情形。...47

Fig. 4-14 各溫度下，計算MSD斜率求得的擴散係數。擴 ...47

Fig. 4-16 Levitt[2]的紀錄離子運動的結果。...49 Fig. 4-17 水分子的速度相關函數。 ...50 Fig. 4-18 Chloride ions的速度相關函數。 ...50 Fig. 4-19 氫鍵與Cl…H鍵結的連續鍵結時間相關函數(Continuous bonded time correlation)。 ...51 Fig. 4-20 氫鍵與Cl…H鍵結的間歇的鍵結時間相關函數(Intermittent

bonded time correlation)。 ...51 Fig. 4-21 全部氫鍵數目統計方法示意圖。...52 Fig. 4-22 全部數目氫鍵數目徑向分佈對Sodium ion-水分子的徑向分佈結 果。...53 Fig. 4-23 全部數目氫鍵數目徑向分佈對Sodium ion-水分子的徑向分佈結 果。...53 Fig. 4-24 系統溫度 298K至 638K時，時間以及分子數目平均之氫鍵數目 統計結果圖。...54 Fig. 4-25 全部數目氫鍵數目徑向分佈對應Chloride ion-水分子的徑向分佈 結果(258K)。 ...55 Fig. 4-26 全部數目氫鍵數目徑向分佈對應Chloride ion-水分子的徑向分佈 結果(298K)。 ...55 Fig. 4-27 離子濃度(Chloride ion)對於平均氫鍵數目的影響。 ...56 Fig. 4-28 Chloride ion的水合作用力幫助水分子排列成四面體空間結構。

...57 Fig. 4-29 模擬溫度從 273K到 638K的O...O徑向分佈函數。 ...59 Fig. 4-30 將Chloride ions水溶液系統之溫度對應O…O徑向分佈函數之第 一波峰(Rv1)值繪製成圖。 ...59 Fig. 4-31 將Chloride ions水溶液系統之不同溫度對應Cl…H徑向分佈函數 值繪製成圖。...60 Fig. 4-32 將Chloride ions水溶液系統之不同溫度對應Cl…O徑向分佈函數 圖。...60 Fig. 4-33 離子濃度變化(0.0091g/cm3至0.273g/cm3)對Cl…O徑向分佈的影 響。似乎影響不大(約 5%以內的變化)。...61 Fig. 4-34 離子濃度變化(0.0091g/cm3至0.273g/cm3)對Cl…H徑向分佈的影 響。似乎也影響不大(約 5%以內的變化)。...61 Fig. 4-35 加入適當數量的Chloride ions 可在高溫時，穩定維持液滴的形 狀。...65 Fig. 4-36 相同溫度與初始化位置的狀態下，沒有加入Chloride ions 的奈

米液滴在高溫時，無法形成液。...65 Fig. 4-37 模擬奈米液滴運動之系統溫度變化情形。...66 Fig. 4-38 奈米液滴平均氫鍵數目分佈統計(273K)。 ...66 Fig. 4-39 奈米液滴平均氫鍵數目分佈統計(330K)。 ...67 Fig. 4-40 奈米液滴平均氫鍵數目分佈統計(360K)。 ...67 Fig. 4-41 奈米液滴平均水分子密度分佈統計(273K)。 ...68 Fig. 4-42 奈米液滴平均水分子密度分佈統計(330K)。 ...68 Fig. 4-43 液滴水分子平均密度分佈(Square)與離子位置分佈機率(長條圖) 的共繪圖。...69 Fig. 4-44 液滴水分子平均密度分佈(Square)與離子位置分佈機率(長條圖) 的共繪圖。...69 Fig. 4-45 液滴水分子平均密度分佈(Square)與離子位置分佈機率(長條圖) 的共繪圖。...70 Fig. 4-46 離子位置分佈機率圖(273K)。 ...70 Fig. 4-47 離子位置分佈機率圖(273K，1320 個水分子)。 ...71 Fig. 4-48 離子位置分佈機率圖(360K，1320 個水分子)。 ...71 Fig. 4-49 水分子平均密度對應離子位置分佈機率圖(330K、1320 個水分 子、13 個Chloride ions)。...72 Fig. 4-50 奈米液滴的運動軌跡瞬照圖(330K)。 ...75 Fig. 4-51 奈米液滴的運動軌跡瞬照圖(330K)。 ...75 Fig. 4-52 奈米液滴的運動軌跡瞬照圖(360K)。 ...76 Fig. 4-53 奈米液滴的運動軌跡瞬照圖(390K)。 ...76 Fig. 4-54 奈米液滴的運動軌跡瞬照圖(330K)。 ...77

Fig. 4-55 4 個Chloride ions與 4 個Sodium ions加入 50 個水分子之中。...77

Fig. 4-56 8 個Chloride ions加入有 50 個水分子的奈米液滴(330K)。 ...78

Fig. 4-57 8 個Chloride ions加入有 50 個水分子的奈米液滴(360K)。 ...78

Fig. 4-58 給予兩液滴特定速度相撞(+5 m/s, -5 m/s)。 ...81 Fig. 4-59 給予兩液滴特定速度相撞(+5 m/s, -5 m/s)。液滴平均動能變化。 ...81 Fig. 4-60 給予兩液滴特定速度相撞(+50 m/s, -50 m/s)。 ...82 Fig. 4-61 給予兩液滴特定速度相撞(+50 m/s, -50 m/s)。液滴平均動能變 化。...82 Fig. 4-62 給予兩液滴特定速度相撞(+150 m/s, -150 m/s)。 ...83 Fig. 4-63 給予兩液滴特定速度相撞(+150 m/s, -150 m/s)。液滴平均動能變 化。...83

Fig. 4-64 給予兩液滴特定速度相撞(+150 m/s, -150 m/s)。形成架橋結構。 ...84 Fig. 4-65 給予兩液滴特定速度相撞(+200 m/s, -200 m/s)。液滴會合併，達 到質心-質心間距離的最小值(約 1.3nm)，形成合併的結果。 ...84 Fig. 4-66 給予兩液滴特定速度相撞(+200 m/s, -200 m/s)。 ...85 Fig. 4-67 建立 100 個水分子(含有 16 個Chloride Ions)的液滴。 ...86

表目錄

Table. 2-1 F3C model在各種性質模擬計算之結果[2]，在與實驗結果有良好 的準度，與其他勢能比較上(文獻值)，兩者結果大約相近，F3C甚至 更接近實驗值(其中gp1 表示徑向分佈函數中第一波峰、gv1 表示第 一波谷、rp1 表示第一波峰的距離、以及rv1 表示第一波谷的距離)。 ...10

Table. 2-2 Ewald sum和截斷勢能法之比較[37]。...10

Table. 3-1 原子質量參數...19 Table. 3-2 鍵結拉伸能量參數...20 Table. 3-3 鍵結彎曲能量參數...20 Table. 3-4 凡得瓦能量參數...20 Table. 3-5 各原子的帶電價數...20 Table. 4-1 溫度與密度改變對水分子結構的影響 ...35 Table. 4-2 模擬巨觀性質之系統設定，不同溫度對應之密度。 ...36 Table. 4-3 鹽類水溶液之徑向分佈函數的文獻結果 ...40 Table. 4-4 有關於鹽類水溶液之配位數的文獻結果 ...40 Table. 4-5 擬合曲線之氣液相密度 ...44 Table. 4-6 不同種類的勢能模擬結果 ...44 Table. 4-7 離子與離子間受到最強作用力的距離 ...49 Table. 4-8 液滴合併前後以及表面積對體積比率之關係 ...85

符號說明

α 水分子對水分子間之∠HOO sc A 截斷函式之補償係數 CN 配位數 v C 速度自相關函數

( )

HB

C t 間歇氫鍵相關函數(Intermittent hydration bond correlation

function)

D 擴散係數

d 氣液界面厚度之參數

t

Δ 步進時間(Timestep)

ε Lennard-Jones (Van der Waal) potential 能量參數

( )

g r , Gr 徑向分佈函數(Radial distribution function)

ξ _{Nose-Hoover thermostat 之摩擦參數} F 作用力場(Force Field) f 真空界電係數 (Vacuum permittivity) g 系統自由度

( )

h t 紀錄時間t 鍵結是否存在的函數

( )

H t 紀錄起始時間至時間t，鍵結是否連續存在的函數 k 波茲曼常數 K_αβ 動能項(計算壓力用)

L 週期性系統(Periodic boundary condition)的邊長

N 系統粒子總數 K N K 時間時的徑向分佈結果 i m _粒子i(原子或分子)的質量 i p 粒子i(原子或分子)的動量向量 Q Nose-Hoover thermostat 之控溫參數 i r _粒子i(原子或分子)的位置向量 i r&, v_i 粒子i(原子或分子)的速度向量 i r &&, ai 粒子i(原子或分子)的加速度向量 ij r 粒子i 與粒子 j 之距離向量 rp1 徑向分佈函數第一波峰值 rv1 徑向分佈函數第一波谷值

S 截斷位移函式(Truncation shift function)

( )

HB

function)

σ Lennard-Jones (Van der Waal) potential 平衡距離參數

D

T 控溫欲達到之溫度

U 勢能場(Potential Field)

W_αβ 維里項

z 吉布氏分離表面之參數(Gibb’s divding surface)

第1章 、緒論

藉由分子動力學 (Molecular dynamics, MD) 模擬研究各種材料的 結構、流變、動態及熱力學等性質的文獻正以爆炸性的速度發表。在 那些連續理論方程式不能良好描述或近似實際結果的部分、實驗的方 法相當困難或不易觀察的部分、以及一些特殊的微特徵的研究，電腦 模擬能提供一個強大的研究工具，深入地瞭解這些現象與性質。分子 動力學將分子視為數個粒子 (All atom model) ，以牛頓與漢米頓力學 的運動方程式，選擇模擬材料的勢能函式 (Potential function) ，搭配適 當的數值方法 (設定積分的步進時間間隔) ，解得粒子的運動軌跡、速 度與受力等。其能夠仔細地顯示運動軌跡，提供上述情形在微觀 (原子 階層) 的觀察與統計熱力學的分析。 本研究室致力於奈米加工模擬領域，開發奈米製程的模擬技術。 目前相當看好奈米液滴應用範圍與學術研究的價值。如以奈米噴流加 工製程(Nanojet printing process) 或奈米噴墨列印製程 (Nano ink-jet printing process) 技術產生奈米液滴，並且用於光電元件的製程，如 Fig. 1-1所示。遂將其研究背景分為以下幾點： 1. 奈米液滴是藉由奈米噴流加工製程所產生 (其製程可大量節 省成本，並且適合溶液製程，特別適合光電高分子材料的應 用) 。 2. 產業界需要相關奈米製程技術之研究與開發，爭取奈米軟性電 子產品的商機。如日本Sony公司以領先全球的製程技術，首先 發表電子紙 (參考 Fig. 1-2) 的高科技產品。此產品廣受矚目， 公司股價隨之上升。 3. 奈米液滴能廣泛應用，例如製造奈米光子材料 (如 Quantum dot) 、奈米金屬導線、奈米焊點與奈米矽基材 (Substrate) 等 等方面的奈米光電元件。另外，其他方面的應用，還包含質譜 儀分析技術 (電灑游離質譜法會產生許多帶電荷的液滴) ，噴

霧塗佈 (Spray coating) 、噴霧冷卻 (Spray cooling) 以及噴霧 成形 (Spray forming) 等等方面。 另外一方面，過去研究所使用的水分子模型 (Water model) [1]在本 研究中重新考慮計算效率與結果準確性的平衡；過去研究中，計算庫 倫靜電力的分子有一困難的挑戰，便是其數學模型在庫倫靜電力的處 理(非收斂性)；過去的演算法是相當複雜且計算量龐大，導致耗費時間 成本高。然而，Flexible three center point model (F3C model) 計算此類 遠距離作用力，具有高效率 (截斷勢能) 與高準確度 (補償方程式) 。 根據 Levitt 等人[2](1997)指出，F3C 是被建立成能夠有效率地模擬巨 分子 (Macromolecules) 在溶液中的運動軌跡 (Trajectories) 的勢能模 型，此最大的不同之處在於能量的收斂性與非鍵結截斷技巧，能夠讓 計算時間大幅減少。目前已經用於模擬許多長時間的巨分子溶液的運 動軌跡，並且在數種溫度與密度下，模擬結果與實驗結果相當符合[2]。 本研究將以 F3C 模擬純水溶液、鹽類水溶液以及帶電水溶液的奈

米液滴 (Charged aqueous nanodroplets) 。在純水溶液與鹽類水溶液的

系統中，控制的變因有溫度、密度 (離子數目) 、加入離子種類 (Na+、

Cl-和 NaCl 離子對) 。統計的性質包括：壓力、擴散係數 (Diffusion

coefficient) 、超臨界溫狀態的模擬、徑向分佈函數 (Radial distribution function) 、平均氫鍵數目、氫鍵時間相關函數 (Hydrogen bond time correlation function) 以及速度自相關函數 (Velocity-Velocity correlation function) 等等。

研究的最後一部份，是模擬帶電水溶液的奈米液滴 (Charged

aqueous nanodroplets) 的運動行為與合併 (Merged) 過程。水溶液是最

廣泛使用的溶液之一。從1970 年開始，微米尺度的水溶液的液滴之性

質或運動過程，已經透過理論分析[3-6]，實驗觀察[7-10]與電腦模擬[11, 12]等方面進行研究。近年來在奈米尺度的液滴 (純水、帶有電荷或高 分子溶液的液滴) ，以理論分析[13, 14]與分子動力學模擬[15-18]的研 究正在發展中。目前在這方面，有許多物理化學性質，運動過程與機

制尚待更深入的瞭解與釐清。模擬帶電水溶液的奈米液滴 (Charged aqueous nanodroplets) 對於研究室自行開發的程式更有兩點貢獻： 1. 建立適合計算高分子溶液的水分子模型之分子動力學模擬的 程式，其具高速的計算能力與性質的準確性 (請參考 Table. 2-1 與Table. 2-2) 。 2. 建立多種統計性質之程式架構：如徑向分佈函數 (Radial distribution function) 、 時 間 相 關 函 數 (Time correlation function) 、氫鍵鍵結的統計、液滴表面積與體積的計算等等 方面 (請參考第三章) 。

帶電水溶液的奈米液滴 (Charged aqueous nanodroplets) 的研究可 分為： z 改變液滴的溫度、尺寸 (Size) 、離子數量以及帶電荷類型分 析液滴結構的改變。 z 觀察液滴-液滴間的吸引或排斥的運動過程，並且提供特定方 向一額外的動能模擬液滴的合併過程。 Fig. 1-1 噴流系統示意圖。藉由微奈米尺度的噴嘴系統製造出奈米液滴，是目前受 到重視的高潛力加工製程系統之一。

Fig. 1-2 左邊圖為日本 Sony 公司 2007 年所發表的電子紙。中間圖為該產品內部結 構示意圖，是以該公司所自行開發的微加工製程技術所製造。 資料來源：http://www.sony.co.jp/sonyinfo/News/Press/200705/07-053/index.html 本研究並非以套裝軟體 (如 Materials Studio 4.2) 或開放源程式碼 (如 DL_POLY 214 program) 進行研究，而是先瞭解分子動力學模擬的 邏輯架構與物理意義，將數學方程式轉譯為C Double Plus 程式碼，並 且設計系統條件的設定方法與建立不同的統計性質所需要的演算法， 處理隨模擬時間輸出資料，如分子的位置、速度與作用力等數值，最 後統計得到所有的數據結果。

第2章 、文獻回顧

2.1 分子動力學模擬

分子動力學理論是 Irving[19]等人(1950)提出，以勢能函數模型 (Potential model) 表示粒子 (包含分子基團或原子) 系統的能量關係， 以數值演算法解出牛頓運動方程式，可決定步進時間的各個粒子軌 跡、速度及加速度等結果，而結果的變化性質可用於計算系統的性質。 Alder [20, 21] 等人以硬球 (Hard spheres) 勢能模擬粒子的運動，硬球 的物理模型是當兩粒子接觸瞬間才有作用力產生，其後發展出軟球 (Soft Spherers) 勢能，又稱為 Lennard Jones potential，而 Rahman[22] 等人 (1964) 用其模擬液態氬的運動與計算氣體的性質。分子動力學模 擬的對象逐漸從原子，理想氣體、簡單液體、牛頓流體、非牛頓流體、 金屬、半導體、以及高分子等等。尤其近數十年來，拜電子計算機的 運算能力大幅提昇，分子模擬逐漸受到重視。以分子動力學理論的數 學模式 (包含分子勢能、數值積分演算法與統計熱力學性質) ，可就實 驗上觀察到的微觀系統的現象與性質，進一步驗證或說明，深入瞭解 實驗結果的意義；對於實驗數據提供內插或外插的結果，尤其是對於 一些特殊的量測區域或難度高的實驗上，更是強力的工具[23]。 就文獻發表數量來看，以分子動力學模擬為關鍵字的文獻成等比 級數上升(1998 到 2007)，如 Fig. 2-2所示，意味分子動力學模擬廣泛地 被研究與受到重視。特別在原子尺度 (Atomic Level) 的現象的觀察 上 ， 提 供 了 許 多 前 瞻 性 與 開 創 性 的 研 究 與 觀 察 結 果 。 例 如 ， 從 Microsecond到Nanosecond的蛋白質的詳細的摺疊路徑[24]、水分子結晶 動態過程[25]、水通過奈米碳管的運動與性質[26]與奈米噴流成形、斷 裂與理論分析[14]等文獻皆發表在Nature或Science等期刊中。 然而，粒子的物理模型一複雜，其數學式也就越難處理。在分子 動力學中，基本勢能的形式如 Fig. 2-1所示，其他勢能如Amber[27]與

Compass[28] model等，數學模型極為複雜且計算繁瑣，在轉譯程式碼 的工作上，面臨極大的困難與挑戰。因此，本研究首先需尋找合適的 勢能，在計算量與結果準確程度取得平衡，將其數學函式轉譯為程式 語言。自行撰寫程式有幾個優點： 1. 程式架構能夠根據研究目的而設計，依據模擬結果進行設計與 改良。 2. 參數的設定自由度高，不受限於軟體開發者所預設的勢能參 數，可依據最新發表的文獻更改參數與設定。 3. 方便加入各種統計性質的程式碼，調整統計細節設定與控制輸 出結果。 本文將就勢能函式、離子水溶液的性質、以及奈米液滴等文獻進 行回顧。

Fig. 2-1最簡化的勢能形式。包含拉伸、彎曲、扭轉與非鍵結形式的凡得瓦和庫倫 勢能(Potential function)[8]

Fig. 2-2 Engineering Village 2 資料庫搜尋以 Molecular dynamics simulation 為關鍵字 的文獻數量與發表年代之長條圖。

2.2 勢能模型

水是自然界中最重要的溶劑，由於其多種特殊 (Anomalies) 的物 理性質，如加入特定的離子會有表面張力-濃度的最小化傾向或是阻礙 液滴 (Bubbles) 的結合 (Coalescence) 。另外，高黏度 (水在 25℃黏度 為0.89cp，而 Pentane 為 0.22cp) 、高表面張力 (水在 20℃時表面張力 為72.75mJ/m2，而CCl4 為 26.6mJ/m2) 與受到氫鍵鍵結影響甚劇等等， 因此水分子的物理性質與結構廣泛地受到分析與研究。Kusalik[29]等人 (1994) 首先以分子動力模擬研究水分子的空間秩序 (Spatial Order) ， 過去液態結構通常會以原子間的徑向分佈函數 (Radial distribution functions) 特徵化，實驗上藉由 X-ray 繞射 (X-ray diffraction) 與中子散 射 (Neutron scattering) 或透過電腦模擬也可以得到，然而這些還不能 十分清楚的證實液態分子的空間秩序。然而，由分子動力學模擬，透 過計算勢能得到運動軌跡，其模擬結果可以仔細地由三維圖形來顯示 局部原子密度與結構的情形，這類結果可以幫助瞭解更多有關於水的 性質與解釋微觀現象。 在分子動力學中，可用於模擬計算水分子之相關性質的模擬有許 多種類，已知約有數十種以上，但其中常見的有：

1. Three Point model : SPC/E model[30]或F3C model[2, 31-34] 2. Multicenter model : TIP4P model[35]

3. Polarizable model : PLO1 model[36]

每一種勢能的計算情形略有不同，如 SPC/E 與 TIP4P model 一般

以固定鍵長與鍵角的演算法計算運動方程式，而 F3C model 鍵長及鍵

角可撓，具有較大的彈性。勢能本身的形式對於計算速度影響甚劇。

本文選擇F3C model[2, 31-34]，其理由如下：

凡得瓦力與庫倫力)的計算方面，經過 Levitt 等人[2, 31-34]的 研究下，建立一組截斷勢能與補償的 (Compensated) 方程式， 後者隨者截斷半徑的改變，調整一比例係數 (Asc) ，目的為減 少凡得瓦與庫倫靜電排斥能量 (Repulsive energy) ，補償因截 斷半徑造成的吸引力損失。而截斷勢能可大幅加快計算速度， 節省計算時間。

2. 在 結 構 性 質 (Structural properties) ， 動 態 性 質 (Dynamics properties) 與熱力學性質 (Thermodynamics properties) 都能 準確地模擬計算而得，換句話說，一般的勢能可能只對某些特 定性質的模擬計算，在其他方面的性質可能就與實驗結果相差 甚遠。然而，Levitt 等人[2, 31-34]所發展的 ENCAD Package

中，所使用的F3C model，卻能夠具有準確的模擬各種性質，

如Table. 2-1 所示，因為此研究團隊先前建立了龐大的分子模

擬與實驗數據資料庫的關係，藉由完整的理論 (量子計算與分 子動力學) 與大量的實驗數據建立了 F3C model。

本研究選用Levitt[2, 31-34]等人所提出Flexible three center point model (簡稱F3C) ，主要因為該模型的數學計算在處理長距離作用力 (作用力為距離的 1 到 3 次方的函數，) 已經截斷處理，且將誤差值考 慮於其補償函數中，而計算216 顆水分子時，時間可節省一倍以上[37] (計算量與粒子數的平方成正比，意味數量越多，計算時間差距將大於 兩倍) 。但是，其勢能的缺點便是在能量計算的準確性略低於Ewald Sum (為一種計算非收斂性作用力的方法，請參考[37-39]) ，如 Table. 2-1與 Table. 2-2所示。文獻[37]也指出，因為兩種勢能所得到運動方程 式之解幾乎相同 (運動軌跡) ，所以即使改變溫度與密度，兩者計算出 來的結構性質都近乎相同。 就一般而言，F3C model 是相當優秀的勢能，目前知道的研究議題

有：Protenis 與 Nucleic acids 在溶液的分子動力學之類的溶液系統模擬

高分子溶液[45]等等，其應用範圍相當廣泛。更重要的是，F3C model 是一種很適合用於模擬分子數量高 (數萬顆) ，系統尺度 (數百奈米) 大的奈米加工系統之勢能。 Table. 2-1 F3C model 在各種性質模擬計算之結果[2]。在與實驗結果有良好的準 度，與其他勢能比較上(文獻值)，兩者結果大約相近，F3C 甚至更接近實驗值(其 中gp1 表示徑向分佈函數中第一波峰、gv1 表示第一波谷、rp1 表示第一波峰的距 離、以及rv1 表示第一波谷的距離)。

Thermodynamic Kinetic Structural

GO…O Heights Positions(Å)

model T(k) U(kcal _mol-1₎ Cv(calmol -1 k-1) D×104(cm2s-1) _{Gp1 Gv1 Gp2 Rp1 Rv1 Rp2} F3C[2] 298 -9.6 26 0.24 3.19 0.83 1.07 2.8 3.3 4.4 TIP4P[46] 298 -10.1 19 0.32 3.0 0.81 1.12 2.8 3.4 4.5 SPCE[30] 300 -10.8 0.25 3.0 0.8 1.10 2.7 3.3 4.7 實驗值 298 -9.9 18 0.23 2.91 0.73 1.14 2.9 3.3 4.5

Table. 2-2 Ewald sum 和截斷勢能法之比較[37]。

方法 總能(kj/mol) 位能(kj/mol) 動能(kj/mol) 計算時間(sec/step) Ewald Sum -34.23598 -41.82354 7.3937 2.44 截斷 勢能法 -32.79025 -40.23228 7.25216 1.21

2.3 鹽類水溶液的性質

水 是 最 常 見 的 溶 劑 之 一 ， 結 構 看 似 簡 單 卻 具 有 許 多 奇 特 (Anomalies) 的性質，尚未有深入與清楚的瞭解。本研究藉由分子動力 學模擬，系統地研究水分子與鹽類水溶液的結構、動態與熱力學性質。 其中，鹽類水溶液的許多性質，例如黏度、擴散性與配位數等等由實 驗與理論的方式證實為濃度的函數，這是因為離子影響水分子結構。 一般認為水分子在短距離 (Short range) 具有不規則的 (Random) 網狀 氫鍵結構，Hao[47]等人 (2007) 而帶電性質 (離子種類的半徑也會影 響，但本研究不做探討) 與濃度改變時，離子會以不同方式改變水分子 的網狀結構。 影響水分子性質最重要因素是改變了氫鍵結構；一般而言，水分 子的氫鍵半徑 (H…O距離1.80Å) 遠比氧原子與氫原子的凡德瓦半徑 合 (1.40 + 1.20 = 2.60Å) 小，但又比一般之氫氧單鍵 (0.97Å) 來得 大。氫鍵的鍵能約在15到20 kJ/mol 之間，而一般的共價鍵的鍵能大約 是400 kJ/mol ， 一 般 室 溫 (25℃) 的 熱 能 大 約 等 於 2.5 kJ/mol 。 Chandrasekhar[46]等人 (1984)模擬Alkai mettal ion的稀薄水溶液，對於 溶質-溶劑、溶劑-溶劑，以及第一與第二水合層的能量關係說明液態水 溶液的結構資訊，離子作用力擾動第一水合層內的氫鍵結構，尤其以 Li+與Na+影響最大。Sivaraja[48]等人 (1992) 以分子動力學模擬液態水 分子的氫鍵缺陷對傳輸 (Transport) 與流動行為的影響。Corongi[49, 50] 等人 (1993) 模擬6種溫度下的水合配位數(Hydration number)隨離子與 水分子間作用力距離的變化。Koishi[51]等人 (2005) 以分子動力學模 擬 奈 米 尺 度 的 雙 層 疏 水 性 薄 膜 之 間 ， 計 算 水 分 子 的 氫 鍵 半 衰 期 (Lifetime) ，而氫鍵結構受到溶液環境影響甚劇。而各種不同系統下， 氫鍵的定義也略有不同，決定於其作用力的第一水合層半徑，如 Koishi[51]等人 (2005) 定義鍵結態生成時，O…H距離小於2.2Å且作用 力能量 (Interaction Energy) 低於-10kj/mol，計算結果約為0.2-0.5ps (一

般情況約為1ps) ；另外，如Caleman[52]等人 (2006) 則以O…O距離小 於3.5Å，以及氫鍵提供者 (Donor) 與接受者 (Acceptor) 夾角小於30° 的幾何條件定義氫鍵的產生；或者如Hanasaki等人 (2006) [53]添加 O…H距離小於2.4Å的限制。本文所使用的定義是參考Ju[15]等人 (2005) 所使用的幾何條件，Ju等人以F3C勢能對純水溶液的奈米液滴的氫鍵結 構進行研究。該幾何限制條件為： 1. O…O的距離小於3.4Å。 2. O…H的距離小於2.4Å。 3. α (α =∠HOO) 小於30° Fig. 2-3 氫鍵形成時的幾何限制。此定義依據 F3C model 的作用力情形所定義的條 件。

通常會依據勢能的徑向分佈函數 (Radial distribution function) 結果設 定氫鍵的截斷距離 (Cut-Off distance) 。另外，絕大部分的文獻都會採 用徑向分佈函數及配位數的統計結果 (其數學式參見第三章) ，分析水 分子的局部結構，方便瞭解水溶液在空間秩序分佈的情形，其細節將 在第三章詳述。 關於離子水和 (Ion hydration) 方面，離子與水分子有長距離作用 力，會形成水合層。Nguygen[54]等人 (1984) 率先以分子動力學模擬 鹼金族離子 (Alkali cations) ，如Li+、Na+與K+以及鹵素離子 (Halide

anions) ，如F-與Cl-的運動。Straatsma[55]等人 (1988) 藉由分子動力學 模擬分析離子水溶液的自由能差距，其研究提供計算自由能的熱力學 積分技巧(Thermodynamic integration technique ) 。Rasaiah[56]等人 (1996) 以SPC/E勢能模擬常溫下，離子 (Alkali和Halide ion) 的移動性 (Mobility) ，證明不同大小與種類的離子，其移動性是掉落在分離的曲 線(Separated curve) 。Obst[57]等人 (1996) 以CHARMM22作用力場， 以分子動力學模擬水合的Alkaline以及Alkaline-Earth 金屬陽離子，研 究顯示，在徑向分佈函數與配位數的統計結果方面與實驗值相當吻 合，計算擴散係數皆略高於實驗結果，氫鍵的殘存時間 (Residence time)

大略隨離子半徑增加而減少，但是Mg2+卻是例外。White[58]等人 (2000)

以First principles molecular dynamics simulations (or Ab initio) 模擬

Na+的溶解情形，第一水合層 (First solvation shell) 與實驗量測結果符

合 (平均有5.2個水分子) ，並且試圖找到任何Na+影響超出水合層之水

分子方位 (Orientation) 的結果，但是沒有找到。在這之前，不曾以模 擬方式預測水合層的結構性質。Bakker[59]等人 (2001) 以實驗的方 式，觀察水分子在水溶液水合層(Aqueous solvation shell) 的化學環境的 動力學，需先加入電解質 (KF、NaCl或NaI等等) 於溶液中，再進一步 以非線性光譜技術 (Nonlinear spectroscopic technology) 觀察選擇的水

分子群組，受到Cl-、Br-或I-對水分子的氫鍵的動力學性質之影響，並 且發表在2001年Science上。Guardia [60]等人 (2005) 分析水溶液的離 子溶液的氫鍵特徵，模擬Alkali metal以及halide 的水溶液系統，藉由 在第一水合層與第二水合層的統計結果，發現離子種類會影響氫鍵的 半衰期 (Lifetime) 。在第一離子的水合層中存在水分子的四面體結 構，且間歇性 (Intermittent) 的氫鍵半衰期會因離子而增加。Rasaiah[47] 等人 (2007) 以SPC/E模型模擬Alkali halide 溶液在結構與動態性質的 不同之處。研究顯示，水分子會有親水性作用力 (Li+) 或類似疏水性 作用力 (Rd+或Cs+) 的行為，種類不同的離子對 (如RdCl與LiCl比較) 及濃度會以相當不同的方式影響系統的平衡狀態與結構性質。

2.4 水溶液的奈米液滴(Aqueous Nanodroplets)

Ju[15]等人 (2006) 研究中指出，目前對於水分子的結構性質還是 瞭解不足，特別在溫度與尺吋 (Size) 的影響。由於微觀尺度並沒有建 立良好的理論幫助我們瞭解這些問題，因此分子動力學模擬能提供有 效的近似結果。而Ju[15]等人 (2006) 分析溫度與尺吋 (Size) 變因對於 奈米純水液滴的影響，研究中指出，液滴可分為兩個區域，其中內部 區域的性質就像是塊材 (Bulk) 一般，沒有特定方位以及尺吋影響。在 外部區域中，液滴的尺寸大大影響水分子分佈情形，在較大尺寸的液 滴表面，因為具有較大的表面張力，氫鍵鍵結數量也較大；反之，較 小尺寸的液滴表面因表面張力較小，氫鍵數量也較少。隨著溫度上升， 兩個區域內氫鍵數目皆下降，但是表面區域的分子數有上升的情形。 電噴霧技術 (產生帶電液滴) 是具有重要實際應用價值的技術，如 用於高分子溶液的電泳分離(Electrophoresis in a concentrated dispersion of polymeric drops) 、蛋白質與核酸的質譜分析 (MS) 、產生帶電墨滴 再加以控制、以及乳化劑等等，都是重要且相當實際的應用。如電噴 霧可使溶液中的分子帶電而離子化，深入電噴霧的巨觀過程，有兩種 機制：離子揮發 (Ion evaporation ) 以及帶電殘存機制 (Charge residue mechanism) ，但在奈米尺度下，學者尚未有定論。而斷裂是電噴霧離 子 化 (Electrospray ionization) 以 及 Electrospray mass spectrometry (ESMS) 的基礎，更可應用於帶電奈米噴流行為。對於奈米液滴而言， Ion evaporation、Coulomb explosion 以及 Charge residue model 等 model 都尚未能解釋奈米尺度的行為[17] 。利用分子動力學模擬研究水溶液 奈 米 液 滴 (Aqueous nanodroplets) 的 斷 裂 機 制 ， 就 我 們 所 知 到 ， Marginean[17]等人 (2003) 的研究，是最早以分子動力學模擬帶電液滴 的斷裂機制，並且統計液滴的應力張量分佈。

在斷裂前的液滴尺寸與帶電價數能夠被紀錄下來，關於這樣臨界斷裂

狀態的研究可以參考Consta[18]等人 (2006) 。他們也指出，液滴能夠

擁 有 (Hold) 多 少 的 價 數 而 不 斷 裂 以 及 庫 倫 爆 炸 極 限 (Coulomb explosion limit) 是目前他們所關心的議題。其研究藉由 Rayleigh＇s charged droplet theory，在相同條件下，除了統計液滴半徑與表面積外， 還統計不同分子數目與溫度下，液滴斷裂次數與半徑變化。研究顯示， 分子動力學模擬的結果大於理論方程式的預測結果，這可能是因為 Rayleigh＇s theory 沒有考慮到離子分佈對於液滴熱力學性質的影 響；與Born＇s theory 不同之處在於，含有一或二個離子的液滴半徑 與液滴初始分子數無關 (也就是液滴半徑是沒有規則) 。其研究也指 出，過去相關模擬文獻，離子分佈位置傾向液滴表面附近 (單一顆離 子) ，在其研究結果中 (多顆離子) ，卻是傾向分佈在液滴的中心部 分。其液滴斷裂機制是以形成類似時架橋 (Bridge) 結構，慢慢形成角 錐形狀在斷裂開來。

第3章 、研究方法

3.1 理論基礎與模擬流程圖

本研究模擬以F3C模型模擬水溶液與奈米液滴 (分別有純水系 統、含有Na+、Cl-或NaCl離子對共存的水溶液系統) 的運動過程，計算 結構、動態與熱力學性質，以溫度、密度以及不同電性的離子對於系 統性質的影響。 分子動力學是以牛頓力學與漢米頓力學計算作用力，預測分子或 粒子在相空間的軌跡 (Phase-Space Trajectories) 變化，其有幾點假設： 1. 所有粒子 (包含分子基團或原子) 皆遵守古典牛頓力學的運動 定律：將粒子間交互的作用力單純的簡化在勢能方程式上，忽 略其他效應的作用現象。 2. 粒子的交互作用滿足疊加原理：意即將所有對粒子作用的各種 效應累加在一起，包括量子效應和多體作用力 (Many body force) ，將其整合於一勢能函式上，因此，當以量子力學或其 他散射實驗測量方法所求得的粒子的勢能函式，即總括了所有 的粒子受力效應，之後再利用古典力學來建立粒子的運動方程 式。 3. 假設每個粒子均為球型，並以質心作計算。 了解上述的基本假設後，首先要先建立研究的數學模型，藉由文 獻所提出的勢能模型 (Potential model) ，以漢米頓力學的方式得到作 用力函式 (Force filed) ，接著建立合適的數值演算法、控溫演算法及 統計方法，將一一說明之。程式的計算流程如Fig. 3-1所示：

Fig. 3-1 本研究的模擬步驟流程圖。在溶液分子的製備，則是達到設定溫度、密度 或壓力 (呈穩定震盪收歛值) 。

3.2 勢能模型

本研究使用Levitt[2, 31-34]等人所發展的Flexible Three-Centered water model (簡稱為F3C model)。各原子質量參數與編號請參考

Table. 3-1，以下分別說明各種作用力的數學形式： 1. 拉伸作用力 (Stretching interaction) 是以式 (3-1) 代表拉伸勢 能(Stretching potential) ，假想一彈簧可描述鍵結間的簡協運 動，用虎克定律來描述其運動行為。式 (3-1) 中Ks為伸縮常 數，d_eq為平衡鍵長。由式 (3-1)微分得到拉伸作用力，如式(3-2) 所示。能量參數參考 2. 3. Table. 3-2。

(

)

2 stretch b eq U =K r−r (3-1) i =m i = −∂U ∂ i F r&& r (3-2) 4. 彎曲勢能 (Bending interaction) 是以式 (3-3) 為鍵角的彎曲勢 能，亦相似虎克定律，不同的是描述鍵結間鍵角的簡協運動， 其中K_θ 為鍵角彎曲常數，θ為三個鄰近原子間的鍵角，θ_eq為平 衡鍵角。經由式 (3-3) 微分得到粒子的作用力函式，i、i+1與i+2 粒子的彎曲作用力分別為式 (3-4) 、式 (3-5) 以及(3-7) 。能 量參數參考Table. 3-3。

(

)

2 Bending eq U =K_θ θ θ− (3-3)

( )

, 1, 2 2 ( - ) i Bending i i i eq F = k_θ θ _{+ +} θ u v (3-4)

(

)(

)

(

)

-3 2 2 , 2, 2 , 2 2, 2 , 2, 2 , 2 , 2, 2 1-i i i i i i i i i i i i i i i i i i i i u d c c c d c c c v c c c + + + + + + + + + + + = − =

( )

2 _{2 ( - )} i Bending i eq F+ = k_θ θ θ w t (3-5)

(

)(

)

(

)

-3 2 , 2, 2 , 2 2 , , 2, 2 , 2 , 2, 2 1-i i i i i i i i i i i i i i i i i i i i w d c c c d c c c t c c c + + + + + + + + + = − = 1 2 0 i i i

Bending Bending Bending

F +F+ +F+ =

Q (3-6)

(

2

)

1 i i i

Bending Bending Bending

F + F F +

∴ = − + (3-7)

其中 di = −ri ri+₁,di+₂ =ri+₂−ri+₁,cii = ⋅d d ci i, i i_{, 1}+ = ⋅d di i+₁,ci+_{2, 2}i+ =di+₂⋅di+₂

5. 非鍵結作用力 (凡得瓦與庫倫作用力) 是以非鍵結力，包含式 (3-8) 凡得瓦力 (Van der Waal force) 以及式 (3-9) 庫倫力勢

能；在F3C model中，Asc為一補償係數 (原理請參考2.2節，數

值 列 於Table. 3-4) ，S函數為截斷勢能項 (Truncation shift function) ，f為電場收斂因子 (Converting factor) 。兩勢能的能

量參數 (如 i 0 , , and qr ε σ 等等) 對應不同原子種類列於Table. 3-4 及Table. 3-5。

( )

(

)

12

(

)

6

( )

1 0 2 0 vdw ij ij ij vdw ij U r =C r r −C r r −S r (3-8)

( )

(

)

12

(

)

6

(

)

(

)

12

(

)

6 1 0 2 0 1 0 2 0 1 2 6 2 ( , 2 ) vdw ij c c ij c c c c sc ij ij S r C r r C r r r r r C r r C r r C A ε C ε ⎡ ⎤ ⎡ ⎤ =_⎣ − _⎦ − − _⎣ − _⎦ = =

( )

i j Columb ij els ij U = f q q⎡_⎣ r −S r ⎤_{⎦ (3-9)}

( )

(

i j

) (

)

(

i j 2

)

els ij c ij c c S r = q q r − r −r q q r Table. 3-1 原子質量參數 原子種類 原子編號 質量(amu) 附註 氫 HW 1.008 水分子的氫 氧 OW 15.999 水分子的氧 鈉 Na 22.989 鈉離子 氯 Cl 35.450 氯離子

Table. 3-2 鍵結拉伸能量參數 拉伸鍵結編號 Kb(kcal/mol-Å2) R0(Å) OW-HW 250 1.000 Table. 3-3 鍵結彎曲能量參數 彎曲鍵結編號 Kθ(kcal/mol-rad2) θ0(Degree) HW-OW-HW 60 109.47 Table. 3-4 凡得瓦能量參數 原子編號 原始 0 r (Å) ε (kcal/mol) 調整(Asc=0.84 for 截斷半徑=6Å)1 0 r (Å) ε (kcal/mol) HW 0.900 0.01 0.8742 0.01190 OW 3.55322 0.18479 3.35260 0.26190 Na 2.7000 0.18479 2.62267 0.2200 Cl 4.7000 0.1000 4.56539 0.11905 Table. 3-5 各原子的帶電價數 原子編號 價數(e) HW 0.410 OW -0.820 Na +1 Cl -1 1 _{其他截斷半徑對應之 A 值請參考[15]}

3.3 計算與分析的方法

在得到作用力函式之後，就可以依照分子動力學模擬的流程圖， 如Fig. 3-1，進行鹽類水溶液的相關模擬研究。鹽類水溶液的模擬分子 有Sodium ion、Chloride ion以及水分子，將其混和於一個系統，可組成

帶負電系統 (Cl-+Water) 、帶正電系統 (Na++Water 以及兩種中性系統

(即純水或NaCl系統) 。由於計算步驟繁雜，此節僅將指出主要的模擬 流程步驟，細節將於後面章節一一說明： Fig. 3-2 水溶液系統初始化結果圖 (216個水分子及4個氯離子)。先設定模擬的水分 子數量，依據不同溫度下的實際密度，計算出模擬盒子的體積，採用FCC排列的方 式，計算出適當晶格長度初始化水分子位置座標；離子則排列於適當位置 (可與水 溶液系統產生作用力) 。實線表示週期性邊界，細節請參考式 (3-11) 與 (3-12) 的 說明。 1. 首先要初始化位置與速度，統計數據前，模擬系統需達到其組 態與溫度的平衡，然而組態取決於粒子 (包括分子基團或原子) 的相對位置，溫度取決於粒子的速度，在系統沒有特殊需求的 情況下，初始化條件只需給定粒子的隨機位置與速度。其中配 置系統粒子的位置之適當性，將會影響到系統的收斂穩定性 (Stability of converge ) 與否，即不適當的初始位置 (不符合真 實環境下的組態) 將會使系統無法在穩定狀態下進行模擬。因

此，我們將對於配置粒子的初始位置與速度之準則，在以下說 明。系統位置初始化，最簡易的初始化方法是結晶型態初始 化，常見三種堆積結構，分別是面心立方晶體結構 (FCC) 、 體心立方晶體結構 (BCC) 、六方最密堆積晶體結構 (HCC) 。 研究是以FCC初始化水溶液位置後再擷取需要的尺度大小，完 成初始化系統位置的工作，如Fig. 3-2所示。至於分子或原子的 速度初始化設定，最常使用的方法是採用亂數給定速度。一般 合乎真實物理狀態，會根據Maxwell分布函數來篩選亂數 (速度 值) 的給定，並將亂數限定在一個範圍內。此外，一般模擬上 會在初始化速度時，將系統的總動量歸零 (相對於質心速 度) 。先將所有粒子的初始動量累加，再平均求得每個粒子的 平均動量，將每個粒子減去此平均動量，則初始總動量即歸 零，如式 (3-10) 所示。 1 1 / N new i i i k k i k m m m N = ⎛ ⎞ =_⎜ − _⎟ ⎝

∑

⎠ v v v (3-10) 2. 系統的控制是根據模擬系統的類型，需針對此模擬系統狀態條 件分別作控制，控制的部分可對溫度 (T) 控制、壓力 (P) 控 制 或 體 積 (V) 控 制 。 本 研 究 為 正 則 系 綜 (Canonical ensemble) ，也就是NVT emsemble，模擬過程中，粒子數 (N) 、 系統體積 (V) 以及系統溫度 (T) 為固定值，其中系統體積與 系統溫度控制方法在下面說明： A. 系統體積控制：可藉由週期性邊界 (Periodic boundary condition) 或虛擬牆 (Virtual wall) 的邊界設定控制。週期 性邊界為分子動力學中經常被採用計算巨觀性質的邊界條 件，消除由於模擬系統邊界造成的表面分子影響效應 (Surface effect) ，Fig. 3-5顯示二維週期邊界示意圖，Virtual Cell之分子是由Main cell之分子映像而得 (位置、速度及加 速度等同於Main cell之粒子) ，其細節可分為：

i. 粒子位置週期性：當粒子運動超出系統，則依據式 (3-11)與式 (3-12) 將粒子平移至系統內部，如Fig. 3-3 以及Fig. 3-4所示。 / 2, / 2, ij ij new ij ij ij new ij r L r r r L r L r r r L > = = − ⎧ ⎨ < − = = + ⎩ (3-11) i new i i new i x L x x L x L x x L > = − ⎧ ⎨ < = + ⎩ (3-12) ii. 最小鏡向作用力：當粒子間距離超過截斷半徑，則考慮 Virtual與Main cell之粒子間作用力，如式 (3-11) ，rij表 示粒子間距離，L為系統邊長，其值設定大於2倍的截 斷半徑。如Fig. 3-5所示。 Fig. 3-3 二維週期性邊界作用示意圖。粒子位置從 R(t’)運動至 R(t”)，該位置超出 系統邊界。

Fig. 3-4 二維週期性邊界作用示意圖。粒子位置從R(t”)透過式 (3-11)或式 (3-12)平 移至最終R(t”)的位置，完成週期性邊界。

Fig. 3-5 二維週期性邊界作用示意圖。Main cell具有一水分子與另一粒子，水分子

與粒子間的作用力距離是最小鏡向距離，其條件參考式(3-11)。

3. 溫度控制：透過分子速度的調整，即可維持系統的溫度，藉由 Nose-Hoover thermostat[61]控制，方程式如式(3-13)至式(3-15)

所示。其中ζ稱為摩擦係數，g為系統自由度，Q為溫控質量 2 D Q=gkT τ i i i m = p v _(3-13) i = − ⋅i ξ i p F p (3-14) 2 1 i D i i d gkT dt Q m ζ ⎡ ⎤ = _⎢ − _⎥ ⎣

∑

⎦ p (3-15) ，τ 為控溫的特徵時間，利用式 (3-13) 至式 (3-15) ，可修正 每個分子的受力關係，進而達到修正系統的速度來控制系統的 溫度，就可將其值帶入數值積分法中 ，就可帶入下個Timestep 的計算。Nose-Hoover thermostat是NVT系統中最常被使用的方 法。此方法的優點在於能夠維持系統能量恆定，較能描述真實 分子運動之行為。上述是用於需控制溫度的系統才進行此步 驟，如果不需控溫，在其計算流程中就不需要此處理。 4. 分子受力的計算為分子動力學的核心，根據牛頓運動定律並配 合適當的分子勢能函式，計算出系統分子在當時所受加速度 a，以便對下個時間分子運動軌跡作預測。然而，分子動力學 的過程相當耗費計算時間，搜尋粒子計算距離是佔主要部分。 另外，模擬系統的大小和模擬時間的長短也使得需要龐大的計 算時間。可以預期，未來必定朝向分子數目龐大的系統，且模 擬的時間間隔數通常也在數千萬Timestep以上 (ns) ，因此適當 的簡化計算過程為分子動力學中重要的一環，否則很難將模擬 時間控制在合理的範圍內，如此一來，即使能正確模擬出重要 的現象，但無法在合理時間內完成模擬也是不合實際的。 5. 加速計算 A. 截斷半徑法：不論粒子間使用的勢能方程式為何種類型， 當粒子與粒子間的距離超過某一段距離後，粒子間的作用 力將會變的極小。因此，可以設定一個截斷距離 (Cut-off distance, Rcut) ，當粒子間的距離大於此距離時，直接假設

此兩粒子間的作用力為零，如此一來，就可節省要計算遠 距離的粒子作用力，且因於截斷距離外之粒子作用力影響 微乎其微，對粒子預測幾乎不造成影響，卻由能有效簡化 龐大的計算量；大部分的分子動力學文獻，都會採用勢能 截斷法來進行模擬。截斷距離長度的選擇，是由勢能函式 的變化程度決定的，如模擬系統僅包含簡單分子，且不考 慮靜電效應與庫倫力時，通常所使用的截斷距離約在 2.5σ~3.0σ之間，σ為兩分子間平衡距離；反之，若考慮具 有電荷的粒子行為時，一般則截斷距離會選擇在5σ以上甚 至超過10σ，不過沒有一定的設定值，主要還是取決於研究 文獻所提供的參數。 B. 加速搜尋法：分子動力學模擬最耗時的部分在於判斷兩分 子間距離，當系統特別龐大時，判斷會更為耗時，運算量 會隨著分子數目的平方成正比關係，即便是分子勢能截斷 法，也需先求得分子間距離。因此一旦系統內分子數過多 時，則會大量拖累運算速度；於是在Verlet (1967) 提出鄰 近列表法來減少判斷分子間距離之運算 (另外，可搭配Cell link法[62]，因本研究系統較小，加速效果有限，故不使用 與說明) 。 i. 鄰近列表法的基本觀念：在系統中的任一分子，其周圍 的分子於短時間內變動 (一般為5fs內) 並不會太大，意 即於此段時間內，此分子僅需和周圍鄰近固定的分子計 算作用力，可以不用對整個系統的分子作計算；當此方 法使用在越大的系統時，加速的效果更為顯著。 ii. 鄰近列表法的使用方法：先將1號粒子列表半徑 (List distance ,rl) 內的粒子都列入清單中 (如Fig. 3-6) ，直到 重新更新列表名單之前，1號粒子只與清單上的粒子計 算作用力。鄰近列表法可控制更新表單的時間間隔，在

計算效率與準確度之間取得適當的平衡。其中列表半徑 大小和更新列表的頻率是相互影響的，一般而言，當粒 子的移動速度越快，更新頻率則要越高或者是列表半徑 要增大。 iii. 鄰近列表法示意圖說明：如Fig. 3-6所示，以1號分子為 中心，Dash-Dash圓圈為1號粒子的截斷半徑範圍， Dash-Dot-Dot圓圈則為1號粒子列表半徑的範圍。在更 新表單前，1號分子要計算的粒子編號為2號以及3號， 這表示2號與3號在更新表單前大概都會在1號粒子的列 表半徑範圍內，而4號粒子可能無法進入截斷半徑內， 因此不考慮列入計算。 Fig. 3-6 鄰近列表法示意圖。此為分子動力學模擬中常用之加速 技巧之一；如果在更大尺寸系統，通常會搭配Cell link 法。 6. 解運動方程式的數值積分法：本研究已知初始條件，粒子位置 和速度 ，再加上由計算粒子的作用力所得的加速度 ，採用 Gear的五階Predictor-Corrector algorithm之數值方法，即可預測 在下一個步階時間下，求得粒子的新位置與速度。其中Gear的 五階Predictor-Corrector algorithm的原理主要是利用五階泰勒

級數之展開式，預測每個粒子的位置和速度，分別對應式 (3-16) 及式 (3-17) 。再根據此二式，推導出Gear Predictor-Corrector algorithm Gear的五階Predictor-Corrector algorithm的處理過程 可 以 分 成 三 個 部 分 ， 分 別 為 預 測 (Prediction) 、 計 算 (Evaluation) 作用力和修正 (Correction) ： A. 預測：首先使用泰勒展開式對位置r( )t、速度v( )t 、加速度r&&( )t 以及加入 3 個時間微分項 ( )iii ( )t r ，r( )iv ( )t ，r( )v ( )t 作下個時間 步階之預測，其計算展開式為矩陣，式 (3-18) 表示之。 B. 計算作用力：由牛頓第二運動定律求出在新預測位置 (t+ Δt) r 上每個粒子所受到的作用力F(t+ Δt)，並轉換為粒子 加速度&&r(t+ Δt)，式 (3-19) 與式 (3-20) 表示作用力與加速 度的計算。 C. 修正：將Predictor步驟計算出的r&&(t+ Δt)與計算作用力的步 驟中的&&r(t+ Δt)相減，求出兩者之間的差值Δ&&r(t+ Δt)，如下 所示，因此我們可以獲得到一個修正項，如式 (3-21) 及 (3-22)所示。 D. 最後，我們可以將第一步驟所得到的預測值加以修正，如 式 (3-23) 所示：其中為α0、α1、α2、α3、α4、以及α5為 Gear Predict-Corrector的參數，會隨著演算法計算階數 (q) 不同而改變參數值，其數值列於Table. 3.1。 2 3 4 5 (iii) (iv) (v) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2! 3! 4! 5! t t t t t+ Δ =t t + t Δ +t t Δ + t Δ + t Δ + t Δ r r r& r&& r r r (3-16) 2 3 4 (iii) (iv) (v) ( ) ( ) ( ) ( ) ( ) ( ) 2! 3! 4! t t t t+ Δ =t t + t Δ +t t Δ + t Δ + t Δ

r& r& r&& r r r (3-17)

( ) ( ) ( ) ( ) ( ) 2 3 4 5 2 3 4 2 3 2 2! 3! 4! 5! 2! 3! 4! 2! 3! iii iii 2! iv iv v ( ) ₁ ( ) ( ) _{0 1} ( ) ( ) _{0 0 1} ( ) ( ) _{0 0 0 1} ( ) 0 0 0 0 1 ( ) ( ) 0 0 0 0 0 1 ( ) t t t t t t t t t t t t _t t t t _t t t t _t t t t _t t t t t t t t Δ Δ Δ Δ Δ Δ Δ Δ Δ Δ + Δ ⎡ _{⎤ ⎡ Δ ⎤} ⎢ _{+ Δ ⎥ ⎢ ⎥} Δ ⎢ _{⎥ ⎢ ⎥} ⎢ + Δ ⎥ _{⎢ Δ ⎥} ⎢ _{⎥ = ⎢ ⎥} + Δ ⎢ ⎥ _⎢ Δ _⎥ ⎢ _{+ Δ ⎥ ⎢ Δ ⎥} ⎢ _{⎥ ⎢ ⎥} ⎢ _{+ Δ ⎥ ⎢⎣ ⎥⎦} ⎣ ⎦ r r v v r r r r r r r r && && ( )v _{( )} t ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ (3-18)

(上標羅馬數字代表微分次數) ( ( )) (t+ Δ = −t) ∂U t+ Δt ∂ r F r (3-19) ( ) ( ) i t t t t m + Δ + Δ =F r && (3-20) ( ) [ ( ) P( )] i i t t t t t t

Δ&&r + Δ = r&& + Δ −&&r + Δ (3-21)

2 2 t Δ Δ ≡R Δr&& (3-22) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 3 3 4 4 5 5 0 1 2 2 ₂ iii iii 3! 3! 3 iv iv 4 4! 4! v v 5 5! 5! t t t t t t t t t t α α α α α α Δ Δ Δ Δ Δ Δ Δ Δ ⎡ ⎤ ⎡ _{⎤ ⎡ ⎤} ⎢ _Δ ⎥ ⎢ _{Δ ⎥ ⎢ ⎥} ⎢ ⎥ ⎢ _{⎥ ⎢ ⎥} ⎢ ⎥ ⎢ _{⎥ ⎢ ⎥} ⎢ ⎥ ⎢₌ ⎥_{+⎢ ⎥Δ} ⎢ ⎥ ⎢ _{⎥ ⎢ ⎥} ⎢ ⎥ ⎢ _{⎥ ⎢ ⎥} ⎢ ⎥ ⎢ _{⎥ ⎢ ⎥} ⎢ ⎥ ⎢ _{⎥ ⎢ ⎥⎣ ⎦} ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ r r v v r r R r r r r r r && && (3-23)

Table 3-1 三～五階 Gear Predict-Corrector Algorithm 參數表

b i α b q = 3 q = 4 q = 5 0 α B 6 1 12019 163 1 α B 6 5 4 3 360251 2 α B 1 1 1 3 α B 3 1 2 1 18 11 4 α B － 121 61 5 α － － ₆₀1 7. 計算性質參數:在計算上預測粒子運動行為時，可以得到每個步 階時間下每個粒子的r、v、a等相關數值，將這些數值經過適 當的統計方法之平均，就可以得到巨觀性質如結構、熱力性質 等…。分子動力學可以準確利用牛頓定律模擬粒子行為狀態， 但所得資訊原先僅有表觀狀態，即粒子運動軌跡狀態而非巨觀 所觀測到的性質，因此需要經過統計力學相關推導才能引入巨 觀的性質行為，其定義如下所示： A. 溫度：對於溫度的計算，已在前面描述 NVT 系統中有所描

述。本研究所模擬的系統，皆在密閉系統中。 B. 壓力：在分子動力學模擬中，自一個給定的作用力場計算 瞬間的壓力張量的過程，是想像一個截斷平面上有許多粒 子通過，計算其動量流量 (Momentum flux) ，如式 (3-24) 所示，其中 1 1 2 N a i i i i K _β m v v_{α β} = =

∑

以及 1 N i i i W_αβ r F_{α β} = =

∑

(V 為體積； i m, r_i, v_i, 和F_i分別為i粒子的質量、位置、速度以及作用 力。W為維里項，α, β落在x, y , z的範圍) ，最後推得式 (3-25) 。

(

)

1 2 P_αβ = V K_αβ +W_αβ (3-24)

(

2

)

1 1 3 N _{i i i i} i P V m v r F = =

∑

+ ⋅ _(3-25)

C. 徑向分佈函數(Radial distribution function,g r

( )

)與配位數

(Hydration number, CN)：式 (3-26)表示一個粒子質心位在

距另一個粒子質心r+ Δr之球殼之內的機率，其中N為粒子

數目，V為球殼體積。而式 (3-27)顯示出鄰近球殼(r+ Δr)

出現幾個粒子，而徑向分佈函數可顯示出一個粒子在時間 平均下，鄰近原子的位置分佈情況，如Fig. 3-9 所示。

Fig. 3-7 徑向分佈函數物理示意圖。若以離子為中心，依據方程式統 計相關位置關係，即可得到溶液結構性質。

( ) (

,

)

1g r V r r N ρ Δ − (3-26)

( )

2

(

,

)

(

,

)

g r = N r Δr N V rρ Δr (3-27) 從 模 擬 資 料 計 算 的 該 函 數 ， 則 將 式 (3-27) 改 寫 為 式 (3-28) ，其中，N_K表示每個Timestep統計的結果，M表示 統計M個Timesteps。而配位數是從g r

( )

積分而得，式 (3-29) 表示配位數。

( )

(

)

(

)

1 1 , , 2 M K k g r N r r M N V rρ r = =

∑

Δ Δ _(3-28)

( )

1 ₂ 0 4 R _io CN = πρ

∫

g r r dr (3-29) (CN) 的計算方法，其中gio代表離子與水分子中的氧。積 分半徑設在第一水合層波谷 (rv1) ，就可以明白計算出離 子與水分子間或水分子與水分子間的水合配位數。CN 與

( )

g r 提供可分析溶液結構的結果，且都能近乎實驗結果 (實 驗的方式是以X-Ray 繞射與中子散射取得材料結構) 。 D. 擴散係數 (Diffusion coefficient) ：以愛因斯坦方程式計算 擴散係數，如式 (3-30) 所示，在微特徵性質上是重要的參 數，如奈米孔洞材料中，流體的擴散係數是重要的參數。 其計算方式列於式 (3-30) ，需要較長的模擬時間，才能取 得足夠的數據，本研究系統大於需180ps以上方能收斂。式

(3-30) 中， 代表累計次數的平均，而且

( ) ( )

0 2 N i i r t −r ⎡ ⎤ ⎣ ⎦

∑

一

式就是方均根位移 (Mean square displacement, MSD) 。

( ) ( )

2 1 lim 0 6 N i t i d D r t r N →∞dt =

∑

⎡_⎣ − ⎤_{⎦ (3-30)} Fig. 3-8 方均根位移統計的方法。是對各分子取平均，另外再對時間 取平均而得。 E. 臨界性質的計算：為了取得氣液相共存的模擬結果，研究 中參考Liew[63]等人 (1998) 模擬的方法，初始化系統如 Fig. 4-8 所 示 ， 將 不 同 溫 度 的 模 擬 數 據 ， 擬 合 到 一 個 Hyperbolic tangent function，如式 (3-31) 所示，ρV 與ρL表

示氣相密度與液相密度，z是Gibb’s dividing surface、以及d 是液氣介面的厚度參數。在得到非等向性 (Orthobaric) 的 密度分佈數據，就能進一步計算水分子的臨界狀態的性 質，最後將密度資料代入到The law of rectilinear diameter 以及 The scaling law of density 方程式中，便是式 (3-32) 及式 (3-33) ，其中Tc為臨界溫度，ρc為臨界密度，C C1, 2,βc 為擬合參數[63]，即可得到臨界狀態的密度與溫度。

( )

z 0.5

(

L V

)

0.5

(

L V

)

tanh[

(

z z0

)

/ ]d ρ = ρ +ρ − ρ −ρ − (3-31)

(

)

1/ 1 c c L V T− =T C ρ −ρ β (3-32)

(

ρL+ρV

)

/ 2=ρc+C2

(

1−T T/ c

)

(3-33) F. 速 度 自 相 關 函 數 (Velocity-Velocity Correlation Function) ：考慮速度在時間t與0時的相關程度，統計V(t) 與V(0)的乘積關係，若系統最後為混亂的，則CV=0 (時間0 與時間t的狀態無關) ，CV = 1 則表示高度相關 (High correlation) 。如式 (3-34) 所示。

( )

_{( ) ( )}

( ) ( )

0 0 0 i i v i i v t v C t v v ⋅ = ⋅ (3-34)

G. 氫 鍵 時 間 相 關 函 數 (Hydrogen bond time correlation function) ：可分為連續氫鍵相關函數 (Continue hydrogen bond correlation function) ，如式 (3-35) ；間歇性氫鍵相關 函數 (Intermittent hydrogen bond correlation function) ，如

式(3-36)。... 表示對全部離子-水分子或水分子-水分子的配 對數目取平均，h t

( )

=1表示在時間t時，氫鍵 (或是Chloride ions對氫原子的鍵結) 存在，反之h t

( )

=0；若H t

( )

=1則表 示該鍵結從初始時間到時間t皆存在，反之則H t

( )

=0。S函 數的收斂時間，可代表鍵結的平均存在時間 (Bond average lifetime) ；C函數的收斂時間，可視為水分子 (或離子) 穿 越可生成鍵結 (與其他水分子) 的區域，到達不會產生鍵 結 的 區 域 ， 所 花 費 的 時 間 ， 亦 為 鍵 結 的 鬆 弛 時 間 (Relaxation time) 。

( )

h 0 H

( ) ( )

_{( ) ( )}

S h 0 h 0 HB t t = ⋅ ⋅ (3-35)

( )

h 0 h

( ) ( )

_{( ) ( )}

C h h HB t t t t ⋅ = ⋅ (3-36) H. 氫鍵與Cl..H鍵結的計算：為了進一步瞭解水溶液的結構， 氫鍵或是Cl…H鍵提供了相當重要的參考依據，因為水溶 液特殊的行為與性質大多是因為氫鍵的存在。氫鍵的幾何 限制定義如2.3所述，分為三個條件：

i. O…O距離小於3.4Å，或者Cl…O距離小於3.95Å。

ii. O…H的距離小於2.4Å，或者Cl…H距離小於2.95Å。

iii. α (α=∠HOO) 小於30°，或者∠HClO小於30°