應用ABEEMσπ極化力場對Zn2+水溶液配位微結構和水交換反應進行分子動力學模擬研究

赫蘭蘭 郭 宇 趙 健 姜新蕊 楊忠志 趙東霞

(遼寧師范大學化學化工學院，遼寧大連116029)

應用ABEEMσπ極化力場對Zn2+水溶液配位微結構和水交換反應進行分子動力學模擬研究

赫蘭蘭 郭 宇 趙 健 姜新蕊 楊忠志*趙東霞*

(遼寧師范大學化學化工學院，遼寧大連116029)

應用ABEEMσπ極化力場，對Zn2+水溶液體系進行分子動力學模擬，探討Zn2+的配位微結構和配體水交換反應。水分子模型采用ABEEM-7P精細水模型。模擬后對體系結構、電荷及動力學性質進行細致分析。結構分析表明，平衡體系中Zn2+的第一層配位數為6，這與實驗值是一致的。水交換反應過程中，溶劑水由O―Zn―O角分線斜上(下)方進攻Zn2+，配位水由O―Zn―O角分線斜下(上)方逐漸遠離。極化力場模擬時Zn2+與交換水間的距離變化波動較大，而固定電荷力場的波動較小。模擬發現，極化力場的徑向分布函數能精細地展示第二、三層配體的配位微結構，第二配位層存在靠近Zn2+的亞殼層，能與第一水合層發生水交換反應，充分體現了Zn2+的極化效應。本文闡明了水交換反應中，Zn2+位點電荷與交換水中氧原子孤對電子位點電荷的規律性變化，從電荷的角度解釋了水交換反應的合理性。ABEEMσπ極化力場模擬Zn2+水溶液獲得第一水合層的平均配位駐留時間為2.0×10－9s，在實驗值范圍內，說明ABEEMσπ極化力場可以合理地模擬Zn2+水溶液體系。

ABEEMσπ極化力場；分子動力學模擬；水交換反應；徑向分布函數；電荷分析；平均配位駐留時間

1 引言

很多二價金屬離子在蛋白中都發揮著重要的作用1－4，例如鎂離子、鈣離子和鋅離子等。Zn2+涵蓋在六大類酶中，在300多種酶中充當協調因子，同時Zn2+在體內反應中作為Lewis酸，起到催化作用。然而，Zn2+參與的這些重要的生命過程都離不開水，少部分含鋅金屬酶的催化過程要利用Zn2+的水交換反應來完成5,6。例如碳酸酐酶II7調節人體內的酸堿平衡就涉及到水交換過程。所以深入理解Zn2+水溶液體系的結構和動態變化至關重要。

在大多數鋅的復合物中，鋅的配位形式主要有四面體構型的4配位形式和八面體構型的6配位形式8。在水溶液中，Zn2+與水分子形成八面體構型。關于Zn2+水溶液體系的性質，近些年來很多實驗組和理論組對其進行了相關研究5,9－16。實驗上的研究方法主要有XRD(X-ray diffration)，X射線吸收光譜和EXAFS(extended X-ray absorption final Structure)，而理論上主要采用QM/MM-MD(quantum mechanical/molecular mechanical-molecular dynamics)和分子力場方法。對于Zn2+水溶液第一水合層配位數和Zn―O平衡距離的問題，實驗和理論上都有明確的報道。例如Kuzmin等10通過EXAFS實驗測得Zn2+水溶液體系中鋅的第一水合層配位數為6.0±0.2，Zn―O平衡距離為(0.206±0.002) nm；Wu等11采用AMOEBA極化力場模擬出Zn2+第一水合層有6個水，同時在徑向分布函數(RDF)分析中Zn―O的第一峰值出現在0.198 nm處；Fatmi等15,17采用QM/MM-MD方法也模擬出Zn2+的第一水合層配位數是6，在RDF中Zn―O的平衡距離為0.211nm；Mohammed等16分別應用QM/MM-MD和經典動力學模擬(CMD)方法同樣也模擬出第一水合層有6個水，他們模擬出Zn―O的RDF第一峰值分別出現在0.216和0.223 nm。從上面研究可以發現Zn2+水溶液中第一水合層配位數為6，Zn―O平衡距離在0.198－0.223 nm之間。對于第二水合層水分子數，實驗和理論上沒有統一的標準。XRD實驗測出了Zn2+的第二層配位數范圍為8－13，并且第二層配體氧距Zn2+的距離為0.421－0.426 nm；Fatmi等15認為第二層配位數范圍為10－18，大多數情況是14，RDF顯示第二水合層的峰值為0.440 nm；Mohammed等16采用QM/MM-MD和CMD的方法模擬出第二層配體數分別為13和16，RDF的第二峰值分別出現在0.454和0.469 nm；Marini等18通過蒙特卡洛模擬得到第二水合層水分子數在20－22之間，水中氧原子距Zn2+的距離在0.430－0.440 nm之間。除此之外，對Zn2+水溶液的光譜信息、水分子的駐留時間和速率常數等也有報道。

