盧 天 陳飛武

(北京科技大學化學與生物工程學院,北京100083)

原子電荷計算方法的對比

盧 天 陳飛武*

(北京科技大學化學與生物工程學院,北京100083)

原子電荷是對化學體系中電荷分布最簡單、最直觀的描述形式之一,在理論和實際應用中都有重要意義.本文介紹了12種重要的原子電荷計算方法的原理和特點,通過大量實例從不同角度比較了它們的優(yōu)缺點.這些方法包括Mulliken、分子環(huán)境中的原子軌道(AOIM)、Hirshfeld、原子偶極矩校正的Hirshfeld布居(ADCH)、自然布居分析(NPA)、Merz-Kollmann(MK)、分子中的原子(AIM)、Merck分子力場94(MMFF94)、AM1-BCC、Gasteiger、電荷模型2(CM2)以及電荷均衡(QEq)方法.最后本文對如何在實際應用中選擇合適的計算方法給出了建議.

原子電荷;計算化學;布居分析;電負性;靜電勢

1 引言

原子電荷,即位于原子中心的點電荷,是對化學體系中電荷分布最簡單、最直觀的描述方式之一.它有很多重要意義,比如幫助化學工作者研究原子在各種化學環(huán)境中的狀態(tài)、1考察分子性質(zhì)、2,3預測反應位點.4另外原子電荷模型在計算化學中也有很多實用價值,如作為分子描述符用于藥物虛擬篩選、5在分子對接和分子動力學/蒙特卡羅模擬中描述靜電作用、6在量子力學(QM)與分子力學(MM)結(jié)合計算中表現(xiàn)QM區(qū)域原子對MM區(qū)域原子的靜電作用勢.7在其它一些理論方法中也需要借助原子電荷,如Pipek-Mezey軌道定域化方法、8電荷自洽的緊束縛密度泛函方法(SCC-DFTB)、9Truhlar課題組開發(fā)的一系列溶劑模型(SM)10-13等.

原子電荷并不是可觀測量,也沒有客觀、唯一的定義,因此獲得原子電荷的方法多種多樣.通過一些實驗數(shù)據(jù)可以間接地考察原子電荷,14如通過分子的多極矩、紅外光譜強度和頻率、配體場分裂能、核磁共振(NMR)位移等,但是數(shù)據(jù)的獲取相對不便,與原子電荷的對應關(guān)系存在較大經(jīng)驗性,也不能分析不穩(wěn)定的體系和狀態(tài).計算化學的發(fā)展使原子電荷方便、快速、可靠的獲得成為了可能.15-18自從1955年Mulliken電荷19-21被提出以來,迄今已有不下60種原子電荷計算方法被相繼提出,近年仍有許多研究者在改進計算方法.

雖然在許多計算方法的原文和部分文獻22-25中都對不同類型原子電荷進行了比較和討論,但是涉及的方法并不全面,測試分子和測試內(nèi)容缺乏統(tǒng)一性.因此,很有必要對較為重要的、文獻中常涉及的原子電荷計算方法進行綜合的比較,使它們的特點、優(yōu)缺點能夠清晰地展現(xiàn),這將有助于在實際問題研究中選擇適合的計算方法.我們選取Mulliken、分子環(huán)境中的原子軌道(AOIM)、26Hirshfeld、27原子偶極矩校正的Hirshfeld布居(ADCH)、28自然布居分析(NPA)、29Merz-Kollmann(MK)、30分子中的原子(AIM)、31Merck分子力場94(MMFF94)、32,33AM1-BCC、34,35Gasteiger、36電荷模型2(CM2)37和電荷均衡(QEq)38共12種方法.在文中將首先介紹這些方法的基本原理和特點,之后通過實際體系測試和比較它們的各方面性能,最后討論計算量的大小,并給出方法選擇上的建議.

2 各種原子電荷計算方法的原理和特點

2.1 Mulliken方法

Mulliken方法19-21是最古老的原子電荷計算方法.它的算法簡單,計算量可忽略不計,幾乎出現(xiàn)在所有量子化學軟件中.首先考慮分子軌道波函數(shù)歸一化條件(假設(shè)為實函數(shù),后同)

r代表空間坐標.將分子軌道?i按以原子為中心的基函數(shù)χm展開

這里C代表系數(shù)矩陣.將式(2)代入式(1),積分后得

其中Sm,n=∫χm(r)χn(r)dr.上式中第一項可視為各個基函數(shù)對軌道獨立貢獻之和,稱定域項;第二項是交叉項,體現(xiàn)了由于每一對基函數(shù)之間的耦合對軌道產(chǎn)生的聯(lián)合貢獻.Mulliken將分子軌道i中基函數(shù)m所占成分定義為

也就是說,定域項被完全劃歸到相應基函數(shù),而交叉項被平分給對應的兩個基函數(shù).將所有軌道中屬于相同原子的基函數(shù)的占據(jù)數(shù)加和就得到了原子占據(jù)數(shù),并直接得到原子電荷

其中Z是原子核電荷數(shù),η是軌道占據(jù)數(shù).

Mulliken電荷存在一些嚴重問題:(1)交叉項平分的做法物理意義不嚴格,沒有考慮到原子間的差異性;(2)基組依賴性非常大.尤其是使用大基組、含彌散函數(shù)基組的情況下問題更為明顯.這主要是因為彌散函數(shù)延展范圍大,容易破壞基組的平衡性,后有一些研究者提出了不同方法試圖改進交叉項劃分的不合理性,39-43其中也包括使用相對較多的L?wdin方法,39但是嚴重的基組依賴性卻并沒消除,也并未普及開來.

2.2 NPA方法

自然布居分析(NPA)29的關(guān)鍵是將帶有隨意性的原始基組(一般為擴展基)描述的波函數(shù)轉(zhuǎn)化到物理意義清晰的正交極小基下描述,使基函數(shù)與原子軌道有明確對應關(guān)系,很大程度避免了基組不平衡性問題對結(jié)果的影響,也避免了劃分交叉項的困難.

NPA的原理如下.首先根據(jù)原始基函數(shù)與原子的對應關(guān)系將單粒子密度矩陣中對應于各個原子的對角塊依次提取并對角化,所得本征向量就是初自然原子軌道(PNAO),其本征值就是PNAO的占據(jù)數(shù).根據(jù)占據(jù)數(shù)可將PNAO分為兩類:(1)極小集軌道.它們擁有較高占據(jù)數(shù),是對原子上的電子密度最主要、最緊湊的描述,與基態(tài)原子的內(nèi)層和價層原子軌道一一對應.(2)里德堡集軌道.這是指PNAO中除極小集軌道之外的、電子占據(jù)數(shù)很少的軌道,它們對原子的電子密度的描述不起主要作用.隨后對這些PNAO進行占據(jù)數(shù)權(quán)重的對稱正交化(OWSO),使它們不僅在原子內(nèi)正交也在原子間正交.相對于其它正交化方法,OWSO方法可以使有意義的極小集軌道變形盡量小,而允許意義不大的里德堡集軌道自由變化,物理意義更明確.PNAO在正交化后被稱為自然原子軌道(NAO).以NAO為基重新構(gòu)建體系單粒子密度矩陣,將對應于相應原子的矩陣對角元加和就得到了原子占據(jù)數(shù).

值得一提的是,由于NPA是自然鍵軌道(NBO)方法44框架內(nèi)的一個組成部分,而且NPA電荷通常由Weinhold等開發(fā)的NBO程序45計算,所以NPA電荷經(jīng)常在文獻中被稱為NBO電荷.

2.3 AOIM方法

和NPA的出發(fā)點相近,黎樂民等26,46提出的分子環(huán)境中的原子軌道(AOIM)方法也是將原始擴展基下的波函數(shù)轉(zhuǎn)化為極小基描述后再做布居分析. AOIM有兩種不同的具體實現(xiàn)方法:

第一種方法46認為分子環(huán)境會使得原子軌道相對于原子在自由狀態(tài)時收縮或擴張,分子環(huán)境的影響通過平均化的球?qū)ΨQ有效外勢來表達其中r,θ和φ是球坐標,球坐標系以原子A的原子核為中心,V是分子勢場.在這樣的有效勢下解原子的徑向薛定諤方程,所得波函數(shù)結(jié)合角度波函數(shù)就是分子環(huán)境下的原子軌道.實際求解時使用變分法,將計算分子時用的擴展基的徑向函數(shù)作為變分函數(shù),就得到了原先擴展基與新的極小基之間的變換關(guān)系.

