水溶液的奈米液滴(Aqueous Nanodroplets)

用於高分子溶液的電泳分離(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章、研究方法

粒子在相空間的軌跡 (Phase-Space Trajectories) 變化，其有幾點假設：

1. 所有粒子 (包含分子基團或原子) 皆遵守古典牛頓力學的運動獻所提出的勢能模型 (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) 中K_s為伸縮常

stretch b eq

U =K r−r (3-1)

Bending i eq

F⁺ = k_θ θ θ w t (3-5)

( )( )

Bending Bending Bending

F +F⁺ +F⁺ =

Q (3-6)

(

)

1 i i

Bending Bending Bending

F ⁺ F F ⁺ function) ，f為電場收斂因子 (Converting factor) 。兩勢能的能量參數 (如ε σ, , and qr₀ ⁱ等等) 對應不同原子種類列於Table. 3-4

Columb ij els ij

U = f q q⎡⎣ r −S r ⎤⎦ (3-9)

Table. 3-2 鍵結拉伸能量參數

拉伸鍵結編號 Kb(kcal/mol-Å²) R0(Å)

OW-HW 250 1.000

Table. 3-3 鍵結彎曲能量參數

彎曲鍵結編號 Kθ(kcal/mol-rad²) θ₀(Degree)

HW-OW-HW 60 109.47

Table. 3-4 凡得瓦能量參數 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 各原子的帶電價數

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 ) 與否，即不適當的初始位置 (不符合真實環境下的組態) 將會使系統無法在穩定狀態下進行模擬。因

此，我們將對於配置粒子的初始位置與速度之準則，在以下說 ensemble) ，也就是NVT emsemble，模擬過程中，粒子數 (N) 、系統體積 (V) 以及系統溫度 (T) 為固定值，其中系統體積與系統溫度控制方法在下面說明：

A. 系統體積控制：可藉由週期性邊界 (Periodic boundary condition) 或虛擬牆 (Virtual wall) 的邊界設定控制。週期性邊界為分子動力學中經常被採用計算巨觀性質的邊界條件，消除由於模擬系統邊界造成的表面分子影響效應 (Surface effect) ，Fig. 3-5顯示二維週期邊界示意圖，Virtual Cell之分子是由Main cell之分子映像而得 (位置、速度及加速度等同於Main cell之粒子) ，其細節可分為：

i. 粒子位置週期性：當粒子運動超出系統，則依據式

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為溫控質量擬的時間間隔數通常也在數千萬Timestep以上 (ns) ，因此適當的簡化計算過程為分子動力學中重要的一環，否則很難將模擬

此兩粒子間的作用力為零，如此一來，就可節省要計算遠時，則會大量拖累運算速度；於是在Verlet (1967) 提出鄰近列表法來減少判斷分子間距離之運算 (另外，可搭配Cell

計算效率與準確度之間取得適當的平衡。其中列表半徑大小和更新列表的頻率是相互影響的，一般而言，當粒子的移動速度越快，更新頻率則要越高或者是列表半徑要增大。

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 Gear Predict-Corrector的參數，會隨著演算法計算階數 (q) 不同而改變參數值，其數值列於Table. 3.1。

(上標羅馬數字代表微分次數)

r&& (3-20)

(t t) [ (_i t t) _i^P(t t)]

Table 3-1 三～五階 Gear Predict-Corrector Algorithm 參數表

bαi^b q = 3 q = 4 q = 5

述。本研究所模擬的系統，皆在密閉系統中。

(Hydration number, CN)：式 (3-26)表示一個粒子質心位在距另一個粒子質心r+ Δr之球殼之內的機率，其中N為粒子數目，V為球殼體積。而式 (3-27)顯示出鄰近球殼(r+ Δr) 出現幾個粒子，而徑向分佈函數可顯示出一個粒子在時間平均下，鄰近原子的位置分佈情況，如Fig. 3-9 所示。

Fig. 3-7 徑向分佈函數物理示意圖。若以離子為中心，依據方程式統 D. 擴散係數 (Diffusion coefficient) ：以愛因斯坦方程式計算

擴散係數，如式 (3-30) 所示，在微特徵性質上是重要的參數，如奈米孔洞材料中，流體的擴散係數是重要的參數。

其計算方式列於式 (3-30) ，需要較長的模擬時間，才能取得足夠的數據，本研究系統大於需180ps以上方能收斂。式

(3-30) 中，代表累計次數的平均，而且 ^N i

( ) ( )

⁰ ² 中參考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)

(

ρ_L+ρ_V

)

/ 2=ρ_c+C2

(

1−T T/ _c

)

(3-33) F. 速度自相關函數 (Velocity-Velocity Correlation

Function) ：考慮速度在時間t與0時的相關程度，統計V(t) 與V(0)的乘積關係，若系統最後為混亂的，則C_V=0 (時間0 function) ：可分為連續氫鍵相關函數 (Continue hydrogen bond correlation function) ，如式 (3-35) ；間歇性氫鍵相關函數 (Intermittent hydrogen bond correlation function) ，如式(3-36)。... 表示對全部離子-水分子或水分子-水分子的配

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

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

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

Fig. 3-9 氫鍵形成的幾何限制條件。參考 Ju[15]等人(2006)。

第4章、結果與討論

4.1 水分子之徑向分佈函數的模擬結果

首先，比較研究模擬所得與F3C model文獻[2, 31-34]的徑向分佈函數結果。分佈結果表示水分子的區域分佈情形，代表溶液系統的結構性質，而系統的空間分佈是由作用力決定。相同的分佈結果能驗證模擬程式與文獻無誤，接者程式才能供後續研究使用。水分子統計結構性質，分別氧-氧原子、氧-氫原子及氫-氫原子的徑向分佈函數。

模擬系統有216顆水分子，密度設定如Table. 4-1，而各溫度下的系統密度，請參考Table. 4-2。以Nose-Hoover thermostat控制溫度，系統為三維週期性邊界。截斷半徑為6-8 Å、A_sc=0.84以及積分的步進時間為1fs至0.5fs。

Fig. 4-1為不同系統狀態下，壓力統計結果，與文獻[2]的結果相當接近。Fig. 4-2、Fig. 4-4以及Fig. 4-6，分別為不同溫度下 (298K、423K 及573K) ，氧對氧、氫對氫以及氫對氧的徑向分佈函數結果。比較研究所得結果，Fig. 4-3、Fig. 4-5以及Fig. 4-7 ，可驗證模擬結果與文獻結果相符合。藉由Gr的結果，積分得到離子水合數(配位數)，並且與文

在文檔中水及鹽類水溶液之熱力學性質、鍵結行為及系統結構之探討:分子動力學模擬的研究 (頁 29-0)

第3章 、研究方法