冰水界面动态结构的分子动力学模拟研究

王瑞，任瑛，陈卫，韩永生，3

（1 中国科学院过程工程研究所，多相 复杂系统国家重点实验室，北京 100190； 2 中国科学院大学化工学院，北京 100049；3 中国科学院粉体材料技术重点实验室，北京 100190）

引 言

冰是水的常见相态，冰的形成和融化过程对水的循环利用和其他的地球化学反应有重要意义，包括生物纳米材料中抗冻蛋白与冰的相互作用[1]、过度积冰对户外设备方面的影响以及抗结冰材料的研发[2]等。水和冰有着诸多反常的特性，如自然界中的液态水在4℃达到密度最大值：温度高于4℃时，随着温度升高，水的密度减小；温度低于4℃时，随着温度升高，水的密度增加。而在水凝固成冰时体积变大，密度变小，不同于一般物质热胀冷缩。在实验条件下，所有结晶冰的表面上的水分子都是不均匀的，表面层的结构、厚度和性质等随着温度和晶面等条件的改变发生变化，从而影响冰晶表面的生长取向，导致了冰晶的多种形态[2-3]。研究冰水相变过程中冰水界面处的动态变化规律对进一步理解冰晶成核和生长过程十分重要，也有利于进一步解决结冰造成的工程应用上的难题，但是对于冰水界面结构变化规律在分子尺度上的实验研究仍然存在一定的困难。

随着计算机技术的发展，模拟计算已经成为科学研究中继理论方法、实验方法之后的第三种研究方法，其可以用来在分子层面上解释一些实验不能揭示的现象和机制[4-5]，例如纳米界面性质与结构的关系[6]和晶体表面自组装生长过程[7]。冰晶在液态水中的生长动力学，特别是冰晶在不同晶体取向间的各向异性，是一个需要在分子尺度上详细研究的重要课题。分子动力学（molecular dynamics，MD）模拟可以在分子尺度上研究冰的成核、不同冰相之间的转化规律以及气液固三相中冰晶动态生长过程，可以实现多元条件下冰晶生长和熔化过程、冰-水以及冰-水蒸气两相平衡态结构变化等方面的模拟研究[8]。

目前已有一些针对冰晶与水的界面和冰晶表面的MD模拟研究，Karim等[9-11]使用分子动力学模拟对六角冰（ice hexagonal,Ih）和液态水的界面进行了模拟，针对水分子的简单点电荷模型[12]（simple point charge, SPC）和四位点模型[13]（the transferable intermolecular potential four-point,TIP4P），模拟了冰水界面的平衡状态，讨论了界面稳定性。Svishchev等[14]利用TIP4P 水分子模型研究了在较大的外部电场条件下过冷的液态水中冰的结晶和生长过程，其冰水界面受到外加电场的影响，进一步影响了冰晶的几何形状和生长速率。近年来，研究者对冰水界面的研究越来越深入，涉及的层次也越发广泛。在冰水界面，密度和自由能等热力学性质发生了显著的变化，如Wang 等[15]在环境压力1 bar（1 bar=0.1 MPa）下计算出不同模型水的固液界面自由能约为37 mJ·m-2（the transferable intermolecular potential four-point with Ewald techniques, TIP4P-Ew）和42 mJ·m-2（TIP5P-Ew），与之前用MD模拟计算的测量范围25～44 mJ·m-2一致。Chung 等[16]研究了极化水分子模型在不同相态之间的电荷转移情况及其对界面性质的影响，如熔化温度、介电特性、能量和密度等，认为电荷转移对水分子内和水分子间的氢键网络作用有重要意义。Louden 等[17]开展了冰水界面摩擦的非平衡分子动力学模拟，认为固液界面的氢键密度与界面摩擦力直接相关。基面{0001}、一次侧棱面{10-1 0}和二次棱柱面{11-2 0}的摩擦力相对于表面特征的剪切速率和剪切方向是不变的，三个晶面的结构和动力界面宽度都是相似的，并且很大程度上独立于剪切速率和方向。Ramírez等[18]从氢键拓扑网络结构的角度出发分析了水分子的气-液和固-液相变过程中的一些界面性质，如平均密度、聚类系数和平均路径长度等。这些性质在接近相变区域时表现出不连续的特征，这与在同一区域中某些热力学性质所表现出的不连续性是相似的，为进一步研究其他含有氢键结构的两相界面提供了更多依据。目前针对冰水界面结构的研究都是在常温常压下进行的，主要是为了在水分子模型的熔点附近得到冰水共存的平衡状态，进而分析界面热力学性质和结构性质。目前得到的普遍认识是界面静态结构的不均匀性和动态分布的不连续性都是区别于体相特征的，也是影响液态水中冰成核、生长等过程的关键因素[19]，但是界面结构在不同温度和晶面条件下的动态变化规律较为复杂，仍需要对其进行深入的研究。

