质子水合结构的振动态密度分析

曾永辉，言天英

（南开大学材料科学与工程学院,新能源材料化学研究所,天津 300350）

质子的传递在化学反应中是一个重要的反应过程，从基本的酸碱化学反应到生物分子中能量的转换和存储、酶的催化以及质子燃料电池中的跨膜传递等一系列的重要理化和生物反应中都与其密切相关［1～6］.质子在水溶液中的迁移速率非常高，室温下，质子在水溶液中的迁移速率要普遍快于其它阳离子［7］.已有研究提出了一种机理用于解释质子在纯水中异常高的迁移速率，即质子可以沿着一条链状的水链从一个水分子跳跃到邻近的一个水分子上而不需要水分子大幅度的移动，这一过程涉及水合氢离子上氢氧共价键连续的断裂和生成，被称为Grotthuss 机理［8］.当一个质子被引入到纯水环境中时，其并不能独立存在于纯水中，而是迅速地与邻近的水分子结合形成一个水合氢离子.水合氢离子在空间上是一个sp3杂化的三角锥结构造型.它的3个氢原子与其第一溶剂化层的3个水分子以一种较强的氢键结合方式构成了一个独特的三配位水合结构，即Eigen［9］提出的H9O4+水合结构，质子上的多余电荷主要离域在这3个水分子上.1976年，Zundel等［10］提出了质子可能存在的另一种水合结构，他们认为质子是与两个水分子结合形成了一种二聚体结构（H5O2+），同理质子的多余电荷也主要分散在这两个水分子上.这两种结构被认为是质子在水溶液中传递所涉及的一系列相互转化结构的形态.

纯水中质子传递的Grotthuss机理主要涉及水合氢离子第一与第二溶剂化层氢键网络的变化，这些微观层面的结构变化往往是在几百个fs或几个ps内完成转换，在实验表征技术方面具有很大的挑战.现今在质子传递的研究方面，应用较多的实验技术主要是红外光谱，其为深入研究质子传递机理提供了许多实验信息.Thämer等［11］利用超快的二维红外光谱技术研究了质子在稀盐酸溶液中的扩散，他们发现Zundel的结构在质子传递过程中起到了重要的作用.Xu等［12］通过分子动力学模拟的方法解析了稀盐酸溶液的红外吸收光谱，针对稀酸溶液红外光谱中的连续宽吸收峰给出了解释，他们发现扭曲的Eigen结构在2000～3000 cm-1范围内存在一个宽吸收峰，而Zundel的结构在此区域内没有出现明显的吸收峰，这说明红外光谱中2000～3000 cm-1范围内的吸收峰主要是来自于扭曲的Eigen 结构的贡献.Bowman等［13］通过高阶的量化模拟，揭示了气相团簇［H+（H2O）4］下Eigen类型结构的光谱信号与实验观测到的气相团簇光谱信号相吻合的结果，其在2000～3000 cm-1范围内存在一个明显的宽吸收峰特征，将有助于对液态下质子的光谱信号的解析.Kulig 和Agmon［14］利用液体中团簇的计算方法对酸性水溶液中Eigen和Zundel结构以及溶剂化层各部分对红外光谱的贡献以及质子发生传递时的特征振动模式进行了相关研究.Fournier 等［15］通过模拟的方法对H+（H2O）n（n=2～28）团簇进行了红外光谱研究，他们发现随着团簇大小的增大，水合氢离子的溶剂化层结构可以扩展到第二溶剂化层以外并且“笼状”的氢键网络结构越加明显，以此来更好地模拟质子在“Bulk”水环境下的传递.Kuo等［16］在对阳离子状态下的气相水团簇（H2O）+n（n=3～11）的红外光谱研究中，发现水分子网络中的氢键强度是确定OH自由基相对于质子化位点（H3O+/H5O2+）的关键因素，该结果显示了水中H3O+与OH自由基接触对的不稳定性.

