使用Multiwfn繪制一維NICS曲線并通過積分衡量芳香性
使用Multiwfn繪制一維NICS曲線并通過積分衡量芳香性
文/Sobereva@北京科音 2023-Aug-6
1 關于NICS曲線及其積分
NICS是一種非常流行的基于磁性質衡量芳香性的方法，在《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）里有專門的詳細介紹，不了解的話務必先閱讀此文的相關部分。NICS大多數情況是考察特定點（如環中心或環中心上/下方1埃處）的磁屏蔽情況。在J. Phys. Chem. A, 123, 3922 (2019)中作者建議從環中心開始在垂直于環方向上對NICS_ZZ（此處Z指垂直于環方向）進行掃描，并將這一維NICS_ZZ曲線進行積分，文中認為此積分值比諸如NICS(0)ZZ、NICS(1)ZZ那樣只考察個別的點在衡量芳香性上更嚴格和可靠。而且通過考察NICS曲線無疑能比只計算個別位置的NICS獲得更多信息，有助于不同體系間橫向討論和對比。在J. Phys. Chem. A, 126, 3433 (2022)中另外的研究者將這種NICS積分的分析思想應用在氨基酸側鏈苯環的芳香性研究上，并將之稱為integral NICS (INICS)指數。類似于常見的NICS(1)ZZ，垂直穿越某環的NICS_ZZ曲線積分值為明顯負值、接近0、明顯正值說明相應的環表現芳香性、非芳香性和反芳香性。
本文將通過例子介紹如何用強大且易用的Multiwfn波函數分析程序結合Gaussian量子化學程序方便地計算、繪制特定直線上的一維NICS曲線并得到其積分值，也即INICS。無論被考察的環是平行于笛卡爾平面還是歪斜、扭曲的，無論是簡單環狀體系還是結構復雜的體系，此功能都能方便地使用。可以還可以把NICS分解為特定軌道貢獻以了解更深入信息，用于比如分離研究sigma和pi電子對芳香性的貢獻，以及分離考察類似18碳環及衍生物體系（見http://www.shanxitv.org/carbon_ring.html）那樣in-plane和out-of-plane兩套pi軌道對芳香性的貢獻。
Multiwfn之前已支持ICSS分析方法并已得到極其廣泛的應用。ICSS相當于在三維立體空間中計算各個位置的NICS（注意相差正負號），見《通過Multiwfn繪制等化學屏蔽表面(ICSS)研究芳香性》（http://www.shanxitv.org/216）里的詳細介紹。此文也示例了在獲得ICSS三維格點數據后就可以輕易用Multiwfn通過插值方式計算和繪制任意方向上的NICS曲線，但它相對于本文介紹的功能有以下不足：
(1)ICSS格點數據計算相當耗時，沒有像樣服務器的話對于中等大小的體系都算不動，而本文介紹的Multiwfn的一維NICS曲線計算功能在筆記本上也能輕易對不太大的體系完成計算
(2)本文介紹的功能支持NICS分解成軌道貢獻，這是Multiwfn的ICSS模塊不支持的（其實原理上ICSS也能這么分解，只是由于技術原因無法實現：當ghost原子數較多的時候無論用Gaussian自身的AICD接口還是借助NBO 7.0的基于正則分子軌道的NCS分解都會卡住算不動）
(3)ICSS模塊不支持考察任意方向NICS分量，比如垂直于某個傾斜的環方向的分量
另外，Multiwfn還有專門繪制二維NICS圖的功能，見《使用Multiwfn巨方便地繪制二維NICS平面圖考察芳香性》（http://www.shanxitv.org/682）。
下面將結合具體例子講解一維NICS曲線的繪制和積分。讀者務必使用2023-Aug-5及以后更新的Multiwfn版本（注意看程序啟動后顯示的更新日期），Multiwfn可以在主頁http://www.shanxitv.org/multiwfn免費下載。不了解Multiwfn者建議參看《Multiwfn入門tips》（http://www.shanxitv.org/167）和《Multiwfn FAQ》（http://www.shanxitv.org/452）。本文用的Gaussian是G16。
2 例1：計算苯、C5H5-和C7H7+垂直于環的NICS曲線
Multiwfn文件包自帶的examples目錄下的NICS_scan子目錄下的benzene.pdb是在B3LYP/6-31G*級別下優化過的苯的結構，體系處在Z=0的XY平面上。此例對它繪制從環中心下方10埃到環中心上方10埃的NICS曲線，掃描路徑垂直穿越環中心，并且考察的是垂直于環方向的NICS分量（即Multiwfn會讀取Gaussian算的各個掃描點的磁屏蔽張量，取投影出的垂直于環方向的分量。本文后面說的NICS一律都是指取這個分量）。
啟動Multiwfn，然后輸入
examples\NICS_scan\benzene.pdb
25 //離域化和芳香性分析
13 //一維NICS掃描和積分
2 //通過一批原子定義環，將掃描的兩個端點置于環中心上方和下方一定距離處，且連線垂直穿越環中心
1-6 //用于1-6號碳原子定義環。之后屏幕上會顯示通過最小二乘法擬合得到的環平面的法矢量
[回車] //用環上的原子的幾何中心作為環中心。此處也可以自己輸入其它方式得到的中心坐標，比如按照http://www.shanxitv.org/108用Multiwfn做電子密度拓撲分析得到的環臨界點（rcp）
10 //一個端點位于環下方10埃處
10 //另一個端點位于環上方10埃處。之后會從屏幕上看到掃描的起點和終點分別為0,0,-10和0,0,10
200 //掃描的點數。此例相當于每隔約0.1埃一個點，足夠精細了。點數越多計算越耗時
1 //產生Gaussian輸入文件
examples\NICS_scan\template_NMR.gjf //這是一個Gaussian做NMR計算的模板文件。打開它便知計算是在B3LYP/6-31+G*下進行的，和普通輸入文件唯一的差異是坐標部分用[geometry]代替
現在Multiwfn在當前目錄下產生了NICS_1D.gjf，可以用GaussView觀看一下，如下所示，各個紫色的Bq原子是要計算NICS的各個位置（當前原子間的成鍵沒顯示出來，想顯示的話，先在File - Preference - Custom Bonding Parameters里加入Bq與C之間的成鍵判斷設置，將所有類型鍵的判斷閾值都改為0。然后再選Edit - rebond判斷成鍵即可）
使用Gaussian對NICS_1D.gjf進行計算，我跑好的文件是examples\NICS_scan\benzene_NICS_1D.out。有個不爛的機子話，此任務轉眼就能算完。
接著在Multiwfn的窗口里輸入2，然后輸入benzene_NICS_1D.out的路徑，Multiwfn就會把其中的各個Bq的磁屏蔽張量都載入了。在之后的界面里可以選1和2分別繪制和保存NICS曲線圖像，選3可以導出數據以便在Origin等程序里更靈活地繪制（導出的文本文件里每一列是什么含義在屏幕上提示得超級清楚，別無視），選4可以計算NICS曲線的積分值，選5可以給出曲線的極小和極大點位置和數值。選-1可以設置考察NICS的哪個分量（更確切來說，是指取各個掃描點的磁屏蔽張量的什么分量），一定要記得默認是考察垂直于掃描路徑上的分量，而用戶也可以改成考察自定義矢量上的分量、X或Y或Z分量、各向同性值、各向異性值。
這里選1繪制NICS曲線圖，看到下圖。圖中X=0的位置是環中心位置。可見在分子上方和下方大約1埃處對垂直于環平面方向的磁場的屏蔽最強。
同時在文本窗口也自動顯示了這個NICS曲線的積分值，為-140.78 ppm*Angstrom。由于從上圖看到曲線兩端的數值已經幾乎為0了，所以當前的積分值是準確的。通常掃平面上下、方10埃的范圍就夠了。如果體系兩側是對稱的，只掃描單方向也可以，相應地積分值也只有一半。J. Phys. Chem. A, 123, 3922 (2019)里算的NICS曲線積分值是對于單側掃描情況而言的。
如果在Multiwfn當前界面里選選項5，會輸出如下信息，說明在環平面上方和下方0.955埃處是NICS最小值位置，數值為-28.7 ppm。
Minimum X (Angstrom): -0.954774 Value: -0.28713400E+02
Minimum X (Angstrom): 0.954774 Value: -0.28713400E+02
Totally found 2 minima, 0 maxima
下面再以類似方式對C5H5-和C7H7+計算一維NICS曲線。這倆物質都是和苯分子一樣是6個pi電子的Huckel芳香性體系，但畢竟結構不同，芳香性會存在一定差異，這可以通過對比NICS曲線來考察。此二者的pdb文件、NICS掃描的輸入和輸出文件見examples\NICS_scan目錄下的C5H5-和C7H7+開頭的文件。按照上文做法產生它們的NICS掃描的輸入文件后，在Gaussian計算之前別忘了修改里面的凈電荷設置成為實際情況。
將苯、C5H5-和C7H7+的NICS曲線數據導出然后在Origin里繪制到一起，得到下圖。為了便于展現差異，只顯示了距離平面5埃范圍內的部分。
NICS積分值為C5H5- = -152.70 ppm*Angstrom，C7H7+ = -142.06 ppm*Angstrom。可見單從NICS掃描的角度來看，C5H5-的芳香性比苯和C7H7+都要強。這個結論其實并不準確，因為苯、C5H5-、C7H7+都是pi芳香性體系，但對于芳香性無貢獻的sigma電子對距離平面較近處的NICS會產生一定影響，所以只有單獨考察pi電子對NICS的貢獻才能最嚴格地對比。下一節就演示只考察pi電子貢獻的NICS曲線的繪制。
3 例2：計算垂直于苯環的NICS-sigma和NICS-pi曲線
在《基于Gaussian的NMR=CSGT任務得到各個軌道對NICS貢獻的方法》（http://www.shanxitv.org/670）中我提到，可以利用Gaussian的AICD接口來實現只考察某些軌道對磁屏蔽值貢獻的方法。對于繪制苯的pi電子貢獻的NICS曲線，只需要把上一節例子里的模板文件改為以下內容即可。此文件是examples\NICS_scan\template_NMR_benzene-pi.gjf。
#p b3lyp/6-31+G* nmr=csgt iop(10/93=2)
[空行]
template file
[空行]
0 1
[geometry]
[空行]
AICD.txt
[空行]
17,20,21
[空行]
[空行]
其中nmr=csgt iop(10/93=2)必須寫，AICD.txt是產生的給AICD程序用的文件，當前情況不用管，文件名隨意，產生后可以刪掉。17,20,21是苯的當前的結構在B3LYP/6-31+G*級別波函數的pi軌道序號。想找軌道序號的話，可以在當前級別計算個單點任務得到fch文件，然后根據《使用Multiwfn觀看分子軌道》（http://www.shanxitv.org/269）介紹的Multiwfn主功能0自行挨個觀看軌道圖形判斷，從HOMO開始依次查看序號更小的軌道，直到找到期望個數的pi軌道。也可以根據《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）里的做法讓Multiwfn自動判斷（注意判斷pi型分子軌道僅限于純平面體系，此文明確說了）。還要注意一點，NMR=CSGT關鍵詞要求Gaussian以CSGT方式計算磁屏蔽張量，和直接寫NMR關鍵詞默認用的GIAO方式的結果有一定定量差異，但定性一致。比如計算苯的NICS(1)ZZ，用CSGT算的是-27.13 ppm，GIAO算的是-28.76 ppm。所以兩種方式算的結果不適合直接橫向對比。
基于以上模板文件由Multiwfn產生的只考慮pi軌道的輸入文件是examples\NICS_scan\benzene-pi_NICS_1D.gjf，Gaussian計算產生的同名的out文件也在此目錄下提供了。注意當Bq多的時候，nmr=csgt iop(10/93=2)任務的計算耗時會增加非常多。不帶Bq的時候在筆者的2*7R32機子上96核并行4秒鐘算完，而帶了當前200個Bq之后花了39秒。
在Multiwfn里載入NICS掃描輸出文件的時候載入examples\NICS_scan\benzene-pi_NICS_1D.out，然后繪制NICS曲線圖，如下所示。可見在上一節的圖中看到的分子平面處的局部極大點沒了，說明那個極大點是sigma軌道所致，當只考慮pi軌道貢獻的話NICS曲線處處為負。
以類似方式，再得到sigma軌道貢獻的NICS曲線，使用的模板文件里軌道序號部分寫1-16,18-19。并且和上面的pi貢獻曲線繪制到一起便于比較，如下所示。可見sigma電子對距離平面很近區域貢獻明顯為正，這充分體現出用NICS(0)ZZ討論pi芳香性體系的芳香性非常不適合，sigma電子的影響會嚴重摻和進去。
以同樣的方式，對C5H5-和C7H7+獲得sigma和pi電子貢獻的NICS曲線，和苯的曲線繪制在一起便于對比，如下所示
從上圖來看，這三個體系的pi芳香性相仿佛。為了能夠更準確地對比，可以看Multiwfn輸出的NICS曲線積分值，如下所示。可見pi芳香性是C7H7+ >= benzene > C5H5-。
benzene sigma: 14.82 pi: -142.45
C5H5- sigma: -0.38 pi: -134.85
C7H7+ sigma: 14.64 pi: -145.02
值得一提的是，本文計算的結果和NICS曲線積分原文J. Phys. Chem. A, 123, 3922 (2019)里的很一致。以苯的NICS-pi曲線積分值作為100%的話，C7H7+的是101.8%，C5H5-的是94.7%，文獻里是102.3%和94.7%，而文獻里是基于GIAO-B3LYP/6-311+G(d)靠Gaussian結合NBO 7.0做的分解（此文獻表3里的絕對數值和本文的不同，一方面是計算級別不同，一方面是文獻里是從0到∞積分，而本文是[-10,10]埃范圍積分。另外，文獻里依靠NBO 7.0實現分解相對于本文做法的一個缺點是需要花錢購買NBO 7.0）。
以上的pi芳香性的結論通過其它方式也可以證明。比如用Multiwfn做ELF-pi的拓撲分析得到二分點數值，C7H7+為0.936，苯為0.897，C5H5-為0.803，也是C7H7+ >= benzene > C5H5-的順序。如果你不了解怎么做這種分析的話，看《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）。
4 例3：計算垂直于無窮烯中某苯環的NICS曲線
在《利用Multiwfn計算傾斜、扭曲環的NICS_ZZ》（http://www.shanxitv.org/261）中介紹過Multiwfn可以對扭曲、傾斜的環計算NICS_ZZ值，Multiwfn計算一維曲線的功能也同樣可以考慮這樣的環。這里以無窮烯中間的環作為例子，結構文件是examples\NICS_scan\infinitene.pdb，是在PBE0/6-31G(d)級別優化過的，如下所示，給了兩個不同的視角。此例要研究的是高亮的原子構成的環，這部分的環明顯是斜著的而且是略微扭曲的。
由上圖可見，不可能在被高亮的環上方和下方對稱地進行掃描，否則其中一頭掃描的點就會接觸到另一側的環了。因此，本例只從這個環的中心垂直向外側掃描10埃距離。
啟動Multiwfn，然后輸入
examples\NICS_scan\infinitene.pdb
25 //離域化和芳香性分析
13 //一維NICS掃描和積分
2 //通過一批原子定義環，將掃描的兩個端點置于環中心上方和下方一定距離處，且連線垂直穿越環中心
35-37,68-69,71 //用于35-37,68-69,71號碳原子定義環，即上圖高亮的那些原子
[回車] //用環上的原子的幾何中心作為環中心
10 //一個端點位于環中心下方10埃處
0 //另一個端點位于環上方0埃處，也即環中心處
[回車] //掃描的點數用推薦的100（推薦值是以間隔0.1埃確定的）
1 //產生Gaussian輸入文件
examples\NICS_scan\template_NMR.gjf //Gaussian做NMR計算的模板文件
用GaussView觀看在當前目錄下產生的NICS_1D.gjf，如下所示。可見Bq確實出現在了期望的位置。如果對某個體系你發現Bq分布在了與預期相反方向，就把環上方和環下方的距離反過來輸入重新產生輸入文件即可。
為節約時間，把NICS_1D.gjf里的6-31+G*改為6-31G*，然后用Gaussian計算之，得到的文件是examples\NICS_scan\infinitene_NICS_1D.out。之后在Multiwfn里輸入2，輸入此文件路徑，再選1繪圖，看到下面的NICS曲線圖（如果想讓此圖左右翻轉，可以選選項-3）
此曲線積分值為-96.46 ppm*Angstrom。僅對環的單側進行積分數值就已經很大了，明顯大于第2節苯分子的兩側NICS積分值的一半（約-70 ppm）。其原因絕非是因為無窮烯中這個環的芳香性比苯還強，而是附近其它的環產生的磁屏蔽效果產生了疊加。使用Multiwfn計算《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）里介紹的多中心鍵級或者《使用Multiwfn計算AV1245指數研究大環的芳香性》（http://www.shanxitv.org/519）里介紹的AV1245分析，就可以明顯看到這個環的芳香性并不如苯。
順帶一提，是否帶彌散函數對NICS曲線影響甚微，不帶彌散函數的話耗時能巨幅降低，尤其是對于無窮烯這樣較大的體系。6-31G*和6-31+G*在2*7R32 96核機子上計算當前任務分別耗時1m6s和4m37s。6-31G*和6-31+G*算的NICS積分值分別是-96.46和-96.19 ppm*Angstrom，相差微乎其微，畫出圖來對比也幾乎沒可察覺的差異。之前的例子帶了彌散函數主要是考慮到要將C5H5-陰離子體系納入對比。只是算中性體系的NICS沒太大必要帶彌散，哪怕要掃描到較遠處。想要比6-31G*更好的結果可以用3-zeta檔次的基組。
5 總結
本文介紹了Multiwfn掃描一維NICS曲線的功能，既可以直接作圖又可以給出積分值，比只討論環中心或上/下方1埃處的值獲得的信息更充分，討論芳香性時更嚴格，而且還能分解成不同軌道的貢獻。從例子可見這種分析在Multiwfn里的操作非常簡單。為了確保計算合理，建議在Multiwfn產生NICS掃描的gjf文件后放到GaussView里檢查一下Bq的位置是否確實正確，如果和預期的不符說明在Multiwfn里輸入的信息有誤。
使用Multiwfn此功能發表文章請記得按Multiwfn啟動后顯示的信息恰當引用Multiwfn的原文。