本文利用分子动力学模拟方法构建冰水两相共存模型，在微正则（NVE）[20]系综下模拟冰水界面结构的动态变化。选择NVE 孤立体系可以避免控温器和控压器施加的力对动力学过程造成影响。NVE 系综下体系的总能量不变，没有外界能量和粒子的输入和输出，完全依靠体系内冰水两相之间的能量交换。研究界面层内水分子在无序和有序状态之间的转变情况，分析运动轨迹和粒子成键情况，总结界面结构的变化规律。利用序参量和密度剖面表征界面结构的动态变化过程，在不同温度和晶面条件下分析界面厚度和性质的变化规律。

1 模拟过程

本工作在超算机群Mole-8.5[21]上通过开源软件LAMMPS（7Aug19）进行算例计算，采用CUDA 版本的GPU 节点，搭载Tesla K80 进行加速计算，一个GPU有12个核，最大计算速度是普通CPU 节点的十倍[22]，可以极大地提高运算速度。将包含2880 个水分子的六角冰Ih 作为初始结构，体系盒子的尺寸大小为4.5 nm×4.7 nm×4.4 nm。选择水分子力场模型为拓展简单点电荷模型[23]（simple point charge extended, SPCE），这是一个三位点刚性全原子模型。该模型可以更好地描述液态水结构，同时三点电荷结构相比于四点电荷的TIP4P 模型能减少计算量，使用SPCE水分子模型描述六角冰Ih时，其在压强为1 bar 时的熔点为215 K[24]。初始六角冰结构的密度为0.91747 g·cm-3，模拟盒子三个方向均采用周期性边界条件，使用SHAKE算法约束水分子空间结构保持。体系为NVE 系综，给定体系不同的初始温度，观察体系结构在孤立系统中的演化。系统的平衡温度和压强可以根据SPCE 水分子模型的相图中水的液相和固相共存曲线确定。

首先对系统进行5000步的能量最小化处理，优化体系结构为更合理的几何结构。用平面z=0.18 nm将初始冰晶结构分成A 和B 两部分，垂直z轴的xy平面为六方晶系的基面{0001}，因此冰水接触的暴露晶面即为基面{0001}。在初始化设置中对A部分的冰晶结构进行固定，即其位置和速度不随时间更新，而B部分的水分子根据麦克斯韦-玻尔兹曼分布的随机数列所产生。弛豫20 ps后解除对A部分位置和速度的固定设置，使体系内两部分进行能量交换，最终体系会重新达到新的平衡态。此时再对整个体系进行温度和压强的弛豫，在等温等压（NPT）系综下弛豫200 ps，采用Nosé-Hoover 恒温器和恒压器来控制体系的温度保持恒定，随着模拟盒子的体积减小，在平衡态时液态水的密度为0.98759 g·cm-3，接近实验值1.0 g·cm-3。最后再进行NVE系综下的10 ns预弛豫过程，体系达到固液两相平衡，体系温度、压强等热力学参量都达到平衡。NVE 孤立体系中存在界面，可以使体系的温度保持在熔点附近达到动态平衡，同时避免形成共存两相平衡态造成的体系过冷问题。对B 部分设置不同的初始温度，经过能量平衡弛豫后，最终体系为冰水共存状态的温度，即为相变温度。若初始温度高于相变温度，B 部分的水分子会将过剩的能量传递给A 部分的冰，导致冰融化；反之，若初始温度低于相变温度，经过能量交换，B 部分水分子运动速度减慢，变成有序冰晶结构，体系逐渐变成全冰相。将轨迹文件导入可视化分子动力学软件VMD (visual molecular dynamics)[25]中可以观察在不同平衡温度时体系结构的变化。在NVE 系综下利用热力学参量来验证体系是否达到平衡态是非常必要的，在温度和压强达到平衡后，统计在NVE 系综下在弛豫的30 ns 内两相体系的各热力学参量。体系的总能量不变，在30 ns的平衡弛豫过程中其波动约为0.028%，体系的温度在熔点215 K附近波动，波动范围为±2 K。

