李曉蕾，孫云嬌，唐 穎，王長生

（遼寧師范大學(xué)化學(xué)化工學(xué)院,大連116029）

生物分子與水分子間的氫鍵協(xié)同效應(yīng)對水環(huán)境中蛋白及DNA的分子結(jié)構(gòu)和性質(zhì)有重要影響.多體極化作用是氫鍵協(xié)同性的主要來源，快速準確計算多體極化作用強度對提升分子模擬方法的準確性十分重要［1~3］.

Xantheas等［4~6］發(fā)現(xiàn)，三體作用對水團簇穩(wěn)定性的影響可高達30%，四體及以上多體作用的貢獻可被忽略.對甲烷水復(fù)合物CH4（H2O）2和CH4（H2O）20的研究表明，三體作用對結(jié)合能的貢獻不可忽略［7］.Paesani等［8，9］通過評估團簇H3O+（H2O）n和OH?（H2O）n的多體作用強度發(fā)現(xiàn)，體系的總相互作用能主要由二體和三體作用貢獻.楊微等［10］使用高精度從頭計算方法計算了團簇（H2O）n（n=8，10，16，20，22，24）的二體、三體和四體作用能，發(fā)現(xiàn)三體作用對總作用能的貢獻可高達25%，四體作用在總作用能中所占的比例不超過3%，五體及以上多體作用在總作用能中所占的比例不超過0.5%.對67個不同尺寸水團簇總?cè)w作用能的計算結(jié)果表明，由30個水分子構(gòu)成的團簇的總?cè)w作用能可高達?209.87 kJ/mol［11］.李書實等［12］對酰胺、二肽、尿嘧啶與水分子形成的氫鍵復(fù)合物的三體效應(yīng)的研究表明，三體作用能占體系總作用能的10%~20%.以上研究結(jié)果表明，含水體系中的三體作用能對體系穩(wěn)定性有重要影響.

使用高精度從頭計算方法可以準確計算得到三體作用能.但是由于三體作用數(shù)目隨體系中分子數(shù)目的增多而急劇增大，致使使用高精度從頭計算方法計算三體作用能變得十分昂貴.以水分子團簇為例，對于含有10個水分子的體系，共計有120個三體作用；當(dāng)體系中水分子的個數(shù)增加到30個時，三體作用數(shù)目急速增加到4060個.因此，發(fā)展一種能夠快速準確計算氫鍵復(fù)合物體系三體作用強度的方法變得十分重要［13］.

在蛋白質(zhì)水體系中，極性氨基酸側(cè)鏈可與水分子形成具有不同氫鍵網(wǎng)絡(luò)結(jié)構(gòu)的團簇.這些團簇中存在大量的三體作用，并對生物分子的穩(wěn)定性產(chǎn)生重要影響.絲氨酸和蘇氨酸的側(cè)鏈均含有一個羥基，與甲醇（CH3OH）和乙醇（CH3CH2OH）的結(jié)構(gòu)相類似.本文使用甲醇、乙醇與水形成的氫鍵復(fù)合物來模擬絲氨酸和蘇氨酸側(cè)鏈與水分子間的作用，使用高精度的量子力學(xué)（QM）方法和可極化偶極-偶極作用模型計算體系中的三體作用能.通過擬合甲醇與水氫鍵復(fù)合物三體作用能隨分子間距離變化的能量曲線確定了可極化偶極-偶極作用模型所需的參數(shù).將可極化偶極-偶極作用模型和所確定的參數(shù)應(yīng)用于快速預(yù)測更多的甲醇、乙醇與水氫鍵復(fù)合物以及脫氧核糖與水氫鍵復(fù)合物的三體作用能，檢驗了模型的準確性和參數(shù)的可轉(zhuǎn)移性.

1 研究對象和計算方法

1.1 研究對象

