應(yīng)用ABEEMσπ極化力場(chǎng)對(duì)Zn2+水溶液配位微結(jié)構(gòu)和水交換反應(yīng)進(jìn)行分子動(dòng)力學(xué)模擬研究

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

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

應(yīng)用ABEEMσπ極化力場(chǎng)對(duì)Zn2+水溶液配位微結(jié)構(gòu)和水交換反應(yīng)進(jìn)行分子動(dòng)力學(xué)模擬研究

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

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

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

ABEEMσπ極化力場(chǎng)；分子動(dòng)力學(xué)模擬；水交換反應(yīng)；徑向分布函數(shù)；電荷分析；平均配位駐留時(shí)間

1 引言

很多二價(jià)金屬離子在蛋白中都發(fā)揮著重要的作用1－4，例如鎂離子、鈣離子和鋅離子等。Zn2+涵蓋在六大類酶中，在300多種酶中充當(dāng)協(xié)調(diào)因子，同時(shí)Zn2+在體內(nèi)反應(yīng)中作為L(zhǎng)ewis酸，起到催化作用。然而，Zn2+參與的這些重要的生命過程都離不開水，少部分含鋅金屬酶的催化過程要利用Zn2+的水交換反應(yīng)來完成5,6。例如碳酸酐酶II7調(diào)節(jié)人體內(nèi)的酸堿平衡就涉及到水交換過程。所以深入理解Zn2+水溶液體系的結(jié)構(gòu)和動(dòng)態(tài)變化至關(guān)重要。

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

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

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

2 研究方法

2.1 ABEEMσπ極化力場(chǎng)

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

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

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

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

EvdW和Eelec是非鍵相互作用勢(shì)能。其中vdW相互作用EvdW的表達(dá)式為：

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

ABEEMσπ模型的獨(dú)特之處在于靜電相互作用項(xiàng)Eelec

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

2.2 ABEEMσπ極化力場(chǎng)參數(shù)和ABEEM-7P水模型

ABEEMσπ極化力場(chǎng)中，分子的電荷位點(diǎn)包括原子位點(diǎn)、鍵(σ/π鍵)位點(diǎn)和孤對(duì)電子位點(diǎn)，各原子周圍鍵和孤對(duì)電子位點(diǎn)的電荷一般是各向異性的，不對(duì)稱的。各位點(diǎn)的電荷是應(yīng)用電負(fù)性均衡方法，求解電負(fù)性均衡方程得到的，具體的能量公式和方程表達(dá)式請(qǐng)見文獻(xiàn)21－23。

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

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

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

表1 Zn2+水溶液體系A(chǔ)BEEMσπ力場(chǎng)參數(shù)Table 1 Parameters of Zn2+aqueous solution forABEEMσπ force field

鍵角參數(shù)的擬合，選取Zn2+-His-2Cys-2(H2O)為模型分子，采用上述方法進(jìn)行優(yōu)化和頻率計(jì)算，獲得穩(wěn)定的結(jié)構(gòu)。分別對(duì)∠O―Zn―O和∠H―O―H每偏離平衡位置1°保存一個(gè)結(jié)構(gòu)，共保存20個(gè)結(jié)構(gòu)，擬合出∠O―Zn―O和∠H―O―H的勢(shì)能曲線，并獲得平衡鍵角和鍵角力常數(shù)。將參數(shù)加入ABEEMσπ力場(chǎng)中，以從頭算獲得的水合鋅離子以及PDB中截取的含Zn2+體系為標(biāo)準(zhǔn)，采用ABEEMσπ力場(chǎng)優(yōu)化上述結(jié)構(gòu)，對(duì)比優(yōu)化前后的結(jié)構(gòu)，反復(fù)調(diào)節(jié)參數(shù)，獲得優(yōu)化后使∠O―Zn―O和∠H―O―H變動(dòng)最小的參數(shù)。

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

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

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

2.3 平均配位駐留時(shí)間

水分子在離子周圍的平均配位駐留時(shí)間36可以通過時(shí)間相關(guān)函數(shù)來得到：

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

2.4 過渡態(tài)理論速率常數(shù)

過渡態(tài)理論的速率常數(shù)20表達(dá)式為：

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

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

2.5 模擬細(xì)節(jié)

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

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

3 結(jié)果與討論

3.1 結(jié)構(gòu)性質(zhì)及電荷分析

3.1.1 捕獲水交換反應(yīng)

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

3.1.2 ABEEMσπ極化力場(chǎng)水交換反應(yīng)過程中電荷的變化

在水交換反應(yīng)的過程中，交換水中氧原子孤對(duì)電子位點(diǎn)電荷也發(fā)生了規(guī)律性的變化，下面具體討論了水溶液體系水交換反應(yīng)過程中電荷的變化。

表2 ABEEMσπ極化力場(chǎng)動(dòng)力學(xué)模擬時(shí)間在1－15 ps內(nèi)，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σπ極化力場(chǎng)Zn2+第一層配體水交換的結(jié)構(gòu)示意圖Fig.4 Structure diagrams of the first shell ligand exchange events near Zn2+in water based on ABEEMσπ polarizable force field

表3 ABEEMσπ極化力場(chǎng)動(dòng)力學(xué)模擬時(shí)間為1－15 ps內(nèi)，Zn(H2O)26+水溶液體系中Zn2+位點(diǎn)電荷及交換水中氧原子的孤對(duì)電子位點(diǎn)電荷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

