謝湖均雷群芳方文軍

(1浙江工商大學(xué)應(yīng)用化學(xué)系 浙江杭州310035;2浙江大學(xué)化學(xué)系 浙江杭州310027)

2013年諾貝爾化學(xué)獎10月9日在瑞典揭曉，法國斯特拉斯堡大學(xué)和美國哈佛大學(xué)教授Martin Karplus、美國斯坦福大學(xué)醫(yī)學(xué)院教授Michael Levitt和美國南加州大學(xué)教授Arieh Warshel因發(fā)展復(fù)雜化學(xué)系統(tǒng)的多尺度模型而共享獎項(xiàng)［1-2］。以前化學(xué)家是用塑料的球和棍來搭建和創(chuàng)造分子模型，而現(xiàn)在則是用計(jì)算機(jī)來輔助建模和計(jì)算。分子和化學(xué)反應(yīng)的精確建模對于化學(xué)的進(jìn)步至關(guān)重要，化學(xué)反應(yīng)的速度非?？?，在幾分之一毫秒間，電子就會從一個原子核跳到另一個原子核，經(jīng)典化學(xué)在這里已無用武之地。Karplus、Levitt和Warshel工作的突破意義在于他們設(shè)法讓量子力學(xué)(quantum mechanics，QM)和分子力學(xué)(molecular mechanics，MM)結(jié)合在化學(xué)過程的建模之中。量子力學(xué)計(jì)算方法可以用來預(yù)測電子結(jié)構(gòu)和化學(xué)反應(yīng)機(jī)理，精確度很高，但只能用來計(jì)算較小的體系。而分子力學(xué)的優(yōu)勢在于計(jì)算簡便，雖然可以用來計(jì)算較大的復(fù)雜體系，但精確度不夠高，而且無法描述化學(xué)鍵生成或斷裂的化學(xué)反應(yīng)過程。他們3人的工作結(jié)合了兩者的長處，發(fā)展出量子力學(xué)和分子力學(xué)組合方法(combined quantum mechanics/molecular mechanics method，QM/MM)。QM/MM組合方法在大分子體系的計(jì)算研究中已展現(xiàn)出越來越強(qiáng)大的功能，已經(jīng)逐步應(yīng)用到化學(xué)、材料、生物學(xué)等各個相關(guān)學(xué)科領(lǐng)域［3-8］。本文就其基本原理和目前的研究進(jìn)展做簡單的介紹。

1 QM/MM組合方法基本原理

1.1 QM/MM組合方法的能量表達(dá)式

如圖1所示，QM/MM組合方法的主要思想是把整個體系分為QM和MM兩部分，其中QM部分用量子力學(xué)方法處理，MM部分用分子力學(xué)方法處理，QM和MM的邊界則用連接原子或凍軌道等方法處理。目前流行的QM/MM能量表達(dá)式有兩種，加和方法(additive schemes)［9-10］和減去方法(subtractive schemes)［11-13］。在一般的計(jì)算中采用的是加和方法，其能量表達(dá)式為:

這里的EQM和EMM分別為單獨(dú)QM和MM部分的能量。EQM-MM為兩個區(qū)域的耦合項(xiàng)。原理上任何的量子力學(xué)方法都可以用來處理QM部分，但是在文獻(xiàn)中報道的一般都是DFT和半經(jīng)驗(yàn)理論計(jì)算的結(jié)果。

這種能量的加和方法在QM/MM組合計(jì)算中廣泛采用，尤其是在生物大分子領(lǐng)域。目前可以使用的軟件包主要有CHARMM［14-15］和AMBER［16-17］。但是此方法也存在著問題，由于存在連接原子［18］，耦合項(xiàng)EQM-MM不易計(jì)算。

圖1 QM/MM組合方法示意圖

1.2 靜電勢項(xiàng)

QM區(qū)域電荷密度和MM區(qū)域電荷模型之間的靜電勢耦合作用，可以在不同的計(jì)算水平下進(jìn)行，主要的區(qū)別是QM和MM區(qū)域相互極化的作用范圍。為此，Bakowies和Thiel定義了3種處理靜電勢相互作用的方法［19］。

1.2.1 機(jī)械嵌入(mechanical embedding)

在這種方法中，QM區(qū)域的計(jì)算本質(zhì)上是在氣相中進(jìn)行的，缺乏與環(huán)境之間的耦合作用。QM與MM區(qū)域靜電勢的相互作用或者被遺漏，或者僅僅在MM水平上計(jì)算。計(jì)算時一般采用剛體的點(diǎn)電荷模型。對于其他的方法，例如鍵偶極方法也在QM部分使用。Morokuma提出的ONIOM模型即采用這種機(jī)械嵌入方案來處理QM與MM區(qū)域間的靜電耦合。

