水分子間氫鍵判別標(biāo)準(zhǔn)的從頭分子動(dòng)力學(xué)模擬研究

陳 靜，潘 章，紀(jì) 青

（1.唐山學(xué)院 基礎(chǔ)部，河北 唐山 063000；2.河北工業(yè)大學(xué) 理學(xué)院，天津 300401）

0 引言

水體系是典型的氫鍵體系，液態(tài)水體系表現(xiàn)出很多特殊的熱力學(xué)性質(zhì)都起源于水分子間的氫鍵網(wǎng)絡(luò)[1-7]．目前還不能通過(guò)物理方法直接觀測(cè)到水中的氫鍵[5]，與氫鍵性質(zhì)相關(guān)的研究主要是采用理論計(jì)算方法和計(jì)算機(jī)模擬方法．在這些研究之前必須明確氫鍵的定義和判別標(biāo)準(zhǔn)[4]．目前有很多種氫鍵的判別方法，如結(jié)合能判據(jù)、幾何判據(jù)、電荷轉(zhuǎn)移判據(jù)、拓?fù)浞椒ㄅ袚?jù)等等．文獻(xiàn)中常用的氫鍵判別方法是結(jié)合能判據(jù)和幾何判據(jù)．結(jié)合能判據(jù)是根據(jù)水分子間結(jié)合能的分布情況人為選擇一個(gè)能量截?cái)嘀底鳛槟芰颗袚?jù)，當(dāng)兩個(gè)水分子的相互作用能量低于截?cái)嘀禃r(shí)就認(rèn)為它們形成了氫鍵，高于截?cái)嘀禃r(shí)就認(rèn)為不形成氫鍵．通常采用的能量截?cái)嘀禐?10kJ/mol[4，8-9]．幾何判據(jù)是通過(guò)人為選擇一個(gè)距離截?cái)嘀祷蛘咭粋€(gè)角度截?cái)嘀底鳛榕袆e是否形成氫鍵的標(biāo)準(zhǔn)．現(xiàn)在廣泛應(yīng)用的幾何判據(jù)要考慮一個(gè)距離（如R或者r）以及一個(gè)角度（如 ， ， ， 和 ）（見圖1，其中 軸為氫鍵受體分子的角平分線， 軸與氫鍵受體分子平面垂直），被稱為“距離—角度”幾何標(biāo)準(zhǔn)[4，7，10-15]．距離和角度的截?cái)嘀荡蠖喔鶕?jù)液態(tài)水的徑向分布函數(shù)來(lái)選取（見圖2）．徑向分布函數(shù)是可以由實(shí)驗(yàn)直接測(cè)量的量，它定義為距離任意一個(gè)給定粒子的一定距離處的粒子密度與平均粒子密度之比．這個(gè)函數(shù)在距離很大時(shí)必然是趨于1的，但是在距離很小時(shí)，它的取值能反映一個(gè)中心粒子周圍其它粒子與中心粒子的相對(duì)位置關(guān)系．圖2給出了液態(tài)水分子的氧原子間的徑向分布函數(shù)．可以看出，在2.5 < <3.3的范圍內(nèi)出現(xiàn)了第一個(gè)主峰，它表明，處在一個(gè)水分子周圍的第一水合層中的水分子都處在這個(gè)范圍內(nèi)．在這個(gè)范圍內(nèi)的水分子通常都會(huì)與中心水分子形成氫鍵，所以有的氫鍵判別標(biāo)準(zhǔn)將氧原子間距落在這個(gè)范圍內(nèi)作為兩個(gè)水分子形成氫鍵的幾何判據(jù)之一．幾何判據(jù)在統(tǒng)計(jì)分析和軟件顯示上存在著很大的優(yōu)勢(shì)，所以現(xiàn)在很多文獻(xiàn)中均采用幾何定義來(lái)分析體系中的氫鍵特點(diǎn)．事實(shí)上，任何氫鍵的定義中的截?cái)嘀刀蓟蚨嗷蛏俚拇嬖谥藶樾裕跓徇\(yùn)動(dòng)環(huán)境下，水分子間的氫鍵不斷的形成和打破．有些水分子會(huì)相距很近，實(shí)際上已經(jīng)處于排斥狀態(tài)，其相互作用能量遠(yuǎn)大于 10kJ/mol．但它們之間的幾何關(guān)系依然可以滿足氫鍵的幾何判據(jù)，還會(huì)被判定為形成了氫鍵．也就是說(shuō)一對(duì)水分子按照能量判據(jù)不能夠滿足氫鍵條件，但是按照幾何標(biāo)準(zhǔn)卻依然可能被判定為形成了氫鍵．在確定氫鍵的判別標(biāo)準(zhǔn)時(shí)需要合理地處理這種明顯然不合理的情況．