第二種方法26利用Sanchez-Portal投影方法47,48將擴展基波函數(shù)投影到Slater極小基軌道(STO).設(shè)L?wdin正交化后的擴展基為φortho,系數(shù)矩陣為C?;所期望的極小基在L?wdin正交化后為χortho,系數(shù)矩陣為C?(mini),則變換關(guān)系為基函數(shù)的約化會導致基組完備性下降,從而使波函數(shù)信息丟失,衡量丟失量的參數(shù)為.它是STO的指數(shù)與方向的函數(shù).為了最大程度避免因約化帶來的基組不完備性增加,通過有限內(nèi)存大尺度限制性優(yōu)化(L-BFGS-B)方法最小化Δ就得到在分子環(huán)境下膨脹、收縮和旋轉(zhuǎn)后的原子軌道.通過上述過程得到極小基下的系數(shù)矩陣后,就可以照常按照Mulliken方法計算原子電荷.下文中涉及的AOIM方法都特指第二種AOIM方法.

2.4 AIM方法

Bader的分子中的原子(AIM)方法31,49是典型的基于實空間劃分的計算原子電荷的方法.AIM方法將電子密度零通量面定義為原子間的分界面,劃分出的每個原子獨立的空間被稱為原子盆.在分界面上沒有電子密度梯度線穿過,即滿足Δρ(r?)·n(r?)=0,這里r?為原子界面上的任意點,n為界面上的單位法矢量,ρ是電子密度函數(shù).這樣的空間劃分從量子力學角度來看理論意義明確,在每個原子盆內(nèi)維里定理得到滿足.對原子盆內(nèi)電子密度積分并與核電荷求差值即得到原子電荷這里ΩA代表A原子盆.注意以AIM的劃分,在特殊條件下可能會出現(xiàn)不含原子核的原子盆,被稱為贗原子.例如鋰金屬的兩個鋰原子間就存在贗原子,這是由金屬鍵所導致的.另外還可能出現(xiàn)一個盆內(nèi)包含不止一個原子核的情況,如KrH+體系中Kr由于電子分布范圍過廣而淹沒了氫.這些情況都無法使用AIM方法來計算原子電荷.

2.5 Hirshfeld方法

Hirshfeld電荷寫為

其中

2.6 ADCH方法

Hirshfeld電荷數(shù)值普遍偏小,51而且偶極矩、靜電勢重現(xiàn)性比較差,25我們認為這主要是忽略了原子偶極矩所引起的.原子偶極矩可定義為

其中rA代表A原子核坐標.原子偶極矩是對原子空間內(nèi)電子密度各向異性分布最重要的描述,它對分子可觀測性質(zhì)有直接貢獻.然而式(8)的積分卻相當于把原子空間內(nèi)電子密度的分布球?qū)ΨQ化了,完全掩蓋了原子偶極矩的效應.因此提出了原子偶極矩校正的Hirshfeld布居(ADCH).28在這個方法中首先計算各個原子的Hirshfeld電荷及原子偶極矩,然后將每個原子偶極矩按下式展開成為周圍原子的校正電荷

其中ΔqAB代表展開A原子偶極矩后對B原子產(chǎn)生的校正電荷.把所有原子偶極矩展開成校正電荷后,將校正電荷累加到原始Hirshfeld電荷上就得到了ADCH電荷.為了保證分子凈電荷在校正過程中不變,展開任一原子A的偶極矩時需滿足這個條件

同時,為了讓原子偶極矩只展開到離當前原子較近的原子上,在實際計算時令式(11)和式(12)以拉格朗日乘子方式作為限制條件,使下式最小化來獲得A原子對其它原子產(chǎn)生的校正電荷

其中rAB是A、B原子間距離,N是總原子數(shù),ν是依賴于原子間距離的函數(shù).它呈倒S形狀,當rAB由0變化到A、B原子范德華半徑之和的過程中νAB由1變化為0.因此,在原子偶極矩展開時,離當前原子越近的原子越容易獲得校正電荷,與當前原子無相互作用的原子(距離超過范德華半徑和)的電荷不會受到影響.

2.7 MK方法

在化學體系中,靜電勢(ESP)由原子核電荷和電子密度兩部分貢獻構(gòu)成對于離原子核較近的區(qū)域,總是由第一項主導,靜電勢數(shù)值為正,這是化學上不感興趣的.而體系范德華表面以外的區(qū)域靜電勢可正可負,具有顯著的化學意義,對于研究受體-配體結(jié)合方式、晶體堆疊方式,預測親核/親電反應位點,通過比較場分析(CoMFA)進行藥物篩選等問題至關(guān)重要.若在分子范德華表面外側(cè)一定范圍內(nèi)選取一批格點,通過最小二乘法令原子電荷在這些點產(chǎn)生的靜電勢

對式(14)計算的靜電勢擬合,就得到了擬合靜電勢電荷.擬合靜電勢電荷計算方法很多,主要差異在于格點位置的選取和擬合過程的細節(jié).最常見的方法有Merz-Kollman(MK)方法、30靜電勢獲得的電荷(CHELP)方法52和基于格點的CHELP(CHELPG)方法.53CHELPG電荷在三者中旋轉(zhuǎn)不變性最好,而通過多極矩重現(xiàn)性分析發(fā)現(xiàn)MK電荷質(zhì)量比CHELP和CHELPG略好,54因此本文選取MK電荷作為擬合靜電勢電荷的代表在后文中與其它類型電荷進行比較.

2.8 AM1-BCC方法

傳統(tǒng)擬合靜電勢方法存在一些弊端,例如:(1)電荷對構(gòu)象依賴性較大;(2)單一構(gòu)象擬合的原子電荷不能體現(xiàn)原子的等價性,例如一氯丙烷中甲基的三個氫是化學等價的,但任何構(gòu)象下擬合出來的這三個氫的電荷都不是全同的;(3)內(nèi)部被包埋的原子電荷擬合不準確,這是由于這些原子距離擬合點過遠所導致的.這三個問題使得擬合靜電勢電荷用于分子動力學模擬等問題中可能導致錯誤的結(jié)果.限制性靜電勢擬合(RESP)方法55可以很大程度解決這些弊端,它對內(nèi)部原子在擬合時施加限制勢以降低數(shù)值不穩(wěn)定性,對等價原子在擬合時強制相等,并且將擬合分為兩步或者多步以提高擬合質(zhì)量.然而計算RESP電荷過程相對復雜,而且需要較為準確的量子化學方法計算靜電勢,用于大批量或者較大分子體系較為耗時.

AM1-BCC34,35是一種簡單的近似計算RESP電荷的方法,其電荷是在AM1電荷(即AM1半經(jīng)驗方法56獲得的波函數(shù)產(chǎn)生的Mulliken電荷)基礎(chǔ)上進行鍵電荷校正(BCC)得到的.校正過程只依賴于原子類型和原子間連接關(guān)系.校正參數(shù)擬合自2755個種類多樣的有機分子的HF/6-31G*下的RESP電荷.例如甲醇中的碳原子形成三個“Csp3-single-H”鍵和一個“Csp3-single-Osp3”鍵,根據(jù)查尋事先擬合的BCC參數(shù)就可得知校正電荷應為3×0.0274+1×0.0835= 0.1657(本文電荷單位均為原子單位),加到碳的AM1電荷上就是AM1-BCC電荷.

2.9 MMFF94方法

Merck分子力場94(MMFF94)32,33獲得電荷的方式在形式上與AM1-BCC十分相似,同樣使用鍵電荷校正過程,即考慮到每個原子形成的化學鍵的極性來對初始電荷校正.與AM1-BCC不同之處是MMFF94的初始電荷由原子類型直接確定,且參數(shù)來自于擬合HF/6-31G*下計算的分子偶極矩及相互作用能.

2.10 CM2方法

電荷模型(CM)系列方法目前包括CM1、57CM2、37,58CM359-61和CM4,62,63它們都是對低級別方法下的Mulliken或L?wdin電荷39進行簡單校正的方法,目的是使原子電荷計算出的偶極矩盡可能精確地重現(xiàn)氣相實驗值或者高精度量子化學計算值.這一類電荷也被用于不同版本的SM系列溶劑模型中.

CM2校正的公式與AM1-BCC和MMFF94的相比略為復雜,CM2電荷可寫為

其中M是原子間Mayer鍵級.64C與D是對實驗氣相偶極矩擬合的參數(shù),它們?nèi)Q于化學鍵的種類,而鍵的種類僅由它相連的兩個原子的元素決定.CM2方法對使用AM1、PM3、HF/MIDI!、HF/6-31G*等級別計算L?wdin電荷的情況分別擬合了校正參數(shù).

2.11 Gasteiger方法

Gasteiger電荷36也被稱為軌道電負性部分均衡(PEOE)電荷,這種方法利用了電負性均衡概念,65即電負性不同的原子成鍵時,電負性較小的原子附近電子密度會流向電負性較大的原子.在這個過程中原先電負性小的原子電負性會增大.當所有原子間電負性相等時,電子密度的分布就是平衡狀態(tài)分布.