但是，在這種方法的應(yīng)用過程中，存在明顯的缺陷和限制。(1)外層的電荷不與QM區(qū)域的密度相互作用，使QM部分不能直接被靜電勢環(huán)境所影響;因此，QM區(qū)域的密度沒有被極化。(2)對于QM區(qū)域的電荷分布，例如在反應(yīng)過程中，電荷模型需要不斷地更新;但是，這又會導(dǎo)致勢能面的不連續(xù)性。(3)采用MM點(diǎn)電荷來處理內(nèi)層區(qū)域所產(chǎn)生的偏差是不能忽略的;一般在程序中，都有很多種力場來處理這些體系，選擇時需要慎重考慮，因?yàn)橐话懔龅陌l(fā)展與這些具體的化合物沒有聯(lián)系。(4)MM的電荷模型依靠其他的力場參數(shù)，這就意味著最后產(chǎn)生的構(gòu)型是平衡態(tài)的描述，而不是重新產(chǎn)生的真實(shí)電荷分布。

1.2.2 靜電場嵌入(electrostatic embedding)

在此方法中，MM區(qū)域電荷分布對QM區(qū)域所產(chǎn)生的極化作用，可以看作是QM區(qū)域電子結(jié)構(gòu)計(jì)算的一部分，因此能夠克服機(jī)械嵌入方法的缺點(diǎn)。其表達(dá)式為:

這里qM為MM區(qū)域點(diǎn)電荷，Zα為QM區(qū)域原子的核電荷，i是電子總數(shù)，M是總的點(diǎn)電荷，α為QM區(qū)域所有的核。

在靜電場嵌入方法中，內(nèi)層區(qū)域的電子結(jié)構(gòu)可以適應(yīng)環(huán)境電荷的變化，并且被環(huán)境所極化。這里的靜電勢是在QM水平下計(jì)算的，相對于機(jī)械嵌入方法，明顯提高了精度，當(dāng)然計(jì)算的代價也會更大。

使用靜電場嵌入方法，需要注意的地方是處理QM/MM的邊界區(qū)域，此時MM的電荷非常接近QM的電子密度，會導(dǎo)致過度極化問題，這在邊界區(qū)域是共價鍵時更加顯著。目前，靜電場嵌入是生物大分子計(jì)算中最流行的方法。

1.2.3 極化嵌入(polarized embedding)

此種方法在靜電場嵌入的基礎(chǔ)之上，又包含了QM區(qū)域電荷分布對MM區(qū)域所產(chǎn)生的極化影響。盡管極化嵌入是最精確的方法，但是對于它的應(yīng)用范圍，仍然是非常有限的。主要的問題在于還沒有建立起合適的生物大分子極化力場。目前，有許多極化的溶劑模型可以采用，最顯著的就是液態(tài)水的模擬。對于蛋白的極化力場，還在進(jìn)一步的發(fā)展中。當(dāng)然，這種方法的應(yīng)用，也會增加計(jì)算的代價，還有可能造成收斂問題。目前常用的極化嵌入方法有polarized point dipoles(PPD)，drude oscillators(DO)和fluctuating charges(FQ)。

1.3 其他的非鍵和成鍵相互作用

除了上面部分討論的靜電勢相互作用，還有范德華和成鍵相互作用貢獻(xiàn)到QM/MM的耦合項(xiàng)。它們的處理比較簡單，一般都在分子力學(xué)的計(jì)算水平下。范德華相互作用一般采用經(jīng)典的Lennard-Jones勢:

Friesner等發(fā)展的QM/MM組合方法［20］，根據(jù)氨基酸模型中氫鍵的幾何結(jié)構(gòu)和鍵能，重新優(yōu)化了QM區(qū)域的范德華相互作用參數(shù)。此時得到的范德華半徑要比OPLS-AA力場大5%～10%，而范德華勢阱深度則沒有改變。增加的范德華排斥用來補(bǔ)償由于MM點(diǎn)電荷引起的QM密度過度極化。最近崔強(qiáng)等［21］提出，在QM/MM組合方法下計(jì)算得到的固相熱力學(xué)數(shù)值，對QM/MM范德華參數(shù)不敏感。

