局部結構對液態水震動的影響：INM分析

國 立 交 通 大 學

物理研究所

碩 士 論 文

局部結構對液態水振動的影響：INM 分析

Local-structure effects on vibrations of liquid water：INM analysis

研 究 生：陳郁竣

指導教授：吳天鳴 教授

局部結構對液態水振動的影響：INM 分析

Local-structure effects on vibrations of liquid water：INM analysis

研 究 生：陳郁竣

Student：Yu-Chun Chen

指導教授：吳天鳴 教授

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

January 2013

Hsinchu, Taiwan, Republic of China

局部結構對液態水振動的影響:：INM 分析

學生：陳郁竣

指導教授：吳天鳴 教授

國立交通大學物理研究所碩士班

摘

要

在這篇論文中，我們將運用分子動力模擬(MD simulation)與瞬間正則模分析研 究水分子的振動行為。此外，我們使用 voronoi 分析及氫鍵狀態分析的方法來探討 局部區域結構對水分子振動行為的影響。 我們的結果顯示，隨著水分子的 voronoi asphericities 增加，在水分子的彎曲運 動頻譜(bending spectra)與伸縮運動頻譜(stretching spectra)分別有藍位移與紅位移 的現象。水分子的氫鍵狀態不同也會影響液態水的振動行為。

Local-structure effects on vibrations of liquid water：INM analysis

Student：Yu-Chun Chen

Advisor：Dr. Ten-Ming Wu

Institute of Physics

National Chiao Tung University

ABSTRACT

In this thesis, we study the intermolecular and intramolecular

vibrations of liquid water through molecular dynamics(MD) simulation and

the instantaneous normal mode (INM) analysis. In addition, we use voronoi

polyhedra(VP) analysis and hydrogen bond configuration analysis to study

the local-structure effects on vibrations of water molecules.

My results indicate that a blue shift and a red shift are found in the

bending and the stretching spectra of liquid water, respectively, when the

voronoi asphericities of water molecules increase. Different H-bond

configurations of water molecules actually influence both intermolecular

and intramolecular vibrations of liquid water.

誌

謝

目

錄

中文摘要 ……… i 英文摘要 ……… ii 誌謝 ……… iii 目錄 ……… iv 圖目錄 ……… v 一、 序論 ……… 1 二、 理論 ……… 2 I. 模型 ……… 2 A. 剛性 SPC/E 水分子模型 ……… 2 B. 非剛性 SPC/E 水分子模型 I ……… 3 C. 非剛性 SPC/E 水分子模型 II ……… 4 II. 振動頻譜 ……… 5 A. 速度頻譜 ……… 6 B. 瞬間正則模分析 ……… 6 III. 局部結構的分析 ……… 8 A. 徑向分佈函數 ……… 8 B. Voronoi 分析 ……… 9 C. 氫鍵狀態分析 ……… 12 三、 分子動力模擬 ……… 17 四、 結果分析 ……… 19 I. 局部結構與液態水的 INM 頻譜 ……… 19 A. 徑向分佈函數 ……… 19 B. 速度頻譜 ……… 20 C. INM 頻譜 ……… 22

II. Subsemble by Voronoi polyhedral analysis ……… 25

A. 徑向分佈函數 ……… 25

B. INM 頻譜 ……… 28

III. Subsemble by Hydrogen bond configuration analysis ……… 31

A. 徑向分佈函數 ……… 31

B. INM 頻譜 ……… 33

五、 結論 ……… 45

附錄 A ……… 46

圖

目

錄

圖(一) 水在 300K 的 asphericity 分佈曲線 ……… 10 圖(二) 兩個水分子間判定氫鍵的示意圖 ……… 12 圖(三) 非剛性水模型 I 四種 VG 的氫鍵狀態分佈圖 ……… 15 圖(四) 非剛性水模型 II 四種 VG 的氫鍵狀態分佈圖 ……… 16 圖(五) 電腦模擬所採用的週期性邊界條件示意圖 ……… 17 圖(六) 水在 1.0gcm3&T 300K的徑向分佈函數 ……… 19 圖(七) 非剛性水模型 II 的歸一化速度相干函數 ……… 20 圖(八) 模型 I 和模型 II 的速度頻譜 ……… 21 圖(九) 非剛性模型 II 的 INM 狀態密度 ……… 23 圖(十) 非剛性模型 I 和非剛性模型 II 的 INM 頻譜比較 ……… 24 圖(十一) 模型 II 分別對氫、氧原子作投影的速度頻譜跟 INM 頻 譜……… 25 圖(十二) 模型 II 四個 VG 的徑向分佈函數 ……… 26 圖(十三) 非剛性模型 I 及非剛性模型 II 的氧-氧徑向分佈函數 … 27 圖(十四) 非剛性模型 I 及非剛性模型 II 的氧-氫徑向分佈函數 … 27 圖(十五) 非剛性模型 I 及非剛性模型 II 的氫-氫徑向分佈函數 … 28 圖(十六) 模型 II 不同 VG 的歸一化 INM 狀態密度 ……… 29 圖(十七) 模型 I 不同 VG 的歸一化 INM 狀態密度 ……… 30 圖(十八) 非剛性模型 I 及非剛性模型 II 的氧-氧徑向分佈函數 … 31 圖(十九) 非剛性模型 I 及非剛性模型 II 的氧-氫徑向分佈函數 … 32 圖(二十) 非剛性模型 I 及非剛性模型 II 的氫-氫徑向分佈函數 … 33 圖(二十一) 非剛性模型 I 氫原子各種 Config. S 的歸一化 INM 狀態密 度……… 34 圖(二十二) 非剛性模型 II 氫原子各種 Config. S 的歸一化 INM 狀態密 度……… 35

圖(二十三) 氧原子各種 Config. S 的 imahinary prat 及平移旋轉運動頻譜─donor H-bond ……… 36

圖(二十四) 氧原子各種 Config. S 的 imahinary prat 及平移旋轉運動頻譜─ acceptor H-bond ……… 37

圖(二十五) 氧原子各種 Config. S 的 imahinary prat 及平移旋轉運動頻譜─total H-bond ……… 38

圖(二十六) 氧原子各種 Config. S 的彎曲運動頻譜─donor H-bond … 39

圖(二十七) 氧原子各種 Config. S 的彎曲運動頻譜─acceptor H-bond … 40

圖(二十八) 氧原子各種 Config. S 的彎曲運動頻譜─total H-bond … 41

圖(三十) 氧原子各種 Config. S 的伸縮運動頻帶─acceptor H-bond … 43

第一章 序論

水是一種無所不在且具特殊性質的液體，最廣為人知的就是水的密度在4C時 有最大值。水的眾多特殊性質與其複雜的氫鍵網絡(Hydrogen-bond network)有極大 的關係。雖然水已經被研究了很長一段時間，我們對它的氫鍵網絡的微觀動力學 仍然所知有限。 近來，由於 Mid-IR pump-probe 技術的發展，提供了一個有效的方法，探討水 的微觀動力學行為。在 2003 年，分析水的 Mid-IR pump-probe 實驗 [1]，了解到溶 解在D₂O液體中的HDO分子，它的OH原子間的瞬間振盪頻率，會受兩個因素影 響。首先，該HDO分子的OH與相鄰D₂O分子的氧原子間是否形成氫鍵的鍵結 (O─HO)。第二，OH與D₂O之間形成氫鍵(O─HO)，當氫鍵的長度越短，OH 原子間的瞬間振盪頻率也越低 [2]。 我們有興趣的是，純水中水分子的OH振動(vibrational motion)如何受氫鍵網絡 的影響。不同的氫鍵鍵結，會導致水分子的局部結構(local structure)不相同。在這 篇論文中，我們使用 voronoi 分析(voronoi polyhedra analysis)及氫鍵狀態分析 (Hydrogen-bond configuration analysis)將水分子依不同的局部結構作分類，然後使 用瞬間正則模分析(instantaneous normal mode analysis)探討不同局部結構水分子的 瞬間振動頻譜，我們就可以知道水分子的振動跟水分子的局部結構間有什麼對應 關係。 本篇的論文架構如下。在第二章會先介紹我們模擬所採用的水分子理論模型 以及計算頻譜的方法，並介紹兩種不同局部結構分析的原理，與這兩種分析方法 間的關係。第三章將交代我們分子動力模擬的條件，以及相關的演算法。第四章 的部分則討論分子動力模擬的結果，包括水的整體性質與不同局部結構對振動頻 譜的影響，並比較不同模型間的頻譜差異。第五章則是本文的總結。

