不同温度下水的分子动力学模拟

1、收修改稿日期:99-05-09江苏省自然科学基金(B K971241资助项目联系人:陆小华1第一作者,周健,男,26岁,博士研究论文不同温度下水的分子动力学模拟周健陆小华王延儒时钧(南京化工大学化学工程系南京210009摘要采用分子动力学(MD 模拟的方法对不同温度下水的氢键结构及扩散性质进行了研究。模拟发现随温度的升高,水分子平均氢键配位数减小,氢键键角分布变宽,键角分布峰值减小,即水分子间氢键作用不断减弱。扩散系数随温度的增加而上升,液体水的扩散系数与温度之间的关系符合Arrhenius 行为。关键词水,氢键,扩散系数,分子动力学(MD ,分子模拟中图号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在化学、生物、地质、海洋等许多学科领域中,都涉及到水溶液问题。水还是化工生产过程中最常用的溶剂以及许多化学反应的介质。但对这一最常用的溶剂

5、的微观结构及其不同寻常的理化特性仍然不是很清楚,而无论开发和设计低能耗、无污染的化工分离工艺,还是开发和设计高新技术含量高的新产品,研究水溶液的氢键结构正起着越来越大的作用1。近年来,采用分子动力学(MD 模拟技术在取得溶液的平衡、结构及动态性质等方面正起着越来越关键的作用2,3。在许多实际化工生产过程中,操作温度范围可以从零度以下的低温一直到几百度的高温,因此,对不同温度下的水进行MD 模拟是很有意义的。1位能对于水分子,其位能模型既可以由量子力学从头计算(Ab initio 获得,如MCY 模型4;又可以由半经验方法关联获得,如ST25、SPC 6、TIP4P 7等模型。MCY 模型由于其

6、形式过于复杂,需要消耗大量的机时,因此在实际模拟过程中研究者大多采用半经验模型。本文选取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 模型假设只有氧原子与氧原子之间存在LJ 作用。u ij 为分子对作用位能,A ,B 分别表示原子A ,B 对分子i ,j 作循环,q A 、q B 为第16卷第4期1999年7月28日计算机与应用化学Computers and Ap

7、plied Chemsitry Vol 116,No 14J uly ,1999A、B原子所带电荷,r AB为原子A,B间距离,r OO为两个分子的氧原子间作用距离,为氧原子LJ 作用参数。r OH为O2H键键长,HOH为H2O2H键角。其参数值见表1。表1SPC水模型参数T able1Potential parameters of SPC w ater#/e=116×10-19C模拟中由于位能的截断,能量将在边界处不连续,会导致Hamiltonian函数不守恒。为克服该问题,本文采用移位的位能函数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为位能截断距离。对位能中静电力的处理,通常采用两种方法。一为Ewald加和法9,另一为反应场法10。由于E2 wald加和法消耗机时较多,本文模拟采用较为省时的反应场法。这时,静电作用位能可写为u AB(r AB=q A q Br AB1+R F-12R F+1r ABr c3(3式(3中,括号后半部分为反应场的作用,R F为环境介电常数,计算中R F取为无穷。模拟过程中,位能所采用的最后表达式为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模拟的细节在对水进行MD模拟时,由于既需要考虑分子间静电作用的处理,又需要考虑复杂的刚体动力学方程的求解。在静电力的处理上,本文采用如前所述的反应场10的方法。在运动方程的求解上,本文在两体有效能近拟、周期性边界条件、时间平均等效于系综平均等基本假设之上,通过求解平动运动方程组和转动运动方程组得到各分子在不同时刻的位置和速度,再通过时间统计平均的方法,求得宏观性质。模拟所用总分子数为256。本文所采用的等温MD模拟基于Nose11所提出的NV T系统的扩展系统法。模拟起始构型为面心立方晶格,水分子的起始取向为随机取向,采用立方周期性边界条件,各分

10、子起始速度按Maxwell分布取样,随后对速度重新标定,以确保体系总动量为零。模拟过程中所取的时间步长为015fs(1fs=10-15s,这是由于水分子的转动惯量较小,会发生高速旋转,因此必须采用较小的时间步以保证运动方程求解的稳定性。模拟过程中位能截断采用球形割去法,截断半径为019nm。截断距离之外分子间相互作用能按平均密度近拟的方法进行校正8。平动数值积分采用4阶G ear预测校正法12,转动运动方程解法采用QUAR TERN ION法13,数值积分采用5阶G ear预测校正法12。每次模拟的总时间为50ps=(1ps=10-12s,其中前25ps使体系达到平衡,后25ps统计各种平衡性

11、质、结构性质和输运性质以供分析之用。3模拟结果和讨论本文分别在193115K、253115K、313115K、373115K、433115K、493115K下对水进行了MD模拟。193115K、253115K时的水为假想的过冷水,密度仍取为110gcm-3。313115K时水的密度为019922gcm-3。373115K、433115K、493115K时的水为饱和液体水,它们的密度分别为019583, 019074,018402gcm-3。模拟结果如下。311结构性质径向分布函数是反映流体微观结构特征的物理量。图1(af给出了不同温度下水分子中各原子对242计算机与应用化学16卷径向分布函数。

12、从图1(a f 中可以看出,较低温度下,氧氧径向分布函数分别在0128nm 、0145nm 和0170nm 处出现峰值,这分别表示中心水分子与最近邻水分子、第二配位圈水分子和第三配位圈水分子间的氧氧距离,随温度的升高,从193115K 到493115K ,氧氧径向分布函数第一峰峰值的变化为4102,3122,2163,2135,2121,2108,其值不断减小,而第一峰峰谷的变化为0154,0174,0192,0193,0194,0195,其值不断升高,同时氧氧第一峰的位置略右移,而氧氧第二峰、第三峰则逐渐消失,这表明随温度的升高,水的有序程度逐渐下降。较低温度下,氧氢径向分布函数在0118n

13、m 和0133nm 处出现峰值,这分别是有氢键作用和无氢键作用的OH 距离。随温度的升高,有氢键作用的氧氢第一峰峰值的变化为2116,1167,1128,1128,1111,0196,0186,其值不断下降,同时氧氢第一峰峰谷值的变化为0108,0116,0126,0133,014,0145,其值不断上升,氧氢第二峰峰值的变化为1172,1162,1150,1145,1136,1130,其值不断下降,这表明温度的升高,水分子间氢键作用不断减弱。表2不同温度下液体水的氢键配位数分布表T able 2T emperature dependence of hydrogencoordination n

14、umber of w ater温度(K 3192317831583139311221 85图1不同温度下液体水的径向分布函数Fig.1Radial distribution functions of water at different temperatures312氢健分析笔者采用赝组分分布函数(Quasi ComponentDistribution Function 14对水分子之间的氢键作用进行了分析。一对水分子若要形成氢键,通常其相互作用能要小于某一能量界限,OH 距离要小于某一位置界限。本文参照文献7以能量界限E HB =-814kJ /mol ,位置界限r OH =0124nm 作

15、为水分子间形成氢键的依据。表2为不同温度下液体水的氢键配位数分布表。从表2中,可以看出在各个温度下,四配位的水分子占绝对多数,但随着温度的升高,四配位所占的比例逐渐减少,而三配位的水分子所占比例则逐渐上升。从总的加权平均氢键配位数来看,随温度的升高,平均氢键配位数从3192降到2185。3424期周健等:不同温度下水的分子动力学模拟氢键键角的定义为A 分子中OH 联线与A 分子中H 和B 分子中O 联线的夹角。图2为不同温度下氢键键角分布。从图中可以看出,随温度的升高,氢键键角分布变宽,键角分布峰值减小。这也表明随温度的升高,氢键作用越来越弱。313扩散性质扩散系数的求取由Einstein 法

16、8在模拟过程中不断统计获得。水的自扩散系数由下式统计获得D s =limt 16N tNi =1r i(t -r i (02。温度对液体的扩散性质的影响通常符合Arrhenius 行为,即有D =D 0exp - E aR T,其中D 0为指前因子,E a 为扩散活化能。笔者将扩散系数模拟数据按上述公式进行拟和得到,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水的扩散系数与温度的关系Fig.3Relationship between diffusion coeffcients of water and temperature图3为温度的倒数与液体水的扩散系数的对数之间的关系图以及模拟数据与实验数据15的比较,从图中可以看出其线性关系良好,符合Arrhenius 行为。与实验值相比,本文由SPC 模型模拟得到的扩散系数结果与实验值吻合基本良好,只是在低温下略有差异,这是由于SPC 水模型参数是拟合常温时液体水的性质而得到的,用于推测其它条件下的物性会有所差异,但我们可以看出扩散系数其随温度的变化关系还是一致的。4结论模拟得到了不同状态条件下

18、水的新的模拟数据。随温度的升高,水分子平均氢键配位数减小,氢键键角分布变宽,键角分布峰值减小,即水分子间氢键作用不断减弱。温度与液体水的扩散系数之间的关系符合Arrhenius 行为。参考文献1房春晖.化学进展,1996,8:3182周健,陆小华,王延儒,时钧.南京化工大学学报,1998,20(3:13周健.流体微观结构及扩散性质的分子动力学模拟研究:博士学位论文.南京:南京化工大学,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