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

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

(大連大學,遼寧省生物有機化學重點實驗室,遼寧大連 116622)

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

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

(大連大學,遼寧省生物有機化學重點實驗室,遼寧大連 116622)

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

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

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

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

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

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

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

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

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

1.1 第一性原理計算靜電勢與靜電場

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

1.2 周期體系的點電荷靜電勢和靜電場

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

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

上式中φmean是一個與位置無關的常數項,φ′(r)是由點電荷模型定義的周期體系中元胞內位置r處的靜電勢:

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

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

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

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

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

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

上式中,A0(ri)是第一性原理計算的網格點ri處的靜電勢或靜電場,Np是元胞內網格點數.χ2由式(7)給出:

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

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

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

2 結果與討論

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結 論

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

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

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

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).

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

?Editorial office of Acta Physico-Chimica Sinica