對于成鍵相互作用，需要采用合理的方案來處理，避免出現(xiàn)雙重計(jì)算。一般的規(guī)則是每一個成鍵項(xiàng)要依靠內(nèi)層和外層的原子。

1.4 QM/MM組合方法對邊界原子的處理

對于邊界原子的處理，一般不可避免會碰到共價鍵斷裂的情況。斷鍵的原理一般是不要涉及到鍵的耦合項(xiàng)，最好是極性的，沒有共軛相互作用。處理的方法具體可以分為3類。連接原子方法(linkatom schemes)［18］，邊界原子方法(boundary-atom schemes)［22-25］和定域軌道方法(localized-orbital schemes)［26-28］。

連接原子方法引進(jìn)了額外的原子中心(通常為H原子)，而這并不是真實(shí)系統(tǒng)的一部分。它的引入，增加了人為的自由度，使得結(jié)構(gòu)優(yōu)化過程更加復(fù)雜。雖然存在缺點(diǎn)，但是此種方法仍舊是最流行最廣泛應(yīng)用的邊界原子處理方法。而我們的QM/MM組合方法計(jì)算中采用的也是這種方法。

在邊界原子方法中，MM邊界的一個原子被一個具有兩性的原子所代替，它可以同時出現(xiàn)在QM和MM的計(jì)算中。大多數(shù)提出的邊界原子方法，都是建立在單價贗勢的基礎(chǔ)上，通過參數(shù)化來得到所期望的性質(zhì)。目前常見的方法有adjusted connection atoms［22］，pseudobonds［23-24］和effective group potentials［25］。

定域軌道方法把有方向性的雜化軌道放在邊界原子處，并使其中的一些軌道凍結(jié)，不參與自洽迭代。主要的方法有l(wèi)ocal self-consistent field(LSCF)［26］，frozen orbitals［27］和generalized hybrid orbitals(GHO)［28］。

1.5 自由能計(jì)算

計(jì)算得到的體系自由能可以與實(shí)驗(yàn)得到的數(shù)據(jù)相比較，因此體系自由能的計(jì)算是QM/MM組合方法中的一個重要環(huán)節(jié)。在描述化學(xué)反應(yīng)的過程中，自由能計(jì)算充分考慮了研究體系的漲落(fluctuation)情況，比靜態(tài)的電子結(jié)構(gòu)計(jì)算獲得的相對能量值更具有物理意義。在QM/MM計(jì)算水平下，為了得到體系的自由能曲線，常用的方法一般有兩種:(1)自由能微擾(free energy perturbation，F(xiàn)EP)［29］;(2)熱力學(xué)積分(thermodynamic integration，TDI)［30］。在處理近過渡態(tài)區(qū)域的分子構(gòu)象時，通過正常的非限制性QM/MM MD方法可能找不到高能區(qū)域的構(gòu)象。因此在研究反應(yīng)的過渡態(tài)時，傘形采樣(umbrella sampling，US)方法［31-32］被廣泛用于高能區(qū)域采樣。一般WHAM分析方法［33］可以與傘形采樣相結(jié)合，最終計(jì)算得到反應(yīng)的自由能曲線。

1.6 QM/MM組合方法的計(jì)算方法架構(gòu)

QM/MM組合方法的使用具有強(qiáng)大的柔性，計(jì)算中能夠采用各種各樣的QM和MM計(jì)算方法。在實(shí)際應(yīng)用中，許多生物大分子的QM/MM組合計(jì)算都采用半經(jīng)驗(yàn)作為QM區(qū)域的計(jì)算方法，用DFT的情況相對較少。對于MM部分，現(xiàn)在采用的都是建立在點(diǎn)電荷模型之上的價力場，到目前為止，還沒有合適的極化力場可以廣泛運(yùn)用。普遍運(yùn)用的生物大分子力場主要有CHARMM［14-15］，AMBER［16-17］，GROMOS［34-35］，OPLS-AA［36-38］;而普通的力場則有MM3［39-40］，MM4［41-42］，MMFF［43-44］和UFF［45］。

目前，流行的QM/MM組合計(jì)算的軟件包有3種:(1)在MM的軟件包中加入QM部分的計(jì)算，常見的軟件包有AMBER和CHARMM;(2)在QM的軟件包中加入MM部分的計(jì)算，常見的軟件包有ADF，GAMESS-UK，Gaussian，NWChem，QSite/Jaguar，Car-Parrinello MD codes with QM/MM capabilities，CPMD，CP-PAW;(3)耦合已經(jīng)存在的QM和MM的程序包，常見的軟件包有ChemShell，QMMM和Q-Chem/Tinker。