關(guān)于液態(tài)水的氫鍵體系的大量研究是采用分子動(dòng)力學(xué)模擬方法進(jìn)行的．這種方法能夠模擬的水體系可以包含多達(dá)幾十萬(wàn)個(gè)水分子，在研究體系的統(tǒng)計(jì)行為時(shí)具有巨大的優(yōu)勢(shì)．在分子動(dòng)力學(xué)模擬中，水分子的力學(xué)狀態(tài)和相互關(guān)系都是按照經(jīng)典力學(xué)方法描述的．然而氫鍵的形成有來(lái)自供體-受體分子間的電荷轉(zhuǎn)移的貢獻(xiàn)．電荷轉(zhuǎn)移是一種典型的量子效應(yīng)，無(wú)法用經(jīng)典力學(xué)方法準(zhǔn)確地描述．而且在常溫環(huán)境下，氫鍵不斷的形成和打破，水分子間的相互作用能也在不斷的變化．因此，對(duì)常溫下液態(tài)水體系的模擬研究采用基于第一性原理的分子動(dòng)力學(xué)方法更加合適[16-17]．基于密度泛函理論的從頭分子動(dòng)力學(xué)方法（Car-Parrinello ab initio molecular dynamics，CPMD）早在20世紀(jì)90年代就開始用于研究液態(tài)水體系，從那時(shí)開始很多應(yīng)用CPMD方法研究水體系的研究工作被相繼報(bào)道[18-19]．目前研究表明，用色散修正后的BLYP交換關(guān)聯(lián)函數(shù)進(jìn)行CPMD模擬，得到的模擬結(jié)果與實(shí)驗(yàn)結(jié)果符合的很好[20-21]．本文用從頭分子分子動(dòng)力學(xué)CPMD方法對(duì)水分子在熱運(yùn)動(dòng)環(huán)境下氫鍵的特征進(jìn)行系統(tǒng)研究，并與分子動(dòng)力學(xué)方法給出的結(jié)果進(jìn)行比較，從而對(duì)現(xiàn)有的氫鍵幾何判據(jù)進(jìn)行合理的改進(jìn)．

1 模擬參數(shù)與計(jì)算細(xì)節(jié)

