静电场中水分子的扩散与液固相变

管星悦，李文飞

(南京大学 物理学院，江苏 南京 210093)

体相水分子系统是人们最为熟悉的分子系统[1].然而，由于分子间的氢键网络结构，水分子系统具有非常丰富的结构与动力学行为，并表现出诸多特殊性质，如高热容、结冰膨胀以及液相密度的非单调温度依赖性等.另外，水分子也是生命体在分子、细胞、组织等多个层次生命运动的关键参与者，因此水分子系统不仅是物理学专业学生学习分子运动与热力学知识的重要实例体系，也是物理、化学、以及生命科学领域的重要研究对象.

体相水的结构与动力学性质一直是人们探讨的重要问题.2006年Laage和Hynes[2]发现水分子存在大角度跳跃转动机制.另外，由于水分子的固有偶极矩，外加静电场可以显著改变水分子系统的结构、动力学与热力学性质.例如，在外加稳恒电场强度大于2.5 V·nm-1时，可以观察到液固相变，即电致结冰[3].Zong等人[4]发现外加静电场能够导致体相水分子系统的黏性产生各向异性.

分子动力学模拟通过求解系统中基本单元的运动方程，从而得到坐标、速度和能量等信息的时间演化，被广泛应用于分子系统结构、动力学与热力学性质的研究.本文基于分子动力学模拟，研究了静电场环境下水分子的大角度跳跃转动运动以及液固相变等性质.结果发现在小电场区域，水分子的大角度跳跃转动运动频次呈现出随静电场强度先增加后减小的非单调依赖性.在大电场区域，水分子系统发生液固相变，即电致结冰.通过对冰相结构的进一步分析，发现电致结冰导致的冰相存在依赖于电场方向的不同冰相的层状混合堆叠结构.本文中的结果加深了人们对水分子系统物理特性的认识，并提供了关于静电场环境下分子扩散与相变等物理现象的新理解.另外，本工作中的分子动力学模拟不仅计算量小、易于实践操作，而且结果具有创新性，因此可作为教学实例培养本科生的实践创新能力和科研兴趣.

1 分子动力学模拟设计

本文采用了SPC/E水模型.其中水分子的键长和键角如图1所示.非共价键相互作用包括原子间的静电势和Lennard-Jones势,即

(1)

上式中，rij为原子i和j之间的距离，rOO为氧原子之间的距离.qi和qj为原子电荷，其中氧原子和氢原子的电荷分别为-0.847 6 e和0.423 8 e.ε0为真空中的介电常数.A=629.4×103kcal·Å12·mol-1和B=625.5 kcal·Å6·mol-1为Lennard-Jones参数.以上SPC/E水模型中的原子电荷与Lennard-Jones参数通过拟合实验测得的液相水密度、扩散常数等静态和输运性质得到.由式(1)给出的能量函数可以计算出每个原子的受力，进而通过求解运动方程得到所有原子的坐标、速度随时间的演化.分子动力学模拟系统由510个水分子组成，初始结构如图1所示.分子动力学模拟采用通用软件Gromacs在等温等压系综下进行，并采用了周期性边界条件.其中，压强设为1个标准大气压，温度设为260 K，高于SPC/E水模型对应的熔点.每次分子动力学模拟时长为5 ns，每个电场强度下独立重复32次模拟进行统计，得到平均值和标准偏差.

水分子氢键网络示意图(上)与SPC/E水分子模型几何结构示意图(下) 分子动力学模拟初始结构

2 静电场下水分子的大角度跳跃转动与扩散

2.1 大角度跳跃转动的电场强度依赖性

图2给出了水分子大角度跳跃转动示意图与代表性轨迹.根据文献[2]的研究，水分子的转动运动呈现出大角度跳跃转动，即某水分子从断开一个氢键到形成新氢键的过程中(图2)，其取向角度存在大幅度的瞬时变化.大角度跳跃转动可由水分子角平分线的取向角度θ来定量表征(其中θ是水分子1的角平分线矢量在2-1-3平面内转过的角度，如图2所示).可以看到在皮秒时间尺度内，水分子取向角度θ从-1.0 rad跳变到1.0 rad，表现出大角度跳跃转动.基于分子动力学模拟轨迹，可以统计出大角度跳跃转动发生的频次.其中，在判定一对水分子之间是否形成氢键时，根据两个水分子的氧原子间距以及氧-氢-氧原子形成的夹角数值，并采用Wernet-Nilson 判据来计算.