2 QM/MM組合方法研究進(jìn)展

自從1976年Warshel和Levitt提出雜化QM/MM的概念［2］，并將其用于研究溶解酵素的反應(yīng)機(jī)理以來，許多課題組相繼提出了不同的QM/MM組合方法，已廣泛應(yīng)用于各類相關(guān)的化學(xué)、生物學(xué)和材料等問題。Kollman小組［46］開發(fā)了AMBER軟件，用來計(jì)算神經(jīng)氨酸苷酶等生物大分子體系。Karplus小組［47］提出了AM1/CHARMM的組合方法思想，應(yīng)用于DNA糖基化酶等生物大分子體系的計(jì)算研究。楊偉濤小組［48］研究了QM/MM組合方法中自由能和靜電勢的計(jì)算，QM與MM邊界的處理等問題，并為組合方法的發(fā)展做出了重要的貢獻(xiàn)。高加力和莫亦榮小組［49］建立了塊定域波函數(shù)方法(BLW)和廣義雜化軌道方法(GHO)，并把它們廣泛應(yīng)用于QM/MM組合方法計(jì)算。Thiel小組［50］發(fā)展了半經(jīng)驗(yàn)方法，并利用MNDO/MM MD結(jié)合熱力學(xué)積分的方法來計(jì)算反應(yīng)自由能。張?jiān)鲚x［51-53］采用極化力場方法，研究了大量蛋白質(zhì)分子的光譜性質(zhì)。曹澤星［54-55］研究了磷酸葡萄糖異構(gòu)酶和鼠李糖異構(gòu)酶的催化機(jī)理，闡明了兩種水解酶在Zn配位結(jié)構(gòu)上的差異性。徐定國［56］研究了透明質(zhì)酸酯裂解酶的催化反應(yīng)，提出了順式消除反應(yīng)的機(jī)理。王永［57］研究了細(xì)胞色素P450蛋白酶活化C-H鍵的反應(yīng)機(jī)理，結(jié)果表明電子轉(zhuǎn)移在催化反應(yīng)中起著重要的作用。馬晶［58］、劉成卜［59-60］等采用QM/MM組合方法，在酶催化反應(yīng)機(jī)理領(lǐng)域也做出了重要貢獻(xiàn)。

3 結(jié)論

美國化學(xué)學(xué)會會長Marinda Li Wu稱今年的諾貝爾獎“令人非常興奮”。她解釋道:“獲獎?wù)咄ㄟ^計(jì)算機(jī)模型，為經(jīng)典實(shí)驗(yàn)科學(xué)與理論科學(xué)的聯(lián)系奠定了基礎(chǔ)。由此得到的見解正在幫助我們開發(fā)新的藥物。比如，他們的成果正在用于決定藥物如何與體內(nèi)蛋白質(zhì)相互作用，從而治療疾病?！盦M/MM組合方法的發(fā)展，并且在研究凝聚態(tài)中的化學(xué)反應(yīng)和生物大分子尤其是酶催化反應(yīng)的機(jī)理等方面有著廣泛的應(yīng)用，為處理龐大而又復(fù)雜的體系提供了一種強(qiáng)有力的工具。

［1］Field M J，Bash P A，Karplus M.J Comput Chem，1990，11:700

［2］Warshel A，Levitt M.J Mol Biol，1976，103:227

［3］Elsasser B，F(xiàn)els G，Weare J H.J Am Chem Soc，2014，136:927

［4］Gotz A W，Clark M A，Walker R C.J Comput Chem，2014，35:95

［5］Pecina A，Lepsik M，Rezac J，et al.J Phys Chem B，2013，117:16096

［6］Thellamurege N M，Si D J，Cui F C，et al.J Comput Chem，2013，34:2816

［7］Nakayama A，Arai G，Yamazaki S，et al.J Chem Phys，2013，139:214304

［8］Golze D，Iannuzzi M，Nguyen M T，et al.J Chem Theory Comput，2013，9:5086

［9］Bakowies D，Thiel W.J Phys Chem，1996，100:10580

［10］Riccardi D，Schaefer P，Cui Q.J Phys Chem B，2005，109:17715

［11］Maseras F，Morokuma K.J Comput Chem，1995，16:1170

［12］Humbel S，Sieber S，Morokuma K.J Chem Phys，1996，105:1959

［13］Vreven T，Morokuma K.Theor Chem Acc，2003，109:125