模擬時(shí)間為2－3 ps的過程中發(fā)生了第二次水交換反應(yīng)(圖4(c)是模擬時(shí)間為3 ps時(shí)的結(jié)構(gòu))。O421處于Zn2+的一、二配位層之間，從O496―Zn―O1177角分線的斜下方進(jìn)攻Zn2+，到達(dá)了Zn2+的第一配位層，占據(jù)O1177的位置，O421孤對(duì)電子位點(diǎn)電荷獲得了－0.059e。同時(shí)O1177從右下方逐漸遠(yuǎn)離Zn2+，到達(dá)了Zn2+的第二配位層，距Zn2+0.411 nm，O1177孤對(duì)電子位點(diǎn)電荷失去了－0.078e。

模擬時(shí)間為5－6 ps過程中發(fā)生了第三次水交換反應(yīng)。圖4(d)和(e)分別代表模擬時(shí)間為5和6 ps時(shí)的結(jié)構(gòu)。5 ps時(shí)，O496、O343和O421與Zn2+距離分別為0.270、0.245、0.261 nm，此時(shí)形成七配位結(jié)構(gòu)，Zn2+的雜化方式為sp3d3雜化，但這個(gè)結(jié)構(gòu)并沒有穩(wěn)定存在。模擬時(shí)間為6 ps時(shí)，O496與Zn2+間的距離為0.446 nm；而O343和O421仍為Zn2+的第一層配體，與Zn2+距離分別為0.213、0.212 nm；此時(shí)Zn2+的雜化方式為sp3d2。在這個(gè)水交換反應(yīng)過程中交換水的孤對(duì)電子位點(diǎn)電荷發(fā)生了規(guī)律性的變化。模擬時(shí)間為4－6 ps，O496由Zn2+的第一配位層運(yùn)動(dòng)到第二配位層，這個(gè)過程中O496的孤對(duì)電子位點(diǎn)電荷由－0.515e變?yōu)椋?.455e。5－6 ps時(shí)，O421和O343孤對(duì)電子位點(diǎn)電荷分別增加了－0.036e和－0.049e，同時(shí)Zn2+上的電荷也由1.036e變?yōu)?.066e，增加了0.03e。這是由于Zn2+與交換水距離變短，它們之間的靜電吸引作用加強(qiáng)。

模擬時(shí)間為11－12 ps時(shí)，發(fā)生了第四次水交換反應(yīng)。仍可得到和上述一致的結(jié)論，從電荷的角度分析，處于Zn2+第二配位層之內(nèi)的水分子受到Zn2+的靜電極化作用較強(qiáng)；靠近Zn2+時(shí)，交換水中O的孤對(duì)電子位點(diǎn)負(fù)電荷電量增大；遠(yuǎn)離Zn2+時(shí)，交換水中O的孤對(duì)電子位點(diǎn)負(fù)電量減少，這樣更容易發(fā)生水交換反應(yīng)。

ABEEMσπ極化力場(chǎng)模擬Zn2+水溶液體系時(shí)，在水交換反應(yīng)過程中Zn2+與交換水孤對(duì)電子位點(diǎn)電荷發(fā)生了規(guī)律性的變化，從電荷的角度解釋了水交換反應(yīng)的發(fā)生，這是固定電荷力場(chǎng)和AMOEBA極化力場(chǎng)不能實(shí)現(xiàn)的。

3.1.3 徑向分布函數(shù)

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

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

3.1.4 角度分布函數(shù)

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

表4 ABEEMσπ固定電荷力場(chǎng)和極化力場(chǎng)模擬Zn2+水溶液體系平衡后，Zn2+各配位層的配位數(shù)及Zn―O的RDF出峰位置與AMOEBA、QM/MM和EXAFS對(duì)比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σπ力場(chǎng)Zn2+水溶液體系O―Zn―O角度分布函數(shù)圖Fig.6 Angle distribution function of O―Zn―O in Zn2+aqueous solution based onABEEMσπ force field

3.2 動(dòng)力學(xué)性質(zhì)分析

3.2.1 平均配位駐留時(shí)間的統(tǒng)計(jì)

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

3.2.2 水交換反應(yīng)速率常數(shù)的統(tǒng)計(jì)

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

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

4 結(jié)論

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

致謝：我們對(duì)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. [錢 萍,楊忠志.物理化學(xué)學(xué)報(bào),2006,22(5),561.] doi:10.3866/PKU.WHXB20060510

(27) Gong,L.D.Sci.Sin.Chim.2013,43(5),519.[宮利東.中國(guó)科學(xué):化學(xué),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.[郭 慈,劉 翠,楊忠志.物理化學(xué)學(xué)報(bào),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.[劉 燕,王芳芳,于春陽,劉 翠,宮利東,楊忠志.物理化學(xué)學(xué)報(bào),2011,27(2), 379.]doi:10.3866/PKU.WHXB20110233

(31) Zhang,Q.;Yang,Z.Z.Chem.J.Chin.Univ.2005,26(12), 2345.[張 強(qiáng),楊忠志.高等學(xué)?；瘜W(xué)學(xué)報(bào),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).

國(guó)家自然科學(xué)基金(21133005,21473083),遼寧省教育廳一般項(xiàng)目(L2014426)和遼寧師范大學(xué)實(shí)驗(yàn)室開放項(xiàng)目(cx20160111)資助