PEOE將原子電負性與原子電荷關(guān)系表達為

其中a、b和c為預設(shè)參數(shù),對于相同元素但價層軌道處于不同雜化狀態(tài)的情況參數(shù)不同.參數(shù)通過實驗測定的相應元素的基態(tài)和電離態(tài)在相應雜化狀態(tài)下的電離勢與電子親合勢來計算.PEOE的迭代過程如同電荷不斷在原子間轉(zhuǎn)移的過程,當?shù)諗繒r電荷分布也就平衡了.在每一步中按照下式計算每一對相連原子間電荷轉(zhuǎn)移量

2.12 QEq方法

電荷均衡方法(QEq)38得到的是在電負性完全均衡下的原子電荷.電負性表達式為電負性均衡條件要求χ1=χ2=…=χN,它提供了N-1個條件,再將所有原子電荷總和等于分子凈電荷作為剩余條件,通過解線性方程即可到得全部原子電荷.式(19)中是廣義化的Mulliken電負性為自庫侖積分,是原子電離勢與電子親和勢的差值.JAB描述了原子間靜電作用,若兩原子距離較遠,則JAB= 1/rAB;若原子間距離較近,電子云靜電屏蔽作用不能忽視,此時JAB被定義為兩原子間庫侖積分.每個原子用單個Slater函數(shù)描述電子密度分布,其中指數(shù)通過擬合堿金屬鹵化物實驗偶極矩數(shù)據(jù)確定,對所有原子都相同.普適型力場UFF66使用的正是QEq電荷模型.

值得一提的是QEq方法可以視為另一種使用廣泛的方法電負性均衡方法(EEM)67的變體.然而EEM的經(jīng)驗性較大,不同研究者以不同的目標數(shù)據(jù)擬合了參數(shù),結(jié)果有很大差異,因此不納入本文的比較.

除了上面涉及的方法外,廣義原子極化張量(GAPT)電荷68,69也常出現(xiàn)在文獻當中.然而GAPT電荷計算過于耗時(和振動分析耗時相近),實際應用意義很有限,因此并不在本文進行介紹和比較.讀者可參看文獻25對它的測試和討論.

3 計算細節(jié)和測試體系

下文將對前面介紹的12種方法得到的原子電荷與化學經(jīng)驗相符程度、對偶極矩和靜電勢的重現(xiàn)性以及基組依賴性進行檢驗,最后討論計算耗時以及計算方法的選擇.測試體系以研究得最為普遍的有機小分子為主,晶體以及包含重元素的體系不在本文討論范圍內(nèi).測試分子皆處于穩(wěn)定構(gòu)型,過渡態(tài)、激發(fā)態(tài)、受外場作用等狀態(tài)也不屬于本文涉及范圍.實際上,Gasteiger、CM2等包含經(jīng)驗參數(shù)的方法也不適用于這些特殊問題.本文所測試的方法分為兩組,第一組不依賴于任何參數(shù),包括Mulliken、AOIM、Hirshfeld、ADCH、NPA、MK和AIM,第二組依賴于擬合的參數(shù),包括MMFF94、AM1-BCC、CM2、Gasteiger和QEq.

除了理論方法和基組依賴性部分的討論外,下文的中性分子、陽離子的結(jié)構(gòu)優(yōu)化和波函數(shù)計算均在HF/6-31G*70-72級別下進行,陰離子的結(jié)構(gòu)優(yōu)化和波函數(shù)計算在HF/6-31+G*73級別下進行.波函數(shù)的計算和結(jié)構(gòu)優(yōu)化,以及Mulliken、NPA、MK、CM2、QEq電荷的計算都通過Gaussian 03程序74完成.計算CM2電荷時的初始電荷由HF/6-31G*級別獲得.為了使擬合質(zhì)量較好,計算MK電荷時通過IOp(6/ 42=6)選項將每單位面積的擬合點由默認的1個增加到6個.Hirshfeld和ADCH電荷通過我們開發(fā)的Multiwfn 2.1.2程序75計算.AIM電荷通過AIMALL 10.05.04程序76計算.AOIM電荷由AOIM 1.1程序77計算(目前沒有其它公開的程序可以計算AOIM電荷).MMFF94電荷通過Avogadro 1.0.3程序78獲得. AM1-BCC和Gasteiger電荷由AmberTools 1.5程序79計算.由于計算AM1-BCC電荷前AmberTools會自動在AM1下進行結(jié)構(gòu)優(yōu)化,我們通過修改代碼將此步驟取消.

4 普通小分子的原子電荷

我們首先檢驗不同方法計算的種類各異的小分子電荷,表1列出了第一組方法計算的結(jié)果.容易看出,不同方法得到的原子電荷的數(shù)值大小在整體上有所不同.Hirshfeld電荷整體偏小,許多原子電荷的大小都不足0.1.ADCH是對Hirshfeld方法的校正,校正后原子電荷普遍增大,電荷絕對值加和比校正前增加了一倍.對于碳、氫、氮、氧構(gòu)成的中性有機分子,Mulliken、AOIM、MK和NPA電荷在數(shù)量級上比較接近,對于個別分子它們也存在著明顯差異,例如ClF3中氯的MK電荷僅為Mulliken、AOIM和NPA方法計算的約一半.很多原子的AIM電荷明顯大于其它方法所得電荷,例如丙酮的氧的電荷達到-1.35,氰化氫的氮的電荷高達-1.48.從電荷絕對值加和來看,AIM計算的電荷值(44.19)達到了Hirshfeld計算的電荷值(9.92)的將近5倍.

表1 第一組方法計算的一些小分子的原子電荷Table 1 Atomic charges of some small molecules calculated by the first group of methods

continued Table 1

Mulliken電荷對于CLi4和BeO這樣高離子性的體系明顯偏小,甚至小于Hirshfeld方法.這和交叉項的平分處理有直接關(guān)系,它導致了本應屬于陰離子的電子被過多地劃歸到了陽離子上.而對這樣的體系A(chǔ)OIM電荷則比較大,這是由于在分子環(huán)境中優(yōu)化STO的指數(shù)后,陰離子的STO較為彌散,一定程度侵入到陽離子空間內(nèi),而陽離子的STO較為緊縮,因此在新的基函數(shù)下做Mulliken分析時陰離子能夠獲得更多的電子.

電荷的大小與其合理性并沒有絕對關(guān)系,考察原子電荷的合理性關(guān)鍵之一是看它能否與一般化學經(jīng)驗、化學觀念相吻合.磺酸基是間位定位基,因此苯磺酸的間位碳的原子電荷應當比鄰對位更負,第一組方法中除了AIM以外都正確地給出了預期的結(jié)果,而AIM方法卻認為間位碳的電荷大于對位的.可見AIM電荷對于研究反應位點并不適合.對于偏二甲肼,由于碳的電負性大于氫,因此氨基氮的電荷應當比中間的氮原子的電荷更負,而只有AIM方法給出的結(jié)論是相反的.磷的Pauling電負性比碳略小,因此在CH2PH中磷的Mulliken、AOIM、NPA、Hirshfeld電荷為不很大的正值并不意外.而ADCH和MK電荷更側(cè)重于等效描述電子密度的實際分布,由于磷的孤對電子產(chǎn)生的極化效應,使得磷的ADCH和MK電荷都為不大的負值,這也是合理的.而AIM給出的磷的電荷高達1.51,明顯缺乏合理性.

在HF/6-31G*波函數(shù)級別下計算乙炔的AIM電荷時出現(xiàn)了前文談到的“贗原子”,它位于兩個碳原子之間.由于乙炔的對稱性,我們將贗原子盆內(nèi)電荷積分值平分給兩個碳原子來獲得碳的AIM電荷,但是當體系缺乏這樣的對稱性時就不能這樣計算AIM電荷了.很多含有炔基的分子也都存在這樣的贗原子,這使得AIM方法的適用性受到一定限制.值得一提的是,將基組替換為與6-31G*級別較為相似的SVP基組80后贗原子隨即消失,并變?yōu)樘?碳間的鍵臨界點,這顯示出在特殊情況下AIM的適用性明顯受基組影響.

