QM/MM方法的研究*

劉景林，曹亞杰，吳云飛

(佳木斯大學(xué))

0 引言

分子動(dòng)力學(xué)方法是較成熟的模擬方法之一，其可以對(duì)由千萬(wàn)個(gè)原子構(gòu)成的生物大分子體系進(jìn)行高效地模擬計(jì)算．但由于其是以經(jīng)典力學(xué)為基礎(chǔ)的，無(wú)法充分地描述電子的運(yùn)動(dòng)，在二十世紀(jì)初誕生的量子力學(xué)能夠充分的地描述電子運(yùn)動(dòng)，現(xiàn)今應(yīng)用量子力學(xué)方法已經(jīng)可以對(duì)電子運(yùn)動(dòng)進(jìn)行精確的計(jì)算，但是由于其需要進(jìn)行大量的并且及其復(fù)雜的積分運(yùn)算，其結(jié)果就是產(chǎn)生了巨大的運(yùn)算量，并且嚴(yán)重的增加了計(jì)算成本．即使在現(xiàn)今最優(yōu)秀的計(jì)算環(huán)境下，也無(wú)法很好地進(jìn)行這種大規(guī)模的運(yùn)算，這種大規(guī)模的運(yùn)算在化學(xué)、生物學(xué)等多門(mén)學(xué)科的發(fā)展中占有重要的地位，為了可以高效地進(jìn)行這種大規(guī)模的運(yùn)算，既能得到精確的結(jié)果，同時(shí)可以有效地降低計(jì)算成本，則應(yīng)運(yùn)而生了被稱作QM/MM方法的混合動(dòng)力學(xué)計(jì)算方法，這是采用量子力學(xué)(Quantum Mechanics，QM)與分子動(dòng)力學(xué)(Molecular Mechanics，MM)相結(jié)合的方法，對(duì)生物大分子體系中需要仔細(xì)觀察的重點(diǎn)部位中所包含的少數(shù)的原子使用QM方法進(jìn)行精確的計(jì)算，而體系中剩余的部分采用MM方法進(jìn)行計(jì)算．這種方法提高計(jì)算結(jié)果的精確性，同時(shí)可以有效地降低計(jì)算成本，因此正被更多的科研人員所采用．

1 量子力學(xué)(Quantum Mechanics，QM)方法概述

量子力學(xué)的方法是以分子中電子的非定域化為基礎(chǔ)，電子的運(yùn)動(dòng)使用函數(shù)描述，依據(jù)海森伯的測(cè)不準(zhǔn)原理，計(jì)算在非定域區(qū)間內(nèi)電子出現(xiàn)的幾率．量子力學(xué)方法中最為普遍的計(jì)算方法為從頭算方法，而所謂的從頭算方法又被稱為分子軌域計(jì)算方法，指僅僅采用了分子軌道理論的三個(gè)基本近似，在此基礎(chǔ)上應(yīng)用變分原理求解分子的Roothaan-Hall方程的方法，把描述體系電子狀態(tài)的波函數(shù)展開(kāi)為原子軌域波函數(shù)的線性組合．原子軌域的波函數(shù)卻為某些特定數(shù)學(xué)函數(shù)(如高斯函數(shù)等)的線性組合．在核運(yùn)動(dòng)和電子運(yùn)動(dòng)分離的近似條件下，可以得到多原子分子體系的電子的哈密頓量表達(dá)式如下:

則其薛定諤方程為:

上式都采用原子單位，角標(biāo)α與β表示的是原子核，i與j表示的是電子．從頭算方法十分精確，但是計(jì)算速度特別緩慢，導(dǎo)致能計(jì)算系統(tǒng)的尺度極其有限，最多不會(huì)超過(guò)100個(gè)原子［1］．為了增加量子力學(xué)方法計(jì)算的效益，自1960年起，引進(jìn)了一些實(shí)驗(yàn)值作為參數(shù)，去替代某些積分項(xiàng)，這種方法被稱作半經(jīng)驗(yàn)分子軌道方法，但還是無(wú)法對(duì)大體系進(jìn)行計(jì)算．

2 分子動(dòng)力學(xué)(Molecular Mechanics，MM)方法概述

分子動(dòng)力學(xué)理論基本的原理是在一個(gè)由N個(gè)粒子組成的系統(tǒng)中的每個(gè)粒子的運(yùn)動(dòng)都遵循Newton方程［2］;通過(guò)對(duì)每一粒子和其周?chē)牧Ｗ娱g相互作用勢(shì)的疊加求得粒子間相互作用勢(shì)［3］．此方法先使用各粒子的位置坐標(biāo)計(jì)算出系統(tǒng)勢(shì)能，再由公式(3)、(4)計(jì)算出系統(tǒng)中各個(gè)粒子受到的力以及其加速度．

