使用Multiwfn做電子密度、ELF、靜電勢、密度差等函數的盆分析
使用Multiwfn做電子密度、ELF、靜電勢、密度差等函數的盆分析
文/Sobereva @北京科音
First release: 2013-Mar-27 Last update: 2023-Jun-16
前言:Multiwfn(http://www.shanxitv.org/multiwfn)主功能17對應于盆分析功能。本文目的在于簡要介紹盆分析的概念、算法,并結合幾個應用實例介紹盆分析功能的基本操作方法并展現盆分析的基本用途。更多的信息參見Multiwfn手冊3.20節和4.17節。筆者在《Multiwfn支持的分析化學鍵的方法一覽》(http://www.shanxitv.org/471)中也專門介紹了盆分析在研究電子結構方面能起到什么價值。不了解Multiwfn的讀者建議閱讀《Multiwfn FAQ》(http://www.shanxitv.org/452),不知道這些分析用到的波函數文件怎么產生的話看《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》(http://www.shanxitv.org/379)。
本文撰寫后Multiwfn又經歷了巨量更新,因此盆分析的數據、極值點序號會和本文例子所示的有些許不同,讀者請以實際在Multiwfn里看到的結果為準。
1 盆、吸引子和盆分析的基本概念
盆(basin)這個概念最早是由Bader在它的Atoms in molecules (AIM)理論中引入的,在AIM理論中盆分析是對于電子密度來進行的。隨后這個概念被Savin、Silvi推廣到用于分析電子定域化函數(ELF)。關于ELF的介紹可參見《ELF綜述和重要文獻小合集》(http://bbs.keinsci.com/thread-2100-1-1.html),里面有大量綜述以及筆者寫的相關博文和論文。實際上,盆這個概念對于任何實空間函數都適用,比如電子密度拉普拉斯函數、動能密度、靜電勢、電子密度差、福井函數等等。
盆是整個空間中的局部空間,而所有的盆所占空間總和就是整個空間。每個盆彼此間通過函數的零通量面(zero-flux surface)來劃分,因而零通量面也叫做盆間面(interbasin surface)。具體來說,零通量面上的各個點上都滿足函數的梯度在垂直于這個面方向上的分量為0的條件。每個盆當中唯一地包含一個“吸引子”(attractor),它是這個盆空間內函數的極大點。吸引子和拓撲分析術語中的(3,-3)型臨界點是一回事。
例如,下圖是電子密度的盆的二維示意圖,這種圖經常出現在AIM分析中
棕色的圓點就是吸引子。由于電子密度的極大點總是很接近原子核(對于H可能偏得多一點。另外還有所謂的非核吸引子,但這里不提及),所以這些吸引子的位置基本就在原子核上。灰色細線是電子密度的梯度線,從各個方向向著吸引子匯集。粗的藍線就是零通量面,將整個分子空間分割為一個個盆。這些盆和原子一一對應,所以AIM理論中將電子密度的盆稱為原子盆。這樣的盆被AIM理論認為描述了原子在分子中所占的空間,所以也被叫做AIM原子空間。上面的圖只是盆在分子平面上的二維示意圖,氫原子盆的三維圖如下所示
注:雖然眾多梯度線組成的面上所有點也都符合“垂直于這個面上函數的梯度分量為0”的條件,但是這樣的面穿越吸引子,所以不特稱為零通量面或盆間面。
將盆劃分好之后,就可以在這些盆的空間里進行各種分析。通常有四類分析,它們都被Multiwfn的盆分析模塊所支持:
(1)在盆里對指定的實空間函數進行積分。例如,積分電子密度,就得到了盆內電子的布居數。如果這個盆是AIM盆,那么得到的就是原子的電子布居數,令原子核電荷減去它就得到了AIM原子電荷。如果這個盆是ELF,且比如對應共價鍵,那么可以通過盆中的電子密度積分值討論有多少電子與共價鍵的構成相關。還可以積分諸如自旋密度、能量密度等等。用于定義盆的實空間函數和在盆中被積分的實空間函數顯然是可以不同的。Multiwfn可以讓任意實空間函數在任意實空間函數定義的盆中被積分。
(2)計算盆內的電多極矩(electric multipole moments),具體來講是指單極矩、偶極矩和四極矩。這對于定量分析電子密度的分布十分有益。比如可以研究孤對電子對應的ELF盆對分子偶極矩的貢獻(如J. Comp. Chem., 29, 1440)、考察形成復合物過程中電子密度被極化的過程、研究盆內電子分布在哪個方向較彌散、偏離球對稱分布程度的多少等等。
(3)計算盆內的定域化指數(localization index, LI)和盆間的離域化指數(delocalization index, DI)。前者衡量平均來說有多少電子定域在盆里;后者衡量的是平均有多少電子離域在兩個盆之間,或者說盆之間共享了多少電子。如果是AIM原子盆,那么它們之間的DI和共價鍵級是可以相對應的。
(4)計算原子對盆布居數的貢獻。比如對應C-N之間成鍵的ELF盆,這個盆有一定電子布居數,我們可以把這個ELF盆空間再用AIM方式進行劃分,來得到C和N各自對這個盆布居數的貢獻。
(2)、(3)、(4)方式的分析本質上也是需要在盆內進行積分才能實現的。所以,要想做盆分析,必須得通過某種算法將盆空間表示出來,并且能夠在其中對特定實空間函數進行積分。這是下一節要討論的。
2 生成盆的算法
由于AIM理論的重要性,好幾十年前就開始有人不斷提出各種各樣的積分AIM原子盆的方法。這些方法都是解析方法,例如J.Comp.Chem.,21,1040和J.Phys.Chem.A,115,13169。遺憾的是,由于盆的邊界復雜,以解析方式積分相當困難,非常耗時,且積分公式都很復雜,不容易實現。這些方法一般需要先搜索到電子密度的(3,-3)及(3,-1)臨界點,然后利用(3,-1)臨界點產生盆間面,然后以(3,-3)臨界點為中心向外發散地做盆內的積分。這類解析方式的積分不僅慢、復雜、通常需要事先找臨界點,而且有一個重要缺點就是基本上只能用于積分電子密度的盆,不適合用于其它實空間函數如ELF。這類方法的唯一好處就是積分精度很不錯。經典的AIM程序如AIMPAC、AIM2000、AIMALL、Morphy之類用的都是這類解析積分方法。值得一提的是在J.Comp.Chem.,30,1082中作者提出一種利用類似Becke的DFT泛函積分的方法積分盆,雖然積分精度有所降低,但是積分速度快了很多,且算法復雜度降了不少。但是這種方法依然只適合于積分電子密度的盆。
Silvi他們組的專用于ELF盆分析的TopMod程序用的是一種基于立方格點的數值方式的盆積分方法。在此程序發表后,逐漸有一些其它的基于立方格點(實際上也完全可以用矩形格點)的盆積分方法和相應程序被提出。這類基于立方格點的數值積分方法速度比較快,最大的好處是可以用于任意類型的實空間函數,而缺點是積分精度遜于解析積分方法。或者說,用的格點數據的格點間距越小則積分精度越高,但是要想達到解析方法那么高的精度,那么計算格點數據的時間以及儲存格點數據所需內存量將甚巨而難以實現。不過,只要對積分精度要求不苛刻,那么這類格點積分的精度一般也夠了。值得一提的是這類方法產生盆之前不需要事先尋找臨界點,吸引子會在產生盆的過程中自動獲得,所以很方便。由于Multiwfn的盆分析模塊強調的重點之一是普適性強,所以采用的是立方格點積分的方法,而非解析方式的積分盆。這類積分方法中目前最好的是near-grid方法(J.Phys.:Condens.Matter,21,084204),這也是Multiwfn默認的方法,而它的前身on-grid方法(Comput.Mat.Sci.,36,354)在Multiwfn中也支持。On-grid方法生成盆的速度更快,但明顯看得出盆間面的劃分準確度不如near-grid方法,因此積分值準確度也差一些。
near-grid方法得到吸引子位置、產生并積分盆的大體過程是:
首先提供要分析的實空間函數的格點數據。然后,依次從每個格點出發,按照梯度最大的方向沿著格點不斷進行移動從而產生出一條軌跡。如果走到一個格點且發現它的數值比相鄰的點都大,那么這個格點就會被當成吸引子,并且這條軌跡上所有點將被指認(assign)到這個吸引子。如果從某個點開始沿著梯度最大方向不斷移動的過程中走到了一個已經被指認了吸引子的點,那么這條軌跡中的所有點也將被指認到這個吸引子。將所有格點循環一遍后,被指認到同一個吸引子的所有點就構成了這個吸引子對應的盆。盆內任意一實空間函數f的積分值可通過此式很容易地得到:∑f(r_i)*dv。其中加和符號循環盆內每個格點,編號以i表示,r_i即是格點的坐標。dv是空間微元,即每個小格子的體積。
積分不同的實空間函數對格點精度要求不同。比如電子密度拉普拉斯函數、源函數(source function)在原子核附近波動極大,所以就很不適合用立方格點方法來積分。即便用了很高精度的格點,仍會發現在盆內積分拉普拉斯函數的結果偏離理想值0比較多。另外,不同區域的盆積分難度也不同。例如電子密度在原子核附近有尖峰,而在原子價層區域變化則較緩,所以在核電子對應的ELF盆里積分電子密度的精度就相對低一些,明顯不如積分價層ELF盆那樣容易得到較高精度。
為了解決AIM盆積分時由于電子密度、源函數、電子密度拉普拉斯函數、能量密度等常見被積函數在較重原子內層區域變化較快而難以準確積分的問題,筆者提出了將均勻格點和原子中心格點相混合的積分方法并在Multiwfn中予以實現。此方法不僅實現容易,精度也高,徹底解決了near-grid方法在AIM盆積分上的軟肋。
3 現有的支持盆分析的程序
那些經典的利用解析方法做電子密度盆分析的AIM程序前面已經提到過了,即AIMPAC、AIM2000、AIMALL、Morphy等,它們沒法用于對其它實空間函數進行盆分析。Bader程序(http://theory.cm.utexas.edu/bader/)是極少數基于立方格點方法進行電子密度盆分析的程序,操作還算方便,但是功能很弱,僅能計算AIM原子盆內的電子密度積分值。還有的程序如InteGriTy,雖然載入的是電子密度格點數據,但是實際上還是用傳統解析方法進行盆積分,只不過任意處的電子密度不是直接計算而是通過格點數據插值得來。
普適型的基于立方格點方式做盆分析的程序最出名的是TopMod(http://www.lct.jussieu.fr/pagesperso/silvi/silvi_english.html),但是它實際上只支持電子密度和ELF的盆分析,也只能在這兩類盆里計算LI/DI和積分電子密度。TopMod的操作超級繁瑣,做個盆分析需要調用多個子程序(top_grid、top_bas、top_pop),得敲一大堆命令。而且程序諸多方面都顯得莫名其妙,頗不好用,而且想可視化盆也非常困難。Multiwfn的盆分析功能的最初開發目標就是完勝TopMod。TopGrid(http://www.lct.jussieu.fr/pagesperso/pilme/topchempage.html)是與TopMod關系甚密的另一個基于立方格點的盆分析程序,只能讀入格點數據來做盆分析,生成盆的算法和TopMod有異,普適性號稱比它要廣,可以支持盆的電多極矩分析。但這程序甚至比TopMod還難用,而且計算結果筆者發現不準,比如觀看它產生的盆會明顯發現盆的形狀和理想形狀往往偏離不小。
相對于上述程序,Multiwfn的盆分析功能無論是普適性、計算速度、積分精度、靈活性還是操作便利性上都有質的飛躍。Multiwfn可以對任意實空間函數劃分盆,對任意實空間函數在盆內進行積分,還可以計算盆的電多極矩、LI和DI,還可以計算原子對盆布居數的貢獻。實空間函數既可以直接由Multiwfn計算也可以從外部文件讀取。分析結果能夠直接可視化而不必再借助于第三方程序。同時,吸引子和盆都可以導出到外部文件,以便使用諸如VMD等可視化程序顯示來得到更好的效果。由于Multiwfn本身代碼的高效,以及作了充分的并行化,做盆分析速度很快,包括C60這樣較大的體系的盆分析也不用花多少時間。均勻格點+原子中心格點積分更是Multiwfn獨有的算法,積分AIM盆精度完勝TopMod和Bader程序。Multiwfn盆分析模塊操作極其方便,只需要按幾個鍵,敲幾下回車即可完成,自動化程度極高。而且輸出信息簡明易懂,各個功能、涉及到的原理在手冊里都能方便地查閱到。可以說,Multiwfn標志著實空間函數的盆分析進入了新紀元。
4 實例
以下六個實例涉及到盆分析模塊的大部分操作,做過一遍后就能對盆分析有個基本的認識。//后面的是注釋。
4.1 HCN分子的電子密度盆(AIM盆)分析
4.1.1 生成盆
啟動Multiwfn并輸入
examples\HCN.wfn
17 //盆分析
1 //生成盆并尋找吸引子
1 //讓Multiwfn計算電子密度格點數據,用來得到它的盆和吸引子
2 //中等質量格點,也就是格點間距為0.1Bohr。格點數據涵蓋的空間范圍程序將根據體系電子密度值自動確定,列表中還有其它很多定義格點設定的方式,這里暫時先不用
此時Multiwfn開始計算格點數據,然后生成盆并確定吸引子的位置。之后吸引子的坐標以及吸引子位置上的函數值會顯示在屏幕上
Attractor X,Y,Z coordinate (Angstrom) Value
1 -0.02645886 -0.02645886 -1.52058663 0.33959944
2 0.02645886 0.02645886 -0.51514986 49.48609717
3 -0.02645886 -0.02645886 0.64904009 76.55832377
由于氫的密度低,很顯然吸引子1就是與氫相對應的。如果之后還想重新看吸引子的信息,以及更詳細的信息,可以選擇選項-3。注意如果開啟了Multiwfn的并行模式(默認是開啟的),每次尋找出的吸引子的編號順序都可能不同,以實際為準,后同。
4.1.2 觀看盆和吸引子
現在選0來可視化剛找出來的吸引子和生成的盆。在彈出來的圖形界面中,紫色的文字對應吸引子的標號。每個吸引子用綠色小球表示,但由于它們對于此例在原子核處,所以得取消選擇show molecule選項才能看到它們。圖形界面中各個控件的含義都很明確就不多提了,試試就明白了。界面右下角(如果是Linux版則在右上角)是盆列表,一個盆對應一個吸引子,盆的編號和吸引子編號是一致的,選擇哪個盆就會顯示出哪個盆。這里我們在盆列表中選1來讓1號盆顯示,并且選show basin interior讓盆內部也顯示,此時看到的圖如下
位于盆邊界的格點通過綠球顯示,使得盆的輪廓被勾勒出來。如果不想顯示盆,在盆列表里選None即可。
在很多AIM分析的文章中,盆只顯示了處于電子密度0.001 a.u.等值面內部區域的部分,如果想只繪制這部分就在圖形界面上方選擇Set basin drawing method - rho>0.001 region only,然后就可看到下面的圖像
在盆列表里還會看到"Unas"(意指Unassigned)和"Boun"(意指travelled to box boundary)。前者指的是產生盆過程中由于某些數值算法原因沒有成功被指認所屬吸引子的格點,后者是指產生盆的過程中按照梯度最大方向移動但是最終移動到格點數據范圍邊界的點。這兩類格點一般沒有什么物理意義,對分析結果影響甚微,通常出現在離體系較遠的地方。值得一提的是,如果為了節省計算時間,沒有讓格點數據范圍囊括整個分子而是只囊括了局部區域,那么可以得到只位于局部區域的盆,與此同時"Boun"類型的點會大量出現,依然無視它即可。在盆列表里選擇"Unas"和"Bound"就可以分別看到這兩類點分布在哪。對于當前體系,在產生盆之后從屏幕上可以看到
The number of unassigned grids: 0
The number of grids travelled to box boundary: 0
也就是說"Unas"和"boun"類型的格點數都為0,這是最普通的情況。也因此對于當前例子在盆里表里選了"Unas"或"Boun"之后什么都不會顯示。
現在點擊RETURN按鈕關閉圖形界面。
順帶一提,如果你希望獲得漂亮的原子盆圖像,并且還能把臨界點和鍵徑也顯示出來,最好的辦法是按照這個視頻的方法操作:《使用Multiwfn和VMD繪制原子盆(AIM盆)》(https://www.bilibili.com/video/av85202089)。主要流程就是用選項-5把指定的盆導出成各個cube文件,然后拖到VMD程序里繪制為等值面。整個過程非常簡單,而圖像效果非常好,例如下圖顯示了丙烯醛的所有鍵徑+臨界點,并且將非氫原子的原子盆都用不同顏色顯示了出來。產生盆時選擇的格點質量越高,盆的邊緣就會越平滑。
4.1.3 積分盆
在我們生成了盆以后,就可以通過盆分析界面里的選項2~5對盆進行各種分析了。先做哪種分析都無所謂。這里先選"2 Integrate real space functions in the basins"做盆積分。我們這里要在盆里積分的是電子密度,所以接下來可以選1,來讓Multiwfn計算各個盆內的格點上的電子密度值以得到盆內的密度積分值。但是,由于前面已經生成過了一遍電子密度格點數據,這套數據仍然在內存里存著,因此就沒必要讓Multiwfn重新算每個點的密度了,直接調用內存里的格點數據的值就行了,這大大節省了時間。所以我們選"0 The values of the grid data stored in memory",來讓內存中存著的格點數據作為被積的函數。立刻結果就輸出出來了
#Basin Integral(a.u.) Volume(a.u.^3)
1 0.7356142812 441.70000000
2 5.3511710147 566.31900000
3 7.9020060745 830.25600000
Sum of above values: 13.98879137
除了每個盆的積分值還有盆的體積。不過這里輸出的盆體積對于當前體系沒什么意義,因為三個盆的總和就是整個格點數據空間范圍,故格點數據空間范圍取得越大盆的體積就越大。輸出中可以看到積分精度不算很理想,總值13.989和理想值14有一定偏離。這主要是因為在前文中提到的,電子密度在重原子的原子核附近衰減速度較快,用均勻格點對這部分沒法積分得很準。為了解決此問題,Multiwfn對于AIM盆積分提供了混合格點積分的選項,即選項7。我們輸入
7 //混合格點進行積分
1 //原子中心+均勻格點,忽略邊界修正
1 //電子密度
結果為
Total result:
#Basin Integral(a.u.) Vol(Bohr^3) Vol(rho>0.001)
1 0.7356461193 441.700 34.884
2 5.3552270393 566.319 102.832
3 7.9086504198 830.256 134.228
Sum of above integrals: 13.99952358
Sum of basin volumes (rho>0.001): 271.944 Bohr^3
這回總積分值13.9995已經十分接近于理想值14了,盆體積Vol(Bohr^3)依然沒什么物理意義。接著還輸出了每個原子的AIM電荷,以及原子體積。
Normalization factor of the integral of electron density is 0.999966
The atomic charges after normalization and atomic volumes:
1 (C ) Charge: 0.644591 Volume: 102.832 Bohr^3
2 (N ) Charge: -0.908920 Volume: 134.228 Bohr^3
3 (H ) Charge: 0.264329 Volume: 34.884 Bohr^3
這里的Volume等同于上面的Vol(rho>0.001),代表的是盆中電子密度大于0.001的體積,按照Bader對范德華表面的定義,這相當于分子范德華表面內原子所占的體積,是有確切物理意義的,可認為是原子的尺寸。
注意,盡管得到的總積分值已經幾乎等同于理想值了,但實際上原子盆的積分值并非很精確,對應的AIM電荷也不是非常精確。為了得到更準確的積分結果,還得對盆邊界進一步修正,這對應于選項7中的2或3,這相比剛才用的1精度有所提升,但會更耗時得多,這三種模式的精度、耗時、內存消耗量在屏幕上已經提示了。這里我們選7再選2然后選1,看看經過邊界修正后的積分結果,結果為
Total result:
#Basin Integral(a.u.) Vol(Bohr^3) Vol(rho>0.001)
1 0.7450375217 442.296 35.480
2 5.2508474961 563.647 100.160
3 8.0036553438 832.332 136.304
Sum of above integrals: 13.99954036
Sum of basin volumes (rho>0.001): 271.944 Bohr^3
Normalization factor of the integral of electron density is 0.999967
The atomic charges after normalization and atomic volumes:
1 (C ) Charge: 0.748980 Volume: 100.160 Bohr^3
2 (N ) Charge: -1.003918 Volume: 136.304 Bohr^3
3 (H ) Charge: 0.254938 Volume: 35.480 Bohr^3
可見現在和剛才得到的AIM電荷有明顯不同,現在得到的是比較準確的結果。如果想進一步提升精度,可以在生成盆的時候用更好的格點,比如High quality grid,乃至Lunatic quality grid。然而,對于大體系,High quality grid已經很耗時、很占內存了。
注意,雖然用選項7比起選項2積分AIM盆精確得多,但是選項7只能積分AIM盆,而選項2基于均勻格點的積分可以用于任何函數的盆。
提示:上面討論得比較多,通常來說,獲得AIM電荷最一般也較準確的做法是進入盆分析模塊,然后選
1 //生成盆
1 //用電子密度劃分盆
2 //中等質量格點
7 //用混合格點積分AIM盆
2 //對邊界進行精確修正
1 //積分電子密度
4.1.4 計算盆的電多極矩
電多極矩的計算公式在手冊3.18.3節里有詳細說明,只不過那些公式里的原子核坐標對于盆分析來說應當被替換為吸引子坐標,而積分范圍應當由模糊空間替換為盆空間。
我們選8來做盆的多極矩計算。注意計算前最好先用選項7中的選項2在精確修正盆邊界的情況下進行一次盆積分(我們上一節已經做過了),因為這樣做盆的邊界會永久性修正,當前計算盆多極矩的時候結果也更準了。另外值得一提的是用選項3也能計算盆的多極矩,但那是基于均勻格點的,雖然對各種類型盆都能用但對AIM盆精度較低,所以這里不用。用選項8得到的碳的多極矩結果如下
Result of atom 1 (C )
Basin electric monopole moment: -5.250847
Basin electric dipole moment:
X= 0.000003 Y= 0.000002 Z= 1.095848 Magnitude= 1.095848
Basin electron contribution to molecular dipole moment:
X= 0.000003 Y= 0.000002 Z= 6.026238 Magnitude= 6.026238
Basin electric quadrupole moment (Cartesian form):
QXX= -0.759559 QXY= -0.000000 QXZ= -0.000002
QYX= -0.000000 QYY= -0.759564 QYZ= -0.000003
QZX= -0.000002 QZY= -0.000003 QZZ= 1.519122
The magnitude of electric quadrupole moment (Cartesian form): 1.519122
Electric quadrupole moments (Spherical harmonic form):
Q_2,0 = 1.519122 Q_2,-1= -0.000003 Q_2,1= -0.000002
Q_2,-2= -0.000000 Q_2,2 = 0.000003
Magnitude: |Q_2|= 1.519122
盆的單極矩就是盆內電子密度積分值的負值。此盆的偶極矩Z分量為明顯的正值,表明盆內的電子密度主要分布于相對于2號吸引子(碳原子核)位置的Z值更負的區域。四極矩的ZZ分量為正,而XX、YY分量為負,表明盆內的電子密度分布相對于球對稱分布來說,在Z方向(分子軸方向)有所延展,而在其余方向有所收縮。四極矩的大小表現盆內電子密度偏離球對稱的程度,此例的數值1.519不小,說明偏離球對稱分布還是挺大的,這顯然是因為形成了化學鍵所引起的。可見,通過盆的電多極矩分析,我們可以以定量的方式明確地討論盆內電子的分布狀態。實際上計算諸如盆內的八極矩、十六極矩也很容易做到,但是用處不是很大,所以Multiwfn最高給到四極矩。
4.1.5 計算定域化和離域化指數
最后,基于AIM盆,我們分析定域化指數(LI)和離域化指數(DI)。選擇功能4,Multiwfn就會輸出盆之間的DI矩陣和盆內的LI值。如下所示
************ Total delocalization index matrix ************
1 2 3
1 0.97290584 0.90037924 0.07252660
2 0.90037924 3.49303265 2.59265341
3 0.07252660 2.59265341 2.66518001
Total localization index:
1: 0.259 2: 3.496 3: 6.658
比如DI(2,3)=2.593,表明2號盆(C的AIM原子空間)和3號盆(N的AIM原子空間)之間平均來說共享了2.593個電子,這和HCN中C與N的形式鍵級3.0挺接近。實際上,DI就是一種共價鍵級的衡量標準,隨著鍵的極性的增強DI值一般會降低。DI矩陣的對角元是相應的列(或行)的矩陣元的加和,如果把DI就當成鍵級看待,那么對于閉殼層體系,DI矩陣的對角元便是原子價。因此,可以說此體系N原子的原子價是2.66。
1號盆(H的AIM原子空間)的定域化指數為0.259,說明平均有0.259個電子定域在這個盆內,這個數值遠小于此盆內電子布居數,表明在HCN體系內,電子很容易從H的AIM空間中離域出去,外界的電子也容易離域進來,換句話說,盆內外的電子容易相互交換。
為了便于考察,免得用戶需要自行對照盆的序號和原子序號,Multiwfn在屏幕上還同時輸出了行和列序號對應于原子序號的DI矩陣和LI值。
4.2 乙炔的ELF盆分析
啟動Multiwfn并輸入
examples\C2H2.wfn
17
1 //生成盆并尋找吸引子
9 //ELF
2 //中等質量格點
0 //觀看結果。圖像如下所示
理論上,在兩個碳之間有一個環形的ELF吸引子,這是由于C-C之間存在的兩個pi鍵所致。在Multiwfn中,吸引子都是以點來描述的。由于數值算法的原因,如圖所見,Multiwfn會在這個環的位置上找出來一大堆數值差不多的吸引子,一起圍成一個環形。顯然不能讓這一圈吸引子一個個地孤立著,否則物理意義就喪失了,所以Multiwfn會自動地將挨得近的且數值相仿佛的吸引子進行“歸簇”(clustering)。這一圈吸引子也因此就被歸為了一個簡并的吸引子(degenerated attractor),即圖中所示的2號吸引子,原本那一圈吸引子就被稱為2號吸引子的“成員吸引子”。而2號盆的空間也就是其成員吸引子對應的盆的總和。有少數時候,在默認設定下Multiwfn自動歸簇的結果可能與預期的有差距,比如這一圈吸引子沒有被歸在一起而被歸成了兩個吸引子,這時應當在生成盆之前先用功能-6調節閾值,包括相對數值差異和距離間隔這兩個條件閾值。如果任意兩個吸引子之間的相對數值差異小于閾值,而且它們的距離也小于閾值的話,則它們會被歸到一起。距離閾值是和格點間距有關的,假設定義的時候輸入了一個值n,那么距離閾值就是n*sqrt(dx^2+dy^2+dz^2),其中dx、dy和dz是三個方向的格點間距。如果將n設為0,那么Multiwfn就不會在尋找完吸引子之后對它們進行歸簇。另外,在生成盆并搜索完吸引子之后,也可以用功能-6的選項3來自行指定將哪些吸引子歸到一起,后面的例子將會用到這個功能。
根據ELF的符號法則,2號盆應該被指認為V(C1,C3),這里V代表這個盆的電子是來自于C1和C3的價層電子。類似地,1號和5號盆應當被分別指認為V(C3,H4)和V(C1,H2)。3號和4號盆應當被分別指認為C(C3)和C(C1),C代表Core,表明內核電子主要都在這個些盆里。在圖形界面中取消Show molecule復選框不讓分子結構顯示出來就能看到這種盆了,如下所示
接下來計算電子密度在ELF盆中的積分值。關閉圖形窗口后選擇功能2,然后選擇選項1(對應電子密度),就會看到如下結果接下來計算電子密度在ELF盆中的積分值。關閉圖形窗口后選擇功能2,然后選擇選項1(對應電子密度),就會看到如下結果
#Basin Integral(a.u.) Volume(a.u.^3)
1 2.2159821115 768.39800000
2 5.3670486307 972.82100000
3 2.0949016050 0.83200000
4 2.0949016050 0.83200000
5 2.2159823917 768.53100000
Sum of above values: 13.98881634
C(C1)和C(C3)的積分值為2.09,很接近2,這也正對應了碳有兩個內核電子。V(C3,H4)和V(C1,H2)的積分值為2.21,也接近2,對應于C-H鍵共享一對兒電子這個經典化學概念。值得一提的是,盡管此C(C1)和C(C3)是對稱的、V(C3,H4)和V(C1,H2)是對稱的,因此原則上它們的積分值應該相同,但實際上因為格點精度是有限,所以積分結果有一些數值差異,使用更高的格點精度能使結果與分子對稱性滿足得更好。按照經典化學鍵理論,乙炔的兩個碳之間應該享有三對兒電子,共6個電子,但是V(C1,C3)的積分值5.37與6偏差略大。不過也沒有理由要求這種基于實空間的分析的結果與經典化學概念完全一致。
然后選功能3再選-1,來得到每個盆的電多極矩分析結果。其中的4號盆,即C(C1)的結果為
Basin 4
Basin electric monopole moment: -2.094902
Basin electric dipole moment:
X= -0.104745 Y= -0.104745 Z= 0.020728 Magnitude= 0.149575
Basin electron contribution to molecular dipole moment:
X= -0.000000 Y= 0.000000 Z= -2.388409 Magnitude= 2.388409
Basin electric quadrupole moment (Cartesian form):
QXX= -0.003034 QXY= -0.007856 QXZ= 0.001555
QYX= -0.007856 QYY= -0.003034 QYZ= 0.001555
QZX= 0.001555 QZY= 0.001555 QZZ= 0.006067
The magnitude of electric quadrupole moment (Cartesian form): 0.006067
Electric quadrupole moments (Spherical harmonic form):
Q_2,0 = 0.006067 Q_2,-1= 0.001795 Q_2,1= 0.001795
Q_2,-2= -0.009071 Q_2,2 = 0.000000
Magnitude: |Q_2|= 0.011205
由于當前體系是軸對稱的,而且軸在Z軸上,所以盆內的偶極矩的X和Y分量都應當為0。不過由于這次只用了medium quality grid,所以這兩個分量不是很接近于0。但即便如此,當前的結果還是有定性意義的,起碼看到偶極矩的大小非常小,而且四極矩的大小也很接近0,這說明原子的內核電子受分子環境影響甚小,基本接近于球對稱分布。
輸入0返回上一級菜單,然后選功能4計算LI和DI,結果如下
************ Total delocalization index matrix ************
1 2 3 4 5
1 1.31231709 1.03496562 0.16085935 0.02291204 0.09358008
2 1.03496562 2.71090549 0.32048710 0.32048709 1.03496568
3 0.16085935 0.32048710 0.51131492 0.00705643 0.02291205
4 0.02291204 0.32048709 0.00705643 0.51131492 0.16085936
5 0.09358008 1.03496568 0.02291205 0.16085936 1.31231716
Total localization index:
1: 1.560 2: 4.011 3: 1.834 4: 1.834 5: 1.560
DI(3,4)對應于C(C1)和C(C3)之間的離域化指數,其數值0.007十分接近于0,這表明在兩個原子的內核區域之間電子很難相互離域。DI(1,3)=0.160和DI(2,3)=0.320雖然不大,但是還是有一定數值的,這說明C(C3)盆內的電子與和它相鄰的僅有的兩個盆V(C3,H4)和V(C1,C3)內的電子相互交換還是有一定幾率的。DI(2,1)和DI(2,5)都接近于1,表明電子在C-C鍵區域和C-H鍵區域間離域很容易。從定域化指數上看,碳的內核的盆的LI為1.83,接近于盆內的電子布居數2.09,表明內核電子整體來說喜歡呆在內核區域,而不容易離域到外頭去,而相應地外界電子也難以進來。而V(C3,H4)的LI為4.011,和它的電子布居數5.37差得挺大,表明V(C3,H4)內的電子容易與外界電子發生交換,這也是價層ELF盆的普遍情況。
4.3 水分子的靜電勢盆分析
雖然很早以前已經有人對靜電勢做過拓撲分析,但是對靜電勢進行盆分析的研究還相當少。
靜電勢同時存在原子核電荷主導的正值的區域,以及電子所主導的負值的區域。對于這種正負值同時存在的實空間函數,Multiwfn的盆分析模塊在生成盆并搜索吸引子之前,會首先將負值區域的數值乘以-1來讓符號顛倒,然后照常生成盆和尋找吸引子,弄完了之后再恢復原先的負值區域的數值。這樣一來,負值區域也會找到吸引子并能產生對應的盆,但這些吸引子實際上是“排斥子”,即函數值的極小點。在Multiwfn程序中以及下文中,不管本質上是吸引子還是排斥子,一律都被叫做吸引子。
啟動Multiwfn并輸入
examples\H2O.fch //當然,也可以用.wfn或.wfx文件作為輸入
17
1 //生成盆并搜索吸引子
12 //計算靜電勢格點數據
2 //中等質量格點
0 //查看結果
從顯示出來的圖中可以看到1、2、3號吸引子對應于原子核電荷產生的靜電勢極大值點。而4、5號吸引子對應于氧的孤對電子導致的靜電勢負值區域中的極小點。諸如這樣的負值區域的極小點以及相應的盆都會用藍色表示,以便與綠色表示的正值區域的極大點和盆進行區分。從盆列表里選擇4號盆,可以看到如下圖像,
注意如命令行窗口中的提示所顯示的,在此例中生成盆的時候有1627個格點最后跑到了格點數據范圍的邊緣,在盆列表里選擇"Boun"就能看到它們,如下圖所示。顯然這些點里分子都很遠,所以這個問題可以無視。
盆分析模塊提供了便利的計算吸引子、原子核間的幾何參數的功能。進入功能-2后可以看到操作方法介紹。例如想測量4號吸引子與1號原子(氧原子)間的距離,就輸入c4 a1,結果為2.238 Bohr。再比如,測量“吸引子4--氧原子核--吸引子5”的角度,就輸入c4 a1 c5,結果為80.84度,這在某種程度上也可以看做是氧的孤對電子的夾角(孤對電子夾角還可以利用ELF、LOL來研究,結果與此有一定差異,見《使用Multiwfn做拓撲分析以及計算孤對電子角度》(http://www.shanxitv.org/108)。然后輸入q退到上一級菜單。
假設我們想手動把兩個孤對電子對應的吸引子合并為一個吸引子,可以輸入
-6 //設定歸簇參數或自行對吸引子進行歸簇
3 //將指定的吸引子歸到一起
4,5 //2和3號吸引子將歸為一個,同時盆也合并到一起
0 //返回
這時再選0打開圖形界面,會看到吸引子的序號都變了,兩個孤對電子對應的吸引子合并成了一個,序號為4。
選功能2,然后選1,來得到電子密度在靜電勢盆中的積分值,結果為
#Basin Integral(a.u.) Volume(a.u.^3)
1 0.8647176869 362.33900000
2 0.8647176862 362.27200000
3 7.5147940814 25.00000000
4 0.7405367228 628.13800000
Sum of above values: 9.98476618
Integral of the grids travelled to box boundary: 0.00000004
雖然氧有兩個孤對電子,但是在兩個孤對電子對應的靜電勢盆,即目前的4號盆當中電子密度積分值僅有0.740。這是因為對任何分子,大部分離分子較近的空間總是由核電荷所主導,故而靜電勢為正。靜電勢為負值的區域都處在分子的較外側,這些區域電子密度是比較有限的,所以電子密度積分值在負的靜電勢盆中不會很大。結果也顯示,那些"Boun"類型的格點內的積分值幾乎為0,這也是為什么前面說這些格點可以忽略。
對靜電勢做盆分析需要獲得靜電勢格點數據,而計算靜電勢格點數據對于稍大的體系都很花時間。如果你想顯著降低這個步驟上的耗時,特別是對于大體系、大基組而言,強烈建議通過令Multiwfn自動借用Gaussian里的cubegen工具來實現(因為cubegen比Multiwfn內部代碼在靜電勢上快不少),只需把settings.ini里的cubegenpath參數設為cubegen的實際路徑,之后Multiwfn里不管哪個分析功能,凡是牽扯到基于波函數計算靜電勢格點數據的時候都會自動調用cubegen來算。詳見http://www.shanxitv.org/435里的說明。
4.4 水分子的電子密度差的盆分析
在進入盆分析之前,首先要利用主功能5把電子密度差的格點數據算出來。電子密度差這里具體指分子密度與組成它的各個原子在孤立狀態下的密度差。這需要涉及到生成原子波函數文件,為了省事,此例就直接用Multiwfn自帶的原子波函數文件。做法也就是把examples目錄下的atomwfn目錄挪到上一級目錄下(假設你是通過雙擊Multiwfn圖標來運行Multiwfn的,否則就將atomwfn放到你調用Multiwfn時所在的目錄下)。關于用Multiwfn做電子密度差的具體內容,見《使用Multiwfn作電子密度差圖》(http://www.shanxitv.org/113)
移動過atomwfn目錄后,啟動Multiwfn,然后輸入
examples\H2O.fch
5 //計算格點數據
-2 //產生變形屬性
1 //電子密度
3 //高質量格點。注意這里的“高質量格點”和前面用到的盆分析模塊里的高質量格點不是一回事。這里定義的是總格點數而非格點間距,故而高質量只是對于較小分子來說的。由于密度差比電子密度、ELF變化都要復雜,為了盆分析可靠,建議用較好的格點質量。
0 //返回主菜單。剛才計算的格點數據會一直保存在內存里。
17 //盆分析
1 //生成盆并尋找吸引子
2 //基于保存在內存里的格點數據來產生盆和尋找吸引子
選擇0進入圖形界面觀看電子密度差的吸引子和相應的盆。顯示和不顯示分子結構時的圖分別如下面兩圖所示
電子密度差的正值和負值部分分別代表原子形成分子后,電子密度增加和減少的區域。上圖中綠球對應正值部分的極大值,藍球對應負值部分的極小值。如果搞不清楚為什么極值點這么分布,對照著電子密度差截面圖就會弄明白。水分子垂直于分子平面和分子平面上的電子密度差圖如下所示
對照平面圖,就很容易理解3號和6號盆對應的是形成O-H鍵造成的電子密度增加區域,而4號和5號盆對應的是氧原子以sp3雜化后由于孤對電子的出現造成的密度增加的區域。
值得一提的是由于格點精度用得仍然不夠高,導致7號吸引子的另一側的等價吸引子沒有被找出來。另外,8號吸引子偏離了理想位置。不過,這樣的吸引子和盆都不是化學上感興趣的,它們對應的是密度差變化的十分細微的結構,沒有太大物理意義,所以可以忽視這些問題、
現在看看O-H鍵形成后造成了多少電子往鍵中間發生了挪動,也就是對3號或6號盆進行積分,被積函數也是電子密度差。于是,進入功能2,然后輸入0,來將內存中保存著的電子密度差格點數據當做被積函數。結果立刻顯示出來,3號和6號盆的積分值都是0.09。而對應于孤對電子的形成造成密度增加的4號或5號盆的積分值是0.21。
如果想要將等值面圖和吸引子進行對比,可以選擇-10退回Multiwfn主菜單,然后進入主功能13,再選-2,當前內存中的格點數據的等值面圖就會顯示出來,之前找出來的吸引子也會同時顯示。下圖是isovalue=0.05的圖
PS:關于保存在內存里的格點數據這里多說幾句。用.cub/.grd文件作為輸入文件的話,其中的格點數據就會一直保存在內存里。另外用主功能5或者13計算或者對格點數據進行操作的話,得到的格點數據也會保存在內存里。當有新的格點數據生成時,舊的格點數據就會被覆蓋。若在盆分析功能里的功能1當中讓Multiwfn計算新的格點數據,那么原先內存里的格點數據就也被覆蓋。
4.5 對乙烷進行基于AIM盆的源函數分析
源函數在Multiwfn手冊2.6節有簡要介紹。在The Quantum Theory of Atoms in Molecules-From Solid State to DNA and Drug Design一書的7.6章有個簡明的綜述建議讀讀。這個函數依賴于兩個坐標,其中一個坐標是參考點。當源函數被用來討論成鍵(包括弱相互作用)問題時,參考點一般取在BCP(鍵臨界點)上。經常通過在AIM原子空間內對源函數進行積分,以研究原子對BCP處密度的貢獻。我們這一節要分析乙烷的各部分對C-H鍵的BCP密度的貢獻值。我們先要進行拓撲分析找出臨界點位置。啟動Multiwfn,依次輸入
examples\ethane.wfn
2 // 拓撲分析
2 // 搜索核臨界點
3 // 搜索BCP
0 // 觀看結果
臨界點分布和編號如下所示
這里我們打算將CP11這個C-H鍵的BCP取做參考點。從文本窗口能找到它的坐標(0.000000000.-1.199262548,-1.909104063),將坐標直接從屏幕上拷貝到系統的剪切板中(如果不知道怎么去做這步,見手冊5.4節)。然后我們看看這個點的電子密度是多少,這一步不是必須的,但是由于對源函數進行全空間積分的結果原理上會精確等于參考點的電子密度值,所以我們現在獲知了它的電子密度,就便于稍后檢驗積分得是否準確。選擇7,然后輸入臨界點編號11,就能得到它的密度0.276277。
現在來設定參考點。一種方法是用settings.ini里的refx、refy和refz參數來設參考點的X,Y,Z坐標,但是這樣設完了還得重啟一下Multiwfn才能生效,這里我們用另一種方法來設。先選-10從拓撲分析界面返回到主菜單,然后輸入1000,就進入了一個秘密界面,選擇1,再把CP13的坐標粘貼進去并按回車,參考點坐標就設好了。
然后,依次輸入以下命令生成AIM盆
17 // 盆分析
1 // 生成盆并尋找吸引子
1 // 電子密度
2 // 中等質量格點
進入功能0,會看到吸引子如下圖分布,我們要考察的是下方的甲基對CP11處密度的貢獻,盆編號可見為1,2,3,4
接下來要對盆進行積分,輸入
7 //用混合格點方法積分AIM盆
1
19 //源函數
結果如下
Atom Basin Integral(a.u.) Vol(Bohr^3) Vol(rho>0.001)
1 (C ) 5 0.00361681 148.640 70.392
2 (H ) 8 0.00289447 465.226 50.058
3 (H ) 6 0.00280055 426.841 50.097
4 (H ) 7 0.00280276 418.343 50.095
5 (C ) 2 0.12224052 149.871 70.396
6 (H ) 1 0.12040142 471.817 50.054
7 (H ) 4 0.01051821 424.105 50.095
8 (H ) 3 0.01051757 432.711 50.097
Sum of above integrals: 0.27579230
Sum of basin volumes (rho>0.001): 441.284 Bohr^3
積分值總和0.27579230與參考點的密度0.276277十分接近,誤差完全在可接受范圍內。下方甲基對應的1~4號盆數值加和為0.12040142+0.12224052+0.01051757+0.01051821=0.2637,占CP11密度的0.2637/0.2763*100%=95.4%,因此甲基是它自身的C-H鍵BCP密度的最主要貢獻者,而其余原子只貢獻4.6%。有人發現這個結論對于不同長度的烷烴都是十分一致的,這被用來說明化學基團的可移植性問題。
4.6 計算成鍵的兩原子對其ELF鍵盆布居數的貢獻
此例我們計算CH3NH2體系中C、N原子對它們之間V(C,N) ELF鍵盆布居數的貢獻,使用的是AIM方式劃分原子空間。用相同的方法還可以計算原子對其它任意類型盆的布居數的貢獻。啟動Multiwfn并輸入
我們先生成記錄ELF盆定義的basin.cub文件,此文件中每個點的數值對應這個點所屬的ELF盆編號。
examples\CH3NH2.wfn
17
1
9 //ELF
2 //中等質量格點
然后進入選項0,檢驗ELF盆編號
我們看到V(C,N)盆的編號是5號。關閉圖形窗口,輸入
-5 //輸出盆
a //將盆信息輸出到當前目錄下的basin.cub
然后照常生成AIM盆,輸入
1 //重新生成盆
1 //重新選擇函數
1 //電子密度
9 //我們必須保證新生成的格點數據的格點設定與basin.cub嚴格一致,最好的做法就是這里選9,直接從basin.cub中讀取格點設定信息
basin.cub
0 //檢查吸引子編號
我們看到N和C對應的吸引子編號分別為2和3。然后關閉圖形窗口輸入
9 //計算原子對外部盆(basin.cub中記錄的盆)的布居數貢獻,選完之后會讀取當前目錄下的basin.cub
2 //要考察的原子對應的吸引子的編號,2對應N
5 //V(C,N) ELF盆對應的編號
結果為1.15866,即N對V(C,N)盆貢獻了1.159個電子。然后輸入
3 //對應C
5 //V(C,N) ELF盆對應的編號
結果為0.46359,這即是C對V(C,N)盆貢獻的電子數,明顯少于N貢獻的,因此C-N鍵可以認為有很顯著的極性。
5 總結&其它
本文介紹了盆分析的意義、算法,結合很多例子對Multiwfn的盆分析操作進行了講解。對于最一般的盆分析,在Multiwfn里其實就這么幾步,簡單得很:
[文件名] // 一般用.mwfn/.wfn/.wfx/.fch/.molden文件
17 // 盆分析模塊
1 // 生成盆并尋找吸引子
1 // 選擇用什么實空間函數定義盆和吸引子,這里假設為電子密度
2 // 產生什么質量的格點數據,這里是中等質量的,適合用作預覽
0 // 觀看結果
然后根據要分析的內容,選擇功能2~9即可。
加入盆分析模塊的Multiwfn 3.0的推出,對于實空間函數的盆分析領域,是具有革新意義的。Multiwfn在普適性、速度、靈活度、簡便易用等各種方面相對以往程序都有著質的飛躍,Multiwfn將ELF盆分析的門檻降低到量化初學者也能輕松完成的程度。本文雖然用的例子都是小分子,但用于較大的分子也沒有問題。除非是分子太大,這樣的話可能會出現內存容量不足以儲存格點數據的問題。對于這種情況,如果只對局部區域感興趣,建議只將格點數據空間范圍定義在相應的局部區域,來減少總格點數從而避免內存不足。