用第一性原理方法獲取周期體系中原子的部分電荷

李亞娜 呂 洋 周立川 陳 理 李慎敏

(大連大學(xué),遼寧省生物有機(jī)化學(xué)重點(diǎn)實(shí)驗(yàn)室,遼寧大連 116622)

用第一性原理方法獲取周期體系中原子的部分電荷

李亞娜 呂 洋 周立川 陳 理 李慎敏*

(大連大學(xué),遼寧省生物有機(jī)化學(xué)重點(diǎn)實(shí)驗(yàn)室,遼寧大連 116622)

利用量子化學(xué)軟件包Crystal計(jì)算了立方周期性邊界條件下液態(tài)水體系的靜電勢(ESP)和靜電場(EF).在此基礎(chǔ)上,提出了一種由第一性原理方法獲取周期體系中原子的部分電荷的快捷方法.該方法把由周期性邊界條件引入的平均靜電勢φmean作為一個擬合參數(shù),通過對第一性原理靜電勢與Ewald加和法靜電勢的最小二乘法擬合而實(shí)現(xiàn).值得說明的是,比較靜電勢與靜電場擬合方法,前者的相對擬合誤差僅為2%-3%,比后者小一個數(shù)量級.考察了四種電荷限制條件下,靜電勢、靜電場擬合的水分子原子部分電荷及偶極矩的分布情況.

原子的部分電荷; 最小二乘擬合; 靜電勢; 靜電場;Ewald加和法

對于化學(xué)工作者來說,原子的部分電荷是一個具有重要意義的不可或缺的概念.許多用于描述體系性質(zhì)的物理、化學(xué)概念或物理量,如分子的極性(偶極矩)、反應(yīng)性(活性)、靜電相互作用等都與體系中電子的分布有關(guān),而描述電子分布的簡便、有效、經(jīng)典的方法是原子的部分電荷(以下簡稱原子電荷).此外,原子電荷也是計(jì)算化學(xué)方法,如分子模擬、蒙特卡洛模擬技術(shù)等所依賴的描述體系電學(xué)性質(zhì)的重要力場參數(shù).然而,原子電荷卻又是一個沒有明確的量子力學(xué)定義,不能被實(shí)驗(yàn)直接測量的人為的物理量.這里的“人為”體現(xiàn)在借助量子力學(xué)或?qū)嶒?yàn)方法獲得的原子電荷往往具有一定的隨意性和不確定性[1].

隨著計(jì)算機(jī)技術(shù)的快速發(fā)展,量子力學(xué)第一性原理方法已成為獲得原子電荷的主要手段.目前常見的方法有:基于分子軌道(Mulliken電荷)[2]、電子密度(Bader電荷)[3]的空間劃分方法;基于原子周圍靜電勢(ESP)[4-8]、靜電場(EF)[9]的數(shù)值擬合方法等. Mulliken電荷的缺點(diǎn)在于其對計(jì)算方法、基組的依賴性較大;Bader電荷的基組依賴性雖小,但它只給出Bader體積內(nèi)的凈電荷,一般不用作模擬力場的電荷參數(shù);靜電勢電荷因其可以較好地描述分子間的相互作用,特別是近年來發(fā)展的參數(shù)限制的靜電勢擬合(RESP)電荷的使用[8,10],使得基于第一性原理靜電勢的擬合電荷已成為目前許多力場的電荷參數(shù),如著名的AMBER力場[11].然而,靜電勢擬合方法不適合于凝聚相周期體系中原子電荷的獲取,原因在于由周期性邊界條件所引入的體系的平均靜電勢是未知的.基于此,Spackman等[9]提出了靜電場擬合電荷的方法.由于靜電場是靜電勢梯度的負(fù)值,因而該方法避開了平均靜電勢的求算.不過,靜電場的第一性原理計(jì)算以及作為矢量的靜電場的擬合,其計(jì)算工作量要遠(yuǎn)大于靜電勢擬合.