2 结果与讨论

2.1 利用序参量表征体系中冰水界面结构动态变化

冰水两相共存的界面是冰晶成核和生长的区域，水分子在冰晶结构表面的动态复杂行为影响冰晶成核和生长过程。因此探究冰晶生长前端的界面微环境首先要对界面结构有清晰的认识。本文采用序参量描述冰水界面结构的变化[26-28]，冰的晶格结构是正四面体的对称结构[29]，但是在液相中由于水分子的无序移动，四面体结构会发生一定程度的扭曲。四面体序参量qi可以用来描述四面体的这种对称性[30]，qi定义为

式中，ψjk为一个给定分子的氧原子与其最近邻的j和k氧原子之间的连线夹角,(°)。如果一个分子位于正四面体的中心，其顶点被最接近的四个水分子中的氧原子所占据，则cosψjk=-1/3。

当qi值接近1 时，表明四面体的对称性较高，此时水分子组成高度有序的冰晶结构；而当qi的值在0附近时，表明四面体的对称性较低，此时水分子组成无序液态水。因此计算四面体序参量可以在一定程度上区别冰水混合结构中类冰水分子和类液水分子，选取合适的序参量的值可以更好地描述冰水界面处水分子的结构转变以及界面结构波动变化规律。

首先分别对冰和水体系进行序参量计算，对比液态水和固相冰晶的序参量q6概率分布曲线。结果发现，在截断半径rc=0.3 nm 时，两者的区分程度最明显，对应的q6阈值为0.41。当序参量q6取值大于0.41时，该对应的水分子即为类冰水分子；当序参量q6取值小于0.41 时，该对应的水分子即为类液水分子，从而可将体系中的类冰水分子和类液水分子区分开来。

图1 为本研究的基本思路，采用序参量来表征液态水中冰晶生长过程中界面结构动态变化规律，如图1(a)所示；图1(b)从左到右依次为液态水、冰水共存和固态冰晶三种状态的结构图及其对应的序参量概率分布曲线。

图1 采用序参量表征冰晶生长过程中冰水界面微观结构的动态变化过程(a)；液态水、冰水共存和固态冰晶体系分子结构及其水分子序参量概率分布曲线(b)Fig.1 The complex dynamic process of the interface microstructure in the process of ice crystal growth in liquid water characterized by order parameters(a);The structure of liquid water,ice-water and ice crystal with probability distribution curve(PDC)of order parameters for water molecules(b)

进一步，利用氧-氧径向分布函数验证不同状态时的体系结构。径向分布函数可用来表征体系内固体或液体分子排列的有序程度，定性判断体系结构。其表达式为：

式中，r为任意粒子与球心(目标粒子)的距离，nm；ρ为体系的密度，g·cm-3；dN是厚度为dr的球壳内的颗粒数量。径向分布函数可以看作整个体系内局域密度与平均密度的比值[31]。

在体系达到平衡态后，统计体系结构的氧-氧径向分布函数的时间平均值，结果如图2 所示。体系结构随温度变化明显，体系温度低于熔点215 K时，体系的特征峰明显，与文献[32]中的冰晶结构一致，此结构长程有序，体系为纯冰状态；体系温度高于熔点时，第一峰值和第二峰值的强度逐渐降低，且远程峰的晶格结构遭到破坏，分子间作用力减弱，特征峰消失，曲线波动趋势减小，逐渐向1靠近，即曲线归一化明显。体系温度在熔点215 K 附近时，氧-氧径向分布函数介于两者之间，可以判断为冰水两相共存状态。

图2 不同温度下体系的氧-氧径向分布函数对比Fig.2 Comparison of the radial distribution function gO—O(r)of the system at different temperatures

2.2 冰晶界面结构动态演变

