冰水界面動態結構的分子動力學模擬研究

王瑞，任瑛，陳衛，韓永生，3

（1 中國科學院過程工程研究所，多相 復雜系統國家重點實驗室，北京 100190； 2 中國科學院大學化工學院，北京 100049；3 中國科學院粉體材料技術重點實驗室，北京 100190）

引 言

冰是水的常見相態，冰的形成和融化過程對水的循環利用和其他的地球化學反應有重要意義，包括生物納米材料中抗凍蛋白與冰的相互作用[1]、過度積冰對戶外設備方面的影響以及抗結冰材料的研發[2]等。水和冰有著諸多反常的特性，如自然界中的液態水在4℃達到密度最大值：溫度高于4℃時，隨著溫度升高，水的密度減小；溫度低于4℃時，隨著溫度升高，水的密度增加。而在水凝固成冰時體積變大，密度變小，不同于一般物質熱脹冷縮。在實驗條件下，所有結晶冰的表面上的水分子都是不均勻的，表面層的結構、厚度和性質等隨著溫度和晶面等條件的改變發生變化，從而影響冰晶表面的生長取向，導致了冰晶的多種形態[2-3]。研究冰水相變過程中冰水界面處的動態變化規律對進一步理解冰晶成核和生長過程十分重要，也有利于進一步解決結冰造成的工程應用上的難題，但是對于冰水界面結構變化規律在分子尺度上的實驗研究仍然存在一定的困難。

隨著計算機技術的發展，模擬計算已經成為科學研究中繼理論方法、實驗方法之后的第三種研究方法，其可以用來在分子層面上解釋一些實驗不能揭示的現象和機制[4-5]，例如納米界面性質與結構的關系[6]和晶體表面自組裝生長過程[7]。冰晶在液態水中的生長動力學，特別是冰晶在不同晶體取向間的各向異性，是一個需要在分子尺度上詳細研究的重要課題。分子動力學（molecular dynamics，MD）模擬可以在分子尺度上研究冰的成核、不同冰相之間的轉化規律以及氣液固三相中冰晶動態生長過程，可以實現多元條件下冰晶生長和熔化過程、冰-水以及冰-水蒸氣兩相平衡態結構變化等方面的模擬研究[8]。

目前已有一些針對冰晶與水的界面和冰晶表面的MD模擬研究，Karim等[9-11]使用分子動力學模擬對六角冰（ice hexagonal,Ih）和液態水的界面進行了模擬，針對水分子的簡單點電荷模型[12]（simple point charge, SPC）和四位點模型[13]（the transferable intermolecular potential four-point,TIP4P），模擬了冰水界面的平衡狀態，討論了界面穩定性。Svishchev等[14]利用TIP4P 水分子模型研究了在較大的外部電場條件下過冷的液態水中冰的結晶和生長過程，其冰水界面受到外加電場的影響，進一步影響了冰晶的幾何形狀和生長速率。近年來，研究者對冰水界面的研究越來越深入，涉及的層次也越發廣泛。在冰水界面，密度和自由能等熱力學性質發生了顯著的變化，如Wang 等[15]在環境壓力1 bar（1 bar=0.1 MPa）下計算出不同模型水的固液界面自由能約為37 mJ·m-2（the transferable intermolecular potential four-point with Ewald techniques, TIP4P-Ew）和42 mJ·m-2（TIP5P-Ew），與之前用MD模擬計算的測量范圍25～44 mJ·m-2一致。Chung 等[16]研究了極化水分子模型在不同相態之間的電荷轉移情況及其對界面性質的影響，如熔化溫度、介電特性、能量和密度等，認為電荷轉移對水分子內和水分子間的氫鍵網絡作用有重要意義。Louden 等[17]開展了冰水界面摩擦的非平衡分子動力學模擬，認為固液界面的氫鍵密度與界面摩擦力直接相關。基面{0001}、一次側棱面{10-1 0}和二次棱柱面{11-2 0}的摩擦力相對于表面特征的剪切速率和剪切方向是不變的，三個晶面的結構和動力界面寬度都是相似的，并且很大程度上獨立于剪切速率和方向。Ramírez等[18]從氫鍵拓撲網絡結構的角度出發分析了水分子的氣-液和固-液相變過程中的一些界面性質，如平均密度、聚類系數和平均路徑長度等。這些性質在接近相變區域時表現出不連續的特征，這與在同一區域中某些熱力學性質所表現出的不連續性是相似的，為進一步研究其他含有氫鍵結構的兩相界面提供了更多依據。目前針對冰水界面結構的研究都是在常溫常壓下進行的，主要是為了在水分子模型的熔點附近得到冰水共存的平衡狀態，進而分析界面熱力學性質和結構性質。目前得到的普遍認識是界面靜態結構的不均勻性和動態分布的不連續性都是區別于體相特征的，也是影響液態水中冰成核、生長等過程的關鍵因素[19]，但是界面結構在不同溫度和晶面條件下的動態變化規律較為復雜，仍需要對其進行深入的研究。