第二章 理論

I. 模型 (Model)

水是地球上很常見的物質，也是構成生命體的主要成分之一，各種關於水的 性質的研究是持續不斷在進行的。其中，電腦模擬在研究水的性質上是很有效的 方法，許多不同的水分子模型已經被提出，並且應用在水的分子動力模擬中，探 討水的結構以及動力學性質，如 ST2, TIP4P, SPC, SPC/E…等等 [3][4][5]。選擇一 適當的水分子模型，我們可以成功的解釋液態水的一些物理特性。 在這篇論文中，我們會介紹三種水分子模型，一種是剛性水分子模型(rigid water model)，該模型將水分子的氫原子跟氧原子固定在一三角形上，原子跟原子 間的距離保持不變，它有固定的鍵長跟鍵角。另外我們討論兩種非剛性水分子模 型(flexible water model)，這兩種非剛性模型原子間的距離不是固定的，使用不同的 分子內位能描述原子跟原子間的相互運動 [6][7]。這三種水分子模型都是建立在 SPC/E 模型(SPC/E model)的基礎上。在本文中，我們主要是想了解水分子內部原 子間的振動行為，對於剛性水分子模型，由於鍵長跟鍵角都是固定的，水分子內 部的原子彼此之間不會振動，只有水分子間的振動。但剛性水分子模型仍然可以 當作一個參考的標準，探討分子內部原子間位能對液態水的結構造成的影響。

A. 剛性 SPC/E 水分子模型 (Rigid SPC/E water model)

剛性 SPC/E 水分子模型沒有分子內部的作用力，系統的總位能只有水分子─ 水分子間位能(Intermolecular potential)，由不同水分子的氫、氧原子之間的庫倫位 能(Coulomb potential)與一描述氧原子─氧原子的 Lennard-Jones potential 所組成， 如下：





                       j i iOjO rigid jO iO rigid j i i j j i Inter r N r M r q q V ₆ , 12 , 3 1 3 1 4 0 ,        (2 - 1)

  j i

r , 表示第i個水分子的 原子與第 j 個水分子的  原子間的距離，qi是第i個分

子的 原子所帶的電荷量，ri_{O j, O}是第i個水分子的氧原子與第 j 個水分子的氧原子

間的距離。M_rigid跟N_rigid 分別代表 Lennard-Jones potential 的強度參數， 是真空₀

中的介電常數。表(一)列出剛性水分子模型所用的各項參數值。 模型參數 數值 模型參數 數值 OH r _1.0A HOH  47 . 109 O

q 0.8476 e M_rigid 103 629.4 kcalÅ12/mol

H

q 0.4328 e N_rigid 625.5 kcalÅ6/mol

表(一) 剛性水分子模型的參數

B. 非剛性 SPC/E 水分子模型 I (Flexible SPC/E water model I) [6]

非剛性水分子模型 I (model I)是在 SPC/E 模型的基礎上，加上分子內部不同原

子間的位能，來描述水分子OH鍵長的伸縮與HOH鍵角的變化。系統的總位能 由水分子─水分子間位能及水分子內部的原子─原子間位能(Intramolecular potential)所構成，如下： Intra Inter Total V V V   (2 - 2) 模型 I 的分子-分子間位能，與剛性水分子模型相同，也是由氫、氧原子之間 的庫倫位能與氧原子─氧原子的 Lennard-Jones potential 所組成，僅所使用之參數 不同，如下：





               j i iOjO I jO iO I j i i j j i Inter r N r M r q q V ₆ , 12 , 3 1 3 1 4 0 ,        (2 - 3) 在這裡，

r

i_,j_代表第i個分子的 原子與第 j 個分子的  原子間的距離，

q

i是第 i個水分子的 原子所帶的電荷量。M 和_I N 是描述模型 I 的 Lennard-Jones potential_I

強度的常數。 水分子內部的位能形式如下：























_                 i _eOH _{i e e}HOH _{i e} e i OH e e i OH e e i OH e Intra k x x k k x x k x x k k x x k V 2 3 2 2 2 3 1 2 1     (2 - 4) 1 i x 是第i個水分子的氧原子與第一個氫原子的距離，x 是第_i2 i個水分子的氧原子與 第二個氫原子的距離，x 是平衡時氧原子與氫原子的鍵長，_e  是第_i i個水分子內部 HOH  所夾的角度， 是平衡時_e HOH的鍵角。表(二)中列出非剛性模型 I 的位 能及系統之參數。 模型參數 數值 模型參數 數值 H q 0.4328 e NI 625.5172kcalÅ /mol 6 O q 0.8476 e keOH 2 / / 06 .

547 kcal mol Angstrom

e   43 . 109 k 1 65 . 1   Angstrom e x 1.0 Angstrom HOH e k 2 / / 9.81

4 kcal mol rad

3 10  I M 630.795 kcalÅ12/mol 表(二) 非剛性水分子模型 I 的參數

在非剛性模型 I 中，HOH鍵角的彎曲運動(Bending motion)只用簡諧運動

(Harmonic motion)描述；氫氧原子間鍵長的伸縮運動(Stretching motion)則是由諧和 項(Harmonic term)以及非諧和項(Anharmonic term)組成。在模型 I 中並沒有考慮鍵 角的彎曲與鍵長的伸縮彼此所造成的相互影響，也就是我們將兩者的運動視為互 相獨立的運動。

C. 非剛性 SPC/E 水分子模型 II (Flexible SPC/E water model II) [7]

非剛性水分子模型 II (model II) 的總位能也是由分子─分子間位能及分子內 部原子─原子間位能所構成。非剛性水模型 II 的分子-分子間位能的形式與剛性模 型及非剛性水模型 I 一樣，也是只有 Lennard-Jones potential 的參數與前兩個模型 的有些微差異：





                       j i iOjO II jO iO II j i i j j i Inter r N r M r q q V ₆ , 12 , 3 1 3 1 4 0 ,        (2 - 5) II M 和N 是模型 II 的 Lennard-Jones potential 強度的常數。 _II 在非剛性模型 II 中，分子內部的位能由四個部分組成：第一項是 Morse potential，是模擬OH鍵長伸縮的位能，第二項是簡諧振盪，描述HOH鍵角的彎 曲行為，第三、四項是鍵長與鍵角變化相互影響的交叉項(Stretch-bend coupling term)：













ae H iH e iOH e iOH i rr r l l e Intra r r r r r r r r r r r k r r r k r k r a D V i i i i i i i i i i                             