表2列出了第二組方法計算的小分子的原子電荷,由于缺少參數(shù),其中刪去了表1中的BeO和CLi4.很明顯Gasteiger電荷相對其它電荷明顯整體偏小,而MMFF94和AM1-BCC電荷通常比CM2和QEq電荷稍大.從所列分子上看,AM1-BCC和CM2電荷都與經(jīng)驗觀念基本一致.由于MMFF94電荷計算方法過于簡單,沒有考慮電子效應的影響,所得苯磺酸的鄰、間、對碳原子電荷相同,均為-0.15,是不合理的.Gasteiger和QEq方法計算的苯磺酸的硫原子電荷都不合理,數(shù)值接近0.由于硫原子與三個電負性明顯比它更大的氧原子相連,故硫原子電荷理應為較大的正值,其它電荷計算方法的結(jié)果也驗證了這一點.另外Gasteiger電荷并沒有正確預測出苯磺酸的間位碳具有比對位的碳更負的電荷,對這種體系有必要改用明確考慮共軛效應的Gasteiger電荷改進方法,即PEPE方法.81,82

表2 第二組方法計算的一些小分子的原子電荷Table 2 Atomic charges of some small molecules calculated by the second group of methods

continued Table 2

圖1 第一組方法計算的甲烷的甲基電荷與取代原子電負性的關(guān)系Fig.1 Relationship between methyl charges calculated by the first group of methods and electronegativity of substituent atoms of methane

5 原子電荷與電負性的相關(guān)性

對于由σ鍵構(gòu)成的典型分子,合理的原子電荷計算方法應當能夠與電負性規(guī)則基本吻合.對甲烷以不同電負性原子進行取代,各種方法計算的甲基的電荷(甲基各原子電荷的加和)如圖1和圖2所示. Mulliken、AOIM、NPA、Hirshfeld、ADCH、Gasteiger、CM2和QEq雖然在曲線整體斜率以及曲線形狀上有別,但是都正確顯示出甲基電荷隨著取代基原子電負性增大而逐漸增加.MMFF94和AIM方法計算的甲烷中甲基電荷基本為0,然而從碳、氫的電負性上來看此時甲基電荷理應為負值.MK方法對于氮、氧、氟取代甲烷體系給出的結(jié)果與電負性關(guān)系不相符,隨著取代原子電負性的增加甲基的MK電荷不增反降.出現(xiàn)這種情況主要原因在于在這三種甲烷取代物中氮、氧、氟的電子極化程度依次減小,如果用原子偶極矩,即式(10)來表示極化程度,則數(shù)值分別為0.296、0.237和0.143 a.u..原子核附近電子極化對分子外側(cè)靜電勢產(chǎn)生的影響會等價地體現(xiàn)在擬合靜電勢方法獲得的電荷中,并且由于在這三種體系中電子極化效應對MK電荷的影響超過了電負性差異的影響,因此MK電荷表現(xiàn)出與電負性規(guī)則相反的變化趨勢.由于AM1-BCC是對擬合靜電勢電荷的近似,故它對于氮、氧取代甲烷體系給出的結(jié)果與電負性關(guān)系相違背也是預料之中的.由于ADCH對Hirshfeld電荷的校正過程中使電子密度極化效應得以等效體現(xiàn),與擬合靜電勢方法有一定共通之處,這導致對于氮、氧、氟取代時ADCH的甲基電荷的變化曲線的斜率略小于Hirshfeld的斜率.然而ADCH并沒有像MK電荷那樣違背電負性規(guī)律,因此在這一點上比MK具有更強的合理性.

圖2 第二組方法計算的甲烷的甲基電荷與取代原子電負性的關(guān)系Fig.2 Relationship between methyl charges calculated by the second group of methods and electronegativity of substituent atoms of methane

對甲烷以不同數(shù)目的氟進行取代時每個氟原子的電荷應有所不同.1至4取代時每個氟原子所連基團分別相當于―CH3、―CH2F、―CHF2和―CF3,基團電負性依次增加,因此合理的方法計算出的氟原子電荷應當依次減小.由圖3、圖4可見,Mulliken、Hirshfeld、ADCH、CM2、Gasteiger、QEq方法計算的結(jié)果都完全符合這個趨勢,而且氟原子電荷的減小是接近線性的.而其它幾種方法則表現(xiàn)出不同程度的不合理性.AOIM電荷與期望的趨勢不符,1取代時氟原子的電荷值比2、3、4取代時的都要小,這說明AOIM在理論上,或者計算程序的數(shù)值算法上有待修正.從1取代變?yōu)?取代時,MK電荷僅減小了0.09,而NPA電荷則基本沒有變化.AIM電荷在1至4取代過程中幾乎都沒有發(fā)生變化,甚至1取代變?yōu)?取代過程中電荷甚至還增大了0.002(可能由于數(shù)值積分誤差導致).MMFF94的氟原子電荷始終為-0.340,這是因為MMFF94原子電荷只能考慮相鄰原子對它的影響,很明顯,這種過于簡單的計算電荷的方式并不能準確描述原子在分子中的狀態(tài). AM1-BCC也未能很好表現(xiàn)氟原子電荷應有的變化趨勢,雖然隨取代數(shù)目增加氟原子電荷確實依次減小,但是曲線卻過于平坦.

圖3 第一組方法計算的不同取代數(shù)目的氟代甲烷的氟原子電荷Fig.3 Atomic charge of fluorine atom calculated by the first group of methods of fluomethane with different numbers of substitutionsTheAIM charge of fluorine atom is too large to be shown on the graph,the values are-0.742,-0.744,-0.744,and-0.737 respectively for 1 to 4 substitutions.

圖4 第二組方法計算的不同取代數(shù)目的氟代甲烷的氟原子電荷Fig.4 Atomic charge of fluorine atom calculated by the second group of methods for fluomethane with different numbers of substitutions

6 偶極矩的重現(xiàn)性

偶極矩是小分子電荷分布的最簡單也是最重要的描述,分子偶極矩重現(xiàn)性的好壞是判斷原子電荷合理性重要的客觀標準.我們選取了20種具有結(jié)構(gòu)特征不同、有代表性,并且已知實驗偶極矩的小分子,對各種方法所得電荷的偶極矩重現(xiàn)性進行了測試,結(jié)果列于表3和表4.注意對于第一類電荷計算方法,只有將原子電荷計算的偶極矩與計算電荷時所用的波函數(shù)的偶極矩期望值相互比較才是有意義的,直接與實驗值相比較是沒有意義的,表4顯示當前所用的HF/6-31G*波函數(shù)偶極矩期望值普遍大于實驗值(平均大0.26 Debye(1 Debye=3.338× 10-30C·m)),這主要是由于忽略了相關(guān)作用導致的. ADCH是偶極矩保守的方法,因此ADCH電荷計算的偶極矩與偶極矩期望值完全相符.由于靜電勢、偶極矩都是分子實際電荷分布的客觀反映,因此以重現(xiàn)靜電勢為目的的MK電荷也顯示出很好的偶極矩重現(xiàn)性.有人提出在擬合靜電勢的過程中將偶極矩能夠精確重現(xiàn)作為限制條件,85但實際上這么做已沒有必要.Mulliken方法的偶極矩重現(xiàn)性從統(tǒng)計結(jié)果上看并不理想,對個別分子其偶極矩誤差很大,例如噻吩偶極矩期望值為0.90 Debye,而Mulliken方法卻認為這是一個幾乎無極性的分子.需要指出的是Mulliken電荷的基組依賴性很大,因此Mulliken電荷的偶極矩重現(xiàn)性也可能會隨基組的不同而有一定變化.AOIM的偶極矩重現(xiàn)性比Mulliken稍有提高.Hirshfeld方法計算的偶極矩從平均無符號誤差(MUE)上看比Mulliken方法更差,而平均含符號誤差(MSE)顯示Hirshfeld的偶極矩是顯著偏低的,這與Hirshfeld電荷偏小有直接關(guān)系.NPA電荷的偶極矩重現(xiàn)性很不理想,有高估偶極矩的傾向,尤其是對PF3的偶極矩高估了4.2倍,而對二甲硫醚的偶極矩卻低估了10倍以上.AIM電荷對分子偶極矩基本沒有重現(xiàn)能力,MUE高達3.00 Debye,這個值已經(jīng)大于大多數(shù)小分子偶極矩.AIM的偶極矩誤差主要來自于高估,這很大程度上是由于AIM電荷數(shù)值過大所導致的.

表3 第一組方法計算的20種小分子的分子偶極矩(單位為Debye)以及與HF/6-31G＊級別下偶極矩算符期望值(〈μ〉)的比較Table 3 Molecular dipole moment(unit in Debye)of twenty small molecules calculated by the first group of methods,compared to the expected value of dipole moment(〈μ〉)operator at the HF/6-31G＊level

