在得到作用力函式之後，就可以依照分子動力學模擬的流程圖，
如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) iP(t t)]
Table 3-1 三～五階 Gear Predict-Corrector Algorithm 參數表
bαib 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
( ) ( )
0 2 中參考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 CorrelationFunction) ：考慮速度在時間t與0時的相關程度，統計V(t) 與V(0)的乘積關係，若系統最後為混亂的，則CV=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)。