只要給出其初始條件，求解體系里各個(gè)粒子運(yùn)動(dòng)方程獲取每一個(gè)粒子在不同時(shí)刻的運(yùn)動(dòng)軌跡，就能夠憑借統(tǒng)計(jì)學(xué)方法去了解體系的狀態(tài)隨時(shí)間變化的規(guī)律．

分子動(dòng)力學(xué)方法由于使用了保守力場(chǎng)，而它是原子位置坐標(biāo)的函數(shù)，這表明了電子的運(yùn)動(dòng)沒(méi)有被考慮，則電子遷移的過(guò)程和電子激發(fā)態(tài)無(wú)法被處理，因此也就不可能去描述化學(xué)鍵的斷裂和生成了．

3 QM/MM方法的實(shí)現(xiàn)

將量子力學(xué)方法與分子動(dòng)力學(xué)方法有機(jī)的結(jié)合起來(lái)，可以有效地避免它們各自缺點(diǎn)，發(fā)揮其優(yōu)點(diǎn)．由此在上世紀(jì)70年代發(fā)展了QM/MM方法，而其實(shí)現(xiàn)的基本原理為ONIOM方法．

3．1 ONIOM方法的基本思想

Morokuma等人提出了ONIOM(Our own nlayered integrated molecular orbital and molecular mechanics)方法［4］，其基本思想是將體系按照其重要性的不同劃分成了若干層，每層需要使用不同等級(jí)的算法，實(shí)現(xiàn)了用多級(jí)算法進(jìn)行模擬計(jì)算大體系，這樣把大體系給層層分割了．下面以三層ONIOM體系來(lái)進(jìn)行分析．

如圖1體系被分成了三層，分別為第一層活躍區(qū)，第二層較活躍區(qū)及第三層不活躍區(qū)．然后再定義三個(gè)新體系．即“內(nèi)核”體系只包含第一層，“媒介”體系包含了第一層和第二層，而“實(shí)際”體系包含了全部三層，可采用三種層級(jí)算法．例如，“實(shí)際”體系可采用經(jīng)典的分子動(dòng)力學(xué)方法(低層級(jí)算法)，“媒介”體系可采用半經(jīng)驗(yàn)量子化學(xué)方法(中層級(jí)算法)，“內(nèi)核”體系可采用最精確的從頭算量子化學(xué)方法(高層級(jí)算法)．

圖1 三層ONIOM方法的示意圖

圖2 三層ONIOM算法的示意圖

圖2上的各點(diǎn)代表采用不同的層級(jí)算法計(jì)算得到對(duì)應(yīng)的體系能量，而F點(diǎn)是采用高層級(jí)算法計(jì)算得到的“實(shí)際”體系能量EF．如果直接使用量子力學(xué)方法去計(jì)算EF是特別困難的，而ONIOM方法是通過(guò)計(jì)算A至E五點(diǎn)的能量后近似地獲得EF，即

可以看出(EC-EB)和(EE-ED)實(shí)際上就是上面定義第二層和第三層能量．計(jì)算相鄰兩個(gè)體系的能量差表示各對(duì)應(yīng)層能量是為了保證體系結(jié)構(gòu)的完整性．因此對(duì)于被分成n層后的體系使用n個(gè)不同層級(jí)算法可以得到體系的總能量如下:

其中E［Level(i)，Model(j)］表示使用第i層級(jí)算法計(jì)算第j個(gè)體系得到的該體系能量．Level(1)為最高層級(jí)的算法，Model(1)為包含了所有粒子的“實(shí)際”體系．E(ONIOMn)的實(shí)質(zhì)為使用n種層級(jí)算法得到了E［Level1(1)，Model(1)］的近似值．

3．2 QM/MM體系的建立

該文中使用雙層的ONIOM方法作為例子介紹了QM/MM方法體系的建立．如圖3所示，將體系劃分為層Ⅰ和層Ⅱ．

圖3 體系分層示意圖

圖3的層Ⅰ與層Ⅱ是通過(guò)加入了鏈接原子(Link Atom，LA)相連．那么對(duì)“實(shí)際”體系使用了次層級(jí)算法分子動(dòng)力學(xué)(MM)去計(jì)算，而對(duì)層Ⅰ構(gòu)成的“內(nèi)核”體系使用高層級(jí)算法量子力學(xué)(QM)去計(jì)算．顯然QM/MM算法得到體系的總能量為［5］

