吳方棣 鄭輝東 吳燕翔

(1.武夷學(xué)院 環(huán)境與建筑工程系，福建 武夷山 354300；2.福州大學(xué) 化學(xué)化工學(xué)院，福建 福州 350108)

水性質(zhì)的分子動(dòng)力學(xué)模擬研究

吳方棣1鄭輝東2吳燕翔2

(1.武夷學(xué)院 環(huán)境與建筑工程系，福建 武夷山 354300；2.福州大學(xué) 化學(xué)化工學(xué)院，福建 福州 350108)

本文采用SPC/E水模型對(duì)水的性質(zhì)進(jìn)行了分子動(dòng)力學(xué)模擬研究。結(jié)果表明，與常態(tài)及亞臨界狀態(tài)下相比，超臨界下水的密度隨溫度壓力的變化更為顯著，這是超臨界水具有較高溶解能力的重要原因；且在超臨界下其自擴(kuò)散系數(shù)明顯增大，表明在超臨界下水有著更高的擴(kuò)散能力。在常態(tài)下水分子主要形成了二至四個(gè)氫鍵的結(jié)構(gòu)，占了70%以上；隨著溫度的升高，每個(gè)水分子平均形成的氫鍵數(shù)減少。水表面張力模擬值與實(shí)驗(yàn)值相比偏小,其原因是未考慮氫鍵等因素對(duì)表面張力的影響。

水；密度；擴(kuò)散系數(shù)；氫鍵度；表面張力；分子動(dòng)力學(xué)模擬

1 引言

分子動(dòng)力學(xué)(Molecular Dynamic，簡(jiǎn)稱MD)是研究在某一條件下分子體系的各種性質(zhì)隨時(shí)間變化，它是以特定粒子(如原子、分子等)為基本研究對(duì)象，將系統(tǒng)看作具有一定特征的粒子集合，運(yùn)用經(jīng)典力學(xué)和統(tǒng)計(jì)力學(xué)方法，通過(guò)研究微觀分子的運(yùn)動(dòng)規(guī)律，得到體系的宏觀特性和基本規(guī)律[1]。分子動(dòng)力學(xué)模擬作為一種理想的計(jì)算機(jī)實(shí)驗(yàn)方法，隨著計(jì)算機(jī)技術(shù)的發(fā)展，近年來(lái)其廣泛應(yīng)用于化學(xué)化工、生命科學(xué)、醫(yī)學(xué)、力學(xué)及材料科學(xué)等各個(gè)領(lǐng)域。

物質(zhì)物性數(shù)據(jù)的獲得通常有實(shí)驗(yàn)測(cè)定、模型計(jì)算以及近年來(lái)迅速發(fā)展起來(lái)的計(jì)算機(jī)模擬等方法。分子動(dòng)力學(xué)模擬不僅可以獲得物質(zhì)的平衡性質(zhì)，還可以獲得物質(zhì)的傳遞性質(zhì)如：擴(kuò)散系數(shù)等。水是我們?nèi)粘Ｉ钪胁豢扇鄙俚奈镔|(zhì)，它在化工生產(chǎn)過(guò)程中被廣泛用作溶劑、冷卻劑、萃取劑等。因此近年來(lái)廣大學(xué)者對(duì)其性質(zhì)也進(jìn)行了較為深入的研究。本文采用分子動(dòng)力學(xué)模擬方法，對(duì)不同條件（從常態(tài)至超臨界狀態(tài)）下水的密度、自擴(kuò)散系數(shù)、氫鍵數(shù)以及表面張力做了較為全面的模擬研究。

2 模擬方法

2.1 勢(shì)能模型選擇

現(xiàn)常用的水的勢(shì)能模型有SPC、SPC/E、TIP4P等，其分子結(jié)構(gòu)及勢(shì)能參數(shù)如表1，其中A、B為L(zhǎng)J參數(shù)。

表1 水模型參數(shù)表Table 1 The parameters of differentwatermodel

為了選擇更為適合的水模型，我們分別計(jì)算了25℃、1atm下這三種模型的密度和擴(kuò)散系數(shù)，結(jié)果如表2所示。

表2 在25℃、1 atm下不同水模型的密度及擴(kuò)散系數(shù)的模擬值Table 2 The simulative density and diffusion coefficient of different water model at 25℃,1 atm

