管星悅，李文飛

(南京大學 物理學院，江蘇 南京 210093)

體相水分子系統(tǒng)是人們最為熟悉的分子系統(tǒng)[1].然而，由于分子間的氫鍵網絡結構，水分子系統(tǒng)具有非常豐富的結構與動力學行為，并表現(xiàn)出諸多特殊性質，如高熱容、結冰膨脹以及液相密度的非單調溫度依賴性等.另外，水分子也是生命體在分子、細胞、組織等多個層次生命運動的關鍵參與者，因此水分子系統(tǒng)不僅是物理學專業(yè)學生學習分子運動與熱力學知識的重要實例體系，也是物理、化學、以及生命科學領域的重要研究對象.

體相水的結構與動力學性質一直是人們探討的重要問題.2006年Laage和Hynes[2]發(fā)現(xiàn)水分子存在大角度跳躍轉動機制.另外，由于水分子的固有偶極矩，外加靜電場可以顯著改變水分子系統(tǒng)的結構、動力學與熱力學性質.例如，在外加穩(wěn)恒電場強度大于2.5 V·nm-1時，可以觀察到液固相變，即電致結冰[3].Zong等人[4]發(fā)現(xiàn)外加靜電場能夠導致體相水分子系統(tǒng)的黏性產生各向異性.

分子動力學模擬通過求解系統(tǒng)中基本單元的運動方程，從而得到坐標、速度和能量等信息的時間演化，被廣泛應用于分子系統(tǒng)結構、動力學與熱力學性質的研究.本文基于分子動力學模擬，研究了靜電場環(huán)境下水分子的大角度跳躍轉動運動以及液固相變等性質.結果發(fā)現(xiàn)在小電場區(qū)域，水分子的大角度跳躍轉動運動頻次呈現(xiàn)出隨靜電場強度先增加后減小的非單調依賴性.在大電場區(qū)域，水分子系統(tǒng)發(fā)生液固相變，即電致結冰.通過對冰相結構的進一步分析，發(fā)現(xiàn)電致結冰導致的冰相存在依賴于電場方向的不同冰相的層狀混合堆疊結構.本文中的結果加深了人們對水分子系統(tǒng)物理特性的認識，并提供了關于靜電場環(huán)境下分子擴散與相變等物理現(xiàn)象的新理解.另外，本工作中的分子動力學模擬不僅計算量小、易于實踐操作，而且結果具有創(chuàng)新性，因此可作為教學實例培養(yǎng)本科生的實踐創(chuàng)新能力和科研興趣.

1 分子動力學模擬設計

本文采用了SPC/E水模型.其中水分子的鍵長和鍵角如圖1所示.非共價鍵相互作用包括原子間的靜電勢和Lennard-Jones勢,即

(1)

上式中，rij為原子i和j之間的距離，rOO為氧原子之間的距離.qi和qj為原子電荷，其中氧原子和氫原子的電荷分別為-0.847 6 e和0.423 8 e.ε0為真空中的介電常數(shù).A=629.4×103kcal·?12·mol-1和B=625.5 kcal·?6·mol-1為Lennard-Jones參數(shù).以上SPC/E水模型中的原子電荷與Lennard-Jones參數(shù)通過擬合實驗測得的液相水密度、擴散常數(shù)等靜態(tài)和輸運性質得到.由式(1)給出的能量函數(shù)可以計算出每個原子的受力，進而通過求解運動方程得到所有原子的坐標、速度隨時間的演化.分子動力學模擬系統(tǒng)由510個水分子組成，初始結構如圖1所示.分子動力學模擬采用通用軟件Gromacs在等溫等壓系綜下進行，并采用了周期性邊界條件.其中，壓強設為1個標準大氣壓，溫度設為260 K，高于SPC/E水模型對應的熔點.每次分子動力學模擬時長為5 ns，每個電場強度下獨立重復32次模擬進行統(tǒng)計，得到平均值和標準偏差.

水分子氫鍵網絡示意圖(上)與SPC/E水分子模型幾何結構示意圖(下) 分子動力學模擬初始結構

2 靜電場下水分子的大角度跳躍轉動與擴散

2.1 大角度跳躍轉動的電場強度依賴性

圖2給出了水分子大角度跳躍轉動示意圖與代表性軌跡.根據(jù)文獻[2]的研究，水分子的轉動運動呈現(xiàn)出大角度跳躍轉動，即某水分子從斷開一個氫鍵到形成新氫鍵的過程中(圖2)，其取向角度存在大幅度的瞬時變化.大角度跳躍轉動可由水分子角平分線的取向角度θ來定量表征(其中θ是水分子1的角平分線矢量在2-1-3平面內轉過的角度，如圖2所示).可以看到在皮秒時間尺度內，水分子取向角度θ從-1.0 rad跳變到1.0 rad，表現(xiàn)出大角度跳躍轉動.基于分子動力學模擬軌跡，可以統(tǒng)計出大角度跳躍轉動發(fā)生的頻次.其中，在判定一對水分子之間是否形成氫鍵時，根據(jù)兩個水分子的氧原子間距以及氧-氫-氧原子形成的夾角數(shù)值，并采用Wernet-Nilson 判據(jù)來計算.

