水、甲醇和乙醇液体微结构性质的Car-Parrinello分子动力学模拟

曾勇平 朱晓敏 杨正华

(扬州大学化学化工学院,江苏扬州225002)

水、甲醇和乙醇液体微结构性质的Car-Parrinello分子动力学模拟

曾勇平*朱晓敏 杨正华

(扬州大学化学化工学院,江苏扬州225002)

采用Car-Parrinello分子动力学(CPMD)方法分别研究了水、甲醇和乙醇的液体微结构性质.研究结果显示:在水、甲醇和乙醇三个体系中O…O径向分布函数曲线的第一个峰位置分别为0.278、0.276和0.275 nm;O…H径向分布函数曲线的第一个峰位置分别为0.178、0.176和0.177 nm.表明基团(氢原子、甲基、乙基)的差异对O…O第一个峰的位置影响很小.但基团的差异对径向分布函数峰高的影响却很显著,由水到乙醇第一个峰的高度逐渐变高.空间分布函数表明氧原子和氢原子在溶剂分子周围有取向地分布,这与径向分布函数所表现出尖锐的第一个峰相一致.氢键分布分析显示,水、甲醇和乙醇的平均氢键数分别为3.62、1.99和1.87,表明水形成了网状氢键结构,而甲醇、乙醇形成链状氢键结构.

水;甲醇;乙醇;CPMD;微结构

1 引言

质子溶剂在有机合成、生物化学及生物技术方面有着重要的应用,因此其结构及动力学性质方面的研究已经引起广泛关注.1-5例如,由于氢键网络的存在,水不仅作为生物、化学反应过程的反应物,而且也作为一种反应媒介大量存在于反应过程中.6-8水的结构及动力学方面的认识,对于了解其存在时的生物、化学反应机理有着重要的意义.甲醇是同时具有疏水基团(甲基)和亲水基团(羟基)的最小有机分子,与水分子相比,结构上的差异使其具有两性特征,例如导致甲醇-水混合物在渗透气化中的不同扩散机理的差异.相对于甲醇分子,乙醇分子的疏水基团更大,其对结构的影响值得探索,而关于这方面的研究也有一些报道.9-11另一方面,甲醇和乙醇作为燃料电池12-17的燃料,深刻认识电池两极放电过程中的质子传递现象,对于开发高性能燃料电池有着重要的理论价值.水、甲醇、乙醇这三种分子在亲、疏水基团上的差异可能导致它们表现出的不同结构及动力学特征.近来,研究者们应用中子和X射线散射,18,19红外和拉曼光谱20,21以及飞秒光谱22来研究溶剂化特征.而许多关于这些溶剂在结构及动力学性质方面的实验报道仍存在较大的争议,这些实验方法研究的内容主要集中在宏观性质方面,对微观结构性质的研究还很有限.而计算和理论方法已经被证明为可以给出见解的强有力工具,并用这些方法已经得到了一些微观方面的性质.经典的分子动力学(MD)计算能够预测凝聚系统的平衡和非平衡性质,但此种方法没有传达任何电子结构方面的信息,另外缺乏普适性的经验力场,导致其研究结果的差异.而从头算分子动力学可以很好地克服这些问题,其仅有的假设是对离子运动作经典力学的近似,在最高水平的电子结构上,完成系统基态电子性质的计算.由于其在模拟精确性和简化模型上的优势而在凝聚态和材料理论设计等方面得到了广泛应用.

本文以水、甲醇、乙醇作为不同亲、疏水基团的研究对象,采用从头算分子动力学研究其结构及氢键性质,系统地考察不同的亲、疏水基团(甲基和乙基)对液体微结构性质的影响.

2 计算方法与体系模型

CPMD是基于第一性原理的分子动力学方法,它将密度泛函和分子动力学结合,并最大限度地发挥了密度泛函的优势.本文采用Car-Parrinello分子动力学方法和局域密度近似赝势对水、甲醇和乙醇的结构性质进行研究.为了取较大的时间步长,水、甲醇、乙醇分子中的氢用同位数氘取代.C原子和O原子选用Martins-Troullier(MT)赝势,23计算赝势非局域部分的贡献采用Kleinman-Bylander方法.24H原子选用Von Barth-Car赝势.25这些赝势的选择在研究纯液体甲醇26,27中得到了较好的结果.平面波截断能量选取为70 Ry(1 Ry=0.5 hartree),采用广义梯度近似(GGA)方法,用BLYP交换关联函数完成密度泛函计算.虚拟电子质量选取800 a.u.以使系统保持在玻恩-奥本海默面,运动方程积分时间步长为5 a.u.(0.12 fs).