利用LAMMPS 构建冰水两相共存模型，冰水接触的界面为{0001}。经过10 ns的预弛豫之后，在熔点附近再进行30 ns的平衡弛豫过程，得到体系中水分子的运动轨迹，在VMD 中观察其三维动态变化过程。

冰水两相达到平衡状态后，选取界面处某一固定位置观察其结构随时间的变化，可以得到不同时刻的界面结构。如图3(a)、(b)所示，以t=1.5、2.3、5.6、7.9 ns四个时刻为例，t=1.5 ns 时，界面该区域全部为冰晶的六元环结构；t=2.3 ns 时，该区域出现了部分无序结构；t=5.6 ns时，该区域全部变为无序的分散结构；t=7.9 ns时，该区域又重新恢复部分有序的冰晶六元环结构。

模拟结果表明，在微观原子尺度上，冰水两相平衡后的界面仍在发生不同程度的移动，界面处无序水分子被吸附在冰面上，与冰层表面的水分子结合生成氢键进而排列成有序的冰晶六元环结构，而同时界面层内冰晶结构的水分子也发生分子间氢键的断裂，脱离晶体表面进入界面层，扩散到液态水分子的体相之中，这两种变化时刻发生。界面层内同时包含无序的液态水分子结构和有序的冰晶结构，两者之间存在可逆动态转变。进一步，通过序参量q6定量统计体系中类冰水分子数目，得到其随时间的变化规律，如图3(c)所示，发现类冰水分子数目在动态平衡过程中也呈现波动变化，证明冰水界面处发生相变过程，氢键的形成和断裂同时进行，并存在局部的能量交换。

图3 不同时间冰水界面处水分子结构(a)；冰水界面示意图(b)；两相平衡态时界面处类冰水分子数目随时间的变化(c)Fig.3 The structures of the ice-water interface at different time(a);Schematic diagram of the ice-water interface(b);The number of ice-like water molecules changes with times at the two-phase equilibrium condition (c)

2.3 温度对冰水界面结构的影响

利用序参量表征冰水界面结构随着温度的变化规律，初始状态为冰水两相共存的平衡态，选择冰水接触面为二次棱柱面{11-2 0}，见图4。

温度为215 K 时，体系中水分子的序参量概率分布曲线为明显的双峰分布，表明体系中冰水两相共存。在降温过程中，右侧分布逐渐增加，说明类冰水分子在不断增加；左侧分布逐渐减小，最后变为靠右侧的单峰，说明类液水分子逐渐转化为类冰水分子，如图4(a)所示。在升温过程中，序参量的分布逐渐变窄且分布整体左移，同时左侧分布增加，右侧分布减少，最后变为靠左侧的单峰，说明此时类冰水分子在逐渐转化为类液水分子，如图4(b)所示。

图4 不同温度下体系水分子序参量概率分布曲线：（a）降温过程；（b）升温过程Fig.4 Probability distribution curve(PDC)of order parameters at different temperatures:(a)cooling process;(b)heating process

进一步，统计沿垂直于界面方向水分子的密度变化，得到如图5 所示不同温度时体系的密度剖面对比图。在过冷度较大时，界面密度波动较大，说明界面结构趋于有序状态。因此界面内水分子密度较大，水分子在界面层内的停留时间较长，形成六角冰结构的概率增大。同时冰相的峰强度更高，其结构高度有序，界面层厚度较薄，冰水界面区分明显。随着温度升高，界面层厚度增加，界面内水分子密度减小，无序性增加。可以预测，随着温度逐渐升高到熔点以上，界面处的峰强度逐渐降低，最终结构密度在1.0 g·cm-3（液态水的密度）处波动。

图5 在不同温度下二次棱柱面{11-2 0}的界面密度剖面图Fig.5 Density profile of the ice-water interface formed by secondary prismatic plane{11-2 0}

2.4 晶面对冰水界面结构的影响