在酸溶液的实验光谱中，常存在一个宽的连续吸收谱带［11，17，18］，但由于在实验层面上很难对这些特征振动峰所涉及的结构动力学过程做出比较具体细致的解析，因此对于理解酸溶液中质子水合结构的动力学信息仍有赖于理论层面的深入考察.此外，质子的传递过程也是一个非常复杂的微观行为，其传递的机制与周围溶剂化层水分子氢键网络结构的动力学耦合作用仍然是一个需要深入探究的课题.为了揭示质子在水溶液中传递的微观结构及相应的动力学特征，本文通过多态经验价键模型（Multistate empirical valence bond model，MS-EVB）的方法，对所模拟的时间轨迹下的质子微观水合结构进行了划分，基于这些划分好的质子水合结构片段，对所需要考察的具体水合结构的特征振动以及所涉及的分子间氢键相互作用进行了比较细致的解析和讨论.

1 分子动力学模拟

计算技术的飞速发展为质子传递的研究提供了有力的辅助，一些宏观实验很难观测到的现象可通过计算机模拟的方法得以实现.在质子传递的研究方面，MS-EVB理论最早由Warshel 和Weiss［19］提出，通过不同的反应态来描述体系所进行的一系列化学反应过程.Vuilleumier 等［20］和Voth 等［21～24］分别发展了各自的多态经验价键模型，使得这一方法被广泛应用于质子在水溶液的传递研究中，本文也是基于此方法.简单来说，该方法主要的思路是将质子传递过程的反应划分成多个价键态，如图1 所示，一个Eigen的结构可以划分成4个不同的态，每个态的H3O+结构不同，且各个态出现的几率也不相同.

Fig.1 Description of the four EVB states in the Eigen structure in MS-EVB simulations

基于第三代的多态价键模型（MS-EVB3）［24］，通过一个自编译的分子动力学模拟程序来实现，模拟的体系包含216个浮动的水分子模型（SPC/Fw）［25］以及一个质子.体系是由一个边长为1.8621 nm，空间三维方向都施加了周期性边界条件的立方盒子构成的，立方盒子中的密度对应于1 g/cm3.分子间的长程静电相互作用通过Smooth particle mesh Ewald（SPME）［26］方法来处理.对于过剩电荷体系的模拟，Ewald的边界处理为导体边界，对应于一个无限大的介电常数，与Voth等［24］和Parrinello 等［27］的研究处理方法相同.体系平衡过程中，采用Nosé-Hoover恒温器［28］使体系的温度保持在300 K，粒子运动的方程通过Velocity Verlet 算法［28］来实现，积分步长设定在1 fs.在体系达到平衡后，去掉Nosé-Hoover恒温器使体系处于一个孤立的微正则（NVE）系综环境下，然后对粒子的运动轨迹进行采样收集.在此次模拟中，一共收集了一个长达10 ns的运动轨迹来对接下来的结构和动力学进行分析处理.在多态价键方法下，绝热态是每个非绝热EVB态的一个线性组合，（其中，ci是第i个态出现的概率幅度，第i个EVB态）.因此，体系中过剩电荷中心（Center of excess charge，CEC）的位置（rCEC）可以通过所有单个EVB态电荷中心的位置（rCOC）与其单个态所占的权重进行求解，即（其中，是第i个态出现的概率）.具有最大概率()的EVB态被看作是一个主态，这个主态结构在数据分析中也被看作是一个最接近水合氢离子的结构形态.

2 结果与讨论

2.1 质子的水合结构