然而，離子水溶液體系存在著一些極其有趣的動態變化，例如水交換反應，這一過程實驗上是很難捕獲到的。水交換反應存在著以下幾種機制，分別是：締合交換(A機制)、解離交換(D機制)、相互交換(I機制)、Ia機制和Id機制。該性質可以通過動力學模擬來實現。水交換過程必然會引起交換水分子的電荷重新分布，產生相應的極化效應，這一點是固定電荷力場不能描述的，可極化力場可以洞悉。目前，模擬水溶液體系的極化力場主要是AMOEBA極化力場11,12,19,20，該力場中考慮了二體相互作用和三體相互作用，雖然可以較精確地模擬離子水溶液體系，但是計算多級相互作用要消耗大量的時間，而且還不能應用到生物大分子體系的研究。Yang等21－25建立和發展的ABEEMσπ極化力場，將體系的電荷劃分到原子、鍵以及孤對電子位點上，可以真實地反映出分子中各位點電荷隨環境的改變，通過電荷分布的變化和轉移反映出極化效應的強弱。目前該力場不僅可以快速計算有機和生物大分子體系的總能量和電荷分布等性質，也成功地應用于水分子團簇26和純水溶液、離子水溶液體系27－30、多肽和蛋白質31,32以及核酸等體系。

本文應用ABEEMσπ極化力場對Zn2+水溶液的配位微結構和水交換反應進行研究，水溶液的水交換反應見式(1)，其中水分子采用ABEEM-7P模型，即將水分子劃分成七個位點，分別是三個原子位點，氧原子的兩個孤對電子位點以及兩個O―H間鍵位點。討論了極化力場模擬平衡后，Zn2+溶液體系的水交換反應過程中顯性水電荷的變化規律，Zn2+各配位層水分子的微結構，第一水合層水分子的平均配位駐留時間以及水交換反應的速率常數。

2 研究方法

2.1 ABEEMσπ極化力場

Zn2+水溶液體系的ABEEMσπ總能量，EABEEMσπ表達式為：

上式中Eb是鍵長伸縮振動能，Zn2+與配位水分子中氧原子之間采用Morse勢能函數，溶劑水分子的鍵長伸縮振動采用Morse勢能函數；公式(3)為該體系鍵長伸縮振動能的表達式。

Eθ是鍵角彎曲振動能，表達式如公式(4)

公式(3－4)中r是實際的鍵長，req是平衡的鍵長，θ是實際的鍵角，θeq是平衡的鍵角，kθ表示力場的鍵角力常數；D是鍵的解離能，α是和鍵長力常數相關的參數(α=(kb/2D)1/2，kb為鍵長伸縮振動力常數)。

EvdW和Eelec是非鍵相互作用勢能。其中vdW相互作用EvdW的表達式為：

公式(5)中，εab、σab和rab分別表示a與b間的勢阱深度、碰撞直徑和距離。此處采用了標準的聯合規則：εab=(εaaεbb)1/2，σab=1/2(σab+σab)，fab表示常數。在分子中，處在不同位置的a－b原子間的fab是不同的，1-2關系和1-3關系的fab=0.0，1-4關系的fab= 0.5，其他情況下fab=1.0。

ABEEMσπ模型的獨特之處在于靜電相互作用項Eelec

公式(6)中，qi和qj分別是位點i和j的電荷，rij是位點i與j間的距離。當i和j之間的最短連接路徑關系(包括σ鍵位點)小于1－6關系時，kij=0；當i和j在氫鍵相互作用區域時，kij=kH-bond22,23；其它所有情況，kij=0.57。

2.2 ABEEMσπ極化力場參數和ABEEM-7P水模型

ABEEMσπ極化力場中，分子的電荷位點包括原子位點、鍵(σ/π鍵)位點和孤對電子位點，各原子周圍鍵和孤對電子位點的電荷一般是各向異性的，不對稱的。各位點的電荷是應用電負性均衡方法，求解電負性均衡方程得到的，具體的能量公式和方程表達式請見文獻21－23。

ABEEMσπ極化分子力場擬合出了一套適用于Zn2+水溶液體系的參數。表1列出了Zn2+水溶液體系的電荷參數(價態電負性χ*和價態硬度2η*)、鍵長參數(平衡鍵長req，鍵的解離能D，以及與鍵長力常數相關的參數α)、鍵角參數(平衡鍵角θeq和鍵角力常數kθ)以及范德華參數(碰撞直徑σ和勢阱深度ε)。溶劑水的電荷、范德華、鍵長和鍵角參數仍采用原來ABEEMσπ力場的參數23,33，其余參數是最新擬合得到的。本文中考慮到第一層配體水分子與溶劑水分子所處的化學環境不同，受到鋅離子的靜電極化作用較強，為了更接近真實情況，第一層配體水分子和溶劑水分子采用不同的參數。當配體水分子和溶劑水分子發生交換時，水分子參數也發生相應的變化。