對于離子或強(qiáng)極性凝聚相周期體系,體系中原子的電荷明顯依賴于分子的幾何構(gòu)象及周圍環(huán)境的變化,換句話說,體系的運(yùn)動過程中存在著顯著的電荷極化和電荷轉(zhuǎn)移現(xiàn)象,而這是傳統(tǒng)的固定電荷力場方法無法正確描述的.完全的第一性原理方法處理工作量巨大,目前還難以實(shí)現(xiàn).對于無限稀溶質(zhì)-溶劑體系、生物質(zhì)-水溶液體系等,變通的方法是通過把系統(tǒng)劃分為量子-經(jīng)典兩個作用區(qū)域,對于所關(guān)心的溶質(zhì)、生物大分子活性區(qū)域等進(jìn)行量子化學(xué)處理，而對其他區(qū)域進(jìn)行經(jīng)典力學(xué)處理，即所謂的量子-經(jīng)典混合方法[12-13].對于不易劃分區(qū)域的研究體系,如水等純?nèi)軇w,類似的應(yīng)用還鮮見報(bào)道.

為了考察上述研究體系的極化作用,近年來,一種發(fā)展較快的浮動電荷力場,即第三代力場受到了越來越多人的關(guān)注[14-15].該方法利用Sanderson的電負(fù)性均衡原理(EEM)[16-17]快速計(jì)算體系中原子、分子的瞬態(tài)電荷,因而可以考察體系中分子間極化作用及電荷轉(zhuǎn)移作用.該方法計(jì)算電荷的準(zhǔn)確性依賴于所研究體系中原子的價(jià)態(tài)電負(fù)性及價(jià)態(tài)硬度參數(shù),這些參數(shù)通常來自于對第一性原理獲取的點(diǎn)電荷的擬合,與擬合訓(xùn)練集組成、大小以及第一性原理所獲取電荷的準(zhǔn)確性密切相關(guān).

水溶劑作為常見的極性溶劑一直備受人們的關(guān)注.目前用于分子模擬的水分子模型多為固定電荷模型,如SPC水,TIPnP(n=3-5)系列水等[18-20].最近Yang等利用原子-鍵電負(fù)性均衡原理(ABEEM)方法構(gòu)建了七點(diǎn)的浮動電荷水模型[21],并將其應(yīng)用于水團(tuán)簇及生物大分子水溶液的模擬中[22-24].Yang的七點(diǎn)水模型在處理水的簇合物中得到了與第一性原理相媲美的計(jì)算結(jié)果[21,25].特別地,當(dāng)將鍵電荷回歸到原子位點(diǎn)時,原子電荷與量子化學(xué)計(jì)算的Mulliken電荷具有相當(dāng)好的線性關(guān)系.

本文中,我們將周期性邊界條件下的液態(tài)水體系作為研究對象,利用第一性原理方法計(jì)算體系的網(wǎng)格靜電勢和靜電場;然后,利用點(diǎn)電荷的直接截?cái)嗉雍头?Rcut)、Ewald加和法靜電勢[26]以及靜電場模型，通過電荷限制條件下的最小二乘擬合法擬合第一性原理靜電勢、靜電場的計(jì)算數(shù)據(jù),進(jìn)而獲取周期水體系的原子部分電荷.我們的目的有二:一是通過第一性原理方法獲取的靜電勢電荷與靜電場電荷的比較,尋找適合于凝聚相體系原子電荷獲取的快速方法;二是在分子動力學(xué)模擬環(huán)境(周期性邊界條件)下,考察水體系中電荷極化和電荷轉(zhuǎn)移作用,為下一步在EEM理論框架內(nèi)進(jìn)行價(jià)態(tài)電負(fù)性及價(jià)態(tài)硬度的參數(shù)化提供可靠的訓(xùn)練集數(shù)據(jù).

1 理論方法與計(jì)算細(xì)節(jié)

1.1 第一性原理計(jì)算靜電勢與靜電場