表1 模拟体系的分子数和尺寸Table 1 Number of molecules and size of cell of the simulated system

体系模型根据分子在300 K温度下的密度确定(盒子尺寸及分子数见表1).首先采用全原子力场(OPLS)进行分子动力学模拟平衡1 ns,最后的平衡构型作为CPMD运行的初始构型.初始构型在300 K温度下以Car-Parrinello动力学热化3 ps,然后平衡2 ps,最后在NVE系综下模拟运行25 ps,取后12 ps作结果分析.经典分子动力学计算在实验室自己搭建的簇计算平台上完成,所有的CPMD计算在ScGrid的deepcomp7000上完成计算.

3 结果与讨论

3.1 二聚体

为了验证计算方法的合理性,本文计算了水、甲醇和乙醇二聚体的结构参数和结合能(二聚体与孤立分子的能量之差)等方面的数据.水、甲醇、乙醇二聚体的几何结构见图1,其二聚体分子间O…H的距离及结合能列于表2.23,28-37

平面波基组的BLYP函数很好地重复了二聚体和络合物的结构参数,但与B3LYP结果比较,呈现低估的趋势;与MP2/TZ2P+dif报道的结构参数较为一致,结果表明本文选用的赝势和截断能量是合理的.由于烷基的吸电子效应,使得水二聚体分子间氢键的距离略大于甲醇和乙醇二聚体.烷基的吸电子效应使得氧原子上的电子密度减小,从而使甲醇和乙醇二聚体分子间氢键距离小于水.然而甲醇分子间氢键的距离又略小于乙醇,分析其原因,这是由于位阻效应所致.由于乙基的体积大于甲基,其空间位阻效应的影响超过了吸电子效应的影响.二聚体的结合能由水到乙醇逐渐变小,即三者稳定性顺序为:水二聚体＞甲醇二聚体＞乙醇二聚体.

图1 优化的二聚体Fig.1 Optimized geometries of the dimer(a)water;(b)methanol;(c)ethanol

表2 水、甲醇和乙醇分子单体和二聚体作用参数比较Table 2 Comparison of structural parameters of isolated molecules and dimers of water,methanol,and ethanol

3.2 径向分布函数

径向分布函数(RDF)是获得体系在原子水平上结构性质的一种重要方法.其定义为:

这里,dV=4πr2dr;dn是dV球壳内的原子数目;＜ρ＞为平均数密度,＜ρ＞=n/V.径向分布函数中第一个峰对应中心原子周围的第一个配位层,第一个峰下面的面积(即g(r)从零积分到第一个最小值处)乘以＜ρ＞即为第一配位层中的原子数目N(即配位数).其表达式为:

图2 水、甲醇和乙醇分子间径向分布函数(a)O…O和(b)O…H(羟基H)及配位数Fig.2 Radial distribution functions of the O…O and O…H(hydroxyl H)and coordination numbers for water, methanol,and ethanol

图2为水、甲醇和乙醇体系中O…O和O…H的径向分布函数及配位数.由图可以看出,径向分布函数第一个峰表现得都很显著,反映出分子的排列呈近程有序性.对比O…O径向分布函数曲线中第一个峰的位置,水中为0.278 nm,甲醇中为0.276 nm,乙醇中为0.275 nm,O…O径向分布函数第一个峰的位置相差约为0.002 nm.在O…H径向分布函数曲线中,水中第一个峰的位置为0.178 nm,甲醇中为0.176 nm,乙醇中为0.177 nm.O…O和O…H的径向分布函数分别约在0.27和0.17 nm处存在波峰,表明氢键的存在.其它的实验和理论方法38-41也得到了类似的结果,由此可以得出,基团(氢原子、甲基、乙基)的差异对峰的位置有较小的影响.但基团大小的差异对径向分布函数峰高的影响却很显著,从水、甲醇到乙醇,第一个峰的高度依次增大.峰高度大处分子出现的概率大,在无穷远处g(r)等于1,表现为远程无序性.峰高说明第一溶剂壳稳定性顺序为:乙醇＞甲醇＞水,其原因由三种溶剂的液体结构性质决定的.在液体水中,形成的氢键结构形式为网状结构.在液体甲醇和乙醇中,形成的氢键结构形式为线性结构,加入一个新分子到线性结构氢键体系中要比加入一个新分子到网状结构的氢键体系中要困难.对于甲醇和乙醇,这是由于乙基的体积大于甲基,位阻效应所致.