本文利用分子動力學模擬方法構建冰水兩相共存模型，在微正則（NVE）[20]系綜下模擬冰水界面結構的動態變化。選擇NVE 孤立體系可以避免控溫器和控壓器施加的力對動力學過程造成影響。NVE 系綜下體系的總能量不變，沒有外界能量和粒子的輸入和輸出，完全依靠體系內冰水兩相之間的能量交換。研究界面層內水分子在無序和有序狀態之間的轉變情況，分析運動軌跡和粒子成鍵情況，總結界面結構的變化規律。利用序參量和密度剖面表征界面結構的動態變化過程，在不同溫度和晶面條件下分析界面厚度和性質的變化規律。

1 模擬過程

本工作在超算機群Mole-8.5[21]上通過開源軟件LAMMPS（7Aug19）進行算例計算，采用CUDA 版本的GPU 節點，搭載Tesla K80 進行加速計算，一個GPU有12個核，最大計算速度是普通CPU 節點的十倍[22]，可以極大地提高運算速度。將包含2880 個水分子的六角冰Ih 作為初始結構，體系盒子的尺寸大小為4.5 nm×4.7 nm×4.4 nm。選擇水分子力場模型為拓展簡單點電荷模型[23]（simple point charge extended, SPCE），這是一個三位點剛性全原子模型。該模型可以更好地描述液態水結構，同時三點電荷結構相比于四點電荷的TIP4P 模型能減少計算量，使用SPCE水分子模型描述六角冰Ih時，其在壓強為1 bar 時的熔點為215 K[24]。初始六角冰結構的密度為0.91747 g·cm-3，模擬盒子三個方向均采用周期性邊界條件，使用SHAKE算法約束水分子空間結構保持。體系為NVE 系綜，給定體系不同的初始溫度，觀察體系結構在孤立系統中的演化。系統的平衡溫度和壓強可以根據SPCE 水分子模型的相圖中水的液相和固相共存曲線確定。

首先對系統進行5000步的能量最小化處理，優化體系結構為更合理的幾何結構。用平面z=0.18 nm將初始冰晶結構分成A 和B 兩部分，垂直z軸的xy平面為六方晶系的基面{0001}，因此冰水接觸的暴露晶面即為基面{0001}。在初始化設置中對A部分的冰晶結構進行固定，即其位置和速度不隨時間更新，而B部分的水分子根據麥克斯韋-玻爾茲曼分布的隨機數列所產生。弛豫20 ps后解除對A部分位置和速度的固定設置，使體系內兩部分進行能量交換，最終體系會重新達到新的平衡態。此時再對整個體系進行溫度和壓強的弛豫，在等溫等壓（NPT）系綜下弛豫200 ps，采用Nosé-Hoover 恒溫器和恒壓器來控制體系的溫度保持恒定，隨著模擬盒子的體積減小，在平衡態時液態水的密度為0.98759 g·cm-3，接近實驗值1.0 g·cm-3。最后再進行NVE系綜下的10 ns預弛豫過程，體系達到固液兩相平衡，體系溫度、壓強等熱力學參量都達到平衡。NVE 孤立體系中存在界面，可以使體系的溫度保持在熔點附近達到動態平衡，同時避免形成共存兩相平衡態造成的體系過冷問題。對B 部分設置不同的初始溫度，經過能量平衡弛豫后，最終體系為冰水共存狀態的溫度，即為相變溫度。若初始溫度高于相變溫度，B 部分的水分子會將過剩的能量傳遞給A 部分的冰，導致冰融化；反之，若初始溫度低于相變溫度，經過能量交換，B 部分水分子運動速度減慢，變成有序冰晶結構，體系逐漸變成全冰相。將軌跡文件導入可視化分子動力學軟件VMD (visual molecular dynamics)[25]中可以觀察在不同平衡溫度時體系結構的變化。在NVE 系綜下利用熱力學參量來驗證體系是否達到平衡態是非常必要的，在溫度和壓強達到平衡后，統計在NVE 系綜下在弛豫的30 ns 內兩相體系的各熱力學參量。體系的總能量不變，在30 ns的平衡弛豫過程中其波動約為0.028%，體系的溫度在熔點215 K附近波動，波動范圍為±2 K。