從表2可以看出，TIP4P計(jì)算的密度與實(shí)驗(yàn)值最接近，但擴(kuò)散系數(shù)相差較大，而SPC/E的密度與擴(kuò)散系數(shù)都與實(shí)驗(yàn)值的誤差很小，因此SPC/E水模型能夠提供更為準(zhǔn)確的數(shù)據(jù)，從而本文采用SPC/E水模型進(jìn)行計(jì)算。

2.2 模擬細(xì)節(jié)

本次的分子動(dòng)力學(xué) (MD)模擬采用的是GROMACS4.0.5[4]軟件包。模擬采用NPT系綜，三維周期邊界條件，壓力耦聯(lián)采用berendsen[5]方法，溫度耦聯(lián)采用v-rescale[6]方法。采用PME[7,8]算法進(jìn)行靜電相互作用力計(jì)算，庫(kù)侖截?cái)喟霃綖?0 A°。范德華截?cái)喟霃綖?2A°。對(duì)所有鍵長(zhǎng)采用LINCS[9]進(jìn)行抑制。模擬的時(shí)間步長(zhǎng)為2 fs，每隔1 ps保存一次坐標(biāo)數(shù)據(jù)，為保證體系達(dá)到平衡，每個(gè)計(jì)算體系模擬時(shí)長(zhǎng)為4-6 ns。

3 結(jié)果與討論

本文模擬過(guò)程水分子數(shù)為4176個(gè)，下面分別對(duì)水的密度，自擴(kuò)散系數(shù)，氫鍵度以及表面張力進(jìn)行研究分析。

3.1 水密度的計(jì)算

采用NPT系綜計(jì)算了不同溫度、壓力下水的密度，并與文獻(xiàn)值進(jìn)行了比較，結(jié)果如表3所示。表3中常態(tài)及亞臨界狀態(tài)密度文獻(xiàn)值由文獻(xiàn)[10]查閱和計(jì)算得到；超臨界下密度的文獻(xiàn)值根據(jù)Saul[11]推導(dǎo)得出的水的38參數(shù)狀態(tài)方程計(jì)算得到。

從表3可以看出，模擬值與文獻(xiàn)值很接近，最大相對(duì)誤差不超過(guò)±10%，其精度基本能夠達(dá)到工程設(shè)計(jì)計(jì)算的要求。比較文獻(xiàn)值與模擬值可以看出，模擬的密度總體趨勢(shì)比文獻(xiàn)值要偏小。從表3數(shù)據(jù)還可以看出，在常態(tài)和亞臨界狀態(tài)下，溫度是影響密度的主要因素，而壓力對(duì)密度的影響較?。辉诔R界條件下，溫度和壓力的變化都會(huì)對(duì)密度產(chǎn)生較大的影響，但也依然保持了常態(tài)及亞臨界條件下的隨著溫度升高密度減小，隨著壓力增大密度也隨之增大的趨勢(shì)。

超臨界水密度是除偶極矩外決定其對(duì)有機(jī)溶劑溶解度大小的另一個(gè)重要因素；通常情況下，其密度增大，溶解度增大；密度減小，溶解度減小[12]。根據(jù)這一特性及超臨界下水的密度隨溫度壓力變化顯著的特點(diǎn)，從而在實(shí)際應(yīng)用過(guò)程中，可以通過(guò)調(diào)節(jié)溫度壓力來(lái)改變超臨界水的溶解能力，進(jìn)而達(dá)到萃取、分離等目的。

表3 水密度模擬值Table 3 The simulative density ofwater

3.2 水的自擴(kuò)散系數(shù)計(jì)算

分子動(dòng)力學(xué)模擬通過(guò)對(duì)研究物系中粒子的運(yùn)動(dòng)方程的求解從而得到粒子的運(yùn)動(dòng)速度和運(yùn)動(dòng)軌跡，通過(guò)統(tǒng)計(jì)計(jì)算即可求得擴(kuò)散系數(shù)，由于是計(jì)算機(jī)實(shí)驗(yàn)，它能夠不受高溫高壓的限制，對(duì)于計(jì)算超臨界、深冷等特殊條件下的擴(kuò)散系數(shù)其優(yōu)勢(shì)更加明顯，因此它是研究擴(kuò)散系數(shù)等傳遞性質(zhì)的有效手段[13]?，F(xiàn)今，基于MD計(jì)算溶液擴(kuò)散系數(shù)主要有兩種方法：一種是Einstein法，基于Einstein方程的，通過(guò)統(tǒng)計(jì)均方位移(MSD)與時(shí)間的變化，當(dāng)時(shí)間足夠長(zhǎng)時(shí)，MSD與時(shí)間基本呈線性關(guān)系，從而計(jì)算斜率即可得到溶液擴(kuò)散系數(shù)，其自擴(kuò)散系數(shù)計(jì)算式如式1所示[14,15]。

