《不同温度下水的分子动力学模拟》由会员分享，可在线阅读，更多相关《不同温度下水的分子动力学模拟（4页珍藏版）》请在金锄头文库上搜索。
1、收修改稿日期: 99 - 05 - 09江苏省自然科学基金(BK971241)资助项目联系人:陆小华 1 第一作者,周健,男, 26岁,博士 研究论文不同温度下水的分子动力学模拟周健 陆小华 王延儒 时钧(南京化工大学化学工程系 南京 210009)摘 要 采用分子动力学(MD)模拟的方法对不同温度下水的氢键结构及扩散性质进行了研究。模拟发现随温度的升高,水分子平均氢键配位数减小,氢键键角分布变宽,键角分布峰值减小,即水分子间氢键作用不断减弱。扩散系数随温度的增加而上升,液体水的扩散系数与温度之间的关系符合Arrhenius行为。关键词 水,氢键,扩散系数,分子动力学(MD) ,分子模拟中图号
2、 O 645116MOLECULAR DYNAMICS SIMULATION OF WATER AT DIFFERENT TEMPERATURESZHOU Jian LU Xiao2Hua WANG Yan2Ru SHIJun( Department of Chemical Engineering , Nanjing University of Chemical Technology ,Nanjing 210009)ABSTRACT The hydrogen bond structure and diffusion behavior of water at different tem pera
3、tures were investigated byMD simulation. With the increase of temperature , the average hydrogen bond coordination number of water decreaseswhile the distribution of hydrogen bond angles broadens and the peak value of hydrogen bond angles reduces , i. e. , thehydrogen bond interaction between water
4、molecules decreases. The relationshi p between the diffusion coefficients of waterand temperature is in good accordance with Arrhenius behavior.KEYWORDS Water Hydrogen bond , Diffusion coefficient , Molecular dynamics (MD) , Molecular simulation在化学、生物、地质、海洋等许多学科领域中,都涉及到水溶液问题。水还是化工生产过程中最常用 的溶剂以及许多化学反
5、应的介质。但对这一最常用的溶剂的微观结构及其不同寻常的理化特性仍然不是很清楚,而无论开发和设计低能耗、无污染的化工分离工艺,还是开发和设计高新技术含量高的新产品,研 究水溶液的氢键结构正起着越来越大的作用1。近年来,采用分子动力学(MD)模拟技术在取得溶液的 平衡、结构及动态性质等方面正起着越来越关键的作用2 ,3。在许多实际化工生产过程中,操作温度范围 可以从零度以下的低温一直到几百度的高温,因此,对不同温度下的水进行MD模拟是很有意义的。1 位能对于水分子,其位能模型既可以由量子力学从头计算(Ab initio)获得,如MCY模型4;又可以由 半经验方法关联获得,如ST25、SPC6、TI
6、P4P7等模型。MCY模型由于其形式过于复杂,需要消耗大量的机时,因此在实际模拟过程中研究者大多采用半经验模型。 本文选取SPC (Simple Point Charge)模型6作为模拟所用位能。其位能形式如下,uij( rij) = Ai BjqAqB rAB+4 rOO12 - rOO6 (1)式中前半部分为静电作用,后半部分为短程Lennard2Jones (LJ)作用,S PC模型假设只有氧原子与氧原 子之间存在LJ作用。uij为分子对作用位能,A ,B分别表示原子A ,B对分子i ,j作循环,qA、qB为第16卷 第4期 1999年7月28日计算机与应用化学 Computers an
7、d Applied ChemsitryVol116 , No14 July , 1999A、B原子所带电荷,rAB为原子A ,B间距离,rOO为两个分子的氧原子间作用距离,为氧原子LJ作用参数。rOH为O2H键键长,HOH为H2O2H键角。其参数值见表1。表1 SPC水模型参数Table 1 Potential parameters of SPC waterrOH(nm)HOH ()/ k( K)(nm)q ( O)( e)#q ( H)( e)SPC0111091477811820131656- 01820141# / e = 11610- 19C模拟中由于位能的截断,能量将在边界处不连续,
8、会导致Hamiltonian函数不守恒。为克服该问题,本文采用移位的位能函数8形式如下,usij( rij) =uij( rij) -uij( rc)rijrc0rijrc(2)式(2)中,usij代表上移的位能,rc为位能截断距离。对位能中静电力的处理,通常采用两种方法。一为Ewald加和法9,另一为反应场法10。由于E2wald加和法消耗机时较多,本文模拟采用较为省时的反应场法。这时,静电作用位能可写为uAB( rAB) =qAqB rAB1+RF-1 2RF+1rAB rc3 (3)式(3)中,括号后半部分为反应场的作用,RF为环境介电常数,计算中RF取为无穷。模拟过程 中,位能所采用的
9、最后表达式为uij( rij) =A rOO12-C rOO6A r12C-C r6C+ Ai BjqAqB rAB1+1 2rAB rc3 -3 2rAB rcrABrc0rABrc(4)2 MD模拟的细节在对水进行MD模拟时,由于既需要考虑分子间静电作用的处理,又需要考虑复杂的刚体动力学方程的求解。在静电力的处理上,本文采用如前所述的反应场10的方法。在运动方程的求解上,本文在两体有效能近拟、周期性边界条件、时间平均等效于系综平均等基本假设之上,通过求解平动运动方程组和 转动运动方程组得到各分子在不同时刻的位置和速度,再通过时间统计平均的方法,求得宏观性质。模拟所用总分子数为256。本文所
10、采用的等温MD模拟基于Nose11所提出的NVT系统的扩展系统法。模拟起始构型为面心立方晶格,水分子的起始取向为随机取向,采用立方周期性边界条件,各分子起 始速度按Maxwell分布取样,随后对速度重新标定,以确保体系总动量为零。模拟过程中所取的时间步长为015fs (1 fs = 10- 15s),这是由于水分子的转动惯量较小,会发生高速旋转,因此必须采用较小的时间 步以保证运动方程求解的稳定性。模拟过程中位能截断采用球形割去法,截断半径为019nm。截断距离之外分子间相互作用能按平均密度近拟的方法进行校正8。平动数值积分采用4阶Gear预测 校正法12,转动运动方程解法采用QUARTERN
11、ION法13,数值积分采用5阶Gear预测 校正法12。每次模拟的总时间为50ps =(1ps = 10- 12s) ,其中前25ps使体系达到平衡,后25ps统计各种平衡性质、结构性质和 输运性质以供分析之用。3 模拟结果和讨论本文分别在193115K、253115K、313115K、373115K、433115K、493115K下对水进行了MD模拟。193115K、253115K时的水为假想的过冷水,密度仍取为110 gcm- 3。313115K时水的密度为019922 gcm- 3。373115K、433115K、493115K时的水为饱和液体水,它们的密度分别为019583 ,0190
12、74 , 018402 gcm- 3。模拟结果如下。311 结构性质 径向分布函数是反映流体微观结构特征的物理量。图1(af)给出了不同温度下水分子中各原子对242计算机与应用化学16卷径向分布函数。从图1(af)中可以看出,较低温度下,氧氧径向分布函数分别在0128nm、0145nm和0170nm处出现峰值,这分别表示中心水分子与最近邻水分子、第二配位圈水分子和第三配位圈水分子间 的氧氧距离,随温度的升高,从193115K到493115K,氧氧径向分布函数第一峰峰值的变化为4102,3122,2163,2135,2121,2108,其值不断减小,而第一峰峰谷的变化为0154,0174,019
13、2,0193,0194,0195,其值不断升高,同时氧氧第一峰的位置略右移,而氧氧第二峰、第三峰则逐渐消失,这表明随温度的升高,水的有序程度逐渐下降。较低温度下,氧氢径向分布函数在0118 nm和0133 nm处出现 峰值,这分别是有氢键作用和无氢键作用的OH距离。随温度的升高,有氢键作用的氧氢第一峰峰值的变 化为2116,1167,1128,1128,1111,0196,0186,其值不断下降,同时氧氢第一峰峰谷值的变化为0108,0116,0126,0133,014,0145,其值不断上升,氧氢第二峰峰值的变化为1172,1162,1150,1145,1136,1130,其值不断下降,这表
14、明温度的升高,水分子间氢键作用不断减弱。表2 不同温度下液体水的氢键配位数分布表Table 2 Temperature dependence of hydrogencoordination number of water温度( K)1931152531153131153731154331154931151010201160193210341788179211043164817413192201932717831319323158321483614338165371604771186312047166381962911821166571729112917381136103319160111013
15、001450152014201267010001000101010101010100加权平均319231783158313931122185图1 不同温度下液体水的径向分布函数Fig. 1 Radial distribution functions of water at different temperatures312 氢健分析笔者采用赝组分分布函数(Quasi ComponentDistribution Function)14对水分子之间的氢键作用进行了分析。一对水分子若要形成氢键,通常其相互作用能要小于某一能量界限, OH距离要小于某一位置界限。本文参照文献7以能量界限EHB= - 8
16、14kJ/mol ,位置界限rOH= 0124 nm作为水分子间形成氢键的依据。表2为不同温度下液体水的氢键配位数分布表。从表2中,可以看出在各个温度下,四配位的水分子占绝对多数,但随着温度的升高,四配位所占的比例逐渐减少,而三配位的水分子所占比例则逐渐上升。从总的加权平均氢键配位数来看,随温度的升高,平均氢键配位数从3192降到2185。3424期周健等:不同温度下水的分子动力学模拟氢键键角的定义为A分子中OH联线与A分子中H和B分子中O联线的夹角。图2为不同温度下氢 键键角分布。从图中可以看出,随温度的升高,氢键键角分布变宽,键角分布峰值减小。这也表明随温度 的升高,氢键作用越来越弱。313 扩散性质 扩散系数的求取由Einstein法8在模拟过程中不断统计获得。水的自扩散系数由下式统计获得Ds=lim t1 6NtNi =1ri( t) -ri(0)2 。温度对液体