_{ }

 2 1 2 1 2 3 1 2 1 2 1 3 2 3 2 1 2 2 1 exp 1 _{ } (2 - 6) 2 1 iOH iOH r r 、 分別是第i個分子氧原子與第一個氫原子、第二個氫原子的距離， 2 1H iH r 是第i個分子兩顆氫原子之間的距離，r 是平衡時氧原子與氫原子的距離，_e r_ae 是平衡時 2 1H iH r 的距離，D_e是 Morse potential 的深度，a,k_,k_r,k_rr是分子內位能的 力常數。在表(三)中列出模型 II 的所有參數。 模型參數 數值 模型參數 數值 e r 1.0Angstrom r _ae 1.6329808 Angstrom e  _109.47 _a 1 567 . 2 Angstrom 3 10  II M 630.437 kcalÅ12/mol e D _{101.9188 kcal /mol} II N 626.003kcalÅ6 /mol  k 2 / / 645606 .

328 kcal mol Angstrom

H q 0.4238 e kr 2 / /

211.4672 kcal mol Angstrom

 O q 0.8476 e krr 2 / /

111.70765kcal mol Angstrom

表(三) 非剛性水分子模型 II 的參數

II. 振動頻譜 (Vibrational spectrum)

在下文中，我們所討論的頻寬，皆是採用 Full Width at Half Maximum(FWHM) 來認定。

A. 速度頻譜 (Power spectrum)

歸一化速度相干函數 (Normalized velocity autocorrelation function) C_v

 

t 是

用來描述單一顆粒子的速度隨時間演化的相干性，其定義如下：

 

 

 

 

0

 

0 0      v v v t v t Cv    (2 - 7)

 

0  v 為 原子在選定時刻的起始速度，v_

 

t 表示在經過時間 t 後 原子的速度，  可為氧原子或氫原子。

將計算得到的氫、氧原子的歸一化速度相干函數作 Fourier cosine transform，

就可以得到氫、氧原子的速度頻譜(Power spectrum) D w ，如下：

 



_



 

 

0

cos

2

dt

wt

t

C

w

D

_{ v}



(2 - 8) 在此，D w 是歸一化的，即



  0 D()d 1。

歸一化速度相干函數也可由速度頻譜作 Inverse Fourier cosine transform 求得：

 



_



 

 

0 D w cos wt dw

t

Cv  (2 - 9)

B. 瞬間正則模分析 (Instantaneous normal mode analysis) [8]

考慮 N 個水分子在一邊長 L 的方型盒子中運動，這個系統在笛卡兒座標系的 總哈密爾頓量(Total Hamiltonian)可表示為：

 

R V r m H N k k k   _ 

_

1  3 1 2 2 1    (2 - 10) 在此，mk為第k個水分子的  原子的質量，rk  表示第k個分子的  原子在笛卡兒 座標系的位置向量，rk_  第k個分子的  原子在笛卡兒座標系的位置向量的一次微 分，V

 

R 是系統的總位能，R是 N 個水分子的系統狀態(Configuration)，總共有 9N 座標，R



r_1Hx,r_1Hy,r_1Hz, ,r_NOx,r_NOy,r_NOz



1 1 1    ，

r

k是rk  在方向上的分量， N k 1 , , ， O,H₁,H₂，  x,y,z。 我們先將座標作變換，由笛卡兒座標系(Cartesian coordinate)轉換到

Mass-weighted coordiante，兩者關係如下：    k k k

m

r

z

_

12 _{(2 - 11)}  k

z

表示在 Mass-weighted coordiante，第k個水分子的  原子的位置向量

z

k



在方 向上的分量。 在 Mass-weighted coordinate，系統的哈密爾頓量可重新寫為：

 

Z V z H N k k   _ 

_

1  3 1 2 2 1   (2 - 12) 定某時刻為時間零點，在時間零點附近的時刻 t，N 個水分子所構成的系統在 Mass-weighted coordinate 之系統狀態可表為：

 

 

 

 

 

 



z

_Hx

t

,

z

_Hy

t

,

z

_Hz

t

,

,

z

_{N x}

t

,

z

_{N y}

t

,

z

_NOz

t



Z

_{t 1 1 1 O O} 1 1 1







(2 - 13) 接著將 t 時刻的系統總位能V

 

Z_t ，對時間零點的系統狀態Z₀作泰勒展開(Taylor’s expansion)：

 



 





 





 





 



               0 2 0 0 0 0 0 2 1 Z Z z z V Z Z Z Z z V Z V Z V t Z l k t t Z k t    (2 - 14) 此主要目的是想探討水在短時間下的動力行為，而在短時間的情況下，高次項的 貢獻很微小可忽略，因此我們只取到二次項。 系統的總哈密爾頓量可寫為

 

0

 

0 0





0

  

0 0



1 3 1 2 2 1 2 1 Z Z Z D Z Z Z Z Z F Z V z H t t t N k k             _{     } 

_

   (2 - 15)

 

 

0 0 0 0 2 0 , 0 ˆ Z l k Z l k l k Z k Z k k z z V V Z D z V V Z F                               (2 - 16) 上式中，F

_{ }

Z₀ 在力學上表示一 9N 維度的力，Dˆ Z

 

0  _{為一 9N×9N 維度的矩陣，此矩}

陣名為 Hessian matrix，有關非剛性水分子模型 II 的 Hessian matrix 的各項 element 請參閱附錄 A。

matrix 的特徵值可由U

 

Z₀ 及Dˆ Z

 