為了使這套參數不僅適用于Zn2+水溶液體系，而且也適用于我們將來要研究的含Zn2+蛋白體系。在確定電荷參數時，我們以B3LYP/6-311++G(d,p)水平下優化后的以及PDB數據庫中截取的三十多種含Zn2+體系作為模型分子，包含二、三、四、五、六個配位數，配體主要為氨基酸側鏈上的原子、水分子和抑制劑分子，Zn2+的基點電荷是1.0。采用從頭算的方法獲得模型分子中各原子的Mulliken電荷，并以此為標準，利用組內編寫的程序，根據線性回歸和最小二乘法擬合出了一套適用于Zn2+水溶液體系以及絕大多數含鋅蛋白電荷分布計算的參數。

Morse勢能參數的擬合，選取B3LYP/6-311++ G(d,p)水平下優化后的Zn2+-His-2Cys-H2O為模型分子，并利用從頭算的方法擬合Zn2+與水分子間的Morse勢能曲線，獲得Zn―O的Morse勢能參數；將參數加入ABEEMσπ力場中，優化1－6)的穩定結構和PDB數據庫中截取的Zn2+配體中含水分子的結構。反復調節參數，獲得優化后使體系Zn―O鍵長變動最小的參數。配體水的O―H鍵解離能仍采用溶劑水的O―H鍵解離能，平衡鍵長做了微調，獲得了適用于第一層配體水分子的平衡鍵長。

表1 Zn2+水溶液體系ABEEMσπ力場參數Table 1 Parameters of Zn2+aqueous solution forABEEMσπ force field

鍵角參數的擬合，選取Zn2+-His-2Cys-2(H2O)為模型分子，采用上述方法進行優化和頻率計算，獲得穩定的結構。分別對∠O―Zn―O和∠H―O―H每偏離平衡位置1°保存一個結構，共保存20個結構，擬合出∠O―Zn―O和∠H―O―H的勢能曲線，并獲得平衡鍵角和鍵角力常數。將參數加入ABEEMσπ力場中，以從頭算獲得的水合鋅離子以及PDB中截取的含Zn2+體系為標準，采用ABEEMσπ力場優化上述結構，對比優化前后的結構，反復調節參數，獲得優化后使∠O―Zn―O和∠H―O―H變動最小的參數。

Li等34報道的Zn2+的vdW參數表明，Zn2+的vdW參數受到金屬離子在力場中的模型、水分子模型和模擬環境等因素的影響，Zn2+的vdW參數可轉移性會受到局限。本文參考了Santosh等35在AMBER力場中進行分子動力學模擬時，研究二肽和Zn2+水溶液體系所用的Zn2+的vdW參數，碰撞直徑σ為0.20471 nm，勢阱深度ε為0.569 kJ·mol－1。以動力學模型平衡時Zn2+的配位結構與實驗數據相近為標準，在Santosh等人的Zn2+的vdW參數附近調節，獲得適用于ABEEMσπ極化力場的Zn2+的vdW參數。配體水中O和H的范德華參數是在溶劑水中O和H的范德華參數的基礎上微調得到的。

水分子七點極化力場非剛體模型(ABEEM-7P)22,23示于圖1。水分子的七點模型包含三個原子位點，兩個O―H化學鍵位點和兩個孤對電子位點。單體水中O―H鍵長和H―O―H鍵角分別設為相應的實驗值0.09572 nm和104.5°，O―H鍵位點設在氫原子和氧原子的共價半徑之比處，孤對電子的位點在距氧原子0.074 nm處。水分子中正電荷集中在氧原子和氫原子上，而負電荷分布在孤對電子區域和化學鍵區域。ABEEM-7P水模型不同于前人的TIP4P、POL5和其他的水分子模型，它具有兩個顯著的特點：(1)ABEEM-7P水模型提供了更大的適用范圍，它不僅可以計算原子、化學鍵和孤對電子的電荷，也提出了氫鍵相互作用區域的概念，通過引入氫鍵擬合函數kH-bond來描述所形成氫鍵的孤對電子和水中氫原子間的靜電相互作用，使氫鍵信息描述得更加精確；(2)在分子力學中，水分子為非剛體結構，允許鍵長和鍵角發生振動，鍵長伸縮振動使用Morse勢能函數來描述，因為它能夠從平衡鍵長到鍵的解離更大范圍描述鍵的伸縮振動。范德華相互作用采用12-6的Lennard-Jones勢能函數。ABEEM-7P模型可以很清晰地描述電荷隨環境的變化和氫鍵信息，所以使用該模型討論Zn2+水溶液的結構和性質。

圖1水分子的ABEEM-7P模型示意圖Fig.1 Diagram ofABEEM-7Pwater molecular model

2.3 平均配位駐留時間

水分子在離子周圍的平均配位駐留時間36可以通過時間相關函數來得到：

公式(7)中τ(r,t)是Heaviside單位函數。如果一個水分子i在半徑為r的區域范圍內，τ(r,t)等于1；如果一個水分子i不在半徑為r的區域范圍內，τ(r,t)等于0。本文中，計算Zn2+第一層配體水的平均配位駐留時間時，r=0.300 nm；第二層時，r在0.345－0.560 nm范圍內；第三層時，r在0.560－0.800 nm范圍內。Nr代表t=0時被統計區域內水分子的個數，同時允許水分子短暫離開區域r的最長時間不超過2 ps。如果水分子在離子周圍的停留時間較長，2 ps的時間不會影響統計計算結果。

