使用Multiwfn計算AV1245指數研究大環的芳香性
2020-Sep-14后記：本文中關于多中心鍵級計算耗時的說法已經完全過時了。因為從2020-Sep-12更新的Multiwfn開始，程序中實現了筆者提出的一種超快的多中心鍵級計算算法，而精度沒有絲毫損失，對于哪怕含有幾十個原子的環也可以非常快速地計算多中心鍵級，這使得多中心鍵級的應用范疇有了革命性的擴展。相比之下，AV1245的實際價值就已經不大了。
使用Multiwfn計算AV1245指數研究大環的芳香性
文/Sobereva@北京科音
First release: 2019-Oct-27 Last update: 2022-Apr-29
0 前言
Multiwfn程序（http://www.shanxitv.org/multiwfn）支持大量分析芳香性的方法，見《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）。其中提到了一個重要指標是多中心鍵級，這個指標對于研究多中心鍵和芳香性非常方便、可靠，如今非常流行。多中心鍵級的一個主要缺點是計算耗時隨考慮的中心數呈指數型增長，雖然算6、7個原子瞬間就能得到結果，但計算12個原子以上的環卻幾乎不可能。
在Phys. Chem. Chem. Phys., 18, 11839 (2016)中提出了AV1245指數，可以視為是多中心鍵級的近似，原文實測發現其數值與多中心鍵級的相關性不錯。由于AV1245的耗時只是隨環中的原子數增加以線性方式增長，因此哪怕對于包含十幾、乃至幾十個原子的大環，AV1245都可以快速地計算出來，這使得定量衡量大環的芳香性成為了可能。Multiwfn的最新版本已經支持了AV1245指數的計算，作為主功能9的子功能11出現。下面交代一些此方法的思想和細節，然后給出兩個具體的計算例子。
值得一提的是，在Carbon, 165, 468 (2020)一文中，筆者將此方法應用于了電子結構非常特殊的18碳環體系用于考察其pi電子全空間離域特征和芳香性，得到了很有價值的結果，充分體現出了AV1245的實際意義。另外，在《深入揭示18碳環的重要衍生物C18-(CO)n的電子結構和光學特性》（http://www.shanxitv.org/640）介紹的筆者的Chem. Eur. J., 28, e202103815 (2022)一文對C18-(CO)n的電子離域的研究中，將AV1245分解為了in-plane和out-of-plane兩類pi電子的貢獻分別討論了它們的離域性。非常建議讀者們看看這兩篇文章里面的相關討論，也推薦在使用AV1245的時候引用這兩篇文章作為應用范例。
1 AV1245的原理
原子序號的輸入順序是可能影響多中心鍵級的計算結果的。一般我們用多中心鍵級研究芳香性的時候，是按照環當中原子順序依次輸入的。還有一種少見的多中心鍵級是在計算過程中考慮所有原子序號的置換(permutation)，表達式見Multiwfn手冊3.11.2節，這種多中心鍵級下面稱為I_perm。I_perm并不適合研究電子沿著某個環的路徑的離域情況，但適合用于考察選定的一批原子彼此間的離域程度（不考慮特定的離域路徑）。
一些研究發現在共軛環當中，1-2和4-5鍵之間通常有下圖所示的這種共振
這種共振特征可以靠1,2,4,5之間的四中心I_perm捕捉到，數值越大說明共振效應越強，也因此暗示這些原子所在的環的相應區域的共軛效應比較強。I_perm和電子共享指數(ESI, electron sharing index)有直接聯系，ESI可以根據此式從I_perm換算過來：ESI=[2/(n-1)!]*I_perm，此處n是中心數。因此，四中心ESI等于四中心I_perm除以3，二者是正比關系。
AV1245將被研究的環上所有四中心ESI取平均作為衡量這個環的芳香性的指標。其思想容易理解，如果環的各個局部位置都存在顯著的上述形式的共振，即環的各個區域都有很強的共軛，那么整個環的共軛就比較強，從而可以說環的芳香性比較大。
AV1245的這種衡量環的芳香性雖然沒有多中心鍵級原理上嚴格，但根據原文的實際測試，發現對于能算得動多中心鍵級的環，AV1245能不錯地重現出多中心鍵級所能展現出的芳香性差異。對于多中心鍵級算不動的大環，經測試AV1245的結果也比較滿足一般化學直覺，應當是合理的。
AV1245絕不是研究大環芳香性的唯一方法，比如往往可以看軌道圖形根據休克爾規則考察、通過ICSS考察（見《通過Multiwfn繪制等化學屏蔽表面(ICSS)研究芳香性》http://www.shanxitv.org/216）、通過ELF-pi或LOL-pi考察（見《在Multiwfn中單獨考察pi電子結構特征》http://www.shanxitv.org/432）、通過AICD環電流圖考察（見《使用AICD 2.0繪制磁感應電流圖》http://www.shanxitv.org/294）。AV1245的一個好處是計算非常簡單和快速，而且可以給出一個簡單直白的數值便于橫向對比。
值得一提的是，原文里計算AV1245用的ESI是基于原子重疊矩陣算的，而且利用AIM方式劃分原子空間，這種方式計算相當耗時、費事，把簡單問題搞復雜。因此筆者將AV1245實現在Multiwfn中的時候沒有用原文的方法，而是以常規的基于密度矩陣和重疊矩陣的方式計算I_perm、折算成四中心ESI，再組成AV1245值。這種方式算AV1245速度很快，結果也合理。只不過由于和原文做法的差異，Multiwfn這樣給出的AV1245數值比原文里的小一些（完全不影響橫向對比），乘上0.635后和原文數值比較接近。這種做法唯一一個缺點就是基組不能帶彌散函數，否則而由于計算出的多中心鍵級不合理，最終導致AV1245的結果也不合理。Multiwfn支持的各種包含基函數信息的格式，如.fch、.molden都可以用于這種情況下計算AV1245，產生方式見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）。
有時候要算陰離子，這種情況往往要帶彌散函數。為了對這種情況也能合理且快速地計算AV1245，Multiwfn還支持基于自然原子軌道(NAO)計算AV1245。此時需要用NBO程序帶著DMNAO關鍵詞的輸出信息作為輸入文件，這樣計算的AV1245不怕彌散函數；而對于本身就沒用彌散函數的時候，這樣計算的結果和常規方式計算的幾乎一樣。
下面看兩個例子。第一個例子是菲，我們用AV1245計算它的中心的六元環和邊緣的六元環，看看芳香性的差異。這個體系用多中心鍵級也可以研究，因此我們可以對比一下AV1245和更嚴格的多中心鍵級的結論是否吻合。第二個例子是卟啉，我們將要通過AV1245研究其中四種可能的大范圍的離域路徑，看看哪種離域特征最強，因而最能體現卟啉的芳香性特征。
2 例1：分析菲的局部芳香性
本例研究的菲的原子序號如下所示
啟動Multiwfn，然后輸入
examples\phenanthrene.fch
9 //鍵級分析
11 //計算AV1245
1,2,3,4,5,6 //邊緣環上的原子序號，按照原子連接關系輸入（順時針輸入還是逆時針輸入都一樣）
結果如下
4-center electron sharing index of 1 2 4 5: 0.01304820
4-center electron sharing index of 2 3 5 6: 0.01245026
4-center electron sharing index of 3 4 6 1: 0.00801537
4-center electron sharing index of 4 5 1 2: 0.01304820
4-center electron sharing index of 5 6 2 3: 0.01245026
4-center electron sharing index of 6 1 3 4: 0.00801537
AV1245 times 1000 for the selected atoms is 11.17127674
由于環上有6個原子，所以如提示所示，程序一共算了6種四中心ESI，然后取平均并乘以1000使得數量級便于考察，最后結果為11.17。
之后我們輸入3,4,8,9,10,7來計算中間環的1000*AV1245，結果為5.01。這體現出，菲的邊緣的環的芳香性強于中間的環。
AV1245給出的這個結論是否可靠？由于六元環比較小，我們可以再用意義更嚴格的六中心鍵級檢驗一下。接著輸入
q //退出AV1245計算界面
0 //返回至主菜單
9 //鍵級分析
11 //計算AV1245
1,2,3,4,5,6 //邊緣環的原子序號
程序輸出的多中心鍵級為0.059。然后再輸入3,4,8,9,10,7，多中心鍵級為0.026，可見這也明顯體現出邊緣環的芳香性顯著強于中間的環。通過其它方法，比如AdNDP、ELF-pi，也同樣可以得到這樣的結論，見《使用AdNDP方法以及ELF/LOL、多中心鍵級研究多中心鍵》（http://www.shanxitv.org/138）。
下面再演示一下基于NAO對菲計算AV1245。要用的輸入文件是examples目錄下的phenanthrene_DMNAO.out，這是一個Gaussian輸出文件，對應的輸入文件是examples目錄下的phenanthrene_DMNAO.gjf。從這個.gjf文件可見，任務要求Gaussian調用NBO模塊進行分析，并且將DMNAO關鍵詞傳遞給NBO模塊來讓它把NAO為基的密度矩陣輸出出來。啟動Multiwfn，依次輸入
examples\phenanthrene_DMNAO.out
9 //鍵級分析
11 //計算AV1245。此時如屏幕提示所示，自動從輸入文件里讀取了NAO信息和以NAO為基的密度矩陣
1,2,3,4,5,6 //邊緣環上的原子序號
計算結果為10.98，和前面基于原始基函數算出來的11.17基本一樣，這是因為當前沒用彌散函數，所以差異必然不會很明顯。這也同時暗示不用彌散函數時，直接基于原始基函數算的AV1245是靠譜的，完全沒必要像原文那樣基于耗時的AIM劃分來計算。
注：基于NAO計算AV1245的時候也可以用獨立的NBO程序的輸出文件作為Multiwfn的輸入文件，不是非得依賴于Gaussian。
3 例2：分析卟啉不同離域路徑的芳香性
這一節用到的卟啉在B3LYP/6-31G*下優化得到的fch文件可以在http://www.shanxitv.org/multiwfn/extrafiles/porphyrin.rar下載。此體系的結構如下所示
計算這個體系的大環的AV1245和上一節其實一樣，按照連接順序依次輸入環上的各個原子序號即可。但是一個一個手動輸入很麻煩，對照著結構看序號還容易眼花看錯，導致結果無意義。為了便于輸入大環中的原子序號，Multiwfn有個很貼心的設計，結合GaussView（版本>=6.0）使用非常方便。首先，進入AV1245計算界面后先輸入d然后按回車進入特殊的輸入模式，然后在GaussView里選擇刷子選擇工具，即Builder - Select Atoms by Brush，然后按住鼠標左鍵滑過要研究的環上每個原子使之被選中成為黃色（如下圖所示），之后進入Tools - Atom Selection，把環上的原子序號從文本框里（見下圖）拷到Multiwfn窗口里，之后Multiwfn就會根據原子連接關系將輸入的序號自動整理成滿足連接關系的序號，然后開始AV1245的計算。
下面我們來實際算一下卟啉的上圖所選中的環的芳香性。啟動Multiwfn后輸入
porphyrin.fch
9 //鍵級分析
11 //計算AV1245
d //在此模式里，環上的原子序號可以按照任意順序輸入
1,3-4,6-8,10-14,16-19,21-22,24 //從GaussView的Atom Selection界面里拷出來的環上的原子序號
輸出信息如下
The order of the atoms in the ring has been successfully identified
Number of selected atoms: 18
Atomic sequence:
1 3 6 8 12 14 17 18 19 22 21 16
13 10 11 7 4 24
4-center electron sharing index of 1 3 8 12: 0.00297431
4-center electron sharing index of 3 6 12 14: 0.00194343
4-center electron sharing index of 6 8 14 17: 0.00287786
...略
4-center electron sharing index of 7 4 1 3: 0.00194343
4-center electron sharing index of 4 24 3 6: 0.00297431
4-center electron sharing index of 24 1 6 8: 0.00259647
AV1245 times 1000 for the selected atoms is 2.75856093
從提示可見，Multiwfn根據輸入的序號和原子連接關系正確識別出了環上的實際原子序列（疑心重的話，可以對照結構圖上的序號去驗證），AV1245乘以1000后的結果為2.758。
類似地，我們對其它幾種可能的共軛路徑也這么算，依次輸入
d
1-2,4-8,10-14,16-20,22-24
d
1,3-4,6,8-10,12-13,15-16,18-19,21-22,24
d
1-2,4-6,8-10,12-13,15-16,18-20,22-24
下圖將上面算過的四種情況的結果匯總在了一起。考察的路徑通過紅色高亮顯示了，黑色文字是在Multiwfn的“d”模式里輸入的原子序號，藍字是AV1245乘以1000后的值。
可見上圖左上角的路徑整體共軛程度最高，因此這個環的AV1245數值可以代表卟啉的芳香性并與其它體系進行對比。
在前述的http://www.shanxitv.org/432一文中，也介紹了怎么繪制卟啉分子平面上方1.2 Bohr 處的LOL-pi圖，結果如下
從LOL-pi圖中的紅色和橙色勾勒出的路徑可以看出，卟啉這個體系中的pi電子的離域路徑傾向于穿越吡咯區域的氮原子，而同時傾向于避開N-H部分的氮原子。這樣的優勢離域路徑正好也是前面我們算出來的AV1245最大的路徑，因此很好地體現出AV1245在研究大環芳香性上是頗有意義的。
需要注意的是，在“d”模式里便利地輸入序號有個限制，也就是被研究的環里每個原子不能同時與超過兩個在環上的原子相連。比如下面的體系，我們不能直接在“d”模式里輸入1-10來讓Multiwfn自動識別出整個10元環上正確的原子連接順序，因為1、9號原子同時連接了三個被考察的環上的原子。
4 總結
AV1245對于考察大環芳香性非常有用，本文的例子已經很大程度體現了其價值和合理性。在原文里有更多例子，感興趣的讀者可以參看。由于AV1245提出的年代較晚，而且之前沒有公開的程序可以實現，因此截止到目前（2019年）還極少在實際研究中被使用，可靠性還需要更多的檢驗。筆者相信，有了Multiwfn這么方便的實現AV1245計算的程序，此方法在未來應該會得到比較多的應用。
從2020-Jun-1更新的Multiwfn開始，Multiwfn又新支持了AVmin指數，這是AV1245的作者在J. Phys. Chem. C, 121, 27118 (2017)中提出的。AVmin對應于計算AV1245過程中涉及的4c-ESI的最小的絕對值。AVmin體現了整個路徑上電子共軛的“瓶頸”，而不是像AV1245那樣體現的是平均共軛程度，因此在衡量大環的芳香性上有其獨特價值。這個特點在Phys. Chem. Chem. Phys., 20, 2787 (2018)中使用各種芳香性指標考察卟啉類體系的離域路徑上有所體現。AVmin會在給出AV1245的時候同時給出，討論例子見Multiwfn手冊4.9.11節。