English

• 1.   引言

  限域水是近年来水科学领域重要的研究热点.在极窄的限域空间内, 水分子间结构与氢键网络会因为各种热力学条件, 例如:限域压强、温度、限域墙提供的波纹状势场等, 组合成极其丰富的物相结构类型.例如五次、六次对称性^[1~7]晶相结构及近期实验观测到的限域方形晶相结构^[8].这些特殊的原子排列结构使得限域水有着异于块体水的特殊性质, 让其在诸多科技领域如生物医药、纳米摩擦、海水淡化、微流控等科技领域中^{[9, 10]}都具有重要的应用潜力.

  2015年, Algara-Siller等^[8]运用透射电镜直接观测到了范德瓦尔斯作用下两片石墨烯中二维限域水的晶体结构—完美的正方形晶格结构, 方冰.这一工作开创了实验直接观测、表征二维限域水结构的先河.除此工作之外, 更多的关于限域水的相结构类型和相图的知识, 来自于第一性原理计算和分子模拟研究.

  在已有的限域水的研究中, 更多的关注集中在限域水的相结构本身^{[11, 12]}, 而不同结构相之间的转换涉及较少. Algara-Siller等的实验工作中报道了电子束的移动会导致方块冰晶界的运动与无序化, 这个有趣的实验发现意味着应用外场控制限域水不同结构相之间的相变过程的可能性.曾经有少数研究讨论过不同限域水结构相之间的转变^[13~18], 但这些研究中没有从一个真正的热力学两相平衡态出发, 也很难获得任何关于限域相变相关的关键界面动力学与热力学性质的定量计算结果.

  本论文中, 我们提出一套基于平衡态分子动力学的分子模拟技术, 可以实现限域单分子层冰和水的恒温、恒压、恒体积下精准的两相平衡, 并用以系统研究限域冰-水相平衡的结构与性质.应用此模拟技术, 我们模拟了SPC/E水模型描述的单分子层厚度的水在0.65 nm限域尺寸、5000 bar限域压强条件下的冰-水(固-液)两相平衡, 测得该限域条件下的相平衡温度.在长时间相平衡状态下(10 ns), 系统表征了该平衡界面(线)的一系列热力学量的面内分布函数.本论文通过平衡态的分子模拟, 展示了粗糙型固-液分界线的热毛细涨落、固-液结构转变的微观机制、以及缺陷在固-液相变区附近的形成与输运, 并且由各种热力学量分布函数呈现了二维限域冰-水共存界面(线)的特殊性质.本论文提出的模拟技术对于研究两相共存的动力学(二维相变)和热力学性质(二维相规律), 及检验相平衡热力学理论在低维体系的适用情况^[19]等具有重要意义.

  我们将在下文中的计算方法部分介绍水模型、限域条件的选取, 阐述分子动力学模拟的具体细节、详细介绍固-液两相平衡建立的方法, 和介绍分析方法.在结果与讨论部分, 我们将报道单分子层固-液平衡态模拟结果, 展示界线处水分子排列结构特性, 并统计分析了界线几个关键热力学量的分布函数.文章的最后, 我们总结本论文的主要结果, 探讨论文所提出的低维限域冰-水相平衡模拟技术的应用前景.

  2.   计算方法

  2.1   水模型

  在Algara-Siller等^[8]的研究工作中, 通过分子模拟对比了SPC/E和TIP4P两种水模型, 得到了近似的模拟计算结果.我们的目的是发展一种分子动力学模拟技术, 水模型的选取并不是核心问题.因此我们在本研究中选用了一个简单的刚性水分子模型, 即SPC/E力场模型^[20].原子所带电荷量和质量都集中在同一个质点上, 氧原子带负电、电荷量大小为0.8476倍电子的电量, 氢原子带正电、电荷量大小为0.4238倍电子的电量.水分子间通过氧原子的Lennard-Jones(LJ)势作用:

  $ {\mathit{\Psi } _{{\rm{oo}}}}(r) = 4{\varepsilon _{{\rm{oo}}}}\left[ {{{\left( {\frac{{{\sigma _{{\rm{oo}}}}}}{r}} \right)}^{12}} - {{\left( {\frac{{{\sigma _{{\rm{oo}}}}}}{r}} \right)}^6}} \right] $

  (1)

  其中ε_oo＝0.15535 kcal/mol, σ_oo＝0.3166 nm, 截断半径2.5σ_oo.

  2.2   受限条件与限域墙

  本工作模拟单分子层固、液相水的受限条件为限域尺寸0.65 nm, 限域压强5000 bar.此受限条件的选择参考了文献^[8]的研究成果, 0.65 nm限域尺寸是能获得单分子层限域冰的最小限域尺寸, 5000 bar的限域压强拟模拟实验中单分子层冰受限于双层石墨烯间范德华力致使的环境压强.

  SPC/E模型的所模拟的限域单分子层方冰结构, 在文献^[8]中被证实不依赖于限域墙模型是简单的刚性墙还是柔性石墨烯.对于刚性墙模型, 受限冰的结构同样独立于墙-水相互作用的强度或疏水性, 以及刚性墙有无原子结构.在本工作中, 水分子被置于两个平行无原子结构的刚性墙内^[21~24].墙与水分子的相互作用采用LJ10-4-3势能函数^[25]:

  $ \begin{array}{l} {\mathit{\Psi } _{{\rm{wo}}}}(z) = 2\pi {\varepsilon _{{\rm{wo}}}}\left( {\frac{4}{5}} \right)\left[ {\frac{2}{5}{{\left( {\frac{{{\sigma _{{\rm{wo}}}}}}{z}} \right)}^{10}} - {{\left( {\frac{{{\sigma _{{\rm{wo}}}}}}{z}} \right)}^4} - } \right.\\ \quad \quad \left. {\frac{{\sqrt 2 }}{{3{{\left( {\frac{z}{{{\sigma _{{\rm{wo}}}}}} + 0.61\sqrt 2 \frac{z}{{{\sigma _{{\rm{wo}}}}}}} \right)}^3}}}} \right] \end{array} $

  (2)

  ε_wo＝0.261856 kcal/mol, σ_wo＝0.326747 nm, 截断半径2.5σ_wo.上式通过积分单层刚性六元碳环中碳原子与近邻水分子中氧原子间的碳-氧相互作用^[26]得到. LJ10-4-3势能函数相比传统的Lennard-Jones势函数, 在描述限域墙-受限水分子间相互作用上更加准确.相比柔性石墨烯墙的模拟, 该方法更能抓住物理本质，且减少计算.

  2.3   分子动力学模拟

  我们构建直角坐标系框架下的模拟盒子.如图 1所示, 模拟盒子在x和y方向上应用周期性边界条件, 单层水分子分布在x、y平面.水分子层上下各有一个理想刚性墙对水提供限域, 因此在z方向上没有应用周期性边界条件. z方向上0.65 nm的限域尺寸提供了单分子层水/冰的存在的空间基础.本论文所采用的限域墙与Algara-Siller等工作中应用刚性单层石墨烯作为限域墙, 有微小区别.通过检测, 我们发现两种刚性限域墙的选取对主要结果没有明显影响.模拟盒子在x、y两个方向上的尺寸分别为42 nm、4.5 nm, z方向尺寸等于限域尺寸0.65 nm.盒子中大概有2000个水分子, 固相、液相各占一半.

  图 1

  

  图 1.  模拟盒子的示意图

  Figure 1.  Schematic diagram of the simulation box

  Single layer of water molecules are confined in a slit pore consists of two parallel rigid walls

  下载: 全尺寸图片 幻灯片

  在当前模拟盒子与周期性边界条件组合下, 为了准确计算长程库伦相互作用, 同时为减少计算耗时, 我们应用Particle-Particle Particle-Mesh算法^[27]处理长程静电相互作用(收敛误差精度取10^－5).由于模拟中z方向为非周期性, 我们遵循Yeh等^[28]去除模拟盒子与镜像之间偶极-偶极相互作用的方法, 引入一个z方向长度为50 Å的真空区域.

  分子动力学模拟的实现来自于美国Sandia国家实验室开发的LAMMPS软件^[29].模拟过程中, 每个“刚性”水分子的运动方程的时间积分应用SHAKE算法^[30], 使每个氢氧键的键长和同一分子两个氢氧键键角分别维持在1.0 Å和109.47°.时间步长取为1 fs (1 fs＝10^－15 s), 模拟温度介于300 K到450 K之间, 控温利用Nosé-Hoover恒温法^[31], 弛豫时间为100.0 fs, 控压利用Anderson恒压法^[31], 弛豫时间为1000.0 fs.我们调控限域水x和y方向上压强维持在5000 bar, 以模拟实验上石墨烯间范德瓦尔斯力致使的限域压强^[8].

  2.4   限域固-液相平衡模拟技术

  基于分子动力学方法的固-液两相平衡模拟技术已成为精确计算固-液界面诸多热力学性质的关键性支撑技术^[32~34], 例如:固-液界面自由能^[35]、过剩应力^{[36, 37]}、溶质偏析^[38]等界面热力学量的计算. Benet、MacDowell和Sanz等利用固-液两相平衡模拟技术研究了, 水-冰固-液界面的毛细波动力学与界面自由能^{[39, 40]}.为了建立高限域压强下的准二维限域冰-水平衡, 我们需要对现有固-液两相平衡模拟技术进行拓展.

  图 2为限域单分子层冰-水两相平衡态模拟技术的流程示意图.第一步, 在两相平衡温度(熔点)附近的试探温度T_i下, 分别制作横截面尺寸A＝L_y(T_i)L_z完全相同的限域单分子层方冰晶体与水液体样本(见后).第二步, 把截面尺寸相同的固、液样品组合同时对组合后的两相共存的样品应用周期性边界条件.第三步, 进行横截面积和粒子数不变的等压-等焓(NP_xAH)分子动力学模拟––基于NPH系综^[41]的分子模拟方法, 这种方法以其自发调整温度的优势较多应用于固液相变的分子动力学模拟中^[42]. P_x＝5000 bar, 持续5 ns平衡固、液两相.在这一步骤中, 体系温度(即粒子平均动能)会自行调节.如果系统温度不在熔点, 固液两相的自由能差别将驱动相变, 但由于焓受到约束, 固-液相物质仅能有限转换, 其引起的体系势能数值的变化会带动动能朝反向变化, 于是体系温度将自发调节, 达到更接近熔点的温度T_j. T_j和T_i的不同随即导致晶体部分晶格尺寸与预设定的A出现了不匹配的情况, 需要依据新温度调节体系截面积A.第四步, 测量第三步中平衡温度T_j, 并依据单分子层方冰热膨胀函数关系计算T_j对应的L_y(T_j), 缩放截面至A＝L_y(T_j)L_z, 再次进行NP_xAH分子动力学模拟, 持续5 ns.第五步, 迭代多次上一步内容, 测量新的平衡温度、调整截面、进行NP_xAH模拟.直到体系中固液两相趋于平衡、体系温度T_j收敛, 至两相平衡温度(熔点)T_m.第六步, 在最终平衡的NP_xAH模拟中计算平均L_x长度, 固定该长度同时固定截面A, 进行10 ns长度的NVT_m模拟, 每ps (1 ps＝10^－12 s)记录一次所有粒子坐标和应力信息用以后续统计分析与计算.

  图 2

  

  图 2.  限域单分子层冰-水两相平衡态模拟技术流程图

  Figure 2.  Flow chart of the MD simulation technique for obtaining the equilibrium coexistence of the monolayer confined ice and water

  下载: 全尺寸图片 幻灯片

  需要说明两点. (1)截面积A在本方法中是至关重要的参量, 寻找熔点的过程中如果忽略晶体晶格对温度微小调整后的变化, 将导致所建立的固-液平衡体系严重偏离的热力学平衡条件.可由公式(7)或者图 3中晶格热膨胀结果, 精确算出体系温度变化后A的需要调整数值. (2)在制作两相平衡初始坐标之前, 应用NP_xAT_i模拟, 分别制作横截面尺寸A＝L_y(T_i)L_z, 限域压强为5000 bar限域单分子层方冰晶体以及限域单分子层水液相样品, 模拟时间长度5 ns.每500 ps输出一个构型结构坐标, 共产生10个固相样品的初始结构和10个液相样品的初始结构.两相平衡态模拟以这10组独立初始坐标出发进行10次重复性模拟.结果显示10次平衡模拟获得的相平衡温度在误差范围是一致的, 不同初始结构对最后固液界面平衡结构的影响非常小. 10次重复性模拟可以提供充分的统计样品和优秀的统计精度.

  2.5   限域冰-水相平衡的热力学参量计算

  为了表征限域单分子层冰-水相平衡体系中, 原子排列结构性质与热力学性质通过相平衡界线的变化行为, 我们计算了一系列热力学参量沿x轴(界线法向)的分布函数^[43].这些分布函数不但能够定量表征界线重要结构与热力学性质的变化规律, 而且可以作为检验2.4节中我们对于两相平衡模拟技术拓展是否成功的重要标准.

  (1) 细粒度密度分布:对xy平面沿x方向上进行细分, 获得宽度为Δx的细粒度“切片”, 本研究中Δx＝0.005 nm, 统计每个“切片”内粒子数面密度, 即单位内氧原子或氢原子数目, 计算公式如下:

  $ \rho (x) = \frac{{\left\langle {{N_x}} \right\rangle }}{{\Delta x{L_y}}} $

  (3)

  其中 ${\left\langle {{N_x}} \right\rangle }$ 为“切片”氧原子或氢原子数目, ΔxL_y细粒度“切片”面积.

  (2) 粗粒度密度分布:利用有限冲激响应(Finite Impulse Response)光滑算法^{[44, 45]}, 把细粒度密度分布转换成粗粒度密度分布.粗粒度分布结果能够更加直观地展现固体与固-液界面区域的热力学量的变化趋势, 辅佐判断平衡态分子模拟的统计抽样是否充分.在这个算法中, 第n个细粒度“切片”中物理量f_n的平均值 $\overline {{f_n}} $ 的计算公式为:

  $ \overline {{f_n}} = \sum\limits_{k = - N}^N {C{e^{( - {k^2}/{\varepsilon ^2})}}} {f_{n + k}} $

  (4)

  C是归一化参数, N和ε是决定加权平均权重因子的参数, 我们取N＝2200, ε＝760.液相界面如液-汽界面, 细粒度/粗粒度密度分布没有区别.固-液界面两种分布有明显的区别.界面性质的细粒度分布对于获得精确的固-液界面动力学理论十分重要, 粗粒度分布则可提供固-液界面连续介质理论所需要的重要数据.

  (3) 压强分量分布:类似于上述密度分布, 细粒度压强张量分布的计算基于统计每个Δx“切片”中所有粒子的理想气体压强分量与维里压强分量之和:

  $ {P_{mn}}(x) = \frac{{\left\langle {\sum\limits_i^{{N_x}} {{M_i}{v_{im}}{v_{in}}} } \right\rangle }}{{3\Delta x{L_y}{L_z}}} + \frac{{\left\langle {\sum\limits_i^{{N_{x'}}} {{r_{im}}{f_{in}}} } \right\rangle }}{{\Delta x{L_y}{L_z}}} $

  (5)

  其中m, n＝x, y, z, 求和号是对每个Δx“切片”中所有原子求和, M_i是第i个原子的质量, ν、r和f是第i个原子的速度、位置坐标和受到的作用力矢量, ΔxL_yL_z是“切片”体积.同样地, 我们利用有限冲激响应光滑算法, 把细粒度分布转换成粗粒度分布.

  (4) 平面内应力分布:从粗粒度P_xx与P_yy压强分量分布, 我们计算平面内应力S(x)的分布:

  $ S(x) = {P_{xx}}(x) - {P_{yy}}(x) $

  (6)

  3.   结果与讨论

  3.1   热膨胀密度回线

  为获得精确熔点, 我们首先测量了侧向压强为5000 bar, 限域冰-水的平衡态密度与温度的变化关系.这个测量基于NP_xyT系综的单相平衡态分子动力学模拟, 模拟盒子在x、y、z三个方向上的尺寸分别为4.5、4.5和4.5 nm, 包含256个水分子.我们遵循逐渐升温和逐渐降温的方式, 模拟温度从350 K到470 K中每隔10 K温度下的平衡态限域水或冰.每个温度的模拟时长为10 ns, 每1 ps记录一次所有粒子坐标用以统计分析.同样地, 每个体系均以不同的初始坐标和初始速度重复模拟10次以获得优秀的统计精度.

  如图 3画出了升温与降温模拟中, 所有温度下冰或水的氧原子面密度结果.在升温过程中, 限域冰的氧原子密度从0.131 Å^－2减少至0.128 Å^－2, 然后发生固-液相变, 密度跳变至0.118 Å^－2, 持续升温至470 K, 密度持续减少至0.117 Å^－2.从470 K开始降温, 体系能够维持单分子层液相至390 K, 密度单调增加至0.121 Å^－2.温度继续降低时, 体系再次经历相变(液-固相变).从380 K降至350 K过程中测量得到的固相密度与之前升温模拟时完全一致.以三次多项式拟合图 3中固相与液相密度的计算结果, 可以获得上述温区内固相与液相单分子层限域水氧原子面密度的数学表达式, 固相密度(单位为Å^－2)为:

  $ \begin{array}{l} \rho _{\rm{O}}^{\rm{s}}\left( T \right) = 0.1456 - 8.2338 \times {10^{ - 5}}T + \\ 2.0430 \times {10^{ - 7}}{T^2} - 2.53850 \times {10^{ - 10}}{T^3} \end{array} $

  (7)

  图 3

  

  图 3.  单分子层限域水和冰的氧原子面密度随温度变化.侧向压强5000 bar, 限域尺寸0.65 nm.实线和虚线分别为升温与降温路径, 方块与圆圈分别代表限域冰与水的密度数据. “热滞回线”型的密度-温度变化规律意味着当前单分子层限域冰-水体系经历了一级(固-液)相变

  Figure 3.  Oxygen area density as the function of temperature for the confined monolayer ice and water. Lateral pressure 5000 bar, pore size 0.65 nm. Solid line and dashed line refer to the heating and cooling path, respectively. The density hysteresis implies that the confined monolayer ice and water undergoes a first order (solid-liquid like) phase transition

  下载: 全尺寸图片 幻灯片

  液相密度为:

  $ \begin{array}{l} \rho _{\rm{O}}^{\rm{l}}\left( T \right) = 0.2343 - 6.6494 \times {10^{ - 4}}T + \\ 1.3488 \times {10^{ - 6}}{T^2} - 9.8883 \times {10^{ - 10}}{T^3} \end{array} $

  (8)

  “热滞回线”型的密度-温度变化规律意味着当前单分子层限域冰-水体系经历了一级(固-液)相变.固相和液相能够在370到410 K的温度区间共存, 意味着固-液两相平衡温度、即熔点就在这个温度区间内.体系相变需要历经熔化或固化形核过程, 导致了过冷液相和过热固相这两种亚稳相在熔点附近以平衡状态出现.需要说明的是, 对于单分子层这一准二维体系, 图 3中展现的一阶相变特征, 否定了目前体系存在KT相变和hexatic相的可能.这其中的原因是, 目前0.65 nm的限域空间内, 水分子仍有较大的自由度可以进行旋转和在z方向上做微小偏移, 仍有区别于理想的二维物质.

  3.2   单分子层限域冰-水相平衡原子结构

  以图 3中回线区域居中的温度, 410 K为试探温度, 我们根据图 2的流程建立起单分子层限域冰-水相平衡态, 最终确定了在5000 bar压强、0.65 nm限域条件下固-液相平衡温度(或熔点)为402.5 K.在这个温度下, 没有连续固化或连续熔化发生, 固相区域与液相区域均没有整体平动发生, 没有任何固相在液相中飘移的现象出现.

  图 4即为NVT_m模拟中一个时刻下模拟体系所有分子的坐标快照.体系中部为固相的单分子层限域冰, 其中氧原子的排列呈现长程有序的平面晶格结构, 大部分氢氧键的方向集中在xy平面内.这些限域冰晶体结构的特点与Algara-Siller等工作中模拟结果非常一致.体系的两侧为液相的单分子层限域水, 通过周期性边界条件在x方向上连接并与中部固相形成左右两个相平衡界线.液相限域水中, 平动与转动型扩散大大高于固相, 总体呈长程无序性结构.与此同时, 我们观察到水分子会形成纳米尺度的局域有序的团簇结构, 在皮秒量级的短时间尺度下, 形成团簇的十到二十个水分子具有明显的集体性, 参与液相限域水的密度涨落与物质输运.为了更加清晰地区分固相和液相并提取出固-液相界线, 我们采用Morris提出的结构序参量^[46], 来定量描述每一个水分子居于结构偏离理想固相晶格的程度(这里我们用二维正方晶格作为理想晶格).理想晶格序参量取1, 完全无序的液相分子序参量趋于0. 图 5中(a)与(b)是同一时刻体系右侧固-液界线区域分子坐标快照, 图 5(b)中展示了以色谱颜色变化标注的水分子(氧原子)的结构序参量数值, 蓝色代表高序参量数值的固相, 红色代表低序参量数值的液相.结构序参量中色彩变化可以更加直观表征出固-液分界线的位置, 从而辅助追踪限域单分子层固-液两相共存的平衡态界线涨落.例如, 图 5(c)和(d)展示出固-液两相在界线区域可以出现大幅的界线涨落和超过10倍晶格尺寸固-液两相交换区间.

  图 4

  

  图 4.  单分子层限域冰-水相平衡分子动力学模拟分子坐标快照.压强5000 bar、温度402.5 K

  Figure 4.  Snapshot of the single layer confined ice-water equilibrium interface obtained from a molecular dynamic simulation. Lateral pressure 5000 bar, under a coexistence temperature of 402.5 K

  下载: 全尺寸图片 幻灯片

  图 5

  

  图 5.  (a) 固-液界线区域分子坐标快照, 红球为氧原子、白球为氢原子; (b)图(a)中氧原子的坐标与结构序, 原子的颜色变化遵循色谱变化, 代表了结构序参量数值改变, 蓝色代表固相, 红色代表液相; (c), (d)固-液界线区域的涨落; (e)~(h)、四个相隔1 ps的连续快照, 黑色的线圈标出了一个点缺陷的扩散

  Figure 5.  (a) Snapshot of the region near the solid-liquid coexistence line; (b) Oxygen atoms in panel (a) is colored coded based on the value of the order parameter, blue for solid and red for liquid; (c), (d) Fluctuation of the coexistence line due to the in-plane solid-liquid phase transition; (e)~(h) Four successive trajectories with a 1 ps time interval, the diffusion of a point defect is marked by black circle

  下载: 全尺寸图片 幻灯片

  从固-液分界线位置出发, 我们进一步统计分析了界线附近的固相冰和液相水的结构性质, 如图 6(a)和图 6(b)所示, 二维的径向分布函数以及结构序参量数值的统计概率密度, 呈现出处于相平衡条件下单分子层限域冰和水的结构差异, 同时反映出相平衡界面区域两相此消彼长、结构共存的结构特性.

  图 6

  

  图 6.  (a) 限域单分子层冰-水两相平衡界线附近固、液相二维氧-氧径向分布函数. (b)固、液相限域水中氧原子结构序参量的统计分布概率

  Figure 6.  (a) A comparison of the 2d O—O radial distribution function between the confined monolayer ice and water phases near the equilibrium crystal-melt interfaces. (b) Probability distributions of the structural order parameter for the oxygen atoms in the confined monolayer ice and water phases

  下载: 全尺寸图片 幻灯片

  图 5(e)到5(h)记录了四个相隔1 ps的连续快照.我们观察发现界线相变, 主要以水分子团簇的形式而非单分子, 在界线处参与固-液转变.例如从5(e)到5(f)界线上侧固相团簇的消逝, 及从5(g)到5(h)界线下侧固相团簇的形成.这些大量高频的微观过程组成了固-液界线涨落和质量输运的动力学机制.上述动力学机制与近期关于异质Al-Pb固-液界面层面内两相平衡界线动力学机制的研究发现相近^[47].两相分界线的涨落和局域结构说明其为粗糙型, 原理上这意味着我们可以应用毛细波涨落方法计算单分子层限域固-液界线的界线自由能, 从而得到更多关于限域相之间形核与相变速率的机理的定量认识, 我们将在后续研究工作中开展探索.

  图 5(e)到5(h)中, 我们可以直接观察到固相中一个点缺陷的传导过程, 用黑色的线圈标出.通过检查, 这个点缺陷产生于固-液界线的涨落, 形成后在固相区域内传导, 经过很长的时间后, 它会再扩散回界线区域.这一结果并不奇怪, 传统固-液界面体系中存在类似的缺陷形成与输运的行为.值得提出的是, 单分子层固-液两相平衡态的模拟实现, 提供了直观研究关于固-液界面空位缺陷动力学的思路.

  3.3   限域冰-水平衡界线热力学量分布函数

  细粒度氧原子的密度分布结果在图 7(a)中画出.坐标0点左侧为固相区域氧原子的排列呈现理想的周期性层状分布, 在小于－40 Å的固相区域内, 密度峰的高度和宽度密度均保持一致, 两个近邻峰之间密度为0, 说明限域冰中的水分子主体处于固体晶格上振动, 不会出现自由扩散.沿着x轴从－40到40 Å, 氧原子的密度波动信号的幅度逐渐减小, 收敛至液相密度.上述密度分布函数在固-液转变的过渡区域的变化行为与粗糙型单质固-液界面的原子密度分布情况相近. 图 7(b)中为氢原子的细粒度密度分布结果.与氧原子密度不同的是, 峰谷的数目比氧原子密度峰的数目多了一倍, 同时观察到两类不同的峰高错落分布.通过与(a)图比较, 较低的氢原子密度峰分布在氧原子两个近邻峰间密度低谷区域, 而较高的氢原子密度峰与氧原子峰的位置一致.这一结果表明氢原子比较规则的分布在氧原子的上下左右, 构成接近方形的氢键网络.

  图 7

  

  图 7.  单分子层限域冰-水相平衡体系平面内细粒密度分布函数, x轴0点对应界线位置. (a)氧原子密度; (b)氢原子密度

  Figure 7.  In-plane fine-grained density profile for the monolayer confined ice-water coexistence, the zero point of x corresponding to the time averaged position of the solid-liquid coexistence line. (a) Area density of oxygen atoms; (b) Area density of hydrogen atoms

  下载: 全尺寸图片 幻灯片

  与Algara-Siller等工作一样, 我们获得的限域冰并非严格理想的正方形晶格.另外, 我们在固相区没有发现氢原子0密度的区域.比较图 7(a)和7(b)中粒子的密度波动信号变化形式, 发现氢原子和氧原子密度分布间的区别.从固相到液相, 氢原子比氧原子更早表现出无序性, 密度峰谷位置的衰减跨越的空间范围更大.这个密度信号上的不完全同步, 对应着在限域冰-水固-液两相平衡体系中, 氧原子(或水分子质心)的结构有序化与水分子转动相关的氢键网络有序化之间存在不同步行为.

  粗粒度氧原子密度分布在图 8(a)中画出.相比而言, 粗粒度密度分布能够更加直观地展现固-液界线区域密度的变化趋势.可以清楚看到分布函数的左右两端, 固相与液相的达到水平的恒定平衡态密度.利用双曲正切函数拟合图中氧原子密度分布函数, 我们可以获得限域冰水固-液界线的宽度大概在7 nm, 大概相当于25~30倍的方形冰晶格常数.这个宽度的数值远远大于金属体系固-液界面的宽度(通常为5~8倍的晶格常数).其中的原因可以归结于当前限域体系两相平衡界线的团簇无序化的微观结构及大幅度的毛细波涨落.这个结果也提醒读者, 在应用本论文的模拟方法模拟二维固-液相平衡体系时, 需要特别注意L_x尺寸的选取.应该避免选择太小的L_x尺寸, 以防止在导致周期性盒子中两个固-液界线间出现强烈相互作用.

  图 8

  

  图 8.  单分子层限域冰-水相平衡体系平面内粗粒分布函数, x轴0点对应相平衡界线位置. (a)氧原子密度; (b)应力; (c)压强分量

  Figure 8.  In-plane coarse-grained profiles for different thermodynamic properties for the monolayer confined ice-water coexistence, the zero point of x corresponding to the time averaged position of the solid-liquid coexistence line. (a) Area density of oxygen atoms; (b) In-plane stress; (c) Three pressure components

  下载: 全尺寸图片 幻灯片

  平衡态固-液界面需要满足块体部分的固体与液体处在流体静压平衡状态.检验一个固-液界面的模拟是否合理, 应力分布函数的计算尤为重要. 图 8(b)展现了体系平面内的应力分布函数, 其中, 固相和液相区域应力为0, 这说明目前建立单分子层限域冰-水相平衡的模拟方法可靠、可行.

  在固-液界线区域, －20 Å到20 Å范围内, 应力分布函数呈现一个巨大的负值峰, 峰值达到－0.6 kbar.负号意味着单分子层限域冰和水的固-液界线处于切向压缩的状态.从压强分量的面内分布函数来看, x压强分量在界线处比y压强分量小, 意味着界线区域相对倾向于在x方向上扩大同时在y方向上缩小.这些力学分量的组合, 对应着相平衡界线上的发生的团簇有序-无序化的微观动力学过程. x和y压强分量在界线以外的区域均等于限域压强5000 bar. z压强分量为1.25 kbar左右, 这个压强由限域尺寸和水分子与限域墙间的相互作用决定.

  除了上述密度、压强分量、应力等热力学量的统计计算之外, 限域两相平衡态模拟还可以提供更多重要热力学性质的分布研究, 如输运、偶极弛豫等.本文将不作具体介绍.

  4.   结论

  本文采用分子动力学方法实现了限域下单分子层水和冰的两相平衡的分子模拟, 精确确定了限域压强5000 bar的两相平衡温度为402.5 K.本文还分别讨论了瞬时时刻和时间平均后得到的固-液界线处结构的特点.通过长时间统计平均, 计算了限域冰-水界线的细粒度粒子密度分布、粗粒度应力分布等界线热力学量的分布函数, 获得了模拟原子动态演化轨迹之外的多个两相平衡界线的热力学性质.此外, 各种粗粒度热力学量的固-液界线分布函数, 印证了限域单分子层冰-水两相精准地处于热力学平衡条件, 形成稳定地共存系统.

  低维限域水存在许多不同结构相, 至今为止, 关于不同限域结构相之间的平衡体系研究较少.本研究中, 拓展固-液两相平衡模拟技术到单分子层物质体系, 我们确实感受到了将该技术从三维拓展到二维体系的困难与挑战性.二维单分子层限域冰的“块体”部分也即其“表面”, 缺陷对晶格结构的破坏性影响相较于块体冰更大, 压强平衡条件很难获得, 二维体系界面(线)纳米尺度的涨落对界面(线)稳定性的影响程度比三维体系高一个数量级, 两相平衡的建立也受固-液界线涨落的严重影响.实现了单分子层体系的固-液平衡, 其它限域条件下的冰相结构与液相水的两相平衡理论上都将能够实现, 并且更加容易.应用本文固-液平衡模拟技术到不同受限条件的限域冰-水体系中, 仅需精确获得限域冰晶格常数与温度的精确变化关系.

  本论文所发展的低维限域固-液平衡模拟技术适用于各种水模型, 包括能够模拟氢键网络畸变行为的高级水分子模型(可极化水分子模型等), 得到科学价值更高的限域相平衡界面性质研究结果.本模拟技术同样适用于及其它限域极性分子体系.

  平衡限域冰-水平衡态界线的原子模拟的成功实现, 奠定了限域相平衡热力学与动力学性质计算的研究基础.例如, 基于两相平衡态的模拟结果, 我们可以进一步追踪界线的热涨落从而计算出两相平衡界线的吉布斯自由能, 将推动限域冰-水固化或熔化相变形核过程的定量理论研究.我们也能从两相平衡态的模拟体系出发, 进行线性响应区间熔化与凝固相变的非平衡态分子动力学模拟, 以获得限域水固-液相变关键动力学性质的定量计算.

• 1. [1]

     Koga, K.; Zeng, X. C.; Tanaka, H. Phys. Rev. Lett. 1997, 79, 5262. doi: 10.1103/PhysRevLett.79.5262

  2. [2]

     Bai, J.; Zeng, X. C. Proc. Nat. Acad. Sci. 2012, 109, 21240. doi: 10.1073/pnas.1213342110

  3. [3]

     Ferguson, A. L.; Giovambattista, N.; Rossky, P. J.; Panagiotopoulos, A. Z.; Debenedetti, P. G. J. Phys. Chem. 2012, 137, 144501. doi: 10.1063/1.4755750

  4. [4]

     Kumar, P.; Buldyrev, S. V.; Starr, F. W.; Giovambat-tista, N.; Stanley, H. E. Phys. Rev. E 2005, 72, 051503. doi: 10.1103/PhysRevE.72.051503

  5. [5]

     Qiu, H.; Guo, W. Phys. Rev. Lett. 2013, 110, 195701. doi: 10.1103/PhysRevLett.110.195701

  6. [6]

     Bai, J. C.; Angell, A.; Zeng, X. C. Proc. Nat. Acad. Sci. 2010, 107, 5718. doi: 10.1073/pnas.0906437107

  7. [7]

     Johnston, J. C.; Kastelowitz, N.; Molinero, V. J. Phys. Chem. 2010, 133, 283101.

  8. [8]

     Algara-Siller, G.; Lehtinen, O.; Wang, F.; Nair, R.; Kaiser, U.; Wu, H.; Geim, A.; Grigorieva, I. Nature 2015, 519, 443. doi: 10.1038/nature14295

  9. [9]

     李海兰, 贾玉香, 胡仰栋, 物理化学学报, 2012, 28, 573. doi: 10.3866/PKU.WHXB201112191Li, H. L.; Jia, Y. X.; Hu, Y. D. Acta Phys.-Chim. Sin. 2012, 28, 573. doi: 10.3866/PKU.WHXB201112191

  10. [10]

      Joshi, R. K.; Carbone, P.; Wang, F. C.; Kravets, V. G.; Su, Y.; Grigorieva, I. V.; Wu, H. A.; Geim, A. K.; Nair, R. R. Science 2014, 343, 752. doi: 10.1126/science.1245711

  11. [11]

      赵梦尧, 杨雪平, 杨晓宁, 物理化学学报, 2015, 31, 1489. doi: 10.3866/PKU.WHXB201506011Zhao, M. Y.; Yang, X. P.; Yang, X. N. Acta Phys.-Chim. Sin. 2015, 31, 1489. doi: 10.3866/PKU.WHXB201506011

  12. [12]

      孙怡然, 于飞, 马杰, 物理化学学报, 2017, 33, 2173. doi: 10.3866/PKU.WHXB201705312Sun, Y. R.; Yu, F.; Ma, J. Acta Phys.-Chim. Sin. 2017, 33, 2173. doi: 10.3866/PKU.WHXB201705312

  13. [13]

      Bai, J.; Angell, C. A.; Zeng, X. C. Proc. Nat. Acad. Sci. 2010, 107, 5718. doi: 10.1073/pnas.0906437107

  14. [14]

      Johnston, J. C.; Kastelowitz, N.; Molinero, V. J. Phys. Chem. 2010, 133, 283101.

  15. [15]

      Chen, J.; Schusteritsch, G.; Pickard, C. J.; Salzmann, C. G.; Michaelides, A. Phys. Rev. Lett. 2016, 116, 025501. doi: 10.1103/PhysRevLett.116.025501

  16. [16]

      Koga, K.; Tanaka, H. J. Phys. Chem. 2005, 122, 104711.

  17. [17]

      Zangi, R.; Mark, A. E. Phys. Rev. Lett. 2003, 91, 025502. doi: 10.1103/PhysRevLett.91.025502

  18. [18]

      Zhao, W. H.; Wang, L.; Bai, J.; Yuan, L. F.; Yang, J.; Zeng, X. C. Acc. Chem. Res. 2014, 47, 2505. doi: 10.1021/ar5001549

  19. [19]

      Frolov, T.; Mishin, Y. J. Phys. Chem. 2015, 143, 044706.

  20. [20]

      Berendsen, H. J. C.; Grigerat, J. R.; Straatsma, T. P. Chem. Soc. 1987, 91, 24.

  21. [21]

      Giovambattista, N.; Rossky, P. J.; Debenedetti, P. G. Phys. Rev. Lett. 2009, 102, 050603. doi: 10.1103/PhysRevLett.102.050603

  22. [22]

      Xia, X. Berkowitz, M. L. Phys. Rev. Lett. 1995, 74, 3193. doi: 10.1103/PhysRevLett.74.3193

  23. [23]

      Kimmel, G. A.; Matthiesen, J.; Baer, M.; Mundy, C. J.; Petrik, N. G.; Smith, R. S.; Dohnalek, Z.; Kay, B. D. J. Am. Chem. Soc. 2009, 131, 12838. doi: 10.1021/ja904708f

  24. [24]

      Yang, J.; Meng, S.; Xu, L.; Wang, E. Phys. Rev. Lett. 2004, 92, 146102. doi: 10.1103/PhysRevLett.92.146102

  25. [25]

      Magda, J. J.; Tirell, M.; Davis, H. T. J. Chem. Phys. 1986, 84, 2901.

  26. [26]

      Werder, T.; Walther, J. H.; Jaffe, R. L.; Halicioglu, T.; Koumoutsakos, P. J. Phys. Chem. B 2003, 107, 1345. doi: 10.1021/jp0268112

  27. [27]

      Hockney, R. W. ; Eastwood, J. W. Computer Simulation Using Particles, CRC Press, U. S., 1988, 55.

  28. [28]

      Yeh, I. C.; Berkowitz, M. L. J. Chem. Phys. 1999, 111, 3155. doi: 10.1063/1.479595

  29. [29]

      Plimpton, S. J. Comput. Phys. 1995, 117, 1. doi: 10.1006/jcph.1995.1039

  30. [30]

      Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. J. Comput. Phys. 1997, 23, 327.

  31. [31]

      Frenkel, D. ; Smit, B. Understanding Molecular Simulation, 2nd ed., Academic Press, New York, 2002.

  32. [32]

      Broughton, J. Q.; Gilmer, G. H. J. Chem. Phys. 1986, 84, 5749. doi: 10.1063/1.449883

  33. [33]

      Broughton, J. Q.; Gilmer, G. H. J. Chem. Phys. 1986, 84, 5759. doi: 10.1063/1.449884

  34. [34]

      Davidchack, R. L.; Laird, B. B. Phys. Rev. Lett. 2000, 85, 4751. doi: 10.1103/PhysRevLett.85.4751

  35. [35]

      Hoyt, J. J.; Asta, M.; Haxhimali, T.; Karma, A.; Napolitano, R. E.; Trivedi, R.; Laird, B. B.; Morris, J. R. MRS Bull. 2004, 29, 935. doi: 10.1557/mrs2004.263

  36. [36]

      Frolov, T.; Mishin, Y. Model. Simul. Mater. Sci. Eng. 2010, 18, 074003. doi: 10.1088/0965-0393/18/7/074003

  37. [37]

      Becker, C. A.; Hoyt, J. J.; Buta, D.; Asta, M. Phys. Rev. E 2007, 75, 061610. doi: 10.1103/PhysRevE.75.061610

  38. [38]

      Beckera, C. A.; Asta, M.; Hoyt, J. J.; Foiles, S. M. J. Chem. Phys. 2006, 124, 164708. doi: 10.1063/1.2185628

  39. [39]

      Benet, J.; MacDowell, L. G.; Sanz, E. J. Chem. Phys. 2014, 141, 034701.

  40. [40]

      Benet, J.; MacDowell, L. G.; Sanz, E. Phys. Chem. Chem. Phys. 2014, 16, 22159. doi: 10.1039/C4CP03398A

  41. [41]

      Andersen, H. C. J. Chem. Phys. 1980, 72, 2384. doi: 10.1063/1.439486

  42. [42]

      Puri, P.; Yang, V. J. Phys. Chem. C 2007, 111, 11776. doi: 10.1021/jp0724774

  43. [43]

      Liang, H. T.; Laird, B. B.; Asta, M.; Yang, Y. Acta Mater. 2018, 143, 329. doi: 10.1016/j.actamat.2017.09.059

  44. [44]

      Davidchack, R. L.; Laird, B. B. J. Chem. Phys. 1998, 108, 9452. doi: 10.1063/1.476396

  45. [45]

      Buta, D.; Asta, M.; Hoyt, J. J. Phys. Rev. E 2008, 78, 031605. doi: 10.1103/PhysRevE.78.031605

  46. [46]

      Morris, J. R. Phys. Rev. B 2002, 66, 144104. doi: 10.1103/PhysRevB.66.144104

  47. [47]

      Yang, Y.; Olmsted, D.; Asta, M.; Laird, B. B. Acta Mater. 2012, 60, 4960. doi: 10.1016/j.actamat.2012.05.016

• 

  图 1  模拟盒子的示意图

  Figure 1  Schematic diagram of the simulation box

  Single layer of water molecules are confined in a slit pore consists of two parallel rigid walls

  下载: 全尺寸图片 幻灯片
  

  图 2  限域单分子层冰-水两相平衡态模拟技术流程图

  Figure 2  Flow chart of the MD simulation technique for obtaining the equilibrium coexistence of the monolayer confined ice and water

  下载: 全尺寸图片 幻灯片
  

  图 3  单分子层限域水和冰的氧原子面密度随温度变化.侧向压强5000 bar, 限域尺寸0.65 nm.实线和虚线分别为升温与降温路径, 方块与圆圈分别代表限域冰与水的密度数据. “热滞回线”型的密度-温度变化规律意味着当前单分子层限域冰-水体系经历了一级(固-液)相变

  Figure 3  Oxygen area density as the function of temperature for the confined monolayer ice and water. Lateral pressure 5000 bar, pore size 0.65 nm. Solid line and dashed line refer to the heating and cooling path, respectively. The density hysteresis implies that the confined monolayer ice and water undergoes a first order (solid-liquid like) phase transition

  下载: 全尺寸图片 幻灯片
  

  图 4  单分子层限域冰-水相平衡分子动力学模拟分子坐标快照.压强5000 bar、温度402.5 K

  Figure 4  Snapshot of the single layer confined ice-water equilibrium interface obtained from a molecular dynamic simulation. Lateral pressure 5000 bar, under a coexistence temperature of 402.5 K

  下载: 全尺寸图片 幻灯片
  

  图 5  (a) 固-液界线区域分子坐标快照, 红球为氧原子、白球为氢原子; (b)图(a)中氧原子的坐标与结构序, 原子的颜色变化遵循色谱变化, 代表了结构序参量数值改变, 蓝色代表固相, 红色代表液相; (c), (d)固-液界线区域的涨落; (e)~(h)、四个相隔1 ps的连续快照, 黑色的线圈标出了一个点缺陷的扩散

  Figure 5  (a) Snapshot of the region near the solid-liquid coexistence line; (b) Oxygen atoms in panel (a) is colored coded based on the value of the order parameter, blue for solid and red for liquid; (c), (d) Fluctuation of the coexistence line due to the in-plane solid-liquid phase transition; (e)~(h) Four successive trajectories with a 1 ps time interval, the diffusion of a point defect is marked by black circle

  下载: 全尺寸图片 幻灯片
  

  图 6  (a) 限域单分子层冰-水两相平衡界线附近固、液相二维氧-氧径向分布函数. (b)固、液相限域水中氧原子结构序参量的统计分布概率

  Figure 6  (a) A comparison of the 2d O—O radial distribution function between the confined monolayer ice and water phases near the equilibrium crystal-melt interfaces. (b) Probability distributions of the structural order parameter for the oxygen atoms in the confined monolayer ice and water phases

  下载: 全尺寸图片 幻灯片
  

  图 7  单分子层限域冰-水相平衡体系平面内细粒密度分布函数, x轴0点对应界线位置. (a)氧原子密度; (b)氢原子密度

  Figure 7  In-plane fine-grained density profile for the monolayer confined ice-water coexistence, the zero point of x corresponding to the time averaged position of the solid-liquid coexistence line. (a) Area density of oxygen atoms; (b) Area density of hydrogen atoms

  下载: 全尺寸图片 幻灯片
  

  图 8  单分子层限域冰-水相平衡体系平面内粗粒分布函数, x轴0点对应相平衡界线位置. (a)氧原子密度; (b)应力; (c)压强分量

  Figure 8  In-plane coarse-grained profiles for different thermodynamic properties for the monolayer confined ice-water coexistence, the zero point of x corresponding to the time averaged position of the solid-liquid coexistence line. (a) Area density of oxygen atoms; (b) In-plane stress; (c) Three pressure components

  下载: 全尺寸图片 幻灯片
•