用CPMD方法中的限定性動(dòng)力學(xué)方法對(duì)水二聚體進(jìn)行模擬．立方體盒子的邊長(zhǎng)為 =18.8973a.u.．模擬采用Becke-Lee-Yang-Parr的交換關(guān)聯(lián)泛函和范數(shù)不變的Troullier Martins贗勢(shì)近似．Kohn-Sham軌道在平面波基組下展開，能量截?cái)酁?0 Ry，虛擬電子質(zhì)量為600 a.u.，時(shí)間步長(zhǎng)設(shè)置為5 a.u.（0.121 fs），溫度為330(±20)K．模擬過(guò)程中限定了二聚體O-O間的距離 和角度 ，每一個(gè)構(gòu)象的模擬產(chǎn)出時(shí)間為2.4 ps．幾何參數(shù)和坐標(biāo)系參見圖1，圖中 軸是氫鍵受體分子的角∠HOH的平分線， 軸與氫鍵受體分子所在平面垂直， 軸與連線O…H的夾角定義為角 ．繞 軸， 軸， 軸旋轉(zhuǎn)的角度分別定義為 ， 和 ．優(yōu)化后的水二聚體幾何參數(shù)為O…O=2.97 ，O…H=1.99 ， =32°以及 =169.42°．為了系統(tǒng)的計(jì)算水二聚體在不同構(gòu)象下的相互作用能，通過(guò)限定性動(dòng)力學(xué)方法[22]將距離 和角度 限定在上述值．同時(shí)，通過(guò)限定供體氫原子和兩個(gè)受體氫原子間的距離實(shí)現(xiàn)受體水分子繞 軸， 軸， 軸的旋轉(zhuǎn)．繞其中的一個(gè)軸旋轉(zhuǎn)時(shí)，該軸被固定在優(yōu)化后的位置，其他兩軸隨著水分子取向變化發(fā)生變化，所有的旋轉(zhuǎn)都起始于最初的優(yōu)化得到的構(gòu)象．文中給出的水二聚體間的相互作用能是二聚體分子在有限溫度下的結(jié)合能平均值．值得注意的是，按上述方法限定距離 和角度 后，水二聚體總是符合通常形成氫鍵的幾何標(biāo)準(zhǔn)的．但是下面會(huì)看到，當(dāng)兩個(gè)水分子的相對(duì)取向發(fā)生變化時(shí)，其相互作用能量有可能超過(guò)氫鍵的能量標(biāo)準(zhǔn)，因而并沒有真正形成氫鍵．

圖1 水二聚體相關(guān)距離和角度示意圖Fig.1 Pictorial representation of thedistances and and angles( , , , and )for water dimer

圖2 應(yīng)用CPMD計(jì)算水分子體系平衡時(shí)的氧-氧徑向分布函數(shù)Fig.2 Oxygen-oxygenradial distributionfunction（RDF）obtained from the CPMD simulation

2 計(jì)算結(jié)果

用CPMD的限定性動(dòng)力學(xué)方法計(jì)算了在2個(gè)水分子的幾何關(guān)系滿足通常的氫鍵判別標(biāo)準(zhǔn)的情況下，其相互作用能量隨二者相對(duì)取向角度的變化．此外，還計(jì)算了二聚體相互作用能量隨受體氧與供體氫之間距離的變化．將計(jì)算結(jié)果與能量判別標(biāo)準(zhǔn)進(jìn)行比較可以明顯地看到通常的氫鍵幾何判別標(biāo)準(zhǔn)存在的問(wèn)題．

2.1 水二聚體間相互作用能與、、的關(guān)系

當(dāng)受體水分子繞 軸旋轉(zhuǎn)時(shí)，旋轉(zhuǎn)角度 是 軸與分子間連線O…H的夾角．順時(shí)針?lè)较蜣D(zhuǎn)動(dòng)時(shí) 角增加，逆時(shí)針?lè)较蜣D(zhuǎn)動(dòng)時(shí) 角減?。畧D3a）中給出了水二聚體間相互作用能隨著 變化的曲線，其中 從 90～90°．在圖3a）中，2個(gè)水分子間相互作用能量的最低點(diǎn)大約在

=30°處，此時(shí)的相互作用能約為 18 kJ/mol．當(dāng)大于 17°時(shí)，水分子間的相互作用能量高于 10kJ/mol，當(dāng) 在 90°～ 40°之間時(shí)，相互作用能甚至變成了正值，說(shuō)明此時(shí)2個(gè)水分子間沒有形成穩(wěn)定的氫鍵二聚體結(jié)構(gòu)，它們甚至處在相互排斥的狀態(tài)．值得指出的是，圖3a）中計(jì)算的水分子間的相互作用能是在O…O=2.97 的條件下得到的．O…O=2.97 是液態(tài)水體系氧-氧間徑向分布函數(shù)的第1個(gè)峰值．因此圖3中的相互作用能曲線可以代表第一水合層中的水分子與中心水分子的相互作用能量與取向角 的關(guān)系．

