电场对水的热力学性质影响的分子动力学模拟

电场对水的热力学性质影响的分子动力学模拟

1水的相互作用水是自然界最常见的液体。它也涉及地质、生物、海洋、化学和其他领域的水的应用。这也是化工生产中最常用的溶剂。随着计算机科学技术的快速发展，计算机模拟方法为研究分子水平提供了一种复杂的系统。近年来，我们对水的微观结构和动力学性质进行了大量研究。corden等人研究了不同类型的水下，并使用原始计算模型模拟了水下的结构。另一方面，通过monterey方法研究了不同温度下水的结构变化。ri使用monterey方法研究了不同温度下水的结构变化。周健等人用分子动力学方法研究了液体水的结构性质。孙伟等人用分子动力学方法研究了江底系统中的水密度和扩散系数。对于水及水溶液在外电场作用下的研究在电化学和生物化学等领域具有重要的意义.Zhu等用分子动力学研究了交流外电场下水的结构性质,Svishchev等研究了外电场下水的电致结冰现象,Jun等用环状结构研究了水在外电场作用下的结构变化.孙炜用分子动力学方法研究了温度在258K,电场强度在0~1×1010V/m时液态水的结构和氢键.但是目前对温度及外加电场对水的热力学性质、微观结构及动力学性质系统研究的却不多见.本文即采用分子动力学方法,系统地研究了温度从273~373K,电场强度从0~2×1010V/m时对水分子团簇的热力学性质、微观结构及动力学性质的研究.2刚体刚体干预的主要机理对于水分子,其势能模型既可以从量子力学从头算获得,如MCY模型,又可由半经验方法关联获得,如SPC、SPCE、TIP3P、TIP4P和TIP5P等.由于MCY模型形式过于复杂,又耗费大量的机时,因此在实际模拟过程中许多研究者大多采用半经验的势能模型.本文采用半经验的SPCE刚体势能模型,分子间的相互作用势函数表达形式为式(1)中,前半部分为长程静电相互作用,后半部分为短程Lennard-Jones作用.Uij(rij)为分子对作用势能,A、B分别表示为A、B原子对i、j做循环,qA为A原子所带电荷,qB为B原子所带电荷,rAB为原子间距离,rOO为两个分子的氧原子间作用距离,σ、ε为氧原子Lennard-Jones作用参数.rOH为O-H键键长,θHOH为H-O-H键角.其具体参数见表1.模拟中由于势能的截断,能量会在边界处发生不连续,将导致Hamiltonian不守恒.为了解决这个问题,本文采用移位的势能函数,其形式如下其中,USij(rij)为上移的位能,rC为势能的截断距离.3u3000模拟方法外电场对原子的作用力表示为Fi=qi⋅E,其中qi为原子i所带电荷,E为电场强度.在温度分别为273K、293K、313K、333K、353K、373K,电场强度E分别为0V/m、0.5×1010V/m、1×1010V/m、1.5×1010V/m和2×1010V/m下对水分子团簇进行分子动力学模拟,所选分子数为256个,系综为正则系综(NVT),压强为1atm,密度为1.0g/cm3,模拟起始构型为面心立方晶格,水分子的起始取向为随机取向,各分子的起始速度按MaxwellBoltzmann分布取样.模拟过程中采用两体有效势近似、最小镜像准则等基本假设,采用立方周期性边界条件,对温度的控制采用Nose-Hoover热浴法,模拟盒子大小由密度确定.模拟过程中,长程静电力的处理采用Ewald加和法,短程作用力的势能截断采用球形割去法(sphericalcutoff),截断半径为半个盒子的边长,截断距离之外的粒子间相互作用能采用长程校正法进行校正,采用SHAKE算法固定几何构型,运动方程的求解采用5阶Gear预估-校正法.模拟过程中所取时间步长为1fs,每次模拟总时间为400ps,前200ps使体系达到平衡,后200ps用于进行统计计算各种物理性质.4结果与讨论4.1最大温度和电场强度的模拟结果本文模拟了水分子团簇在温度分别为273K、293、313K、333K、353K、373K,以及电场强度分别为0V/m、0.5×1010V/m、1×1010V/m、1.5×1010V/m和2×1010V/m下的内能.其模拟结果见表2.从表2中看出,同一电场强度下,随温度的升高,内能一直在增加,这是因为温度升高,加剧了分子的热运动以及分子之间的相互碰撞,使内能增加.同一温度下,随电场强度的增加,内能也减小.因此,我们得出,内能与温度成正比与电场强度成反比.4.2氧氢和氧氢径向分布函数中水的峰谷分析径向分布函数(RadialDistributionFunction)gij(r)是与时间无关的量度,它的物理意义是指距离某中心粒子r处出现另一粒子的概率密度与体系的平均密度之比.可表示为式中,ρ(r)为距离中心粒子r处的概率密度,ρ为体系的平均密度.当某一粒子距离中心粒子附近(r值很小)时该粒子出现的概率密度不同于体系的平均密度,但当与中心粒子很远(r值很大)时,该粒子的概率密度趋于体系的平均密度,即r值很大时,径向分布函数的值应接近与1,它反映了流体的微观结构.在不同温度和外加电场下,我们可以通过径向分布函数来了解水的微观结构的变化.本文计算了温度为273K、293K、313K、333K、353K、373K,以及电场强度为0V/m、0.5×1010V/m、1×1010V/m、1.5×1010V/m和2×1010V/m下的径向分布函数鉴于篇幅有限,图1只列出了温度在273K、313K、373K,及电场强度在0V/m、0.5×1010V/m、1×1010V/m、1.5×1010V/m和2×1010V/m下的径向分布函数.从图1中看出同一电场强度下,随温度的升高,氧氧和氧氢径向分布函数第一峰和第二峰的峰位不变,峰值降低,峰谷抬升.对于氧氢径向分布函数,分别在0.17nm、0.32nm出现第一峰和第二峰,这分别是有氢键与无氢键作用的OH距离,随温度的升高,有氢键作用的第一峰峰值不断下降,第一峰谷不断上升,这与中子衍射实验结果吻合,表明温度升高,水分子间的氢键作用不断减弱.以上分析表明,随温度的升高,水的有序程度逐渐下降.同一温度下,随电场强度的增加,氧氧和氧氢径向分布函数的第一峰和第二峰峰值不断上升,峰谷下降峰位其本保持不变,这与Kiselev等的模拟结果相同.这一结果表明,水的第一配位圈结构得到了增强,且水的第二第三配位圈的结构受电场作用的影响更加明显,第三第四峰明显出现,甚至第五第六峰都有可能出现这表明,外加电场使水的有序程度增加,氢键作用增强.但是随温度的升高,第三第四峰的峰值逐渐减小,到373K时,只在E=20V/nm时,峰位比较明显,在0~15V/nm时,几乎消失.对于图1中的(b)和(d)图,我们注意到随电场强度增加,在0.3~0.6nm之间,氧氢径向分布函数出现了峰值的波动,我们称之为次峰,且随电场强度增加而明显,这可能是由于外电场的作用,使水分子的空间取向发生了变化,受电场作用水分子发生旋转改变了氢原子的位置,使得在0.3~0.6nm时,氧氢在此区域出现了小的聚集,但是温度的增加,加剧了分子的碰撞使这种次峰减弱.4.3温度和电场强度对自扩散系数的影响为了了解温度和外电场对水分子团簇的动力学性质的影响,本文计算了不同温度和电场强度下水的自扩散系数.自扩散系数可通过两种等价的方法求得.一种是对速度自相关函数求积分的Green-Kubo法,另一种方法是对均方位移求斜率的Einstein法.Green-Kubo公式为式中,m为水分子的质量,KB为Boltzmann常数,T为体系的温度,vi(0)、vi(t)分别为0时刻和t时刻粒子i的速度Einstein公式为式中,N为模拟体系的分子数,t为时间,ri(t)和ri(0)分别对应粒子在t时刻和0时刻的位置.我们知道,对于速度自相关函数在某一时间很快从1降到0附近,此后一直在0附近振荡,表明这一时间之后,粒子的速度与时间已不再相关,这一时间称为松弛时间.由于松弛时间很难准确确定.本文则采用Einstein法计算不同温度和电场强度下水的自扩散系数.为了说明本文模拟时所选SPCE模型的正确性,表3给出了电场强度为0V/nm时不同温度下水的自扩散系数以及其他文献值.表3可以看出随着温度的升高自扩散系数增加,温度升高,加剧了分子的运动和碰撞,使分子偏离初始位置增大,相互作用力减弱,有利于扩散,因此自扩散系数随温度的增加而增加.文献的扩散系数由其所给的均方位移求得,从表3看出本文的计算值与文献接近,而相比文献用第一性原理分子动力学方法得到的扩散系数要小一些.而且模拟结果的变化趋势与文献所得水的自扩散系数随温度升高而升高的结论是一致的这说明本文选用的模型是合理的.表4给出了不同温度和电场强度下水的自扩散系数.由表4看出,同一温度下,随电场强度的增加,水的自扩散系数降低,表明在外电场作用下,电场对原子产生的作用力即Fi=qi⋅E阻碍了粒子的运动.在低温时,有场时的自扩散系数是无场时的十分之一,这与文献所得的结论一致.综上所述,我们得出结论,温度升高,自扩散系数增加,电场强度增加,自扩散系数较小.5电场强度对氧氢径向分布函数的影响(1)同一电场强度下,随温度的升高,内能增加,温度升高,加剧了分子的热运动以及分子之间的相互碰撞使内能增加.同一温度下,随电场强度的增加,内能也减小.我们得出,水的内能与温度成正比与电场强度成反比.(2)随着温度的升高,水分子间的氢键作用不断减弱,水的有序程度逐渐下降.随电场强度的增加,水的第一配位圈结构得到了增强,第二第三配位圈的结构受电场作用的影响更加明显,第三第四峰明显出现,甚至第五第六峰都有可能出现,这表明,外加电场使水的有序程度增加,氢键作用增强.在0.3~0.6nm之间,氧氢径向分布函数出现了峰值的波动即次峰,且随电场强度增加而明显,这可能是由于外电场的作用,使水分子的空间取向发生了变