在擬合參數(shù)過程中,MMFF94的參數(shù)主要來自于使MMFF94電荷產(chǎn)生的偶極矩逼近HF/6-31G*密度產(chǎn)生的分子偶極矩,AM1-BCC的參數(shù)來自于使AM1電荷在校正后逼近HF/6-31G*級別下計算的RESP電荷.它們都使用HF/6-31G*下的數(shù)據(jù)作為目標數(shù)據(jù)是因為這個波函數(shù)級別對偶極矩的高估被認為可以等效反映出分子在凝聚相中存在的極化效應,因此MMFF94和AM1-BCC的偶極矩應當與HF/6-31G*下的偶極矩期望值相比較.從表4可見這兩種電荷的偶極矩重現(xiàn)性對大部分分子基本可以令人滿意,MMFF94的MUE為0.36 Debye,雖然AM1-BCC和MMFF94在鍵電荷校正過程中很相似,但AM1-BCC的MUE更小,為0.25 Debye,其改善在一定程度上可能是由于AM1-BCC使用了更有意義的初始電荷.Gasteiger的偶極矩重現(xiàn)性無論以氣相實驗值為參考還是以HF/6-31G*偶極矩為參考都顯得一般,并且對多數(shù)分子有低估偶極矩的傾向.CM2電荷能夠很好地重現(xiàn)氣相實驗偶極矩, MUE僅為0.15 Debye,對所有分子都沒有出現(xiàn)過大誤差,這與CM2提出的目標完全相一致.QEq對實驗偶極矩和HF/6-31G*下的偶極矩重現(xiàn)性都很不理想,對許多分子偶極矩高估了近一倍甚至一倍以上,嚴重低估偶極矩的情況同樣存在,如2-甲基嘧啶.

表4 第二組方法計算的20種小分子的分子偶極矩(單位為Debye)以及與HF/6-31G＊級別下偶極矩算符期望值和氣相實驗值的比較Table 4 Molecular dipole moment(unit in Debye)of twenty small molecules calculated by the second group of methods, compared to the expected value of dipole moment operator at the HF/6-31G＊level and gas-phase experimental value

綜上所述,Mulliken、AOIM、Hirshfeld、NPA、Gasteiger、QEq電荷,尤其是AIM電荷,對偶極矩重現(xiàn)性都不好.ADCH和MK都能很準確重現(xiàn)波函數(shù)的偶極矩期望值,若波函數(shù)完全精確,則這兩種電荷計算的偶極矩將與氣相的實際值相一致,而利用CM2方法則可以在較低波函數(shù)級別下很準確地預測氣相偶極矩.MMFF94和AM1-BCC能夠較好地重現(xiàn)HF/6-31G*的偶極矩.

7 靜電勢的重現(xiàn)性

當研究的問題涉及分子間較近距離的靜電相互作用時偶極模型顯得過于簡單,只有分子范德華表面外側(cè)精確靜電勢也能被較好重現(xiàn)時,才說明這種原子電荷適合研究分子間近距離靜電作用.又由于靜電勢是可觀測量,因此靜電勢的重現(xiàn)性與偶極矩一樣也是檢驗原子電荷合理性的重要、客觀的指標.其重現(xiàn)性常以原子電荷在分子范德華表面外側(cè)所取的格點位置上產(chǎn)生的靜電勢VModel與波函數(shù)計算的精確靜電勢VReal之間的相對均方根偏差(RRMSE)來反映,計算公式為

我們選取了21種中性分子和4種離子用于測試不同原子電荷計算方法的靜電勢重現(xiàn)性,見表5和表6.其中既包含一部分測試偶極矩重現(xiàn)性所用的小分子,也包含結(jié)構(gòu)更復雜的生物、藥物分子.計算RRMSE所用格點設(shè)定與計算MK電荷時的一致,參考靜電勢在HF/6-31G*級別下獲得.

由于MK電荷正是通過最小化靜電勢重現(xiàn)性誤差得到的,因此它的RRMSE擁有理論最低值.從RRMSE的平均值上看,其余11種方法的靜電勢重現(xiàn)能力可以大致分為四個檔次:(1)重現(xiàn)性很好,平均RRMSE在0.25上下.這一類包括ADCH、MMFF94、AM1-BCC和CM2.如果改用能夠比較準確地重現(xiàn)氣相偶極矩的波函數(shù),即使用B3LYP/ cc-pVTZ波函數(shù)代替HF/6-31G*去計算參考靜電勢,則CM2的平均RRMSE值會從0.268降低到0.228.(2)重現(xiàn)性不錯.AOIM屬于此類,平均RRMSE為0.337.(3)重現(xiàn)性一般,RRMSE平均值在0.4-0.6.包括Mulliken、Hirshfeld、NPA、Gasteiger和QEq.(4)重現(xiàn)性很差,這一類只包括AIM,平均RRMSE高達1.53.

表5 第一組方法產(chǎn)生的靜電勢相對于HF/6-31G＊級別下精確靜電勢的相對均方根偏差Table 5 Relative root mean square error(RRMSE)of the ESPproduced by the first group of methods relative to the exact ESPat the HF/6-31G＊level

通過檢驗MK電荷的數(shù)據(jù),可以看出原子電荷模型本身對于一些體系存在明顯局限性.對于丙烷和CH3PH2分子,即便是MK電荷,其RRMSE也分別高達0.72和0.51,因此我們未將這兩個分子納入RRMSE平均值的統(tǒng)計中.若想更準確地描述靜電相互作用,即降低RRMSE,就必須增加一些非原子中心的點電荷,或者在原子中心引入點多極矩.25若在計算MK電荷時在每個原子中心都引入一個點偶極矩(點偶極矩向量與原子電荷一起被擬合),則丙烷和CH3PH2分子的RRMSE會分別降至0.22和 0.13,靜電勢重現(xiàn)性大大提高.而對于一些體系用原子電荷模型表現(xiàn)靜電勢是足夠準確的,如MK電荷對Aspirin的RRMSE僅有0.043.

從表中可見所有方法對離子的靜電勢重現(xiàn)性都顯著好于中性分子,RRMSE降低了約一個數(shù)量級.出現(xiàn)這種情況是由于帶電分子外圍靜電勢主要由分子單極矩所主導,而所測試的所有原子電荷計算方法都是單極矩保守的,即原子電荷加和等于分子凈電荷.不同方法對帶電分子靜電勢重現(xiàn)性的優(yōu)劣關(guān)系和中性分子情況是相似的.

8 波函數(shù)質(zhì)量的依賴性

圖5比較了第一類方法計算的乙酸的甲基碳原子電荷的基組依賴性,從左到右基組完備性逐漸增加.從STO-3G86提升至3-21G87的過程中各種方法計算的結(jié)果都有不同程度的變化,通常被認為具有很好基組穩(wěn)定性的NPA電荷的變化值甚至大于公認基組穩(wěn)定性差的Mulliken電荷.注意由于極小基的完備性很低,所以提升至分裂價基后結(jié)果的明顯變化并不能說明方法的合理性差.從6-311G**88開始除了Mulliken方法外所有方法的結(jié)果都已經(jīng)收斂,提升至aug-cc-pVTZ89,90的過程中結(jié)果變化都不顯著.其中Hirshfeld和ADCH電荷收斂得最早,從3-21G開始就基本不發(fā)生變化.同時也說明ADCH對Hirshfeld電荷的校正并沒有破壞Hirshfeld電荷良好的基組穩(wěn)定性.MK和AOIM從6-31G*開始收斂,而NPA和AIM收斂得相對稍遲.唯一不能隨基組增加而收斂的方法是Mulliken,曲線在基組增大過程中始終強烈波動,對6-311G**基組下重原子增加彌散函數(shù)成為6-311+G**73后,其它方法得到的甲基碳電荷變化均不超過0.01,而Mulliken電荷數(shù)值則降低了0.15,這是因為碳的彌散函數(shù)嚴重侵入到了氫的原子空間內(nèi),原本在6-311G**基組下屬于氫的電子密度現(xiàn)在有一部分被碳的基函數(shù)所描述,因此氫的電子占據(jù)數(shù)一部分被歸屬到了碳.而進一步提升到6-311++G(2df,2p)73后,由于氫的彌散函數(shù)也明顯侵入到了碳的原子空間內(nèi),碳的電子占據(jù)數(shù)有不少被劃歸給了氫,因此碳的Mulliken電荷又增加了許多.6-311++G(2df,2p)已經(jīng)是完備性很高的基組,繼續(xù)提升質(zhì)量后原子電荷應當不出現(xiàn)明顯改變,然而提升至aug-cc-pVTZ后Mulliken電荷變化卻仍高達0.43.因此從基組依賴性的角度上講, Mulliken方法是存在嚴重缺陷的,沒有辦法“唯一”地獲得Mulliken電荷,基組的輕微改變就可能使結(jié)果有定性的變化.同時也說明提升基組質(zhì)量并不會提升Mulliken電荷的合理性,用大基組計算Mulliken電荷是沒有意義的.實際上大基組,尤其是包含彌散函數(shù)的情況,基函數(shù)的化學意義往往比較差,與原子軌道缺乏對應性,而極小基則能夠與原子軌道相互對應,會使Mulliken方法的計算過程更有物理意義.以上的討論對于乙酸的其它原子的電荷也同樣適用.