使用二階微擾理論的MP2/6-31+G（d，p）方法優(yōu)化得到了1個甲醇分子與2個水分子構(gòu)成的氫鍵三聚體CH3OH·（H2O）2、2個甲醇分子與1個水分子構(gòu)成的氫鍵三聚體（CH3OH）2·H2O、3個甲醇分子構(gòu)成的氫鍵三聚體（CH3OH）3、1個乙醇分子與2個水分子構(gòu)成的氫鍵三聚體CH3CH2OH·（H2O）2、2個乙醇分子與1個水分子構(gòu)成的氫鍵三聚體（CH3CH2OH）2·H2O、3個乙醇分子構(gòu)成的氫鍵三聚體（CH3CH2OH）3和1個脫氧核糖分子與2個水分子構(gòu)成的氫鍵三聚體C4H9O3CHO·（H2O）2的穩(wěn)定幾何構(gòu)型.采用常溫常壓（NPT，p=101.325 kPa，T=298 K）系綜進行分子動力學(xué)模擬，得到了1個甲醇分子與20個水分子構(gòu)成的團簇CH3OH·（H2O）20、1個乙醇分子與20個水分子構(gòu)成的團簇CH3CH2OH·（H2O）20和1個脫氧核糖分子與12個水分子構(gòu)成的團簇C4H9O3CHO·（H2O）12的隨機結(jié)構(gòu).團簇獲取過程如下：將CH3OH，CH3CH2OH和C4H9O3CHO對應(yīng)放入邊長為1.5，1.7和2.1 nm的立方水盒子中，分別在NPT（p=101.325 kPa，T=298 K）系綜下進行2 ns的分子動力學(xué)模擬，步長均設(shè)置為1 fs，分別隨機選取了2 ns時CH3OH及其距離最近的20個水分子、1 ns時CH3CH2OH及其距離最近的20個水分子和800 ps時C4H9O3CHO及其距離最近的12個水分子，依次得到了團簇CH3OH·（H2O）20，CH3CH2OH·（H2O）20和C4H9O3CHO·（H2O）12的結(jié)構(gòu).

1.2 QM方法

在n個分子所構(gòu)成的多聚體中，由第i，j，k個分子所構(gòu)成的三聚體的三體作用能（V3B，kJ/mol）定義如下［14］：

總?cè)w作用能（IE3B，kJ/mol）為體系中所有三體作用能之和［15］，即

式中：E（i，j，k）（kJ/mol），E（i，j）（kJ/mol）和E（i）（kJ/mol）分別為三聚體、二聚體和單體的能量.所有甲醇、乙醇與水氫鍵復(fù)合物的三聚體、二聚體和單體的能量均由包含基組重疊誤差校正［16］的MP2/aug-ccpVTZ方法計算獲得，脫氧核糖與水氫鍵復(fù)合物的三聚體、二聚體和單體的能量均由包含基組重疊誤差校正的MP2/aug-cc-pVDZ方法計算獲得.全部計算均使用Gaussian 09［17］程序包完成.

1.3 基于可極化偶極-偶極作用模型計算三體作用能

根據(jù)可極化偶極-偶極作用模型（PBFF）［18，19］，對于由n個單體構(gòu)成的體系，極化作用對體系總能量的貢獻［Epol（1，···，n），kJ/mol］可表示為

式中：δμ=μ?μ0為某化學(xué)鍵的誘導(dǎo)偶極矩；μ（C·m）為團簇或復(fù)合物體系中單體的某個極性化學(xué)鍵的鍵偶極矩；μ0（C·m）為該單體處于孤立狀態(tài)時的這個極性化學(xué)鍵的鍵偶極矩；μ0，δμ和μ′0，δμ′分屬于不同的單體分子.θ，θ'，ζ，r是與兩個鍵偶極相對位置有關(guān)的幾何參量［20］.式（3）的求和遍及不同單體間所有偶極-偶極作用.

對于孤立分子，不存在環(huán)境或其它分子對它的影響，式（3）變?yōu)?/p>

對于由單體i和j兩個分子構(gòu)成的體系，可以將一個分子看成為另一個分子的環(huán)境，對i，j兩個單體間的所有偶極-偶極作用求和后，將得到兩個分子間的極化作用Epol(i，j).同理，對i，j，k3個單體間的所有偶極-偶極作用求和后，將得到3個分子間的極化作用Epol(i，j，k).

所以，可極化偶極-偶極作用模型下的三體作用能的計算如下：