另一種是Green-Kubo法，通過(guò)計(jì)算體系中組分的速度相關(guān)函數(shù)再積分得到，其計(jì)算式如式2所示[14,15]。

由于Einstein法計(jì)算方法相對(duì)簡(jiǎn)便且準(zhǔn)確性與Green-Kubo法相當(dāng)，因此本文中的擴(kuò)散系數(shù)采用Einstein方法進(jìn)行計(jì)算。

本文用Einstein方法計(jì)算了不同溫度、壓力下水的自擴(kuò)散系數(shù)，結(jié)果如表4所示。表4中常壓及亞臨界狀態(tài)下擴(kuò)散系數(shù)文獻(xiàn)值由Krynicki[16]關(guān)聯(lián)式計(jì)算得出；超臨界條件下擴(kuò)散系數(shù)文獻(xiàn)值由Lamb[17]的實(shí)驗(yàn)數(shù)據(jù)公式擬合得到。

比較模擬值與文獻(xiàn)值可知，水自擴(kuò)散系數(shù)的模擬值與文獻(xiàn)值的最大相對(duì)誤差為-8.98%，其平均相對(duì)誤差在±10%以內(nèi)。因此在實(shí)驗(yàn)數(shù)據(jù)不足的情況下，用分子動(dòng)力學(xué)模擬的方法來(lái)預(yù)測(cè)水的擴(kuò)散系數(shù)是可行的[18]。從表4可以看出，從常態(tài)至超臨界狀態(tài)，水的自擴(kuò)散系數(shù)都是隨著溫度的增加而增大，隨著壓力的增大而減小。與常態(tài)及亞臨界條件相比，在超臨界狀態(tài)下，水有著更高的擴(kuò)散系數(shù)，并且擴(kuò)散系數(shù)隨溫度、壓力的變化更加敏感。

表4 水的自擴(kuò)散系數(shù)Table 4 The self-diffusion coefficient ofwater

3.3 水中氫鍵成鍵度

為了統(tǒng)計(jì)不同條件下的氫鍵變化情況，我們首先要給出合適的氫鍵成鍵標(biāo)準(zhǔn)。氫鍵的定義有分別按能量標(biāo)準(zhǔn)和按幾何標(biāo)準(zhǔn)定義[19,20]，同時(shí)也有采用能量-幾何混合的標(biāo)準(zhǔn)來(lái)定義[21]。本文由于是用坐標(biāo)文件來(lái)統(tǒng)計(jì)氫鍵的，因此擬采用幾何的標(biāo)準(zhǔn)來(lái)定義氫鍵。Teixeira[22]的中子衍射實(shí)驗(yàn)表明，當(dāng)氫鍵偏離O-H鍵軸30°時(shí)，即圖1中θ小于150°時(shí)，氫鍵將被破壞；通過(guò)水分子的徑向分布函數(shù)研究表明，水分子間氧－氫間第一殼層的距離為2.4 A°，水分子內(nèi)的氧－氫鍵距離為1.1 A°。為了編程統(tǒng)計(jì)過(guò)程的方便，從而我們定義氫鍵如下：如圖1所示，當(dāng)圖1中∠ α≤30°并且R≤3.5 A°時(shí)氫鍵形成，當(dāng)∠ α>30°，或 R>3.5 A°時(shí)氫鍵斷裂，這一標(biāo)準(zhǔn)結(jié)合了本文體系的實(shí)際情況，同時(shí)也被眾多國(guó)內(nèi)外學(xué)者所采用[20,23]。

圖1 水中氫鍵的定義Fig.1 The definition of H-bonds in water box