六角冰的生长晶面主要有三个，基面{0001}、一次棱柱面{10-1 0}和二次棱柱面{11-2 0}。其中一次棱柱面{10-1 0}和二次棱柱面{11-2 0}的暴露界面具有几何粗糙结构，而基面{0001}具有分子平面结构。基面和棱柱面的双层高度分别为0.37 nm 和0.39 nm。在冰水共存体系中，界面处存在大量吸附的水分子，其与周围的水分子没有形成理想的六方晶系结构，但其分子间氢键作用力远大于体相中液态水分子之间的相互作用力，界面冰晶水分子的排列方式不同，即暴露的氢键数目和角度不同，会对界面处液态水分子产生不同取向和大小的作用力。

针对以上三个晶面分别建立冰水两相共存体系，利用序参量q6定量表征不同晶面形成的冰水界面结构变化规律。三个体系中以基面为接触面形成初始两相结构时的类冰水分子和类液水分子数目分别为1650和1230，以一次棱柱面为接触面形成初始两相结构时的类冰水分子和类液水分子数目分别为1480和1400，以二次棱柱面为界面的初始两相结构时的类冰水分子和类液水分子数目分别为1536 和1344。在经过弛豫达到两相共存的平衡态时，分别统计三个体系中类冰和类液水分子的序参量，得到图6 的概率分布曲线。其中体系的冰水接触面分别是(a)基面、(b)一次棱柱面、(c)二次棱柱面。经过时间平均和系综平均后得到两相平衡状态时类冰水分子和类液水分子数目。分别对比初始两相结构，发现一次棱柱面和二次棱柱面形成的两相体系中类冰水分子和类液水分子数目更接近，在孤立体系的模拟环境中也更稳定。

图6 不同晶面形成的冰水两相体系的序参量概率分布曲线，数字分别代表类冰水分子（右）和类液水分子（左）数目Fig.6 Probability distribution curve(PDC)of the order parameters in ice-water system formed by three different crystal planes,and the figures in the histogram represent the number of water molecules at liquid (left)and solid phases(right),respectively

进一步，沿着垂直于界面方向计算密度剖面观察界面结构，如图7所示。其中图7(a)为基面形成的冰水界面结构。其表面趋向于平坦，这与实验[33]观察结果一致。图7(b)为一次棱柱面{10-1 0}对应的冰水界面结构，与基面相同，一次棱柱面的每一层都有一个双峰，第一个双峰含有六角形环上的三个氧，第二个双峰与其他三个氧对应。图7(c)为二次棱柱面{11-2 0}对应的界面结构，其密度剖面为单峰，这是由界面水分子的排列结构决定的。

图7 不同晶面形成的冰水界面的密度剖面图Fig.7 Density profile of the ice-water interface formed at different crystal planes

通过计算三个不同晶面从冰晶结构过渡到液态水分子的界面层厚度，发现在相同的温度下，一次棱柱面和二次棱柱面的界面层厚度分别为0.35 nm 和0.20 nm，小于基面的界面层厚度0.40 nm。基面的边缘结构平坦，界面层较厚，氢键在晶面上发生二维重组；棱柱面形成的界面层比基面窄，界面水分子的有序性强，氢键在晶面上发生三维重组，边缘结构接近六元环。

本研究在分子尺度上阐明了体系间生长过程中界面结构的各向异性，即不同晶面的生长速度不同会导致冰晶生长形态的多样性。因此不同于气相雪花形成过程中基面生长速度比棱柱面快，液态水中冰晶在过冷度相同的条件下其棱柱面生长速度更快，倾向生成侧枝较长的枝晶。

3 结 论

本文利用分子动力学方法模拟研究了冰水两相共存体系中冰水界面结构的动态变化过程，揭示了不同温度和晶面条件下界面结构的变化规律，通过序参量、径向分布函数和密度剖面等对该过程进行了定量表征，得到主要结论如下。

（1）在NVE系综下，当冰水两相达到平衡态时，冰水界面结构在有序冰晶六元环和无序水分子之间可逆动态转变。

（2）晶面类型影响冰水界面动态结构，基面结构较为平坦，形成的冰水界面层较厚，分子间氢键在晶面上发生二维重组，水分子排列的有序性较差；而一次和二次棱柱面结构粗糙，分子间氢键在晶面上发生三维重组，有序性更强。

（3）过冷度对冰水界面动态结构及其演化趋势有重要影响，过冷度越大，界面结构越趋于冰晶六元环的有序状态，界面内水分子密度越大，水分子在界面层内停留时间越长，晶面生长的概率越大。