誘導(dǎo)偶極矩δμ使用下式計算：

式中：q0（C）為孤立分子的原子分數(shù)電荷；q（C）為復(fù)合物中相對應(yīng)的原子分數(shù)電荷；d（nm）為化學(xué)鍵的鍵長；由于不同理論計算方法計算的原子分數(shù)電荷與其真實值間存在差異，為了得到與真實值更為接近的理論計算結(jié)果，在式（6）中加入了電荷修正系數(shù)c，其數(shù)值依賴于計算原子分數(shù)電荷時所使用的方法，c的值應(yīng)≥0.c＞1表示計算原子分數(shù)電荷所用的方法低估了復(fù)合物與孤立分子中原子分數(shù)電荷的差異，c<1則表示方法高估了復(fù)合物與孤立分子中原子分數(shù)電荷的差異.

1.4 參數(shù)的確定

使用可極化偶極-偶極作用模型計算體系的三體作用能需已知θ，θ'，ζ，r，d，q0，q，μ0和c等9個量，其中，θ，θ'，ζ，r，d的值可通過體系中原子的坐標計算得到，q和q0由AM1方法計算獲得.對于本文所研究的醇-水體系，水分子的O—H鍵偶極的固有偶極矩μ0和電荷修正系數(shù)c直接沿用我們前期工作確定的值，即μ0（O—H）=5.04×10?30C·m［21］，c=1.86［22］.甲醇分子含有3種化學(xué)鍵，分別是極性較強的C—O和O—H鍵以及弱極性化學(xué)鍵C—H鍵，乙醇分子較甲醇分子多一個非極性化學(xué)鍵C—C鍵.本文中忽略非極性化學(xué)鍵C—C鍵和弱極性化學(xué)鍵C—H鍵的鍵偶極矩，醇分子中的O—H鍵的固有偶極矩參數(shù)和電荷修正系數(shù)取值與水分子的相同.因此只需要確定極性化學(xué)鍵C—O的固有偶極矩參數(shù)和電荷修正系數(shù).

C—O鍵的固有偶極矩參數(shù)μ0和電荷修正系數(shù)c的確定方法如下：選取如圖1所示的（CH3OH）2·H2O氫鍵復(fù)合物為模型分子，在參與形成氫鍵的兩個O原子中心設(shè)置一個虛原子，改變虛原子與另一個O原子間的距離R，分別使用包含基組重疊誤差校正的MP2/aug-cc-pVTZ方法和本文方法計算不同距離R所對應(yīng)結(jié)構(gòu)的三體作用能，可獲得三體作用能隨分子間距離變化的能量曲線.調(diào)整C—O鍵的參數(shù)μ0和c的值，直至采用本文方法計算的三體作用能隨分子間距離變化的能量曲線與MP2/aug-cc-pVTZ方法計算結(jié)果的誤差最小，從而確定了醇分子C—O鍵的固有偶極矩參數(shù)μ0（C—O）=6.67×10?30C·m和電荷修正系數(shù)c（C—O）=0.根據(jù)式（6），c的值為0意味著不考慮C—O偶極的極化.

Fig.1 Structures of the trimers in training set

2 結(jié)果與討論

2.1 三體作用能及其隨分子間距離變化的計算

為了檢測可極化偶極-偶極作用模型的準確性和參數(shù)的可轉(zhuǎn)移性，將其應(yīng)用于計算甲醇、乙醇和脫氧核糖與水氫鍵復(fù)合物的三體作用能.圖2為由MP2/6-31+G（d，p）方法優(yōu)化得到的甲醇、乙醇和脫氧核糖與水分子形成的三聚體的穩(wěn)定結(jié)構(gòu)，其中結(jié)構(gòu)A~C為CH3OH·（H2O）2，結(jié)構(gòu)D和E為（CH3OH）2·H2O，結(jié)構(gòu)F為（CH3OH）3，結(jié)構(gòu)G為CH3CH2OH·（H2O）2，結(jié)構(gòu)H為（CH3CH2OH）2·H2O，結(jié)構(gòu)I為（CH3CH2OH）3，結(jié)構(gòu)J和K為C4H9O3CHO·（H2O）2.