［14］MacKerell A D Jr，Brooks B，Karplus M，et al.CHARMM:The Energy Function and Its Parameterization with an Overview of the Program∥von Rague Schleyer P.Encyclopedia of Computational Chemistry.Vol 1.Chichester:Wiley，1998:271

［15］http:∥www.charmm.org/

［16］Case D A，Cheatham T E，Darden T，et al.J Comput Chem，2005，26:1668

［17］http:∥amber.scripps.edu/

［18］Singh U C，Kollman P A.J Comput Chem，1986，7:718

［19］Bakowies D，Thiel W.J Phys Chem，1996，100:10580

［20］Murphy R B，Philipp D M，F(xiàn)riesner R A.J Comput Chem，2000，21:1442

［21］Riccardi D，Li G，Cui Q.J Phys Chem B，2004，108:6467

［22］Antes I，Thiel W.J Phys Chem A，1999，103:9290

［23］Zhang Y，Lee T S，Yang W.J Chem Phys，1999，110:46

［24］Zhang Y.J Chem Phys，2005，122:024114

［25］Alary F，Poteau R，Heully J L，et al.Theor Chem Acc，2000，104:174

［26］Assfeld X，Rivail J L.Chem Phys Lett，1996，263:100

［27］Philipp D M，F(xiàn)riesner R A.J Comput Chem，1999，20:1468

［28］Gao J，Amara P，Alhambra C，et al.J Phys Chem A，1998，102:4714

［29］Zwanzig R W.J Chem Phys，1954，22:1420

［30］Senn H M，Thiel S，Thiel W.J Chem Theory Comput，2005，1:494

［31］Torrie G M，Valleau J P.J Comput Phys，1977，23:187

［32］Bartels C，Karplus M J.J Comput Chem，1997，18:1450

［33］Bouzida D，Swendsen R H，Kollman P A，et al.J Comput Chem，1992，13:1011

［34］Scott W R P，Hünenberger P H，Tironi I G，et al.J Phys Chem A，1999，103:3596

［35］Jorgensen W L，Maxwell D S，Tirado-Rives J.J Am Chem Soc，1996，118:11225

［36］Jorgensen W L.OPLS Force Fields∥von RaguéSchleyer P.Encyclopedia of Computational Chemistry.Vol 3.Chichester:Wiley，1986:1998

［37］Kaminski G A，F(xiàn)riesner R A，Tirado-Rives J，et al.J Phys Chem B，2001，105:6474

［38］Lii J H，Allinger N L.J Am Chem Soc，1989，111:8566

［39］Lii J H，Allinger N L.J Comput Chem，1991，12:186

［40］Allinger N L，Chen K S，Lii J H.J Comput Chem，1996，17:642

［41］Langley C H，Lii J H，Allinger N L.J Comput Chem，2001，22:1396

［42］Halgren T A.J Comput Chem，1996，17:490

［43］Halgren T A.J Comput Chem，1996，17:520

［44］Rappe A K，Casewit C J，Colwell K S，et al.J Am Chem Soc，1992，114:10024

［45］Casewit C J，Colwell K S，Rappe A K.J Am Chem Soc，1992，114:10035

［46］Case D A，Seetin M G，Sagui C，et al.AMBER 10.University of California，San Francisco，2008

［47］Bash P A，F(xiàn)ield M J，Davenport R C，et al.Biochemistry，1991，30:5826

［48］Hao H，Yang W T.Annu Rev Phys Chem，2008，59:573

［49］Mo Y R，Zhang Y，Gao J L.J Am Chem Soc，1999，121:5737

［50］Bakowies D，Thiel W.J Phys Chem，1996，100:10580

［51］Ji C G，Zhang J Z H.J Am Chem Soc，2011，133:17727

［52］Wang X W，He X，Zhang J Z H.J Phys Chem A，2013，117:6015

［53］Xu J，Zhang J Z H，Xiang Y.J Phys Chem A，2013，117:6373

［54］Wu R B，Xie H J，Cao Z X，et al.J Am Chem Soc，2008，130:7022

［55］Wu R B，Gong W J，Ting L，et al.J Phys Chem B，2012，117:1984

［56］Zheng M，Xu D G.J Phys Chem B，2013，117:10161

［57］Shaik S，Cohen S，Wang Y，et al.Chem Rev，2010，110:949

［58］Jiang N，Ma J.J Phys Chem B，2010，114:11241

［59］Wang J H，Hou Q Q，Dong L H，et al.J Mol Graph Model，2011，30:148

［60］Sheng X，Gao J，Liu Y J，et al.J Mol Graph Model，2013，44:17