对O…O和O…H径向分布函数第一个最小值及原子的配位数汇于表3.液体水第一溶剂壳的配位数为4.01,液体甲醇和乙醇第一溶剂壳的配位数分别为2.00和1.97.由第一溶剂壳的配位数可以得到,甲醇分子和乙醇分子的配位数相近,水分子配位数大约是甲醇分子和乙醇分子配位数的两倍.在液体水、甲醇、乙醇相邻的两个分子之间较容易形成氢键,配位数在一定程度上也反映了分子周围的氢键数.

3.3 空间分布函数

径向分布函数(RDF)仅仅给出的是整个空间距离上平均结构信息,不能给出溶液结构详细的局部信息,以及同等距离上的粒子在不同位置分布密度的高低.为了获得局域分布情况,本文计算了原子的空间分布函数.

图3为第一溶剂壳中氧原子和氢原子分别在水、甲醇和乙醇分子周围的空间分布图.氧原子和氢原子在溶液分子周围有取向地分布,这与径向分布函数表现出的尖锐的第一个峰一致.对于水分子,在分布强度较大时,靠近H原子侧的氧原子分布范围更广.这是因为分布在氧原子上的孤对电子直接与参考水分子(即水分子平均结构)中的氢原子形成氢键,而氢原子均匀分布在参考水分子中氧原子孤对电子的两侧,直接与水分子中的氧原子形成氢键.

表3 水、甲醇和乙醇中分子间O…O和O…H径向分布函数第一个最小值及配位数Table 3 The first minima and coordination number in radial distribution functions of O…O and O…H for water, methanol,and ethanol

图3 第一溶剂壳原子空间分布函数Fig.3 Spatial distribution of atoms in the first solvation shellO atoms relative to a water(a),methanol(b),and ethanol(c) molecules and H atoms relative to a water(d),methanol(e), and ethanol(f)molecules

与水中氧原子的分布类似,甲醇的氧原子分布的两个区域有明显的不同,作为质子接受体的甲醇分子中的氧原子分布更为集中.而氢原子均匀分布在参考甲醇分子中氧原子孤对电子的两侧,这种分布形式与水中的分布一致.

乙醇与甲醇的分布有较大的区别,乙醇的分布区域范围更广,而甲醇的分布比较集中.分析其原因是甲醇分子的稳定形式为反式构象,且所有的反式构象形式都是相同的,因此其分布也较为集中.但是,乙醇分子的稳定形式为邻位交叉构象,乙醇分子的邻位交叉构象有两种形式,因此乙醇中氧原子和氢原子的分布会随羟基的转动(也就是两种稳定构象的变化)而平均集中分布在两个区域.由于羟基的指向在两个稳定构象之间变化,所以氧原子和氢原子的分布范围较广.

3.4 氢键分布分析

为了进一步说明水、甲醇和乙醇分子周围分子的分布形式,这里统计了三种分子周围的平均氢键数目.分子间是否存在氢键不仅在距离上有严格的限制,在角度上也有明确的限定,本文判定分子间是否存在氢键的标准为:对水,相连的两个分子间O…O的距离小于0.34 nm,θ(H－O…O)小于30°;对甲醇和乙醇,相连的两个分子间O…O的距离小于0.35 nm,θ(H－O…O)小于30°.42定义fn为氢键所占的百分数,n(n=0,1,2,3…)为氢键的个数,fn和每个分子周围氢键平均数目nHB见表4.

表4 分子所包含不同氢键数的比例及每个分子周围氢键的平均数目(nHB)Table 4 Percentages of molecules with different numbers of hydrogen bonds and mean number(nHB)of H-bonds per molecule