按照上述標(biāo)準(zhǔn)，本文采用VMD軟件并結(jié)合Tcl/Tk編程，統(tǒng)計(jì)了不同溫度、壓力下水分子的成鍵度，結(jié)果如表5所示。從表5可以看出，298 K，0.1 MPa時(shí)，成鍵度為2.73，與Mas[24]用勢(shì)函數(shù)SAPT-5s的計(jì)算值2.77很接近，但與Soper[25]測(cè)得該條件下的實(shí)驗(yàn)值3.58相比偏小，其主要原因還是現(xiàn)有水勢(shì)能模型的準(zhǔn)確性不足引起的。雖然結(jié)果與實(shí)驗(yàn)值有一定的偏差，但用于氫鍵估算和了解氫鍵變化規(guī)律是可行的。不論在何種條件下，隨著溫度的升高，壓力減小，水的平均成鍵度減?。环粗?，成鍵度增加?？梢钥闯?，在超臨界條件下，即使到了873 K，50 MPa下，水中仍然有少量的氫鍵存在，Gorbaty[26]也用實(shí)驗(yàn)證明了在超臨界條件下仍有氫鍵的存在。現(xiàn)今，對(duì)于具體在什么情況下水中氫鍵會(huì)消失仍沒(méi)有定論。

通常，每個(gè)水子可以形成1到5個(gè)氫鍵[27]。本文對(duì)不同條件下水分子各種氫鍵的概率做了統(tǒng)計(jì)，如表6所示。

表5 水中氫鍵的平均成鍵度Table5 The averagenumberofH-bonds perwatermolecule

從表6可以看出，在常態(tài)下，未形成氫鍵的水分子所占的比例很?。凰肿又饕纬闪硕了膫€(gè)氫鍵的結(jié)構(gòu)，占了70%以上。隨著溫度的升高，形成三、四、五個(gè)氫鍵的水分子數(shù)逐漸減少，而未成鍵和形成一、二個(gè)氫鍵的水分子數(shù)所占比例增加。與常態(tài)下相比，亞臨界下未形成氫鍵的水分子數(shù)所占比例增幅明顯，在接近超臨界(633 K,22.05 MPa)時(shí)達(dá)到了33.05%。在亞臨界相同壓力下，溫度升高使形成二、三、四、五個(gè)氫鍵的水分子數(shù)都減少，而未成鍵和形成一個(gè)氫鍵的水分子數(shù)增加；在亞臨界相同溫度下，隨著壓力的增加，未成鍵和形成一、二個(gè)氫鍵的水分子數(shù)有所減小，同時(shí)形成三、四、五個(gè)氫鍵的水分子數(shù)有所增加。到了超臨界條件，形成5個(gè)氫鍵的水分子消失，與亞臨界條件下相比未形成氫鍵的水分子數(shù)大幅增加，在873 K，50 MPa時(shí)達(dá)到了73%。在超臨界相同壓力下，溫度升高使得形成各種類型氫鍵的水分子數(shù)都減少；而相同溫度下，壓力升高未成鍵水分子數(shù)減少，形成一、二、三、四個(gè)氫鍵的水分子數(shù)都有所增加。

表6 不同條件下各種氫鍵的概率Table 6 The probability of the different kinds of H-bonds in water

3.4 水表面張力的計(jì)算

表面張力通過(guò)壓力張量的分量計(jì)算得到[28,29]，如式3：Pxx，Pyy，Pzz分別是系統(tǒng)壓力張量在各個(gè)軸的分量，bar，Lz是模擬體系Z軸的長(zhǎng)度，nm。通過(guò)單位換算即可得到水的表面張力。計(jì)算得到不同溫度下水的表面張力如圖2。圖2中實(shí)線為Vargaftik[30]等人的實(shí)驗(yàn)值。從圖2可以看出，模擬值與實(shí)驗(yàn)值相比都呈偏小的趨勢(shì)，其原因可能是在計(jì)算過(guò)程中的簡(jiǎn)化計(jì)算，沒(méi)有考慮氫鍵等因素對(duì)表面張力的影響，而水中卻含有大量的氫鍵，從而使得模擬值比實(shí)驗(yàn)值偏小。

圖2 水的表面張力Fig.2 The surface tension ofwater

4 結(jié)論