2 結果與討論

2.1 利用序參量表征體系中冰水界面結構動態變化

冰水兩相共存的界面是冰晶成核和生長的區域，水分子在冰晶結構表面的動態復雜行為影響冰晶成核和生長過程。因此探究冰晶生長前端的界面微環境首先要對界面結構有清晰的認識。本文采用序參量描述冰水界面結構的變化[26-28]，冰的晶格結構是正四面體的對稱結構[29]，但是在液相中由于水分子的無序移動，四面體結構會發生一定程度的扭曲。四面體序參量qi可以用來描述四面體的這種對稱性[30]，qi定義為

式中，ψjk為一個給定分子的氧原子與其最近鄰的j和k氧原子之間的連線夾角,(°)。如果一個分子位于正四面體的中心，其頂點被最接近的四個水分子中的氧原子所占據，則cosψjk=-1/3。

當qi值接近1 時，表明四面體的對稱性較高，此時水分子組成高度有序的冰晶結構；而當qi的值在0附近時，表明四面體的對稱性較低，此時水分子組成無序液態水。因此計算四面體序參量可以在一定程度上區別冰水混合結構中類冰水分子和類液水分子，選取合適的序參量的值可以更好地描述冰水界面處水分子的結構轉變以及界面結構波動變化規律。

首先分別對冰和水體系進行序參量計算，對比液態水和固相冰晶的序參量q6概率分布曲線。結果發現，在截斷半徑rc=0.3 nm 時，兩者的區分程度最明顯，對應的q6閾值為0.41。當序參量q6取值大于0.41時，該對應的水分子即為類冰水分子；當序參量q6取值小于0.41 時，該對應的水分子即為類液水分子，從而可將體系中的類冰水分子和類液水分子區分開來。

圖1 為本研究的基本思路，采用序參量來表征液態水中冰晶生長過程中界面結構動態變化規律，如圖1(a)所示；圖1(b)從左到右依次為液態水、冰水共存和固態冰晶三種狀態的結構圖及其對應的序參量概率分布曲線。

圖1 采用序參量表征冰晶生長過程中冰水界面微觀結構的動態變化過程(a)；液態水、冰水共存和固態冰晶體系分子結構及其水分子序參量概率分布曲線(b)Fig.1 The complex dynamic process of the interface microstructure in the process of ice crystal growth in liquid water characterized by order parameters(a);The structure of liquid water,ice-water and ice crystal with probability distribution curve(PDC)of order parameters for water molecules(b)

進一步，利用氧-氧徑向分布函數驗證不同狀態時的體系結構。徑向分布函數可用來表征體系內固體或液體分子排列的有序程度，定性判斷體系結構。其表達式為：

式中，r為任意粒子與球心(目標粒子)的距離，nm；ρ為體系的密度，g·cm-3；dN是厚度為dr的球殼內的顆粒數量。徑向分布函數可以看作整個體系內局域密度與平均密度的比值[31]。

在體系達到平衡態后，統計體系結構的氧-氧徑向分布函數的時間平均值，結果如圖2 所示。體系結構隨溫度變化明顯，體系溫度低于熔點215 K時，體系的特征峰明顯，與文獻[32]中的冰晶結構一致，此結構長程有序，體系為純冰狀態；體系溫度高于熔點時，第一峰值和第二峰值的強度逐漸降低，且遠程峰的晶格結構遭到破壞，分子間作用力減弱，特征峰消失，曲線波動趨勢減小，逐漸向1靠近，即曲線歸一化明顯。體系溫度在熔點215 K 附近時，氧-氧徑向分布函數介于兩者之間，可以判斷為冰水兩相共存狀態。

圖2 不同溫度下體系的氧-氧徑向分布函數對比Fig.2 Comparison of the radial distribution function gO—O(r)of the system at different temperatures

2.2 冰晶界面結構動態演變