圖2 水分子大角度轉動事件示意圖(上)與典型軌跡(下)

圖3給出了大角度跳躍轉動頻次以及平均每個水分子形成的氫鍵數(shù)目隨電場強度的變化.結果顯示，大角度跳躍轉動頻次和氫鍵數(shù)目隨電場強度呈非單調變化.隨著電場強度的增加，大角度跳躍轉動頻次先增加，后減小，在0.4 V·nm-1附近時達到最大.同時，水分子平均氫鍵數(shù)呈現(xiàn)出相反的趨勢，即隨著電場強度先減小，后增加，且平均氫鍵數(shù)較小的區(qū)域轉動頻次較大.值得注意的是，在文獻[5]中，作者基于反應平衡與玻耳茲曼分布所建立的理論模型計算得到的結論認為水分子平均氫鍵數(shù)隨電場強度單調遞增.本文基于分子動力學模擬觀察到了氫鍵數(shù)隨電場強度的非單調變化行為，展現(xiàn)了微觀分子模擬對研究水分子結構與動力學的重要性.水分子傾向于形成四配位的氫鍵結構，但是對于液相水，熱漲落使得實際形成的氫鍵數(shù)往往小于4，從而產生氫鍵網絡缺陷.根據(jù)文獻[2]的結果，這種氫鍵網絡缺陷與大角度跳躍轉動密切相關.以上結果說明，隨著電場強度的增加，氫鍵網絡結構偏離標準的正四面體配位結構,導致更多的氫鍵網絡缺陷.隨著電場強度的進一步增加，外電場引入的靜電相互作用傾向于導致和電場方向匹配的新四面體序，降低氫鍵網絡缺陷，進而減小大角度跳躍轉動頻次.特別是隨著靜電場的增強并形成新的四面體氫鍵網絡結構，會觀察到液固相變(見下文討論).

可以由下式定義的正四面體序參量Q定量描述正四面體的形成程度[6]為

(2)

四面體序參量Q的平均值隨電場強度變化圖線 偶極矩對關聯(lián)的徑向分布

2.2 擴散系數(shù)的電場依賴性

根據(jù)以往研究結果，體相水分子的大角度跳躍轉動與水分子的擴散過程密切相關[2].圖5給出了不同電場強度下水分子的均方位移(MSD=〈|x(t)-x(0)|2〉，其中x(t)表示t時刻某個水分子的坐標，而x(0)表示該水分子在初始時間點的坐標；〈…〉表示對水分子以及時間的平均.)以及通過直線擬合得到的擴散系數(shù)D.其中，均方位移MSD隨時間t的變化采用了三維空間正常擴散形式，即MSD=6Dt.可以看到，MSD的斜率以及對應的擴散系數(shù)D隨著電場強度先增加，在0.4 V·nm-1附近達到極大值，然后隨電場強度降低.值得一提的是，在Zong等人的理論與分子動力學模擬研究中[4]，觀察到了外加靜電場能夠導致體相水分子系統(tǒng)黏性的各向異性以及非單調變化，這與本文得到的結果相符合.另外，比較圖3與圖5，我們可以看到，這種擴散系數(shù)隨電場強度的變化趨勢與大角度跳躍轉變的變化趨勢相一致，進一步支持大角度跳躍轉動在水分子擴散過程中的關鍵作用[2].

不同電場強度下水分子的均方位移(MSD)隨時間的演化 擴散系數(shù)D隨電場強度的變化圖線

3 靜電場環(huán)境下的液固相變

3.1 電致結冰

以上關于水分子動力學的討論主要集中在小電場區(qū)間.當進一步增強電場時，會觀察到外電場導致的液固相變，即電致結冰[3].圖6給出了平均氫鍵數(shù)以及正四面體序參量Q隨電場的變化.可以看到大電場下平均氫鍵數(shù)目單調遞增的趨勢，并且在電場強度處于2～3 V·nm-1之間存在一個轉變點.經過此轉變點, 水分子快速逼近對應氫鍵飽和構型的四配位氫鍵結構，從而出現(xiàn)了有序的冰相結構.相應地，平均四面體序參量也在轉變點附近發(fā)生了快速增加.以上結果說明在電場強度處于2～3 V·nm-1之間，體相水系統(tǒng)發(fā)生了液固相變，即電致結冰現(xiàn)象.這種電致結冰現(xiàn)象已在以往的實驗和分子動力學模擬中觀察到[3].

大電場下平均氫鍵數(shù) 平均四面體序參量Q隨電場強度的變化

3.2 不同冰相的層狀混合堆疊結構