2.4 過渡態理論速率常數

過渡態理論的速率常數20表達式為：

其中300 K時校正因子k0≈(kBT/h)≈0.6×1013s－1，R=8.314×10－3kJ·mol－1·K－1，ΔG≠是發生水交換反應的活化能，通過分析平均力勢(PMF)36獲得，平均力勢及有效平均力勢(EPMF)可以通過徑向分布函數獲得，如公式(9－10)：

其中β=1/(kBT)，kB為玻爾茲曼因子，T為熱力學溫度，r代表離子與水分子中氧原子間的距離，r*代表活化勢壘的位置，W(r)和Weff(r)分別表示平均力勢和有效平均力勢。

2.5 模擬細節

將Zn2+融入Tinker程序中具有500個水分子的waterbox立方體水盒子中，模擬體系中水分子的個數為497個，加入兩個抗衡離子Cl－，保持體系凈電荷為0。體系盒長為2.4662 nm。ABEEMσπ力場中各區域劃分情況見圖2，包括原子、鍵(Zn―O/O―H鍵)和孤對電子區域。分別采用ABEEMσπ固定電荷和極化力場的成鍵模型，對建好的模擬體系進行優化和動力學模擬。采用正則系綜(NVT)，利用Berendsen熱浴控制每個分子的溫度，并利用Verlet跳蛙法積分運動方程，選取積分步長為1 fs，極化力場每隔1 ps重新計算體系中各位點的電荷，各體系每隔1 ps保存一次坐標文件，模擬時間為2 ns?？紤]周期性邊界條件和最近鏡像，采用截斷半徑方法(截斷半徑均為1.2331 nm)，并且利用開關函數調節非鍵相互作用。在整個模擬過程中，對Zn2+和抗衡離子進行固定，其余原子可以自由移動。模擬初速度通過Maxwell-Boltzmann分布隨機給出，模擬過程中原子的運動遵守牛頓運動方程，通過調節速度使體系達到設定的溫度300 K，整個體系保持局域守恒。

圖2 ABEEMσπ力場中Zn(H2O)62+各區域示意圖Fig.2 Diagram of Zn(H2O)62+divided into regions in ABEEMσπ force field

3 結果與討論

3.1 結構性質及電荷分析

3.1.1 捕獲水交換反應

圖3 ABEEMσπ力場模擬Zn2+水溶液體系過程直接觀測到水交換反應Fig.3 Direct observation of the water exchange events near Zn2+in water based onABEEMσπ force field

3.1.2 ABEEMσπ極化力場水交換反應過程中電荷的變化

在水交換反應的過程中，交換水中氧原子孤對電子位點電荷也發生了規律性的變化，下面具體討論了水溶液體系水交換反應過程中電荷的變化。

表2 ABEEMσπ極化力場動力學模擬時間在1－15 ps內，Zn(H2O)62+水溶液體系水中氧原子與Zn2+間的距離變化情況Table 2 Distance between Zn2+and oxygen atom of water,during dynamic simulation time 1－15 ps, based onABEEMσπ polarizable force field

圖4 ABEEMσπ極化力場Zn2+第一層配體水交換的結構示意圖Fig.4 Structure diagrams of the first shell ligand exchange events near Zn2+in water based on ABEEMσπ polarizable force field

表3 ABEEMσπ極化力場動力學模擬時間為1－15 ps內，Zn(H2O)26+水溶液體系中Zn2+位點電荷及交換水中氧原子的孤對電子位點電荷Table 3 Charge of Zn2+site and lone pair site of oxygen atom of exchange water,during dynamic simulation time 1－15 ps,based onABEEMσπ polarizable force field

模擬時間為2－3 ps的過程中發生了第二次水交換反應(圖4(c)是模擬時間為3 ps時的結構)。O421處于Zn2+的一、二配位層之間，從O496―Zn―O1177角分線的斜下方進攻Zn2+，到達了Zn2+的第一配位層，占據O1177的位置，O421孤對電子位點電荷獲得了－0.059e。同時O1177從右下方逐漸遠離Zn2+，到達了Zn2+的第二配位層，距Zn2+0.411 nm，O1177孤對電子位點電荷失去了－0.078e。