统计显示水分子周围有4个氢键的占57.2%,其次为3个氢键,占30.9%,没有出现孤立的水分子和有6个氢键的水分子,极少数水分子周围存在一个氢键,少部分的水分子周围有2个或者5个氢键,平均氢键数为3.62.在水中,水分子中的两个氢原子都可以参与形成氢键,这点与甲醇和乙醇分子不同.甲醇和乙醇分子周围有两个氢键的居多,分别占76.4%和76.6%,少部分有1个和3个氢键,没有出现4个以上氢键的甲醇和乙醇分子,但出现了孤立的甲醇和乙醇分子,甲醇和乙醇的平均氢键数分别为1.99和1.87.这与甲醇和乙醇分子本身的结构有关,在这两个分子中,只有一个羟基氢原子参与氢键的形成.这里得到的分子周围氢键数与上面第一溶剂壳的配位数相一致.其它文献43-45也得到了类似的氢键数及配位数.为了直观地说明分子周围的氢键数,图4给出了水、甲醇和乙醇体系中具有不同氢键数目的分子数随模拟时间变化图.

由图4可以看出,在本文的模拟时间中,没有出现孤立的水分子和分子周围有6个氢键的水分子.结果显示48个水分子中有28个左右的水分子周围有4个氢键,15个左右的水分子有3个氢键,具有2个和5个氢键的水分子较少,基本在1到5个之间变化,偶尔会有1到2个具有一个氢键的水分子出现.大部分的水分子周围有4个氢键也说明了液体水中的氢键形式为四面体的网状结构.

甲醇和乙醇周围的氢键数与水体系有着较大的不同.在模拟期间,没有出现4个氢键以上的甲醇和乙醇分子,大多数甲醇和乙醇分子(大约20个左右)周围具有2个氢键,这也表明了液体甲醇和乙醇中的氢键结构为链状结构.甲醇和乙醇中有1个氢键的分子为氢键链两端的分子,有3个氢键的分子是在链上出现的分支.同时出现了孤立的甲醇和乙醇分子,在甲醇中,孤立分子不超过2个,在乙醇中孤立分子不超过3个.为了更直观地了解甲醇和乙醇的氢键链状结构,图5给出了水、甲醇和乙醇的氢键结构图.

图4 水(a)、甲醇(b)和乙醇(c)中分别具有0-6个氢键的分子数随时间的变化Fig.4 Number of molecules with 0-6 hydrogen bonds for water(a),methanal(b),and ethanol(c)along the simulation run

水分子中的氧原子作为质子接受体可以接受不同水分子中的两个氢原子形成两个氢键,同样,同一个水分子中的两个氢原子作为质子给体可以与不同水分子中的氧原子形成两个氢键,其氢键形式为四面体网络结构.甲醇和乙醇分子中的氧原子作为质子接受体只接受其它一个分子中羟基上的氢原子形成一个氢键,羟基上的氢原子作为质子给体与其它分子中的氧原子形成一个氢键,因此其局部为链状氢键结构.由此也可以得出,甲醇和乙醇中不同的疏水基团(甲基、乙基)对氢键的链状结构影响较小.

图5 水(a)、甲醇(b)和乙醇(c)第一溶剂壳的氢键构型Fig.5 Hydrogen bonding around solvent molecule in water(a),methanol(b),and ethanol(c)

4 结论

采用CPMD方法计算了水、甲醇和乙醇的径向分布函数、配位数和氢键分布等性质.考察了基团极性对流体结构的影响.计算结果表明在径向分布函数中,基团(氢原子、甲基、乙基)的差异对O…O第一个峰的位置影响很小.但基团体积的差异对径向分布函数峰高的影响却很显著,由水、甲醇到乙醇,第一个峰的高度逐渐变高.对于水分子,在分布强度比较大时,靠近H原子侧的氧原子分布范围更广.甲醇的氧原子分布的两个区域有明显的不同,作为质子接受体的甲醇分子中的氧原子分布更为集中.而氢原子均匀分布在参考甲醇分子中氧原子孤对电子的两侧.乙醇的氧原子分布区域范围更广,产生此种现象的原因为乙醇分子的邻位交叉构象有两种形式,由于羟基的指向在两个稳定构象之间变化,乙醇中氧原子和氢原子的分布会随羟基的转动而平均集中分布在两个区域.氢键分析显示水形成了网状氢键结构,而甲醇、乙醇形成链状氢键结构.