本文采用分子動(dòng)力學(xué)方法對(duì)水的物性進(jìn)行了模擬研究，其中水密度與自擴(kuò)散系數(shù)的模擬值與文獻(xiàn)值基本一致，可見(jiàn)采用分子動(dòng)力學(xué)模擬方法來(lái)模擬預(yù)測(cè)水的性質(zhì)是可行的。尤其當(dāng)在超臨界條件下實(shí)驗(yàn)數(shù)據(jù)難以測(cè)定時(shí)，分子動(dòng)力學(xué)模擬的優(yōu)勢(shì)就更為明顯。模擬結(jié)果表明：超臨界下水的密度隨溫度壓力的改變而變化顯著，這一特性是超臨界水用于萃取、分離有機(jī)溶劑的重要依據(jù)；水在超臨界下比常態(tài)下有著更高的擴(kuò)散能力。對(duì)水中的氫鍵計(jì)算結(jié)果表明：在常態(tài)下水中每個(gè)水分子平均形成的氫鍵數(shù)以二至四個(gè)氫鍵為主，占?xì)滏I總數(shù)的70%以上；然而隨著溫度的升高，每個(gè)水分子形成的氫鍵數(shù)減少。水表面張力模擬值比實(shí)驗(yàn)值小，其主要原因是計(jì)算過(guò)程未考慮氫鍵等因素對(duì)表面張力的影響。

[1]任華.分子模擬在界面相互作用計(jì)算中的應(yīng)用 [D].西安：西北工業(yè)大學(xué)，2007.

[2]W.L.Jorgensen,J.Chandrasekhar,J.D.Madura,R.W.Impey,M.L.Klein.Comparison of simple potential functions for simulating liquid water[J].J.Chem.Phys.,1983,79,926-935.

[3]H.J.C.Berendsen,J.R.Grigera,and T.P.Straatsma.The Missing Term in Effective Pair Potentials[J].J.Phys.Chem.,1987,91:6269-6271.

[4]B.Hess,C.Kutzner,D.van der Spoel.GROMACS 4:Algorithms for Highly Efficient,Load-Balanced,and Scalable Molecular Simulation[J].J.Chem.Theory Comput.,2008,4:435-447.

[5]Berendsen,H.J.C.,Postma,J.P.M.,DiNola,A.,Haak,J.R.Molecular dynamicswith coupling to an external bath [J].J.Chem.Phys.1984,81:3684-3690.

[6]Bussi,G.,Donadio,D.,Parrinello,M.Canonical sampling through velocity rescaling [J].J.Chem.Phys.2007,126:014101.

[7]Essmann,U.,L.Perera,M.L.Berkowitz,T.Darden,H.Lee,and L.Pedersen.A smooth particlemesh ewald potential[J].J.Chem.Phys.,1995,103:8577-8592.

[8]T.Darden,D.York,L.Pedersen.Particle Mesh Ewald:An N-log(N)method for Ewald sums in large systems[J].J.Chem.Phys.,1993,98:10089-10092.

[9]Hess,B.,H.Bekker,H.Berendsen,and J.Fraaije.LINCS:A Linear Constraint Solver for molecular simulations [J].J.Comp.Chem.,1997.18:1463-1472.

[10]馬沛生.化工熱力學(xué) [M].北京:化學(xué)工業(yè)出版社，2005.

[11]A.Saul,W.Wagner.A fundamental equation for water covering the range from the melting line to 1273 K at pressures up to 25000 MPa[J].J.Phys.Chem.Ref.Data,1989,18(4):1537-1564.

[12]楊潤(rùn)田，王志鋒，張海峰，等.超臨界水的應(yīng)用及其環(huán)境中材料腐蝕的研究 [J].表面技術(shù)，2007，36(5):84-87.

[13]肖吉，陸九芳，陳健，李以圭.超臨界水中氣體擴(kuò)散系數(shù)的分子動(dòng)力學(xué)模擬 [J].高?；瘜W(xué)工程學(xué)報(bào)，2001，15(1):6-10.

[14]H.Inomata,S.Saito,P.G.Debenedetti.Molecular Dynamics Simulation of Infinitely Dilute Solutions of Benzene in Supercritical CO2 [J].Fluid Phase Equilibria,1996,116:282-288.

[15]石劍.超臨界CO2體系擴(kuò)散系數(shù)的實(shí)驗(yàn)研究和分子動(dòng)力學(xué)模擬 [D].天津:天津大學(xué)，2006.

