不同溫度下水的分子動(dòng)力學(xué)模擬

1、收修改稿日期:99-05-09江蘇省自然科學(xué)基金(B K971241資助項(xiàng)目聯(lián)系人:陸小華1第一作者,周健,男,26歲,博士研究論文不同溫度下水的分子動(dòng)力學(xué)模擬周健陸小華王延儒時(shí)鈞(南京化工大學(xué)化學(xué)工程系南京210009摘要采用分子動(dòng)力學(xué)(MD 模擬的方法對(duì)不同溫度下水的氫鍵結(jié)構(gòu)及擴(kuò)散性質(zhì)進(jìn)行了研究。模擬發(fā)現(xiàn)隨溫度的升高,水分子平均氫鍵配位數(shù)減小,氫鍵鍵角分布變寬,鍵角分布峰值減小,即水分子間氫鍵作用不斷減弱。擴(kuò)散系數(shù)隨溫度的增加而上升,液體水的擴(kuò)散系數(shù)與溫度之間的關(guān)系符合Arrhenius 行為。關(guān)鍵詞水,氫鍵,擴(kuò)散系數(shù),分子動(dòng)力學(xué)(MD ,分子模擬中圖號(hào)O 645116MOLECU LAR

2、 DY NAMICS SIMU LATIONOF WATER AT DIFFERENT TEMPERATURESZHOU Jian LU X iao 2Hua WANG Y an 2Ru SHI J un(Depart ment of Chemical Engineering ,N anjing U niversity of Chemical Technology ,Nanjing 210009ABSTRACT The hydrogen bond structure and diffusion behavior of water at different tem peratures were

3、investigated by MD simulation.With the increase of tem perature ,the average hydrogen bond coordination number of water decreases while the distribution of hydrogen bond angles broadens and the peak value of hydrogen bond angles reduces ,i. e.,the hydrogen bond interaction between water molecules de

4、creases.The relationshi p between the diffusion coefficients of water and temperature is in good accordance with Arrhenius behavior.KE YWOR DS Water Hydrogen bond ,Diffusion coefficient ,Molecular dynamics (MD ,Molecular simulation在化學(xué)、生物、地質(zhì)、海洋等許多學(xué)科領(lǐng)域中,都涉及到水溶液?jiǎn)栴}。水還是化工生產(chǎn)過(guò)程中最常用的溶劑以及許多化學(xué)反應(yīng)的介質(zhì)。但對(duì)這一最常用的溶劑

5、的微觀結(jié)構(gòu)及其不同尋常的理化特性仍然不是很清楚,而無(wú)論開(kāi)發(fā)和設(shè)計(jì)低能耗、無(wú)污染的化工分離工藝,還是開(kāi)發(fā)和設(shè)計(jì)高新技術(shù)含量高的新產(chǎn)品,研究水溶液的氫鍵結(jié)構(gòu)正起著越來(lái)越大的作用1。近年來(lái),采用分子動(dòng)力學(xué)(MD 模擬技術(shù)在取得溶液的平衡、結(jié)構(gòu)及動(dòng)態(tài)性質(zhì)等方面正起著越來(lái)越關(guān)鍵的作用2,3。在許多實(shí)際化工生產(chǎn)過(guò)程中,操作溫度范圍可以從零度以下的低溫一直到幾百度的高溫,因此,對(duì)不同溫度下的水進(jìn)行MD 模擬是很有意義的。1位能對(duì)于水分子,其位能模型既可以由量子力學(xué)從頭計(jì)算(Ab initio 獲得,如MCY 模型4;又可以由半經(jīng)驗(yàn)方法關(guān)聯(lián)獲得,如ST25、SPC 6、TIP4P 7等模型。MCY 模型由于其

6、形式過(guò)于復(fù)雜,需要消耗大量的機(jī)時(shí),因此在實(shí)際模擬過(guò)程中研究者大多采用半經(jīng)驗(yàn)?zāi)Ｐ?。本文選取SPC (Simple Point Charge 模型6作為模擬所用位能。其位能形式如下,u ij (r ij =A i B jq A q B r AB +4r OO 12-r OO6(1式中前半部分為靜電作用,后半部分為短程Lennard 2Jones (LJ 作用,S PC 模型假設(shè)只有氧原子與氧原子之間存在LJ 作用。u ij 為分子對(duì)作用位能,A ,B 分別表示原子A ,B 對(duì)分子i ,j 作循環(huán),q A 、q B 為第16卷第4期1999年7月28日計(jì)算機(jī)與應(yīng)用化學(xué)Computers and Ap

7、plied Chemsitry Vol 116,No 14J uly ,1999A、B原子所帶電荷,r AB為原子A,B間距離,r OO為兩個(gè)分子的氧原子間作用距離,為氧原子LJ 作用參數(shù)。r OH為O2H鍵鍵長(zhǎng),HOH為H2O2H鍵角。其參數(shù)值見(jiàn)表1。表1SPC水模型參數(shù)T able1Potential parameters of SPC w ater#/e=116×10-19C模擬中由于位能的截?cái)?能量將在邊界處不連續(xù),會(huì)導(dǎo)致Hamiltonian函數(shù)不守恒。為克服該問(wèn)題,本文采用移位的位能函數(shù)8形式如下,u s ij(r ij=u ij(r ij-u ij(r cr ijr c

8、0r ijr c(2式(2中,u s ij代表上移的位能,r c為位能截?cái)嗑嚯x。對(duì)位能中靜電力的處理,通常采用兩種方法。一為Ewald加和法9,另一為反應(yīng)場(chǎng)法10。由于E2 wald加和法消耗機(jī)時(shí)較多,本文模擬采用較為省時(shí)的反應(yīng)場(chǎng)法。這時(shí),靜電作用位能可寫(xiě)為u AB(r AB=q A q Br AB1+R F-12R F+1r ABr c3(3式(3中,括號(hào)后半部分為反應(yīng)場(chǎng)的作用,R F為環(huán)境介電常數(shù),計(jì)算中R F取為無(wú)窮。模擬過(guò)程中,位能所采用的最后表達(dá)式為u ij(r ij=Ar OO12-Cr OO6Ar12C-Cr6C+AiBjq A q Br AB1+12r ABr c3-32r A

9、Br cr ABr c0r ABr c(42MD模擬的細(xì)節(jié)在對(duì)水進(jìn)行MD模擬時(shí),由于既需要考慮分子間靜電作用的處理,又需要考慮復(fù)雜的剛體動(dòng)力學(xué)方程的求解。在靜電力的處理上,本文采用如前所述的反應(yīng)場(chǎng)10的方法。在運(yùn)動(dòng)方程的求解上,本文在兩體有效能近擬、周期性邊界條件、時(shí)間平均等效于系綜平均等基本假設(shè)之上,通過(guò)求解平動(dòng)運(yùn)動(dòng)方程組和轉(zhuǎn)動(dòng)運(yùn)動(dòng)方程組得到各分子在不同時(shí)刻的位置和速度,再通過(guò)時(shí)間統(tǒng)計(jì)平均的方法,求得宏觀性質(zhì)。模擬所用總分子數(shù)為256。本文所采用的等溫MD模擬基于Nose11所提出的NV T系統(tǒng)的擴(kuò)展系統(tǒng)法。模擬起始構(gòu)型為面心立方晶格,水分子的起始取向?yàn)殡S機(jī)取向,采用立方周期性邊界條件,各分

10、子起始速度按Maxwell分布取樣,隨后對(duì)速度重新標(biāo)定,以確保體系總動(dòng)量為零。模擬過(guò)程中所取的時(shí)間步長(zhǎng)為015fs(1fs=10-15s,這是由于水分子的轉(zhuǎn)動(dòng)慣量較小,會(huì)發(fā)生高速旋轉(zhuǎn),因此必須采用較小的時(shí)間步以保證運(yùn)動(dòng)方程求解的穩(wěn)定性。模擬過(guò)程中位能截?cái)嗖捎们蛐胃钊シ?截?cái)喟霃綖?19nm。截?cái)嗑嚯x之外分子間相互作用能按平均密度近擬的方法進(jìn)行校正8。平動(dòng)數(shù)值積分采用4階G ear預(yù)測(cè)校正法12,轉(zhuǎn)動(dòng)運(yùn)動(dòng)方程解法采用QUAR TERN ION法13,數(shù)值積分采用5階G ear預(yù)測(cè)校正法12。每次模擬的總時(shí)間為50ps=(1ps=10-12s,其中前25ps使體系達(dá)到平衡,后25ps統(tǒng)計(jì)各種平衡性

11、質(zhì)、結(jié)構(gòu)性質(zhì)和輸運(yùn)性質(zhì)以供分析之用。3模擬結(jié)果和討論本文分別在193115K、253115K、313115K、373115K、433115K、493115K下對(duì)水進(jìn)行了MD模擬。193115K、253115K時(shí)的水為假想的過(guò)冷水,密度仍取為110gcm-3。313115K時(shí)水的密度為019922gcm-3。373115K、433115K、493115K時(shí)的水為飽和液體水,它們的密度分別為019583, 019074,018402gcm-3。模擬結(jié)果如下。311結(jié)構(gòu)性質(zhì)徑向分布函數(shù)是反映流體微觀結(jié)構(gòu)特征的物理量。圖1(af給出了不同溫度下水分子中各原子對(duì)242計(jì)算機(jī)與應(yīng)用化學(xué)16卷徑向分布函數(shù)。

12、從圖1(a f 中可以看出,較低溫度下,氧氧徑向分布函數(shù)分別在0128nm 、0145nm 和0170nm 處出現(xiàn)峰值,這分別表示中心水分子與最近鄰水分子、第二配位圈水分子和第三配位圈水分子間的氧氧距離,隨溫度的升高,從193115K 到493115K ,氧氧徑向分布函數(shù)第一峰峰值的變化為4102,3122,2163,2135,2121,2108,其值不斷減小,而第一峰峰谷的變化為0154,0174,0192,0193,0194,0195,其值不斷升高,同時(shí)氧氧第一峰的位置略右移,而氧氧第二峰、第三峰則逐漸消失,這表明隨溫度的升高,水的有序程度逐漸下降。較低溫度下,氧氫徑向分布函數(shù)在0118n

13、m 和0133nm 處出現(xiàn)峰值,這分別是有氫鍵作用和無(wú)氫鍵作用的OH 距離。隨溫度的升高,有氫鍵作用的氧氫第一峰峰值的變化為2116,1167,1128,1128,1111,0196,0186,其值不斷下降,同時(shí)氧氫第一峰峰谷值的變化為0108,0116,0126,0133,014,0145,其值不斷上升,氧氫第二峰峰值的變化為1172,1162,1150,1145,1136,1130,其值不斷下降,這表明溫度的升高,水分子間氫鍵作用不斷減弱。表2不同溫度下液體水的氫鍵配位數(shù)分布表T able 2T emperature dependence of hydrogencoordination n

14、umber of w ater溫度(K 3192317831583139311221 85圖1不同溫度下液體水的徑向分布函數(shù)Fig.1Radial distribution functions of water at different temperatures312氫健分析筆者采用贗組分分布函數(shù)(Quasi ComponentDistribution Function 14對(duì)水分子之間的氫鍵作用進(jìn)行了分析。一對(duì)水分子若要形成氫鍵,通常其相互作用能要小于某一能量界限,OH 距離要小于某一位置界限。本文參照文獻(xiàn)7以能量界限E HB =-814kJ /mol ,位置界限r(nóng) OH =0124nm 作

15、為水分子間形成氫鍵的依據(jù)。表2為不同溫度下液體水的氫鍵配位數(shù)分布表。從表2中,可以看出在各個(gè)溫度下,四配位的水分子占絕對(duì)多數(shù),但隨著溫度的升高,四配位所占的比例逐漸減少,而三配位的水分子所占比例則逐漸上升。從總的加權(quán)平均氫鍵配位數(shù)來(lái)看,隨溫度的升高,平均氫鍵配位數(shù)從3192降到2185。3424期周健等:不同溫度下水的分子動(dòng)力學(xué)模擬氫鍵鍵角的定義為A 分子中OH 聯(lián)線與A 分子中H 和B 分子中O 聯(lián)線的夾角。圖2為不同溫度下氫鍵鍵角分布。從圖中可以看出,隨溫度的升高,氫鍵鍵角分布變寬,鍵角分布峰值減小。這也表明隨溫度的升高,氫鍵作用越來(lái)越弱。313擴(kuò)散性質(zhì)擴(kuò)散系數(shù)的求取由Einstein 法

16、8在模擬過(guò)程中不斷統(tǒng)計(jì)獲得。水的自擴(kuò)散系數(shù)由下式統(tǒng)計(jì)獲得D s =limt 16N tNi =1r i(t -r i (02。溫度對(duì)液體的擴(kuò)散性質(zhì)的影響通常符合Arrhenius 行為,即有D =D 0exp - E aR T,其中D 0為指前因子,E a 為擴(kuò)散活化能。筆者將擴(kuò)散系數(shù)模擬數(shù)據(jù)按上述公式進(jìn)行擬和得到,D 0=4132×10-3cm 2s -1,E a =12 102kJ /mol 。圖2不同溫度下水的氫鍵鍵角分布Fig.2Distribution of hydrogen bond angle of liquid water at different temperatu

17、res圖3水的擴(kuò)散系數(shù)與溫度的關(guān)系Fig.3Relationship between diffusion coeffcients of water and temperature圖3為溫度的倒數(shù)與液體水的擴(kuò)散系數(shù)的對(duì)數(shù)之間的關(guān)系圖以及模擬數(shù)據(jù)與實(shí)驗(yàn)數(shù)據(jù)15的比較,從圖中可以看出其線性關(guān)系良好,符合Arrhenius 行為。與實(shí)驗(yàn)值相比,本文由SPC 模型模擬得到的擴(kuò)散系數(shù)結(jié)果與實(shí)驗(yàn)值吻合基本良好,只是在低溫下略有差異,這是由于SPC 水模型參數(shù)是擬合常溫時(shí)液體水的性質(zhì)而得到的,用于推測(cè)其它條件下的物性會(huì)有所差異,但我們可以看出擴(kuò)散系數(shù)其隨溫度的變化關(guān)系還是一致的。4結(jié)論模擬得到了不同狀態(tài)條件下

18、水的新的模擬數(shù)據(jù)。隨溫度的升高,水分子平均氫鍵配位數(shù)減小,氫鍵鍵角分布變寬,鍵角分布峰值減小,即水分子間氫鍵作用不斷減弱。溫度與液體水的擴(kuò)散系數(shù)之間的關(guān)系符合Arrhenius 行為。參考文獻(xiàn)1房春暉.化學(xué)進(jìn)展,1996,8:3182周健,陸小華,王延儒,時(shí)鈞.南京化工大學(xué)學(xué)報(bào),1998,20(3:13周健.流體微觀結(jié)構(gòu)及擴(kuò)散性質(zhì)的分子動(dòng)力學(xué)模擬研究:博士學(xué)位論文.南京:南京化工大學(xué),19984Matsuoka O ,Clementi E ,Y oshimine M.J Chem Phys ,1976,64:13515Stillinger F H ,Rahman A.J Chem Phys ,1974,60:15456Berendsen H J C ,Postma J P M ,Gunsteren