在MS-EVB 的模型下，一个Eigen 的结构可以由前两个最大的EVB 态出现的几率来表示，；同理Zundel 的结构对应于这两个态的几率为≈0.45.换言之，当0.47 时，此时质子的水合结构可以看作一个Eigen，同理时，此时质子的水合结构是一个Zundel，因此可以通过计算的分布来更直观地显示质子在这两种水合结构的概率分布［图2（A）］.很显然，在MS-EVB的模型下，质子的水合结构更亲近于Eigen的结构形态.为了更全面了解质子的水合结构中水合氢离子和周围配位水分子的结构关系，一个包含氧氧距离的二维结构分布如图2（B）所示，其中氧氧的距离是指两个EVB 态的两个氧原子间的距离.从图2（B）可见，当0.42以及处于前两个最大态的状态下，H3O+的氧氧原子间的距离rO*O≈0.25 nm时，质子水合结构的分布最为显著，表明在Eigen 的结构中，主态水合氢离子的氧原子与周围配位水分子的氧原子的距离在0.25 nm左右.图2（C）给出了质子沿反应坐标δ的一个概率分布，反应坐标的定义为δ=rO1xH*-rO*H*，其中O*表示水合氢离子上的氧原子，O1x是最邻近水合氢离子的水分子上的氧原子，H*表示传递的质子［图2（C）插图］.从反应坐标的分布来看，其分布的峰值处于|δ|≈0.03 nm，根据之前的研究［7］，质子的水合结构还可通过反应坐标的几何方式来划分，当|δ|＞0.02 nm时可以看作一个Eigen的结构；同理，当|δ|＜0.01 nm 时，此时质子的水合结构表现为一个Zundel 类型，很明显δ的分布也表明在MSEVB方法下，质子的水合结构也是以Eigen 结构的分布为主.图2（D）给出了质子反应坐标δ与氧氧原子距离rO*O1x的二维分布，需要说明的是，此处的氧氧距离不同于图2（B）中的rO*O，rO*O1x为水合氢离子上的氧原子与最邻近它的水分子上氧原子之间的距离.可见，当|δ|≈0.03 nm，以及rO*O1x≈0.242 nm时，质子水合结构的分布最为显著，与图2（B）相比也可以发现，rO*O1x的值要小于rO*O≈0.25 nm，因为此时第二主态所具有的水分子并不总是最邻近水合氢离子，但两者所显示的分布趋势一致.

Fig.2 Distribution of the difference between and of the two EVB states of the highest probabilities(A),two-dimensional probability distribution P(-rO*O)(B),distribution of P(δ) for the proton transport along the reaction coordinate δ(inset)(C)and two-dimentional probability distribution function P(δ,rO*O1x)(D)

首先，根据几何标准将水合氢离子第一溶剂化层的水分子进行数据分类［图3（A）］.按照3个最近配位水分子离其它氧原子的距离，分别标记为O1x（最近邻水分子）、O1y（稍远的水分子）和O1z（最远的水分子），为了方便叙述，O代表一个水分子标识，而不纯粹表示氧原子.图3（B）展示了3个配位水分子与水合氢离子的径向分布函数（RDF）以及总体上水合氢离子与周围配位水分子平均的RDF.可见，水合氢离子与3个水分子的平均RDF分布峰值在0.25 nm处，与图2（B）的结果一致，此峰值的出现可能会认为是由H3O+第一溶剂化层中的3个等效的水分子引起的，但事实并非如此.通过单独对这3个水分子的RDF 分析发现，在这个平均的RDF 曲线下会裂解为3 个小的RDF 子集，对应的峰值分别位于0.241，0.251 和0.261 nm 处，即周围的3 个水分子以0.01 nm 的距离增长远离中心的H3O+.很显然有一个水分子更接近于H3O+，即O1x，这个水分子与H3O+的组合被称为特殊分子对（Special pair，SP）［7］.由此可知，Eigen 的结构是一个由中心H3O+与距离互不相等的3 个配位水分子构成的一个不对称扭曲结构，其中包含了一个“SP”，可能使得中心的H3O+结构状态更趋向于稳定；而对于另一种质子水合结构Zundel，可以被视为在质子的传递过程中的过渡态［7］.图3（C）给出了水合氢离子周围溶剂水分子的空间分布情况，在一个氢键长度（0＜rOO＜0.35 nm）的距离范围内，3个水分子以水合氢离子为中心有序地分布在其周围，表明水合氢离子上的3个氢原子与其周围的水分子较为紧密地结合在一起，而在其氧原子的孤对电子区域几乎没有水分子的分布，这种独特的属性使水合氢离子成为一种最小的两性分子［29］.

Fig.3 Configuration of the three coordination waters in the first shell of hydronium(A),partial and total radial distribution functions of the oxygen on the three coordinating waters of the hydronium(B)and spatial density distribution of oxygen(red) and hydrogen(white) of the coordinating water in around hydronium(C)

