醇及脫氧核糖與水分子間三體作用強(qiáng)度的快速準(zhǔn)確計(jì)算

李曉蕾，孫云嬌，唐 穎，王長(zhǎng)生

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

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

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

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

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

1 研究對(duì)象和計(jì)算方法

1.1 研究對(duì)象

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

1.2 QM方法

在n個(gè)分子所構(gòu)成的多聚體中，由第i，j，k個(gè)分子所構(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方法計(jì)算獲得，脫氧核糖與水氫鍵復(fù)合物的三聚體、二聚體和單體的能量均由包含基組重疊誤差校正的MP2/aug-cc-pVDZ方法計(jì)算獲得.全部計(jì)算均使用Gaussian 09［17］程序包完成.

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

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

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

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

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

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

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

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

1.4 參數(shù)的確定

使用可極化偶極-偶極作用模型計(jì)算體系的三體作用能需已知θ，θ'，ζ，r，d，q0，q，μ0和c等9個(gè)量，其中，θ，θ'，ζ，r，d的值可通過(guò)體系中原子的坐標(biāo)計(jì)算得到，q和q0由AM1方法計(jì)算獲得.對(duì)于本文所研究的醇-水體系，水分子的O—H鍵偶極的固有偶極矩μ0和電荷修正系數(shù)c直接沿用我們前期工作確定的值，即μ0（O—H）=5.04×10?30C·m［21］，c=1.86［22］.甲醇分子含有3種化學(xué)鍵，分別是極性較強(qiáng)的C—O和O—H鍵以及弱極性化學(xué)鍵C—H鍵，乙醇分子較甲醇分子多一個(gè)非極性化學(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ù)合物為模型分子，在參與形成氫鍵的兩個(gè)O原子中心設(shè)置一個(gè)虛原子，改變虛原子與另一個(gè)O原子間的距離R，分別使用包含基組重疊誤差校正的MP2/aug-cc-pVTZ方法和本文方法計(jì)算不同距離R所對(duì)應(yīng)結(jié)構(gòu)的三體作用能，可獲得三體作用能隨分子間距離變化的能量曲線.調(diào)整C—O鍵的參數(shù)μ0和c的值，直至采用本文方法計(jì)算的三體作用能隨分子間距離變化的能量曲線與MP2/aug-cc-pVTZ方法計(jì)算結(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 三體作用能及其隨分子間距離變化的計(jì)算

為了檢測(cè)可極化偶極-偶極作用模型的準(zhǔn)確性和參數(shù)的可轉(zhuǎn)移性，將其應(yīng)用于計(jì)算甲醇、乙醇和脫氧核糖與水氫鍵復(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è)異構(gòu)體之間和（CH3OH）2·H2O的2個(gè)異構(gòu)體之間的結(jié)構(gòu)既相似又略有不同，相似點(diǎn)在于，這些三聚體均形成3條氫鍵，且氫鍵鍵長(zhǎng)差別不大；不同點(diǎn)在于，甲醇分子中的—CH3和水分子中未參與形成氫鍵的H原子取向存在一定差異.從圖2還可見(jiàn)，這11個(gè)三聚體中，每個(gè)三聚體均形成3條O—H…O型氫鍵，且每個(gè)分子既作質(zhì)子供體，又做質(zhì)子受體，因此推測(cè)體系的三體作用將表現(xiàn)為吸引，即三體作用能小于零.

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

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

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

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

對(duì)甲醇、乙醇與水分子形成的三聚體，在參與形成氫鍵的兩個(gè)O原子中心放置一個(gè)虛原子，改變虛原子與另一個(gè)O原子間的距離（R），計(jì)算三體作用能，可得到三體作用能隨分子間距離變化的能量曲線（圖3）.對(duì)脫氧核糖與兩個(gè)水分子形成的三聚體，在兩個(gè)水分子的O原子的中心設(shè)置一個(gè)虛原子，分別改變脫氧核糖分子中參與形成氫鍵的O原子與兩個(gè)水分子中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可知，對(duì)于甲醇-水三聚體和乙醇-水三聚體而言，三體作用能隨分子間距離變化的能量曲線存在一定的相似性，即，當(dāng)分子間距離R介于0.21~0.48 nm之間時(shí)，三體作用能介于?16.00~?0.50 kJ/mol之間.綜合圖3和圖4可知，對(duì)甲醇、乙醇和脫氧核糖與水分子形成的三聚體而言，隨著分子間距離R的變大，三體作用強(qiáng)度逐漸減弱，由此可推測(cè)，三體作用為近程作用；PBFF模型計(jì)算得到的能量曲線不僅較好地描述了三體作用能隨分子間距離變化的趨勢(shì)，而且模型得到的三體作用能與MP2方法的計(jì)算結(jié)果符合得很好，說(shuō)明PBFF模型可準(zhǔn)確預(yù)測(cè)三體作用能隨三聚體結(jié)構(gòu)的變化，進(jìn)一步驗(yàn)證了PBFF模型的準(zhǔn)確性.

2.2 團(tuán)簇體系總?cè)w作用能的計(jì)算

相較于三聚體，甲醇、乙醇、脫氧核糖與多個(gè)水分子形成的團(tuán)簇的結(jié)構(gòu)和性質(zhì)與分子存在于水溶液中的情況更為接近.為了進(jìn)一步驗(yàn)證模型的準(zhǔn)確性，將PBFF模型應(yīng)用于計(jì)算CH3OH·（H2O）20，CH3CH2OH·（H2O）20和C4H9O3CHO·（H2O）12的總?cè)w作用能，為了比較，還給出了AMOEBAbio18方法的計(jì)算結(jié)果，其中AMOEBAbio18方法的結(jié)果使用Tinker程序包計(jì)算得到.圖5為這3個(gè)團(tuán)簇的結(jié)構(gòu)圖.表2列出了采用MP2方法、AMOEBAbio18方法和PBFF模型計(jì)算的總?cè)w作用能結(jié)果，其中，第4列為以MP2結(jié)果作為基準(zhǔn)，AMOEBAbio18方法計(jì)算結(jié)果的誤差（Δ），第6列為以MP2結(jié)果作為基準(zhǔn)PBFF模型計(jì)算結(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可知，對(duì)于團(tuán)簇CH3OH·（H2O）20和CH3CH2OH·（H2O）20，使用MP2/aug-cc-pVTZ方法計(jì)算的總?cè)w作用能分別為?59.25和?35.82 kJ/mol，使用AMOEBAbio18方法計(jì)算的總?cè)w作用能分別為?55.52和?27.51 kJ/mol，誤差分別為?3.73和?8.31 kJ/mol，使用PBFF模型計(jì)算的總?cè)w作用能分別為?59.41和?35.86 kJ/mol，誤差僅為0.16和0.04 kJ/mol；對(duì)于團(tuán)簇C4H9O3CHO·（H2O）12，使用MP2/aug-cc-pVDZ方法計(jì)算的總?cè)w作用能為?7.15 kJ/mol，使用AMOEBAbio18方法計(jì)算的三體作用能為?8.42 kJ/mol，誤差為1.27 kJ/mol，使用PBFF方法計(jì)算的總?cè)w作用能為?6.57 kJ/mol，誤差僅為?0.58 kJ/mol.上述結(jié)果說(shuō)明，對(duì)于這3個(gè)大尺寸團(tuán)簇，PBFF模型計(jì)算的總?cè)w作用能與MP2方法的計(jì)算結(jié)果符合得很好，模型的計(jì)算精度優(yōu)于AMOEBAbio18方法，進(jìn)一步證明了模型的準(zhǔn)確性.

3 結(jié) 論

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