模擬時間為5－6 ps過程中發生了第三次水交換反應。圖4(d)和(e)分別代表模擬時間為5和6 ps時的結構。5 ps時，O496、O343和O421與Zn2+距離分別為0.270、0.245、0.261 nm，此時形成七配位結構，Zn2+的雜化方式為sp3d3雜化，但這個結構并沒有穩定存在。模擬時間為6 ps時，O496與Zn2+間的距離為0.446 nm；而O343和O421仍為Zn2+的第一層配體，與Zn2+距離分別為0.213、0.212 nm；此時Zn2+的雜化方式為sp3d2。在這個水交換反應過程中交換水的孤對電子位點電荷發生了規律性的變化。模擬時間為4－6 ps，O496由Zn2+的第一配位層運動到第二配位層，這個過程中O496的孤對電子位點電荷由－0.515e變為－0.455e。5－6 ps時，O421和O343孤對電子位點電荷分別增加了－0.036e和－0.049e，同時Zn2+上的電荷也由1.036e變為1.066e，增加了0.03e。這是由于Zn2+與交換水距離變短，它們之間的靜電吸引作用加強。

模擬時間為11－12 ps時，發生了第四次水交換反應。仍可得到和上述一致的結論，從電荷的角度分析，處于Zn2+第二配位層之內的水分子受到Zn2+的靜電極化作用較強；靠近Zn2+時，交換水中O的孤對電子位點負電荷電量增大；遠離Zn2+時，交換水中O的孤對電子位點負電量減少，這樣更容易發生水交換反應。

ABEEMσπ極化力場模擬Zn2+水溶液體系時，在水交換反應過程中Zn2+與交換水孤對電子位點電荷發生了規律性的變化，從電荷的角度解釋了水交換反應的發生，這是固定電荷力場和AMOEBA極化力場不能實現的。

3.1.3 徑向分布函數

為了描述ABEEMσπ力場中Zn2+周圍水分子的微觀結構，我們分析了動力學模擬平衡后Zn―O之間的徑向分布函數(RDF)。圖5中黑線和紅線分別表示ABEEMσπ固定電荷和極化力場模擬Zn2+水溶液體系的Zn―O徑向分布函數圖。表437中列出了ABEEMσπ固定電荷及極化力場模擬Zn2+水溶液體系平衡后，Zn2+各配位層的配位數及RDF的峰值位置。從圖5中可以看出固定電荷力場與極化力場獲得的RDF有很大區別。極化力場RDF的第一峰要比固定電荷力場的高出一倍，說明極化力場Zn2+與第一層配體之間的相互作用強于固定電荷力場的。固定電荷和極化力場獲得的RDF第一峰的積分數都是6，說明Zn2+在水溶液中以六配位的形式存在較為穩定，這與實驗值10是一致的。ABEEMσπ極化力場的RDF的第二層配體的分布情況不同于ABEEMσπ固定電荷力場和AMOEBA極化力場11；ABEEMσπ固定電荷力場以及AMOEBA極化力場的RDF第二峰是一個較為平緩的峰，而ABEEMσπ極化力場第二層配體存在幾個亞殼層，并且是很陡的尖峰，充分展現了距Zn2+一定距離的很小殼層內水分子的分布有一定秩序；這種微結構體現了Zn2+對水溶液體系的極化效應。ABEEMσπ固定電荷力場第二層配體氧與Zn2+距離的范圍為0.320－0.510 nm，第二層配位數約為18.2個；極化力場第二層配體與Zn2+間距離的范圍為0.348－0.560 nm，共28個水分子，與Zn2+距離范圍為0.345－0.380 nm的亞殼層內有7個水分子，這些水分子可與第一層配體發生水交換反應。相比于ABEEMσπ固定電荷力場，ABEEMσπ極化力場配體與Zn2+距離0.628 nm處出現了第三峰，第三層配體的配位數約為50.0個。以上分析表明ABEEMσπ極化力場模擬Zn2+溶液與固定電荷力場以及AMOEBA極化力場相比更精確。

圖5 ABEEMσπ力場Zn2+水溶液體系Zn―O徑向分布函數(RDF)以及積分線Fig.5 RDFs of Zn―O and integration lines for Zn2+aqueous solution based onABEEMσπ force field color online

3.1.4 角度分布函數

圖6是ABEEMσπ固定電荷力場和極化力場模擬Zn2+水溶液體系平衡后O―Zn―O的角度分布函數圖。從圖中可以看出固定電荷力場和極化力場模擬Zn2+水溶液體系平衡后，O―Zn―O的角度分布函數圖中有兩個明顯的峰，固定電荷力場得到的兩個峰值分別出現在88.5°和171.5°，峰的范圍較寬且存在小的肩峰，是較為扭曲的八面體構型。而極化力場得到的兩個峰分別出現在89.5°和178.5°，并且峰的寬度較小，是比較規則的八面體構型。另外，極化力場獲得的O―Zn―O角度分布函數的峰值遠大于固定電荷力場的峰值。說明極化力場所模擬的Zn2+水溶液體系形成6配位的八面體構型的幾率更大。

表4 ABEEMσπ固定電荷力場和極化力場模擬Zn2+水溶液體系平衡后，Zn2+各配位層的配位數及Zn―O的RDF出峰位置與AMOEBA、QM/MM和EXAFS對比Table 4 Coordination number of different hydration shells and the peaks of Zn―O RDFs in Zn2+aqueous solution equilibrium system based onABEEMσπ fixed force field and polarizable force field,AMOEBA,QM/MM and EXAFS