通過對以上分子動力學模擬數(shù)據(jù)進行詳細分析，可以研究電致結冰得到的冰相結構特征.在通常條件下，液態(tài)水結冰得到的冰相有六方相冰Ice h 和四方相冰Cubic Ice.這兩種主要冰相均呈現(xiàn)正四面體結構，即四面體序參量Q為1，因此無法用Q值分辨這兩種冰相.四面體序參量Q作為標量序，僅能表征構型接近正四面體的程度，無法表征該四面體的空間取向.而判定冰相則需要鄰近四面體的空間取向排列信息，也即需要定義一個包含四面體取向信息的矢量序，并對其統(tǒng)計近鄰關聯(lián).基于這樣的思路，Moore等人[7]提出了CHILL算法來解決冰相結構的定量表征問題，將四面體取向投影到球諧函數(shù)作為單個四面體的矢量序，并由此定義近鄰四面體矢量序關聯(lián)系數(shù)aij：

(3)

其中

(4)

對于室溫液態(tài)水、無定形冰、Ice h和Cubic Ice冰相，四面體矢量序關聯(lián)aij的概率分布具有明顯的區(qū)別.例如，對于Cubic Ice冰相中的水分子，其aij的分布表現(xiàn)為-1.0處的單峰；而Ice h冰相中水分子的aij分布同時在-1.0和-0.11處出現(xiàn)兩個峰.圖7給出的是在電場強度為0, 2, 3 V·nm-1時水分子的aij分布.在電場強度為0和2 V·nm-1時，aij在-1到+1的較大范圍內均具有顯著分布，表現(xiàn)出液態(tài)水的特性.而當電場強度增加為3 V·nm-1時，不僅可以觀察到-1.0處對應Cubic Ice冰相的峰，同時觀察到-0.5附近的另一個峰.這說明電致結冰所得到的冰相并非是單純的Ice h或Cubic Ice冰相，而是存在不同峰形表征出的混合冰相.圖7中灰色曲線表示3 V·nm-1下，不同的單條分子動力學軌跡得到的aij分布.結果顯示不同軌跡得到兩種冰相的比例具有較大的漲落，但雙峰分布特征是一致的.

圖8給出的是電場強度為3 V·nm-1時的電致結冰示意圖.其中圖8(右)展示了分子動力學模擬得到的兩相混合堆疊冰相代表性結構，其中左側區(qū)域的 Cubic Ice冰相，對應于圖7中-1處的峰.右側區(qū)域為較少報道的Polar Ice B冰相[3].可以驗證，這種Polar Ice B對應圖7中-0.5處的峰.

圖7 不同電場強度下序參量關聯(lián)aij分布

圖8 電場強度為3 V·nm-1時的電致結冰示意圖[在電場作用下，體相水由無序的初始構型(左)形成沿電場方向堆疊的混合冰相結構(右)]

人們通常認為具有高平移對稱性的均一冰相自由能最低.最近，Lupi等學者[8]發(fā)表在Nature的工作報道了不一樣的結論：由于熵效應，多種冰相混合堆疊的結構具有更低的自由能.本文中電場環(huán)境下觀察到的Cubic Ice與Polar Ice B的混合堆疊結構與Lupi等人的結論相符合.值得強調的是，文獻[8]中的兩種冰相的堆疊是各向同性的，而本文電致結冰得到的冰相堆疊方向存在明顯的各向異性.圖9給出了電場強度為3 V·nm-1下分子動力學模擬得到的沿電場方向X不同位置處水分子的aij散點圖，以及沿X、Y、Z方向的平均值.可以看到，沿電場方向X，aij平均值在-1.0和-0.5兩種取值交替出現(xiàn).而沿垂直于電場的Y、Z方向，aij平均值為處于-0.5與-0.1之間的恒定值.這說明沿著電場方向存在兩種冰相的堆疊，而垂直于電場方向上不存在堆疊.以上結果展示了電致結冰過程中存在外電場導致的不同冰相堆疊方向的對稱性破缺現(xiàn)象.

圖9 冰相序參量關聯(lián)〈aij〉沿X、Y、Z方向不同位置的平均值(沿電場方向X存在不同冰相的混合堆積，灰色數(shù)據(jù)點表示沿X方向aij的散點圖)

4 結語

本文基于分子動力學模擬，研究了不同電場強度下，體相水分子的大角度跳躍轉動、擴散以及電致結冰.相關結果首次展現(xiàn)了小電場區(qū)域大角度跳躍轉動頻次存在非單調電場依賴性.這種大角度跳躍轉動的反常電場依賴性與水分子擴散系數(shù)隨電場強度變化的非單調性存在較強相關性.在大電場區(qū)域，分子動力學模擬重現(xiàn)出前人研究中得到的電致結冰現(xiàn)象，并發(fā)現(xiàn)了電致結冰導致的冰相體系中，存在兩種不同冰相沿電場方向堆疊的新現(xiàn)象.該分子動力學模擬涵蓋水分子系統(tǒng)的擴散、轉動運動、相變等多個物理過程，有助于學生從微觀層次理解分子擴散與相變等物理現(xiàn)象的分子機制[9,10]，可作為教學實例培養(yǎng)本科生實踐創(chuàng)新能力和科研興趣.

致謝：作者感謝南京大學王煒教授的指導和討論.