[16]K.Krynicki, C.D.Green, D.W.Sawyer. Pressure and temperature dependence of self-diffusion in water[J].Faraday Discuss.Chem.Soc.,1978,66:199-208.

[17]W.J.Lamb,G.A.Hoffman,J.Jonas.Self-diffusion in compressed supercritical water[J].J.Chem.Phys.,1981,74(12):6875-6880.

[18]孫煒，黃素逸，王存文，池汝安.超臨界水密度和自擴(kuò)散系數(shù)預(yù)測(cè)的分子動(dòng)力學(xué)模擬 [J].華中科技大學(xué)學(xué)報(bào)，2008，36(5):103-105.

[19]T.Lazaridis,M.Karplus.Orientational Correlations and Entropy in Liquid Water[J].J.Chem.Phys.,1996,105(10)：4294-4316.

[20]F.W.Starr,J.K.Nielsen,H.E.Stanley.Fast and Slow Dynamics of Hydrogen Bonds in Liquid Water [J].Phys.Rev.Lett.,1999,82(11):2294-2297.

[21]A.G.Kalinichev,S.V.Churakov.Size and topology of molecular clusters in supercritical water:a molecular dynamics simulation [J].Chemical Physics Letters,1999,302:411-417.

[22]J.Teixeira,M.C.Bellissent-Funel,S.H.Chen,Dynamics of water studied by neutron scattering [J].J.Phys.:Condens.Matter,1990,2:105-108.

[23]A.Luzar,D.Chandler.Hydrogen-bond kinetics in liquid water[J].Nature,1996,379(4):55-57.

[24]E.M.Mas,R.Bukowski,K.Szalewicz.Ab initio three-body interactions for water.II.Effects on structure and energetics of liquid[J].J.Chem.Phys.,2003,118(10):4404-4413.

[25]A.K.Soper,F.Bruni,M.A.Ricci.Site-site pair correlation functions of water from 25 to 400°C:Revised analysis of new and old diffraction data [J].J.Chem.Phys.,1997,106(1):247-254.

[26]Yu.E.Gorbaty, A.G.Kalinichev. Hydrogen Bonding in Supercritical Water.1.Experimental Results [J].J.Phys.Chem.,1995,99(15):5336-5340.

[27]張秀欣.液態(tài)水分子間氫鍵相互作用的從頭分子動(dòng)力學(xué)研究 [D].天津：河北工業(yè)大學(xué)，2007.

[28]C.A.Croxton,Statistical Mechanics of the Liquid Surface[M].John Wiley&Sons,Chichester,1980.

[29]N.J.P.Nijmeijer,A.F.Bakker,C.Bruin,et al.A molecular dynamics simulation of the Lennard-Jones liquid-vapor interface[J].J.Chem.Phys.,1988,89:3789-3792.

[30]N.B.Vargaftik,B.N.Volkov,L.D.Voljak.International Tables of the Surface Tension of Water[J].J.Phys.Chem.Ref.Data,1983,12(3):817-820.

Molecular Dynamics Simulation on Physical Properties of Water

WU Fangdi1ZHENG Huidong2WU Yanxiang2

(1.Department of Environment and Architecture Engineering，Wuyi University，Wuyishan，F(xiàn)ujian 354300;2.College of Chemistry and Chemical Engineering,Fuzhou University,Fuzhou，F(xiàn)ujian 350108)

The physical properties of water under different conditions were studied with SPC/Emodel by MD simulation.The density of supercritical water changed more significantly with the variation of temperature and pressure than that under subcritical and normal conditions,which is an important basis for supercritical water being used for extraction and separation of organic solvents.Moreover，supercritical water has higher diffusion ability than normal water.The studies showed that each water molecules mainly formed by two to four hydrogen bonds under normal condition.The number of hydrogen bonds in water reduced as the temperature increased.The simulation value of surface tension of water was smaller than experimental value for do not take into account the effects of H-bonds on the surface tension.

water;density;diffusion coefficient;degree of H-bond;surface tension;MD simulation

Q7

A

1674－2109(2011)05－0061－06

2011-09-07

吳方棣（1985-），男，漢族，助教，碩士，主要研究方向：傳質(zhì)與分離。