(1) Susan,B.R.;Lawrence,R.P.;Gerhard,H.;Joel,D.K.;Richard, L.M.;Antonio,R.J.Am.Chem.Soc.2000,122,966.

(2)Ayala,R.;Martinez,J.M.;Pappalardo,R.R.;Muňoz-Paez,A.; Marcos,E.S.J.Phys.Chem.B 2008,112,5416.

(3) Rudolph,W.W.;Fischer,D.;Irmer,G.;Pye,C.C.Dalton Trans. 2009,33,6513.

(4)Yu,X.C.;Lin,K.;Hu,N.Y.;Zhou,S.G.;Liu,S.L.Acta Phys.-Chim.Sin.2010,26,2473. [余小春,林 珂,胡乃银,周晓国,刘世林.物理化学学报,2010,26,2473.]

(5) Liu,Y.;Wang,F.F.;Yu,C.Y.;Liu,C.;Gong,L.D.;Yang,Z.Z. Acta Phys.-Chim.Sin.2011,27,379.[刘 燕,王芳芳,于春阳,刘 翠,宫利东,杨忠志.物理化学学报,2011,27,379.]

(6) Savage,P.E.Chem.Rev.1999,99,603.

(7) Hass,K.C.;Schneider,W.F.;Curioni,A.;Andreoni,W.Science 1998,282,265.

(8) Staehelin,J.;Hoigne,J.Environ.Sci.Technol.1985,79,1206.

(9) Danten,Y.;Tassaing,T.;Besnard,M.J.Phys.Chem.A 2006, 110,8986.

(10) Saiz,L.;Guardia,E.;Padro,J.A.J.Chem.Phys.2000,113, 2814.

(11) van Erp,T.S.;Meijer,E.J.J.Chem.Phys.2003,118,8831.

(12) Chen,R.;Zhao,T.S.Electrochem.Commun.2007,9,718.

(13) Neburchilov,V.;Martin,J.;Wang,H.;Zhang,J.J.Power Sources 2007,169,221.

(14)Zhou,W.J.;Zhou,B.;Li,W.Z.;Zhou,Z.H.;Song,S.Q.;Sun, G.Q.;Xin,Q.;Douvartzides,S.;Goula,M.;Tsiakaras,P. J.Power Sources 2004,126,16.

(15) Li,W.;Liang,C.;Zhou,W.;Qiu,J.;Zhou,Z.;Sun,G.;Xin,Q. J.Phys.Chem.B 2003,107,6292.

(16) Zhou,W.;Zhou,Z.;Song,S.;Li,W.;Sun,G.;Tsiakaras,P.;Xin, Q.Appl.Catal.B:Environ.2003,46,273.

(17) Song,S.;Tsiakaras,P.Appl.Catal.B:Environ.2006,63,187.

(18) Megyes,T.;Radnai,T.;Grósz,T.;Pálinkás,G.J.Mol.Liq. 2002,101,3.

(19) Herdman,G.J.;Salmon,P.S.J.Am.Chem.Soc.1991,113, 2930.

(20) Guillot,B.;Marteau,P.;Obriot,J.J.Chem.Phys.1990,93, 6148.

(21) Pye,C.C.;Rudolph,W.W.J.Phys.Chem.1998,102,9933.

(22)Omta,A.W.;Kropman,M.F.;Woutersen,S.;Bakker,H.J. Science 2003,301,347.

(23) Sprik,M.;Hutter,J.;Parrinello,M.J.Chem.Phys.1996,105, 1142.

(24) Wernet,P.;Nordlund,D.;Bergmann,U.;Cavalleri,M.;Odelius, M.;Ogasawara,H.;Näslund,L.Å.;Hirsch,T.K.;Ojamäe,L.; Glatzel,P.;Pettersson,L.G.M.;Nilsson,A.Science 2004,204, 995.

(25)Haughney,M.;Ferrario,M.;McDonald,R.J.Phys.Chem. 1987,91,4934.

(26) Saiz,L.;Padró,J.A.;Guàrdia,E.J.Phys.Chem.1997,101,78.

(27) Pagliai,M.;Cardini,G.;Righini,R.;Schettino,V.J.Chem. Phys.2003,119,6655.

(28) Handgraaf,J.W.;Erp,T.S.;Meijer,E.J.Chem.Phys.Lett. 2003,367,617.