表6 第二組方法產(chǎn)生的靜電勢相對于HF/6-31G＊級別下精確靜電勢的相對均方根偏差Table 6 RRMSE of the ESPproduced by the second group of methods relative to the exact ESPat the HF/6-31G＊level

圖5 不同方法結(jié)合不同基組在Hartree-Fock級別下獲得的乙酸的甲基碳原子電荷Fig.5 Atomic charge of methyl carbon of acetic acid calculated by different methods under various basis-sets at the Hartree-Fock level

第一類方法對計算波函數(shù)所用理論方法的依賴性如圖6所示,其中包括不含相關(guān)效應的Hartree-Fock方法、廣義梯度近似(GGA)泛函BP86、91,92雜化泛函B3LYP,93以及MP2和CCSD70方法,基組均使用cc-pVDZ.89由于相關(guān)效應對多重鍵影響較大,所以這里選取乙酸的羰基氧原子用于測試.由于AOIM 1.1程序只能讀取SCF波函數(shù),因此MP2和CCSD的結(jié)果未在圖中顯示,但實際上利用自然軌道形式不難將AOIM的應用范圍擴展到后HF波函數(shù)上.從圖6可見對于任何一種計算原子電荷的方法,無論通過何種方式引入相關(guān)效應,所得結(jié)果都很相近,數(shù)值差異最大不超過0.05.而Hartree-Fock波函數(shù)和各種相關(guān)波函數(shù)下的計算結(jié)果卻差異明顯,引入相關(guān)效應后所有方法給出的羰基氧原子電荷均明顯減小,這反映出相關(guān)效應會減小化學鍵的極性這一普遍現(xiàn)象.通過檢驗乙酸各個原子的電荷會發(fā)現(xiàn)Hirshfeld和ADCH電荷受電子相關(guān)效應影響相對于其它方法更弱,比如BP86下Hirshfeld和ADCH的羰基氧原子電荷相對于Hartree-Fock下的結(jié)果差異分別為0.07、0.09,而Mulliken、AOIM、NPA、MK和AIM電荷變化分別為0.16、0.13、0.13、0.12和0.18.

圖6 不同方法結(jié)合cc-pVDZ基組和不同哈密頓得到的乙酸的羰基氧原子電荷Fig.6 Atomic charge of carbonyl oxygen of acetic acid calculated by different methods in combination with cc-pVDZ basis-set and various Hamiltonians

由上可見,至少對于普通有機分子,獲得原子電荷并不需要高級別后HF方法計算波函數(shù),選擇適當?shù)拿芏确汉椒ㄓ嬎愕牟ê瘮?shù)就可以滿足要求.對于NPA、MK、AIM,尤其是Hirshfeld和ADCH方法,只需要中等質(zhì)量的基組就已經(jīng)足夠.在此之上繼續(xù)提高波函數(shù)質(zhì)量意義不大.

9 計算的耗時

計算方法的耗時大小直接關(guān)系到它的適用領(lǐng)域.計算耗時不僅取決于方法本身的復雜度,也明顯依賴于計算程序的代碼效率、數(shù)值方法和編譯、運行時的軟硬件環(huán)境.另外,很多方法的可變參數(shù)也會影響耗時,比如AIM、Hirshfeld方法中對積分的期望精度、AOIM的L-BFGS-B步驟中的收斂閾值、靜電勢擬合時擬合點的密度等.本節(jié)對各種方法計算耗時的討論只是定性的,計算耗時以當前普通實驗室的計算能力作為參照.

MMFF94、Gasteiger和QEq方法由于不依賴于波函數(shù),方法本身也比較簡單,所以計算這三種電荷總耗時極小,即使對于幾千個原子的體系也可以在數(shù)秒內(nèi)給出結(jié)果,它們也能夠用于需要高效計算大批量小分子電荷的情況,比如藥物高通量篩選.由于QEq計算量小,而且是依賴于幾何結(jié)構(gòu)的,能反映出原子電荷對周圍環(huán)境變化的響應,因此QEq可以作為浮動電荷力場94的電荷模型,即在動力學/蒙特卡羅模擬過程中每隔一定步數(shù)或構(gòu)象改變一定程度后自動更新原子電荷.

AM1-BCC電荷的計算總耗時完全來自于AM1單點計算,BCC校正過程幾乎不需要耗時.對于幾十個原子體系的AM1計算,通常在數(shù)秒內(nèi)能夠完成.對于數(shù)百個原子的體系,AM1計算就較為耗時,而上千個原子的體系則十分困難.但是如果使用線性標度方法,如MOZYME方法95來加速AM1步驟,則AM1-BCC電荷的適用范圍可以被大大擴展.

得到波函數(shù)后,計算Mulliken電荷不需要額外耗時,而計算NPA、MK、Hirshfeld和ADCH電荷的計算耗時也只占波函數(shù)計算耗時的很小比例.因此可以近似認為計算這幾種電荷的總耗時就是計算波函數(shù)的耗時.

使用AOIM 1.1程序計算基組較大或原子數(shù)稍多體系的AOIM電荷比較慢.例如在HF/6-31G*下計算Ephedrine(C10H15NO)分子的AOIM電荷耗費約8 min,而使用默認參數(shù)、單線程下使用Gaussian 03計算這個分子波函數(shù)耗時40 s.我們認為計算速度緩慢不是AOIM理論本身的原因,而是因為AOIM 1.1代碼中的L-BFGS-B最小化模塊效率有待改善.對代碼和算法進行優(yōu)化和并行化后有望使AOIM電荷計算速度顯著加快.

AIM電荷計算耗時極長.我們使用計算效率較高的AIM程序AIMALL,并采用默認的積分精度和積分算法時,計算Ephedrine分子全部AIM電荷總耗時為10.5 min,是計算波函數(shù)時間的近16倍.AIM電荷計算費時是因為AIM的原子邊界復雜,導致原子盆不容易被準確積分,盡管已有不少研究者探索出更有效的原子盆積分方法,96-98但AIM電荷計算緩慢的問題始終沒有被徹底解決.

CM2的校正步驟幾乎不需要額外時間,總耗時取決于使用何種理論級別計算初始電荷.如果要求計算快速,可以用半經(jīng)驗方法,如果要求對實際氣相偶極矩有更好的重現(xiàn)能力,則可以用從頭算方法.實際上從原文獻的統(tǒng)計數(shù)據(jù)來看,半經(jīng)驗與從頭算方法所得CM2電荷的偶極矩重現(xiàn)能力對于多數(shù)體系差距并不很大.

10 計算方法的選擇

結(jié)合上文對各種原子電荷計算方法的對比結(jié)果、方法特點以及目前應用的現(xiàn)狀,我們對原子電荷計算方法的選擇給出以下建議.

對于分子體系的量子化學研究,建議使用NPA和ADCH電荷.NPA電荷目前應用十分廣泛,合理性總體表現(xiàn)較好.但注意不宜使用NPA電荷分析靜電相互作用,同時注意NPA用于過渡金屬、鑭系錒系金屬時存在一定問題.23,99問題主要是由于這些元素的極小集/里德堡集劃分標準比較含糊,而不同的劃分又會影響OWSO過程,進而影響電荷數(shù)值.最近提出的ADCH方法無論在電荷合理性、基組依賴性、偶極矩重現(xiàn)性和靜電勢的重現(xiàn)性上都表現(xiàn)不錯,電荷計算時間也遠小于波函數(shù)計算時間,是十分有前景的原子電荷計算方法,其普適性、穩(wěn)定性有待在更多的實際應用中進行檢驗.

基于原子電荷的分子模擬任務需要原子電荷能夠較好地描述靜電相互作用.通常來說以MK為代表的擬合靜電勢電荷是最適合用于分子模擬的.但如果模擬的是有機分子,更適宜使用RESP電荷,因為它解決了普通擬合靜電勢方法存在的構(gòu)象依賴性、內(nèi)部原子電荷數(shù)值不穩(wěn)定性和原子等價性問題,若這些問題存留則可能會影響模擬過程的合理性.注意由于在擬合力場參數(shù)時范德華作用、1-4位相互作用和靜電作用是相互耦合的,故范德華參數(shù)和二面角參數(shù)對電荷模型有依賴性,因此若模擬任務依賴于已有力場參數(shù),則原子電荷的計算方法應當盡量與力場原文中的方法一致.