图2 水分子大角度转动事件示意图(上)与典型轨迹(下)

图3给出了大角度跳跃转动频次以及平均每个水分子形成的氢键数目随电场强度的变化.结果显示，大角度跳跃转动频次和氢键数目随电场强度呈非单调变化.随着电场强度的增加，大角度跳跃转动频次先增加，后减小，在0.4 V·nm-1附近时达到最大.同时，水分子平均氢键数呈现出相反的趋势，即随着电场强度先减小，后增加，且平均氢键数较小的区域转动频次较大.值得注意的是，在文献[5]中，作者基于反应平衡与玻耳兹曼分布所建立的理论模型计算得到的结论认为水分子平均氢键数随电场强度单调递增.本文基于分子动力学模拟观察到了氢键数随电场强度的非单调变化行为，展现了微观分子模拟对研究水分子结构与动力学的重要性.水分子倾向于形成四配位的氢键结构，但是对于液相水，热涨落使得实际形成的氢键数往往小于4，从而产生氢键网络缺陷.根据文献[2]的结果，这种氢键网络缺陷与大角度跳跃转动密切相关.以上结果说明，随着电场强度的增加，氢键网络结构偏离标准的正四面体配位结构,导致更多的氢键网络缺陷.随着电场强度的进一步增加，外电场引入的静电相互作用倾向于导致和电场方向匹配的新四面体序，降低氢键网络缺陷，进而减小大角度跳跃转动频次.特别是随着静电场的增强并形成新的四面体氢键网络结构，会观察到液固相变(见下文讨论).

可以由下式定义的正四面体序参量Q定量描述正四面体的形成程度[6]为

(2)

四面体序参量Q的平均值随电场强度变化图线 偶极矩对关联的径向分布

2.2 扩散系数的电场依赖性

根据以往研究结果，体相水分子的大角度跳跃转动与水分子的扩散过程密切相关[2].图5给出了不同电场强度下水分子的均方位移(MSD=〈|x(t)-x(0)|2〉，其中x(t)表示t时刻某个水分子的坐标，而x(0)表示该水分子在初始时间点的坐标；〈…〉表示对水分子以及时间的平均.)以及通过直线拟合得到的扩散系数D.其中，均方位移MSD随时间t的变化采用了三维空间正常扩散形式，即MSD=6Dt.可以看到，MSD的斜率以及对应的扩散系数D随着电场强度先增加，在0.4 V·nm-1附近达到极大值，然后随电场强度降低.值得一提的是，在Zong等人的理论与分子动力学模拟研究中[4]，观察到了外加静电场能够导致体相水分子系统黏性的各向异性以及非单调变化，这与本文得到的结果相符合.另外，比较图3与图5，我们可以看到，这种扩散系数随电场强度的变化趋势与大角度跳跃转变的变化趋势相一致，进一步支持大角度跳跃转动在水分子扩散过程中的关键作用[2].

不同电场强度下水分子的均方位移(MSD)随时间的演化 扩散系数D随电场强度的变化图线

3 静电场环境下的液固相变

3.1 电致结冰

以上关于水分子动力学的讨论主要集中在小电场区间.当进一步增强电场时，会观察到外电场导致的液固相变，即电致结冰[3].图6给出了平均氢键数以及正四面体序参量Q随电场的变化.可以看到大电场下平均氢键数目单调递增的趋势，并且在电场强度处于2～3 V·nm-1之间存在一个转变点.经过此转变点, 水分子快速逼近对应氢键饱和构型的四配位氢键结构，从而出现了有序的冰相结构.相应地，平均四面体序参量也在转变点附近发生了快速增加.以上结果说明在电场强度处于2～3 V·nm-1之间，体相水系统发生了液固相变，即电致结冰现象.这种电致结冰现象已在以往的实验和分子动力学模拟中观察到[3].

大电场下平均氢键数 平均四面体序参量Q随电场强度的变化

3.2 不同冰相的层状混合堆叠结构