(29) Kim,K.;Jordan,K.D.J.Phys.Chem.1994,98,10089.

(30) Provencal,R.A.;Casaes,R.N.;Roth,K.;Paul,J.B.;Chapo,C. N.;Saykally,R.J.J.Phys.Chem.A 2000,104,1423.

(31) Xu,X.;Goddard,W.A.,III.J.Phys.Chem.A 2004,108,2305.

(32) Benedict,W.S.;Gailan,N.;Plyler,E.K.J.Chem.Phys.1956, 24,1139.

(33) Boyd,S.L.;Boyd,R.J.J.Chem.Theory Comput.2007,3,54.

(34) Sasada,Y.;Takano,M.;Satoh,T.J.Mol.Spectrosc.1971,38,33.

(35) Culot,J.P.Symposium on Gas Phase Molecular Structure,4th ed.;Austin:Texas,1972;paper T8.

(36) Han,G.Z.;Zhang,C.;Gao,J.G.;Qian,P.Acta Phys.-Chim. Sin.2011,27,1361.[韩光占,张 超,高吉刚,钱 萍.物理化学学报,2011,27,1361.]

(37) Coussan,S.;Bouteiller,Y.;Perchard,J.P.;Zheng,W.Q. J.Phys.Chem.A 1998,102,5789.

(38) Narten,A.H.;Levy,H.A.J.Chem.Phys.1971,55,2263.

(39) Narten,A.H.;Nabenschuss,A.J.Chem.Phys.1984,80,3387.

(40) Svishchev,I.M.;Kusalik,P.G.J.Chem.Phys.1994,100,5165.

(41) Kuo,W.;Mundy,J.;McGrath,J.;Siepmann,J.I.;Van de Vondele,J.;Sprik,M.;Hutter,J.;Chen,B.;Klein,M.L.; Mohamed,F.;Krack,M.;Parrinello,M.J.Phys.Chem.B 2004, 108,12990.

(42) Padró,J.A.;Saiz,L.;Guàrdia,E.J.Mol.Struct.1997,416,243.

(43) Jorgensen,W.L.J.Phys.Chem.1986,90,1276.

(44) Boese,A.D.;Doltsinis,N.L.;Handy,N.C.;Sprik,M.J.Chem. Phys.2000,112,1670.

(45) Soper,A.K.;Bruni,F.;Ricci,M.A.J.Chem.Phys.1997,106, 247.

August 2,2011;Revised:September 2,2011;Published on Web:September 6,2011.

Car-Parrinello Molecular Dynamics Simulations of Microstructure Properties of Liquid Water,Methanol and Ethanol

ZENG Yong-Ping*ZHU Xiao-Min YANG Zheng-Hua
(College of Chemistry and Chemical Engineering,Yangzhou University,Yangzhou 225002,Jiangsu Province,P.R.China)

Car-Parrinello molecular dynamics(CPMD)calculations were performed on the solvent structure properties of water,methanol,and ethanol.The results show that the first O…O peaks of the radial distribution functions(RDFs)in the three systems are 0.278 nm for water,0.276 nm for methanol, and 0.275 nm for ethanol.The first O…H peaks of the radial distribution functions(RDFs)in the three systems are at 0.178 nm for water,0.176 nm for methanol,and 0.177 nm for ethanol.This indicates that the hydrophobic groupings(hydrogen,methyl,and ethyl)have little influence on the first peak position. However,the intensity of the RDFs increases from water to methanol and ethanol.The spatial distribution functions show that the oxygen and hydrogen atoms of other solvent molecules have characteristic orientations on the reference molecules in these systems.The results are in agreement with the first sharp peak of the radial distribution functions.We analyzed the hydrogen bonds using a statistical method.The results show that the average hydrogen bond numbers are 3.62 for water,1.99 for methanol,and 1.87 for ethanol.Therefore,different hydrogen-bonded network structures are formed for liquid water,methanol, and ethanol via hydrogen bonds.

Water;Methanol;Ethanol;CPMD;Microstructure

10.3866/PKU.WHXB20112779

∗Corresponding author.Email:ypzeng@yzu.edu.cn;Tel:+86-514-87975590-9109.

The project was supported by the National Natural Science Foundation of China(20806064).

国家自然科学基金(20806064)资助项目

O645;O641