Fig.2 Structures of trimers in testing set

由圖2可知，CH3OH·（H2O）2的3個異構(gòu)體之間和（CH3OH）2·H2O的2個異構(gòu)體之間的結(jié)構(gòu)既相似又略有不同，相似點在于，這些三聚體均形成3條氫鍵，且氫鍵鍵長差別不大；不同點在于，甲醇分子中的—CH3和水分子中未參與形成氫鍵的H原子取向存在一定差異.從圖2還可見，這11個三聚體中，每個三聚體均形成3條O—H…O型氫鍵，且每個分子既作質(zhì)子供體，又做質(zhì)子受體，因此推測體系的三體作用將表現(xiàn)為吸引，即三體作用能小于零.

表1 列出了使用MP2/aug-cc-pVTZ方法計算的甲醇-水三聚體和乙醇-水三聚體的三體作用能以及使用MP2/aug-cc-pVDZ方法計算的脫氧核糖-水三聚體的三體作用能.為了比較，表1還給出了PBFF的計算結(jié)果.

從表1可見，這11個三聚體的三體作用能均為負值，即三體作用表現(xiàn)為吸引，與我們的推斷一致，即，當(dāng)所有分子既作氫鍵供體又作氫鍵受體時，體系表現(xiàn)協(xié)同性；對比6個甲醇-水三聚體（A~F）的三體作用能，發(fā)現(xiàn)三體作用強度隨三聚體中所含甲醇數(shù)目的增多而增強，對比3個乙醇-水三聚體（G~I）的三體作用能也可得到類似的結(jié)論.

比較MP2方法和PBFF模型計算的三體作用能，可見PBFF模型計算的三體作用能與MP2方法結(jié)果的最大絕對誤差（Maximum absolute error）出現(xiàn)在三聚體H，即（CH3CH2OH）2·H2O，使用MP2/aug-ccpVTZ方法計算該三聚體的三體作用能為?10.25 kJ/mol，PBFF模型計算的結(jié)果為?8.87 kJ/mol，兩種方法結(jié)果間的誤差僅為1.38 kJ/mol；對于這11個三聚體的三體作用能，兩種方法計算結(jié)果的均方根誤差（Root-mean square error）為0.77 kJ/mol.以上結(jié)果說明，PBFF模型較好地預(yù)測了這11個三聚體的三體作用能.

為了驗證PBFF模型對甲醇、乙醇和脫氧核糖與水氫鍵復(fù)合物不穩(wěn)定結(jié)構(gòu)的適用性，將其應(yīng)用于計算甲醇-水、乙醇-水和脫氧核糖-水三聚體的三體作用能隨分子間距離變化的能量曲線.

對甲醇、乙醇與水分子形成的三聚體，在參與形成氫鍵的兩個O原子中心放置一個虛原子，改變虛原子與另一個O原子間的距離（R），計算三體作用能，可得到三體作用能隨分子間距離變化的能量曲線（圖3）.對脫氧核糖與兩個水分子形成的三聚體，在兩個水分子的O原子的中心設(shè)置一個虛原子，分別改變脫氧核糖分子中參與形成氫鍵的O原子與兩個水分子中O原子的距離以及與虛原子間的距離，可得到三組三體作用能隨分子間距離變化的能量曲線（圖4）.

Table 1 Three-body interaction energies(IE3B)calculated by different methods for the testing trimers A—K in Fig.2

Fig.3 Three?body interaction energy(V3B)vs.the intermolecular distances(R)for the testing trimers A—I

Fig.4 Three?body interaction energy vs.the intermolecular distances for the testing trimers J and K

由圖3可知，對于甲醇-水三聚體和乙醇-水三聚體而言，三體作用能隨分子間距離變化的能量曲線存在一定的相似性，即，當(dāng)分子間距離R介于0.21~0.48 nm之間時，三體作用能介于?16.00~?0.50 kJ/mol之間.綜合圖3和圖4可知，對甲醇、乙醇和脫氧核糖與水分子形成的三聚體而言，隨著分子間距離R的變大，三體作用強度逐漸減弱，由此可推測，三體作用為近程作用；PBFF模型計算得到的能量曲線不僅較好地描述了三體作用能隨分子間距離變化的趨勢，而且模型得到的三體作用能與MP2方法的計算結(jié)果符合得很好，說明PBFF模型可準確預(yù)測三體作用能隨三聚體結(jié)構(gòu)的變化，進一步驗證了PBFF模型的準確性.