2.2 Eigen和Zundel的时间演化

为了考察不同状态下质子周围的溶剂化环境，整个10 ns的轨迹主要被划分为Eigen和Zundel的时间片段.在模拟中，记录轨迹每个时刻下最具有水合氢离子形态的氧编号O*以及其最近邻水分子的氧编号O1x，如图4（A）和（B）所示.在质子的水合结构为一个Eigen 的状态时，此时水合氢离子上氧原子的编号是不随时间变化的，而最邻近它的水分子O1x，根据Agmon等［7］的研究可知，这一阶段是质子传递的初期，“SP dance”，水合氢离子第一溶剂化层的3 个配位水会快速地交换其特殊分子对（SP）的身份，即O1x的编号会以2～3个不同编号的变化出现，将大于100 fs的这段时间归类为质子所处的Eigen时间片段［图4（A）］.同理，当质子的水合结构处于一个Zundel的状态下时，此时质子是在毗邻的两个水分子之间快速地来回穿梭，对于主态下的H3O+，其O*的编号是在两个原子编号下随时间不断地切换，同样，将大于100 fs的这段时间归类为一个Zundel的时间片段［图4（A）］.

Fig.4 Evolution of the identity of O*(A)and O1x(B),as atomic index in the simulation,from a trajectory segment of 7 ps and evolution of the reaction coordinate(δ)(C)

通过反应坐标δ可将质子的水合结构划分为Eigen和Zundel的结构，即|δ|＞0.02 nm时为一个Eigen结构，如图4（C）中主要集中在黑色虚线上方的片段，在|δ|＜0.01 nm时为Zundel结构，即图4（C）中位于两条绿色虚线之间的片段.通过反应坐标可以检验上述划分的结构片段是否合理，从整个时间段来看，这些划分的Eigen和Zundel结构片段大部分符合δ结构的划分标准，当然划分的这些结构片段并不都是包含各自理想的Eigen或Zundel结构，在划分的Eigen片段中也有一部分Zundel结构，同样Zundel的片段中也包含一些Eigen结构.

纯水中质子传递的Grotthuss机理主要涉及水合氢离子溶剂化层氢键网络的变化，这些微观层面的结构变化往往是在几百个fs或几个ps内转换完成，在实验表征方面具有很大的挑战.采用振动的态密度（Density of states，DOS）来分析质子水溶液体系的振动响应，图5（D）的黑色实线为包含一个质子和216个水分子体系中所有氢原子的DOS频谱，与实验下1 mol/L稀盐酸以及通过MS-EVB3方法模拟的红外光谱（Infrared spectrum，IR）相比［14］，体系的主要特征峰得到了很好的再现，但对于稀酸溶液的宽吸收特征，这部分信息有所掩盖，因为仅仅存在一个质子，其对整个体系的扰动十分微弱，所以整个体系的频谱主要是来自于水分子的贡献，这也是实验中所探测的稀酸溶液的红外光谱与纯水谱图相近的原因.因此，为了更深入了解水溶液中质子对周围溶剂化层水分子的扰动行为，首先对不同的质子水合环境做了划分.图5（A）～（C）展示了不同质子水合环境下的态密度谱，在分析各时间片段下的响应信息时，如对质子扩散行为的研究，同样选取过剩电荷中心（CEC）来计算各片段下的态密度.将图5（B）和（C）的谱线对比可见，CEC 在整个时间段下的态密度与在Eigen 片段下的态密度一致，这是因为在MS-EVB 中，整个体系主要以Eigen 的结构为主.从各自DOS 峰的主要特征来看，整体以及Eigen 片段的DOS中，在2000～3000 cm-1范围内有一个很明显的连续宽吸收谱带特征，这一响应特征也被归为来自于一个扭曲的Eigen结构的贡献［12］，而对于Zundel的结构片段在这一波数区间没有明显的吸收峰出现［图5（A）］.Bowman等［30］对气相团簇下的水合质子的研究揭示了红外光谱的特征振动峰与质子水合结构的状态密切相关，随着反应坐标δ值的减小，质子的水合结构也就越接近一个Zundel的结构，此时2000～3000 cm-1范围内的宽吸收特征峰信号也越来越弱.虽然体系CEC的DOS总体上反应出了质子所处环境的一个动态响应行为，但对于质子水合结构较为全面的动态结构信息特别是质子传递的一个特征响应却仍然不是很清晰，为此将在以Eigen和Zundel的时间片段下，分别只对它们各自的Eigen结构（H9）和Zundel结构（H5）进行了DOS分析.