圖6 ABEEMσπ力場Zn2+水溶液體系O―Zn―O角度分布函數圖Fig.6 Angle distribution function of O―Zn―O in Zn2+aqueous solution based onABEEMσπ force field

3.2 動力學性質分析

3.2.1 平均配位駐留時間的統計

通過平均配位駐留時間的相關函數可獲得Zn2+第一水合層水分子的駐留時間(見公式(7))。經統計ABEEMσπ可極化力場中Zn2+水溶液的第一層配體水分子的平均配位駐留時間為2.0×10－9s，實驗值38在10－10－10－9s之間。AMOEBA極化力場11Zn2+溶液第一層配體水的平均駐留時間為2.2× 10－9s，Classical固定電荷力場11分子動力學模擬得到第一層配體的平均駐留時間約為1.5×10－10s。ABEEMσπ可極化力場獲得的Zn2+第一水合層的平均配位駐留時間與其他力場相差不大，且在實驗值范圍內。ABEEMσπ可極化力場模擬Zn2+水溶液體系得到第二、三層配體水分子的平均配位駐留時間分別為4.3×10－10和4.0×10－10s。從上述數據可以看出，ABEEMσπ極化力場外層水的平均配位駐留時間要比內層水的短，說明外層水分子受到Zn2+的極化作用小。

3.2.2 水交換反應速率常數的統計

圖7 ABEEMσπ極化力場Zn2+與水分子的平均力勢(W(r))和有效平均力勢(Weff(r))隨Zn2+和水分子間的距離變化情況Fig.7 Potential of mean force(W(r))and effective potential of mean force(Weff(r))between Zn2+and water with the distance change between Zn2+and water based on ABEEMσπ polarizable force field

圖7展示了Zn2+與水分子間的平均力勢和有效平均力勢隨Zn2+與水分子間的距離變化情況。從圖中可以獲悉ABEEMσπ極化力場模擬Zn2+水溶液體系，獲得Zn2+與第一水合層水分子的解離能壘約為7.8kBT(1kBT約為2.5 kJ·mol－1)。與其他固定電荷力場相比，ABEEMσπ極化力場的優勢是不僅可以獲得第一層配體水交換反應的活化能壘，而且可以獲得第二、三層配體水交換的活化能壘，它們分別約為0.5kBT和1.0kBT。ABEEMσπ極化力場的PMF活化能壘中，第一水合層的水分子與Zn2+之間的解離能壘最大，這是由于第一水合層的水分子與Zn2+間的相互作用較強；第二層配體水分子發生水交換反應的能壘最小，說明第二層配體很容易與第一層配體和第三層配體發生水交換反應，自由移動性強；第三層配體水交換反應能壘比第二層配體的大，這是由于第三層配體間水分子的平均氫鍵數目比第二層配體的多，水分子移動時要克服的氫鍵相互作用較大。ABEEMσπ極化力場可以精細地反映各層配體水發生水交換反應的難易情況，與其特有的極化力場和ABEEM-7P水模型密切相關。

由平均力勢分析獲得ABEEMσπ極化力場第一水合層水交換反應的活化能為19.5 kJ·mol－1，Michael等39采用從頭算的方法獲得氣相下Zn2+第一水合層水交換反應的活化能壘在17.6－19.2 kJ·mol－1之間，ABEEMσπ極化電荷力場獲得的第一水合層水交換反應的活化能與從頭算結果非常接近，說明ABEEMσπ極化電荷力場模擬Zn2+水溶液體系的結果是合理的。通過公式(8)獲得極化力場過渡態理論的速率常數為2.4×109s－1；其倒數與我們獲得的第一層配體水分子的平均配位駐留時間吻合。實驗報道40，298 K時Zn2+水溶液中第一層配體水交換反應的速率常數約為2×107s－1。溫度升高水交換反應的速率會加快，這說明ABEEMσπ力場獲得300 K時，Zn2+第一水合層水交換反應的速率常數的數量級是合理的。

4 結論

應用ABEEMσπ固定電荷力場和極化力場模擬Zn2+水溶液體系，獲得的Zn2+第一水合層的配位數與實驗值一致。對RDF進行分析，獨特的ABEEMσπ極化力場和ABEEM-7P精細水模型模擬Zn2+水溶液發現了第一、二和三層配體水分子的微結構，Zn2+的第二層配體存在亞殼層，可與第一層配體發生水交換反應。ABEEMσπ極化力場模擬Zn2+溶液體系獲得了水交換反應過程中Zn2+位點與交換水分子孤對電子位點電荷的規律性變化，解釋了水交換反應的合理性。ABEEMσπ極化力場Zn2+水溶液體系第一水合層水分子平均配位駐留時間在實驗值范圍內，并且第一水合層水分子的解離能壘與從頭算很接近，說明了ABEEMσπ極化力場模擬Zn2+水溶液體系的結果是合理的。本文對Zn2+水溶液體系的配位微結構以及水交換反應的探討，為進一步研究其它過渡金屬離子水溶液中的水交換反應，以及一些酶催化過程中的水交換反應奠定基礎。