利用LAMMPS 構建冰水兩相共存模型，冰水接觸的界面為{0001}。經過10 ns的預弛豫之后，在熔點附近再進行30 ns的平衡弛豫過程，得到體系中水分子的運動軌跡，在VMD 中觀察其三維動態變化過程。

冰水兩相達到平衡狀態后，選取界面處某一固定位置觀察其結構隨時間的變化，可以得到不同時刻的界面結構。如圖3(a)、(b)所示，以t=1.5、2.3、5.6、7.9 ns四個時刻為例，t=1.5 ns 時，界面該區域全部為冰晶的六元環結構；t=2.3 ns 時，該區域出現了部分無序結構；t=5.6 ns時，該區域全部變為無序的分散結構；t=7.9 ns時，該區域又重新恢復部分有序的冰晶六元環結構。

模擬結果表明，在微觀原子尺度上，冰水兩相平衡后的界面仍在發生不同程度的移動，界面處無序水分子被吸附在冰面上，與冰層表面的水分子結合生成氫鍵進而排列成有序的冰晶六元環結構，而同時界面層內冰晶結構的水分子也發生分子間氫鍵的斷裂，脫離晶體表面進入界面層，擴散到液態水分子的體相之中，這兩種變化時刻發生。界面層內同時包含無序的液態水分子結構和有序的冰晶結構，兩者之間存在可逆動態轉變。進一步，通過序參量q6定量統計體系中類冰水分子數目，得到其隨時間的變化規律，如圖3(c)所示，發現類冰水分子數目在動態平衡過程中也呈現波動變化，證明冰水界面處發生相變過程，氫鍵的形成和斷裂同時進行，并存在局部的能量交換。

圖3 不同時間冰水界面處水分子結構(a)；冰水界面示意圖(b)；兩相平衡態時界面處類冰水分子數目隨時間的變化(c)Fig.3 The structures of the ice-water interface at different time(a);Schematic diagram of the ice-water interface(b);The number of ice-like water molecules changes with times at the two-phase equilibrium condition (c)

2.3 溫度對冰水界面結構的影響

利用序參量表征冰水界面結構隨著溫度的變化規律，初始狀態為冰水兩相共存的平衡態，選擇冰水接觸面為二次棱柱面{11-2 0}，見圖4。

溫度為215 K 時，體系中水分子的序參量概率分布曲線為明顯的雙峰分布，表明體系中冰水兩相共存。在降溫過程中，右側分布逐漸增加，說明類冰水分子在不斷增加；左側分布逐漸減小，最后變為靠右側的單峰，說明類液水分子逐漸轉化為類冰水分子，如圖4(a)所示。在升溫過程中，序參量的分布逐漸變窄且分布整體左移，同時左側分布增加，右側分布減少，最后變為靠左側的單峰，說明此時類冰水分子在逐漸轉化為類液水分子，如圖4(b)所示。

圖4 不同溫度下體系水分子序參量概率分布曲線：（a）降溫過程；（b）升溫過程Fig.4 Probability distribution curve(PDC)of order parameters at different temperatures:(a)cooling process;(b)heating process

進一步，統計沿垂直于界面方向水分子的密度變化，得到如圖5 所示不同溫度時體系的密度剖面對比圖。在過冷度較大時，界面密度波動較大，說明界面結構趨于有序狀態。因此界面內水分子密度較大，水分子在界面層內的停留時間較長，形成六角冰結構的概率增大。同時冰相的峰強度更高，其結構高度有序，界面層厚度較薄，冰水界面區分明顯。隨著溫度升高，界面層厚度增加，界面內水分子密度減小，無序性增加。可以預測，隨著溫度逐漸升高到熔點以上，界面處的峰強度逐漸降低，最終結構密度在1.0 g·cm-3（液態水的密度）處波動。

圖5 在不同溫度下二次棱柱面{11-2 0}的界面密度剖面圖Fig.5 Density profile of the ice-water interface formed by secondary prismatic plane{11-2 0}

2.4 晶面對冰水界面結構的影響