2.2 團簇體系總?cè)w作用能的計算

相較于三聚體，甲醇、乙醇、脫氧核糖與多個水分子形成的團簇的結(jié)構(gòu)和性質(zhì)與分子存在于水溶液中的情況更為接近.為了進一步驗證模型的準確性，將PBFF模型應(yīng)用于計算CH3OH·（H2O）20，CH3CH2OH·（H2O）20和C4H9O3CHO·（H2O）12的總?cè)w作用能，為了比較，還給出了AMOEBAbio18方法的計算結(jié)果，其中AMOEBAbio18方法的結(jié)果使用Tinker程序包計算得到.圖5為這3個團簇的結(jié)構(gòu)圖.表2列出了采用MP2方法、AMOEBAbio18方法和PBFF模型計算的總?cè)w作用能結(jié)果，其中，第4列為以MP2結(jié)果作為基準，AMOEBAbio18方法計算結(jié)果的誤差（Δ），第6列為以MP2結(jié)果作為基準PBFF模型計算結(jié)果的誤差（Δ）.

Fig.5 Structures of CH3OH·(H2O)20(A),CH3CH2OH·(H2O)20(B)and C4H9O3CHO·(H2O)12(C)

Table 2 Total three-body interaction energies(IE3B)calculated by different methods

由表2可知，對于團簇CH3OH·（H2O）20和CH3CH2OH·（H2O）20，使用MP2/aug-cc-pVTZ方法計算的總?cè)w作用能分別為?59.25和?35.82 kJ/mol，使用AMOEBAbio18方法計算的總?cè)w作用能分別為?55.52和?27.51 kJ/mol，誤差分別為?3.73和?8.31 kJ/mol，使用PBFF模型計算的總?cè)w作用能分別為?59.41和?35.86 kJ/mol，誤差僅為0.16和0.04 kJ/mol；對于團簇C4H9O3CHO·（H2O）12，使用MP2/aug-cc-pVDZ方法計算的總?cè)w作用能為?7.15 kJ/mol，使用AMOEBAbio18方法計算的三體作用能為?8.42 kJ/mol，誤差為1.27 kJ/mol，使用PBFF方法計算的總?cè)w作用能為?6.57 kJ/mol，誤差僅為?0.58 kJ/mol.上述結(jié)果說明，對于這3個大尺寸團簇，PBFF模型計算的總?cè)w作用能與MP2方法的計算結(jié)果符合得很好，模型的計算精度優(yōu)于AMOEBAbio18方法，進一步證明了模型的準確性.

3 結(jié) 論

對于醇分子，僅考慮極性較強的C—O和O—H鍵的鍵偶極矩，忽略非極性化學(xué)鍵C—C鍵和弱極性化學(xué)鍵C—H鍵的鍵偶極矩，通過擬合甲醇水三聚體的三體作用能，確定了可極化偶極-偶極作用模型所需參數(shù).將模型及參數(shù)應(yīng)用于計算更多的甲醇、乙醇和脫氧核糖與水氫鍵復(fù)合物的三體作用能及其隨分子間距離變化的能量曲線.計算結(jié)果表明，可極化偶極-偶極作用模型和本文所確定的參數(shù)可準確預(yù)測具有不同結(jié)構(gòu)三聚體的三體作用能以及較大團簇的總?cè)w作用能，計算精度與MP2方法的精度相當(dāng).目前，可極化偶極-偶極作用模型已經(jīng)可以準確地預(yù)測水團簇的三體作用強度、醇及脫氧核糖與水分子間的三體作用強度，期望為準確快速模擬水環(huán)境中蛋白質(zhì)和核酸的結(jié)構(gòu)和性質(zhì)提供有力的工具.