對于較大分子體系或者大量有機小分子體系的原子電荷的計算建議使用AM1-BCC方法,此方法表現(xiàn)靜電作用較好且計算速度較快.如果需要進一步降低計算耗時,可以改用MMFF94方法.另外,如果所研究的大分子體系結(jié)構(gòu)具有規(guī)律性,如蛋白質(zhì)、多糖、核酸、高分子,也可以使用能夠更好重現(xiàn)靜電勢的MK方法獲得每個片段的電荷再進行拼接.

估算氣相分子偶極矩最適合使用CM2方法,可以以較低的計算花費獲得很精確的氣相偶極矩.

最后,本文對目前使用廣泛的Mulliken、AIM和Gasteiger電荷進行簡要評述.

Mulliken方法盡管存在原理不嚴格、低估離子性化合物的鍵的極性、基組依賴性大、偶極矩和靜電勢重現(xiàn)性差等諸多缺陷,但是從本文的測試看到,Mulliken電荷對于多數(shù)小分子并沒有出現(xiàn)過于明顯的不合理性,又由于Mulliken電荷幾乎被所有量子化學軟件支持,不耗費額外計算時間,因此Mulliken電荷雖然不推薦使用,但仍可以作為定性的參考.在極小基下產(chǎn)生的Mulliken電荷比在更高級別基組下產(chǎn)生結(jié)果在原理上更合理.如果所用基組包含彌散函數(shù),則應當舍棄Mulliken電荷.另外,由于Mulliken電荷受基組影響明顯,比較Mulliken電荷時應注意必須處于在相同基組下,否則結(jié)果沒有意義.

AIM方法一般不建議使用.盡管它擁有嚴格的物理解釋,但是在測試中它的結(jié)果往往明顯違背化學觀念,對偶極矩和靜電勢重現(xiàn)性都很差,當體系存在贗原子時無法使用,而且計算耗時極長.

Gasteiger電荷由于具有計算方便、快速的特點,被很多藥物設(shè)計、分子對接軟件所采用,常被用于描述受體-配體靜電相互作用、生成CoMFA任務所需要的靜電勢場.然而從前文的對比結(jié)果可見它的靜電勢重現(xiàn)性并不理想,因此并不推薦Gasteiger電荷用于這些場合.相比之下使用AM1-BCC和MMFF94電荷更為合適.

11 總 結(jié)

原子電荷是十分古老、簡單的模型,可以追溯到化合價概念的提出.隨著理論化學的發(fā)展以及人們對電子結(jié)構(gòu)認識的加深,原子電荷模型不僅沒有被遺棄,新的計算方法還在不斷涌現(xiàn),這是由于它有著重要的理論和實用意義.本文介紹了目前使用最為廣泛的12種計算原子電荷的方法,通過大量分析對它們的特點進行了全面的比較和評述.對比分析中顯示,許多計算方法的結(jié)果間存在很大差異,錯誤地選用計算方法可能得到不可靠甚至沒有意義的原子電荷.只有充分掌握了眾多計算方法各自的原理和特點,才能根據(jù)實際問題選擇最為適用的一種.

(1)Qian,B.H.;Ma,W.X.;Lu,L.D.;Yang,X.J.;Wang,X.Acta Phys.-Chim.Sin.2010,26,610.[錢保華,馬衛(wèi)興,陸路德,楊緒杰,汪 信.物理化學學報,2010,26,610.]

(2) Zheng,W.R.;Xu,J.L.;Xiong,R.Acta Phys.-Chim.Sin.2010, 26,2535.[鄭文銳,徐菁利,熊 瑞.物理化學學報,2010,26, 2535.]

(3)Shen,T.;Du,F.P.;Liu,T.;Yao,G.W.;Wu,Z.;Fang,M.M.; Xu,X.J.;Lu,H.Z.Acta Phys.-Chim.Sin.2011,27,1831. [申 濤,杜鳳沛,劉 婷,姚廣偉,吳 崢,方萌萌,徐筱杰,路慧哲.物理化學學報,2011,27,1831.]

(4) Zhou,J.J.;Chen,H.M.;Xie,G.R.;Ren,T.R.;Xu,Z.H.Prog. Chem.1998,10,55.[周家駒,陳紅明,謝桂榮,任天瑞,許志宏.化學進展,1998,10,55.]

(5) Ji,G.D.;Zhao,Y.H.;Yuan,X.J.Northeast Normal Univ. (Natural Science Edition)1998,47. [籍國東,趙元慧,袁 星.東北師范大學學報(自然科學版),1998,47.]

(6) Ding,Y.F.;Zhang,Y.;Zhang,D.H.;Li,Z.P.Acta Phys.-Chim. Sin.2010,26,1651.[丁元法,張 躍,張大海,李仲平.物理化學學報,2010,26,1651.]

(7) Laio,A.;VandeVondele,J.;Rothlisberger,U.J.Phys.Chem.B 2002,106,7300.

(8) Pipek,J.;Mezey,P.G.J.Chem.Phys.1989,90,4916.

(9) Elstner,M.;Porezag,D.;Jungnickel,G.;Elsner,J.;Haugk,M.; Frauenheim,T.;Suhai,S.;Seifert,G.Phys.Rev.B 1998,58, 7260.

(10) Giesen,D.J.;Cramer,C.J.;Truhlar,D.G.J.Phys.Chem.1995, 99,7137.

(11) Giesen,D.J.;Hawkins,G.D.;Liotard,D.A.;Cramer,C.J.; Truhlar,D.G.Theor.Chem.Acc.1997,98,85.

(12) Li,J.B.;Hawkins,G.D.;Cramer,C.J.;Truhlar,D.G.Chem. Phys.Lett.1998,288,293.

(13) Thompson,J.D.;Cramer,C.J.;Truhlar,D.G.J.Phys.Chem.A 2004,108,6532.

(14) Meister,J.;Schwarz,W.H.E.J.Phys.Chem.1994,98,8245.

(15) Cramer,C.J.Essentials of Computational Chemistry,2nd ed.; John Wiley&Sons:West Sussex,2004;pp 309-324.

(16) Jensen,F.Introduction to Computational Chemistry,2nd ed.; John Wiley&Sons:West Sussex,2007;pp 293-304.

(17)Young,D.C.Computational Chemistry;John Wiley&Sons: New York,2001;pp 99-105.

(18) Cioslowski,J.Electronic WavefunctionAnalysis.In Encyclopedia of Computational Chemistry;Schleyer,P.v.R. Ed.;John Wiley&Sons:West Sussex,1998;Vol.2,pp 892-905.

(19) Mulliken,R.S.J.Chem.Phys.1955,23,1841.

(20) Mulliken,R.S.J.Chem.Phys.1955,23,1833.

(21) Mulliken,R.S.J.Chem.Phys.1955,23,2338.

(22) Bachrach,S.M.PopulationAnalysis and Electron Densities from Quantum Mechanics.In Reviews in Computational Chemistry;Lipkowitz,K.B.,Boyd,D.B.Eds.;VCH Publishers:New York,1994;Vol.5,pp 171-227.

(23) Clark,A.E.;Sonnenberg,J.L.;Hay,P.J.;Martin,R.L. J.Chem.Phys.2004,121,2563.

(24) Martin,F.;Zipse,H.J.Comput.Chem.2005,26,97.

(25) Wiberg,K.B.;Rablen,P.R.J.Comput.Chem.1993,14,1504.

(26)Lu,H.G.;Dai,D.D.;Yang,P.;Li,L.M.Phys.Chem.Chem. Phys.2006,8,340.

(27) Hirshfeld,F.L.Theor.Chem.Acc.1977,44,129.

(28)Lu,T.;Chen,F.W.J.Theor.Comput.Chem.accepted.

(29) Reed,A.E.;Weinstock,R.B.;Weinhold,F.J.Chem.Phys. 1985,83,735.

(30)Besler,B.H.;Merz,K.M.,Jr.;Kollman,P.A.J.Comput.Chem. 1990,11,431.

(31)Bader,R.F.W.;Beddall,P.M.J.Chem.Phys.1972,56,3320.

(32)Halgren,T.A.J.Comput.Chem.1996,17,520.

(33) Halgren,T.A.J.Comput.Chem.1996,17,616.

(34) Jakalian,A.;Bush,B.L.;Jack,D.B.;Bayly,C.I.J.Comput. Chem.2000,21,132.

(35) Jakalian,A.;Jack,D.B.;Bayly,C.I.J.Comput.Chem.2002, 23,1623.

(36) Gasteiger,J.;Marsili,M.Tetrahedron 1980,36,3219.

(37) Li,J.B.;Zhu,T.H.;Cramer,C.J.;Truhlar,D.G.J.Phys. Chem.A 1998,102,1820.

(38)Rappe,A.K.;Goddard,W.A.J.Phys.Chem.1991,95,3358.

(39) Cusachs,L.C.;Politzer,P.Chem.Phys.Lett.1968,1,529.