式(7)中的Eqm［6］ 和 Emm

［7］分別為

Pμν為密度矩陣元，Hμν為單電子對(duì)應(yīng)原子軌道的哈密頓量矩陣元，F(xiàn)μ，ν為Fock矩陣元，為核的排斥能．且

用MOPAC軟件［8］中的AM1半經(jīng)驗(yàn)量子化學(xué)方法［9］，選用的是STO-3G基組計(jì)算QM區(qū)域量子化學(xué)能量Eqm(qm+L)，使用GROMACS軟件［10］選用OPLS全原子力場(chǎng)去計(jì)算QM區(qū)域的能量Emm(qm+L)及QM和MM兩區(qū)域的總能量Emm(total)．

3．3 層邊界的處理

如何處理層邊界，不同的模擬軟件采用不同的方法．在GROMACS軟件中是把兩層邊界的化學(xué)鍵人為斷開(kāi)，將兩端組成該化學(xué)鍵的原子分別劃歸為層Ⅰ與層Ⅱ后再加入一個(gè)通常都是虛擬的氫原子的鏈接原子．鏈接原子的質(zhì)量和所帶電量都為零，因此對(duì)計(jì)算的最終結(jié)果沒(méi)影響，而在計(jì)算模擬的時(shí)候，把它放到使用量子力學(xué)計(jì)算的層Ⅰ中，鏈接原子放在斷開(kāi)的成鍵的兩個(gè)原子之間，并且這三個(gè)原子在一條直線上．例如在OPLASS力場(chǎng)中，斷開(kāi)了鍵長(zhǎng)為0．153 nm的C—C鍵，用“虛擬的氫原子”來(lái)作鏈接原子，C—H鍵長(zhǎng)為0．108 nm，則鏈接原子與層Ⅰ的碳原子之間的距離就是0．108 nm，而鏈接原子與層Ⅱ的碳原子之間距離為0．045 nm，通常令鏈接原子靠近層Ⅱ那端的原子，兩個(gè)碳原子的坐標(biāo)分 別 為 (x1，y1，z1) 和 (x2，y2，z2)，C—H 鍵 和C—C 鍵的鍵長(zhǎng)之比是0．108/0．153=0．706．那么鏈接原子坐標(biāo)是(x0，y0，z0)，分別為:x0=x1+|x1-x2|× 0．706，y0=y1+|y1-y2|× 0．706，z0=z1+|z1-z2|×0．706，其中角標(biāo)為1的是層Ⅰ的碳原子，角標(biāo)為2的是層Ⅱ的碳原子．

4 結(jié)束語(yǔ)

目前，QM/MM方法已經(jīng)用于溶液體系中的化學(xué)變化、生物大分子的處理、酶催化原理等方面，得到了很好的發(fā)展．在未來(lái)的科研工作中，該方法會(huì)越來(lái)越受到重視．

［1］ 陳光巨，黃元河．量子化學(xué)．上海:華東理工大學(xué)出版社．

［2］ Alder B J，Wainwright T E．Phase transition for a hardsphere system［J］．J Chem Phys，1957，27:1208-1209．

［3］ 徐筱杰，侯廷軍，喬學(xué)斌，等．計(jì)算機(jī)輔助藥物分子設(shè)計(jì)［M］．北京:化學(xué)工業(yè)出版社．

［4］ Morokuma K，Svensson M，Humbel S，et al．ONIOM:A Multilayered Integrated MO+MM Method for Geometry Optimizations and Single Point Energy Predictions［J］．J Phys Chem，1996，100:19357-19363．

［5］ Obara S，Saika A．Efficient recursive computation of molecular integrals over Cartesian Gaussian functions［J］．J Chem Phys，1986，84，7(1):3963-3974．

［6］ 徐光憲，黎樂(lè)民，王德民．量子化學(xué):中冊(cè)［M］．北京:科學(xué)出版社，2009．

［7］ William L，Jorgensen，David S．Maxwell，Julian Tirado-Rives．Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids［J］．J Am Chem Soc，1996，118，11225-11236．

［8］ Dewar M J S．Development and status of MINDO/3 and MNDO［J］．J Mol Struct，1983，100:41．

［9］ Hess B，Bekker H，Berendsen H J C，et al．LINCS:A linear constraint solver for molecular simulations［J］．J Comp Chem，1997，18:1463-1472．

［10］ Bussi G，Donadio D，Parrinello M．Canonical sampling through velocity rescaling［J］．J Chem Phys，2007，126:014101-014101-7．