繞 軸轉(zhuǎn)動(dòng)的角度 是指 軸固定不動(dòng)， 軸繞著軸順時(shí)針轉(zhuǎn)動(dòng)時(shí)轉(zhuǎn)過(guò)的角度．圖3b）是水分子間相互作用能與 的關(guān)系．根據(jù)構(gòu)象的對(duì)稱性， 的變化范圍是0～90°．可以看到，當(dāng) 大于72°時(shí)，水分子間的相互作用能量高于 10 kJ/mol．

繞 軸旋轉(zhuǎn)的角度 是指 軸不動(dòng)， 軸繞著 軸順時(shí)針轉(zhuǎn)動(dòng)時(shí)轉(zhuǎn)過(guò)的角度．在繞 軸轉(zhuǎn)動(dòng)過(guò)程中，角度 始終等于32°．圖3c）展示了水分子間相互作用能與 的關(guān)系，由于構(gòu)象的對(duì)稱性， 的變化范圍是0～180°．從圖3c）可以看出當(dāng) 大于80°時(shí)，水分子間的相互作用能已經(jīng)大于 10 kJ/mol，甚至是變成了正值．

2.2 水二聚體間相互作用能量與rO…H的關(guān)系

圖3 水二聚體在限定距離 和角度 時(shí)不同構(gòu)象下相互作用能量，其中a），b），c）分別為水二聚體間相互作用能量隨角度 ， ， 的變化關(guān)系。Fig.3 Theinteraction energy of thewater dimer asa function of theangle , , with thedistance and theangle constrained

水二聚體中受體氧和供體氫之間的距離O…H（圖1中的 ）是判別氫鍵的最常用幾何量．關(guān)于O…H的判別標(biāo)準(zhǔn)存在許多差異．有些文獻(xiàn)中給出的氫鍵的幾何判據(jù)中考慮了O…H的上下限，例如“ -”判據(jù)中， 取值范圍被限定為1.59 < <2.27[11]．也有一些幾何判據(jù)中沒有考慮O…H的下限，如 “ - ”判據(jù)[4]．為了確定O…H對(duì)水分子間相互作用能量的影響，本文對(duì)水分子間相互作用能量與O…H的關(guān)系進(jìn)行了系統(tǒng)計(jì)算，計(jì)算時(shí) 和 保持不變．圖4展示了水分子間相互作用能與 的關(guān)系．當(dāng) <1.62時(shí)，水分子間的相互作用能已經(jīng)超過(guò)了 10kJ/mol，并且隨著 減小，相互作用能迅速增加．

圖4 水二聚體間相互作用能與 O…H的變化關(guān)系Fig.4 Theinteraction energy of thewater dimer as afunction of thedistance O…H

3 分析與討論

