趙崇理， 趙東霞， 楊忠志

(遼寧師范大學(xué) 化學(xué)化工學(xué)院，遼寧 大連 116029)

水作為自然界最常見的物質(zhì)被廣泛地研究[1].作為水固態(tài)形式的冰，人們研究的相對較少[2].計算機模擬經(jīng)常用來幫助解釋實驗中分子層面上的運動過程，如冰的成核、冰的相間轉(zhuǎn)化以及過冷液-液相平衡等.利用任何一種水分子模型模擬相圖都是具有挑戰(zhàn)性的.由于冰中的水分子間強的氫鍵作用，常見冰的密度通常小于1.0 g/cm3.冰密度的實驗值經(jīng)常被用來作為調(diào)節(jié)水分子力場參數(shù)和檢驗水分子模型優(yōu)劣的標(biāo)準(zhǔn).TIP7P水模型[3]采用了7個電荷位點，對水分子的電荷分布劃分更細致，能準(zhǔn)確模擬液相水密度、汽化熱、介電常數(shù)等性質(zhì)，同時作為固定電荷力場計算效率比極化力場高.

目前，用于模擬估算晶體熔點的方法通常有4種，分別是“絕對”自由能法[4]、Gibbs-Duhem積分法[5]、固-液相共存法[6]以及自由表面法[7].其中，前兩種是基于Gibbs自由能的方法，平均誤差為±6.0 K；后兩種是直接模擬相表面的方法，平均誤差為±2.0 K[8].自由表面法準(zhǔn)確度高，而且便于觀察融化過程，所以本文選擇了自由表面法估算TIP7P模型下冰的熔點.

1 計算方法

1.1 TIP7P水分子模型

TIP7P水分子為剛性模型，分子的鍵長鍵角保持不變，分子內(nèi)相互作用能為0.TIP7P水分子體系中，分子i和j間的勢能表達式如下：

(1)

其中，a和b為原子編號，c和d為電荷位點編號，εab和σab分別為范德華作用勢阱深度和距離參數(shù)，qc和qd為位點電荷，rab(或rcd)為原子(或位點)間距離.范德華參數(shù)εab和σab的混合規(guī)則均為幾何平均.TIP7P水分子力場的詳細參數(shù)見文獻[3].

1.2 冰的建模

按照Buch等人[11]的方法構(gòu)建具有Ih群對稱性冰.首先在六面體晶胞(a=b=0.451 8 nm,c=0.735 6 nm,α=β=90°,γ=120°)中，按晶胞分?jǐn)?shù)坐標(biāo)±(1/3,2/3,1/16)和±(2/3,1/3,9/16)放置4個氧原子，這樣建出的冰晶密度約為0.920 g/cm3.然后將晶胞復(fù)制平移得到(8,8,4)晶胞陣列，整個冰塊中有1 024個氧原子，1個氧原子與其他4個氧原子相鄰且呈正四面體，2個相鄰氧原子間距約0.276 nm.接下來要在周期性邊界條件下添加氫原子，滿足兩個原則：每2個相鄰氧原子之間有1個氫原子；1個氧原子與2個氫原子成共價鍵.由于建模中氫原子是無序排列的，所以加上第三個原則：整個冰塊總偶極矩為0.具體做法：先將每2個相鄰氧原子之間放置1個氫原子隨機靠近其中1個氧原子，然后利用蒙特卡羅方法隨機生成水分子以滿足第二個原則.其中，蒙特卡羅方法嘗試選取1個氫原子從氧原子a切換到與另1個氧原子b成鍵，計算2個氧原子a,b的成共價鍵數(shù)之差|na-nb|.當(dāng)|na-nb|增加時，該嘗試被禁止；當(dāng)|na-nb|減少時，該嘗試被接受；當(dāng)|na-nb|不變時，接受的概率為1/2.當(dāng)檢測到所有氧原子成2根共價鍵則停止嘗試.不斷變化初始隨機數(shù)種子，計算整個冰塊的偶極矩，當(dāng)結(jié)果接近于0時即滿足了第三個原則.

1.3 計算細節(jié)