致謝：我們對Ponder教授提供的Tinker源程序表示衷心的感謝。

(1) Cowan,J.A.Chem.Rev.1998,98(3),1067.doi:10.1021/ cr960436q

(2) Silverman,D.N.;McKenna,R.Accounts Chem.Res.2007,40 (8),669.doi:10.1021/ar7000588

(3) Imanishi,M.;Matsumura,K.;Tsuji,S.;Nakaya,T.;Negi,S.; Futaki,S.;Sugiura,Y.Biochemistry 2012,51(16),3342. doi:10.1021/bi300236m

(4) Kepp,K.P.Chem.Rev.2012,112(10),5193.doi:10.1021/ cr300009x

(5) Ohtaki,H.;Radnai,T.Chem.Rev.1993,93(3),1157. doi:10.1021/cr00019a014

(6) Helm,L.;Merbach,A.E.Coord.Chem.Rev.1999,187(1), 151.doi:10.1016/s0010-8545(99)90232-1

(7) Nielsen,P.M.;Fago,A.J.Inorg.Biochem.2015,149,6. doi:10.1016/j.jinorgbio.2015.04.013

(8) Maynard,A.T.;Covell,D.G.J.Am.Chem.Soc.2001,123(6), 1047.doi:10.1021/ja0011616

(9) Asthagiri,D.;Pratt,L.R.;Paulaitis,M.E.;Rempe,S.B.J.Am. Chem.Soc.2004,126,1285.doi:10.1021/ja0382967

(10) Kuzmin,A.;Obst,S.;Purans,J.J.Phys.:Condens.Matter 1997,9,10065.

(11) Wu,J.C.;Piquemal,J.P.;Chaudret,R.;Reinhardt,P.;Ren,P.Y. J.Chem.Theory Comput.2010,6(7),2059.doi:10.1021/ ct100091j

(12) Xiang,J.Y.;Ponder,J.W.J.Comput.Chem.2013,34(9),739. doi:10.1002/jcc.23190

(13) Babu,C.S.;Lim,C.J.Phys.Chem.A 2006,110(2),691. doi:10.1021/jp054177x

(14) Fatmi,M.Q.;Hofer,T.S.;Randolf,B.R.;Rode,B.M.J.Phys. Chem.B 2008,112(18),5788.doi:10.1021/jp710270z

(15) Fatmi,M.Q.;Hofer,T.S.;Randolf,B.R.;Rode,B.M.J.Phys. Chem.B 2006,110(1),616.doi:10.1021/jp0546655

(16) Mohammed,A.M.;Loeffler,H.H.;Inada,Y.;Tanada,K.; Funahashi,S.J.Mol.Liq.2005,119(1－3),55.doi:10.1016/j. molliq.2004.10.008

(17) Fatmi,M.Q.;Hofer,T.S.;Randolf,B.R.;Rode,B.M.J.Phys. Chem.B 2007,111(1),151.doi:10.1021/jp0654213

(18) Marini,G.W.;Texler,N.R.;Rode,B.M.J.Phys.Chem.1996, 100(16),6808.doi:10.1021/jp953375t

(19) Soniat,M.;Hartman,L.;Rick,S.W.J.Chem.Theory Comput. 2015,11(4),1658.doi:10.1021/ct501173n

(20) Kurnikov,I.V.;Kurnikova,M.J.Phys.Chem.B 2015,119(32), 10275.doi:10.1021/acs.jpcb.5b01295

(21) Zhao,D.X.;Liu,C.;Wang,F.F.;Yu,C.Y.;Gong,L.D.;Liu,S. B.;Yang,Z.Z.J.Chem.Theory Comput.2010,6(3),795. doi:10.1021/ct9006647

(22) Yang,Z.Z.;Wu,Y.;Zhao,D.X.J.Chem.Phys.2004,120(6), 2541.doi:10.1063/1.1640345

(23) Wu,Y.;Yang,Z.Z.J.Phys.Chem.A 2004,108(37),7563. doi:10.1021/jp0493881

(24) Li,X.;Yang,Z.Z.J.Phys.Chem.A 2005,109(18),4102. doi:10.1021/jp0458093

(25) Yang,Z.Z.;Cui,B.Q.J.Chem.Theory Comput.2007,3(4), 1561.doi:10.1021/ct600379n

(26) Qian,P.;Yang,Z.Z.Acta Phys.-Chim.Sin.2006,22(5),561. [錢 萍,楊忠志.物理化學學報,2006,22(5),561.] doi:10.3866/PKU.WHXB20060510

(27) Gong,L.D.Sci.Sin.Chim.2013,43(5),519.[宮利東.中國科學:化學,2013,43(5),519.]doi:10.1007/s11426-012-4787-3

(28) Li,X.;Yang,Z.Z.J.Chem.Phys.2005,122(8),084514. doi:10.1063/1.1853372

(29) Guo,C.;Liu,C.;Yang,Z.Z.Acta Phys.-Chim.Sin.2010,26(2),478.[郭 慈,劉 翠,楊忠志.物理化學學報,2010,26 (2),478.]doi:10.3866/PKU.WHXB20100219