Fig.5 Density of states(DOS) of the center of excess charge(CEC) from the Zundel(A),Eigen(B) and total segments(C),DOS of the H+(H2O)216 aqueous systems(D)

2.3 Eigen和Zundel的振动态密度

对于酸溶液，实验光谱的研究结果显示在2000～3000 cm-1范围内常可以见到一个宽的连续吸收谱带［11，17，18］.但是对于这部分特征振动峰的来源，从实验层面上还未给出具体的解析.通过此次划分出的Eigen时间片段，质子的水合结构大部分为一个扭曲的Eigen（），因此可以在此片段下以为对象，对以Eigen的质子水合结构进行更为细致的动力学研究.图6（A）给出了H9O4+以及中心H3O+上的所有氢原子的DOS，从两者谱图的特征来看，与实验观察到的宽的连续吸收谱带的结果相吻合，它们在2000～3000 cm-1范围内均具有一个明显的宽吸收特征峰.通过对比，可以推断出这部分的贡献应该是来自于Eigen中心的H3O+上氢原子的伸缩振动，更为确切的说是H3O+上两个与O1y和O1z氢键相结合的氢原子的伸缩振动特征.因为对于H3O+第一溶剂化层的水分子，均以一个比较强的氢键和H3O+上的氢键合，强氢键效应必然会削弱它的OH振动，使它的伸缩振动峰出现一个红移（最大峰值在2840 cm-1处），可从的DOS 上对比看出H3O+周围的配位水分子的OH 振动，如图6（A）插图中H2O 上所标注的Stretch（Str.），其特征峰在3560 cm-1附近非常显著.在1680 cm-1处，两者也各有一个显著的特征峰，归属于来自H3O+上可能出现质子传递行为的贡献，即图6（A）插图中H3O+上所标注的Str.，相比于O1y和O1z，它的OH振动红移现象更为明显，因为质子与最近邻的水分子O1x具有一个更强的氢键作用.在此需要说明的是，此时Eigen中心的H3O+，其质子并没有和周围的水分子进行实际的传递，所以这个峰值不能简单看作为质子传递的特征峰.在1520 cm-1处，的DOS出现了一个特征峰，而H3O+的DOS并未出现，因为其是水合氢离子周围配位水分子的弯曲振动［图6（A）插图中H2O上所标注的Bend］所引起的.

Fig.6 DOS of the Eigen structure(red) and H3O+ cation(black) taken from the Eigen segements of the trajectory(A),similar spectra for the Zundel structure(red) and H3O+ cation(black) taken from the Zundel part of the trajectory(B)

3 结 论

以质子的水溶液体系为研究对象，采用多态经验价键（MS-EVB）方法，对质子的水合结构以及溶剂化层结构进行了相应的分析研究.在基于MS-EVB 的模拟方法下，水溶液中质子的水合结构主要以Eigen和Zundel的结构为主，其中作为一个过渡态的Zundel结构在模拟中所占的比率并不高，更多的质子水合结构以Eigen的形式出现.为了更细致地探讨在不同质子化的水合环境下质子水合结构的动态响应特征，整个模拟轨迹被划分为以Eigen 和Zundel 为主的两个时间片段.对质子主要的Eigen 和Zundel 水合结构所具有的特征振动响应进行了分析，从Eigen 结构的振动特征峰来看，在2000～3000 cm-1范围内具有一个明显的连续宽吸收谱带，这些特征峰的出现与水合氢离子第一溶剂化层内的强氢键作用密切相关.相反，对于Zundel结构，并未出现明显的特征峰，但在1760 cm-1处的肩峰归属为质子传递的特征振动，与之前相关研究结果相符.