上面的二聚體相互作用能量是在2個(gè)水分子間的距離和角度 或角度 被限定在優(yōu)化后的數(shù)值上進(jìn)行計(jì)算的．按照傳統(tǒng)的氫鍵幾何標(biāo)準(zhǔn)，如“ - ”判據(jù)，判斷是否水分子間形成氫鍵，無(wú)論這2個(gè)水分子的相對(duì)位置怎樣變化，它們都被判斷為形成了氫鍵．但是上面對(duì)2個(gè)分子的相互作用能量的計(jì)算結(jié)果表明，隨著2個(gè)水分子間相對(duì)取向的變化，有些構(gòu)象下二者并沒有形成有效的氫鍵．例如，當(dāng)1） < 17°（圖3a））時(shí)，2）>72°（圖3b））時(shí)，3） >80°（圖3c））時(shí)，2個(gè)水分子間的相互作用能已經(jīng)大于 10 kJ/mol，甚至是正值．按照氫鍵的結(jié)合能標(biāo)準(zhǔn)，在這些構(gòu)象下二者并沒有形成氫鍵．通過(guò)仔細(xì)的分析這些構(gòu)象下的特點(diǎn)發(fā)現(xiàn)，這些構(gòu)象下兩個(gè)水分子間的氫原子間的距離H…H都明顯小于2.4 （這里2.4 是2個(gè)氫原子的范德華半徑之和）．圖5中顯示了當(dāng)轉(zhuǎn)動(dòng)角度處在前面指出的范圍內(nèi)的幾個(gè)特殊值時(shí)，相應(yīng)的H…H的取值，其中，當(dāng) = 17°時(shí)，H…H=2.03 （圖5a）），當(dāng) =72°時(shí)，H…H=2.05 （圖5b）），當(dāng) =80°時(shí)，H…H=2.02（圖5c））．在這些情況下，氫原子間會(huì)產(chǎn)生較強(qiáng)的靜電排斥作用．本文相信是氫原子間的靜電排斥作用導(dǎo)致了分子間氫鍵的破壞．圖3a）、圖3b）、圖3c）中曲線的變化趨勢(shì)都明顯的說(shuō)明了這一點(diǎn)．例如在圖3a）中，當(dāng) 90°< < 17°時(shí)，分子間的相互作用能逐漸的減弱甚至變成正值，與之相應(yīng)的是H…H從2.03減小到1.50 ．此時(shí)氫原子間的排斥作用會(huì)增加的非常顯著，排斥的作用導(dǎo)致了水分子間結(jié)合能的減小，這與圖3a）中能量曲線是一致的．圖3b）和圖3c）中的結(jié)合能減弱，甚至變?yōu)檎狄彩怯捎陬愃频脑蛟斐傻模肿娱g相互作用能量與O…H的關(guān)系也表明，當(dāng)O…H減小時(shí)靜電排斥作用大大削弱了分子間的氫鍵作用．所以在判斷是否形成氫鍵的幾何判據(jù)中必須要加入O…H的下限．

4 結(jié)論

本文應(yīng)用從頭分子動(dòng)力學(xué)軟件CPMD結(jié)合限定性動(dòng)力學(xué)的方法，對(duì)熱環(huán)境下水二聚體不同空間取向下的相互作用能量進(jìn)行了系統(tǒng)的計(jì)算．將計(jì)算結(jié)果與傳統(tǒng)的氫鍵幾何標(biāo)準(zhǔn)的判別結(jié)果進(jìn)行比較，發(fā)現(xiàn)水分子間氫鍵的幾何判斷標(biāo)準(zhǔn)需要考慮兩個(gè)新的幾何限制，即兩個(gè)水分子間氫原子的間距H…H需要有下限，同時(shí)受體氧原子和供體氫原子間的距離O…H也需要有下限．根據(jù)本文的計(jì)算結(jié)果，將H…H>2.0 ，O…H>1.62加入氫鍵幾何標(biāo)準(zhǔn)中可以使幾何標(biāo)準(zhǔn)與結(jié)合能標(biāo)準(zhǔn)（ 10kJ/mol）的判別結(jié)果符合的更好，從而可以有效的減少幾何判別方法在統(tǒng)計(jì)水體系氫鍵時(shí)產(chǎn)生的統(tǒng)計(jì)偏差．

圖5 旋轉(zhuǎn)受體水分子時(shí)兩個(gè)水分子氫原子間的最短距離（a）對(duì)應(yīng)于= 17°，b）對(duì)應(yīng)于 =72°，c）對(duì)應(yīng)于 =80°．為了顯示旋轉(zhuǎn)過(guò)程在d）中給出了水二聚體優(yōu)化后的初始構(gòu)象）Fig.5 Theshortest intermolecular H…H distances(between centersof mass)of thewater dimer when a） = 17°,b） =72°and c） =80 °.Theinitial stateof thewater dimer and theoriginal axes 0,0 and 0 of theacceptor water aredrawn in d）．