六角冰的生長晶面主要有三個，基面{0001}、一次棱柱面{10-1 0}和二次棱柱面{11-2 0}。其中一次棱柱面{10-1 0}和二次棱柱面{11-2 0}的暴露界面具有幾何粗糙結構，而基面{0001}具有分子平面結構。基面和棱柱面的雙層高度分別為0.37 nm 和0.39 nm。在冰水共存體系中，界面處存在大量吸附的水分子，其與周圍的水分子沒有形成理想的六方晶系結構，但其分子間氫鍵作用力遠大于體相中液態水分子之間的相互作用力，界面冰晶水分子的排列方式不同，即暴露的氫鍵數目和角度不同，會對界面處液態水分子產生不同取向和大小的作用力。

針對以上三個晶面分別建立冰水兩相共存體系，利用序參量q6定量表征不同晶面形成的冰水界面結構變化規律。三個體系中以基面為接觸面形成初始兩相結構時的類冰水分子和類液水分子數目分別為1650和1230，以一次棱柱面為接觸面形成初始兩相結構時的類冰水分子和類液水分子數目分別為1480和1400，以二次棱柱面為界面的初始兩相結構時的類冰水分子和類液水分子數目分別為1536 和1344。在經過弛豫達到兩相共存的平衡態時，分別統計三個體系中類冰和類液水分子的序參量，得到圖6 的概率分布曲線。其中體系的冰水接觸面分別是(a)基面、(b)一次棱柱面、(c)二次棱柱面。經過時間平均和系綜平均后得到兩相平衡狀態時類冰水分子和類液水分子數目。分別對比初始兩相結構，發現一次棱柱面和二次棱柱面形成的兩相體系中類冰水分子和類液水分子數目更接近，在孤立體系的模擬環境中也更穩定。

圖6 不同晶面形成的冰水兩相體系的序參量概率分布曲線，數字分別代表類冰水分子（右）和類液水分子（左）數目Fig.6 Probability distribution curve(PDC)of the order parameters in ice-water system formed by three different crystal planes,and the figures in the histogram represent the number of water molecules at liquid (left)and solid phases(right),respectively

進一步，沿著垂直于界面方向計算密度剖面觀察界面結構，如圖7所示。其中圖7(a)為基面形成的冰水界面結構。其表面趨向于平坦，這與實驗[33]觀察結果一致。圖7(b)為一次棱柱面{10-1 0}對應的冰水界面結構，與基面相同，一次棱柱面的每一層都有一個雙峰，第一個雙峰含有六角形環上的三個氧，第二個雙峰與其他三個氧對應。圖7(c)為二次棱柱面{11-2 0}對應的界面結構，其密度剖面為單峰，這是由界面水分子的排列結構決定的。

圖7 不同晶面形成的冰水界面的密度剖面圖Fig.7 Density profile of the ice-water interface formed at different crystal planes

通過計算三個不同晶面從冰晶結構過渡到液態水分子的界面層厚度，發現在相同的溫度下，一次棱柱面和二次棱柱面的界面層厚度分別為0.35 nm 和0.20 nm，小于基面的界面層厚度0.40 nm。基面的邊緣結構平坦，界面層較厚，氫鍵在晶面上發生二維重組；棱柱面形成的界面層比基面窄，界面水分子的有序性強，氫鍵在晶面上發生三維重組，邊緣結構接近六元環。

本研究在分子尺度上闡明了體系間生長過程中界面結構的各向異性，即不同晶面的生長速度不同會導致冰晶生長形態的多樣性。因此不同于氣相雪花形成過程中基面生長速度比棱柱面快，液態水中冰晶在過冷度相同的條件下其棱柱面生長速度更快，傾向生成側枝較長的枝晶。

3 結 論

本文利用分子動力學方法模擬研究了冰水兩相共存體系中冰水界面結構的動態變化過程，揭示了不同溫度和晶面條件下界面結構的變化規律，通過序參量、徑向分布函數和密度剖面等對該過程進行了定量表征，得到主要結論如下。

（1）在NVE系綜下，當冰水兩相達到平衡態時，冰水界面結構在有序冰晶六元環和無序水分子之間可逆動態轉變。

（2）晶面類型影響冰水界面動態結構，基面結構較為平坦，形成的冰水界面層較厚，分子間氫鍵在晶面上發生二維重組，水分子排列的有序性較差；而一次和二次棱柱面結構粗糙，分子間氫鍵在晶面上發生三維重組，有序性更強。

（3）過冷度對冰水界面動態結構及其演化趨勢有重要影響，過冷度越大，界面結構越趨于冰晶六元環的有序狀態，界面內水分子密度越大，水分子在界面層內停留時間越長，晶面生長的概率越大。