(30) Liu,Y.;Wang,F.F.;Yu,C.Y.;Liu,C.;Gong,L.D.;Yang,Z.Z. Acta Phys.-Chim.Sin.2011,27(2),379.[劉 燕,王芳芳,于春陽,劉 翠,宮利東,楊忠志.物理化學學報,2011,27(2), 379.]doi:10.3866/PKU.WHXB20110233

(31) Zhang,Q.;Yang,Z.Z.Chem.J.Chin.Univ.2005,26(12), 2345.[張 強,楊忠志.高等學?；瘜W學報,2005,26(12), 2345.]doi:10.3321/j.issn:0251-0790.2005.12.029

(32) Guan,Q.M.;Cui,B.Q.;Zhao,D.X.;Gong,L.D.;Yang,Z.Z. Chin.Sci.Bull.2008,53(8),1171.doi:10.1007/s11434-007-0493-5

(33) Yang,Z.Z.;Wu,Y.;Zhao,D.X.J.Chem.Phys.2004,120(6), 2541.doi:10.1063/1.1640345

(34) Li,P.F.;Roberts,B.P.;Chakravorty,D.K.;Merz,K.M.,Jr.J. Chem.Theory Comput.2013,9(6),2733.doi:10.1021/ ct400751u

(35) Santosh,M.S.;Lyubartsev,A.P.;Mirzoev,A.A.;Bhat,D.K.J. Phys.Chem.B 2010,114(49),16632.doi:10.1021/jp108376j

(36) Sp?ngberg,D.;Rey,R.;Hynes,J.T.;Hermansson,K.J.Phys. Chem.B 2003,107(18),4470.doi:10.1021/jp027230f

(37) Zhang,Q.;Yang,Z.Z.Chem.Phys.Lett.2005,403,242. doi:10.1016/j.cplett.2005.01.011

(38) Salmon,P.S.;Bellissent-Funel,M.C.;Herdman,G.J.J.Phys.: Condens.Matter 1990,2,4297.doi:10.1088/0953-8984/2/18/ 027

(39) Michael,H.;Clark,T.;Eldik,R.V.J.Am.Chem.Soc.1997,119 (33),7843.doi:10.1021/ja970483f

(40) Akesson,R.;Pettersson,L.G.M.;Sandstroem,M.;Siegbahn,P. E.M.;Wahlgren,U.J.Phys.Chem.1993,97(15),3765. doi:10.1021/j100117a023

Study on Coordination Microstructure and Water Exchange Reaction of Zn2+Aqueous Solutions through Molecular Dynamics Simulations Using the ABEEMσπ Polarizable Force Field

HE Lan-Lan GUO Yu ZHAO Jian JIANG Xin-Rui YANG Zhong-Zhi*ZHAO Dong-Xia*
(School of Chemistry and Chemical Engineering,Liaoning Normal University,Dalian 116029,Liaoning Province,P.R.China)

Coordination microstructure and water exchange of Zn2+aqueous solution were studied through molecular dynamics simulations based on the ABEEMσπ polarizable force field with the precise ABEEM-7P model.Structural and dynamical properties were investigated in detail.We show that the first-shell water coordination number is six,which is in agreement with the experimental value.During the water exchange process,H2O,which attacks or leaves the first solvation shell of Zn2+,is orientated on the upper or lower inclined side of the∠O―Zn―O bisector.The distance change between Zn2+and the oxygen atom of the exchange water for the polarizable force field fluctuates more than that observed for the fixed charge force field.The radial distribution function(RDF)of the polarizable force field showed clearly the microstructure of the second and third solvation shells.The subshell of the second solvation shell was found to exchange with the first solvation shell.The polarizable effect of Zn2+has been fully expressed.The charges of the Zn2+site and the lone-pair sitesin exchange water regularly change,demonstrating the rationality of water exchange.The mean ligand residence time(2.0×10－9s)of the first-shell water produced by theABEEMσπ polarizable force field is within the range of the experimental value.Zn2+aqueous solution can be reasonably simulated through theABEEMσπ polarizable force field.

ABEEMσπ polarizable force field;Molecular dynamic simulation;Water exchange reaction; Radial distribution function;Charge analysis;Mean ligand residence time

O641

10.3866/PKU.WHXB201609193

Received:August 9,2016;Revised:September 19,2016;Published online:September 19,2016.

*Corresponding authors.YANG Zhong-Zhi,Email:zzyang@lnnu.edu.cn;Tel:+86-411-82159607.

ZHAO Dong-Xia,Email:zhaodxchem@lnnu.edu.cn.

The project was supported by the National Natural Science Foundation of China(21133005,21473083),General Project of Education Department of Liaoning Province,China(L2014426),and Laboratory Opening Program of Liaoning Normal University,China(cx20160111).

國家自然科學基金(21133005,21473083),遼寧省教育廳一般項目(L2014426)和遼寧師范大學實驗室開放項目(cx20160111)資助