使用Multiwfn考察周期性體系的芳香性
使用Multiwfn考察周期性體系的芳香性
Using Multiwfn to study aromaticity for periodic systems
文/Sobereva@北京科音 2024-Jul-31
0 前言
衡量化學體系的芳香性的方法非常多，見《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）。研究芳香性的文章數目甚巨，但大多數都是對分子、原子團簇這樣的孤立體系研究的，主要在于很少有程序能夠支持將衡量芳香性的方法用于周期性體系。如今Multiwfn已經可以將很多方法用于周期性體系了，包括多中心鍵級、AV1245、AVmin、PDI、FLU、FLU-pi、PLR、HOMA、Bird、ELF二分值、Shannon芳香性指數、芳香環的環臨界點屬性，等等，還可以通過繪制LOL-pi函數圖像直觀考察共軛情況。在本文中，將對其中大部分方法在Multiwfn中的操作進行演示。作為例子的體系是一個共價有機框架化合物（COF）的一層，此體系在《使用Multiwfn對周期性體系做鍵級分析和NAdO分析考察成鍵特征》（http://www.shanxitv.org/719）中已經作為例子分析過，結構如下所示，本文要通過不同方法對比綠色虛線里三個環的芳香性的大小。可見1號環是C3N3，原子序號為1,2,21,22,11,12。2號環和3號環相當于萘片段中的兩個六元環，2號環序號為19,20,50,49,18,48，3號環序號為14,48,18,44,46,16。這里序號是按照成鍵關系排的。
筆者假定讀者已經閱讀過《衡量芳香性的方法以及在Multiwfn中的計算》充分了解了各種衡量芳香性方法的特征，也假定讀者有了Multiwfn的基本使用常識，不了解者建議閱讀《Multiwfn入門tips》（http://www.shanxitv.org/167）和《Multiwfn FAQ》（http://www.shanxitv.org/452）。Multiwfn可以從官網http://www.shanxitv.org/multiwfn免費下載，注意務必使用2024-Jul-1及以后更新的版本（注意看Multiwfn啟動時的更新日期的提示），否則情況和本文所述不符。
本文涉及的多數芳香性分析方法都是基于波函數的，上面COF這個體系的由CP2K程序在PBE/DZVP-MOLOPT-SR-GTH級別計算得到的molden格式的波函數文件和http://www.shanxitv.org/719里用的是一致的，請從此文的鏈接里下載，是其中的COF-MOS-1_0.molden。此文件的產生方法在此文里也詳細說了，不會用CP2K者必須仔細看。本文涉及到的諸如HOMA這樣的芳香性指數是基于幾何結構的，只要給Multiwfn提供的文件里包含結構信息就可以，如cif、pdb、xyz、gjf、mol2等常見格式都可以。如果被分析的區域是跨晶胞的，則輸入文件必須能給Multiwfn提供晶胞信息，哪些格式能提供晶胞信息見《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）里的說明。
2 計算多中心鍵級
多中心鍵級是非常好、很嚴格的考察芳香性的指標。這一節計算一下當前體系中三個六元環的六中心鍵級，以此考察它們六中心共軛的強弱，這正比于它們的芳香性。多中心鍵級的概念參看《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）和《使用AdNDP方法以及ELF/LOL、多中心鍵級研究多中心鍵》（http://www.shanxitv.org/138）。
啟動Multiwfn，載入COF-MOS-1_0.molden，然后輸入
25 //離域性與芳香性分析
1 //多中心鍵級
1,2,21,22,11,12 //計算1號環。原子序號按原子連接關系輸入
此時得到多中心鍵級結果為0.0399。類似地，再分別輸入19,20,50,49,18,48 [回車]和14,48,18,44,46,16 [回車]計算2、3號環，多中心鍵級分別為0.0212和0.0281。
根據六中心鍵級大小可見此體系中的C3N3環（1號環）的芳香性顯著強于六元碳環，而不與C3N3環連接的六元碳環（3號環）的芳香性比與C3N3環連接的六元碳環（2號環）更強。
3 計算AV1245和AVmin
AV1245的概念和計算方法在《使用Multiwfn計算AV1245指數研究大環的芳香性》（http://www.shanxitv.org/519）中介紹過，這里就不多說了。這里用AV1245算一下前述的三個環的芳香性。
啟動Multiwfn，載入COF-MOS-1_0.molden，然后輸入
25 //離域性與芳香性分析
2 //AV1245
1,2,21,22,11,12 //計算1號環。原子序號按連接關系輸入
看到以下結果
AV1245 times 1000 for the selected atoms is 7.35835761
AVmin times 1000 for the selected atoms is 7.358354 ( 11 12 2 21)
即AV1245乘以1000和AVmin乘以1000都為7.358。
再輸入19,20,50,49,18,48計算2號環，結果為
AV1245 times 1000 for the selected atoms is 4.46380673
AVmin times 1000 for the selected atoms is 3.441219 ( 20 50 18 48)
再輸入14,48,18,44,46,16計算3號環，結果為
AV1245 times 1000 for the selected atoms is 6.30980750
AVmin times 1000 for the selected atoms is 4.014836 ( 46 16 48 18)
可見，無論是從AV1245還是AVmin上來看，芳香性都是1號環>3號環>2號環，這和多中心鍵級的結論完全一致。
4 計算PDI
這一節用PDI衡量芳香性。注意PDI只能用于六元環。Multiwfn中PDI是基于模糊空間定義的離域化指數算的，對周期性體系默認用的是Hirshfeld原子空間劃分。
啟動Multiwfn，載入COF-MOS-1_0.molden，然后輸入
15 //模糊空間分析
5 //PDI
[回車] //用默認的格點間距計算原子重疊矩陣（AOM）。對較大體系可以用更大的比如0.35 Bohr格點間距明顯節約這一步的計算時間，對結果影響甚微
1,2,21,22,11,12 //第一個環里的原子序號
結果如下
Delocalization index of 1(C ) -- 22(N ): 0.104520
Delocalization index of 2(N ) -- 11(C ): 0.104520
Delocalization index of 21(C ) -- 12(N ): 0.104520
PDI value is 0.104520
再依次輸入19,20,50,49,18,48 [回車]和14,48,18,44,46,16 [回車]分別計算第2、3個環的PDI，結果分別為0.083658和0.096507。PDI越大芳香性越強，可見PDI給出的芳香性大小的結論也是1號環>3號環>2號環。
5 計算FLU和FLU-pi
先計算FLU芳香性指數。啟動Multiwfn，載入COF-MOS-1_0.molden，然后輸入
15 //模糊空間分析
6 //FLU
[回車] //用默認的格點間距計算原子重疊矩陣（AOM）
現在從屏幕上可以看到做FLU計算用的參數。輸入第一個環里的原子序號1,2,21,22,11,12，得到如下結果
Atom pair Contribution DI
1(C ) -- 2(N ): 0.000800 1.474298
2(N ) -- 21(C ): 0.000312 1.508760
21(C ) -- 22(N ): 0.000800 1.474298
22(N ) -- 11(C ): 0.000312 1.508761
11(C ) -- 12(N ): 0.000800 1.474298
12(N ) -- 1(C ): 0.000312 1.508761
FLU value is 0.003335
再依次輸入19,20,50,49,18,48 [回車]和14,48,18,44,46,16 [回車]分別計算第2、3個環的FLU，結果分別為0.017882和0.012199。由于FLU越小芳香性越強，因此FLU的芳香性大小的結論是1號環>3號環>2號環，和前述芳香性指標結論相同。
下面再來計算FLU-pi。算這個需要告訴Multiwfn當前體系的所有pi占據軌道序號。可以按照《使用Multiwfn觀看分子軌道》（http://www.shanxitv.org/269）肉眼一個個看來記錄序號，但太麻煩。建議按照《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）介紹的做法自動指認。啟動Multiwfn，載入COF-MOS-1_0.molden，然后輸入
100 //其它功能（Part 1）
22 //自動檢測pi軌道
0 //當前軌道是離域的（分子/晶體軌道屬于這類）
馬上Multiwfn就顯示了占據的pi軌道序號，為47,50,60,61,64,70-73,77,81,82,84-89,92-94。下面計算FLU-pi，接著輸入
0 //對識別出的pi軌道什么都不做
0 //返回主菜單
15 //模糊原子空間分析
7 //FLU-pi
[回車] //用默認的格點間距計算原子重疊矩陣（AOM）
47,50,60,61,64,70-73,77,81,82,84-89,92-94 //pi占據軌道的序號
之后輸入第一個環里的原子序號1,2,21,22,11,12，得到如下結果
Average of DI-pi is 0.387267
Atom pair Contribution DI
1(C ) -- 2(N ): 0.000191 0.375470
2(N ) -- 21(C ): 0.000191 0.399063
21(C ) -- 22(N ): 0.000191 0.375471
22(N ) -- 11(C ): 0.000191 0.399063
11(C ) -- 12(N ): 0.000191 0.375470
12(N ) -- 1(C ): 0.000191 0.399063
FLU-pi value is 0.001148
再依次輸入19,20,50,49,18,48 [回車]和14,48,18,44,46,16 [回車]分別計算第2、3個環的FLU-pi，結果分別為0.028167和0.026396。由于FLU-pi越小芳香性越強，因此FLU-pi的芳香性大小的結論是1號環>3號環>2號環，和前述芳香性指標結論相同。
6 計算HOMA
HOMA是流行的基于環上的鍵長特征衡量芳香性的方法。給Multiwfn提供的輸入文件里有結構信息就行了。由于COF-MOS-1_0.molden里也包含結構信息，所以此例還是用這個作為輸入文件。啟動Multiwfn，載入COF-MOS-1_0.molden，然后輸入
25 //芳香性分析
6 //HOMA和Bird芳香性指數
0 //開始計算HOMA（基于默認參數）
現在屏幕上顯示了計算HOMA用的參數。輸入第一個環里的原子序號1,2,21,22,11,12，得到如下結果
Atom pair Contribution Bond length(Angstrom)
1(C ) -- 2(N ): -0.003592 1.349180
2(N ) -- 21(C ): -0.001191 1.342742
21(C ) -- 22(N ): -0.003592 1.349181
22(N ) -- 11(C ): -0.001191 1.342741
11(C ) -- 12(N ): -0.003592 1.349181
12(N ) -- 1(C ): -0.001191 1.342741
HOMA value is 0.985651
可見HOMA為0.985651，并且環上各個鍵的鍵長以及對HOMA的貢獻量都給出了。再依次輸入19,20,50,49,18,48 [回車]和14,48,18,44,46,16 [回車]分別計算第2、3個環的HOMA，結果分別為0.598160和0.781579。由于HOMA越接近1芳香性越強，因此HOMA的芳香性大小的結論是1號環>3號環>2號環，和前述芳香性指標結論相同。
Bird芳香性指數的計算方法和HOMA非常類似，只不過是在主功能25里的子功能6里選擇2而非此例的0而已，這里就不演示了。
7 計算Shannon芳香性指數
Shannon芳香性指數是基于被考察的環上的各個鍵的鍵臨界點的電子密度定義的，因此計算它之前必須先用Multiwfn搜索臨界點。這方面的操作在《使用Multiwfn結合CP2K做周期性體系的atom-in-molecules (AIM)拓撲分析》（http://www.shanxitv.org/717）里有詳細說明，這里就不對細節做具體介紹了。
啟動Multiwfn，載入COF-MOS-1_0.molden，然后輸入
2 //拓撲分析
3 //用每一對原子連線的中點作為初猜點搜索臨界點
0 //在圖形窗口里查看結果
在圖形界面里要求只顯示(3,-1)臨界點，也就是鍵臨界點，并且要求把臨界點序號顯示出來，看到下圖。可見第一個環上的臨界點序號為72,76,67,59,55,61。
點擊圖形界面右上角的Return按鈕，然后輸入
20 //計算Shannon芳香性
72,76,67,59,55,61 //第1個環上的BCP序號
馬上看到結果：
Electron density at CP 72: 0.3352748688 Local entropy: 0.2993061920
Electron density at CP 76: 0.3318277479 Local entropy: 0.2979424307
Electron density at CP 67: 0.3352750629 Local entropy: 0.2993062683
Electron density at CP 59: 0.3318276569 Local entropy: 0.2979423945
Electron density at CP 55: 0.3352752224 Local entropy: 0.2993063310
Electron density at CP 61: 0.3318279246 Local entropy: 0.2979425011
Total electron density: 2.0013084835
Total Shannon entropy: 1.7917461175
Expected maximum Shannon entropy: 1.7917594692
Shannon aromaticity index: 0.0000133518
即Shannon芳香性指數為0.0000133518，每個鍵臨界點的電子密度和局部熵也都給出了，這些是計算Shannon芳香性指數的中間量。
類似地也輸入第2、3個環上的鍵臨界點的序號，分別為62,64,54,45,43,51和64,71,82,84,75,65，結果分別為0.0013759152和0.0010309007。由于Shannon芳香性指數越小芳香性越強，因此芳香性大小的結論是1號環>3號環>2號環，和前述芳香性指標結論相同。
8 基于環臨界點屬性判斷芳香性
在某個環中央區域的環臨界點位置上，若垂直于環的方向上電子密度曲率越負，說明這個環的芳香性越強。這種判斷芳香性的方法需要利用AIM拓撲分析得到的環臨界點的屬性。為了用這種分析，和上一節一樣先用主功能2做拓撲分析，然后在選項0里只顯示(3,+1)臨界點，即環臨界點，看到的圖如下。可見第1、2、3號環的環臨界點序號分別為66、53、73。
關閉圖形窗口后輸入
21 //計算順著某個方向的電子密度梯度和曲率
66 //第1個環的環臨界點序號
1 //指定某個方向
0,0,1 //當前體系平行于XY平面，因此計算Z方向的電子密度的梯度和曲率
結果如下
Electron density is 0.0287379670 a.u.
Electron density gradient is 0.0000000000 a.u.
Electron density curvature is -0.0261250099 a.u.
可見，在第一個環的環臨界點位置上垂直于這個環的電子密度曲率為-0.0261250099 a.u.。類似地再使用這個功能考察第2、3個環分別對應的第53、73號環臨界點垂直于環平面的電子密度曲率，結果分別為-0.0154770974和-0.0159388407 a.u.。由于數值越負芳香性越強，因此芳香性大小的結論是1號環>3號環>2號環，和前述芳香性指標結論相同。
9 ELF-pi二分點考察芳香性
這一節通過各個環上ELF-pi函數的等值面首次發生二分的ELF-pi數值（ELF-pi二分值）來衡量芳香性的大小。這個值既可以通過反復微調等值面數值來獲得，也可以通過對ELF-pi做拓撲分析來獲得，前者更直觀，后者更精確，相關信息見《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）。這一節演示前者的做法。
首先需要產生ELF-pi的格點數據。啟動Multiwfn，載入COF-MOS-1_0.molden，然后輸入
100 //其它功能（Part 1）
22 //檢測pi軌道
0 //軌道都是離域形式
2 //設置其它軌道占據數為0
0 //返回主菜單
5 //計算格點數據
9 //ELF
9 //利用晶胞平移矢量定義格點信息
[回車] //坐標原點為(0,0,0)
[回車] //盒子三個方向尺寸和相應晶胞邊長一致
0.15 //用較精細的0.15 Bohr格點間距，使得通過觀看等值面獲得二分值能盡可能準確
2 //將格點數據導出為cube文件
現在當前目錄下就有了ELF.cub，對應ELF-pi的格點數據。用VMD顯示其等值面（不會操作的話用《在VMD里將cube文件瞬間繪制成效果極佳的等值面圖的方法》http://www.shanxitv.org/483里提供的腳本）。選Display - Orthographic用正交視角，在Graphics - Representation里將等值面數值（isovalue）由小到大調節，會發現在等值面數值為0.761的時候正好1號環上等值面首次發生二分，位置如下圖箭頭所示。因此1號環的二分點數值為0.761。
由上圖可見，在1號環ELF-pi等值面二分的時候，2、3號環早就發生了二分、等值面間間隔已經很大了，這說明1號環的芳香性遠強于2、3號環。
再將等值面數值調節到0.623，會看到2、3號環上首次出現了ELF-pi的二分，如下圖紅色箭頭所示。這個鍵是被兩個六元環共享的，怎么區分哪個芳香性更強？這需要再看其它的鍵。在下圖紫色箭頭所示的2號環上的位置，只要等值面數值再增加一點點就會二分，而3號環就沒有這個情況，說明3號環的芳香性比2號環更強、環上的電子共軛遇到的瓶頸更少。因此ELF-pi的分析結論和其它方法完全一致。
10 LOL-pi圖形分析
作為前述分析的擴展和補充，這一節繪制LOL-pi圖直觀展示一下單層COF體系的pi共軛情況，這能十分清楚地讓大家理解三個環上電子離域特征的差異。這種分析在前述的《在Multiwfn中單獨考察pi電子結構特征》文中有大量介紹，在筆者的Theor. Chem. Acc, 139, 25 (2020)中也有很多例子，歡迎閱讀和引用。
啟動Multiwfn并載入COF-MOS-1_0.molden后，依次輸入
100 //其它功能（Part 1）
22 //檢測pi軌道
0 //軌道都是離域形式
2 //設置其它軌道占據數為0
0 //返回
4 //繪制平面圖
10 //LOL
1 //填色圖
[回車] //用默認的格點數
1 //XY平面
4.25 //當前體系中的原子都在Z=3.25位置，繪制平面上方1 Bohr處的圖像顯然應該設Z=4.25 Bohr
利用后處理菜單對作圖效果的一番調節，得到下圖。如果不會調的話，仔細把Multiwfn手冊4.4節的所有例子都仔細領會了，結合后處理菜單中提示得很清楚的選項名字去理解，自然就明白了。當前用的色彩刻度是0到0.65，超過上限的區域顯示為白色。
由上圖可見，在1號環上pi共軛很顯著，LOL-pi在環上分布得比較均勻。而如圖上我標注的箭頭所示，2號環上有三個鍵的pi共軛顯著低于其它地方，因此這個環上的六中心pi共軛作用明顯受到了很大遏制。而3號環上的pi共軛整體相對好點，但也不是很均勻，而是在不同鍵上有大有小，所以芳香性只是介于1和2號環之間。
11 總結
Multiwfn支持極為豐富的考察化學體系芳香性的方法，本文對其中已支持周期性體系的大部分方法的操作過程進行了簡明扼要的演示。從結果可見盡管這些方法思想差異很大，但對于當前研究的這個體系中的三個六元環，它們給出的芳香性順序完全一致，互相印證了彼此的可靠性。當然對于一些特殊情況，由于一些方法原理上的局限性、穩健性和普適性的不足也有可能導致結果存在差異。本文示例的只是一個很標準、簡單的體系，希望大家能充分舉一反三，將Multiwfn應用于廣泛體系的芳香性的研究中。
使用Multiwfn做任何分析在發表文章時都請務必記得按照程序啟動時的提示恰當引用Multiwfn原文，給別人代算時也必須明確告知對方這一點。