盧 天 陳飛武

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

原子電荷計(jì)算方法的對(duì)比

盧 天 陳飛武*

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

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

原子電荷;計(jì)算化學(xué);布居分析;電負(fù)性;靜電勢(shì)

1 引言

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

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

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

2 各種原子電荷計(jì)算方法的原理和特點(diǎn)

2.1 Mulliken方法

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

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

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

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

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

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

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

2.2 NPA方法

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

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

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

2.3 AOIM方法

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

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

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

2.4 AIM方法

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

2.5 Hirshfeld方法

Hirshfeld電荷寫(xiě)為

其中

2.6 ADCH方法

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

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

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

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

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

2.7 MK方法

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

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

2.8 AM1-BCC方法

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

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

2.9 MMFF94方法

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

2.10 CM2方法

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

CM2校正的公式與AM1-BCC和MMFF94的相比略為復(fù)雜,CM2電荷可寫(xiě)為

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

2.11 Gasteiger方法

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

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

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

2.12 QEq方法

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

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

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

3 計(jì)算細(xì)節(jié)和測(cè)試體系

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

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

4 普通小分子的原子電荷

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

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

continued Table 1

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

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

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

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

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

continued Table 2

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

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

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

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

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

圖3 第一組方法計(jì)算的不同取代數(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 第二組方法計(jì)算的不同取代數(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)性

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

表3 第一組方法計(jì)算的20種小分子的分子偶極矩(單位為Debye)以及與HF/6-31G＊級(jí)別下偶極矩算符期望值(〈μ〉)的比較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ù)過(guò)程中,MMFF94的參數(shù)主要來(lái)自于使MMFF94電荷產(chǎn)生的偶極矩逼近HF/6-31G*密度產(chǎn)生的分子偶極矩,AM1-BCC的參數(shù)來(lái)自于使AM1電荷在校正后逼近HF/6-31G*級(jí)別下計(jì)算的RESP電荷.它們都使用HF/6-31G*下的數(shù)據(jù)作為目標(biāo)數(shù)據(jù)是因?yàn)檫@個(gè)波函數(shù)級(jí)別對(duì)偶極矩的高估被認(rèn)為可以等效反映出分子在凝聚相中存在的極化效應(yīng),因此MMFF94和AM1-BCC的偶極矩應(yīng)當(dāng)與HF/6-31G*下的偶極矩期望值相比較.從表4可見(jiàn)這兩種電荷的偶極矩重現(xiàn)性對(duì)大部分分子基本可以令人滿意,MMFF94的MUE為0.36 Debye,雖然AM1-BCC和MMFF94在鍵電荷校正過(guò)程中很相似,但AM1-BCC的MUE更小,為0.25 Debye,其改善在一定程度上可能是由于AM1-BCC使用了更有意義的初始電荷.Gasteiger的偶極矩重現(xiàn)性無(wú)論以氣相實(shí)驗(yàn)值為參考還是以HF/6-31G*偶極矩為參考都顯得一般,并且對(duì)多數(shù)分子有低估偶極矩的傾向.CM2電荷能夠很好地重現(xiàn)氣相實(shí)驗(yàn)偶極矩, MUE僅為0.15 Debye,對(duì)所有分子都沒(méi)有出現(xiàn)過(guò)大誤差,這與CM2提出的目標(biāo)完全相一致.QEq對(duì)實(shí)驗(yàn)偶極矩和HF/6-31G*下的偶極矩重現(xiàn)性都很不理想,對(duì)許多分子偶極矩高估了近一倍甚至一倍以上,嚴(yán)重低估偶極矩的情況同樣存在,如2-甲基嘧啶.

表4 第二組方法計(jì)算的20種小分子的分子偶極矩(單位為Debye)以及與HF/6-31G＊級(jí)別下偶極矩算符期望值和氣相實(shí)驗(yàn)值的比較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電荷,對(duì)偶極矩重現(xiàn)性都不好.ADCH和MK都能很準(zhǔn)確重現(xiàn)波函數(shù)的偶極矩期望值,若波函數(shù)完全精確,則這兩種電荷計(jì)算的偶極矩將與氣相的實(shí)際值相一致,而利用CM2方法則可以在較低波函數(shù)級(jí)別下很準(zhǔn)確地預(yù)測(cè)氣相偶極矩.MMFF94和AM1-BCC能夠較好地重現(xiàn)HF/6-31G*的偶極矩.

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

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

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

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

表5 第一組方法產(chǎn)生的靜電勢(shì)相對(duì)于HF/6-31G＊級(jí)別下精確靜電勢(shì)的相對(duì)均方根偏差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

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

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

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

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

表6 第二組方法產(chǎn)生的靜電勢(shì)相對(duì)于HF/6-31G＊級(jí)別下精確靜電勢(shì)的相對(duì)均方根偏差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級(jí)別下獲得的乙酸的甲基碳原子電荷Fig.5 Atomic charge of methyl carbon of acetic acid calculated by different methods under various basis-sets at the Hartree-Fock level

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

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

9 計(jì)算的耗時(shí)

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

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

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

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

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

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

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

10 計(jì)算方法的選擇

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

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

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

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

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

最后,本文對(duì)目前使用廣泛的Mulliken、AIM和Gasteiger電荷進(jìn)行簡(jiǎn)要評(píng)述.

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

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

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

11 總 結(jié)

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

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

(2) Zheng,W.R.;Xu,J.L.;Xiong,R.Acta Phys.-Chim.Sin.2010, 26,2535.[鄭文銳,徐菁利,熊 瑞.物理化學(xué)學(xué)報(bào),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. [申 濤,杜鳳沛,劉 婷,姚廣偉,吳 崢,方萌萌,徐筱杰,路慧哲.物理化學(xué)學(xué)報(bào),2011,27,1831.]

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

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

(6) Ding,Y.F.;Zhang,Y.;Zhang,D.H.;Li,Z.P.Acta Phys.-Chim. Sin.2010,26,1651.[丁元法,張 躍,張大海,李仲平.物理化學(xué)學(xué)報(bào),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).

國(guó)家自然科學(xué)基金(20773011)資助項(xiàng)目

O641