目前適用于周期體系第一性原理計(jì)算的商用軟件包并不多,而直接包含靜電場計(jì)算的則更少, Dovesi等[27]發(fā)展的Crystal軟件包是其中的一個.本文中,使用Crystal 06軟件包在HF/6-21G基組下,計(jì)算了立方體周期性邊界條件下液態(tài)水體系的網(wǎng)格靜電勢和靜電場.具體步驟如下:首先,利用常規(guī)的MD軟件包GROMACS[28],選取固定電荷的SPC水模型,在NVT系綜下(T=298 K)分別對兩個尺寸的液態(tài)水盒子進(jìn)行了500 ps的動力學(xué)模擬.模擬中,水的密度 ρ=1.0 g·cm-3,水盒子邊長 a為 1.0,1.1 nm,與其對應(yīng)的水分子數(shù)分別為34,45.然后,針對模擬平衡后的兩個體系的坐標(biāo)構(gòu)象進(jìn)行數(shù)據(jù)采集,通過隨機(jī)取樣方法分別選取10個水盒子構(gòu)象為樣本數(shù)據(jù).最后,利用Crystal程序計(jì)算每個樣本構(gòu)象的單點(diǎn)能,水分子的Mulliken電荷,以及每個原子周圍給定網(wǎng)格點(diǎn)的靜電勢和靜電場.其中,第一性原理靜電勢和靜電場的采集坐標(biāo)(網(wǎng)格點(diǎn))來自于水盒子(元胞)內(nèi)每個原子范德華半徑以外的均勻分布三維網(wǎng)格點(diǎn),網(wǎng)格點(diǎn)間距為0.05 nm.對于邊長為1.0,1.1 nm的兩個體系,元胞內(nèi)網(wǎng)格點(diǎn)數(shù) Np分別是3388和3483.圖1給出了a=1.0 nm體系的計(jì)算網(wǎng)格點(diǎn)示意圖.

1.2 周期體系的點(diǎn)電荷靜電勢和靜電場

對于孤立分子或分子簇,無窮遠(yuǎn)處的靜電勢為零.然而,對于無限重復(fù)的周期體系(如凝聚相體系的常用處理方法),由于電子密度分布的周期性,無窮遠(yuǎn)處的靜電勢不為零.不過,當(dāng)元胞內(nèi)體系中總的凈電荷為零時,元胞內(nèi)靜電勢φ(r)通常采用如下邊界條件[9],

式(1)也是Crystal程序中元胞邊界條件的選擇方式[29],這與Ewald加和法處理長程靜電相互作用的邊界條件是一致的[30].為滿足式(1)的邊界條件,周期體系的靜電勢可以寫成如下兩部分:

上式中φmean是一個與位置無關(guān)的常數(shù)項(xiàng),φ′(r)是由點(diǎn)電荷模型定義的周期體系中元胞內(nèi)位置r處的靜電勢:

其中,ri是元胞內(nèi)原子i的位置矢量,qi是原子i的部分電荷,l是水盒子的尺寸,N是元胞中原子的個數(shù),k表示重復(fù)單元的序號.當(dāng)元胞內(nèi)凈電荷為零時,式(3)是收斂的,我們稱其為靜電勢的截?cái)嗉雍头?但其收斂速度較慢.本文中,以元胞為中心,沿坐標(biāo)軸三個方向分別取k為-10到10的21×21×21個重復(fù)單元.靜電勢的Ewald加和法可以很好地克服式(3)收斂慢的缺點(diǎn),近年來已被廣泛地應(yīng)用于長程靜電相互作用能的計(jì)算.根據(jù)靜電勢的Ewald加和法,元胞內(nèi)靜電勢可改寫為:

上式中,V是元胞體積,α是高斯函數(shù)的寬度系數(shù),用以描述原子位點(diǎn)電荷密度的分布,k是r在倒易空間中的對應(yīng)量,erfc是誤差函數(shù)的補(bǔ)函數(shù).值得注意的是,與靜電勢的Ewald加和法不同,這里的φ′(r)ES只包含F(xiàn)ourier變換項(xiàng)和短程校正項(xiàng)兩項(xiàng),不包含靜電能中的自能項(xiàng)部分.

需要說明的是,由于式(2)中φmean的存在,給利用靜電勢法精確地?cái)M合體系的原子電荷帶來了一定的困難.最近,Spackman等[9]給出了靜電場擬合電荷的方案,其中,E(r)是靜電勢φ(r)梯度的負(fù)值,如(5)式:

也就是說,E(r)與φmean的選取無關(guān),從而避開了φmean的求算,減少了擬合的不確定性.不過,由于E(r)是矢量,相應(yīng)的第一性原理計(jì)算與擬合工作量較靜電勢的計(jì)算與擬合要大許多.