₀ 求得： N U D U l k l l k k (Z ) ˆ (Z ) (Z ) 1,2, ,9 ) Z ( , 0 , 0 , 0 , 0 2     _  

_

           ， _{(2 - 17)} ) Z (₀ 

 是第th INM 的頻率， U_{ k,}(Z₀)是th INM 在 Mass-weighted

coordinate zk_方向上對應的特徵向量的分量。

因為_(Z₀)和 U_{ k,}(Z₀)會隨著系統狀態Z0 

的選擇不同而改變，我們對不同時 間的系統狀態作系綜的平均(Ensemble average)，即可求得 INM 狀態密度函數(Total INM density of state)如下：

 

_



 



   N INM Z w w N w D 9 1 0 9 1     (2 - 18) 在此，角括號( )表示對不同的系統狀態作系綜平均，且

_

DINM_{ }wdw 1。原則上，

 

w DINM 包含了每一個分子的氫原子跟氧原子的貢獻，為了區分不同原子對 INM 狀 態密度的貢獻，我們定義一個投影算符(Projection operator)p_, 如下：

 

   N k xyz k U p 1 , , 2 , ,      (2 - 19) 當  為 H 表示氫原子，  為 O 表示氧原子。由投影算符 p_, ，氫原子與氧原子對 INM 狀態密度的貢獻可表示為：

 

_



 



   N INM Z w w p N w D 9 1 0 , 9 1        (2 - 20)

III. 局部結構的分析

A. 徑向分佈函數 (Radial distribution function)

徑向分佈函數g

 

r 是描述系統中，所有粒子平均相對位置的分佈情況，其定義 如下：

 

 

N r d r g k k ab ab 



 (2 - 21)

為系統在空間中水分子的平均密度，N為系統中水分子的總個數，dabk

 

r 是當我 們將 k 分子的

a

原子置於中心，距離a原子在

r

與r dr的球殼間，

b

原子分佈的 密度。角括號表示對不同的系統狀態作系綜平均。在我們計算中，水分子的系統 組態是由我們模擬所產生出來的，可以藉由徑向分佈函數g_ab

 

r 來描述液態水的微 觀結構。在此a 跟

b

可以分別是氫原子或氧原子。我們計算氧─氧、氧─氫、氫 ─氫的徑向分佈函數，分別以g_OO r , g_OH r , gHH

 

r 表示。

B. Voronoi 分析 (Voronoi analysis)

Voronoi 分析是一種空間分割的演算法，其概念是在一個 N 維空間中，若有 k 個粒子分佈在其中，我們可以將這個空間分割成 k 個 Voronoi 單位晶胞(Voronoi unit cell)，每個 Voronoi 單位晶胞內部只有一個粒子。對空間內所有的 Voronoi 單位晶 胞，在晶胞內部的任意一點與該點距離最相近的粒子，即為該晶胞內的粒子。 水分子是由三個原子所組成，因此我們將這個概念應用在水分子的質心或水 分子的氧原子上。在本文中，我們選擇對氧原子作分析。 Voronoi 分析是在水分子跟相鄰水分子的氧原子連線上作一中垂面，將所有的 分子對都依照上述的方式作分割，然後連結這些中垂面，則每一個水分子都可以 得到一個多面體。對每一個多面體，我們定義 Asphericity  ： 2 3 36 V A    (2 - 22) A跟V分別代表多面體的總表面積跟體積 [9]。不同的，水分子的局部結構也會 不同。當多面體為球體時，為最小值 1；當多面體為正四面體，值為2.25。隨 著值逐漸變大，多面體的幾何性質會逐漸趨向非球型化(Aspherical)。 在圖(一)中，展示我們的三種水分子模型的 Asphericity 分佈曲線。我們將 Asphericity 分佈劃分為四個區間，每一區間表示為一 Voronoi group (VG)。表(四)、 表(五)分別列出非剛性水分子模型 I 跟模型 II 四個 VG 所對應的區間以及該區間所 佔的平均分子數目的百分率。由表(四)、表(五)四個 VG 的分佈情形，在我們的模

擬中，大多數水分子分佈在 VG II 和 III，而 VG I 和 IV 所佔比例很少。

圖(一) 水在 300K 的 asphericity 分佈曲線，三種模型模擬的條件皆為 1.0gcm 3&T 300K



 

 ，

solid, dash, dot-dash line 分別表示非剛性水模型 I, 非剛性水模型 II, 剛性水模型的 Asphericity 分 佈曲線。 VG I VG II VG III VG IV η below-1.46 1.46-1.72 1.72-1.98 1.98-above L  0.005948 0.505424 0.444793 0.043835 表(四) 非剛性模型 I 四種 VG 所對應的 Asphericity 區間 L



是 VG L 的 averaged number fraction， _{N N}

L / L   ， NL 是 VG L 的平均分子數，N是分子 的總數。 VG I VG II VG III VG IV η below-1.46 1.46-1.72 1.72-1.98 1.98-above L  0.006531 0.502044 0.447872 0.043553 表(五) 非剛性模型 II 四種 VG 所對應的 asphericity 區間 接著我們定義一個選擇算符(selection operator)k

 

L ：

 

       otherwise L of region the in is k molecules of y asphericit the if L k 0 VG 1 在此，L 代表四種 VG 的標記，。

選擇算符有下列的恆等式：     L k k L k N L L    





1 (2 - 23) L N 是在某一系統狀態時， VG L 區間的水分子的數目。 由選擇算符，我們可以得到每一個 VG 的氧─氧、氧─氫、氫─氫的徑向分 佈函數，由下式求得：

 

 

 

L k k k ab L ab N L r d r g 



  (2 - 24) 在此，上標 L 代表 VG 的標記， a 跟b可是氫原子或氧原子，  為系統在空間中水 的平均密度。 同樣的，我們對不同的 VG L 定義次系綜投影算符(Subensemble projection operator) L p ,：

 

 

    N k xyz k k L L U p 1 , , 2 , ,      (2 - 25)   k U _, 是 Hessian matrix 對角化之轉換矩陣。 不同 VG L 的次系綜投影算符 L p ,滿足以下的關係式， IV III II I p p p p p_,  _,  _,  _,  _,。 在某一系統狀態時，VG L 的分子數目，可由

_

_{ }

   N k k L L N 1 求得。因為每一個 VG 次系綜(VG subensemble)的平均分子數目不相同，對 INM 頻譜上的貢獻也會有差 異。為了方便比較四個 VG 次系綜之間的關係，我們對 VG L 定義一個原子的歸

一化 INM 狀態密度(Normalized  atom INM density of state) DINM(L)

 

w

 ： ：

 

_



 



   N L L L INM P Z w w N w D 9 1 , 0 ) ( 1        (2 - 26)

 

w DINM(L)  的物理意義即在 原子的 INM 頻譜D

 

w INM  中，VG L 內每一個水分子的 原子之平均貢獻。

C. 氫鍵狀態分析 (Hydrogen bond configuration analysis) Voronoi 分析在解析液態水的結構跟力學性質是很有效的方法，然而，液態水 的力學行為與水的氫鍵網絡有極高的相關性。在氫鍵狀態分析中，我們依據水分 子所鍵結的氫鍵的數目，以及氫鍵鍵結的種類來將水分子作分類。 在本文中，我們採用能量來判斷水分子間是否有氫鍵的鍵結。設定一個能量 上限E_HB，如果兩個水分子的分子間位能E 小於此能量上限



E EHB



，我們就認定 這兩個水分子之間有氫鍵的鍵結。很明顯的，兩個水分子是否有氫鍵的鍵結與能 量上限E_HB的選擇有很大的關聯 [10]，通常選擇能量介於10kj /mol到20kj /mol 之間，在本文中，我們選定能量上限為E_HB 12kj/mol。 確定兩個水分子之間有氫鍵後，我們需要知道氫鍵是分別接在兩個水分子的 哪個原子上。兩個不同水分子間的氫─氧對有四種組合，我們選擇距離最短的氫─ 氧對作為此二分子氫鍵鍵結的位置，則每一氫鍵可表示為(O─HO)，參考圖 (二)。在氫鍵的位置確定後，我們將提供氫鍵中氫原子的水分子稱為 donor，接受 此氫原子的水分子稱為 acceptor。 圖(二) 兩個水分子間判定氫鍵的示意圖 依據上述的判定方法，所有的水分子都可以依其所鍵結的氫鍵分類，我們利 用DnAm這種標記法來描述一個水分子的氫鍵狀態。例如，D2A1代表此水分子與

其相鄰水分子共有三個氫鍵，兩個氫鍵鍵結對此水分子是 donor H-bond，一個鍵結 對此水分子為 acceptor H-bond。在我們的模擬系統中，可以找到很多種不同的氫 鍵狀態(H-bond configuration)，在本文中，挑選出十組可能的氫鍵狀態來探討氫鍵 跟水分子振動的關係。我們以 Config. S 來表示這十組氫鍵狀態，S 是不同氫鍵狀 態的標記。非剛性模型 I 跟非剛性模型 II 的十組 Config. S 及其平均分子數目的百 分率，參考表(六)、表(七)。

Config. S D2A3 D3A2 D2A2 D1A2 D2A1

S

 0.0506 0.0160 0.5282 0.1323 0.1878

Config. S D0A2 D1A1 D2A0 D0A1 D1A0

S

 0.0071 0.0658 0.0030 0.0052 0.0011

表(六) 非剛性模型 I 十種氫鍵狀態及其分子數目的平均百分率

S

 是 Config. S 的 averaged number fraction。

Config. S D2A3 D3A2 D2A2 D1A2 D2A1

S

 0.0514 0.0161 0.5264 0.1313 0.1877

Config. S D0A2 D1A1 D2A0 D0A1 D1A0

S

 0.0074 0.0670 0.0031 0.0055 0.0012

表(七) 非剛性模型 II 十種氫鍵狀態及其分子數目的平均百分率

與 Voronoi 分析相似，我們對每一種 Config. S 定義原子的歸一化 INM 狀態

密度(Normalized  atom INM density of state)DINM(S)

 

w

 ：

 

_



 



   N S S S INM _{w w Z P} N w D 9 1 , 0 ) ( 1        (2 - 27) S N 是在某一系統狀態時，擁有氫鍵狀態 Config. S 的水分子數目，

_

_{ }

   N k k S S N 1 。 S

p_{ ,} 是每一種 Config. S 的次系綜投影算符(Subensemble projection operator)，

 

 

    N k xyz k k S S U p 1 , , 2 , ,      ， U k,是 Hessian matrix 對角化之轉換矩陣，k

 

S 是

 

       otherwise S Config as same the is k molecules of ion Configurat HB the if S k 0 . 1

 

S k  滿足以下兩個恆等式：

 

 

S k k S k N S S    





1 (2 - 28) 以上我們介紹了 Voronoi 分析與氫鍵狀態分析兩種方法，用來討論液態水的 局部結構，接下來我們要探討這兩種方法分析的結果之間的相互關係。表(八)與表 (九)為模型 I 及模型 II 對四個 VG 的分子個別計算的氫鍵狀態分佈，從這個分佈我 們可以知道每一種 VG 內水分子氫鍵可能的鍵結情況。圖(三)及圖(四)是模型 I 及 模型 II 四個 VG 的氫鍵狀態分佈圖。 Voronoi group I

Config. S D2A3 D3A2 D2A2 D1A2 D2A1

S

 0.0048 0.0063 0.1093 0.1663 0.2157

Config. S D0A2 D1A1 D2A0 D0A1 D1A0

S

 0.0426 0.3387 0.0103 0.0859 0.0136

Voronoi group II

Config. S D2A3 D3A2 D2A2 D1A2 D2A1

S

 0.0369 0.0175 0.4309 0.1606 0.2281

Config. S D0A2 D1A1 D2A0 D0A1 D1A0

S

 0.0106 0.0984 0.0041 0.0079 0.0018

Voronoi group III

Config. S D2A3 D3A2 D2A2 D1A2 D2A1

S

 0.0669 0.0150 0.6340 0.1031 0.1457

Config. S D0A2 D1A1 D2A0 D0A1 D1A0

S

 0.0029 0.0268 0.0018 0.0010 0.0003

Voronoi group IV

Config. S D2A3 D3A2 D2A2 D1A2 D2A1

S

 0.0612 0.0095 0.7261 0.0733 0.1113

Config. S D0A2 D1A1 D2A0 D0A1 D1A0

S

 0.0015 0.0144 0.0012 0.0004 0.0002

圖(三) 非剛性水模型 I 四種 VG 的氫鍵狀態分佈圖 從圖(三)的分佈我們可以很清楚的知道，模型 I 在 VG IV 的分子中，D2A2是 出現機率最高的，超過七成的水分子是這種氫鍵狀態，而且 VG IV 的分子大部分 鍵結的氫鍵數目都在兩個以上。VG I 的情況是另一種極端，與 VG IV 相較有很大 的不一樣。VG I 分子的分佈是比較多樣的，D1A1雖是出現機會最高的，但其他幾 種氫鍵狀態所占的比例也有一定分量，且相差並不懸殊。如果單看各個 VG 的 2 A 2 D 所佔比例，我們發現隨著 VG 的 Asphericity 增加，找到一個擁有D2A2的水 分子的機會也逐漸增加。但是對D1A1狀態，找到水分子的機率是隨著 Asphericity 增加而減少。 從圖(四)觀察到模型 II 不同 VG 的氫鍵分佈情形，與模型 I 所計算出的分佈是 很相似的。D1A1與D2A2分別是 VG I 與 VG IV 中出現比例最高的氫鍵類別，且 隨 Asphericity 增加，D2A2在各個 VG 中的比例逐漸變大。

Voronoi group I

Config. S D2A3 D3A2 D2A2 D1A2 D2A1

S

 0.0045 0.0050 0.1069 0.1586 0.2076

Config. S D0A2 D1A1 D2A0 D0A1 D1A0

S

 0.0418 0.3479 0.0121 0.0922 0.0160

Voronoi group II

Config. S D2A3 D3A2 D2A2 D1A2 D2A1

S

 0.0370 0.0176 0.4266 0.1597 0.2292

Config. S D0A2 D1A1 D2A0 D0A1 D1A0

S

 0.0111 0.1009 0.0042 0.0085 0.0018

Voronoi group III

Config. S D2A3 D3A2 D2A2 D1A2 D2A1

S

 0.0682 0.0151 0.6336 0.1023 0.1453

Config. S D0A2 D1A1 D2A0 D0A1 D1A0

S

 0.0030 0.0267 0.0018 0.0010 0.0003

Voronoi group IV

Config. S D2A3 D3A2 D2A2 D1A2 D2A1

S

 0.0632 0.0092 0.7232 0.0746 0.1099

Config. S D0A2 D1A1 D2A0 D0A1 D1A0

S

 0.0014 0.0155 0.0012 0.0006 0.0001

表(九) 模型 II 四種 voronoi group 的分子各種 Config. S 及其水分子的平均百分率

第三章 分子動力模擬

分子動力學模擬是很廣泛地應用在研究液體的力學與結構性質。在本文中， 我們模擬在溫度T 300K時，256 個水分子在一正方形盒子中的運動，盒子內水分 子的密度為 3 / 0 . 1 g cm   ，且這 256 個水分子的運動滿足週期性的邊界條件，參考 圖(五)。 圖(五) 電腦模擬所採用的週期性邊界條件示意圖。 由於庫倫作用屬於長距離的作用力，每個水分子除了會跟模擬盒子內部的水 分子之間有庫倫交互作用，在盒子外的其它水分子也會跟這個水分子有作用力， 且這個作用力是不可被忽略的。在本文中我們使用 Ewald summation method 計算 盒子外分子對每一水分子所受的庫倫力總合。對兩個氧原子間的 Lennard-Jones

在截斷的地方維持連續。在我們的模擬計算中，所採用的演算法是 leap-frog algorithm[11]，每一步的時間間隔為 0.2 飛秒(femtosecond)。 在開始進行模擬前，我們讓水分子以面心立方晶格結構(face-centered cubic lattice structure)分布在方型盒子內，且有一隨機分佈的速度，這個速度分佈必須滿 足 Maxwell-Boltzmann distribution。為了確定系統已達平衡狀態，我們讓系統先經 過 1,500,000 步的模擬過程，才開始進行我們的模擬實驗，並收取數據，每五個時 間間隔收取一組 configuration，總共收取 1,000,000 組。

第四章 結果分析

在本章中，我們將非剛性模型 II 所計算的結果與剛性模型[13]以及非剛性模 型 I 的結果作比較。

I.

局部結構與液態水的 INM 頻譜

A. 徑向分佈函數 在圖(六)中，展示著三種模型的水分子在溫度 300K 時所計算的徑向分佈函 數。首先，我們先比較兩種非剛性水模型的徑向分佈函數，由圖中可以看到兩種 模型所計算的g_OO r , g_OH r , gHH

 

r 是非常相近的，也就是兩種水分子模型所模擬出 的液態水，結構極為相似。這表示我們在模型 II 中所增加的鍵長與鍵角變化的交 叉項，不會對液態水的結構有很明顯的影響。

圖(六) 水在 1.0gcm3&T 300K的徑向分佈函數(a) g_{OO }r , (b) g_OH r , (c) gHH r ，solid, dash, dot-dash line 分別表示模型 II, 模型 I, 剛性模型的徑向分佈函數。

接著是剛性模型與非剛性模型的比較，對三種徑向分佈函數從圖(六)可看出兩 者所計算的結果，在徑向上的變化趨勢是相似的，但非剛性模型的 first maximum 位置比較接近中心原子且強度比剛性模型稍高一點。first minimum 的位置非剛性 模型仍然比較靠近中心原子，但強度則是剛性模型比較高。 B. 速度頻譜 非剛性模型 II 氫原子跟氧原子的歸一化速度相干函數展示在圖(七)。在圖中， 我們可以很明確的看到氫原子的歸一化速度相干函數 CvH

 

t 有高頻率振盪行為， 而相同的振盪行為在氧原子的歸一化速度相干函數 CvO

 

t 中也有出現但強度相對 是很微弱的。

圖(七) 非剛性水模型 II 的歸一化速度相干函數，solid, dash line 分別表示模型 II 的氫、氧原子的 歸一化速度相干函數。

圖(八)是模型 I 和模型 II 的歸一化速度相干函經過 Fourier cosine transform 所 得到的速度頻譜。我們所計算的模型 II 的速度頻譜顯示液態水有三個頻帶

(frequency band)。這三個頻帶的頻率由低到高，分別是由水分子之間的平移旋轉運 動( 1 1 1200 ~ 0cm cm )、分子內部的彎曲振動( 1 1 2100 ~ 1600cm cm )以及OH伸縮 運動(_3000cm1 _~3900cm1_{)所產生的。氫原子對高頻的彎曲運動和伸縮運動兩個頻} 帶有很顯著貢獻，這符合我們在C_vH

 

t 中所觀察到的高頻率振盪行為，而氧原子主 要的貢獻集中於低頻的水分子平移和旋轉運動，高頻部分的貢獻只佔速度頻譜中 很小的一部分。

圖(八)模型 I 和模型 II 的速度頻譜 (a) Oxygen, (b) Hydrogen。

在兩種非剛性模型的比較上，我們先討論氫原子的速度頻譜，在頻率_1200cm1

以下的部分，因為主要是受水分子間的平移旋轉運動所影響，而兩種模型分子間

位能的形式相同，因此兩者在這個部分所計算的結果大致吻合。在頻率_{1200cm 以}1

差異。由彎曲運動所產生的頻帶，模型 II 的分佈大概以頻率_{1880cm 為中心，比}1 模型 I 的 1 1650cm 略為往高頻區域移動。強度則是模型 I



D_H

 

 0.00115



較模型 II



DH

 

 0.00095



為高。頻寬的部份，模型 I 大約在 1 1 1720 ~ 1620cm cm ，模型 II 則是在1810cm1 ~1940cm1。有關伸OH縮運動的頻帶，兩種模型分佈的中心位 置差不多都在_{3450cm 附近，但模型 I 的頻寬較窄，約為}1 _3340cm1 _~3540cm1_， 模型 II 的頻寬較寬，大約 1 1 3680 ~ 3260cm cm 。但在模型 II 的分佈中，可以觀察 到有一個 shoulder 的情形，位置大約在_{3550cm 。其原因與水分子內兩對 OH 原}1 子間的伸縮運動是否對稱有關，當兩對 OH 原子的距離是同時伸長或縮短(對稱)， OH 的振盪頻率會比 OH 一伸一縮(非對稱)的情況低[12]。由圖(八)可知，模型 II 所模擬的水分子，兩對 OH 原子的振盪是以對稱的情況較多。 C. INM 頻譜

INM 分析中，hessian matrix 的特徵值並不一定為正數。當特徵值為正數時， 所求得的 INM 我們稱之為 stable mode；當特徵值為負數時，我們會得到虛數的頻 率，這種虛數頻率的 INM 我們稱之為 unstable mode。在計算 INM 的狀態密度時， 我們將虛數的頻率所對應的部分劃於頻率座標的負端。 圖(九)為我們利用 INM 分析，所求得的模型 II 的 INM 狀態密度。在實數頻率 (real frequency)的部分，狀態密度有三個頻帶，頻率 1 1200 ~ 0 cm 的部份對應水分

子的平移旋轉運動，在這個區域的 Total INM 狀態密度有兩個 peak，位置在_50cm1

及 1 670cm 。藉由氫原子跟氧原子在 Total INM 狀態密度上的貢獻，我們可以了解 在_{50cm 附近的 peak，主要是由氧原子的運動所產生的，而}1 _{670cm 的 peak 則是}1 與氫原子的運動有關。頻率在 1 1 2200 ~ 1400cm cm 和 1 1 4600 ~ 2300cm cm 兩個頻 帶，分別對應水分子內部的彎曲以及伸縮運動，在這兩個頻帶氫原子的貢獻較顯 著。在虛數頻率(imagine frequency)的部分，狀態密度也有兩個 peak，位置在

1

50 

圖(九) 非剛性模型 II 的 INM 狀態密度，solid line 是 Total INM 狀態密度，dash, dot-das line 分別表 示 INM 狀態密投影在氫原子跟氧原子上的分量。 在圖(十)中，我們比較模型 I 跟模型 II 的 INM 狀態密度。在虛數頻率的部份 以及實數頻率的平移旋轉運動兩部分，兩種非剛性模型的 INM 狀態密度沒有明顯 的差異。但在彎曲運動的頻帶，模型 II 的 peak 位置在 1 1850cm ，相較於模型 I 的 peak 位置在_{1650cm ，模型 II 是向高頻區域移動。在頻寬的部分，模型 II 大約在}1 1 1 1730 ~ 1590cm cm ，而模型 I 約為 1 1 2020 ~ 1720cm cm ，模型 II 的頻寬約是模 型 I 的兩倍。但模型 II 的峰值大約為模型 I 的一半。在伸縮運動的頻帶，兩種模型 的 peak 位置相差不大，大約都在 1 3500cm 附近，但模型 II 的頻寬 ( 1 1 4000 ~ 3000cm cm )，比模型 I( 1 1 3780 ~ 3160cm cm )更寬，且峰值方面，模型 II 比模型 I 低。模型 II 所觀察到的頻譜變寬的結果，可能是因為鍵長與鍵角變化的 交叉項(stretch-bend coupling term)的增加，讓水分子內的彎曲及伸縮運動變得更為 複雜所導致的。

圖(十) 非剛性模型 I 和非剛性模型 II 的 INM 頻譜比較，黑線表示模型 II 的 INM 狀態密度，紅線 表示模型 I 的 INM 狀態密度，Solid line 是 Total INM 狀態密度，dash, dot-dash line 分別表示 INM 狀態密度投影在氫原子跟氧原子上的分量。

我們由速度頻譜及瞬間正則模分析皆可求得液態水的狀態密度。圖(十一)為模 型 II 使用兩種方法所求得的狀態密度的比較。由於速度頻譜的狀態密度沒有 unstable mode，我們取 INM 頻譜中實數頻率的部分作歸一化，比較兩種計算方法

的區別。在彎曲跟OH伸縮運動的頻帶，INM 頻譜及速度頻譜的 peak 位置相近，

可是 INM 頻譜的頻寬較寬。但是在平移旋轉運動的頻帶，不論是氧原子或氫原子 的狀態密度，我們都可以明顯看出兩種方法計算的頻譜有明顯的差異。造成這種 差異的原因，可能是因為瞬間正則模分析使用 harmonic approximation，忽略高次 項的貢獻所導致。

圖(十一) 模型 II 分別對氫、氧原子作投影的速度頻譜跟 INM 頻譜 (a)氧原子, (b)氫原子，solid,dash line 分別表示 real INM 頻譜跟速度頻譜，real INM 頻譜代表 INM 狀態密度中實數頻率的部分。

II.

Subsemble by Voronoi polyhedral analysis

A. 徑向分佈函數

圖(十二)中展示四個 VG 的徑向分佈函數g_OOL

 

r 、g_OHL

 

r 、gHHL

 

r 。在g

 

r

L

OO 中，

水分子的結構可以分類成兩種極端的情形：一種是沿著徑向方向上，粒子分佈具 有一個相似的有序結構，如 VG II, III ,IV 所展示的。另一情況是，VG I 水分子的

局部結構，在徑向方向上，gI

 

r

OO 的高度會上下起伏，但變化的幅度不大。也就是

說 VG I 水分子的分佈結構比較像是無序的分佈。而且g_OOI

 

r 的 First maximum 及

Second maximum 的高度很接近，與 VG II, III ,IV 所觀察到的情形很明顯的不同。 我們從圖(十二)也可以觀察到，三種徑向分佈函數的 Frist maximum，隨著

圖(十二) 模型 II 四個 VG 的徑向分佈函數 (a) _gL _{ r}

OO , (b) g  r L

OH , (c) g  r L

HH ， solid, dot, dash,

dot-dash line 分別對應 VG I, II, III, IV。

接下來我們將比較非剛性模型 I 跟非剛性模型 II 的gL

_{ }

r OO 、g

 

r L OH 、g

 

r L HH ，分 別展示在圖(十三)、圖(十四)、圖(十五)中。由這三張圖我們可以看到兩種模型所 計算的不同 VG 的徑向分佈函數幾乎是一致的，結合前面所計算的整體水的徑向 分佈函數及 Asphericity 分佈圖，我們可以知道兩種模型模擬的液態水的結構幾乎 是一樣。這再一次說明模型 II 中增加的鍵長與鍵角變化的交叉項對液態水結構的 影響確實是很微弱的。

圖(十三) 非剛性模型 I 及非剛性模型 II 的氧-氧徑向分佈函數 (a) gI

_{ }

r OO , (b) g

 

r II OO , (c)

 

r gIII OO , (d) g

 

r IV

OO ，solid, dash line 分別表示模型 II 跟模型 I。

圖(十四) 非剛性模型 I 及非剛性模型 II 的氧-氫徑向分佈函數 (a) gI

_{ }

r OH , (b) g

 

r II OH , (c)

 

r gIII _{, (d)}

_{ }

r

圖(十五) 非剛性模型 I 及非剛性模型 II 的氫-氫徑向分佈函數 (a) g_HHI

 

r , (b) g_HHII

 

r , (c)

 

r

g_HHIII , (d) g_HHIV

 

r ，solid, dash line 分別表示模型 II 跟模型 I。

B. INM 頻譜 圖(十六)是模型 II 對氫、氧原子在不同 VG 貢獻的歸一化 INM 狀態密度。由 於氧原子的 INM 狀態密度主要集中於在平移旋轉運動的頻帶，彎曲跟伸縮運動的 貢獻佔整體比例相對很低，為了便於觀察，我們將氧原子的 INM 狀態密度以頻率 1 00 12 cm 為分水嶺，分成兩部分來看，分別置於圖(十六)中的(b)、(c)。  首先，我們從氫原子的貢獻開始討論，在虛數頻率的部分，隨 VG 的 Asphericity 增加，unstable mode 所佔的比例會逐漸減少。實數頻率的部分，我們可以觀察到 當 Asphericity 增大，平移旋轉運動以及彎曲運動的頻帶有藍位移(blue shift)的現 象，但伸縮運動的頻帶則是出現紅位移(red shift)的情形。在平移旋轉運動的頻帶，

不單只有 peak 位置的移動，峰值也有變化，在頻率_{70cm 的 peak 其峰值會隨}1

圖(十六) 模型 II 不同 VG 的歸一化 INM 狀態密度 (a)氫原子, (b)氧原子, (c)氧原子，solid, dot, dash, dot-dash line 分別對應 VG I, II, III, IV。

氧原子在不同 VG 的 INM 狀態密度所觀察到的變化，大致上與氫原子對 INM 狀態密度的貢獻所看到的隨 Asphericity 變化的趨勢很相似，只有平移旋轉運動的 頻帶與氫原子所觀察到的趨勢不一樣。氧原子平移旋轉運動的頻帶在 VG I 只有一

個 peak，位置大約是_{70cm ，但隨著 Asphericity 增加，在 300}1 _{cm 附近 INM 狀態}1

密度逐漸增加，在 VG III 時，我們已經可以觀察到有一個 shoulder 的產生。在氧 原子平移旋轉運動的頻帶依然能看到頻譜藍位移的現象。

(十七)中不論是氫原子或氧原子，我們都可以觀察到當 Aspheticity 增加，unstable mode 的比例減少、平移旋轉運動和彎曲運動頻帶的藍位移、以及伸縮運動頻帶的 紅位移，這個結果與模型 II 所計算的是一致的。

圖(十七) 模型 I 不同 VG 的歸一化 INM 狀態密度 (a)氫原子, (b)氧原子, (c)氧原子，solid, dot, dash, dot-dash line 分別對應 VG I, II, III, IV。

由兩種模型的結果可知，我們使用幾何方法的 Voronoi 分析，依照水分子不同 局部結構將水分子作分類，來探討結構上的差異與水分子動力學性質的關係是有 效的。雖然兩種模型所模擬的水分子在力學性質上有一些不同的地方，但兩種模 型都可以觀察到不同局部結構對 INM 頻譜的影響，且局部結構的改變與 INM 頻譜

變化之間的關係，對兩種非剛性水模型皆是相同的。由這個對應關係，當我們知 道一個水分子周遭其他水分子的分佈情形，我們就可以大概了解這個水分子 INM 頻譜改變的趨勢。

III.

Subsemble by Hydrogen bond configuration analysis

A. 徑向分佈函數 非剛性模型 I 與非剛性模型 II 三種徑向分佈函數gS

 

r OO 、g

 

r S OH 、g

 

r S HH ，分 別展示於圖(十八)、圖(十九)、圖(二十)。由這三張圖兩種模型結果的比較，我們 可以知道這兩種水分子模型所模擬的液態水的結構是很相近的，與我們在 Voronoi 分析所得到的結果是一致的。 圖(十八) 非剛性模型 I 及非剛性模型 II 的氧-氧徑向分佈函數 (a) g_OOD2A3

 

r , (b) g_OOD3A2

 

r , (c)

 

r g_OOD2A2 , (d) g_OOD1A2

 

r , (e) g_OOD2A1

 

r , (f) g_OOD0A2

 

r , (g) g_OOD1A1

 

r , (h) g_OOD2A0

 

r , (i) g_OOD0A1

 

r , (j) g_OOD1A0

 

r ，solid, dash line 分別表示模型 I 跟模型 II。

在模型 I 與模型 II 的gS

 

r

OO 中，水分子的徑向分佈情形可大致分為兩類：第一

種，水分子的 total H-bond 數目在三個以上，其結構分佈是比較有序的；第二種，

水分子的 total H-bond 數目小於三個，其分佈沒有明顯的規律性。gS

 

r

OO 的 first

maximum 高度受水分子 total H-bond 的數目影響，水分子的 total H-bond 數目越多，

first maximum 的高度越高。在gS

 

r

OH ，水分子的徑向分佈函數與水分子 acceptor

H-bond 數目有關，first maximum 的高度隨水分子 acceptor H-bond 數目增加而升

高。gHHS

 

r 的部份，當水分子的 total H-bond 數目越多，徑向分佈函數的 first

maximum 會提高，然後 first minimum 會降低。

圖(十九) 非剛性模型 I 及非剛性模型 II 的氧-氫徑向分佈函數 (a) g_OHD2A3

 

r , (b) g_OHD3A2

 

r , (c)

 

r

g_OHD2A2 , (d) g_OHD1A2

 

r , (e) g_OHD2A1

 

r , (f) g_OHD0A2

 

r , (g) g_OHD1A1

 

r , (h) g_OHD2A0

 

r , (i) g_OHD0A1

 

r , (j) g_OHD1A0

 

r ，solid, dash line 分別表示模型 I 跟模型 II。

圖(二十) 非剛性模型 I 及非剛性模型 II 的氫-氫徑向分佈函數 (a) g_HHD2A3

 

r , (b) g_HHD3A2

 

r , (c)

 

r

g_HHD2A2 , (d) g_HHD1A2

 

r , (e) g_HHD2A1

 

r , (f) g_HHD0A2

 

r , (g) g_HHD1A1

 

r , (h) g_HHD2A0

 

r , (i) g_HHD0A1

 

r , (j) g_HHD1A0

 

r ，solid, dash line 分別表示模型 I 跟模型 II。

B. INM 頻譜

圖(二十一)、圖(二十二)中，分別展示模型 I 與模型 II 氫原子各種 Config. S 的 歸一化 INM 狀態密度。在圖(二十一)中，unstable mode 的比例隨著 donor H-bond 的數目增加而減少，在 donor H-bond 數目相同時，水分子的 Acceptor H-bond 數目 越多，unstable mode 的比例越少。在平移旋轉運動、彎曲運動及伸縮運動的頻帶 中，donor H-bond 的數量都是影響氫原子 INM 頻譜的主要因素，其次是 aceeptor H-bond 的數目，但表現在氫原子 INM 頻譜上的行為卻是不同的。Donor H-bond 和 acceptor H-bond 數量的增加在平移旋轉運動、彎曲運動頻帶上，是讓 INM 頻譜

向高頻移動，但在伸縮運動的頻帶所看到的變化是 INM 頻譜會往低頻移動。模型 II 所計算的結果與模型 I 是相同的，不過 acceptor H-bond 數目對彎曲運動及伸縮 運動的頻帶所產生的影響，模型 II 是較模型 I 明顯。

圖(二十一) 非剛性模型 I 氫原子各種 Config. S 的歸一化 INM 狀態密度 (a)three donor H-bond, (b)two donor H-bond, (c)one donor H-bond, (d)no donor H-bond。

圖(二十二) 非剛性模型 II 氫原子各種 Config. S 的歸一化 INM 狀態密度 (a)three donor H-bond, (b)two donor H-bond, (c)one donor H-bond, (d)no donor H-bond。

氧原子的 INM 頻譜因為彎曲運動及伸縮運動的貢獻較少，我們將整個 INM 頻 譜分成三個部分個別討論，而且我們除了討論 donor H-bond 數目對 INM 頻譜的影 響，也會分析 acceptor H-bond 數目以及 Total H-bond 數目與 INM 頻譜之間的關係。 圖(二十三)、圖(二十四)、圖(二十五)分別展示兩種非剛性水模型依照不同 donor

轉運動頻帶(translational and rotational spectrum)。從圖(二十三)、圖(二十四) 的結 果，我們可以發現氧原子 INM 頻譜的變化與水分子 donor H-bond 及 acceptor H-bond 的數量沒有明確的關係。若是由水分子 Total H-bond 數目來看，圖(二十五)顯示隨 著 Total H-bond 數目增加，模型 I 與模型 II 都有 unstable mode 減少且平移旋轉運 動的頻率帶向高頻偏移的趨勢。

圖(二十三) 氧原子各種 Config. S 的 imahinary prat 及平移旋轉運動頻譜 (a)three donor H-bond, (b)two donor H-bond, (c)one donor H-bond, (d)no donor H-bond。

圖(二十四) 氧原子各種 Config. S 的 imahinary prat 及平移旋轉運動頻譜 (a)three acceptor H-bond, (b)two acceptor H-bond, (c)one acceptor H-bond, (d)no acceptor H-bond。

圖(二十五) 氧原子各種 Config. S 的 imahinary prat 及平移旋轉運動頻譜 (a)five total H-bond, (b)four total H-bond, (c)three total H-bond, (d)two total H-bond, (e)one total H-bond。

圖(二十六)、圖(二十七)、圖(二十八) 分別是兩種非剛性水分子模型依照不同 donor H-bond、acceptor H-bond、Total H-bond 數目數目分類的氧原子彎曲運動頻帶 (bending spectrum)。圖(二十七)的結果顯示水分子 acceptor H-bond 的數量不是造成 INM 頻譜偏移的主要因素。比較圖(二十六)、圖(二十八)的頻譜，我們發現模型 I 的 INM 頻譜變化與水分子的 donor H-bond 數量有顯著的關係，當 donor H-bond 的

數目變多，INM 頻譜會有藍位移的現象，且峰值會逐漸下降；但在模型 II 的 INM 頻譜，隨著水分子 total H-bond 及 donor H-bond 數目的不同，INM 頻譜都可觀察到 藍位移現象，且峰值不會隨 total H-bond 及 donor H-bond 數目不同而變化。

圖(二十六) 氧原子各種 Config. S 的彎曲運動頻譜 (a)three donor H-bond, (b)two donor H-bond, (c)one donor H-bond, (d)no donor H-bond。

圖(二十七) 氧原子各種 Config. S 的彎曲運動頻譜 (a)three acceptor H-bond, (b)two acceptor H-bond, (c)one acceptor H-bond, (d)no acceptor H-bond。

圖(二十八) 氧原子各種 Config. S 的彎曲運動頻譜 (a)five total H-bond, (b)four total H-bond, (c)three total H-bond, (d)two total H-bond, (e)one total H-bond。

最後的伸縮運動頻帶(stretching spectrum)我們一樣將 INM 頻譜依照不同數量 的 donor H-bond、acceptor H-bond、Total H-bond 作分類，觀察頻譜的變化，結果 分別展示在圖(二十九)、圖(三十)、圖(三十一)。由圖(二十九)、圖(三十一)的結果， 我們發現當水分子 donor H-bond 及 total H-bond 的數量改變，模型 I 與模型 II 都可 觀察到 INM 頻譜往低頻移動的現象。從圖(三十)中，模型 I 與模型 II 的 INM 頻譜

變化與水分子 acceptor H-bond 的數目沒有明確的對應關係。

圖(二十九) 氧原子各種 Config. S 的伸縮運動頻帶 (a)three donor H-bond, (b)two donor H-bond, (c)one donor H-bond, (d)no donor H-bond。

圖(三十) 氧原子各種 Config. S 的伸縮運動頻帶 (a)three acceptor H-bond, (b)two acceptor H-bond, (c)one acceptor H-bond, (d)no acceptor H-bond。

圖(三十一) 氧原子各種 Config. S 的伸縮運動頻帶 (a)five total H-bond, (b)four total H-bond, (c)three total H-bond, (d)two total H-bond, (e)one total H-bond。

第五章 結論

雖然非剛性模型 I 與非剛性模型 II 採用不同的分子內位能描述分子內部原子 間的運動，而且我們使用速度頻譜或瞬間正則模分析所計算的頻譜，也都顯示兩 種非剛性水模型的分子內部的微觀力學行為是不一樣的，但在結構性質上，兩者 表現的行為是很相似的。 在液態水局部結構分析的部分，非剛性模型 I 與非剛性模型 II 隨水分子的 Asphericity 增加，所觀察到的 INM 頻譜變化是一致的，在彎曲運動頻帶，頻譜是 往高頻偏移，而伸縮運動頻帶則是往低頻區移動。若是以水分子的氫鍵狀態來看， 兩種非剛性模型的氫原子 INM 頻譜都與水分子的 donor H-bond 數目有關，且不同 donor H-bond 數目在 INM 頻譜上的變化趨勢是相同的；但氧原子的 INM 頻譜，在 不同的頻帶造成變化 INM 頻譜變化的因素，非剛性模型 I 與非剛性模型 II 並不完 全是一樣的。