通过对以上分子动力学模拟数据进行详细分析，可以研究电致结冰得到的冰相结构特征.在通常条件下，液态水结冰得到的冰相有六方相冰Ice h 和四方相冰Cubic Ice.这两种主要冰相均呈现正四面体结构，即四面体序参量Q为1，因此无法用Q值分辨这两种冰相.四面体序参量Q作为标量序，仅能表征构型接近正四面体的程度，无法表征该四面体的空间取向.而判定冰相则需要邻近四面体的空间取向排列信息，也即需要定义一个包含四面体取向信息的矢量序，并对其统计近邻关联.基于这样的思路，Moore等人[7]提出了CHILL算法来解决冰相结构的定量表征问题，将四面体取向投影到球谐函数作为单个四面体的矢量序，并由此定义近邻四面体矢量序关联系数aij：

(3)

其中

(4)

对于室温液态水、无定形冰、Ice h和Cubic Ice冰相，四面体矢量序关联aij的概率分布具有明显的区别.例如，对于Cubic Ice冰相中的水分子，其aij的分布表现为-1.0处的单峰；而Ice h冰相中水分子的aij分布同时在-1.0和-0.11处出现两个峰.图7给出的是在电场强度为0, 2, 3 V·nm-1时水分子的aij分布.在电场强度为0和2 V·nm-1时，aij在-1到+1的较大范围内均具有显著分布，表现出液态水的特性.而当电场强度增加为3 V·nm-1时，不仅可以观察到-1.0处对应Cubic Ice冰相的峰，同时观察到-0.5附近的另一个峰.这说明电致结冰所得到的冰相并非是单纯的Ice h或Cubic Ice冰相，而是存在不同峰形表征出的混合冰相.图7中灰色曲线表示3 V·nm-1下，不同的单条分子动力学轨迹得到的aij分布.结果显示不同轨迹得到两种冰相的比例具有较大的涨落，但双峰分布特征是一致的.

图8给出的是电场强度为3 V·nm-1时的电致结冰示意图.其中图8(右)展示了分子动力学模拟得到的两相混合堆叠冰相代表性结构，其中左侧区域的 Cubic Ice冰相，对应于图7中-1处的峰.右侧区域为较少报道的Polar Ice B冰相[3].可以验证，这种Polar Ice B对应图7中-0.5处的峰.

图7 不同电场强度下序参量关联aij分布

图8 电场强度为3 V·nm-1时的电致结冰示意图[在电场作用下，体相水由无序的初始构型(左)形成沿电场方向堆叠的混合冰相结构(右)]

人们通常认为具有高平移对称性的均一冰相自由能最低.最近，Lupi等学者[8]发表在Nature的工作报道了不一样的结论：由于熵效应，多种冰相混合堆叠的结构具有更低的自由能.本文中电场环境下观察到的Cubic Ice与Polar Ice B的混合堆叠结构与Lupi等人的结论相符合.值得强调的是，文献[8]中的两种冰相的堆叠是各向同性的，而本文电致结冰得到的冰相堆叠方向存在明显的各向异性.图9给出了电场强度为3 V·nm-1下分子动力学模拟得到的沿电场方向X不同位置处水分子的aij散点图，以及沿X、Y、Z方向的平均值.可以看到，沿电场方向X，aij平均值在-1.0和-0.5两种取值交替出现.而沿垂直于电场的Y、Z方向，aij平均值为处于-0.5与-0.1之间的恒定值.这说明沿着电场方向存在两种冰相的堆叠，而垂直于电场方向上不存在堆叠.以上结果展示了电致结冰过程中存在外电场导致的不同冰相堆叠方向的对称性破缺现象.

图9 冰相序参量关联〈aij〉沿X、Y、Z方向不同位置的平均值(沿电场方向X存在不同冰相的混合堆积，灰色数据点表示沿X方向aij的散点图)

4 结语

本文基于分子动力学模拟，研究了不同电场强度下，体相水分子的大角度跳跃转动、扩散以及电致结冰.相关结果首次展现了小电场区域大角度跳跃转动频次存在非单调电场依赖性.这种大角度跳跃转动的反常电场依赖性与水分子扩散系数随电场强度变化的非单调性存在较强相关性.在大电场区域，分子动力学模拟重现出前人研究中得到的电致结冰现象，并发现了电致结冰导致的冰相体系中，存在两种不同冰相沿电场方向堆叠的新现象.该分子动力学模拟涵盖水分子系统的扩散、转动运动、相变等多个物理过程，有助于学生从微观层次理解分子扩散与相变等物理现象的分子机制[9,10]，可作为教学实例培养本科生实践创新能力和科研兴趣.

致谢：作者感谢南京大学王炜教授的指导和讨论.