1.3 電荷限制的最小二乘法擬合

通過Lapack軟件包中DGGLSE程序,利用周期體系點(diǎn)電荷靜電勢、靜電場模型,針對Crystal 06計(jì)算給出的水盒子體系10個構(gòu)象的網(wǎng)格靜電勢、靜電場數(shù)據(jù),分別進(jìn)行了四個電荷限制條件下的最小二乘法擬合.擬合評價(jià)判據(jù)由下式定義的相對擬合誤差的百分比給出,

上式中,A0(ri)是第一性原理計(jì)算的網(wǎng)格點(diǎn)ri處的靜電勢或靜電場,Np是元胞內(nèi)網(wǎng)格點(diǎn)數(shù).χ2由式(7)給出:

其中,A(ri)是式(2)-(5)定義的ri處靜電勢或靜電場. Q是元胞內(nèi)中的凈電荷(本文中Q=0),qi是原子i的部分電荷,λ是拉格朗日乘因子.需要說明的是,對于周期體系,由于φmean是未知的,擬合中為獲得靜電勢擬合電荷,我們設(shè)其為一新的擬合參數(shù)qN+1,并構(gòu)造新的擬合判據(jù):

實(shí)際計(jì)算中我們還通過引入多個拉格朗日乘因子的方式嘗試了其它電荷限制條件下的擬合情況.例如,若只允許分子電荷的極化而不允許分子間電荷轉(zhuǎn)移,即限制每個水分子的凈電荷為零.此外,我們還仿照Kollman引入限制函數(shù)χ200rstr,考察了Mulliken電荷限制下的電荷擬合(RESP)[10],這時,

式(9)中χ20就是式(8)中的χ2;式(10)中,qMui00為Crystal 06計(jì)算的原子i的Mulliken電荷.

2 結(jié)果與討論

2.1 不同方法擬合誤差的比較

利用Crystal量子化學(xué)軟件包,首先計(jì)算了立方體周期性邊界條件下液態(tài)水體系的第一性原理網(wǎng)格靜電勢和靜電場.然后,利用式(2)-(5),通過電荷限制的最小二乘擬合獲得了元胞中每個水分子的原子部分電荷.依據(jù)擬合數(shù)據(jù)來源劃分,本文的工作可歸納為靜電勢擬合和靜電場擬合兩種;而依據(jù)對周期體系點(diǎn)電荷靜電勢計(jì)算模型的劃分,靜電勢擬合又可歸納為簡單的截?cái)嗉雍头?Rcut)和Ewald加和法兩種.此外,為了合理地描述體系中電荷極化與轉(zhuǎn)移作用,針對以上三種擬合方法,分別考察了四種電荷限制條件下的擬合情況:(1)元胞中總的凈電荷為零;(2)滿足條件(1),同時加入Mulliken電荷限制; (3)滿足條件(1),同時加入分子的凈電荷為零限制; (4)條件(1)、(2)和(3)同時滿足的情況.不同方法的擬合判據(jù)由10個水盒子構(gòu)象的平均擬合誤差及均方差給出(見表1).較大的體系,即45個水分子元胞中氧、氫原子電荷的擬合結(jié)果見表2,擬合電荷的分布情況分別見圖2和圖3.