(40) Stout,E.W.;Politzer,P.Theor.Chem.Acc.1968,12,379.

(41) Doggett,G.J.Chem.Soc.A 1969,229.

(42) Christoffersena,R.E.;Baker,K.A.Chem.Phys.Lett.1971,8,4.

(43) Bickelhaupt,F.M.;van Eikema Hommes,N.J.R.;Fonseca Guerra,C.;Baerends,E.J.Organometallics 1996,15,2923.

(44) Weinhold,F.Natural Bond Orbital Methods.In Encyclopedia of Computational Chemistry;Schleyer,P.v.R.Ed.;John Wiley& Sons:West Sussex,1998;Vol.2,pp 1792-1811.

(45) Glendening,E.D.;Badenhoop,J.K.;Reed,A.E.;Carpenter,J. E.;Bohmann,J.A.;Morales,C.M.;Weinhold,F.NBO,Version 5.0,2001.http://www.chem.wisc.edu/～nbo5/.

(46) Liu,W.;Li,L.Theor.Chem.Acc.1997,95,81.

(47) Sanchez-Portal,D.;Artacho,E.;Soler,J.M.Solid State Commun.1995,95,685.

(48) Sanchez-Portal,D.;Artacho,E.;Soler,J.M.J.Phys.:Condens. Matter 1996,8,3859.

(49) Bader,F.W.Atoms in Molecules:A Quantum Theory;Oxford University Press:New York,1994.

(50) Nalewajski,R.F.;Parr,R.G.Proc.Natl.Acad.Sci.U.S.A. 2000,97,8879.

(51) Davidson,E.R.;Chakravorty,S.Theor.Chem.Acc.1992,83, 319.

(52) Chirlian,L.E.;Francl,M.M.J.Comput.Chem.1987,8,894.

(53)Breneman,C.M.;Wiberg,K.B.J.Comput.Chem.1990,11, 361.

(54) Sigfridsson,E.;Ryde,U.J.Comput.Chem.1998,19,377.

(55) Bayly,C.I.;Cieplak,P.;Cornell,W.;Kollman,P.A.J.Phys. Chem.1993,97,10269.

(56) Dewar,M.J.S.;Zoebisch,E.G.;Healy,E.F.;Stewart,J.J.P. J.Am.Chem.Soc.1985,107,3902.

(57) Storer,J.W.;Giesen,D.J.;Cramer,C.J.;Truhlar,D.G. J.Comput.-Aided Mol.Des.1995,9,87.

(58) Li,J.B.;Williams,B.;Cramer,C.J.;Truhlar,D.G.J.Chem. Phys.1999,110,724.

(59) Thompson,J.D.;Cramer,C.J.;Truhlar,D.G.J.Comput. Chem.2003,24,1291.

(60) Winget,P.;Thompson,J.D.;Xidos,J.D.;Cramer,C.J.; Truhlar,D.G.J.Phys.Chem.A 2002,106,10707.

(61) Kalinowski,J.A.;Lesyng,B.;Thompson,J.D.;Cramer,C.J.; Truhlar,D.G.J.Phys.Chem.A 2004,108,2545.

(62) Olson,R.M.;Marenich,A.V.;Cramer,C.J.;Truhlar,D.G. J.Chem.Theory Comput.2007,3,2046.

(63) Kelly,C.P.;Cramer,C.J.;Truhlar,D.G.J.Chem.Theory Comput.2005,1,1133.

(64) Mayer,I.Chem.Phys.Lett.1983,97,270.

(65) Sanderson,R.T.Science 1951,114,670.

(66) Rappe,A.K.;Casewit,C.J.;Colwell,K.S.;Goddard,W.A.; Skiff,W.M.J.Am.Chem.Soc.1992,114,10024.

(67)Mortier,W.J.;Ghosh,S.K.;Shankar,S.J.Am.Chem.Soc. 1986,108,4315.

(68) Cioslowski,J.Phys.Rev.Lett.1989,62,1469.

(69) Cioslowski,J.J.Am.Chem.Soc.1989,111,8333.

(70) Szabo,A.;Ostlund,N.S.Modern Quantum Chemistry,1st rev ed.;Dover Publications:New York,1989.

(71) Hariharan,P.C.;Pople,J.A.Theor.Chem.Acc.1973,28,213.

(72) Hehre,W.J.;Ditchfield,R.;Pople,J.A.J.Chem.Phys.1972, 56,2257.

(73) Frisch,M.J.;Pople,J.A.;Binkley,J.S.J.Chem.Phys.1984, 80,3265.

(74) Frisch,M.J.;Trucks,G.W.;Schlegel,H.B.;et al.Gaussian 03, Revison E.01;Gaussian Inc.:Wallingford,CT,2004.

(75) Lu,T.Multiwfn,Version 2.1.2;2011.http://Multiwfn.codeplex. com.

(76) Keith,T.A.AIMALL,Version 10.05.04,2010.

(77) Lu,H.G.AOIM,Version 1.1,2006;http://faculty.sxu.cn/luhg/ aoim.html.

(78)Avogadro:an Open-Source Molecular Builder and Visualization Tool,Version 1.0.3,2011.

(79) Case,D.A.;Darden,T.A.;Cheatham,T.E.C.,III;et al. AmberTools,Version 1.5;2011.

(80) Sch?fer,A.;Horn,H.;Ahlrichs,R.J.Chem.Phys.1992,97, 2571.

(81) PETRA Manual.http://www2.ccc.uni-erlangen.de/software/ petra/manual(accessed Sep 12,2011).

(82) Marsili,M.;Gasteiger,J.Croat.Chem.Acta 1980,53,601.

(83) The Open Babel Package,Version 2.3.0;2010;http://openbabel. sourceforge.net.

(84) Sanner,M.F.J.Mol.Graph.Model.1999,17,57.

(85)Woods,R.J.;Khalil,M.;Pell,W.;Moffat,S.H.;Smith,V.H. J.Comput.Chem.1990,11,297.

(86) Hehre,W.J.;Stewart,R.F.;Pople,J.A.J.Chem.Phys.1969, 51,2657.

(87) Binkley,J.S.;Pople,J.A.;Hehre,W.J.J.Am.Chem.Soc.1980, 102,939.

(88) Krishnan,R.;Binkley,J.S.;Seeger,R.;Pople,J.A.J.Chern. Phys.1980,72,650.

(89) Dunning,J.T.H.J.Chem.Phys.1989,90,1007.

(90) Kendall,R.A.;Dunning,T.H.;Harrison,R.J.J.Chem.Phys. 1992,96,6796.

(91) Becke,A.D.Phys.Rev.A 1988,38,3098.

(92) Perdew,J.P.Phys.Rev.B 1986,33,8822.

(93) Becke,A.D.J.Chem.Phys.1993,98,1372.

(94) Patel,S.;Brooks,C.L.Mol.Simul.2006,32,231.

(95) Stewart,J.J.P.Int.J.Quantum Chem.1996,58,133.

(96) Biegler-K?nig,F.W.J.Comput.Chem.2000,21,1040.

(97) Biegler-K?nig,F.W.;Bader,R.F.W.;Tang,T.H.J.Comput. Chem.1982,3,317.

(98)Sanville,E.;Kenny,S.D.;Smith,R.;Henkelman,G.J.Comput. Chem.2007,28,899.

(99) Maseras,F.;Morokuma,K.Chem.Phys.Lett.1992,195,500.

September 13,2011;Revised:October 25,2011;Published on Web:October 31,2011.*

.Email:chenfeiwu@ustb.edu.cn;Tel:+86-10-62332689.

Comparison of Computational Methods for Atomic Charges

LU Tian CHEN Fei-Wu*
(School of Chemical and Biological Engineering,University of Science and Technology Beijing,Beijing 100083,P.R.China)

Atomic charge is one of the simplest and the most intuitive description of charge distribution in chemical systems.It has great significance in theory and in practical applications.In this article we introduce the basic principles and special characteristics of twelve important computational methods for the determination of atomic charges and compare their pros and cons from various aspects by considering a large number of instances.These methods include Mulliken,atomic orbitals in molecules(AOIM), Hirshfeld,atomic dipole moment corrected Hirshfeld population(ADCH),natural population analysis(NPA), Merz-Kollmann(MK),atom in molecules(AIM),Merck molecular force field 94(MMFF94),AM1-BCC, Gasteiger,charge model 2(CM2),and charge equilibration(QEq).Finally some general suggestions on how to choose a proper method for practical applications are given.

Atomic charge;Computational chemistry;Population analysis;Electronegativity; Electrostatic potential

10.3866/PKU.WHXB2012281

The project was supported by the National Natural Science Foundation of China(20773011).

國家自然科學基金(20773011)資助項目

O641