[1]Buckingham AD，Bene JED，McDowell SAC．Thehydrogen bond[J]．Chem Phys Lett，2008，463（1）：1-10．

[2]Li Jingyuan，Gong Xiaojing，Lu Hangjun，etal．Electrostatic gatingof ananometer water channel[J]．Proc Natl Acad Sci，2007，104（10）：3687-3692．

[3]Tanaka H．Simplephysical model of liquid water[J]．JChem Phys，2000，112（2）：799-809．

[4]Kumar R，Schmidt JR，Skinner L．Hydrogen bonding definitionsand dynamicsin liquid water[J]． JChem Phys，2007，126：204107．

[5]Fernández-Serra MV，Artacho E．Electronsand Hydrogen-Bond Connectivity in Liquid Water[J]．Phys Rev Lett，2006（96）：016404．

[6]Smith JD，Cappa CD，Wilson KR，etal．Energeticsof Hydrogen Bond Network Rearrangementsin Liquid Water[J]．Science，2004，306：851-853．

[7]Wernet P，Nordlund D，Bergmann U，et al．The Structureof the First Coordination Shell in Liquid Water[J]．Science，2004，304：995-999．

[8]Stillinger FH．Water revisited[J]．Science，1980，209：451-457．

[9]Jorgensen WL，Chandrasekhar J，Madura JD，et al．Comparison of simple potential functions for simulating liquid water[J]．JPhys Chem B，1983，79（2）：926-935．

[10]Luzara A．Resolving thehydrogen bond dynamicsconundrum[J]．JChem Phys，2000，113（23）：10663-10675．

[11]Kuo IFW，Mundy CJ．An ab Initio Molecular Dynamics Study of the Aqueous Liquid-Vapor Interface[J]．Science，2004，303：658-660．

[12]Todorova T，Seitsonen AP，Hutter J，et al．Molecular Dynamics Simulation of Liquid Water:Hybrid Density Functionals[J]．JPhys Chem B，2006，110：3685-3691．

[13]Luzar A，Chandler D．Effect of Environment on Hydrogen Bond Dynamicsin Liquid Water[J]．Phys Rev Lett，1996，76（6）：928-931．

[14]Luzar A，Chandler D．Hydrogen-bond kineticsin liquid water[J]．Nature，996，379：55-57．

[15]Zielkiewicz J．Structural properties of water:Comparison of the SPC,SPCE,TIP4Pand TIP5Pmodels of water[J]．JChem Phys，2005，123（10）：104501．

[16]Morozov AV，Kortemme T，Tsemekhman K，et al．Close agreement between the orientation dependence of hydrogen bonds observed in protein structuresand quantum mechanical calculations[J]．Proc Natl Acad Sci，2004，101（18）：6946-6951．

[17]Nilsson A，Ogasawara H，Cavalleri M，et al．The hydrogen bond in ice probed by soft x-ray spectroscopy and density functional theory[J]．J Chem Phys，2005，122（15）：154505．

[18]Car R，Parrinello M．Unified approach for molecular dynamicsand theory[J]．Phys Rev Lett，985，55（22）：2471-2474．

[19]Tuckerman ME，Parrinello M，Integratingthe Car-Parrinelloequations．I．Basic integrationtechniques[J]．JChem Phys，1994，101（2）：1302-1315．

[20]Grimme S．Density functional theory with London dispersion corrections[J]．WIREs Comput Mol Sci，2011，1：211-228．

[21]Zimmerli U，Parrinello M，Koumoutsakos P．Dispersioncorrectionstodensity functionalsfor water aromatic interactions[J]．JChem Phys，2004，120（6）：2693-2699．

[22]Schmidt J，VandeVondele V，Kuo IFW，et al．Isobaric-Isothermal Molecular Dynamics Simulations Utilizing Density Functional Theory:An Assessment of the Structureand Density of Water at Near-Ambient Conditions[J]．JPhys Chem B，2009，113（35）：11959-11964．