由表1我們看到,在靜電勢擬合法中Ewald加和法的相對誤差及標(biāo)準(zhǔn)偏差明顯低于截?cái)嗉雍头?在較小的34個水分子體系中,Ewald加和法最大誤差出現(xiàn)在限制條件4的情況,只有3.75%,比截?cái)嗉雍头ㄗ钚≌`差5.05%(限制條件1)還小.另一方面,靜電場擬合法雖然不涉及平均靜電勢φmean,但由靜電場擬合電荷給出的網(wǎng)格靜電場與第一性原理方法計(jì)算得到的網(wǎng)格數(shù)值相比,相對誤差達(dá)到了30%-40%,遠(yuǎn)大于靜電勢擬合相對誤差.我們認(rèn)為在相同的網(wǎng)格點(diǎn)處,靜電勢與靜電場數(shù)值值域分布的差異是導(dǎo)致上述情況的主要原因.在當(dāng)前的水體系中,由于φmean的存在,靜電勢值域在-0.25--0.05 a.u.,其峰值在φmean附近;而作為矢量的靜電場，其值域落在以原點(diǎn)對稱的-0.10-0.10 a.u.區(qū)間內(nèi),其峰值在0.0 a.u.處.在較多的網(wǎng)格點(diǎn)上靜電場數(shù)值趨近于零的事實(shí)使得式(6)中分母項(xiàng)數(shù)值較小,表現(xiàn)為靜電勢擬合出現(xiàn)較大的相對誤差.事實(shí)上,Spackman等[9]在利用靜電場法擬合尿素晶體時也發(fā)現(xiàn)該方法存在較大的相對擬合誤差.值得一提的是,通過截?cái)嗉雍头o電勢擬合的電荷,再利用式(5)而得到的網(wǎng)格靜電場數(shù)值與第一性原理計(jì)算的靜電場比較,相對誤差僅比靜電場擬合法略大.這暗示著即使是對靜電場計(jì)算有較高要求的模擬體系,考慮到靜電勢擬合計(jì)算效率較靜電場擬合高很多,通過靜電勢擬合獲取電荷的方法也是值得考慮的.這一點(diǎn)也可以從靜電場擬合電荷通過式(3)所計(jì)算的靜電勢與第一性原理比較具有較大的偏差而進(jìn)一步被確認(rèn).

表1 兩個水體系四種電荷限制條件下擬合評價(jià)Table 1 Fitting evaluation of the four charge restrained fits on the two water systems

此外,無論是靜電勢擬合還是靜電場擬合,限制條件的加入會使擬合誤差相應(yīng)地增加.綜合考慮各種擬合情況容易發(fā)現(xiàn),加入限制條件3,即不考慮分子間電荷轉(zhuǎn)移的情況,擬合誤差變化較小,而加入Mulliken電荷限制(RESP限制)后,相對的擬合誤差變化較大.

表1還給出了由靜電勢擬合得到的45個水分子的平均靜電勢.可以看出,不同限制條件下, Ewald加和法得到的樣本均方根偏差均很小,說明平均靜電勢與體系的構(gòu)象相關(guān)性不大,這與我們預(yù)期的結(jié)果相一致.

需要說明的是,利用Ewald加和法擬合電荷,不僅靜電勢的相對誤差和標(biāo)準(zhǔn)偏差小、計(jì)算效率高、收斂快;而且,比較兩個水分子體系的計(jì)算結(jié)果可知,擬合誤差與體系尺寸的依賴性幾乎可以忽略不計(jì).可以說,對于周期性體系,利用Ewald加和法通過擬合第一性原理靜電勢獲取力場中原子部分電荷的方法是首選方法.

最近,Woo等[31]通過引入改進(jìn)的誤差函數(shù)方法,提出了一種利用靜電勢計(jì)算周期體系原子電荷的方法,并成功應(yīng)用于金屬有機(jī)物多孔材料體系中.

2.2 不同方法電荷分布的比較

接下來,我們討論不同方法下,擬合得到的水分子中氧原子、氫原子的部分電荷及其分布情況.對于我們考察的兩個水盒子體系,由于體系尺寸對擬合的結(jié)果影響不大,為節(jié)省篇幅,這里我們只給出較大的45個水分子體系的計(jì)算結(jié)果(見表2及圖2和圖3).容易發(fā)現(xiàn),無論哪種擬合方法,在引入Mulliken電荷限制后(限制條件(2)與(4)),氧、氫擬合電荷的平均值,相對于限制條件(1)與(3),其絕對值均呈減小的趨勢,同時,電荷分布范圍也明顯變窄.而通過限制條件(3)擬合的氧、氫電荷,不考慮電荷轉(zhuǎn)移,具有最大的漲落,以截?cái)喾o電勢擬合為例,氧的最大原子電荷為0.327e,而最小值為-1.423e.由圖2和圖3可見,在限制條件(1)與(3)下,靜電勢截?cái)喾〝M合的電荷呈現(xiàn)最大的分布范圍,而Ewald加和法與靜電場擬合的電荷呈現(xiàn)相似強(qiáng)度的分布;與限制條件(2)與(4)比較,無論是靜電勢擬合還是靜電場擬合,氧原子的電荷具有較大的負(fù)電荷,而氫原子具有較大的正電荷.值得強(qiáng)調(diào)的是引入 Mulliken限制的Ewald加和法擬合的水分子電荷與Crystal程序計(jì)算的Mulliken電荷分布具有良好的一致性.

表2 45個水分子體系中不同擬合方法獲得的原子電荷和分子電荷漲落Table 2 Atomic charges by different fitting methods and its charge fluctuation for each molecule of the 45 water systems

下面,我們討論一下分子的電荷轉(zhuǎn)移問題.在限制條件(1)、(2)下,利用對Crystal程序計(jì)算的靜電勢或靜電場擬合所得到的水分子原子電荷,計(jì)算了體系中水分子的平均凈電荷漲落,ΔQ(見表2)用以衡量水分子中電荷轉(zhuǎn)移能力.我們發(fā)現(xiàn),靜電勢截?cái)嗉雍头〝M合的電荷具有最大的電荷轉(zhuǎn)移,單分子中最大的電荷漲落為0.410e;Ewald加和法擬合的電荷轉(zhuǎn)移最小,單分子中最大的電荷漲落僅為0.084e.

2.3 不同方法水分子偶極矩及其分布的比較

由于分子偶極矩是一個實(shí)驗(yàn)可觀測量,因此,計(jì)算分子偶極矩及其分布已成為評價(jià)凝聚相力場好壞的重要一環(huán).本文中,由于所選體系的Crystal幾何全優(yōu)化工作量巨大,難以實(shí)施,但考慮到我們所研究的水體系中原子的部分電荷可能是決定偶極矩大小關(guān)鍵因素,因此,對于不同擬合方法得到的原子電荷,分別計(jì)算了45個水分子體系的平均偶極矩并給出了動力學(xué)模擬中10幀水構(gòu)象450個水分子的偶極矩分布,計(jì)算結(jié)果分別見表3和圖4.

表3列出了元胞中45個水分子10個構(gòu)象的平均偶極矩、均方根偏差及樣本中單個水分子偶極矩的最大、最小值.可以看到,在Mulliken電荷限制條件(2)和(4)下,靜電勢和靜電場擬合電荷均給出較小的分子偶極矩.偶極矩減小的這種變化趨勢,使得其更趨向于固定電荷的SPC剛性水的偶極矩值,7.84× 10-30C·m[31],以及第一性原理Mulliken電荷的偶極矩值,7.71×10-30C·m,而與體相水的實(shí)驗(yàn)觀測值, 8.34×10-30C·m[32],有一定的偏差.我們認(rèn)為導(dǎo)致計(jì)算與實(shí)驗(yàn)觀測差異的主要原因是本文中水的幾何結(jié)構(gòu)為SPC剛性水構(gòu)型,同時在第一性原理計(jì)算時只進(jìn)行了單點(diǎn)能計(jì)算,沒有進(jìn)行幾何結(jié)構(gòu)的優(yōu)化.

需要說明的是由10幀構(gòu)象450個水的偶極矩的標(biāo)準(zhǔn)偏差來看,與上節(jié)討論的氧、氫原子的電荷分布一致,不考慮電荷轉(zhuǎn)移體系的(限制條件(4))偶極矩,其漲落最小,這一點(diǎn)可以從圖3偶極矩分布圖中清晰看到.有趣的是,在Mulliken電荷限制條件下,盡管靜電場擬合電荷與Ewald加和法靜電勢擬合電荷存在一定的差異,但偶極矩的分布,包括峰寬和峰強(qiáng),卻呈現(xiàn)出很好的一致性.由于偶極矩作用決定了許多體系的靜電性質(zhì),這表明兩種方法擬合的電荷所產(chǎn)生的靜電性質(zhì)的相似性.

表3 不同擬合方法獲得電荷計(jì)算的45個水分子體系中水分子的平均偶極矩Table 3 Average dipole moment(μ)of the 45 water system calculated by the derived charges from different fitting methods

3 結(jié) 論

針對凝聚相體系,提出了一種利用第一性原理靜電勢獲取原子部分電荷的快捷方法.該方法把周期性邊界條件引入的平均靜電勢φmean作為一個擬合參數(shù),通過對第一性原理靜電勢與截?cái)嗉雍头ɑ駿wald加和法點(diǎn)電荷靜電勢進(jìn)行數(shù)值擬合,進(jìn)而獲得周期體系力場中原子的部分電荷.該方法相對擬合誤差與相應(yīng)的靜電場擬合電荷所給出的相對誤差要小得多.特別地,Ewald加和法的靜電勢擬合相對誤差僅為2%-3%,比靜電場相對誤差小一個數(shù)量級.此外,與我們期望的一致,Ewald加和法擬合的φmean隨體系幾何構(gòu)象、尺寸的選取變化不大.考慮到分子間的靜電相互作用主要與靜電勢有關(guān),并且靜電勢擬合的計(jì)算效率要遠(yuǎn)高于靜電場擬合,我們認(rèn)為Ewald加和法是一種通過第一性原理計(jì)算獲取周期性邊界條件下原子電荷的有效、快捷方法.

在不明顯增加擬合相對誤差的前提下,我們討論了四種電荷限制條件下的靜電勢和靜電場電荷擬合情況.結(jié)果表明,無論是靜電勢還是靜電場擬合,引入Mulliken電荷限制的擬合電荷,其平均值較小,分布較窄.另一方面,在只限制元胞凈電荷為零時(限制條件(1)),分子間不僅存在電荷極化還存在電荷轉(zhuǎn)移.其中,電荷轉(zhuǎn)移的大小與擬合的方法存在較大的依賴性.靜電勢截?cái)嗉雍头〝M合的電荷,分子平均的電荷轉(zhuǎn)移最大,為0.410e,而Ewald加和法電荷轉(zhuǎn)移最小,僅為0.084e.為了避免電荷轉(zhuǎn)移被人為地夸大,我們還考察了加入分子電中性的兩個電荷限制條件的擬合(限制條件(3)和(4)).需要指出的是,在分子電中性限制下,分子中原子部分電荷的分布范圍明顯增大,顯然,這種擬合電荷呈現(xiàn)了較強(qiáng)的電荷極化作用.

利用獲得的水體系中原子的部分電荷,我們計(jì)算了水分子的平均偶極矩,與第一性原理計(jì)算的偶極矩比較,引入Mulliken電荷限制的體系呈現(xiàn)了更好的趨勢.此外,利用第一性原理快速獲得體系的原子部分電荷,也將為EEM方法中對于原子價(jià)態(tài)電負(fù)性、價(jià)態(tài)硬度參數(shù)的擬合提供可靠的訓(xùn)練集.

1 Verstraelen,T.;Speybroeck,V.V.;Waroquier,M.J.Chem.Phys., 2009,131:044127

2 Mulliken,R.S.J.Chem.Phys.,1955,23:1833

3 Bader,R.F.W.;Matta,C.F.J.Phys.Chem.A,2004,108:8385

4 Breneman,C.M.;Wiberg,K.B.J.Comput.Chem.,1990,11:361

5 Arroyo,S.T.;Martin,J.A.S.;Carcia,A.H.Chem.Phys.Lett., 2002,357:279

6 Besler,B.H.;Merz,K.M.;Kollman Jr.,P.A.J.Comput.Chem., 1990,11:431

7 Singh,U.C.;Kollman,P.A.J.Comput.Chem.,1984,5:129

8 Wang,J.;Wolf,R.M.;Caklwell,J.W.;Kollman,P.A.;Case,D. A.J.Comput.Chem.,2004,25:1157

9 Whitten,A.E.;Mckinnon,J.J.;Spackman,M.A.J.Comput. Chem.,2006,27:1063

10 Wang,J.;Cieplak,P.;Kollman,P.A.J.Comput.Chem.,2000,21: 1049

11 Case,D.A.;Pearlman,D.A.;Caldwell,J.W.;Cheatham III,T.E.; Wang,J.;Ross,W.S.;Simmerling,C.;Darden,T.;Merz,K.M.; Stanton,R.V.;Cheng,A.;Vincent,J.J.;Crowley,M.;Tsui,V.; Gohlke,H.;Radmer,R.;Duan,Y.;Pitera,J.;Massova,I.;Seibel,G. L.;Singh,U.C.;Weiner,P.;Kollman,P.A.AMBER 7 user′s Manual.California:University of California,2002

12 Gao,J.;Xia,X.Science,1992,258:631

13 Field,M.J.;Bash,P.A.;Karplus,M.J.Comput.Chem.,1990,11: 700

14 Patel,S.;Brooks,C.L.J.Comput.Chem.,2004,25:1

15 Varekova,R.S.;Koca,J.J.Comput.Chem.,2006,27:396

16 Sanderson,R.T.Chemical bond and bond energies.New York: Academic Press,1976

17 Sanderson,R.T.Polar covalence.New York:Academic Press, 1983

18 Berendsen,H.J.C.;Grigera,J.R.;Straatsma,T.P.J.Phys.Chem., 1987,91:6269

19 Jorgensen,W.L.;Madura,J.D.J.Am.Chem.Soc.,1983,105: 1407

20 Mahoney,M.W.;Jorgensen,W.L.J.Chem.Phys.,2000,112: 8910

21 Yang,Z.Z.;Wu,Y.;Zhao,D.X.J.Chem.Phys.,2004,120:2541

22 Zhang,Q.;Yang,Z.Z.Chem.Phys.Lett.,2005,403:242

23 Yang,Z.Z.;Zhang,Q.J.Comput.Chem.,2006,27:1

24 Yang,Z.Z.;Qian,P.J.Chem.Phys.,2006,125:064311

25 Wu,Y.;Yang,Z.Z.J.Phys.Chem.A,2004,108:7563

26 Ewald,P.P.Ann.Phys.,1921,64:253

27 Dovesi,R.;Saunders,V.R.;Roetti,C.;Orlando,R.;Zicovich-Wilson,C.M.;Pascale,F.;Civalleri,B.;Doll,K.;Harrison,N.M.; Bush,I.J.;D′Arco,P.;Llunell,M.Crystal 06 user′s manual. Torino:University of Torino,2006

28 Spoel,D.v.d.;Lindahl,E.;Hess,B.;Groenhof,G.;Mark,A.E.; Berendsen,H.J.C.J.Comput.Chem.,2005,26:1701

29 Saunders,V.R.;Freyria-Fava,C.;Dovesi,R.;Salasco,L.;Roetti, C.Mol.Phys.,1992,77:629

30 Stewart,R.F.God Jugosl Cent Kristalogr,1982,17:1

31 Campa?á,C.;Mussard,B.;Woo,T.K.J.Chem.Theory Comput., 2009,5:2866

32 Shirono,K.;Daiguji,H.Chem.Phys.Lett.,2006,417:251

Atomic Partial Charges for Periodic Systems from First-Principles Calculations

LI Ya-Na Lü Yang ZHOU Li-Chuan CHEN Li LI Shen-Min*
(Liaoning Key Laboratory of Bio-Organic Chemistry,Dalian University,Dalian 116622,Liaoning Province,P.R.China)

We calculated the electrostatic potential(ESP)and electric field(EF)of periodic liquid water systems using the quantum chemistry software package,Crystal.We propose a method to obtain atomic partial charges rapidly for periodic systems based on first-principles calculations.In this method,the average electrostatic potential φmean, which is introduced to meet the periodic boundary condition,is taken as a parameter during the least squares fitting of the ESP from first-principles calculations and used in the Ewald summation.A comparison of the two methods,i.e., ESP and EF fitting,reveals that the relative root mean-square deviation(RMS)of the former method is only 2%-3%, which is one order of magnitude smaller than that of the latter method.In addition,the distribution of the derived atomic partial charges and dipole moments for the water system are discussed using four charge restrained fits.

Atomic partial charge; Least square fitting; Electrostatic potential; Electrostatic field; Ewald summation

O641

Received:April 24,2010;Revised:June 21,2010;Published on Web:August 27,2010.

*Corresponding author.Email:shenmin@dl.cn;Tel:+86-411-87402384.

The project was supported by the National Natural Science Foundation of China(20973027,20633050)and Program for Liaoning Excellent Talents in University,China(2007R02).

國家自然科學(xué)基金(20973027,20633050)和遼寧省高等學(xué)校優(yōu)秀人才支持計(jì)劃項(xiàng)目(2007R02)資助

?Editorial office of Acta Physico-Chimica Sinica