模擬計算在GROMACS[12]分子動力學(xué)軟件包中實現(xiàn).將包含有1 024個水分子的冰塊進行恒溫恒壓模擬.SETTLE算法限制每個水分子幾何構(gòu)型變化而成為剛體.范德華和靜電相互作用的截斷半徑均為1 nm，長程相互作用均用SPME方法來修正.動力學(xué)采用蛙跳積分，時間步長為1 fs.溫度控制采用Nosé-Hoover熱浴，溫度弛豫時間為0.5 ps.壓強控制采用Parrinello-Rahman壓浴，3個維度上分別維持在1 bar，壓強弛豫時間為1.0 ps.通過4 ns的預(yù)模擬來確定冰塊的三維尺寸.將原有冰塊在z的上下兩方向分別擴大1倍真空空間，x和y方向尺寸保持不變，構(gòu)成冰表面模型體系，然后進行10 ns恒溫恒容模擬.本文模擬了232、236、238和240 K 4個溫度，并觀察其內(nèi)能和晶格對稱性的變化.

2 結(jié)果與討論

圖1中顯示溫度為232 K和236 K冰始終晶格排列規(guī)則(見圖2(A))，其內(nèi)能未有明顯突躍，冰未融化；而溫度為238 K大約8 ns時和溫度為240 K大約5 ns時內(nèi)能發(fā)生突躍，晶格排列被打亂(見圖2(B))，冰由表面至內(nèi)部全部融化.觀察到的最低的融化溫度為238 K，最高的未融化溫度為236 K，取平均值即為冰的熔點237 K，誤差為±1 K.

由恒壓模擬得到的體相冰的密度0.908 g/cm3接近于實驗值0.922 g/cm3.如果只是模擬體相冰，由于水分子受到作用力比較對稱，所以比較難融化，即容易過熱.但在筆者構(gòu)建的冰的表面體系中，冰表面一側(cè)是真空沒有作用力,另一側(cè)是冰內(nèi)部的水分子，其對表面水分子進行拖拽來維持冰表面.冰表面總有一些水分子由于分子熱運動，擺脫了與冰內(nèi)部水分子之間的氫鍵作用而自由旋轉(zhuǎn)，但仍未遠離冰的表面.冰表面具有的這一薄液相層的現(xiàn)象，稱為預(yù)融化，可以有效降低冰表面的張力，已有相關(guān)的實驗測定和模擬計算報道[13].在圖2(A)中可以看到冰內(nèi)部結(jié)構(gòu)保持完好，但兩側(cè)表面已有部分水分子融化，這些水分子像在液相中那樣運動.

圖1 平均每個分子的內(nèi)能隨時間變化的關(guān)系

值得注意的是，在模擬的過程中需要保持模擬體系的周期性尺寸不變.一方面，可以維持z方向(圖2中的水平方向)真空區(qū)域的大小；另一方面，冰表面融化的薄液層在x和y方向自由運動會引起這兩個方向壓強波動.如果調(diào)節(jié)這兩個方向的尺寸來維持壓強的話，冰內(nèi)部會受影響，從而使冰相從內(nèi)部融化.所以為了準(zhǔn)確地觀察冰的融化情況，先通過體相的冰恒溫恒壓模擬來確定兩個方向的尺寸，之后的模擬一定要保持體系尺寸固定不變.

圖2 冰薄層經(jīng)10 ns模擬后出現(xiàn)兩種情況 (A)未融化和(B)完全融化

3 結(jié)論與展望

使用TIP7P水分子模型模擬冰表面的融化過程，得到的冰的熔點(237±1) K具有較小的誤差.所有溫度下冰的表面都有預(yù)融化層.當(dāng)升高溫度時，預(yù)融化層厚度增加.在熔點附近，融化層迅速向冰的內(nèi)部擴張.TIP7P水模型能較準(zhǔn)確地模擬冰的密度，而在構(gòu)建TIP7P模型時未考慮模擬冰的性質(zhì)，這說明TIP7P水模型具有一定的可轉(zhuǎn)移性.為TIP7P模型冰的模擬繪制相圖積累了經(jīng)驗.