QM/MM方法的研究*

劉景林，曹亞杰，吳云飛

(佳木斯大學)

0 引言

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

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

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

則其薛定諤方程為:

上式都采用原子單位，角標α與β表示的是原子核，i與j表示的是電子．從頭算方法十分精確，但是計算速度特別緩慢，導致能計算系統的尺度極其有限，最多不會超過100個原子［1］．為了增加量子力學方法計算的效益，自1960年起，引進了一些實驗值作為參數，去替代某些積分項，這種方法被稱作半經驗分子軌道方法，但還是無法對大體系進行計算．

2 分子動力學(Molecular Mechanics，MM)方法概述

分子動力學理論基本的原理是在一個由N個粒子組成的系統中的每個粒子的運動都遵循Newton方程［2］;通過對每一粒子和其周圍的粒子間相互作用勢的疊加求得粒子間相互作用勢［3］．此方法先使用各粒子的位置坐標計算出系統勢能，再由公式(3)、(4)計算出系統中各個粒子受到的力以及其加速度．

只要給出其初始條件，求解體系里各個粒子運動方程獲取每一個粒子在不同時刻的運動軌跡，就能夠憑借統計學方法去了解體系的狀態隨時間變化的規律．

分子動力學方法由于使用了保守力場，而它是原子位置坐標的函數，這表明了電子的運動沒有被考慮，則電子遷移的過程和電子激發態無法被處理，因此也就不可能去描述化學鍵的斷裂和生成了．

3 QM/MM方法的實現

將量子力學方法與分子動力學方法有機的結合起來，可以有效地避免它們各自缺點，發揮其優點．由此在上世紀70年代發展了QM/MM方法，而其實現的基本原理為ONIOM方法．

3．1 ONIOM方法的基本思想

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

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

圖1 三層ONIOM方法的示意圖

圖2 三層ONIOM算法的示意圖

圖2上的各點代表采用不同的層級算法計算得到對應的體系能量，而F點是采用高層級算法計算得到的“實際”體系能量EF．如果直接使用量子力學方法去計算EF是特別困難的，而ONIOM方法是通過計算A至E五點的能量后近似地獲得EF，即

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

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

3．2 QM/MM體系的建立

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

圖3 體系分層示意圖

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

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

［7］分別為

Pμν為密度矩陣元，Hμν為單電子對應原子軌道的哈密頓量矩陣元，Fμ，ν為Fock矩陣元，為核的排斥能．且

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

3．3 層邊界的處理

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

4 結束語

目前，QM/MM方法已經用于溶液體系中的化學變化、生物大分子的處理、酶催化原理等方面，得到了很好的發展．在未來的科研工作中，該方法會越來越受到重視．

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

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

［3］ 徐筱杰，侯廷軍，喬學斌，等．計算機輔助藥物分子設計［M］．北京:化學工業出版社．

［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］ 徐光憲，黎樂民，王德民．量子化學:中冊［M］．北京:科學出版社，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．